图（graph）是计算机科学中一类重要的数据结构，由节点和连接这些节点的边构成. 在图上进行随机游走（random walks on graphs）是一种常用的图信息挖掘方式，所述随机游走本质上是一条马尔可夫链（Markov chain），随机游走在各步走到的节点对应马尔可夫链在各时刻的状态.

具体来说，给定随机游走的源节点$s$，假设该游走在当前时刻走到节点$u$（在起始时刻$u$即为$s$），在下一时刻，该游走会从节点$u$的出邻居集合中随机选取一个节点$v$，从节点$u$走到节点$v$. 随机游走选取节点$v$的概率与边$(u,v)$的权重有关. 在无权图上，图上所有边的权重为1. 因此，随机游走会从节点$u$的出邻居集合中均匀随机地选取一个节点作为节点$v$. 若在有权图上，图上各边带有边权，随机游走则会按照节点$u$出边的边权分布选择节点$v$. 具体而言，如果用${d}_{u}$表示节点$u$的出度，用${{ A}}_{uv}$表示边$(u,v)$的权重，则若一条随机游走在当前时刻走到节点$u$，其在下一时刻会以$\dfrac{{{ A}}_{uv}}{{d}_{u}}$的概率从节点$u$的所有出邻居中选出节点$v$，从节点$u$走到节点$v$.

计算图上的一条随机游走在其游走各步走到各节点的概率被称为随机游走概率的计算（random-walk probability computation）问题，在可控误差下进行随机游走概率的近似计算是图论、数据挖掘、网络分析等领域的核心研究问题，相关技术在图节点邻近度计算^[1-7]、图神经网络^[8-10]、社区发现^[11-12]等领域均有广泛应用.

本文关注动态有权图上随机游走概率的近似计算问题. 此处，动态有权图指图上各边带有权重、图结构会随时间动态变化的图数据类型. 动态有权图在实际应用场景中极为常见. 例如，在社交网络中，社交用户与图节点对应，用户之间的聊天关系与图节点间的连边相对应，用户之间联系的频繁程度可被量化为边权，用户的注册与注销行为可由图结构上节点数的增减反映. 因此，真实应用场景中的社交网络可由动态有权图来刻画. 类似地，道路交通网络也可用动态有权图刻画，其中站点对应图节点，站点之间的通达关系对应节点间的连边，站点间的距离反映为边权，站点和道路的增减反映为图结构的动态变化.

尽管动态有权图可以更好地建模真实网络结构，但由于其上随机游走概率计算的复杂性，现有方法大多面向静态无权图，而在动态有权图上的执行效率不高. 具体而言，现有计算随机游走概率的方法大多是蒙特卡洛方法（Monte-Carlo method），即通过在图上产生大量的随机游走来完成图上随机游走概率的近似估计. 因此，随机游走概率的计算效率和各条随机游走的生成效率直接相关，而生成一条随机游走的时间取决于在游走各步从当前节点的出邻居集合中采样节点的时间. 在静态无权图上，借助简单的均匀采样（uniform sampling）技术即可在$O\left(1\right)$时间内完成1步随机游走采样. 但在动态有权图上，随机游走采样需借助特定的采样技术方能高效完成. 现有方法普遍调用权重采样（weighted sampling）技术进行随机游走采样，但受限于权重采样的采样模式设定，其在动态有权图上的采样效率不高，且普遍依赖于复杂的采样结构. 这使得现有计算随机游走概率的方法在动态有权图上或计算效率低下，或受制于复杂的采样结构而难以实时处理动态图的结构变化. 如何在动态有权图上高效完成随机游走概率的近似计算，虽经多年研究但仍未被解决.

针对上述问题，本文提出将硬币翻转采样（coin flip sampling）技术应用至动态有权图上的随机游走概率问题中. 本文证明，相较于既往借助权重采样技术在动态有权图上生成随机游走的方法，基于硬币翻转采样的方法更具灵活性，无需复杂的采样结构即可在动态有权图上同时实现随机游走概率的高效计算和实时更新.

具体而言，本文首先为现有硬币翻转采样方法设计了一种支持动态图上实时更新的采样结构，该采样结构可以在不影响现有硬币翻转采样方法采样效率的前提下，实时处理动态有权图上的边增删操作.

进一步地，本文将硬币翻转采样方法应用至动态有权图上的随机游走概率计算中，提出了一个新颖的算法CoinFlipWalk. 本文证明，CoinFlipWalk算法可以在保证随机游走概率计算结果无偏性的前提下，同时达到近似最优的随机游走概率计算复杂度和最优的采样结构更新复杂度. 特别地，CoinFlipWalk算法同时支持无权图上的随机游走概率计算，可以在最优的时间复杂度下完成计算.

最后，本文在多个大规模动态有权图上进行了模拟实验. 实验结果显示，在达到同级别近似计算误差的前提下，CoinFlipWalk算法的运行时间和更新时间均小于既往方法，这证实了CoinFlipWalk算法相较于既往方法的优越性.

1.   预备知识

1.1   符号和计算模型

本文使用$G=(V,E)$表示所关注的有权动态图结构，其中$V$指代图$G$的节点集，${E}$指代边集. 本文用$n=\left|V\right|$表示图$G$的节点数，用$m=\left|E\right|$表示图$G$的边数. 对于图上任意节点$u$，本文用${N}_{\rm in}\left(u\right)$和${N}_{\rm out}\left(u\right)$分别表示节点$u$的入邻居节点集合和出邻居节点集合；用${n}_{u}=\left|{N}_{\rm out}\right(u\left)\right|$表示节点$u$出邻居的个数. 此外，本文用$\boldsymbol{A}$表示图$G$的邻接矩阵，用${{ A}}_{uv}$表示边$(u,v)$的权重，用${d}_{u}=\displaystyle\sum _{v\in {N}_{\rm out}\left(u\right)}^{}{{ A}}_{uv}$表示节点$u$的出度，用${{ A}}_{u}^{\mathbf{*}}$表示节点$u$出边边权的最大值.

本文使用word RAM 模型分析算法的时间复杂度. Word RAM模型是一个常用的计算模型，用以模拟真实的计算机运行情况. Word RAM模型的基本假设是：一个可随机存取的机器（random access machine）可以对一个由$w$位组成的字进行位操作，且机器寻址大小（即${2}^{w}$）与问题规模大小相匹配. 需要注意的是，随机游走概率计算问题的问题规模大小为$O(n+m)= O({n}^{2})$. 因此，为遵从Word RAM模型的基本假设，本文亦假设所用数值（例如，任一边权${{ A}}_{uv}$）都可以用$O({\rm log}\,n)$比特来表示. 此外，为了表述与分析方便，本文假设图$G$上各边权均满足${{ A}}_{uv}\ge 1$. 对于任意给定图$G$，在读图时对其上各边的边权${{ A}}_{uv}$进行统一处理${{ A}}_{uv}\leftarrow {{ A}}_{uv}/{w}_{\min}$，在随机游走概率计算全部结束后再进行逆处理${{ A}}_{uv}\leftarrow {w}_{\min}\cdot {{ A}}_{uv}$，即可保证${{A}}_{uv}\ge 1$假设成立且不改变给定图$G$的任何性质. 此处${w}_{\min}$指图$G$上各边边权的一个下界.

1.2   随机游走概率的定义与性质

给定图$G$上的节点$s$作为随机游走的源节点，给定正整数$L$作为随机游走的长度，从节点$s$出发的随机游走在其第$L$步走到图上各节点的随机游走概率向量${{\boldsymbol \pi }}_{s}^{\left(L\right)}$的定义为：

其中，$\boldsymbol{P}=\boldsymbol{A}^{\rm{T}}\boldsymbol{D}^{-1}$被称为随机游走的转移概率矩阵, ${{\boldsymbol e}}_{s}$是一个$n$维的独热向量（one-hot vector），其只有第$s$维是1，其余维为0. 对于图上任意节点$u$，${{ \pi }}_{s}^{\left(L\right)}\left(u\right)$即表示从源点$s$处出发的随机游走在其第$L$步走到节点$u$的概率. 为了表述的简洁，本文将在不引起歧义的前提下默认随机游走的源节点为$s$，并省略${{\boldsymbol \pi }}_{s}^{\left(L\right)}$和${{ \pi }}_{s}^{\left(L\right)}\left(u\right)$的下脚标$s$，将其简化为${{\boldsymbol \pi }}_{}^{\left(L\right)}$和${{ \pi }}_{}^{\left(L\right)}\left(u\right)$.

由式（1）中随机游走概率向量可得：${{\boldsymbol \pi }}^{\left(0\right)}={{\boldsymbol e}}_{s}$，即从源点$s$出发的随机游走在其第0步只能走到节点$s$. 对于任意$\mathscr{l}\in \{1, 2, … \}$，从源节点$s$出发的随机游走在其第$\mathscr{l}$步走到任意图节点$v$的概率${{ \pi }}^{\left(\mathscr{l}\right)}\left(v\right)$满足关系式：

式（2）表示从节点$s$出发的随机游走在其第$\mathscr{l}$步走到节点$v$的概率可被拆分为2项：该游走在其第$\mathscr{l}-1$步走到图节点$u$的概率，即${{ \pi }}^{\left(\mathscr{l}-1\right)}\left(u\right)$和该游走在其最后一步由节点$u$走到节点$v$的概率，即${{ A}}_{uv}/d_u$. 式（2）可被进一步拓展为关系式：

其中，$\mathscr{l}\in \left\{\mathrm{0,1}, … \right\}$，$i\in \left\{\mathrm{0,1}, … ,\mathscr{l}\right\}$.

1.3   动态有权图上的随机游走概率计算

本文关注动态有权图上随机游走概率的近似计算问题，该问题的具体定义有：

定义1. 给定动态有权图$G=(V,E)$、源节点$s\in V$、随机游走步长$L$、相对误差参数$c$、相对误差阈值$\delta$、失败概率参数${p}_{f}$，本文期望对图上所有节点$u$的随机游走概率${{ \pi }}^{\left(L\right)}\left(u\right)$进行近似估计，所得估计值${\hat{{ \pi }}}^{\left(L\right)}\left(u\right)$以不小于$1-{p}_{f}$的概率满足：

为了表述和分析的方便，本文假设相对误差$c$和失败概率参数${p}_{f}$均为常数，既往多个工作中亦作此假设.

本文期望设计近似算法以较小的计算复杂度完成定义1所要求的近似计算问题. 同时，本文关注所提近似算法所依赖的采样结构在图结构动态变化后的更新时间复杂度.

2.   现有工作

图上随机游走由于其刻画节点相关性的能力，在图算法、图机器学习领域有广泛的研究和应用. 近年来，随着研究的深入，有权图和动态图逐渐引起研究者的重视，对有权图和动态图上随机游走应用的工作也有了初步的探究. 例如，郭佳雯等人^[13]提出利用随机游走计算动态异质图的网络嵌入，张小驰等人^[14]针对子图查询问题提出基于半Markov随机游走的算法，Wu等人^[15]提出利用随机游走来辅助计算有权图的嵌入. 除了本文提出的近似算法之外，如何利用矩阵运算精确计算随机游走概率也是一个重要的问题. 例如，Jung等人^[16]提出利用分块求解的方法计算带重启的随机游走（random walk with restart）概率. 总的来说，关于有权图和动态图随机游走的研究还处于初步阶段，仍存在大量有价值的问题等待探索.

本节其余部分将对随机游走概率计算的蒙特卡洛方法作以介绍. 此外，本节对权重采样和硬币翻转采样的采样模式和代表性算法作以概述. 表1展示了采样算法在有权图、无权图上单次采样和更新所需的时间开销.

表  1  各采样方法在采样节点$\mathit{u}$的邻居的时间复杂度

Table  1.  Time Complexities for Each Sampling Method at Neighbor of Sampling Node $\mathit{u}$

采样方法	有权图上的单次 采样时间	无权图上的单次 采样时间	更新时间开销
前缀和采样[17]	$O\left(1+\mathrm{log}{n}_{u}\right)$	$O\left(1+\mathrm{log}{n}_{u}\right)$	$O\left(1+\mathrm{log}{n}_{u}\right)$
别名采样[18-19]	$O\left(1\right)$	$O\left(1\right)$	$O\left({n}_{u}\right)$
拒绝采样[20]	$O\left({n}_{u}\cdot \dfrac{{{ A}}_{u}^{*}}{{d}_{u}}\right)\star$	$O\left(1\right)$	$O\left(1+\mathrm{log}\mathrm{log}{n}_{u}\right)$
硬币翻转采样	$O\left(1+\mathrm{log}{n}_{u}\right)\star$[21]	$O\left(1\right)\star$[22-23]	$O\left(1\right)$（本文）
注：⋆ 表示期望时间复杂度.

下载: 导出CSV 

| 显示表格

2.1   随机游走概率计算的蒙特卡洛方法

现有计算随机游走概率的方法大致可以分为两类：蒙特卡洛方法（Monte Carlo method）和确定性概率传播方法（deterministic push method）. 但是，确定性概率传播方法无法保证计算结果满足估计误差. 故本文仅关注随机游走概率计算的蒙特卡洛方法，对于确定性概率传播方法暂不多作讨论.

蒙特卡洛方法^[24]是解决随机问题的经典方法，其基本思想是重复采样大量样本，将各样本中目标结果出现的次数占总采样次数的比例作为目标结果的估计值. 随机游走概率计算的蒙特卡洛方法会重复地在给定图结构上的给定源点$s$处生成随机游走. 对于图上任意节点$u$，蒙特卡洛方法将记录所有随机游走中在第$L$步走到图上节点$u$的随机游走条数，用该条数占随机游走总条数的比例作为${\pi }^{\left(L\right)}\left(u\right)$的估计值.

算法1是用于随机游走概率计算的蒙特卡洛方法的伪代码. 给定动态有权图$G=(V,E)$、源节点$s$、步长$L、$随机游走条数${n}_{r}$，算法1会在图$G$上重复产生${n}_{r}$条$L$步随机游走，每条随机游走均从节点$s$出发，算法1记录所有${n}_{r}$条随机游走中在第$L$步走到图上各节点（例如节点$u$）的条数（例如${n}_{r}\left(u\right)$），用${n}_{r}\left(u\right)/{n}_{r}$作为从源节点出发的随机游走在第$L$步走到节点$u$的概率估计. 具体而言，算法1在每次生成随机游走前先定义第0步的概率向量$\boldsymbol{\hat{\pi}}_{ }^{\left(0\right)}=\boldsymbol{e}_s$，即限制第0步从节点$s$出发. 对于每一步$l\in [0,L-1]$，首先设置下一步的概率向量${\hat{{ \boldsymbol\pi }}}_{}^{(\mathscr{l}+1)}$为零向量，再从当前概率向量${\hat{{\boldsymbol \pi }}}_{}^{\left(\mathscr{l}\right)}$中取出有非零概率的节点，从该节点出发采样该节点的邻居，被采样到的邻居即是下一步随机游走的方向. 采样结束后，更新${\hat{{\boldsymbol \pi }}}_{}^{(\mathscr{l}+1)}$中被采样邻居的概率值. 可以注意到，蒙特卡洛算法在生成随机游走时，每步均会调用采样技术从当前节点的出邻居中采样节点. 因此，随机游走各步的采样效率是影响随机游走概率计算效率的关键.

算法1. 随机游走概率计算的蒙特卡洛方法.

输入：动态有权图$G=(V,E)$，源节点$s$，步长$L$，随机游走条数${n}_{r}$；

输出：${\boldsymbol \pi }^{\left(L\right)}$的估计值.

① for $w=\mathrm{1,2},… ,{n}_{r}$ do

②　定义n维向量${\hat{{\boldsymbol \pi }}}_{}^{\left(0\right)}$并作初始化：${\hat{{\boldsymbol \pi }}}_{}^{\left(0\right)}\leftarrow {{\boldsymbol e}}_{s}$；

③　for $l=\mathrm{0,1},… ,L-1$ do

④　　定义n维向量${\hat{{ \boldsymbol \pi }}}_{}^{(\mathscr{l})}$并作初始化：${\hat{{ \boldsymbol \pi }}}_{}^{(\mathscr{l}+1)}\leftarrow {\boldsymbol{0}}$；

⑤　　for each $u\in V$ with nonzero ${\hat{{\pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)$ do

⑥　　　调用一个采样方法，从集合${N}_{\text{out}}\left(u\right)$中采 　　　样节点并更新向量${\hat{{\boldsymbol \pi }}}_{}^{(\mathscr{l}+1)}$向量对应维；

⑦　　end for

⑧　　Clear ${\hat{{ \boldsymbol \pi }}}_{}^{\left(\mathscr{l}\right)}$；

⑨　end for

⑩　${{\boldsymbol q}}_{w}\leftarrow {\hat{{\boldsymbol \pi }}}_{}^{\left(L\right)}$；

⑪ end for

⑫ ${\hat{{\boldsymbol \pi }}}^{\left(L\right)}\leftarrow \dfrac{1}{{n}_{r}}\cdot \displaystyle\sum _{w=1}^{{n}_{r}} {{\boldsymbol q}}_{w}$；

⑬ return ${\hat{{\boldsymbol \pi }}}^{\left(L\right)}$.

2.2   权重采样方法概述

现有工作普遍调用权重采样技术完成随机游走各步的采样. 权重采样技术致力于从一个给定概率分布中采样1个元素，具体地，以从图节点$u$的出邻居集${N}_{\rm out}\left(u\right)$中采样出邻居节点为例，应用于随机游走概率计算的权重采样方法会从概率分布$\left\{{{ A}}_{uv}/{d}_{u}\right\}$中选取一个出邻居节点$v$. 此处，${{ A}}_{uv}/{d}_{u}$表示节点$v$的采样概率. 因为$\displaystyle\sum _{v\in {N}_{{\rm out}\left(u\right)}}{{ A}}_{uv}/{d}_{u}=1$，故节点$u$所有出邻居的采样概率形成了一个概率分布，权重采样技术即希望从该概率分布中选出一个节点，使得任意出邻居节点被选中的概率为${{ A}}_{uv}/{d}_{u}$. 将权重采样方法运用到算法1中实现随机游走概率计算，则每次采样都会有且仅有1个邻居被采样，也即每步随机游走都只会走到1个节点.

单次权重采样方法的采样时间复杂度和动态更新时间复杂度的理论下界均为$\mathrm{\Omega}\left(1\right)$，但现有权重采样方法的采样或更新复杂度较$\mathrm{\Omega}\left(1\right)$差距较大. 这是由于，权重采样的采样结果对于各采样节点是不独立的，具有互斥性. 如果节点$v$被选中，则说明节点$u$其它出邻居都未被选中. 这种采样互斥性为权重采样的采样结构设计带来了麻烦，其要求采样结构必须足够复杂方可支持这种采样依赖. 这使得现有的权重采样方法或采样效率不高；或依赖于复杂的采样结构难以在动态图上即时更新. 上述问题成为实现动态有权图上随机游走概率高效计算的阻碍.

下文将对3种代表性的权重采样技术及其各自存在的问题作以概述.

2.2.1   别名采样

别名采样（alias sampling）是一种经典的权重采样方法^[18-19]，其构建了一个名为别名表（alias table）的采样结构来实现权重采样. 借助别名表，别名采样方法即可在静态有权图上实现$O\left(1\right)$时间的权重采样. 但是，由于别名表结构非常复杂，不支持动态图上的实时更新，一旦图结构发生变化（例如出现一条边的增删），别名采样方法都需要花费$O\left({n}_{u}\right)$的时间重建别名表，这严重阻碍了别名采样方法在动态有权图上的应用.

2.2.2   前缀和采样

前缀和采样（prefix sum sampling）方法构建了一个二叉搜索树来存储采样概率的前缀和. 借助该二叉搜索树结构，前缀和采样方法可以在动态有权图上以$O\left(1+\mathrm{log}{n}_{u}\right)$的时间复杂度完成1次采样. 但是，对于动态图上1条边的增删，前缀和采样方法需要花费$O(1+\mathrm{log}{n}_{u})$的时间更新该二叉搜索树结构，无法做到实时更新. 特别地，对于各元素采样概率近似相同的情况（无权图上随机游走的各步采样是该情况的一个极端例子），前缀和采样方法仍需要花费$O(\mathrm{log}{n}_{u})$的时间完成1次采样，这说明二叉搜索树结构的灵活性较差，无法根据元素的采样概率分布进行采样策略调整，从而造成时间浪费.

2.2.3   拒绝采样

拒绝采样（rejection sampling） ^[20]也是一种经典的权重采样方法. 以有权图上节点$u$处的随机游走采样为例，拒绝采样方法首先从${N}_{\rm out}\left(u\right)$中均匀随机地采样一个候选节点$v$，再以概率${{ A}}_{uv}/{{ A}}_{u}^{*}$接受节点$v$，以$1- {{ A}}_{uv}/ {{ A}}_{u}^{*}$的概率拒绝该候选节点，其中${{ A}}_{u}^{*}$为节点$u$出边边权的最大值. 拒绝采样方法若拒绝了该候选节点，则继续重复上述过程直至接受了1个候选节点. 可以注意到，拒绝采样方法不依赖于任何采样结构，其更新时间消耗仅来自于图结构变化后对值${{ A}}_{u}^{*}$的维护，更新时间复杂度为$O\left(1+\mathrm{log}\mathrm{log}{n}_{u}\right)$. 但是，拒绝采样完成单次采样的期望时间复杂度为$O\left(1+{n}_{u}\cdot {{ A}}_{u}^{*}/{d}_{u}\right)$，这一复杂度在边权分布高度偏斜的图上可能接近$O\left({n}_{u}\right)$，该复杂度量级较高，会造成不必要的时间浪费.

2.3   硬币翻转采样方法概述

本节将介绍硬币翻转采样定义和现有的硬币翻转采样方法.

硬币翻转采样是一种经典的采样模式，其关注如何以指定概率对各元素进行独立地采样. 具体地，以图节点$u$处的随机游走采样为例，硬币翻转采样会对节点$u$的每一个出邻居节点进行独立采样，节点$u$的每个出邻居节点都将以采样概率被采出. 如果将硬币翻转采样运用到算法1中实现随机游走概率计算，那么每次采样可能会得到0个或多个邻居，也即每步随机游走可能没有走到任何节点或走到了多个节点. 如果有1步随机游走没有走到任何节点，那么这次随机游走模拟就结束在了中间步而没有走到要求的步；如果走到了多个节点，那么下一步就需要从这些节点出发，分别采样它们的邻居并游走到所有被采样的邻居.

更为形象地，在节点$u$处的硬币翻转采样过程可以被理解为一个硬币投掷过程，硬币翻转采样在节点$u$的每个出邻居节点处都独立地抛掷一个硬币. 对于节点$u$的某一出邻居$v$，硬币翻转采样在$v$处抛掷的硬币正面向上的概率为${{ A}}_{uv}/{d}_{u}$，反面向上的概率为$1-{{ A}}_{uv}/{d}_{u}$. 硬币翻转采样将返回所有硬币抛掷时正面向上的节点作为采样结果.

对比硬币翻转采样和权重采样，硬币翻转采样的采样结果对各采样元素是独立的：节点$u$的多个出邻居节点可能被同时采样到. 但是，权重采样的采样结果是互斥的，每次权重采样只会从出邻居集合中采出1个节点. 硬币翻转采样的采样独立性为其采样结构的设计增添了很大的灵活度. 如本文第4节所示，仅需构造一个非常松散的采样结构，现有硬币翻转方法就可以同时做到高效的采样和更新.

2.3.1   硬币翻转采样的研究现状

现有的硬币翻转采样方法大多针对静态采样场景，即元素个数和采样概率始终保持不变的情况，早期Tsai等人^[25]提出了时间复杂度为$O(1+n\sqrt{\mathrm{min}\{\tilde{p},1-\tilde{p}\}})$的算法，之后Bringmann等人^[26]将单次采样的时间复杂度改进到最优的$O(1+\mu )$，达到理论下界. Guo等人^[21]针对静态影响力最大化问题提出了预先对元素概率排序的方法以获得更好的实际运行效率，相关算法的理论时间复杂度为$O(1+\mu +\mathrm{log}{n}_{u})$.

动态场景下的硬币翻转采样近年来才逐渐受到关注，所述动态场景即指元素数量和采样概率会随时间变化. Yi等人^[22]提出了动态硬币翻转采样算法ODSS，同期，Bhattacharya等人^[23]提出了SetSampler，这2个方法均同时达到了最优采样时间复杂度$O(1+\mu )$和最优更新时间复杂度$O\left(1\right)$.

但是，ODSS和SetSampler算法只能在$O\left(1\right)$时间内处理单个元素的变动（增加或删除1个元素，修改1个元素的概率），而在有权动态图上的随机游走中，节点$u$的任一边权发生变化或邻居数量发生变化，${d}_{u}$就会发生变化，从而每个邻居的采样概率${{ A}}_{uv}/{d}_{u}$都发生变化，ODSS和SetSampler在这种情况下需要花费$O\left({n}_{u}\right)$的时间才能完成一次更新.

2.3.2   BioSampling算法

BioSampling是硬币翻转采样的代表算法，对于含${n}_{u}$个元素的硬币翻转采样问题，BioSampling算法可以在$O(1+\mu+\mathrm{log}n_u)$的采样时间复杂度下完成1次采样，其中$\mu$指所有元素采样概率之和，$\mathrm{\Omega }(1+\mu )$即为该问题采样时间复杂度的理论下界. 值得注意的是，BioSampling算法目前不支持采样元素的动态增删，本文第4节将为其设计支持动态图上实时更新的采样结构.

下文将简述BioSampling算法的具体流程. 以图节点$u$处的随机游走采样为例，BioSampling算法在读图阶段提前将节点$u$的出邻居集合划分为$\lceil\mathrm{log}{n}_{u}\rceil +1$个组，其中第$i$组（$i\in [1,\lceil\mathrm{log}{n}_{u}\rceil$）包含所有满足${p}_{v}\in \left(\dfrac{1}{{2}^{i}},\left. \dfrac{1}{{2}^{i-1}}\right]\right.$的节点$v$，此处$i\in \{1, 2,… ,\lceil \mathrm{log}{n}_{u}\rceil \}$，${p}_{v}={{ A}}_{uv}/{d}_{u}$指节点$v$的采样概率；第$\lceil\mathrm{log}{n}_{u}\rceil +1$组包含所有满足${p}_{v}\in \left(0\right. ,1/{2}^{\lceil\mathrm{log}{n}_{u}\rceil }]$的节点$v$. BioSampling算法依次扫描这$\lceil\mathrm{log}{n}_{u}\rceil +1$个组，并在各组中分别进行硬币翻转采样.

在第$i$组（$i\in [1,\lceil\mathrm{log}{n}_{u}\rceil$）的硬币翻转采样中，BioSampling算法首先设置一个采样概率上界${\bar{p}}_{}=1/{2}^{i-1}$，将该采样概率上界${\bar{p}}_{}$视为该组内各节点的采样概率并进行硬币翻转采样. 具体地，先产生一个二项分布随机数$b\thicksim B({n}_{u},{\bar{p}}_{})$，再从该组内所有节点中均匀随机地采样$c$个节点. 对于每个被采样到的节点，称为候选节点$v$，以${p}_{v}/\bar p$的概率采出候选节点$v$，以$1-{p}_{v}/\bar p$的概率拒绝采出节点$v$. 可以注意到，节点$v$最终被BioSampling算法采出的概率为$\bar{p}\cdot {p}_{v}/\bar p={p}_{v}$，故BioSampling算法在各组内进行的硬币翻转采样是无偏的.

值得注意的是，在前$\lceil\mathrm{log}{n}_{u}\rceil$组中，组内采样概率的上界$\bar p$与该组内任一节点$v$的采样概率${p}_{v}$最多相差2倍. 因此在期望上，BioSampling算法在该组内采出的候选节点的数量不会超过最终被采样的节点数量的2倍. 因此，BioSampling算法在前$\lceil\mathrm{log}{n}_{u}\rceil$组的采样时间开销不超过$1+2\mu +\lceil\mathrm{log}{n}_{u}\rceil =O\left(1+\mu +\mathrm{log}{n}_{u}\right)$.

对于第$\lceil\mathrm{log}{n}_{u}\rceil +1$组的硬币翻转采样，BioSampling算法设置采样概率上界${\bar{p}}_{}=1/{n}_{u}$并产生二项分布随机数$c\thicksim B({n}_{u},{\bar{p}}_{})$，若$c < 1$则直接结束该组内的采样，否则在该组内均匀随机地采样$c$个候选节点，并以${p}_{v}/\bar{p}$的概率采出候选节点$v$. 由于${\bar{p}}_{}=1/{n}_{u}$，故该组采出候选节点的期望数不超过$1$，不影响BioSampling算法单次采样的时间复杂度.

值得注意的是，上述产生二项分布随机数$c\thicksim B({n}_{u},{\bar{p}}_{})$并作$c$次均匀随机采样的过程可以通过重复产生几何分布随机数的方法在$O\left({n}_{u}{\bar{p}}_{}\right)$的期望时间内实现. 具体地，BioSampling算法首先令$k'=0$，产生一个几何分布随机数$k\thicksim G\left({\bar{p}}_{}\right)$并采样节点$u$的第$k$个邻居作为候选节点；再令$k'=k$，产生一个几何分布随机数$k\thicksim G\left({\bar{p}}_{}\right)$ 并采样节点$u$的第$k+k'$个邻居作为候选节点；重复上述操作直至$k+k' > {n}_{u}$.

3.   动态图上可实时更新的硬币翻转采样结构

本节提出了一个支持动态更新的硬币翻转采样结构，该结构可以使BioSampling算法在$O\left(1\right)$时间内处理完动态有权图上1条边的增删. 值得注意的是，动态图上的节点增删操作可以转化为多个边增删操作（即逐一增加或删除该节点的所有邻边），边权变化操作可以转化为1次边删除操作和1次边增加操作. 因此，为了表述方便，本文只关注动态有权图上的边增删操作.

3.1   采样结构

本文在图1上各节点处构建了支持实时动态更新的硬币翻转采样结构. 以图上某一节点$u$为例，节点$u$处构建的采样结构包含出邻居节点索引表、出邻居的分组情况表、出邻居分组情况汇总表3个部分.

图  1  硬币翻转采样技术的采样结构

Figure  1.  Sampling structure of coin flip sampling technique

下载: 全尺寸图片 幻灯片

1）出邻居节点索引表

出邻居节点索引表（index table）包含节点$u$的所有出邻居的索引. 该索引表使用哈希表维护，哈希表中存储的元素为指向节点$u$各出邻居信息的指针，元素键值（key）为出邻居节点的ID. 由此，给定节点$u$某一出邻居节点$v$的ID信息，借助该索引表即可在$O\left(1\right)$时间内找到出邻居节点$v$的相关信息. 值得注意的是，该索引表结构是动态图维护的基本结构. 所有处理动态图的算法都需要建立类似的结构以实时定位到增删的边.

2）出邻居的分组情况表

如第3.3节所述，BioSampling算法会提前将待采样元素按照其采样概率分成多个组，并在各组内分别执行采样操作. 为此，本文为节点$u$的所有出邻居建立分组情况表，该表在形式上包括$\lceil\mathrm{log}{{ A}}_{u}^{\mathbf{*}}\rceil =O\left(\mathrm{log}{d}_{u}\right)$个组（group），各组均由数组（array）维护，其中第$i$组容纳所有${{ A}}_{uv}\in ({2}^{i-1},{2}^{i}]$的节点$v$及${{ A}}_{uv}$的值，此处$v$指代节点$u$的一个出邻居，$i\in [1,\lceil\mathrm{log}{{ A}}_{u}^{\mathbf{*}}\rceil ]$. 为了表述的方便，本文用${N}_{\rm out}^{\left(i\right)}\left(u\right)$指代节点$u$分组情况表中的第$i$组. 上述节点$u$的出邻居节点索引表中以$v$为键值的哈希元素，其存储的指针即指向$v$在节点$u$的分组情况表中所在的位置. 由此，给定节点$u$某一出邻居节点$v$的ID信息，借助节点$u$的出邻居节点索引表中存储的指针即可在$O\left(1\right)$时间内找到节点$v$在其分组情况细节表中存储的位置，进而读取到边${{ A}}_{uv}$的权重.

3）出邻居的分组情况汇总表

本文额外为每个图节点$u$建立一个分组情况汇总表，该汇总表由一个数组（array）和一个二进制数组（bit array）维护，用于存储节点$u$分组情况表的一些统计信息. 具体而言，该汇总表的数组结构长度为$\lceil\mathrm{log}{{ A}}_{u}^{\mathbf{*}}\rceil$，该表第$i$项存储${N}_{\rm out}^{\left(i\right)}\left(u\right)$中的元素个数和一个指向${N}_{\rm out}^{\left(i\right)}\left(u\right)$数组结构的指针. 特别地，对于所有非空的${N}_{\rm out}^{\left(i\right)}\left(u\right)$，其组的ID（即$i$）记录在该汇总表的二维数组结构中，即非空组ID对应位置标记为1，其余组为0. 由此，借助该汇总表便可完成BioSampling中所要求的操作，即将采样元素分组并在各组中依次进行采样. 具体而言，给定节点$u$，逐一访问节点$u$分组情况汇总表的二进制数组便可找到所有指向非空${N}_{\rm out}^{\left(i\right)}\left(u\right)$的指针，进而找到节点$u$分组情况表中的所有组，从而得以在各组中完成BioSampling算法所要求的采样操作.

3.2   采样结构对图更新操作的处理步骤

本节将具体说明所提采样结构处理动态图上一条边增删操作的方法.

1）边增（edge insertion）操作. 以节点$u$为例，当在图$G$上插入边$(u,v)$时，节点$u$的出邻居数增加. 因此，在边增操作后，该采样结构首先需要在节点$u$的分组情况汇总表中更新出邻居$v$所在组的信息. 为表述方便，此处用${N}_{\rm out}^{\left(i\right)}\left(u\right)$指代$v$所在的组；其次，根据汇总表中存储的指针，找到${N}_{\rm out}^{\left(i\right)}\left(u\right)$对应的数组结构，在数组的最后加入节点$v$及${{ A}}_{uv}$的信息. 进一步地，创建一个指向${N}_{\rm out}^{\left(i\right)}\left(u\right)$中节点$v$存储位置的指针，并在节点$u$的邻居节点索引表新增一个哈希元素存储该指针.

2）边删（edge deletion）操作. 以节点$u$为例，当从图$G$上删除边$(u,v)$时，节点$u$的出邻居数减少. 因此，在边删操作后，该采样结构首先需要在节点$u$的邻居节点索引表找到其存储的指针，通过该指针找到节点$v$的相关信息在${N}_{\rm out}^{\left(i\right)}\left(u\right)$的存储位置，此处用${N}_{\rm out}^{\left(i\right)}\left(u\right)$指代$v$所在的组. 然后将$v$存储在${N}_{\rm out}^{\left(i\right)}\left(u\right)$中的信息删除，再将${N}_{\rm out}^{\left(i\right)}\left(u\right)$尾项中的信息补入$v$在${N}_{\rm out}^{\left(i\right)}\left(u\right)$中的原位置，保证${N}_{\rm out}^{\left(i\right)}\left(u\right)$中元素存储的连续性，同时也需要修改尾项元素在邻居节点索引表中的键值，使其指向更新后的位置. 最后，在节点$u$的分组情况汇总表中更新${N}_{\rm out}^{\left(i\right)}\left(u\right)$组的信息.

3.3   采样结构的动态更新时间分析

本节将对所提采样结构处理动态图上边增删操作的时间复杂度，以及所提采样结构的空间复杂度进行分析.

引理1. 在处理图结构边增删的时间开销之外，所提采样结构处理单条边增删操作的时间开销为$O\left(1\right)$.

证明. 引理1的证明较为直观. 注意，所有处理动态图的算法都需要建立类似邻居节点索引表的结构以实时定位到需要增删的边. 因此，由邻居节点索引表所引起的边增删时间不在本文分析范围之内. 由此，根据第4.2节所述更新步骤，所提采样结构处理动态图上边增删的各步操作所消耗的均为$O\left(1\right)$. 证毕.

4.   动态有权图上的随机游走概率计算方法

本节将硬币翻转采样方法应用至动态有权图上的随机游走概率计算中，提出了一种新颖的高效算法CoinFlipWalk. CoinFlipWalk也是一种蒙特卡洛算法，其通过在给定图结构上重复地生成随机游走实现随机游走概率的近似估计. 但与既往方法不同的是，CoinFlipWalk在每条随机游走中都调用BioSampling这一硬币翻转采样方法进行随机游走采样. 下文将证明CoinFlipWalk可以在近似最优（忽略log因子）的时间复杂度内完成定义1中所要求的随机游走概率近似计算. 结合第4节提出的支持实时动态更新的硬币翻转采样结构，CoinFlipWalk即可同时做到近似最优的随机游走概率计算复杂度和最优的采样结构更新复杂度.

算法2. CoinFlipWalk的单步采样操作.

输入：满足${\hat{\pi }}_{}^{\left(\mathscr{l}\right)}\left(u\right) > 0$的节点$u$；

输出：更新后${\boldsymbol \pi }_{}^{(\mathscr{l}+1)}$的估计值.

① for each non-empty group ${N}_{\text{out}}^{\left(i\right)}\left(u\right)$ of $u$ do

②　Let ${p}^{\mathrm{*}}=\dfrac{{\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {2}^{i}}{{d}_{u}}$ for all $v\in {N}_{\text{out}}^{\left(i\right)}\left(u\right)$；

③　if ${p}^{\mathrm{*}}\ge 1$ then

④　　for each $v\in {N}_{\text{out}}^{\left(i\right)}\left(u\right)$ do

⑤　　　${\hat{{ \pi }}}_{}^{(\mathscr{l}+1)}\left(v\right)\leftarrow {\hat{{ \pi }}}_{}^{(\mathscr{l}+1)}\left(v\right)+\dfrac{{\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {{ A}}_{uv}}{{d}_{u}}$；

⑥　　end for

⑦　else

⑧　　定义${n}_{u}^{\left(i\right)}$ 以记录集合 in ${N}_{\rm out}^{\left(i\right)}\left(u\right)$中的节点数；

⑨　　$idx=0；$

⑩　　while $idx < {n}_{u}^{\left(i\right)}$ do

⑪　　　生成一个几何分布随机数 $rg\sim G\left({p}^{\mathrm{*}}\right)$；

⑫ $idx\leftarrow idx+rg$；

⑬　　　定义 $v$ 以表示集合 ${N}_{\rm out}^{\left(i\right)}\left(u\right)$中的第$idx$个 　　　节点数；

⑭　　　生成一个均匀随机数 $r\sim U\left(\mathrm{0,1}\right)$；

⑮　　　if $r < \dfrac{{{ A}}_{uv}}{{2}^{i}}$ then ${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)\leftarrow {\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)+1$；

⑯　　　end if

⑰　　end while

⑱　end if

⑲ end for

⑳ return ${\hat{\boldsymbol \pi }}_{}^{(\mathscr{l}+1)}$.

4.1   CoinFlipWalk算法说明

CoinFlipWalk算法的主体框架仍为蒙特卡洛方法框架，该框架已在算法1中展示. CoinFlipWalk中使用的硬币翻转采样方法具体如算法2所示，该采样方法以BioSampling为基础，且针对随机游走概率计算做了适应性调整. 将算法2嵌入算法1的行⑥即得到完整的CoinFlipWalk算法.

具体而言，给定动态权重图$G=(V,E)$、源节点$s$、随机游走步长$L$和随机游走条数${n}_{r}$，CoinFlipWalk算法从源节点$s$处生成${n}_{r}$条“广义”的随机游走，其中各条随机游走在其游走第$\mathscr{l}$步的结果都由一个$n$维向量${\hat{{ \boldsymbol\pi }}}_{}^{\left(\mathscr{l}\right)}$记录，$\mathscr{l}\in \{0,1,… \}$. $\hat{\boldsymbol{\pi}}_{ }^{\left(0\right)}$被初始化为${{\boldsymbol e}}_{s}$，其余所有$\mathscr{l}\ge 1$的$\hat{\boldsymbol{\pi\boldsymbol{ }}}_{ }^{\left(0\right)}$向量均被初始化为$n$维的全0向量. 算法2展示了CoinFlipWalk中“广义”随机游走第$\mathscr{l}$步的生成过程，其本质上是基于向量${\hat{{ \boldsymbol\pi }}}_{}^{\left(\mathscr{l}\right)}$的值对向量${\hat{{ \boldsymbol\pi }}}_{}^{\left(\mathscr{l}+1\right)}$的更新过程. 具体说明为：

对于图上各节点$u$，若${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)$非零，则依据下述步骤使用${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)$更新节点$u$的出邻居节点$v$所对应的${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$.

将节点$u$的所有出邻居节点分成多个组，其中第$i$组容纳所有${{ A}}_{uv}\in ({2}^{i-1},{2}^{i}]$的出邻居节点$v$；设第$i$组中所有节点的采样概率为${p}^{\mathrm{*}}=\left({\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {2}^{i}/{d}_{u}\right)$. 若${p}^{\mathrm{*}}\ge 1$，则更新第$i$组中所有节点$v$所对应的${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$：

反之，若${p}^{\mathrm{*}} < 1$，则从第$i$组中进行硬币翻转采样，对应算法2中的行⑧~⑭. 首先设定第$i$组中每个节点$v$的采样概率均为${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {2}^{i}/{d}_{u}$来采样候选节点. 因为此时各节点的采样概率相同，故第$i$组中的硬币翻转采样退化为二项分布采样，如第3.3节所述，二项分布采样可以用等价的几何分布采样代替，即从第$i$组中的第1个节点开始，产生几何分布随机数$k\thicksim G\left({p}_{}^{*}\right)$，其中${p}_{}^{*}={\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {2}^{i}/{d}_{u}$，返回第$k$个邻居作为第一个候选节点. 之后，继续从第$k+1$个邻居开始并重复上述操作，直至跳过第$i$组中的所有节点. 对于所有采样到的候选节点，需要保证其采样概率正确，所以还需执行算法2中的行⑫~⑭. 对每个采样到的候选节点$v$，以${{ A}}_{uv}/{2}^{i}$的概率对节点$v$所对应的${\hat{{ \pi }}}_{w}^{\left(\mathscr{l}+1\right)}\left(v\right)$进行更新：

以$1-{{ A}}_{uv}/{2}^{i}$的概率不更新${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$. 至此，CoinFlipWalk算法的一步采样结束，完成了基于向量${\hat{{ \boldsymbol\pi }}}_{}^{\left(\mathscr{l}\right)}$的值对向量${\hat{{ \boldsymbol\pi }}}_{}^{\left(\mathscr{l}+1\right)}$的更新过程. 根据算法1的框架，从$\mathscr{l}=0$开始重新上述过程$L$次，即可完成一次随机游走模拟对概率向量$\hat{\boldsymbol{\pi\boldsymbol{ }}}_{ }^{\left(L\right)}$的更新. 重复模拟${n}_{r}$次，即可得到4.2节所分析的有理论保证的估计结果.

4.2   理论分析

本节将对 CoinFlipWalk 算法的理论性质进行分析并给出详细证明.

4.2.1   无偏性

本节首先证明CoinFlipWalk算法计算得到的随机游走概率是无偏的. 具体如引理1所示.

引理1. 对于$\forall u\in V$、$\forall\mathscr{l}\in [0,L]$，将算法1的行⑥替换为算法2后，执行算法1所返回的${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}$满足：

证明. 引理1可通过数学归纳法证明. 首先，考虑起始状态：${\hat{{\boldsymbol \pi }}}_{}^{\left(0\right)}={{\boldsymbol e}}_{s}$，故对于$\forall u\in V$，有：$\mathbb{E}\left[{\hat{{ \pi }}}_{}^{\left(0\right)}\left(u\right)\right]= {{ \pi }}^{\left(0\right)}\left(u\right)$. 进一步地，假设对某$\mathscr{l}\in [0,L-1]$，$\mathbb{E}\left[{\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\right]= {{ \pi }}^{\left(\mathscr{l}\right)}\left(u\right)$对$\forall u\in V$成立，由于：

因此，$\mathbb{E}\left[{\hat{{ \pi }}}^{\left(\mathscr{l}+1\right)}\left(v\right)\right]=\mathbb{E}\left[\mathbb{E}\left[{\hat{{ \pi }}}^{\left(\mathscr{l}+1\right)}\left(v\right)\mid {\hat{{\boldsymbol\pi }}}^{\left(\mathscr{l}\right)}\right]\right]$. 进而有，

代入$\mathbb{E}\left[{\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\right]={{ \pi }}^{\left(\mathscr{l}\right)}\left(u\right)$的假设和式（2），可得：

下文将具体说明式（4）成立的原因. 回顾算法2，对于第$i$组中所有节点$v$，若$\left({\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right){2}^{i}/{d}_{u}\right)\ge 1$，则算法1将${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$确定地增加${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {{ A}}_{uv}/{d}_{u}$. 若$\big({\hat{{ \pi }}}_{}^{(\mathscr{l})}(u)\cdot {2}^{i}/ {d}_{u}\big) <1$，则算法1以${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {2}^{i}/{d}_{u}$的采样概率从第$i$组中采样出节点$v$，再以${{ A}}_{uv}/{2}^{i}$的概率将${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$增加$1$. 因此，当$\left({\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {2}^{i}/{d}_{u}\right) < 1$时，${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$的期望增加值仍为${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot \dfrac{{2}^{i}}{{d}_{u}}\cdot \dfrac{{{ A}}_{uv}}{{2}^{i}}={\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot \dfrac{{{ A}}_{uv}}{{d}_{u}}$，故式$（4）$可得，从而引理1得证. 证毕.

4.2.2   单次采样时间

本节进一步分析CoinFlipWalk完成单次硬币翻转采样的时间开销.

引理2. 对于$\forall u\in V$、$\forall\mathscr{l}\in [0,L-1]$，算法2的期望时间复杂度上界为$O\left(1+\mathrm{log}{n}_{u}\right)$.

证明. 令$Cost^{\left(\mathscr{l}\right)}\left(u\right)$表示在节点$u$处执行1次算法2的期望时间复杂度，令$Cost^{\left(\mathscr{l}\right)}(u,v)$表示算法2中更新${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$所花费的时间，根据算法2有：

其中，$\mathrm{log}{{ A}}_{u}^{*}$项与本文第4节所提采样结构有关，即要完成节点$u$处的硬币翻转采样，需要至少花费$O\left(1\right)$时间遍历节点$u$分组情况细节表中的每个非空组，所以算法2最多需要花费和组数相同的时间遍历非空组，故算法2的时间复杂度中包含$O\left(\mathrm{log}{{ A}}_{u}^{*}\right)$. 式（5）中$\displaystyle\sum _{v\in {N}_{{\rm out}\left(u\right)}}Cost^{\left(\mathscr{l}\right)}\left(u,v\right)$表示算法2更新节点$u$的所有邻居节点$v$的${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$所花费的时间，该时间即为算法2在各个非空组中的采样时间之和.

值得注意的是式（5）中的$\mathrm{log}{{ A}}_{u}^{*}$项可被进一步减小为$\mathrm{log}{n}_{u}$. 具体而言，设第$i$组中所有节点的采样概率为${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {2}^{i}/{d}_{u}$，算法2在扫描节点$u$的分组情况细节表时，可以$1-1/{n}_{u}$的概率跳过所有${2}^{i}/{d}_{u}$小于$1/{n}_{u}^{2}$的组. 这是因为，节点$u$最多有${n}_{u}$个邻居满足${2}^{i}/{d}_{u}\le 1/{n}_{u}^{2}$，这些邻居中至少有1个被采样到的概率不超过$1-{\left(1-\dfrac{1}{{n}_{u}^{2}}\right)}^{{n}_{u}}$，进而不超过$\dfrac{1}{{n}_{u}}$，也即有大于$1-1/{n}_{u}$的概率是不需要采样所有这些满足条件的邻居节点的. 所以，算法2在运行之初可以产生一个$\left[\mathrm{0,1}\right]$之间的均匀随机数$r$，若$r\le 1-1/{n}_{u}$，则跳过所有满足$\dfrac{{2}^{i}}{{d}_{u}}\le 1/{n}_{u}^{2}$的分组；否则逐一扫描这些分组中的邻居节点以判断是否需要采样，该过程最多消耗$O\left({n}_{u}\right)$的时间. 注意，满足$\dfrac{{2}^{i}}{{d}_{u}}\ge 1/{n}_{u}^{2}$的分组最多有$2\mathrm{log}{n}_{u}$组. 因此，$\mathrm{log}{{ A}}_{u}^{*}$项从期望上可被减小为：$2\mathrm{log}{n}_{u}+{n}_{u}\cdot \dfrac{1}{{n}_{u}}= O(1+\mathrm{log}{n}_{u})$. 实现时，可以始终维护首个采样概率下界超过$1/{n}_{u}^{2}$的组ID，记为$j$，若产生的随机数$r\le 1-1/{n}_{u}$，直接从第$j$组开始扫描.

另一方面，根据算法2，若${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right){2}^{i}/{d}_{u}\ge 1$，则有$Cost^{\left(\mathscr{l}\right)}\left(u,v\right)=1\le \left({\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right){2}^{i}/{d}_{u}\right)$；若$\left({\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {2}^{i}/{d}_{u}\right) < 1$，则以${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {2}^{i}/{d}_{u}$的概率有$Cost^{\left(\mathscr{l}\right)}\left(u,v\right)=1$. 因此，

综上，$Cost^{\left(\mathscr{l}\right)}\left(u\right)$的条件期望可被表示为：

由于$\mathbb{E}\left[Cost^{\left(\mathscr{l}\right)}\left(u\right)\right]=\mathbb{E}\left[\mathbb{E}\left[Cost^{\left(\mathscr{l}\right)}\left(u\right)\mid {\hat{{\boldsymbol\pi }}}^{\left(\mathscr{l}\right)}\right]\right]$成立，结合$Cost^{(\mathscr{l})}(u)$的条件期望，并由引理1可进一步得到式（5）的上界：

注意$\displaystyle\sum _{v\in {N}_{{\rm out}\left(u\right)}}{{ \pi }}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot \dfrac{2{{ A}}_{uv}}{{d}_{u}}\le \displaystyle\sum _{v\in {N}_{{\rm out}\left(u\right)}}{{ \pi }}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)\le 1$始终成立，故引理2得证. 证毕.

4.2.3   CoinFlipWalk时间复杂度

最后，本节证明CoinFlipWalk算法可以在近似最优的时间复杂度下完成定义1中所要求的随机游走概率计算问题. 定理3给出了满足定义1所需的随机游走次数. 结合定理3和引理2可得CoinFlipWalk的计算复杂度为$\tilde{O}({L}^{2}/\delta )$. 此处$\tilde{O}$为省略log项的大O表示法. 值得注意的是，该问题的平凡下界为$\Omega \left(1/\delta \right)$，$L$常被设置为常数. 因此，如果忽略log系数，CoinFlipWalk算法的期望时间复杂度达到了理论下界. 以下为引理3及其证明.

引理3. 当$n_r\ge\dfrac{L}{p_f\cdot c_{ }^2\delta}$时，CoinFlipWalk所返回的$\hat{\boldsymbol{\pi}}^{\left(L\right)}$对于任意节点$u\in V$以不小于$1-{p}_{f}$的概率满足：

证明. 为证明引理3，先对$\forall u\in V$求出${\hat{{ \pi }}}^{\left(L\right)}\left(u\right)$方差的上界，再借助切比雪夫不等式得到随机游走条数${n}_{r}$的上界.

回顾算法2，对于采样概率小于1的组$i$内的邻居$v$，${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$以${\hat{{ \pi }}}_{}^{\left(\mathscr{l}\right)}\left(u\right)\cdot {{ A}}_{uv}/{d}_{u}$的概率增加1，否则不增加. 令${{r}_{}^{\left(\mathscr{l}\right)}}_{u}\left(v\right)$为算法2在节点$u$处进行硬币翻转采样时${\hat{{ \pi }}}_{}^{\left(\mathscr{l}+1\right)}\left(v\right)$的增量，则有：

进而有，

因为${\hat{{ \pi }}}^{\left(\mathscr{l}+1\right)}\left(v\right)=\displaystyle\sum _{u\in {N}_{\rm in}\left(v\right)}{r}_{u}^{\left(\mathscr{l}\right)}\left(v\right)$，以及硬币翻转采样的独立性，有：

进而有：

进一步地，由于方差满足全方差公式（total variance law），故有：

因此，${Var}\left[{\hat{{ \pi }}}^{\left(L\right)}\left(u\right)\right]$可被表示为：$\mathbb{E}\left[Var\left[{\hat{{ \pi }}}^{\left(L\right)}\left(u\right)\mid \right.\right. \left. \left. {\hat{{ \pi }}}^{(L-1)}\right]\right]+{Var} \left[\displaystyle\sum _{x\in {N}_{\rm in}\left(u\right)}\dfrac{{{ A}}_{xu}\cdot {\hat{{ \pi }}}^{(L-1)}\left(x\right)}{{d}_{x}} \right]$. 注意$\dfrac{{{ A}}_{xu}}{{d}_{x}}={\hat{{ \pi }}}_{x}^{\left(1\right)}\left(u\right)$，所以${Var} \left[\displaystyle\sum _{x\in {N}_{\rm in}\left(u\right)} \dfrac{{{ A}}_{xu}\cdot {\hat{{ \pi }}}^{(L-1)}\left(x\right)}{{d}_{x}}\right]$=${Var}\left[\displaystyle\sum _{x\in {N}_{\rm in}\left(u\right)} {{ \pi }}_{x}^{\left(1\right)}\left(u\right)\cdot {\hat{{ \pi }}}^{(L-1)}\left(x\right)\right]$，这也可以写作${Var}\left[\displaystyle\sum _{x\in V}{{ \pi }}_{x}^{\left(1\right)}\left(u\right)\cdot {\hat{{ \pi }}}^{(L-1)}\left(x\right)\right]$. 根据全方差公式，上述推导可以得到更一般的结论：对于$\mathscr{l}\in \left[1,L\right]$，有

由式（6）推导可知，调用一次全方差公式，可得：

反复调用式（6） ，$Var\left[\hat{\boldsymbol{\mathit{\pi}}}^{\left(L\right)}\left(u\right)\right]$最终可被写为：

注意，对于$\forall u\in V$， ${\hat{{ \pi }}}^{\left(\mathscr{l}\right)}\left(u\right)$的取值在给定${\hat{{\boldsymbol \pi }}}^{(\mathscr{l}-1)}$后是独立的. 因此，

将式（8）代入式（7）即可得到：

由此，$Var\left[{\hat{{ \pi }}}^{\left(L\right)}\left(u\right)\right]$可被写为：

代入 ${\hat{{ \pi }}}^{(\mathscr{l}-1)}\left(y\right)$ 的无偏性，即有：

由此得到${\hat{{ \pi }}}^{\left(L\right)}\left(u\right)$方差的上界.

进一步地，将该方差界代入切比雪夫不等式（Chebyshev’s inequality）中：

并将其代入所得方差界，即有：

因此，当${n}_{r}\ge \dfrac{L}{{p}_{f}\cdot \delta \cdot {c}^{2}}$，CoinFlipWalk计算得到的随机游走概率近似向量$\hat{\boldsymbol{\pi}}^{\left(L\right)}$满足定义1所示的误差要求. 本引理由此得证. 证毕.

5.   实验结果分析

本节展示CoinFlipWalk在大规模真实图和合成图上的实验结果. 实验均用C++实现，使用g++编译器且开启了O3优化进行编译. 本节所涉实验均在一台配备有Intel 2.20 GHz CPU和754 GB内存的Linux机器上完成.

5.1   实验设置

5.1.1   基线方法

本文将CoinFlipWalk与2种基线方法PrefixWalk和RejectionWalk进行了比较，这2种方法也是蒙特卡洛方法，其在随机游走各步分别调用前缀和采样方法以及拒绝采样方法. 具体而言，CoinFlipWalk算法在算法1的行⑥调用算法2，而PrefixWalk和 RejectionWalk算法在算法1的行⑥分别调用了前缀和采样方法以及拒绝方法. 由于别名采样方法仅支持静态加权图且每次边插入/删除操作都必须完全重建索引结构，因此本文在实验中未包括别名方法. 此外，本文采用幂方法计算有权图上随机游走概率的真值. 幂方法初始化${{\boldsymbol \pi }}^{\left(0\right)}={\boldsymbol e}_{s}$，并重复计算$\boldsymbol{\pi}^{(i+1)}\leftarrow\boldsymbol{P}\boldsymbol{\pi}^{\left(i\right)}$，直到得出$\boldsymbol{\pi}^{\left(L\right)}$. 在本文中，$L$为常数整数，在实验中除非另有说明，设$L$=10.

5.1.2   数据集

本文在4个真实有权图和2个合成图上进行了实验. 为了方便参考，表2总结了6个真实图数据集的统计信息. 其中，cos²φ反映有权图的不平衡性，cos²φ的值越小，图权重分布的不平衡性越大. 真实图分别为：Tags (TAG)^[27]，Threads (TH)^[27]，Reddit (RR)^[28-29]和Colisten (CS)^[30]. 数据集TAG和TH都源自问答网站Stack Overflow. 具体来说，TAG和TH数据集中的节点代表问题标签，边权表示2个标签在同一个问题系列中同时出现的频率. 数据集RR是根据Reddit网站上的评论构建出的网络，其中的节点代表用户，边的权重代表用户评论用户的帖子的次数. 数据集CS是一个歌曲数据集，歌曲与图节点对应，边权是2首歌曲出现相似音乐片段的次数. 合成图包括了2个基于子图模式（motif）生成权重的有权图，即IndoChina (IC)^[31-32]和Twitter (TW)^[33-34]. 这2个数据集本身是真实世界中的无权图. 本文参考已有文献[35-36]将这2个数据集转化为基于子图模式的有权图，对于每一条边，设置其参与的特定子图模式，如一个三角形子图的数量作为权重.

表  2  数据集参数

Table  2.  Parameters of Datasets

数据集	节点数$n$	边数$m$	$\dfrac{m}{n}$	${\mathrm{c}\mathrm{o}\mathrm{s}}^{2}\varphi$
TAG	49 945	$\mathrm{8\;294\;604}$	166	0.27
RR	$\mathrm{8\;396\;162}$	$\mathrm{870\;580\;840}$	103	0.90
TH	$\mathrm{2\;321\;767}$	$\mathrm{42\;012\;344}$	18	0.97
CS	$\mathrm{3\;604\;308}$	$\mathrm{3\;854\;964\;026}$	1070	0.29
IC	$\mathrm{7\;414\;768}$	$\mathrm{295\;191\;370}$	40	0.31
TW	$\mathrm{41\;652\;230}$	$\mathrm{1\;687\;583\;684}$	40	0.34

下载: 导出CSV 

| 显示表格

5.1.3   测试方法和评价指标

本文根据${n}_{u}/m$的分布随机生成10个查询节点. 所有实验结果均是每种方法以随机生成的10个查询节点为源节点进行随机游走的平均结果. 对于更新时间的测试，我们从原始图中随机删除100 000条边，然后将这100 000条边重新插入图中. 对于每种方法, 我们计算200 000次边更新操作中的平均更新时间开销作为该方法的更新开销. 对于真实图上各方法随机游走的近似误差，我们的评价指标包括Precision@50、MaxError. 其中，Precision@k是衡量近似精度的一个流行指标，定义为这些近似方法返回的结果中概率估计值最大的前$k$个节点与通过基准真值得到的实际概率最大的前$k$个节点相符的百分比. 实验时，可以通过调整每种方法中的误差参数$\delta$值来达到所需的Precision@k. $\delta$值越小，现有方法模拟的随机游走次数就越多. 因此$\delta$值越小, 每种方法在随机游走模拟过程中花费的查询时间就越长，达到的Precision@k值也就越高. 本文设置$k=50$. MaxError是另一个常用的衡量指标，该指标量化了所求的近似结果与真值相比的近似误差，定义为$MaxError=\max\limits_{u\in V} \left|{\pi }^{\left(L\right)}\left(u\right)-{\hat{\pi }}^{\left(L\right)}\left(u\right)\right|$.

5.2   查询时间与更新时间开销的相关分析

本节对各方法在动态有权图上计算和更新随机游走概率的相对大小关系进行了对比，图2具体展示了CoinFlipWalk及其他方法在达到相近Precision@50结果时的随机游走概率计算时间（查询时间）和动态更新时间，各方法名旁显示的数字为各方法所达到的Precision@50结果. 根据图2所示结果，PrefixWalk算法在各数据集上的随机游走概率动态更新时间普遍较高，RejectionWalk算法的随机游走概率动态更新时间较低但计算时间较高，而本文所提的CoinFlipWalk算法的随机游走概率计算和更新时间均较低. 特别地，由此体现出CoinFlipWalk的优越性.

图  2  在真实有权图上各方法的更新时间随查询时间的变化

Figure  2.  Update time varies with query time on real-world weighted graphs in different methods

下载: 全尺寸图片 幻灯片

5.3   查询时间与近似精度的相关分析

本节对各方法的随机游走概率计算效率进行了对比，图3具体展示了真实有权图上CoinFlipWalk及其他算法在随机游走概率计算中的Precision@50与查询时间的变化关系，图中各算法对应查询时间从左到右依次对应取0.1，0.01，…，10^–5，查询时间越小，各算法所产生的随机数量越多，算法运行时间越长. 值得注意的是，相比于PrefixWalk和RejectionWalk，CoinFlipWalk在各数据集上都可以在更短的计算时间内达到更高的Precision@50，这体现出CoinFlipWalk的优越性. 此外，图4展示了真实有权图上CoinFlipWalk及其他方法在随机游走概率计算中所得结果的MaxError与查询时间的变化关系. 类似地，相比于PrefixWalk和RejectionWalk，CoinFlipWalk在各数据集上都可以在更短的查询时间内达到更小的MaxError，这再次体现出CoinFlipWalk的计算高效性.另外值得注意的是，从图3和图4中可以看出，PrefixWalk在边权分布不均（即cos²$\varphi$较小的图数据集）的图上较在边权分布较为均匀（即cos²$\varphi$较大的图数据集）的图上的计算效率更高.与之相对地，RejectionWalk在边权分布均匀的图上较在边权分布较为不均的图上的计算效率更高.进一步地，本节还对CoinFlipWalk及其他算法在合成图上的随机游走计算效率进行了对比. 其中，图5和图6分别展示了CoinFlipWalk，PrefixWalk和Reject-ionWalk算法在2个合成图数据集上进行随机游走概率计算时所得Precision@50和MaxError与查询时间的变化关系.与前文所述类似，CoinFlipWalk在各数据集上都可以在更短的计算时间内达到更大的Precision@50和更小的MaxError.进一步地，本文还分别将CoinFlipWalk，PrefixWalk和RejectionWalk算法计算所得的随机游走概率近似结果应用于局部社区发现（local clustering）的场景中，按照sweep方法^[35]基于随机游走概率近似结果发现给定源节点周围的社区，并计算所得社区的导度（一种社区质量衡量指标，导度值越低，所得社区质量越好）. 图7展示了CoinFlipWalk，PrefixWalk和RejectionWalk三个算法所得社区的导度与算法计算随机游走概率时间之间的相对变化关系，可以看出，CoinFlipWalk可以在各数据集上用更短的查询时间取得更小的导度，这再次证明CoinFlipWalk的高效性，也体现出CoinFlipWalk的实际应用价值.

图  3  在真实有权图上各方法的Presicion@50随查询时间的变化

Figure  3.  Presicion@50 of different methods varies with query time on real-world weighted graphs

下载: 全尺寸图片 幻灯片

图  4  在真实有权图上各方法的MaxError随查询时间的变化

Figure  4.  MaxError of different methods varies with query time on real-world weighted graphs

下载: 全尺寸图片 幻灯片

图  5  在合成有权图上各方法的Precision@50随查询时间的变化

Figure  5.  Precision@50 of different methods varies with query time on synthetic weighted graphs

下载: 全尺寸图片 幻灯片

图  6  在合成有权图上各方法的MaxError随查询时间的变化

Figure  6.  MaxError of different methods varies with query time on synthetic weighted graphs

下载: 全尺寸图片 幻灯片

图  7  在合成有权图上各方法的导度随查询时间的变化

Figure  7.  Conductance of different methods varies with query time on synthetic weighted graphs

下载: 全尺寸图片 幻灯片

5.4   更新时间分析

本节对各方法的随机游走概率更新时间进行了对比. 其中，图8展示了各方法在真实有权图上的随机游走概率更新时间，图9展示了各方法在合成有权图上的随机游走概率更新时间. 由图8和图9可以看出，本文所提CoinFlipWalk算法在各数据集上所需的更新时间均小于PrefixWalk和RejectionWalk算法，这与本文理论分析所得结论一致，证明了CoinFlipWalk算法的高效性. 此外，由图8和图9还可看出，PrefixWalk 的更新时间开销显著大于 CoinFlipWalk 和 RejectionWalk，这一点亦与本文分析吻合，即前缀和方法每次边插入/删除操作都需要O(logn)的更新时间复杂度.相比之下，CoinFlipWalk 和 RejectionWalk 的更新时间复杂度仅为O(1)和O(loglogn).

图  8  在真实有权图上各方法的更新时间

Figure  8.  Update time of different methods on real-world weighted graphs

下载: 全尺寸图片 幻灯片

图  9  在合成有权图上各方法的更新时间

Figure  9.  Update time of different methods on synthetic weighted graphs

下载: 全尺寸图片 幻灯片

6.   结　　论

本文研究了动态有权图上随机游走概率的近似计算问题，提出了一种基于硬币翻转采样的随机游走概率算法. 相比于传统的基于权重采样的随机游走概率计算方法，本文所提方法在随机游走概率计算复杂度和采样结构的更新复杂度两方面都具有优势. 未来或可基于2个问题展开研究：1）硬币翻转采样技术是否可以进一步优化？2）单点随机游走计算也是相关领域的核心问题之一，所提方法是否可以应用于单点随机游走概率上的高效计算？

作者贡献声明： 王涵之负责算法设计；易璐负责实验验证；魏哲巍和甘骏豪负责理论优化；袁野、文继荣和杜小勇提出论文修改意见. 王涵之和易璐是共同第一作者.