社会学、信息科学与技术、生物学和统计物理学等领域的各种系统中存在大量的关系信息, 例如, 学术论文间的引用关系, 万维网中网页间的超链接关系, 生物细胞系统中蛋白质的物理交互关系, 航空交通网络中城市间航班往来, 社交网络中人和人的交互等.人们常常采用网络对这些复杂系统中的关系信息进行建模, 称为复杂网络.随着移动互联网、普适计算^[1]、生物计算、社会感知计算^[2]、城市计算等应用技术的发展, 对包含复杂关系信息的海量数据进行分析和挖掘已经成为很多行业和领域最迫切的需求, 为复杂网络分析带来了新的机遇和挑战.

贝叶斯模型作为一种重要的统计建模方法因其具有完备的理论基础和可解释性, 已经被广泛用于对网络数据进行建模和预测.可交换性假设是采用贝叶斯模型对网络数据进行建模的重要前提, 网络数据的可交换性是指节点的出现顺序不影响统计网络模型的定义.de Finetti^[3]提出了随机变量序列的可交换性, 称作de Finetti表示理论; Aldous^[4]和Hoover^[5]提出了随机数组的可交换性, 称作Aldous-Hoover表示理论; Kallenberg^[6]提出了随机测度(采用连续时间中的随机点过程表示网络数据)的联合可交换性和分部可交换性, 称作Kallenberg表示理论. Aldous-Hoover表示理论推动了贝叶斯统计模型用于网络数据建模的发展, 很多著名的随机图模型如ER模型^[7]、SBM模型^[8]、MMSB模型^[9]、IRM模型^[10]等都可以纳入这一框架(Aldous-Hoover可交换图, 简称AHEG^[11]), 对这一类模型的性质和推理方法的研究也取得了很大的进展.然而, AHEG模型要么只能生成空图, 要么只能生成稠密图, 不能产生稀疏图^[12], 即AHEG模型不能同时满足可交换性和稀疏性.

实证结果表明, 真实世界中的很多复杂网络都具有节点度幂律分布的稀疏特征, Barabási提出的PA (preferential attachment)模型^[13]为了生成具有幂律特征的稀疏网络, 放弃了可交换性, 因而PA模型只适用于解释网络整体结构属性(如节点度分布), 不能用于进行社区发现和链路预测^[14].Bollobás等人、Wolfe等人、Borgs等人提出一个折衷的方向, 基于重新调节标度的方法来产生有限可交换的稀疏图序列, 然而这类方法生成的可交换稀疏网络却不满足可投影特性, 也就不保证随机图随节点数增加产生的图序列具有一致的分布(可投影特性是采用流方式处理网络数据的重要前提), 因而基于Aldous-Hoover表示理论的统计生成模型具有先天缺陷^[11].

Caron等人^[15]利用随机测度和随机图间的联系提出了采用可交换随机测度来对网络数据进行建模, 将网络表示成随机点过程而不是邻接矩阵, 基于Kallenberg表示理论, Caron-Fox模型既可以生成稠密可交换图也可以生成稀疏可交换图, 解决了可交换性和稀疏性二者不可兼得的矛盾.稀疏可交换图建模由此成为学术界的研究热点, 伴随大量的相关研究工作产生了很多重要的理论和方法.

本文对稀疏可交换图建模的研究进展进行综述.第1节讨论随机图模型、Graphon模型、Graphex模型、贝叶斯非参数混合模型、可交换性表示理论、泊松点过程、完全随机测度、带独立增量的归一化随机测度、可交换分区概率函数等理论基础.第2节介绍有向多图和带自边的无向图所对应的Caron-Fox模型.第3节总结进行稀疏可交换图模型的随机模拟所涉及的截断采样和边缘化采样方法的研究发展, 重点描述通过条件采样构造CRM的Adaptive Thinning算法和通过边缘化采样构建稀疏可交换图的Pólya Ura Scheme方法.第4节对Caron-Fox模型后验推理过程采用的MCMC近似推理算法进行综述, 详细阐述哈密顿MCMC算法, 并对切片采样算法和指数倾斜稳定分布的采样进行说明.第5节对稀疏可交换图建模的最新研究进展和研究前景进行介绍.

1 相关理论 1.1 随机图模型(random graph model)

随机图是利用观测数据理解复杂网络结构的重要工具, 其基本原理是:观测到的图结构是由一个潜在数据生成过程(即随机图模型)产生的, 从最大似然性思想出发, 在模型空间中寻找使观测数据具有最大出现可能性的模型, 其独特之处在于可以用概率的方法来证明所需要的图的存在性, 而不必把图真实构造出来.

定义1.1. 一个随机图模型是一族编了索引的图值(graph-valued)随机变量G_{s, ϕ}, s定义了图的尺寸(取值自一个完全有序集合S), ϕ是模型参数, 决定了图G_{s, ϕ}的一些分布特性, G_{s, ϕ}的分布记为g_{s, ϕ}^[11].

例1.1(ER模型):经典随机图模型由埃尔德什和莱利^[7], 是n个节点上的一族简单随机图G_{n, p}, n∈N, 节点间连边的概率为p∈(0, 1), 边的产生是相互独立的.可以将G_{n, p}直观地表示成一个n×n随机邻接矩阵(一个对称的存放0, 1值的n×n数组, 其对角线全部是0), 统计分析的目标之一就是利用给定观测数据推理参数p.

为了更准确地对真实网络进行建模, 随机图模型的一个重要属性就是节点数和边数之间的关系.随机图模型G_{s, ϕ}, 给定参数ϕ, 记s_n↑∞表示一个递增趋向于无穷大的尺寸参数序列.给定一个图G, ν=|v(G)|表示节点数, e=|e(G)|表示边数, 假设n→∞时ν→∞.若当n→∞时$\frac{{\sqrt e }}{\nu } \to 0$以概率1成立, 则称随机图序列(G_{s, ϕ})是稀疏的, 节点数为v的稀疏图其边数为e=o(v²).

例1.2(KEG):Kallenberg可交换图G_{u, ϕ}, 尺寸参数s不是自然数n(这是一种常用的随机图族编索引的方法), 而是一个非负实数v, $\nu \propto \sqrt e , $e是图的边数的期望值; KEG有3个可能的组成部分:隔离边I、无限星形结构S、包含图结构主要信息的有限部分Θ, 参数ϕ是一个三元组$(I, S, \mathit{\Theta } ), I \in {{\bf{R}}_ + }, S:{{\bf{R}}_ + } \to {{\bf{R}}_ + }, \mathit{\Theta } :{\bf{R}}_ + ^2 \to [0, 1]{\rm{, }}$称作graphex^[11].通常只讨论没有隔离边和无限星形结构的KEG, 所以I=S=0, 不加区分地称Θ为graphex.

统计网络分析的固定模式是:观测到的网络g_s被建模成G_{s, ϕ}的一个实例(参数s已知, ϕ未知), 网络分析的目的就是推理出ϕ.某些随机图模型中, 序列G_{s1, f}, G_{s2, ϕ}, ...是一个网络演化动力学模型, 尺寸s与采样大小有关, 随着采样到的观测数据的增多, 图的尺寸也在增大, 此时自然就必须考虑不同尺寸下的网络应该在分布上满足一致性(consistency).一致性可以通过要求随机图模型分布是可投影的(projective)来表述.可投影性定义在可投影系统(projective system)中, 即一族可测映射(f_{s, t}; s≤t∈S), f_{s, t}将尺寸为t的随机图映射到尺寸为s≤t的随机图, f_{t, t}是同等映射, f_{r, t}=f_{r, s}of_{s, t}, (r≤s≤t)^[11].

定义1.2. 若存在可投影系统f_{s, t}; s≤t∈S, 使得对于任意s < t和参数ϕ, 都有${G_{s, \phi }}\underline{\underline d} {f_{s, t}}({G_{t, \phi }}){f_{s, t}}({G_{t, \phi }})$, 则称随机图模型是可投影的, $\underline{\underline d} $表示在分布上相等^[11].

定义1.3. Graphon是定义在概率空间上的, 在[0, 1]上取值的可测对称函数Θ:[0, 1]²→[0, 1].Graphon是无限随机图序列(G_{s, ϕ})的极限^[11].

例1.3(AHEG模型):Aldous-Hoover可交换图是n个节点上的一族简单随机图G_{n, Θ}, n仍然表示节点数, 但是参数ϕ不再是概率p, 而是一个对称可测随机函数Θ:[0, 1]²→[0, 1], 称作graphon^[11].ER模型可以看做是AHEG的一个特例, 即Θ(x, y)=p是一个常graphon.

对于无限图, Aldous-Hoover表示理论断言:一个随机图具有可交换性当且仅当其分布是定义在特定的一族各态历经测度(ergodic measures)上的一个混合分布, 每一个各态历经测度就是一个graphon, 由此导致的一个必然结果是n个节点间的期望连边数为$\left( \begin{array}{l} n\\ 2 \end{array} \right){\left\| \mathit{\Theta } \right\|_1}$(${\left\| \mathit{\Theta } \right\|_1}$表示Θ的1范数), 因此AHEG生成的要么是空图, 要么是稠密图^[11].

随机图模型的一般结构:给定一个Graphon, 必定存在一个对应的具有可数无限可交换性的随机图G(n, Θ), 相应的随机邻接矩阵(G_ij)_ij∈N定义为

$ \left\{ \begin{array}{l} \quad \quad \quad \;\;\mathit{\Theta } \sim \mu \quad \quad \quad \, \\ \quad \quad \quad \;{U_i}{ \sim _{iid}}Uniform[0, 1]\quad {\rm{for}}\, i \in {\bf{N}}\\ {G_{ij}}|\mathit{\Theta } , {U_i}, {U_j}{ \sim _{ind}}Bernoulli(\mathit{\Theta } ({U_i}, {U_j})) \end{array} \right. $

(1.1)

相应的随机邻接测度${\{ ({\theta _i}, {\theta _j})\} _{{\theta _i}, {\theta _j} \in {{\bf{R}}_ + }}}$定义为

$ \left\{ \begin{array}{l} \quad \quad \quad \;\;\, \, {N_\alpha } \sim Poisson(c\alpha )\\ \quad \quad \;\{ {\theta _i}\} |{N_\alpha } \sim {\;_{ind}}Uniform[0, \alpha ]\\ \quad \quad \;\{ {\vartheta _i}\} |{N_\alpha } \sim {\;_{ind}}\;Uniform[0, 1]\\ ({\theta _i}, {\theta _j})|\mathit{\Theta } , {\vartheta _i}, {\vartheta _j} \sim {\;_{ind}}{\rm{ }}Bernoulli(\mathit{\Theta } ({\vartheta _i}, {\vartheta _j})) \end{array} \right. $

(1.2)

$ \mathit{\Theta } ({\vartheta _i}, {\vartheta _j}) = \left\{ {\begin{array}{*{20}{l}} {1 - \exp ( - 2{{\bar \rho }^{ - 1}}({\vartheta _i}){{\bar \rho }^{ - 1}}({\vartheta _j}), }&{{\rm{if}}\;{\vartheta _i} \ne {\vartheta _j}}\\ {1 - \exp ( - {{\bar \rho }^{ - 1}}{{({\vartheta _i})}^2}), }&{{\rm{if}}\;{\vartheta _i} = {\vartheta _j}} \end{array}} \right. $

(1.3)

将一个分布建模成由很多简单分布混合而成, 既是一种有用的非参数密度估计方法, 又是一种对解释变量间依赖关系的潜在类进行识别的重要方法.可以采用以某个先验分布作为混合比例的层次贝叶斯框架来处理由可数无限个成份组成的混合分布, 混合比例有助于找到起决定作用的混合成份.

非参数贝叶斯估计的灵活性能够带来更好的预测性能, 原因在于其学习能力不会饱和, 从而使得预测性能可以随着观测数据的增多而持续提升^[16].最大后验估计方法进行推理的目标只是找到某个特定的使后验最大的参数, 进行预测的时候只用一个模型; 而NPB估计方法进行推理的目标的则是学习参数的分布(即考虑所有可能的参数), 进行预测的时候把不同的模型都考虑进来, 无穷多种模型按相应的重要性权重一起发挥作用.令D表示观测样本集, Θ表示模型的所有参数, x^*表示做预测的新样本, $\hat y$表示x^*的预测值, MAP和NPB的对比见表 1.

1.2 可交换性表示理论

模型空间中存在大量可选的随机图模型, 然而大部分模型都只能建模出真实网络的一部分特性, 换个角度可能就变成了病态模型, 因此很难评估这些模型的统计可用性.想要找到既易于处理, 又足够灵活(可以较准确地解释真实世界中各种现象)的随机图模型, 就必须做出一定假设, 可交换性假设是非常重要的一种假设, 是指生成模型不依赖于观测样本出现的顺序, 或者说变换节点的标签不会改变图模型的概率分布(即节点的标签不提供图结构的任何信息)^[11].

可交换性是最常见的一种统计不变性(invariant), 也称为概率对称(probability symmetry), 具有统计不变性的分布族, 其结构可以通过各态历经分解(ergodic decomposition)来解释, 或者说采用表示理论来更直观地刻画^[11].以下首先介绍最简单的随机序列的可交换性表示理论, 然后将可交换性推广到更复杂的随机结构, 包括随机数组的可交换性表示理论、随机测度的可交换性表示理论, 并且对与可交换随机测度密切相关的可交换随机分区进行阐述.各种表示理论见表 2.

定义1.4. 一个随机变量序列, 若对序列元素下标进行任意排列不改变序列的联合分布, 称序列是可交换的(X₁, X₂, ...)$\underline{\underline d} $(X_π(1), X_π(2), ...), 对于任意的下标置换π^[11].

定理1.1(de Finetti表示理论). 随机序列(X_i)_i∈N具有可交换性当且仅当存在Χ上的一个随机概率测度Θ, 满足${X_1},{X_2} \ldots |\mathit{\Theta }\mathop :\limits^{{\rm{iid}}} \mathit{\Theta }$, 即以Θ为条件, 观测变量独立并且服从Θ-分布^[11].

从概率角度来看, Θ表示了观测数据的公共结构, 而统计推理出的条件概率Pr(X_i|Θ)则捕获了每个观测变量独有的随机性.

定义1.5. 二维随机数组X=(X_ij)_ij∈N称作可交换数组, 若满足(X_ij)(X_π(i)σ(j))(X_π(i)σ(j)), 对于任意下标置换π和σ.当π=σ时, 称作联合可交换, 否则称为分部可交换^[11].

定理1.2(Aldous-Hoover表示理论). 二维随机数组(X_ij)_ij∈N具有可交换性当且仅当存在一个随机可测函数F:[0, 1]³→$\mathbb{X}$, 满足(X_ij)(F(U_i, U_j, U_ij))(F(U_i, U_j, U_ij)), 其中, (U_i)_i∈N, (U_ij)_ij∈N是取自均匀分布Uniform[0, 1]的随机变量构成的序列和矩阵, 并且有U_ij=U_ji^[11].

定义1.6. 空间$\mathit{\mathbb{X}}$={0, 1}时, 二维随机数组X就是一个以N为节点集的随机图G的邻接矩阵, 对于无向图, X是一个对称矩阵, 若X满足定理1.2, 则称随机图G是可交换的.单一图对应联合可交换, 二部图对应分部可交换^[11].

此时, Aldous-Hoover表示理论可以进一步简化为:无向图是可交换的, 当且仅当存在一个参数对称的随机函数Θ:[0, 1]²→[0, 1], 使定理1.2满足$F\left( {{U_i}, {U_j}, {U_{ij}}} \right){\rm{ = }}\left\{ {\begin{array}{*{20}{l}} {1, }&{{\rm{if}}\, {U_{ij}} < \mathit{\Theta }({U_i}, {U_j})}\\ {0, }&{{\rm{otherwise}}} \end{array}} \right..$U_i与节点相关联, U_ij与边相关联, 将函数F分解成Θ:[0, 1]²→[0, 1]和H:[0, 1]²→$\mathbb{X}$, 使得(X_ij)$\underline{\underline d} $(F(U_i, U_j, U_ij))(H(U_ij, Θ(U_i, U_j)))(H(U_ij, Θ(U_i, U_j)))^[11].

Caron等人建立了随机图和对称点过程的对应关系:将随机图中的节点看作任意实数x, 将边看作${\bf{R}}_ + ^2$上的偶对(x, y), 就可以将随机图与${\bf{R}}_ + ^2$上的简单点过程对应起来, 若节点i和j有边相连则平面上的位置(θ_i, θ_j)或(θ_j, θ_i)处有一个点, 如图 1所示.

定义1.7. 一个简单无向图可以采用实平面${\bf{R}}_ + ^2$上的一个对称的点过程$Z = \sum\limits_{i, j} {{z_{ij}}\delta ({\theta _i}, {\theta _j})} $来表示, 节点i按出现的顺序嵌在实数轴的正半轴R₊的位置θ_i上, 并且关联了一个参数ω_i(解释为该节点的社交能力(sociability)), 若节点i和j有边相连则z_ij=1, δ(q_i, θ_j)是一个Dirac测度, 把单位质量集中在平面上(q_i, θ_j)处^[15].

定义1.8. ${\bf{R}}_ + ^2$上的点过程Z是可交换的, 当且仅当对于任意的h > 0, 在任意下标置换π和σ下, 都有:

$ (Z({A_i}, {A_j}))\underline{\underline d} (Z({A_{\pi (i)}}, {A_{\sigma (j)}}))(Z({A_{\pi (i)}}, {A_{\sigma (j)}})) $

(1.4)

这里, 区间A_i=[h(i–1), hi], i=N, 考虑任意小的区间A_i以确保节点i和j不会落入同一个区间^[15].

采用一个纯原子的简单随机测度来表示简单点过程更便于进行处理, 将点过程的每个点用随机测度中的每个原子来表示, 就得到了随机图的邻接测度(adjacent measure)表示, 即KEG(Kallenberg可交换图).不同于邻接矩阵表示的AHEG, KEG在连续空间取节点, 因此无限图G在某个有限区域截断得到受限图G_t=G(.I[0, t]²), 其节点个数是个随机量, 此外, KEG中节点至少要与1条边关联^[11].

定义1.9. ${\bf{R}}_ + ^2$上的随机测度W, 若对于R₊上任意的测度保持变换f都有$W \circ {(f \otimes f)^{ - 1}}\underline{\underline d} $W, 则称W是可交换的(⊗表示张量的内积运算)^[11].

定理1.3(Kallenberg表示理论). ${\bf{R}}_ + ^2$上的随机测度ζ是联合可交换的当且仅当几乎确定满足:

$ \xi = \sum\limits_{i, j} {f({\alpha _{}}, {\vartheta _i}, {\vartheta _j}, {\zeta _{\left\{ {i, j} \right\}}}){\delta _{{\theta _i}, {\theta _j}}}} + $

(1.5)

$ \sum\limits_{j, k} {(g(\alpha , {\vartheta _j}, {{\cal X}_{jk}}){\delta _{{\theta _j}, {\sigma _{jk}}}}\, + g'(\alpha , {\vartheta _j}, {{\cal X}_{jk}}){\delta _{{\sigma _{jk}}, {\theta _j}}})} \, + $

(1.6)

$\sum\limits_k {(l(\alpha , {\eta _k}){\delta _{{\rho _k}, {{\rho '}_k}}}\, + l'(\alpha , \, {\eta _k}){\delta _{{{\rho '}_k}, {\rho _k}}})} \, + $

(1.7)

$\sum\limits_j {(h(\alpha , \, {\vartheta _j})({\delta _{{\theta _j}}} \otimes \Lambda )\, + h'(\alpha , \, {\vartheta _j})({\delta _{{\theta _j}}} \otimes \Lambda ))} + \beta {\Lambda _D} + \gamma {\Lambda ^2}$

(1.8)

$f:{\bf{R}}_ + ^4 \to {{\bf{R}}_ + }$是可测函数, α≥0, ${{\mathit{\mathbb{N}}}_{2}}=\{\{i, j\}|(i, j)\in {{\mathbb{N}}^{2}}, $ζ_{(i, j)}是$\{i, j\}\in {{\mathit{\mathbb{N}}}_{2}}$构成的U-array(独立均匀分布随机变量组成的数组), $\{({{\theta }_{i}}, {{\vartheta }_{j}})\}$是$\bf{R}_{+}^{2}$上的独立单位率(unit-rate)泊松过程^[11].

因为邻接测度是纯原子的, 所有带Lebesgue成分的项必须测度为0, 因此式(1.8)取0, 式(1.7)对应了graphex的S部分, 式(1.6)对应了graphex的I部分, 式(1.5)对应了graphex的Θ部分.

$ f({\alpha _{}}, {\vartheta _i}, {\vartheta _j}, {\zeta _{\left\{ {i, j} \right\}}}) = \left\{ {\begin{array}{*{20}{l}} {1{\rm{, }}}&{{\zeta _{\left\{ {i, j} \right\}}} \le \mathit{\Theta} ({\vartheta _i}, {\vartheta _j})}\\ {0{\rm{, }}}&{{\rm{otherwise}}} \end{array}} \right. $

(1.9)

$ \mathit{\Theta} ({\vartheta _i}, {\vartheta _j}) = \left\{ {\begin{array}{*{20}{l}} {1 - \exp ( - 2{{\bar \rho }^{ - 1}}({\vartheta _i}){{\bar \rho }^{ - 1}}({\vartheta _j}){\rm{, }}}&{{\rm{if}}\;{\vartheta _i} \ne {\vartheta _j}}\\ {1 - \exp ( - {{\bar \rho }^{ - 1}}{{({\vartheta _i})}^2}){\rm{, }}}&{{\rm{if}}\;{\vartheta _i} = {\vartheta _j}} \end{array}} \right. $

(1.10)

1.3 泊松点过程(Poisson point process)

随机变量X是定义在样本空间Ω上的函数, 当x是X的观测值时, 存在Ω中的ω使得x=X(w).随机向量(X₁, ...X_n)是定义在样本空间Ω上的n元函数, 同时要研究更多的随机变量时, 就要引入随机过程的概念, 设T是(–∞, +∞)的子集, 则称随机变量的集合{X_t|t∈T}是随机过程, 称T为该随机过程的指标集(index set), 称{x_t|x_t=X_t(w), t∈T}为{X_t|t∈T}的一次观测, 当T=[0, +∞)或T=(–∞, +∞)时, {X_t|t∈T}的一次观测是一条曲线, 称作随机过程的一条轨迹(trajectory).泊松过程是一种时间连续、状态离散的随机过程^[17].

定义1.10(独立增量过程(dependent incremental process)). 用随机变量N(t)表示时间段[0, t]内某类事件发生的个数, 则{N(t); t≥0}是计数过程(counting process), 满足如下条件:①对t≥0, N(t)是取非负整数值的随机变量; ②对t > s≥0, N(t)≥N(s); ③对t > s≥0, N(t)–N(s)是时间段(s, t]内的事件发生数; ④ {N(t)}的轨迹是单调不减右连续阶梯函数.若对于任意n和0≤t₁ < t₂ < ... < t_n都有随机变量N(0), N(0, t₁], N(t₁, t₂), ..., N(t_n–1, t_n)相互独立, 则这个计数过程称为独立增量过程^[17].

定义1.11(平稳增量过程(stable incremental process)). 若在长度相等的时间段内, 事件发生个数的概率分布是相同的, 即对于任意s > 0, t₂ > t₁≥0随机变量N(t₁, t₂)和N(t₁+s, t₂+s)同分布, 则称计数过程为平稳增量过程^[17].

定义1.12(泊松过程(Poisson process)). 满足条件:① N(0)=0;② {N(t)}是独立增量过程; ③对任意t, s≥0都有N(s, t+s]服从参数为λt的泊松分布, 即${\rm{Pr}}(N(s, t + s] = k) = \frac{{{{(\lambda t)}^k}}}{{k!}}{{\rm{e}}^{ - \lambda t}}, {\rm{ }}k = 0, 1, ...{\rm{, }}$称计数过程为泊松过程^[17].

常数λ取正值, 称为泊松过程的强度(intensity), 若λ不随时间变化, 则泊松过程是齐次的, 若λ随时间变化(或者说与时间有关), 则泊松过程是非齐次的^[17].

定理1.4(叠加原理). 设Π₁, Π₂, ...是可数个相互独立的泊松点过程, 对应的计数过程分别为N₁, N₂, ..., 强度分别为μ₁, μ₂, ..., 则$\mathit{\Pi } = \bigcup\nolimits_{i = 1}^\infty {{\mathit{\Pi } _i}} $也是一个泊松过程, 对应的计数过程为$N = \sum\limits_{i = 1}^\infty {{N_i}} , $强度为$\mu = \sum\limits_{i = 1}^\infty {{\mu _i}} $^[17].

定理1.5. 令Φ₁, Φ₂, ...是一个率为λ的泊松过程中事件发生的时间(或者位置), 给定N(t)=n, 随机变量Φ₁, Φ₂, ..., Φ_n的联合概率密度函数为${f_{{\mathit{\Phi} _{{{\kern 1pt} _1}}}, {\mathit{\Phi} _{_{{\kern 1pt} 2}}}, ..., {\mathit{\Phi} _{_{{\kern 1pt} n}}}|N(t) = n}}({\mathit{\Phi} _{{\rm{ }}1}}, {\mathit{\Phi} _{{\rm{ }}2}}, ..., {\mathit{\Phi} _{{\rm{ }}n}}) = n!{t^{ - n}}$^[17].

定理1.5可以通俗地解释为给定区间内事件的总数量, 事件发生的位置是某种方式的均匀分布.

定理1.6(映射原理(mapping proposition)). $(\mathit{\mathbb{S}}, \mathcal{S})$是一个可测空间, $\mathit{\mathbb{S}}$上的随机点集Π及其计数过程N是一个强度为μ的泊松点过程, $(\mathit{\mathbb{T}}, \mathcal{T})$是另外一个可测空间, $f:\mathit{\mathbb{S}}\to \mathit{\mathbb{T}}$是一个可测函数, 若集值函数μ与f的反函数复合$\mu \, \circ \, {f^{ - 1}}$是非原子的(non-atomic), 则f(Π)={f(x):x∈Π}也是一个泊松点过程, 强度为$\mu \, \circ \, {f^{ - 1}}$^[17].

1.4 几类常用的离散非参数先验

任何贝叶斯模型都有一个重要的组成部分, 被表示成一个随机可测函数Θ, Θ可以用于构建贝叶斯模型的自然参数, 实现了对数据的解耦(或者说建立起数据间的联系)^[11].de Finetti理论刻画了与随机序列相关的模型参数Θ, Aldous-Hoover表示理论刻画了与随机数组相关的模型参数Θ, Kallenberg表示理论刻画了与随机测度相关的模型参数Θ^[18].

通过为可测函数指定先验Q, 就可以为可交换图模型指定先验.Lijoi等人2013年综述了几类非常适合作为混合模型的非参数先验Q的概率模型^[19], 包括:完全随机测度(completely random measure, 简称CRM)、带独立增量的规一化随机测度(normalized random measures with independent increment, 简称NRMI)、吉布斯型先验(Gibbs-Type prior)等.很多先验都是基于CRM对进行适当的变换得到的, CRM的一个重要特性是其几乎确定的离散性, 对CRM进行变换得到的随机测度也继承了这一特性, 因此它们以概率1选择了离散分布.应用这些先验结构的方式主要有两种:①直接用于对观测数据进行建模, 当已知这些数据是由一个离散分布产生的; ②若数据是由连续分布(或者更复杂的混合分布)产生的, 先验结构就是层次混合模型的某个的构造块.

1.4.1 CRMs(完全随机测度)

定义1.13. $(\mathit{\mathbb{T}}, \mathcal{T})$是一个完全并可分的度量空间(complete and separable metric space), $\mathcal{T}$是相应的Borel σ-代数, W是$\mathit{\mathbb{T}}$上的随机测度, A₁, A₂, ...是$\mathcal{T}$中任意可数个互不相交的可测集, 若随机变量W(A₁), W(A₂), ...相互独立并且满足W($W\left( {\bigcup\limits_i {{A_j}} \;W\left( {\bigcup\limits_i {{A_j}} } \right)} \right. = \sum\limits_j {W\left( {{A_j}} \right)} $则称W是一个CRM^[20].

Kingman^[20]证明了一个CRM W可以分解成3个独立部分之和:不含原子的非随机的测度W_d、不含原子的非负随机质量全体W_f, 含固定原子的非负随机质量全体W_r, 并且对于不包含固定原子部分W_r和确定部分W_d的任意CRM W, 存在一个泊松过程(随机点集Π=(w_i, θ_i)_i∈N及其计数过程N(dw, dθ))使得$W(\text{d}\theta )=\int_{{{\bf R}_{+}}}{wN(\text{d}w, \text{d}\theta )}$(一个CRM可以表示成泊松随机测度的线性函数)^[17], W可以等价地表示为$W=\sum\limits_{i=1}^{\infty }{{{w}_{i}}{{\delta }_{{{\theta }_{i}}}}}\text{, }$随机质量ω_i表示某种含义的权重(如, 节点的社交能力), Dirac测度${{\delta }_{{{\theta }_{i}}}}$将单位质量集中在θ_i处.

依据Campbell定理^[20], 一个CRM由其Laplace函数唯一刻画, W(A)的Laplace变换由以下公式给出:$\mathbb{E}\left[ \exp (-tW(A)) \right]=\text{exp}\left( -\int_{{{\bf{R}}_{+}}\times \mathbb{A}}{[1-\text{exp}(-tw)]\nu (\text{d}\omega , \text{d}\theta )} \right)\text{, }t\ge 0\text{, }$$\nu\left( {{\rm{d}}\omega ,{\rm{d}}\theta } \right)$是泊松过程的列维强度(Lévy intensity); 拉普拉斯指数(Laplace exponent)定义为$\psi (t)=\int_{0}^{\infty }{(1-{{\text{e}}^{-wt}})\rho (w)\text{d}w}.$$\nu$包含了点过程的所有位置(location)和跳跃(jump)的全部信息, 若$\nu \left( {{\rm{d}}\omega ,{\rm{d}}\theta } \right) = \rho \left( {{\rm{d}}\omega } \right){\mu _0}\left( {{\rm{d}}\theta } \right)$, 随机质量与随机位置相互独立, 则这个泊松过程是齐次的; 若$\nu \left( {{\rm{d}}\omega ,{\rm{d}}\theta } \right) = \rho \left( {{\rm{d}}\omega {\rm{|d}}\theta } \right){\mu _0}\left( {{\rm{d}}\theta } \right)$, 随机质量依赖于随机位置, 则这个泊松过程是非齐次的(inhomogeneous)^[20].

1.4.2 NRMIs(带独立增量的规一化随机测度)

定义1.14. W是空间$(\mathit{\mathbb{T}}, \mathcal{T})$上的CRM, $T=W(\mathit{\mathbb{T}})$表示总的随机质量.并且几乎确定满足0 < T < ∞, 则随机概率测度$\widetilde {W}=\frac{W}{T}=\sum\limits_{i=1}^{\infty }{{{{\widetilde{w}}}_{i}}}{{\delta }_{{{\theta }_{i}}}}$被定义为带独立增量的规一化随机测度^[20].

为确保定义1.14中的规一化是良定的, T必须几乎确定是有限正值, 这一点可以由列维测度ρ的性质来保证

$ \int_{{{\mathbb{R}}^{+}}}{\rho (\text{d}w)=+\infty }, \int_{{{\bf R}^{+}}}{(1-{{\text{e}}^{-\omega }})\rho (\text{d}w) < +\infty } $

(1.11)

满足(1.11)表明CRM在任意区间[0, S]有无限多个跳(对于任意的T < ∞), 称这样的CRM为无限活动CRM(infinite activity CRM)^[20].

例1.4:GP(伽马过程)的列维强度有如下形式ρ_α(dω)μ₀(dθ)=aω^–1e^–ωdωμ₀(dθ), α > 0, 记GP为W_α, 其总质量为T_α, 列维测度ρ_α满足(1.11), 从而有${{\widetilde{W}}_{\alpha }}=\frac{{{W}_{\alpha }}}{{{T}_{\alpha }}}$是一个良定的随机概率测度, 我们称这个NRMI为一个以α作为集中参数(concentration parameter), 以μ₀作为基分布的DP(狄利赫里过程)^[11].

例1.5:广义伽马过程的列维强度如下形式${{\rho }_{\alpha , \sigma , \tau }}\left( \text{d}\omega \right){{\mu }_{0}}(\text{d}\theta )=\frac{\alpha }{\Gamma (1-\sigma )}{{\omega }^{-\sigma -1}}{{\text{e}}^{-\tau \omega }}\text{d}\omega {{\mu }_{0}}(\text{d}\theta ), $α > 0, σ∈(0, 1), τ≥0, 记GPP为W_{α, σ, τ}, 其总质量为T_{α, σ, τ}, 列维测度ρ_{α, σ, τ}满足(1.11), 从而有${{\widetilde{W}}_{_{\alpha , \sigma , \tau }}}=\frac{{{W}_{\alpha , \sigma , \tau }}}{{{T}_{\alpha , \sigma , \tau }}}$是一个良定的随机概率测度, 我们称这个NRMI是一个具有参数(α, s, τ), 以μ₀作为基分布的归一化广义伽马过程(NGPP)^[15].

1.4.3 EPPF(可交换分区概率函数)

NRMI的离散分布特性很自然地引发了对生成数据的划分结构的分析, 事实上, 给定取值于某个度量空间$(\mathit{\mathbb{T}}, \mathcal{T})$上的n个观测量X=(X₁, ..., X_n), 数据的离散分布意味着数据间应该是有联系的(ties), 因此Pr(X_i=X_j) > 0, (i≠j), 即X取了k≤n个不同的值, 由此, 可以用[n]:={1, 2, ..., n}的随机分区π来作为X的一个等价的表示方式^[10], 可交换分区可以帮助我们更有效地构建关于NRMI的后验分析方法.

定义1.15. p是[n]的一族随机子集, 下标i和j同属于某个子集c当且仅当X_i=X_j, 则称π是[n]的随机分区(random partition)^[11].

如将中任意子集A_i对应的唯一取值记作$X_{{{A}_{i}}}^{*}$(例如, 当且仅当${{X}_{2}}={{X}_{3}}={{X}_{6}}={{X}_{9}}=X_{{{A}_{1}}}^{*}, $有A₁={2, 3, 6, 9}), 在混合模型中, 产生随机分区的过程可以看做是将观测变量指派到所隶属的成分, 即聚类或者社区划分随机变量X采样自可交换序列(X₁, X₂, ...), 相应的随机分区也是可交换的, 其概率密度函数只依赖于分区数量|π|=k和每个子集的大小n_i=|A_i|$(1\le i\le k, {{n}_{i}}\ge 1, \sum\nolimits_{i=1}^{k}{{{n}_{i}}}=n).$

定义1.16. π={A₁, ..., A_k}是[n]的某个随机分区的概率$\Pr (\pi = \{ {A_1}, ..., {A_k}\} ) = \prod\nolimits_n^k {({n_1}, ..., {n_k})} , $$\prod\nolimits_n^k {\mathit{\prod} _n^k({n_1}, ..., {n_k})} $是一个对称函数且满足加性规则(addition rule)$\prod\nolimits_n^k {({n_1}, ..., {n_k})} = \prod\nolimits_{n + 1}^{k + 1} {({n_1}, ..., {n_k}, 1)} + \sum\limits_{i = 1}^k {\prod\nolimits_{n + 1}^{k + 1} {({n_1}, ..., {n_i} + 1, ..., {n_k})} } , $这样的函数称为EPPF^[11].

2 稀疏可交换图建模

Caron和Fox利用随机测度和随机图间的联系提出了采用可交换随机测度来对网络数据进行建模, 将网络表示成随机点过程而不是邻接矩阵, 基于Kallenberg表示理论, Caron-Fox模型既可以生成稠密可交换网络也可以生成稀疏可交换图.以下介绍Caron和Fox提出的有向多图、带自边的无向图所对应的稀疏可交换图模型^[15].

2.1 有向多图(directed multigraphs)

可数无限多个节点的集合V=(θ₁, θ₂, ...), θ_i∈R₊, 其上的有向多图(两个节点θ_i, θ_j间可以有n_ij≥1条边, 并且允许节点θ_i自己到自己有边, 边是有方向的)可以表示成R₊上原子随机测度D^[15].

$D = \sum\limits_{i = 1}^\infty {} \sum\limits_{j = 1}^\infty {} {n_{ij}}\, \delta ({\theta _i}, {\theta _j})$

(2.1)

与每个节点θ_i相关联的参数ω_i > 0(sociability)由一个原子随机测度W定义, W的分布是一个列维强度为ρ(dω)μ₀(dθ)的齐次CRM^[15].

$ W = \sum\limits_{i = 1}^\infty {{w_i}{\delta _{{\theta _i}}}} , {\rm{ }}W \sim CRM(\rho , {\mu _0}) $

(2.2)

给定W, 则D由一个强度为乘积测度(product measure)$\widetilde{W} = W \times W$的泊松过程(Poisson process, 简称PP)生成

$D|W \sim PP(W \times W)$

(2.3)

本质上, PP是随机点集Π及其计数过程N, 式(2.1)可以非正式地解释为:两个节点θ_i, θ_j间有边, 则位置(q_i, θ_j)或(q_j, θ_i)是Π中的点, 每个点的计数n_ij由泊松分布Poisson(w_iw_j)得到.

2.2 带自边的无向图(undirected graphs with self edge)

对有向图$D = \sum\limits_{i = 1}^\infty {} \sum\limits_{j = 1}^\infty {} {n_{ij}}\, \delta ({\theta _i}, {\theta _j})$做一个简单的变换就可以得到定义1.7中的无向图$Z = \sum\limits_{i, j} {{z_{ij}}\delta ({\theta _i}, {\theta _j})} , $可以令z_ij=z_ji=1若有n_ij+n_ji > 0, 否则令z_ij=z_ji=0, 也就是说, 如果有向图两个节点θ_i, θ_j间只要有1条有向边, 那么对应的无向图中θ_i, θ_j间就会有1条无向边, 因此有向多图对应的无向图中允许节点自己和自己连边(一个现实的例子是在社交网络中, 一个用户对自己的博文发表评论).无向图的构建过程可用层次结构模型(式2.4)来表示^[15].

$ \left\{ \begin{array}{l} W = \sum\limits_{i = 1}^\infty {{w_i}{\delta _{{\theta _i}}}} \quad \quad \quad \quad \quad \quad \;\quad W \sim CRM(\rho , {\mu _0})\\ D = \sum\limits_{i = 1}^\infty {} \sum\limits_{j = 1}^\infty {} {n_{ij}}\, \delta ({\theta _i}, {\theta _j})\quad \quad \quad D|W \sim PP(W \times W)\\ Z = \sum\limits_{i, j} {\min ({n_{ij}} + {n_{ji}}, 1)\, \delta ({\theta _i}, {\theta _j})} \end{array} \right. $

(2.4)

因为Pr(z_ij=1|w)=Pr(n_ij+n_ji > 0|w)并且n_ij和n_ji是相互独立的随机变量(给定W), 由泊松过程的叠加原理(superposition principle), 两个强度是w_iw_j的泊松随机变量相加, 得到强度是2w_iw_j的泊松随机变量.又因为Pr(n_ij+n_ji > 0|w)=1–Pr(n_ij+n_ji=0|w), 所以当给出节点的社交能力w=(w_i)后, 也可以采用式(2.5)直接定义出无向图, 可以看到式(2.5)和式(2.4)的构建过程是等价的^[15].

$ {\rm{Pr}}({z_{ij}} = 1|w) = \left\{ {\begin{array}{*{20}{l}} {1 - \exp ( - 2{w_i}{w_j}), }&{i \ne j}\\ {1 - \exp ( - w_i^2), }&{i = j} \end{array}} \right. $

(2.5)

2.3 GPP(广义伽马过程)

Caron等人证明了当W是一个GPP, 并且当σ≥0时, 可以得到稀疏网络.GGP具很多有非常好的特性, 比如其参数具有非常好的可解释性, GGP是条件共轭的.Hougaard^[21]和Caron等人^[22]对GGP进行了研究.

GGP的定义在例2.5给出, σ≥0和σ < 0的情况下, GGP具有显著不同的特性, 当σ < 0时, GGP是一个有限活动CRM(finite activity CRM), 在区间[0, α]上的跳跃的个数J依概率1有限, 并且J的分布是一个率为$ - \frac{\alpha }{\sigma }{\tau ^\sigma }$的泊松分布, 而跳跃w_i独立同分布于伽马分布Gamma(–σ, t).当σ≥0时, GGP是一个无限活动CRM, GGP在任何区间[s, t]上的跳跃的个数J是无限的, 这里有几个常见的特例:s=0, τ > 0时对应GP(伽马过程); σ∈(0, 1), τ=0时对应稳定过程(stable process); $\sigma = \frac{1}{2}, \tau > 0$时对应逆高斯过程(inverse-Gaussian process).GGP的尾列维强度定义为$\bar \rho (w) = \int_{{\kern 1pt} x}^{{\kern 1pt} \infty } {\frac{1}{{\Gamma (1 - \sigma )}}{\omega ^{ - \sigma - 1}}{{\rm{e}}^{ - \tau \omega }}{\rm{d}}\omega } = \left\{ {\begin{array}{*{20}{l}} {\frac{{{\tau ^\sigma }\Gamma ( - \sigma , \tau x)}}{{\Gamma (1 - \sigma )}}, }&{\tau > 0}\\ {\frac{{{x^{ - \sigma }}}}{{\Gamma (1 - \sigma )\sigma }}, }&{\tau = 0} \end{array}} \right., $$\Gamma (\alpha , x)$是一个不完全伽马函数^[15].

3 稀疏可交换图的随机模拟

稀疏可交换图的生成模型是一个非参数贝叶斯模型, 非参数贝叶斯有无限多个参数, 使得模型的复杂性随着采样尺寸的增大而增大, 由于基于随机模拟的方法需要在有限维参数空间进行采样, 因而需要通过边缘化或者截断方式对无限维参数空间进行采样.Caron等人提出了采用截断方式对GPP进行模拟, 采用边缘化对稀疏可交换图进行模拟^[15].

3.1 对GGP进行随机模拟

截断方法也称条件采样方法(conditional sampler), 通过合适的方法采样W的有限足够多个原子, 采用一个有限维近似来替代无限维先验.这一类方法最早由Muliere等人^[23]提出, 他们证明了截断近似与无限维先验间的误差(采用总方差范式total varation norm表示)可以被选择小于某个特定值.Ishwaran等人在采用截棍方式构建NRMI $\widetilde{u}$方面做了大量工作, 他们^[24]提出了适用于DP先验模型的截断方法.Ishwaran等人^[25]研究了一类更广泛的截棍先验模型的截断采样方法, 用总方差范式表示截断误差, 并提出了一种简单的块Gibbs采样后验推理方法.Papaspiliopoulos等人^[26]提出的MH采样通过移动相互改变一对原子来加速混合; Walker^[27]提出了DP先验模型的切片采样方法; Kalli等人^[28]提出了NRMI先验模型的切片采样方法, 通过使用自然无序表示避免了参数的弱可识; Favaro等人^[29]提出了σ-稳定Poisson-Kingman先验模型的切片采样方法.在切片采样中, 辅助变量被引入到后验分布中, 使得Gibbs采样的所有完全条件具有有限维分布.重要的是, 这一类随机截断方法可以对后验分布进行精确采样, 避免了截断误差.然而, 尚不清楚怎样将这类方法推广到更广泛的非参数先验模型.Orbanz等人提出了采用单位率泊松过程构建NRMI的逆列维方法(inverse levy method).

3.1.1 通过条件采样构造CRM的Adaptive Thinning算法

NRMI作为先验时, 对应的CRM W通常只有不含原子的非负随机质量全体W_f, 只有当给出观测数据后才会将含固定原子的非负随机质量全体W_r引入到后验中.对于W_f, 因为$\widetilde{W}$是齐次NRMI, 所以随机位置独立同分布采样自μ₀, 随机质量与[S, ∞)上的具有指数倾斜强度测度ρ′(ds)=e^–Usρ(ds)的泊松随机测度同分布, 可以采用Adaptive Thinning方法进行采样^[30], 见算法1.

算法1. Adaptive Thinning.

输入:强度为$\nu '$的泊松随机测度, 截断级别S.

输出:对强度为$\nu '$的泊松随机测度在[S, ∞)进行的有限采样N.

算法步骤:

1.  令N:=Ø, t:=S

2.  迭代一些操作一直到结束

    a)  从参数是1的指数分布采样得到r

    b)  若r > W_t(∞), 结束; 否则令$t': = W_t^{ - 1}(r)$

    c)  以概率$\nu '$(t′)/w_t(t′)接受采样, 即令N:=N∪(t′)

    d)  令t:=t′, 进行下一轮迭代

3.  返回N(N就是对强度为[S,∞)的泊松随机测度在[S, ∞)进行的有限采样)

Thinning是从一种泊松随机测度进行采样的方法, 分为两步:首先从一个提议分布(一个比目标分布强度更高的泊松随机测度)采样出一些点, 然后以提议分布和目标分布的强度之比作为概率接受或拒绝每个采样^[31].

如图 2所示, 在Adaptive Thinning算法中, 从提议分布中采样点时, 从截断级别S出发从左向右迭代进行, 令$\nu '(s)$是ρ′(ds)在Lebesgue测度下的密度, 对于任意的t∈R₊, 存在一个函数w_t(s)满足${w_t}(t) = \nu '(t)$和w_t(s)≥w_t′ (s)≥$\nu '(s)$(对于任意的s, t′≥t), 对于NGGP, $\nu '(s) = \frac{\alpha }{{\Gamma (1 - \sigma )}}{s^{ - 1 - \sigma }}{{\rm{e}}^{ - s(\tau + U)}}, $每次采样到一个点t时, 引入一族渐进界限$\frac{\alpha }{{\Gamma (1 - \sigma )}}{t^{ - 1 - \sigma }}{{\rm{e}}^{ - s(\tau + U)}}.$令得到的界限w_t作为采样下一个点提议分布的强度, 随着t′的增加, 界限被收紧, 因此拒绝率不断被减小, 注意w_t(s)和${W_t}(s) = \int_{{\kern 1pt} t}^{{\kern 1pt} s} {{w_t}(s'){\rm{d}}} s'$的逆函数都可解析得到, 通过对w_t(s)其求积分然后求逆得到$W_t^{ - 1}(r) = t - \frac{1}{{\tau + U}}\log \left( {1 - \frac{{r(\tau + U)\Gamma (1 - \sigma )}}{{a{t^{ - 1 - \sigma }}{{\rm{e}}^{ - t(\tau + U)}}}}} \right), $并且有$\int_{{\kern 1pt} t}^{{\kern 1pt} \infty } {{w_t}(s'){\rm{d}}} s' < \infty , $所以, 当有限数量的点被采样后迭代将结束^[30].算法的描述如下.

3.1.2 通过单位率泊松过程构造CRM的逆列维方法

随机变量序列{ω_i}和{q_i}可以由截棍(stick breaking)过程构造^[32], 也可以由单位率泊松过程构造^[33].逆列维方法利用泊松随机测度的映射理论按从大到小的顺序对随机质量进行采样.对于R₊上的任意单调函数f(x), 记f^–1表示f(x)的右连续逆函数, 因此${f_{^{ - 1}}}(y) = \left\{ \begin{array}{l} Inf\{ x|f(x) \ge y\} {\rm{, }}\quad {\rm{若}}f{\rm{是个非减函数}}\\ Inf\{ x|f(x) \le y\} {\rm{, }}\quad {\rm{若}}f{\rm{是个非增函数}} \end{array} \right.{\rm{.}}$

定理3.1(CRMs的泊松采样). W是一个CRM, 列维强度为$\nu $(dω, dθ), 尾列维强度(tail Lévy intensity)定义为$\bar \rho (w) = \int_{{\kern 1pt} x}^{{\kern 1pt} \infty } {\rho (w){\rm{d}}w} {\rm{, }}$则W=W_f+W_r$\sum\limits_i {{J_i}\delta ({\theta _i})} + \sum\limits_k {{{\bar \rho }^{ - 1}}({\vartheta _k})\delta ({\theta _k})} \sum\limits_i {{J_i}\delta ({\theta _i})} + \sum\limits_k {{{\bar \rho }^{ - 1}}({\vartheta _k})\delta ({\theta _k})} {\rm{, }}$这里, {θ₁, θ₂, ...}是固定的跳跃位置^[34].

给定${\bf{R}}_ + ^2$上一个单位率泊松过程(θ_i, ${\vartheta _i}$), 构造一个列维强度为ρ(dω)μ₀(dθ)的CRM的逆列维方法是直接令${w_i} = {\bar \rho ^{ - 1}}({\vartheta _i}), $$\bar \rho (x)$是尾列维强度, ${\bar \rho ^{ - 1}}$是一个单调函数, 称逆列维强度(inverse Lévy intensity)^[34].

3.2 对稀疏可交换图进行模拟

NRMI作为先验时, 若对应的CRM采用截棍过程表示时, 网络的生成模型为(2.4), 若对应的CRM采用unit-rate poisson表示时, 网络的生成模型为(1.2).以下介绍Caron和Fox提出的对生成模型为(2.4)的稀疏可交换图所进行的模拟^[15].

3.2.1 将图限制在某个区域(graph restriction)

采用邻接测度表示的随机图是一个无限图(因为W(R₊)=∞), 而实际应用中所研究的网络都是有有限数量的边, 但是如何表达或者界定这个有限量是不知道的, 所以只能考虑将D或Z限制在某个区域[0, α]²得到相应的有限图D_a或Z_α, a作为一个模型参数通过推理得到, 相应的有CRM W_α和$\mu _0^\alpha $, 记$Z_\alpha ^* = {Z_\alpha }([0, \alpha ])$为[0, α]²上的总质量, 相应的有$D_\alpha ^* = {D_\alpha }([0, \alpha ])$和$W_\alpha ^* = {W_\alpha }([0, \alpha ]), $令${\bar W_\alpha } = \frac{{{W_\alpha }}}{{W_\alpha ^*}}, $则有限有向图的模型为

$ \text{For}~k=1,...,D_{\alpha }^{*}~\text{and}~j=1,2\\\left\{ \begin{align} & {{W}_{\alpha }}\tilde{\ }CRM(\rho ,\mu _{0}^{\alpha }) \\ & D_{\alpha }^{*}|W_{\alpha }^{*}\tilde{\ }Poisson(W{{_{\alpha }^{*}}^{2}})\left( \text{这里用到了泊松过程的叠加原理,定理1.4} \right) \\ & {{U}_{kj}}|W_{\alpha }^{*}{{\tilde{\ }}_{ind}}{{{\bar{W}}}_{\alpha }}\left( \text{这里用到了定理1.5} \right) \\ & {{D}_{\alpha }}=\sum\nolimits_{k=1}^{D_{\alpha }^{*}}{\delta ({{U}_{k1}},{{U}_{k2}})} \\ \end{align} \right. $

(3.1)

随机变量U_kj∈R₊对应图中的节点, 节点偶对(U_k1, U_k2)对应U_k1到U_k2的一条边, 有向边的数量$D_\alpha ^*$依赖于CRM的质量$W_\alpha ^*$, 对于每一条有向边, 所关联的节点U_kj采样自${\bar W_\alpha }$, 因为${\bar W_\alpha }$是NRMI(以概率1离散分布), 所以U_kj可以取${N_\alpha } \le 2D_\alpha ^*$个不同的值(N_α对应网络中度至少是1的节点的数量).

3.2.2 有限维生成过程(finite dimensional generative process)

由于W是一个无限活动, 所以不能直接对${W_\alpha } \sim CRM(\rho , \mu _0^\alpha )$进行采样, 通过引入罐子模型(Pólya ura scheme, 简称PUS)^[35]得到一个有限维过程, 方法如下:令$({U'_1}, ..., {U'_{2D_\alpha ^*}}) = ({U_{11}}, {U_{12}}, ..., {U_{D_\alpha ^*1}}, {U_{D_\alpha ^*2}})$, $t = W_\alpha ^*$, 对于GPP, 可以积分得到${\bar W_\alpha }$并且能够推理出${U'_{n + 1}}$在给定$(W_\alpha ^*, {U'_1}, ..., {U'_n})$后的条件分布, 因为${\bar W_\alpha }$以概率1离散分布, 所以${U'_1}, ..., {U'_n}$取k < n个不同的值${\bar U'_i}\;(1 \le i \le k), $相应的重数为$1 \le {m_i} \le n$(即有${U'_1}, ..., {U'_n}$中有m_i个取了同一个值${\bar U'_i}\;$), 利用可交换分区表示, 可以得到[n]的某个随机分区π={A₁, ..., A_k}, 其EPPF为$\prod\nolimits_n^k {({m_1}, ..., {m_k}|t)} , $给定$(t, {U'_1}, ..., {U'_n})$后${U'_{n + 1}}$的预测分布可以由EPPF得到${U'_{n + 1}}|(t, {U'_1}, ..., {U'_n}) \sim \frac{{\prod\nolimits_{n + 1}^{k + 1} {({m_1}, ..., {m_k}, 1|t)} }}{{\prod\nolimits_n^k {({m_1}, ..., {m_k}|t)} }}\frac{1}{\alpha }\mu _0^\alpha + \sum\limits_{i = 1}^k {\frac{{\prod\nolimits_{n + 1}^k {({m_1}, ..., {m_i} + 1, ..., {m_k}, 1|t)} }}{{\prod\nolimits_n^k {({m_1}, ..., {m_k}|t)} }}} {\delta _{{{\bar U'}_i}}}, $采用这种PUS表示, 我们可以将生成模型(3.1)重写为

$ \text{For}\ k=1,...,D_{\alpha }^{*}\ \text{and}\ j=1,2\left\{ \begin{align} & {{W}_{\alpha }}\tilde{\ }{{P}_{W_{\alpha }^{*}}} \\ & D_{\alpha }^{*}|W_{\alpha }^{*}\tilde{\ }Poisson(W{{_{\alpha }^{*}}^{2}}) \\ & {{U}_{kj}}|W_{\alpha }^{*}{{\tilde{\ }}_{ind}}\text{Urn process}. \\ & {{D}_{\alpha }}=\sum\nolimits_{k=1}^{D_{\alpha }^{*}}{\delta ({{U}_{k1}},{{U}_{k2}})} \\ \end{align} \right. $

${P_{W_\alpha ^*}}$是$W_\alpha ^*$的分布, 在σ≥0的情况下, ${P_{W_\alpha ^*}}$是一个指数倾斜的稳定分布, 可以对${P_{W_\alpha ^*}}$进行精确采样并且能够估计出EPPF, 就可以对图模型进行准确的采样.当$W_\alpha ^*$不能被精确采样时, 在某些列维强度可以借助于Adaptive Thinning, 其他情况下采用逆列维方法.

边缘化采样(marginal sampler)通过求W的边缘分布移除了由于W是无限活动带来的麻烦, 适用于预测分布体系显式给定的情形.边缘化得到的是先验的PUS表示, 这种表示可以用来定义有效的MCMC算法.MacEachern和Neal对在DPM模型上基于PUS表示的后验推理算法进行了总结, Favaro等人将这一类算法推广到了NRMI作为先验的混合模型, 这些方法的不足在于只有某些先验存在合适的PUS, 很多先验的PUS难以得到^[36].

4 稀疏可交换图模型的推理

稀疏可交换图的生成模型是一个以NRMI作为先验的非参数贝叶斯混合模型, 其后验分布的计算涉及到很复杂的积分, 大多不可能精确计算得到, 需要采用近似计算方法^[37], 比如马尔可夫链蒙特卡洛(Markov Chain Monte Carlo, 简称MCMC)采样方法.常用的MCMC方法包括:Gibbs采样、Metropolis-Hastings采样(MH)、切片采样、哈密顿蒙特卡罗(Hamiltonian Monte Carlo, 简称HMC)、朗之万动力学方法(Langevin dynamics).

4.1 切片采样算法

无限维先验不允许直接使用基于随机模拟的方法来进行后验分布推理, 因此需要对W进行截断, 仅对有限个原子进行处理.对给出X和U时CRM W的后验进行推理时, 切片采样算法引入一个辅助变量S_i, 称为切片变量, 其条件分布为S_i|X_i, W:Uniform(0, W({X_i}), W({X_i}表示原子X_i在W中的质量, 以S_i为条件, X_i的取值对应到在W中的原子的质量肯定是大于等于S_i的, 因此如果只需要有限个质量大于等于S_i的原子就可以更新X_i时, S_i实际上就作为了对W的截断级别^[36].也就是说, 在整个数据集上, 采样的状态空间包括X, 切片变量S, CRM W和辅助变量U, 采样方法其实就是一个Gibbs采样, 迭代更新X, 然后更新U, 然后同时更新W和S.

这里只讨论对W和S的更新, W的后验既包含W_f, 又包含由观测数据引入的W_r部分, 切片变量只依赖于W中固定原子的质量, 另外我们只需要W中有限个质量大于全局截断级别S=min_i∈[n]S_i的位置随机的原子, 因此对W和S进行充分采样的方法是先采样W_r, 然后采样S, 最后采样W_f.对于W_r, 定理4.2表明每个固定原子对应于X中的一个唯一值$\{ X_k^*:k \in \pi \} $, 其质量相互独立且与位置无关, 条件分布为${\rm{Pr}}({J'_k} \in {\rm{d}}s|U, X) \propto {s^{\left| k \right|}}{{\rm{e}}^{ - Us}}\rho ({\rm{d}}s), $当$\widetilde{W}$是NGGP时${\rm{Pr}}({J'_k} \in {\rm{d}}s|U, X) \propto {s^{\left| k \right| - \sigma - 1}}{{\rm{e}}^{ - (U + \tau )s}}, $对固定原子的位置的更新可以采样Bush和MacEachern提出的加速步骤进行.当W_r更新后, 对S的更新可以通过从其分布进行独立采样得到每个S_i.

4.2 哈密顿蒙特卡罗采样方法 4.2.1 哈密顿动力学方法

MH算法的一个主要的局限是它具有随机游走的行为, 而在状态空间中遍历的距离与步骤数量只是平方根的关系, 仅仅通过增加步长的方式是无法解决这个问题的, 因为这会使得拒绝率变高.HMC采样方法弥补了MH算法的缺陷, 可以更有效地对状态空间进行搜索.HMC起源于对哈密顿动力学(Hamiltonian dynamics)下进行变化的物理系统的行为的模拟, 通过将概率仿真转化为哈密顿系统的形式, 利用哈密顿动力学的框架能够让系统状态发生较大的改变, 同时让拒绝概率较低^[38].

哈密顿动力学对应于在连续时刻T下的状态变量z={z_i}的演化.经典的动力学由牛顿第二定律描述(物体的加速度正比于所受的力), 是关于时间的二阶微分方程.通过引入中间的动量变量r可以将二阶微分方程分解为两个相互耦合的一阶方程, r对应于状态变量z的变化率, 元素为${r_i} = \frac{{{\rm{d}}{z_i}}}{{{\rm{d}}{\cal T}}}{\rm{, }}$从动力学的角度, z_i可以被看做位置变量, 对于每个位置变量, 都存在一个对应的动量变量, 位置和动量组成的联合空间被称为相空间.

不失一般性, 将概率分布Pr(z)写成下面的形式:${\rm{Pr}}(z) = \frac{1}{{{Z_p}}}\exp ( - E(z)), $其中E(z)可以看做状态z处的势能, 系统的总能量是势能和动能之和H(z, r)=E(z)+K(r), 其中H是哈密顿函数, 加速度是动量的变化率, 通过施加力的方式确定, 是势能的负梯度, 即$\frac{{{\rm{d}}{r_i}}}{{{\rm{d}}{\cal T}}} = - \frac{{\partial E(z)}}{{\partial {z_i}}}, $将动能定义为$K(r) = \frac{1}{2}{\left\| r \right\|^2} = \frac{1}{2}\sum\limits_i {r_i^2} , $将系统的动力学用哈密顿方程的形式表示出来, 形式为$\left\{ \begin{array}{l} \;\frac{{{\rm{d}}{z_i}}}{{{\rm{d}}{\cal T}}} = \frac{{\partial H}}{{\partial {r_i}}}\\ \;\frac{{{\rm{d}}{r_i}}}{{{\rm{d}}{\cal T}}} = - \frac{{\partial H}}{{\partial {z_i}}} \end{array} \right..$

哈密顿动态系统的一个重要性质是在动态系统的变化过程中, 哈密顿函数H的值是一个常数, 这一点通过求微分的方式很容易看出来, 即$\frac{{{\rm{d}}H}}{{{\rm{d}}{\cal T}}} = \sum\limits_i {\left\{ {\frac{{\partial H}}{{\partial {z_i}}}\frac{{{\rm{d}}{z_i}}}{{{\rm{d}}{\cal T}}} + \frac{{\partial H}}{{\partial {r_i}}}\frac{{{\rm{d}}{r_i}}}{{{\rm{d}}{\cal T}}}} \right\}} = \sum\limits_i {\left\{ {\frac{{\partial H}}{{\partial {z_i}}}\frac{{\partial H}}{{\partial {r_i}}} - \frac{{\partial H}}{{\partial {r_i}}}\frac{{\partial H}}{{\partial {z_i}}}} \right\}} = 0.$哈密顿动态系统的第2个重要性质是动态系统在相空间中体积不变, 这被称为Liouville定理, 换句话说, 如果我们考虑变量(z, r)空间中的某个区域, 那么当这个区域在哈密顿动态过程下的变化时, 它的形状可能会改变, 但是它的体积不会改变.

考虑相空间上的联合概率分布, 其总能量是哈密顿函数, 概率分布的形式为${\rm{Pr}}(z, r) = \frac{1}{{{Z_H}}}\exp ( - H(z, r)), $利用体系的不变性和H的守恒性, 可以看到哈密顿动态系统会使得Pr(z, r)保持不变.虽然H是不变的, 但是z和r会发生变化因此通过在某个有限的时间间隔上对哈密顿动态系统积分, 就可以让z以某种系统化的方式发生较大的变化, 避免了随机游走的行为.

4.2.2 混合蒙特卡罗算法

HMC也称混合蒙特卡罗(hybrid Monte Carlo), 是一种对连续分布进行近似的MCMC方法:给定独立观测变量集D, 后验分布Pr(z|D)∞exp(–U(z)), 势能函数$U = - \sum\limits_{x \in D} {\log {\rm{Pr}}(x|z)} - \log {\rm{Pr}}(z), $引入一组辅助动量变量r对参数空间进行扩展(r与z具有相同的维度), 先对联合分布$\pi (z, r) \propto \exp \left( { - U(z) - \frac{1}{2}{r^T}{M^{ - 1}}r} \right)$进行采样, 然后丢弃r的样本, 就得了后验分布的样本, M是一个质量矩阵(常被设置成单位矩阵I), 与r一起定义了动能.哈密尔顿函数$H(z, r) = \frac{1}{2}{r^T}{M^{ - 1}}r + U(z)$, 系统的动力学方程式为$\left\{ \begin{array}{l} \;\frac{{{\rm{d}}{z_i}}}{{{\rm{d}}{\cal T}}} = \frac{{\partial H}}{{\partial {r_i}}} = {M^{ - 1}}r{\rm{d}}{\cal T}\\ \;\frac{{{\rm{d}}{r_i}}}{{{\rm{d}}{\cal T}}} = - \frac{{\partial H}}{{\partial {z_i}}} = - \nabla U(z){\rm{d}}t \end{array} \right., $不能直接对连续系统进行采样, 需要先将其离散化.常用的离散化方法是蛙跳积分(leapfrog integrator), 利用下列公式对位置变量和动量变量的离散时间近似$\hat z$和$\hat r$进行交替的更新.

${\hat r_i}\left( {{\cal T} + \frac{\varepsilon }{2}} \right) = {\hat r_i}({\cal T}) - \frac{\varepsilon }{2}\frac{{\partial E}}{{\partial {z_i}}}(\hat z(r)), $

${\hat z_i}({\cal T} + \varepsilon ) = {\hat z_i}({\cal T}) + \varepsilon {\hat r_i}\left( {{\cal T} + \frac{\varepsilon }{2}} \right), $

${\hat r_i}({\cal T} + \varepsilon ) = {\hat r_i}\left( {{\cal T} + \frac{\varepsilon }{2}} \right) - \frac{\varepsilon }{2}\frac{{\partial E}}{{\partial {z_i}}}(\hat z({\cal T} + \varepsilon )).$

蛙跳积分对动量变量的更新形式是半步更新, 步长为ε/2, 接着是对位置变量的整步更新, 步长为ε, 然后是对动量变量的第2个半步更新.连续地使用几次蛙跳, 对动量变量的半步更新就可以结合到步长为ε的整步更新中.于是, 位置变量的更新和动量变量的更新互相之间以蛙跳的形式结合.为了将动态系统归进一个时间间隔${\cal T}$, 我们需要进${\cal T}/\varepsilon $个步骤.因为蛙跳积分中的每一步对z_i或者r_i的更新都只是另一个变量的函数, 即将相空间的一个区域进行形变而不改变它的体积, 所以蛙跳积分精确地保持了相空间的体积不变性.

对于一个非零的步长ε, 蛙跳算法的离散化会在哈密顿动力学方程的积分过程中引入误差, 为了消除与离散化过程关联的数值误差, HMC将Hamiltonian动态系统框架与Metropolis算法结合在一起, 构造了一个马尔可夫链, 它由对动量r的随机更新和使用蛙跳算法对哈密顿动力系统的更新交替组成.每次应用蛙跳算法之后, 基于哈密顿函数H的值, 确定Metropolis准则, 确定生成的候选状态被接受或者拒绝.HMC完美地模拟了哈密顿动态系统, 因此每个这种候选状态都会高概率地被接受(因为H的值会保持不变)^[39].HMC采样过程描述见算法2.

算法2. Hamiltonian Monte Carlo采样.

输入:起始位置z⁽¹⁾, 步长ε.

输出:后验分布Pr(z|D)的有限采样N.

算法步骤:

for t=1, 2, ...do

1.  对动量r进行重新采样:${r^{({\cal T})}} \sim Norm(0, M), \;({z_0}, {r_0}) = ({z^{({\cal T})}}, {r^{({\cal T})}})$

2.  对离散化的对哈密顿动力系统进行采样

    a)   ${z_0}: = {z_0} - \frac{\varepsilon }{2}\nabla U({z_0})$

    b)   for i=1 to m do

     $\begin{array}{l} {z_i}: = {z_{i - 1}} + \varepsilon {M^{ - 1}}{r_{i - 1}}\\ {r_i}: = {r_{i - 1}} - \varepsilon \nabla U({z_i}) \end{array}$

    end

    c)  ${r_m}:\; = {r_m} - \frac{\varepsilon }{2}\nabla U({z_m})$

    d)  $(\hat z, r): = ({z_m}, {r_m})$

  3.进行Metroplis-Hastings修正, 消除离散化带来的误差

    a)  计算接受率$\rho : = {{\rm{e}}^{H(\hat z, r) - H({z^{({\cal T})}}, {r^{({\cal T})}})}}$

    b)  从[0, 1]上的均匀分布采样一个μ, 若min(1, ρ) > μ, 则${z^{(t + 1)}}: = \hat z$

    End

与基本的MH方法不同, HMC能够利用到对数概率分布的梯度信息和概率分布本身的信息^[40].在函数最优化领域有一个类似的情形, 大多数可以得到梯度信息的情况下, 使用哈密顿动力学方法是很有优势的.可以直观地解释为, 这种现象是由于下列的事实造成的:在n维空间中, 与计算函数本身的代价相比, 计算梯度所带来的额外计算代价通常是一个与n无关的固定因子.而与函数本身只能传递一条信息相比, n维梯度向量可以传递n条信息.

4.3 稀疏可交换图模型的后验分布的推理

首先对受限CRM W_φ的条件分布W_φ|D_φ进行分析, W_φ|D_φ可以被分解成两部分之和, 一部分是全体带随机权值w_i的固定位置θ_i的测度, 对应着观测到的网络数据中至少关联了1条边的节点全体, 另一部分是权值随机且位置随机的测度, 对应了其他节点, 用w^*表示这些节点的总的随机质量.

定理4.1. 令$({\theta _1}, ..., {\theta _{{N_\alpha }}})$是有向图D_α的支撑点, ${D_\alpha } = \sum\nolimits_{1 \le i, j \le {N_\alpha }} {{n_{ij}}\delta ({\theta _i}, {\theta _j})} $, 令${m_i} = \sum\nolimits_{j = 1}^{{N_\alpha }} {({n_{ij}} + {n_{ji}}) > 0} $, i=1, ..., N_α, 则有${W_\alpha }{\rm{|}}{D_\alpha }$${w_*}\sum\limits_{i = 1}^\infty {{{\widetilde{P}}_i}{\delta _{{{\mathit \theta }_i}}} + \sum\limits_{i = 1}^{{N_\alpha }} {{w_i}{\delta _{{\theta _i}}}} } {w_*}\sum\limits_{i = 1}^\infty {{{\widetilde{P}}_i}{\delta _{{{\mathit \theta }_i}}} + \sum\limits_{i = 1}^{{N_\alpha }} {{w_i}{\delta _{{\theta _i}}}} } {\rm{, }}$这里, ${\theta _i} \sim Unif([0, \alpha ]){\rm{, }}$权值序列${({\widetilde{P}_i})_{i = 1, 2, ...}}$满足${\widetilde{P}_1} > {\widetilde{P}_2} > ...$且$\sum\nolimits_{i = 1}^\infty {{{\widetilde{P}}_i}} = 1{\rm{, }}$给定${w_*}$, $({\widetilde{P}_i})$服从一个列维强度为ρ的Poisson-Kingman分布, 即$({\widetilde{P}_i})|{w_*} \sim PK(\rho |{w_*}).$最终, 给定D_α, 权值$(w, {w_*})$的联合依赖条件分布为$({w_1}, ..., {w_{{N_\varphi }}}, {w_*})|{D_\alpha } \propto \left[ {\coprod\limits_{i = 1}^{{N_\alpha }} {w_i^{{m_i}}} } \right]{{\rm{e}}^{ - {{\left( {\sum\nolimits_{i = 1}^{{N_\alpha }} {{w_i} + {w_*}} } \right)}^2}}}\left[ {\coprod\limits_{i = 1}^{{N_\alpha }} {\rho ({w_i})} } \right] \times g_\alpha ^*({w_*}), $这里, $g_\alpha ^*$是随机变量$W_\alpha ^* = {W_\alpha }(\left[ {0, \alpha } \right])$的概率密度函数, 其拉普拉斯变换为$\mathit{\mathbb{E}}\left[ {{\text{e}}^{-tW_{\alpha }^{*}}} \right]={{\text{e}}^{-\alpha \psi (t)}}$^[15].

需要注意的是, 规一化的权值${{({{\widetilde{P}}_{i}})}_{i=1, 2, ...}}$和位置${{({{\theta }_{i}})}_{i=1, 2, ...}}$不是似然可识别的(likelihood identifiable), 因为观测到的网络节点数据仅包含节点权值和${{w}_{*}}$的信息, 还需要注意的是$({{w}_{1}}, ..., {{w}_{{{N}_{\alpha }}}}, {{w}_{*}})\text{ }\!\!|\!\!\text{ }{{D}_{\alpha }}$不依赖于节点位置$({{\theta }_{1}}, ..., {{\theta }_{{{N}_{\alpha }}}}), $我们考虑的是齐次CRM, 这一点非常重要, 因为节点位置$c$通常无法观测到的, 推理过程不会考虑节点位置.

James等人^[41]通过引入一个辅助随机变量U刻画了NRMI所导出的EFFP, 令Γ_n是一个Gamma随机变量, 其位置参数是1, 形状参数是n, 并且独立于总质量T, 取值为正的随机变量U_n=Γ_n/T, 容易证明对于任意的n≥1可以给出的U_n密度函数${{f}_{{{U}_{n}}}}(W)=\frac{{{W}^{n-1}}}{\Gamma (n)}\int_{{{\bf R}_{+}}}{{{t}^{n}}{{\text{e}}^{-Wt}}{{f}_{T}}(t)\text{d}t}.$

特别地, $\text{Pr}(\pi =\{{{A}_{1}}, ..., {{A}_{k}}\}, \{X_{k}^{*}\in \text{d}{{x}_{k}}:k\in \pi \}, U\in \text{d}W|W)=\frac{1}{\Gamma (n)}{{W}^{n-1}}{{\text{e}}^{-TW}}\text{d}W\prod\limits_{k\in \pi }{W{{(\text{d}{{x}_{k}})}^{\left| k \right|}}}$(给定W, X和U的联合条件分布), 利用Palm公式, 可以由(w, w_*)的联合依赖条件分布推导出对EPPF和预测分布体系的刻画.

定理4.2. 令$\widetilde{W}$是一个齐次NRM, 其列维测度为ρ, 基分布为μ₀, 拉普拉斯指数为ψ(W), 边缘化$\widetilde{W}$后, X和U的联合分布为$\text{Pr}(\pi =\{{{A}_{1}}, ..., {{A}_{k}}\}, \{X_{k}^{*}\in \text{d}{{x}_{k}}:k\in \pi \}, U\in \text{d}W)=\frac{1}{\Gamma (n)}{{W}^{n-1}}{{\text{e}}^{-\psi (W)}}\text{d}W\prod\limits_{k\in \pi }{{{\mathcal{C}}_{\left| k \right|}}(W){{\mu }_{0}}(\text{d}{{x}_{k}})}, $${{\mathcal{C}}_{m}}(W)$是指数倾斜列维测度${{\text{e}}^{-us}}\rho (\text{d}s)$的第m阶矩${{\mathcal{C}}_{m}}\left( W \right)=\int_{0}^{\infty }{{{w}^{m}}{{\text{e}}^{-wt}}\rho (w)\text{d}w}.$特定地, 当将辅助随机变量U边缘化后, EPPF的表达式为$\text{Pr}(\pi =\{{{A}_{1}}, ..., {{A}_{k}}\})=\int_{{{R}_{+}}}{\frac{1}{\Gamma (n)}{{u}^{n-1}}{{\text{e}}^{-\psi (W)}}\text{d}W\prod\limits_{k\in \pi }{{{\mathcal{C}}_{\left| k \right|}}(W)\text{d}W}}, $因为序列$\{X_{k}^{*}\in \text{d}{{x}_{k}}:k\in \pi \}$独立同分布于μ₀, 相应的预测分布体系为$\text{Pr}({{X}_{n+1}}\in \text{d}x|U, X)\propto {{\mathcal{C}}_{\left| k \right|}}(W){{\mu }_{0}}(\text{d}x)+\sum\limits_{k\in \pi }{\frac{{{\mathcal{C}}_{\left| k \right|+1}}(U)}{{{\mathcal{C}}_{\left| k \right|}}(U)}}{{\delta }_{X_{k}^{*}(\text{d}x)}}\text{, }$注意, U在给出X时的后验分布为$\text{Pr}(U\in \text{d}W|X)\propto {{W}^{n-1}}{{\text{e}}^{-\psi (W)}}\text{d}W\prod\limits_{k\in \pi }{{{\mathcal{C}}_{\left| k \right|}}(W)}$^[15].

定理4.2表明, 在边缘采样中, CRM W被边缘后, 就可以采样得到代表X的分布的分区结构和聚类参数.

定理4.3. 令$\widetilde{W}$是一个齐次NRMI, 列维测度为ρ, 基分布为μ₀, 对应的CRM W在给出X和U时的后验分布为$W|U, X\sim {W}'+\sum\limits_{k\in \pi }{{{{{J}'}}_{k}}{{\delta }_{X_{k}^{*}}}}\text{, }$${W}'$是一个齐次CRM, 具有指数倾斜列维强度${\upsilon }'(\text{d}s, \text{d}y)={{\text{e}}^{-Us}}\rho (\text{d}s){{\mu }_{0}}(\text{d}y)\text{, }$随机质量$\{{{{J}'}_{k}}:k\in \pi \}$相互独立并且与${W}'$独立, 其条件分布为$\text{Pr}({{{J}'}_{k}}\in \text{d}s|U, X)=\frac{1}{{{\mathcal{C}}_{\left| k \right|}}(U)}{{s}^{\left| k \right|}}{{\text{e}}^{-Us}}\rho (\text{d}s)\text{, }$$\widetilde{W}$在给出X和U时的后验分布就是对W, U, X进行规一化^[15].

定理4.3表明一个齐次NRMI $\widetilde{W}$在给出X和U时的后验分布仍然是一个NRMI.

以下介绍Caron和Fox提出的对生成模型为(2.4)的稀疏可交换图所进行的后验推理^[15].针对CRM是一个GGP的情况, 其后验推理的MCMC采样过程, 令$\phi =(\alpha , \sigma , \tau )$是超参数集合, 这些超参数也需要推理得到, 为这些超参数假设合适的先验$\text{Pr}(\alpha )\propto \frac{1}{\alpha }, \ \text{Pr}(\sigma )\propto \frac{1}{1-\sigma }, \ \text{Pr}(\tau )\propto \frac{1}{\tau }$, 为了强调列维测度和w_*的概率密度函数对超参数的依赖, 记列维强度为ρ(w|σ, τ)和$g_{\alpha , \sigma , \tau }^{*}({{w}_{*}})\text{.}$对于有向多图, 需要近似推理的是后验分布$\Pr ({{w}_{1}}, ..., {{w}_{{{N}_{\varphi }}}}, {{w}_{*}}, \phi |$${{({{n}_{ij}})}_{1\le i, j\le {{N}_{\varphi }}}}), $对于无向图, 需要要近似推理的是后验分布$\text{Pr}({{w}_{1}}, ..., {{w}_{{{N}_{\varphi }}}}, {{w}_{*}}, \phi |{{({{z}_{ij}})}_{1\le i, j\le {{N}_{\varphi }}}}).$

在观测数据是一个无向简单图的情况下, 需要推算未知的有向边数n_ij, 观测到z_ij=1(i≤j)时, 引入一个潜变量${{\bar{n}}_{ij}}={{n}_{ij}}+{{n}_{ji}}$, 其条件分布为

${{\bar{n}}_{ij}}|z, w\sim \left\{ \begin{array}{*{35}{l}} {{\delta }_{0}}\text{, } & \text{if}\ {{z}_{ij}}=0 \\ tPoisson(2{{w}_{i}}{{w}_{j}})\text{, } & \text{if}\ {{z}_{ij}}=0, i\ne j \\ tPoisson({{w}_{i}}^{2})\text{, } & \text{if}\ {{z}_{ij}}=0, i=j \\ \end{array} \right.$

(4.1)

此处, tPoisson(λ)是一个零截断泊松分布, 其概率分布密度函数为$\frac{{{k^\lambda }\exp ( - \lambda )}}{{(1 - \exp ( - \lambda ))k!}}\, , \, {\rm{for}}\, k = 1, 2, ...$, 方便起见, 令${\bar n_{ij}} = {\bar n_{ji}}, $${m_i} = \sum\nolimits_{j = 1}^{{N_\varphi }} {{{\bar n}_{ij}}} .$为了更有效地得到所要推理的后验分布, 使用HMC与Gibbs采样相结合对$({w_1}, {w_2}, ..., {w_{{N_\varphi }}})$进行更新, HMC需要计算对数后验的梯度${\left[ {{\nabla _{{\omega _{1:{N_\alpha }}}}}\log (\Pr ({w_{1:{N_\alpha }}}, {w_*}|{D_\alpha })} \right]_i} = {m_i} - \sigma - {w_i}\left( {\tau + 2\sum\limits_{j = 1}^{{N_\alpha }} {{w_j} + 2{w_*}} } \right), $${\omega _i} = \log {w_i}$.采用Metropolis-Hastings算法对总质量w_*和超参数ϕ进行更新.需要注意的是, 除了特定的当$\sigma = 0, \frac{1}{2}$时, 其他情况下w_*的概率密度函数$g_{\alpha , \sigma , \tau }^*({w_*})$都不能被解析表示, 需要通过对$g_{\alpha , \sigma , \tau }^*({w_*})$的指数倾斜进行采样得到w_*^[15].

综上所述, 当稀疏可交换图是一个有向图时, 后验分布${\rm{Pr}}({w_1}, ..., {w_{{N_\alpha }}}, {w_*}, \phi |{({n_{ij}})_{1 \le i, j \le {N_\alpha }}})$的MCMC推理过程包括以下两个步骤, 若是一个无向图是, 后验分布${\rm{Pr}}({w_1}, ..., {w_{{N_\alpha }}}, {w_*}, \phi |{({z_{ij}})_{1 \le i, j \le {N_\alpha }}})$的MCMC推理还需要加上第3个步骤.

1.给定${w_*}, \phi , ({n_{ij}})$, 采用HMC方法对社交能力参数$w = ({w_1}, {w_2}, ..., {w_{{N_\alpha }}})$进行更新;

2.给定$({w_1}, {w_2}, ..., {w_{{N_\alpha }}})$和(n_ij), 采用MH算法对总质量w_*和超参数ϕ进行更新;

3.给定$({w_1}, {w_2}, ..., {w_{{N_\alpha }}})$, w_*, ϕ, 采用条件分布或MH算法对${\bar n_{ij}}$进行更新.

步骤1:采用HMC对$({w_1}, {w_2}, ..., {w_{{N_\alpha }}})$进行更新, 令L≥1是蛙跳步数, ε > 0是步长, 对数后验的梯度记作$U'(w, {w_*}, \phi ) = {\nabla _{{\omega _{1:{N_\alpha }}}}}\log (Pr({w_{1:{N_\alpha }}}, {w_*}|{D_\alpha }){\rm{, }}$HMC算法过程如下.

(1) 首先对辅助动量变量$r:N(0, {I_{{N_\alpha }}})$进行采样.

(2) (w, r)为哈密顿动力系统的初始状态, 通过蛙跳离散化得到候选状态$(\widetilde{w}, \widetilde{r})$, 步骤如下.

① ${\widetilde{r}^{(0)}} = r + \frac{\varepsilon }{2}U'(w, {w_*}, \phi )$${\widetilde{w}^{(0)}} = w$;

② ${\rm{for}}\;l = 1, ..., L - 1$

    $\begin{array}{l} \log {{\widetilde{w}}^{(l)}} = \log {{\widetilde{w}}^{(l - 1)}} + \varepsilon {{\widetilde{r}}^{(l - 1)}}, \\ {{\widetilde{r}}^{(l)}} = \;{{\widetilde{r}}^{(l)}} + \varepsilon U'({{\widetilde{w}}^{(l)}}, {w_*}, \phi ); \end{array}$

③ $\log \widetilde{w} = \log {\widetilde{w}^{(L - 1)}} + \varepsilon {\widetilde{r}^{(L - 1)}}, $

    $\begin{array}{l} \widetilde{r} = \; - \left[ {{{\widetilde{r}}^{(L - 1)}} + \frac{\varepsilon }{2}U'(\widetilde{w}, {w_*}, \phi )} \right], \\ \widetilde{w} = {{\widetilde{w}}^{(L)}}. \end{array} $

(3) 以概率accep=min(1, ρ)接受$(\widetilde{w}, \widetilde{r})$, ρ为接受率, 以下省略关于ρ的计算表达式, 可参阅文献[15].

步骤2:采用MH算法对${w_*}, \alpha , \sigma , \tau $进行更新, 从提议分布$q(\widetilde \alpha , \widetilde \sigma , \widetilde \tau , {\widetilde{w}_*}|\alpha , \sigma , \tau , {w_*})$采样得到$(\widetilde \alpha , \widetilde \sigma , \widetilde \tau , {\widetilde{w}_*})$以概率accep=min(1, mhr)接受它.可以对提议分布进行因子分解, 得到:

$ q(\widetilde \alpha , \widetilde \sigma , \widetilde \tau , {\widetilde{w}_*}|\alpha , \sigma , \tau , {w_*}) = q(\widetilde \tau |\tau )q(\widetilde \sigma |\sigma )q(\widetilde \alpha |\widetilde \sigma , \widetilde \tau , {w_*})q({\widetilde{w}_*}|\widetilde \alpha , \widetilde \sigma , \widetilde \tau , {w_*}). $

这里,

$ \begin{array}{c} q(\widetilde \tau |\tau ) = \log normal(\widetilde \tau ;\log (\tau ), \sigma _\tau ^2),\\ q(\widetilde \sigma |\sigma ) = \log normal(1 - \widetilde \sigma ;\log (1 - \sigma ), \sigma _\tau ^2), \\ q(\widetilde \alpha |\widetilde \sigma , \widetilde \tau , {w_*}) = Gamma\left( {\widetilde \alpha ;{N_\varphi }, \frac{{{{\left( {\tau + 2\sum {{w_i} + {w_*}} } \right)}^{\widetilde \sigma }} - {\tau ^{\widetilde \sigma }}}}{{\widetilde \sigma }}} \right), \\ q({{\widetilde{w}}_*}|\widetilde \alpha , \widetilde \sigma , \widetilde \tau , {w_*}) = g_{\widetilde \alpha , \widetilde \sigma , \widetilde \tau + 2\sum {{w_i}} + {w_*}}^*({{\widetilde{w}}_*}) = \frac{{\exp \left( { - 2\sum {{w_i}} - {w_*}} \right)g_{\widetilde \alpha , \widetilde \sigma , \widetilde \tau }^*({{\widetilde{w}}_*})}}{{\exp ( - {\psi _{\widetilde \alpha , \widetilde \sigma , \widetilde \tau }}({{\widetilde{w}}_*}))}}. \end{array} $

${\widetilde{w}_*}$的提议分布的选择基于${\widetilde{w}_*}$的分布$g_{\widetilde \alpha , \widetilde \sigma , \widetilde \tau }^*({\widetilde{w}_*})$的指数倾斜, 这样就可以简化计算接受率ρ的表达式.指数倾斜也称埃舍尔变换, 是对随机测度施加的一个指数式改变.若函数$\varphi :[0, \infty ] \to {\bf{R}}$连续且各阶导数存在并满足(–1)^kφ^(k)(t)≥0(对于任意t∈(0, ∞), k∈N), 则称φ完全单调.令Ψ_∞表示满足φ(0)=1的全体完全单调函数φ:[0, ∞]→[0, 1]构成的集合, 由Bernstein定理:[0, ∞]上的分布函数F的Laplace-Stieltjes变换是一个完全单调函数φ:[0, ∞]→[0, 1]当且仅当φ∈Ψ_∞, 即$\varphi = {\cal L}{\cal S}(F)$或$F = {\cal L}{{\cal S}^{ - 1}}(\varphi )\quad {\rm{iff}}\;\varphi \in {\Psi _\infty }$, 容易证明$\widetilde \varphi (t): = \varphi (\lambda {\rm{ + }}t)/\varphi (t)$也在Ψ_∞中, 对应的分布函数$\widetilde F = {\cal L}{{\cal S}^{ - 1}}(\widetilde \varphi )$称作是分布函数$F = {\cal L}{{\cal S}^{ - 1}}(\varphi )$的指数倾斜(λ称作倾斜参数), 因为若F的密度函数为f, 则$\widetilde F$的密度函数$\widetilde f(x) = \exp ( - \lambda x)f(x)/\varphi (\lambda ), x \in [0, \infty )$是对f做了指数倾斜^[42].

步骤3:对潜变量${\bar n_{ij}}$的更新, 可以通过直接对条件分布(4.1)进行采样得到.当网络中的边数量众多的时候, 一个更有效的更新方法是使用一个Metropolis-Hastings提议$q({{\tilde{\bar{n}}}_{ij}}|{{\bar{n}}_{ij}})=\left\{ \begin{array}{*{35}{l}} \frac{1}{2}\text{,} & 若{{{\tilde{n}}}_{ij}}={{n}_{ij}}+1\ ,\ {{n}_{ij}}>1 \\ \frac{1}{2}\text{,} & 若{{{\tilde{n}}}_{ij}}={{n}_{ij}}-1\ ,\ {{n}_{ij}}>1 \\ 1\text{,} & 若{{{\tilde{n}}}_{ij}}={{n}_{ij}}+1\ ,\ {{n}_{ij}}=1 \\ 0\text{,} & 其他 \\ \end{array} \right.,$并且以概率$\min \left( 1,\frac{{{{\bar{n}}}_{ij}}!}{{{{\tilde{\bar{n}}}}_{ij}}!}{{((1+{{\delta }_{ij}}){{w}_{i}}{{w}_{j}})}^{{{{\tilde{\tilde{n}}}}_{ij}}-{{{\tilde{n}}}_{ij}}}}\frac{q({{{\bar{n}}}_{ij}}|{{{\tilde{\bar{n}}}}_{ij}})}{q({{{\tilde{\bar{n}}}}_{ij}}|{{{\bar{n}}}_{ij}})} \right)$接受从提议分布中得到的采样.

4.4 对指数倾斜稳定分布进行采样(sampling exponentially tilted stable distribution)

在σ≥0的情况下, $W_\varphi ^*$是一个服从指数倾斜稳定分布(ETSD)的随机变量.稳定分布S(α, β, γ, d; 1)又称为非高斯稳定分布、重尾分布, 由列维提出, 是唯一满足广义中心极限定理的分布(即无限多个可能方差无限大的独立分布的随机变量之和, 其极限分布是稳定分布)^[43].高斯分布(α=2)、柯西分布(α=1, β=0)、列维分布(α=1/2, β=1)都仅是稳定分布的特例, 很多不满足经典的中心极限定理的数据都可以用稳定分布来描述, 因此具有更普遍的应用范围.同高斯分布具有指数衰减的拖尾不同, 稳定分布的拖尾以平方律衰减, 衰减的速度与指数α有关, α越小, 分布的衰减就越缓慢, 分布的拖尾越重.

指数倾斜稳定分布是严格稳定分布天然的指数族分布, 采用指数倾斜可以在求解贝叶斯估计的近似公式时不涉及似然函数的条件最大值, 求解过程更稳定, 并且显著减少了所需要的计算时间.将严格稳定分布与指数倾斜相结合已经被证明在很多应用(例如, 重要性采样、罕见事件的模拟、保险精算等)中非常有效^[44].

Brix、Rosinski、Ridout、Devroye、Hofert对指数倾斜稳定分布采样算法进行了研究, 针对分布为$\widetilde S \sim \widetilde S(\alpha , 1, $$\cos {(\alpha {\rm{ \mathsf{ π} }}/2)^{1/a}}, {1_{\{ \alpha = 1\} }}, \lambda {1_{\{ \alpha \ne 1\} }};1), $$\alpha \in (0, 1], \lambda \in [0, \infty )$的指数倾斜稳定分布, Hofert^[44]提出了一种快速拒绝采样算法, Devroye^[45]提出了一种双拒绝采样算法, 这两种算法都属于精确采样算法.快速拒绝采样算法简单容易实现, 当参数范围很大时采样速度也很快; 双拒绝采样算法比较复杂, 但是在所有参数组合上的复杂性界是均等的, 非常适用于那些即使在线性复杂度λ^α下也很耗时的参数的采样问题.Hofert^[42]修正了双拒绝采样算法中的两个小问题, 并且解决了在α较小时存在的一些临界数值计算过程, 使得双拒绝采样算法成为一种非常重要的指数倾斜稳定分布采样算法.

5 实验验证

Caron等人^[15]采用他们提出的生成模型模拟生成了一系列随机图(s和τ取不同值情况下), 并与采用ER模型^[7]、PA模型^[13]、Lloyd^[18]模拟生成的随机图进行了对比, 结果表明Caron-Fox模型能够模拟生成节点度分布呈现幂律分布而且具有重尾特性的稀疏图, 并且能够处理节点度分布尾部的指数截断(这一性质是进行网络数据分析时非常有用的一个重要性质), 但是采用3个作对比的生成模型模拟生成的随机图则不能表现出这些性质.

如图 3(a)所示, 在log-log坐标系下, 采用GGP作为先验的Caron-Fox模型在σ分别取0.2, 0.5, 0.8的情况下模拟生成的随机图, 其节点度分布呈现出明显的幂律分布, 而p=0.05的ER模型G(n, p)和Lloyd模型模拟生成的随机图, 其节点度分布则不能呈现出幂律分布, PA模型模拟生成的随机图的节点度分布虽然呈现出了幂律分布, 但不是可交换图, 图 3(b)所示的是当τ分别取0.1, 1, 5时的情况, 可以得到和图 3(a)类似的结论(图例中的BA表示PA模型, 因为PA模型是由Barabási和Albert共同提出的, 所以也常常被称作BA模型).

如图 4所示, 采用GGP作为先验的Caron-Fox模型模拟生成的随机图表现出了这样一种特性:随着节点个数的增加, 图中度是1的节点的个数是线性增加, PA模型模拟生成的随机图有类似的性质, 但是ER模型模拟生成的随机图中度是1的节点的个数是随着节点个数的增加呈现指数减少的, Lloyd模型模拟生成的随机图则几乎没有度是1的节点.

如图 5所示, 采用GGP作为先验的Caron-Fox模型模拟生成的随机图表现出了这样一种特性:随着节点个数n的增加, 图中边的个数呈现o(n²)增加(稀疏图的特征), PA模型模拟生成的随机图有类似的性质, 但是ER模型和Lloyd模型模拟生成的随机图中, 边的个数呈现Θ(n²)增加(稠密图的特征).

实证结果表明, 当σ≥0时, 采用GGP作为先验的Caron-Fox模型可以生成稀疏可交换图, 而基于Aldous- Hoover表示理论的ER模型和Lloyd模型只能生成稠密可交换图, PA模型虽然模拟生成稀疏图, 但不是可交换的.关于实证分析的更详细过程, 可参阅文献[15].

6 最新进展与研究展望

在随机图模型方面, 通过对经典Graphon的定义进行推广, Borgs等人^[46]提出了Graphon process模型, Veitch等人^[11]提出了Graphex模型, 这两种更广义的稀疏可交换图模型都可以对一大类可交换图进行建模, 并且将Caron-Fox模型和传统可交换稠密图作为特例纳入了一个统一的模型框架.Crane等人^[47]认为, 尽管Caron-Fox模型能够生成满足可交换性的稀疏网络, 但是尚不清楚可交换性是否在各种采样策略下都能得到保证, 他们提出了一种边可交换的随机图模型(edge exchangeable model); Jacobs和Clauset^[48]对各种网络生成模型从哲学思想、表示和统一性这3个方面进行了总结, 提出了网络生成模型在可解释性、同质性、优化、模型可实现、模型选择、模型检验等方面所面临的挑战和机遇, 并且指出Caron-Fox模型只能生成超线性密度(superlinear density)O(n^α), 1 < α < 2的稀疏可交换图, 还不能生成更稀疏的O(n)数量级的节点数对应O(n)数量级的边数的稀疏可交换图, 还需要在理论创新和推理算法方面进行大量的研究才能使得这种基于连续空间的网络生成模型被应用到更广泛的真实网络数据分析中.

社区结构作为网络模型的重要特征刻画了网络的多种属性和功能, 对人们准确理解复杂系统的特性有十分重要的意义.从社区结构这一中观层面对复杂网络进行研究对人们准确理解复杂系统的特性有十分重要的意义, 既能够弥补宏观层面粒度过粗所造成的很多网络特性无法观察到的缺陷, 又避免了微观层面粒度过细所带来的丢失共性并且计算复杂度高等问题.社区发现作为网络数据分析的重要内容已经得到了学术界的广泛关注和深入研究, Herlau等人^[14]将SBM与Caron-Fox模型相结合, 为稀疏可交换图引入了社区结构, 并提出了稀疏可交换图的非重叠静态社区发现方法; Todeschini和Caron^[49]提出了稀疏可交换图的重叠静态社区发现方法, 该方法受到了Zhou^[50]提出的无限边分区模型(Infinite edge partition model)的启发, 假设网络中的节点i关联了一组权重参数w_ik, k=1, ..., p(可以解释为节点的不同方面的社交能力), 节点i关联的权重参数是有依赖关系的混合随机测度(compound random measure)^[51], 网络的生成模型为${z_{ij}} \sim Bern\left( {1 - \exp \left( { - \sum\nolimits_{k = 1}^p {\sum\nolimits_{l = 1}^p {{\eta _{kl}}{w_{ik}}{w_{jl}}} } } \right)} \right), $这样就可以同时发现网络中的同配结构和异配结构.Borgs等人^[46]提出了将MMSB与Graphon process模型相结合的重叠静态社区发现方法.

在稀疏可交换图的动态演化方面, Palla等人^[52]提出了采用动态泊松点过程对节点社交能力随时间演化的特性进行建模, 可以同时捕获节点间边的平滑演化(smooth evolution)和长程演化(long term evolution).Matias^[53]指出采用离散模型对网络进行建模, 在研究网络和社区结构的动态演化时, 需要将数据聚集到预先指导的时间段以得到网络快照序列, 数据聚集会带来信息损失, 并且时间段的选择会直接影响预测结果, 但是采用随机过程作为模型参数的连续模型, 则能够很自然地将时间信息引入网络模型, 避免了上述问题.

在近似推理方法方面, 现有的关于稀疏可交换图的研究工作几乎都采用的是MCMC推理方法.相对于随机变分推理(stochastic variational inference, 简称SVI)方法, MCMC方法因为可扩展性和可并行性差从而导致其很难应用于大数据环境下的数据分析, 因此如何采用SVI对稀疏可交换图进行推理是一个非常重要的研究课题. Roychowdhury^[54]提出了伽马过程(GP)的变分推理方法; Tank等人^[55]提出了贝叶斯非参数混合模型的流式变分推理算法; 缘于MCMC与VI在能量函数和信息熵层面颇有渊源, Salimans等人^[56]提出通过桥接二者之间的嫌隙将MCMC和VI结合起来进行推理; Wolf等人^[57]提出了一种将HMC和VI相结合的算法, 把HMC步骤包含到变分下界, 以加速后验推理过程的收敛.

Caron-Fox模型中提出的“节点社交能力”这一概念可以为网络分析提供非常合理的社会学解释——节点间建立边的概率取决于两个节点潜在的社交能力, 并且正是因为具有不同方面的社交能力(或者说不同方面的潜在兴趣), 决定了一个节点可以同时隶属于不同的社区.“节点社交能力”与度纠正的随机块模型(degree- corrected SBM)^[58]中的“度偏好”异曲同工, 都能使网络生成模型很好地拟合具有各种特定度分布的真实网络, 使得Caron-Fox模型与SBM模型相结合时, 不需要再进行度纠正.

稀疏可交换图模型是对传统可交换稠密图模型的扩展.关于传统可交换稠密图的研究, 从网络建模、网络演化到社区发现、社区演化, 再到链路预测、影响力分析、网络传播等各个方面的研究工作, 都已经历时弥久并且成果丰硕.尽管传统可交换稠密图的先天不足导致其无法再继续成为对真实复杂网络进行数据分析的利器, 但是很多已经被提出的研究问题和研究成果可以在进行稀疏可交换图研究时被广泛继承和发展.

综上所述, Caron等人的开创性工作为复杂网络数据分析带来了新的转机, 稀疏可交换图模型不仅解决了可交换性与稀疏性二者不可兼顾的矛盾, 并且具有很多独特的优良基因.随着对稀疏可交换图模型以及贝叶斯非参数混合模型的深入研究, 必将可以推动复杂网络数据分析, 尤其是网络社区结构动态演化研究工作的进一步发展.

致谢 文中关于哈密顿动力学方法和混合蒙特卡罗算法的部分内容引自马毅未公开发表的对文献[39]的翻译手稿, 在此向马毅表示感谢.