介数中心度（betweenness centrality, BC）通过计算经过顶点的最短路径数来确定顶点在图结构中的重要程度，是图顶点重要性计算的一种流行方式，最早由Freeman^[1]提出. 顶点的介数中心度越大，则顶点对图中其他顶点的控制能力越强，介数中心度算法被广泛应用于社会网络分析^[2-3]、蛋白质交互网络影响预测^[4]、社区发现^[5]等.

现有的介数中心度算法主要聚焦于普通图，然而现实场景中建模的图常为时态图，边上带有时态信息. 例如电子邮件网络^[6]边上附带有邮件发送/接收时间；社交网络^[7]边上附带有人与人的接触时间.图1是一个电子邮件网络示例，顶点表示用户，边表示用户之间发送/接收邮件的关系，边上时间戳表示邮件发送/接收或转发/接收的时间. 具体地，用户1在3个不同时间点向用户2共发送了3份邮件，而后用户2将邮件1和邮件2转发给用户3和用户4，并将邮件3转发给用户5.电子邮件网络上介数中心度计算有助于精确推断网络结构^[8]、用户活跃程度等. 图2是一个农贸市场接触网络示例，顶点表示供货商、商户、消费者；边表示他们之间的接触时刻. 具体地，消费者1和消费者2在商户1消费；消费者3在商户1和商户2消费. 蔬菜供货商为商户1送货，水果供货商为商户1和商户2送货. 在传染病疫情暴发时，接触网络上介数中心度计算有助于识别超级传播者^[9]，控制传染病的传播.

图  1  电子邮件网络示例

Figure  1.  An example of email network

下载: 全尺寸图片 幻灯片

图  2  农贸市场接触网络示例

Figure  2.  An example of agricultural market contact network

下载: 全尺寸图片 幻灯片

普通图上介数中心度计算忽略了重要的时态信息，而时态信息对于信息流传播扩散^[10]具有重要作用. 鉴于此，本文研究的时态图顶点介数中心度计算方法与已有的普通图介数中心度计算方法相比，时态图包含时态信息，且1对顶点之间存在多条被不同时间戳标记的边，因此时态图上介数中心度计算难度更大. 具体原因可分为2方面：

1）时态图顶点介数中心度需要根据时态最短路径计算，而时态最短路径的定义方法多样，且计算时需要考虑时态边之间的时序依赖关系. 例如图3（a）给出的时态图示例中，边上的值为时间戳信息. 如果不考虑时态信息，则a到d的最短路径为$a \to c \to d$，但实际上a经c到d是不可达的，因为a到c的时间为4 （> 3）. 所以计算时态最短路径时需要考虑时态边与边之间的时序依赖关系，记a到d的时态最短路径为$a \to b \to c \to$d.

图  3  时态图和普通图示例

Figure  3.  Examples of temporal graph and general graph

下载: 全尺寸图片 幻灯片

2）针对时态图，普通图上介数中心度的计算理论与方法已不再适用，需要设计全新的理论与方法. 这是因为普通图中的介数中心度计算方法主要根据Brandes算法^[11]设计，Brandes算法有效的关键理论是最短路径的子路径依然是最短路径，即最优子结构特性. 然而时态最短路径并不满足此特性. 例如图3（b）给出的普通图中，$a \to c \to d$是1条最短路径，则子路径$a \to c$一定是最短路径. 而在时态图中（如图3（a）所示），时态最短路径$a \to b \to c \to d$的子路径$a \to b \to c$很明显不是时态最短路径，a到c的时态最短路径为$a \to c$.

为了解决这2个难点，本文根据时态边之间的时序依赖关系，定义了严格（时态递增）和非严格（时态非递减）2种时态路径类型，提出了基于消息传播的2阶段迭代计算框架以高效计算时态图顶点介数中心度. 其中，第1阶段采用自顶向下的广度优先遍历方式计算时态最短路径；第2阶段采用自底向上的方式计算顶点的后继节点和孩子节点对其介数中心度的贡献值，并设计了基于消息传播机制的迭代累积计算方法. 为了提高效率和可扩展性，实现了基于OpenMP（open multiprocessing）框架的多线程并发算法FTBC（fast temporal betweenness centrality）. 概括而言，本文的主要贡献有3点：

1）提出了自顶向下和自底向上结合的2阶段计算框架，并设计了基于消息传播机制的迭代累积计算方法以高效计算时态图顶点介数中心度.

2）提出了基于OpenMP框架的多线程并发算法FTBC，以提高时态图顶点介数中心度计算的效率和可扩展性，并理论分析了算法的复杂度.

3）基于8个真实的时态图数据集进行了大量的实验评估，验证了多线程FTBC算法相比目前流行的方法计算性能更优，可扩展性更强.

1.   相关工作

本节分别概述已有的普通图和时态图介数中心度计算方法.

1.1   普通图介数中心度计算方法

针对普通图，Brandes^[11]推导了经典的成对依赖、迭代计算理论，并基于此提出了精确介数中心度计算方法. 算法空间复杂度为$O(n + m)$，时间复杂度为$O(nm)$（无权图）或O（nm+n²log₂ n）（有权图），其中$n$和$m$分别表示顶点数和边数.Erdős等人^[12]提出了Brandes++算法，该算法采用分治策略，利用图基础结构加速计算.Sariyüce等人^[13]提出了BADIOS算法，该算法针对无向无环图，先将一些特殊顶点压缩，并利用顶点、桥边将图划分为多个子图，而后在子图上计算顶点的介数中心度.Baglioni等人^[14]提出利用图的拓扑特征加速计算. 但文献[12-14]所提算法的最坏时间复杂度仍与Brandes算法^[11]相同.

为了提高介数中心度计算效率，研究人员提出了近似算法. 为了避免计算所有顶点对之间的最短路径，近似算法的总体思想是基于采样方法计算部分顶点之间的最短路径，并基于此估算所有顶点的介数中心度. 为了保证结果质量，近似算法通常理论推导采样数或迭代更新介数中心度值，直到满足预置的终止条件.Brandes等人^[15]提出了基于顶点采样的方法，其利用霍夫丁不等式^[16]估计误差概率.Riondato等人^[17]提出了基于最短路径采样的方法RK，其利用顶点直径（vertex diameter, VD）和VC（vapnik-chervonenkis）维^[18]估算达到要求精度所需的最少样本数.RK首先采样一对顶点，然后再基于采样顶点对的一条最短路径而不是所有最短路径来近似计算顶点的介数中心度. 此外，由于计算精确的顶点直径需要所有顶点对之间的最短路径，此操作非常耗时，因此RK随机采样一个顶点，并根据该顶点到图中其他顶点的单源最短路径距离估计顶点直径. 基于RK，Borassi等人^[19]提出了KADABRA算法，其采用双向广度优先搜索方式来减少最短路径的采样时间，在介数中心度估计时允许为每个顶点设置不同的概率置信度.Riondato等人^[20]提出了ABRA算法，其利用渐进式随机抽样，基于拉德马赫平均值和伪维估计采样数. Cousins等人^[21]提出了Bavarain算法，其采用蒙特卡洛经验拉德马赫平均值^[22]估计更加严格的样本数上界，进而保证结果精度.Pellegrina等人^[23]提出了SILVAN算法，与Bavarain算法相同的是，SILVAN同样采用蒙特卡洛经验拉德马赫平均值估计样本数上界；与Bavarain不同的是，Bavarain仅适用于均匀边界，而SILVAN可以用于均匀和非均匀边界.

然而，无论是普通图介数中心度精确计算方法还是近似计算方法，方法有效的关键理论是最优子结构性质，而时态图不满足最优子结构特性，因此无法直接扩展到时态图.

1.2   时态图介数中心度计算方法

一些工作将时态图看作是一系列图的镜像，即动态图，而后研究动态图介数中心度计算.Lee等人^[24]将动态图分解为连通分量后再进行BC值更新计算以减少搜索空间. Green等人^[25]提出最短路径树结构加速动态图BC值的更新.Kourtellis等人^[26]将最短路径树结构压缩存储以减少内存占用.Kas等人^[27]拓展了动态APSP（all pairs of shortest path）算法^[28]以实现BC值动态更新.Bergamini等人^[29]提出了半动态的近似计算方法以支持顶点和边的插入操作.Hayashi等人^[30]提出了全动态的近似算法以支持边和顶点的插入和删除. 然而文献[25-30]算法均将时态图视为图镜像集合，没有考虑时态边之间的时序依赖关系，本质上还是基于普通图的最优子结构性质设计的方法.

与本文研究问题最相似的工作为Buß等人^[31]提出的时态图顶点介数中心度精确计算算法，其首先构建顶点的前置图，使得前置图上满足最优子结构性质，而后将Brandes算法理论扩展，运用在前置图上进行计算. 另一篇较相似的工作为Tsalouchidou等人^[32]提出的算法，其将路径长度与持续时间结合起来作为时态最短路径定义标准，并在时态图定长静态窗口上计算介数中心度精确值. 本文实验测试了Buß等人^[31]提出的算法，发现当时态图边数为3万时，其在内存16GB的机器上已无法运行.Tsalouchidou等人^[32]提出的算法对于时间离散程度较大的时态图而言，需要设置较大的静态窗口值，导致复制静态窗口时耗费了大量的时间. 可见，文献[31-32]工作存在计算效率较低、可扩展性差的问题.

2.   问题定义

本节主要介绍时态图、时态路径、时态最短路径、时态介数中心度的相关概念，并给出问题定义.

定义1. 时态图. 本文定义的时态图既可以是有向的，又可以是无向的，表示为$G = (V,E,T)$. 其中，$V$表示顶点集合；$E$表示时态边集合；$T$为时间戳集合. 时态图中2点之间可以有多条时态边. 具体地，时态边e_i = （u_i, v_i,, t_i）表示顶点u_i与v_i之间的事件发生时间为t_i，t_i∈T.

定义2. 严格与非严格时态路径. 从顶点u到v_k的时态路径表示为p = u$\xrightarrow{{{t_1}}}$v₁$\xrightarrow{{{t_2}}}$v₂$\xrightarrow{{{t_3}}}$v₃…v_k-1$\xrightarrow{{{t_k}}}$v_k. 其中，p的第1条边e₁ = （u, v₁, t₁），第i条边e_i = （v_i−1, v_i, t_i） （1 < i < k），第k条边e_k = （v_k−1, v_k, t_k）使得对于任意的$1 \leqslant i < k$，t_i $\leqslant$ t_i+1. 特别地，当$k = 1$，p也是1条时态路径. 进一步地，如果对于任意的$1 \leqslant i < k$，t_i $<$ t_i+1，则p称为严格时态路径；否则p称为非严格时态路径. 令|p|=k表示时态路径的长度.

定义3. 时态最短路径. 给定从顶点u到v的时态路径p_s，如果不存在从u到v的其他时态路径p满足路径长度$|p| \lt |p_s|$，则p_s是时态最短路径. 特别地，当|p_s| = 1时，p_s也是一条时态最短路径.

定义4. 时态介数中心度. 给定时态图$G =$ $(V, E, T)$，令σ_sf表示顶点s到顶点f的时态最短路径数目；σ_sf （v）表示由顶点s经顶点v到顶点f的时态最短路径数目，则对于所有顶点$v \in V$，v的时态介数中心度TBC（v）定义为TBC（v）$= \displaystyle\sum\limits_{s,v,f \in V,s \ne v \ne f} {\dfrac{{{\sigma _{sf}}(v)}}{{{\sigma _{sf}}}}}$.

3.   算法介绍

本节详细阐述基于消息传播的2阶段迭代计算框架以及基于OpenMP框架的多线程并发算法FTBC.

3.1   2阶段迭代计算框架

为了清晰地说明框架的整体思路和2阶段计算过程，首先定义了分裂点集合.

定义5. 分裂点集合. 给定时态图$G = (V,E,T)$，对于任意的顶点$v \in V$，$v$的分裂点集合表示为S（v） = {（v, t_m）|1$\leqslant$m$\leqslant$h}，其中t_m是$v$入边中到达v的时间实例，$h$表示不同的到达时间实例数.

图4（a）给出了一个时态图G示例.G由6个顶点和15条时态边组成. 以顶点b为例：b的分裂点集合$S(b) = \{ (b,0),(b,2),(b,4)\}$.

图  4  时态图示例和FTBC执行中间结果

Figure  4.  An example of temporal graph and intermediate results of FTBC execution

下载: 全尺寸图片 幻灯片

2阶段迭代计算框架中，第1阶段采用自顶向下的广度优先遍历方式计算时态最短路径. 具体地，将每个顶点u作为源点，自顶向下计算源点u到其他顶点和分裂点的最短路径数. 这个阶段需要保存的主要数据结构为：

1）σ_uv和σ_u（v,t）分别记录源点u到v和源点u到分裂点$(v,t)$的时态最短路径数.

2）D_uv和D_u（v,t）分别记录源点u到v和源点u到分裂点（v, t）的时态最短路径长度.

3）flag（v, t）标记（v, t）是否是源点u到v的时态最短路径的终点，如flag（v, t）=1，则表示（v, t）是源点u到v的时态最短路径的终点；否则flag（v, t）=0.

4）P（v, t）记录（v, t）的前驱分裂点集合. 对于时态最短路径$s$$\xrightarrow{{{t_0}}} … \xrightarrow{{{t_k}}}$$u$$\xrightarrow{{{t_{k + 1}}}}$$v$而言，（u, t_k）为（v, t_k+1）的一个前驱分裂点.

第2阶段基于消息传播的机制，采用自底向上的方式迭代计算每个顶点的所有分裂点的时态介数中心度. 这个阶段需要保存的主要数据结构为：

TBC（v）和δ_u（v, t）分别记录v的时态介数中心度和源点u经过分裂点（v, t）的时态最短路径中，（v, t）的所有后续节点对（v, t）的贡献值大小，其中${\delta _{u}}(v,t) =$ $\displaystyle\sum\limits_{u,v,f \in V,u \ne v \ne f} {\dfrac{{{\sigma _{uf}}(v,t)}}{{{\sigma _{uf}}}}}$.

由于时态最短路径不满足子结构特性，因此需要分别计算分裂点的后继节点和孩子节点对其时态介数中心度的贡献值. 具体地，对于时态最短路径$s$$\xrightarrow{{{t_0}}} … \xrightarrow{{{t_k}}}$$\,u \,$$\xrightarrow{{{t_{k + 1}}}}$$v$$\xrightarrow{{{t_{k + 2}}}}$$w$$\xrightarrow{{{t_{k + 3}}}} … \xrightarrow{{{t_n}}}$$z$而言，（v, t_k+1）称为分裂点（u, t_k）的孩子节点；（w, t_k+2）…（z, t_n）称为（u, t_k）的后继节点. 则对于分裂点（u, t_k），需要计算2部分贡献值：第1部分为其后继节点（w, t_k+2）…（z, t_n）对其时态介数中心度的贡献值，需要迭代计算得到；第2部分为其孩子节点（v, t_k+1）对其时态介数中心度的贡献值.2部分贡献值均通过消息传播给分裂点（u, t_k）. 以图4（a）为例：当根据时态最短路径$a$$\xrightarrow{0}$$b$$\xrightarrow{1}$$c$$\xrightarrow{2}$$e$ 更新分裂点（b, 0）的时态介数中心度时，需要考虑其孩子节点（c, 1）和其后继节点（e, 2）的贡献值. 由于flag（c, 1）=0，即（c, 1）不是a到c的时态最短路径的终点，因此（c, 1）对（b, 0）的贡献值为0；只需计算第1部分贡献值，即后继节点（e, 2）对（b, 0）的贡献值.

引理1.　可推导顶点介数中心度公式TBC$(v) =$ $\displaystyle\sum\limits_{s,v,f \in V,s \ne v \ne f} {\dfrac{{{\sigma _{sf}}(v)}}{{{\sigma _{sf}}}}} = \displaystyle\sum\limits_{(v,t) \in S(v)} {{\delta _{s}}(v,t) = } \displaystyle\sum\limits_{(v,t) \in S(v)}\; {\displaystyle\sum\limits_{(w,t'):(v,t) \in P(w,t')} {\left(\dfrac{{{\sigma _{sw}}(v,t)}}{{{\sigma _{sw}}}}\right.}\times} { flag(w,t') +\left. \dfrac{{{\sigma _{sw}}(v,t)}}{{{\sigma _{sw}}(w,t')}} \times {\delta _{s}}(w,t')\right) }$.

证明. 由${\delta _{s}}(v,t) = \displaystyle\sum\limits_{s,v,f \in V,s \ne v \ne f} {\frac{{{\sigma _{sf}}(v,t)}}{{{\sigma _{sf}}}}}$以及定义4可得TBC$(v) = \displaystyle\sum\limits_{s,v,f \in V,s \ne v \ne f} {\frac{{{\sigma _{sf}}(v)}}{{{\sigma _{sf}}}}} = \displaystyle\sum\limits_{(v,t) \in S(v)} {{\delta _{s}}(v,t)}$. 进一步地，可推导出计算公式${\delta _{s}}(v,t) = \displaystyle\sum\limits_{v,f \in V,s \ne v \ne f} {{\delta _{sf}}(v,t)} = \displaystyle\sum\limits_{(w,t'):(v,t) \in P(w,t')} {\displaystyle\sum\limits_{f \in V,s \ne v \ne f} {{\delta _{sf}}((v,t),(v,t) \to (w,t'))} }$：

1）当$f = w$，即$s$$\xrightarrow{{{t_0}}}… \xrightarrow{t}$$v$$\xrightarrow{{t'}}$$w$时，有${\delta _{sf}}((v,t), ((v,t) \to (w,t'))) = \dfrac{{{\sigma _{sw}}(v,t) \times flag(w,t')}}{{{\sigma _{sw}}}}$；

2）当$f \ne w$，即$s$$\xrightarrow{{{t_0}}} …$$\xrightarrow{t}$$v$$\xrightarrow{{t'}}$$w$$\xrightarrow{{{t_k}}} …$$\xrightarrow{{{t_n}}}$$f$时，有${\delta _{sf}}((v,t), ((v,t) \to (w,t'))) = \dfrac{{{\sigma _{sw}}(v,t)}}{{{\sigma _{sw}}(w,t')}} \times \dfrac{{{\sigma _{sf}}(w,t')}}{{{\sigma _{sf}}}} = \dfrac{{{\sigma _{sw}}(v,t)}}{{{\sigma _{sw}}(w,t')}} \times {\delta _{s}}(w,t')$.

因此最终TBC$(v) = \displaystyle\sum\limits_{(v,t) \in S(v)}\; \displaystyle\sum\limits_{(w,t'):(v,t) \in P(w,t')} {\left(\dfrac{{{\sigma _{sw}}(v,t)}}{{{\sigma _{sw}}}} \times \right.} {\left. flag\left(w,t'\right) + \dfrac{{{\sigma _{sw}}(v,t)}}{{{\sigma _{sw}}\left(w,t'\right)}} \times {\delta _{s}}\left(w,t'\right) \right)}$.证毕.

3.2   FTBC算法

基于2阶段迭代计算框架，本节提出了算法1.

算法1.　时态图介数中心度计算算法FTBC.

输入：时态图$G = (V,E,T)$，线程数$\# threadnum$；

输出：所有顶点的介数中心度{TBC$(u),u \in V$}.

① 创建$\# threadnum$个线程，TBC$(u) \leftarrow 0$；

② for each $u \in V$

③ 　线程$i$（${\text{1}} \leqslant i \leqslant \# threadnum$）执行方法

Compute（u,TBC）；

④ end for

⑤ if 所有线程终止

⑥ 　输出TBC$\{ (u),u \in V\}$；

⑦ end if

线程：Compute（u,TBC）；

⑧ 初始化栈S；

⑨ 从u出发BFS方式遍历严格/非严格时态路径；

⑩ 首次访问$(w,t')$则入栈S，计算${\sigma _{u(w,t')}}$，${D_{u(w,t')}}$， ${\sigma _{uw}}$，${D_{uw}}$，$P(w,t')$；

⑪ if $(w,t')$是u到w的时态最短路径的终点

⑫ 　$flag(w,t') \leftarrow 1$；

⑬ end if

⑭ while $(w,t') \leftarrow pop(S)$

⑮ 　if $flag(w,t') = 1$/* 计算孩子节点的贡献值*/

⑯ 　　for each $(v,t) \in P(w,t')$

⑰ 　　　${\delta _{u}}(v,t) \leftarrow {\delta _{u}}(v,t) + \dfrac{{{\sigma _{u(v,t)}}}}{{{\sigma _{uw}}}}$；

⑱ 　　end for

⑲ 　end if

⑳ 　for each $(v,t) \in P(w,t')$/* 计算后继节点的贡 献值*/

㉑ 　　${\delta _{u}}(v,t) \leftarrow {\delta _{u}}(v,t) + \dfrac{{{\sigma _{u(v,t)}}}}{{{\sigma _{u(w,t')}}}} \times {\delta _{u}}(w,t')$；

㉒ 　end for

㉓ end while

㉔ TBC加写锁；

㉕ TBC$(w) \leftarrow$TBC$(w) + {\delta _{u}}(w,t')$；/* 根据引理1累加 计算顶点介数中心度*/

㉖ TBC释放锁.

FTBC算法的输入为时态图G和自定义的线程数，输出为G中所有顶点的时态介数中心度. 首先，FTBC创建线程，初始化TBC数组（行①）. 然后，对于时态图中的每一个顶点u，FTBC委派空闲线程执行Compute方法计算顶点的TBC值（行②~④）. 最后，如果所有线程终止，FTBC返回所有顶点的TBC值（行⑤~⑦）.

Compute方法是一个2阶段迭代计算的过程. 阶段1通过广度优先搜索的方式遍历时态图以完成距离、前驱分裂点和最短路径数的计算，并确定分裂点的flag值（行⑧~⑬）. 具体地，Compute首先初始化栈S（行⑧）. 然后从当前源点u出发遍历严格或非严格时态路径，将首次遍历到的分裂点（w, $t'$）加入S中，并计算源点u到分裂点（w, $t'$）的最短路径数σ_u（w,t'）、源点u到分裂点（w, $t'$）的距离D_u（w,t'）、源点u到顶点w的最短路径数${\sigma _{uw}}$、源点u到w的距离D_uw以及分裂点（w, $t'$）的前驱分裂点集合P（w, $t'$）（行⑨~⑩）. 如果分裂点（w, $t'$）确定是u到w的时态最短路径的终点，则令flag（w, $t'$） = 1（行⑪~⑬）. 阶段2根据引理1自底向上迭代计算分裂点的时态介数中心度，进而累加得到最终顶点的时态介数中心度（行⑭~㉖）. 具体地，当栈S不为空时，从S中弹出分裂点（w, $t'$），如果flag（w, $t'$） = 1，则计算（w, $t'$）对其前驱分裂点（v, t）的贡献值（行⑭~⑲）；接着，累加计算（v, t）经（w, $t'$）到达的所有后继节点对（v, t）的贡献值（行⑳~㉒）. 迭代计算这2部分贡献值直至S为空. 最后，Compute累加计算顶点的时态介数中心度（行㉔~㉖）.

以图4（a）所示时态图为FTBC算法输入，图4（b）给出了源点为a时FTBC算法第1阶段自顶向下计算得到的时态最短路径数σ_uv和σ_u（v,t）、时态最短路径长度D_uv和D_u（v,t）、$flag(v,t)$以及前驱分裂点集合P（v, t）. 以顶点c为例：a到分裂点（c, 5）的时态最短径为$a \xrightarrow{5} c$，则D_a（c,5） = 1，σ_a（c,5） = 1；a到分裂点（c, 1）的时态最短路径为$a$$\xrightarrow{0}$$b$$\xrightarrow{1}$$c$，则D_a（c,1） = 2，σ_a（c,1）= 1；a到分裂点（c, 3）的时态最短路径为$a$$\xrightarrow{{0/2}}$ $b$$\xrightarrow{3}$$c$，D_a（c,3） = 2，σ_a（c,3） = 2.因为D_a（c,5） = 1最小，所以可以得出：flag（c,5）$= 1$；flag（c,1）$= 0$；flag（c,3）$= 0$；D_ac = 1；σ_ac = 1.

图4（c）给出了源点为a时FTBC算法第2阶段自底向上计算得到的贡献值δ_u（v, t）和δ_u（v）. 首先，因为（f, 4）和（f, 6）的flag = 1，因此计算（f, 4）对其前驱分裂点（d, 3）的贡献值以及（f, 6）对其前驱分裂点（d, 3）和（d, 5）的贡献值，明显地，有δ_a（f, 4） = δ_a（f, 6） = 0，则根据引理1可得：δ_a（d, 3）$= \dfrac{2}{7} + \dfrac{2}{7} = \dfrac{4}{7}$；δ_a（d, 5）$= \dfrac{3}{7}$；δ_a（d）$= 1$. 然后，计算flag = 1时的e的分裂点对其前驱分裂点的贡献值，即计算（e, 2）对（c, 1）的贡献值以及（e, 4）对（c, 1）和（c, 3）的贡献值. 明显地，有δ_a（e, 2） = δ_a（e, 4） = 0，则根据引理1可得：δ_a（c,1）$= \dfrac{1}{4} + \dfrac{1}{4} = \dfrac{1}{2}$；δ_a（c,3）$= \dfrac{2}{4} = \dfrac{1}{2}$；δ_a（c）$= 1$. 接着，计算flag = 1的d的分裂点对其前驱分裂点的贡献值，以及经d的分裂点到达的所有后继节点对其前驱分裂点的贡献值，即计算（d, 3），（f, 4），（f, 6）对（b, 0）和（b, 2）以及（d, 5）和（f, 6）对（b, 0），（b, 2），（b, 4）的贡献值. 则：${\delta _{a}}(b,0) = \dfrac{1}{5} + \dfrac{4}{7} \times \dfrac{1}{2} + \dfrac{1}{5} + \dfrac{3}{7} \times \dfrac{1}{3} = \dfrac{{29}}{{35}}$；${\delta _{a}}(b,2) =$ $\dfrac{1}{5} + \dfrac{4}{7} \times \dfrac{1}{2} + \dfrac{1}{5} + \dfrac{3}{7} \times \dfrac{1}{3} = \dfrac{{29}}{{35}}$；${\delta _{a}}(b,4) =$$\dfrac{1}{5} + \dfrac{3}{7} \times \dfrac{1}{3} =$$\dfrac{{12}}{{35}}$. 而后，因为S（c）中分裂点（c, 1）和（c, 3）的flag = 0，所以只需计算经（c, 1）和（c, 3）到达的所有后继节点对（c, 1）和（c, 3）前驱分裂点的贡献值，即计算（e, 2），（e, 4）和（f, 3）对（b, 0）和（e, 4）对（b, 2）的贡献值. 则：${\delta _{a}}(b,0) =$ $\dfrac{{29}}{{35}} + \dfrac{1}{2} \times \dfrac{1}{1} + \dfrac{1}{2} \times \dfrac{1}{2} = \dfrac{{221}}{{140}}$；${\delta _{a}}(b,2) = \dfrac{{29}}{{35}} + \dfrac{1}{2} \times \dfrac{1}{2} = \dfrac{{151}}{{140}}$；δ_a（b）=3.

复杂度分析. 算法1中的主要时间开销为Compute方法的计算开销，Compute方法自顶向下计算时态最短路径花费O（|V|²|T|²），|V|和|T|分别代表时态图顶点数和时间戳集合的大小. 因此FTBC算法的时间复杂度为$O(\dfrac{{{{\left| V \right|}^3}{{\left| T \right|}^2}}}{{\# threadnum}})$. 每个线程保存源点到分裂点的时态最短路径数，需要$O(|V||V|)$空间，边信息需要$O(|E|)$空间，|E|为时态图的边数，FTBC的空间复杂度为$O(\# threadnum|V||T|+|E|)$.

4.   实验分析

本节在真实的数据集中对FTBC算法进行实验测试，并与2种流行算法进行对比，以验证FTBC的效率.

4.1   实验环境和数据集介绍

本文采用了8个真实的数据集进行实验测试. email^[8]数据集是一家中型制造企业员工之间的内部电子邮件通信网络，hypertext^[33]是参会者面对面的接触网络；以hs为前缀的3个数据集^[34-36]是由高中生与朋友构成的联络网络；hospital^[37]，school^[36]，infectious^[33]分别为患者和护工、老师和学生、参展人之间构成的接触网络. 其中email数据集来自KONECT^[38]，其他7个数据集均来自ScocialPattern^[39].表1给出了数据集的统计信息，其中|V|表示顶点数，|E|表示边数，|D|表示时态区间. 时态区间为整个时态图中最大时间戳和最小时间戳的差值.

表  1  数据集信息

Table  1.  Dataset Information

数据集	$\left| V \right|$	$\left| E \right|$	$\left| { {D} } \right|$
hypertext	113	20818	212340
hs2011	126	28560	272330
hs2012	180	45047	729500
hs2013	327	188508	363560
hospital	75	32424	347500
email	167	82927	23430482
school	242	125773	116900
infectious	10972	415912	123837267

下载: 导出CSV 

| 显示表格

本文将FTBC与2个流行算法SBT算法^[31]和SWTBC算法^[32]进行对比. 相关实验代码分别从文献[40-41]中获取. 在实验测试过程中，为取得较好的计算性能，除infectious外的所有数据集默认线程数均设置为8，infectious数据集由于内存限制，线程数设置为1.本文所有实验程序均使用C++语言编写，实验测试环境为一台配置为英特尔至强CPU处理器E5-2640 v4 2.40 GHz，128 GB内存，Linux系统版本为CentOS 7.9的服务器.

4.2   实验结果分析

4.2.1   时态介数中心度分布实验

根据时态介数中心度值，将顶点分散在8个桶中，如果满足$\mathop {\max }\limits_j \dfrac{{T BC(v)}}{{T B{C_{{\text{max}}}}}} \leqslant \dfrac{j}{8},j \in \left\{ {1, 2, …. 8} \right\}$，则顶点$v$被分到第$j$个桶中. 其中TBC_max为时态图中顶点的最大时态介数中心度. 图5（a）~（f）分别展示6个数据集hypertext，hs2012，hs2013，school，hs2011，hospital上的时态介数中心度分布结果. 从图5看出，顶点的介数中心度呈幂律分布，基于非严格和严格时态路径下的顶点介数中心度分布差距不大.

图  5  不同数据集的顶点介数中心度分布

Figure  5.  TBC distribution of different datasets

下载: 全尺寸图片 幻灯片

4.2.2   时间效率对比实验

表2给出了FTBC，SBT，SWTBC算法的时态介数中心度计算时间. 其中，FTBC和SBT算法名称前加上前缀“N”为基于非严格时态路径计算的时态介数中心度，加上前缀“S”为基于严格时态路径计算的时态介数中心度. SWTBC算法仅支持非严格的时态路径. 首先可以看出无论基于严格的还是非严格的时态最短路径计算方式，FTBC算法的计算效率最高. 具体地，FTBC计算时间比SBT快0.7~3倍，比SWTBC算法最多快4个数量级. 这是因为，FTBC算法采用2阶段迭代计算框架，基于引理1，运用并发机制高效迭代计算介数中心度. SBT算法为了降低内存使用，需要不断花费时间清空数据结构中的值；SWTBC算法由于数据集中时间离散程度大，需要设置静态窗口值较大，在复制静态窗口中的图数据方面花费了大量的时间，导致其效率最低.

表  2  顶点介数中心度计算时间

Table  2.  TBC Computation Time s

数据集	NFTBC	SFTBC	NSBT	SSBT	SWTBC
hypertext	3.9	4.1	12.4	12.3	−1
hs2011	2.3	2.3	4.0	4.0	29101.9
hs2012	11.0	11.2	32.1	32.1	−1
hs2013	65.9	68.3	244.3	243.5	−1
hospital	8.0	8.3	17.4	17.4	28612.8
email	218.4	220.7	1001.8	999.7	−1
school	44.1	42.7	178.6	179.1	40872.0
infectious	1280.1	1278.0	1689.4	1646.1	−1
注: −1表示程序12 h内未运行完成；加粗数值表示最快计算速度.

下载: 导出CSV 

| 显示表格

4.2.3   线程数对效率的影响实验

实验验证了线程数对时态介数中心度计算时间的影响. 图6给出了6个数据集分别在线程数为1，8，16，24，32，40进行实验的结果. 从图6可以看出时态介数中心度的计算时间随着线程数的增加呈先减后增的趋势. 这是因为随着线程数的增加，顶点并发计算数增多，计算效率提高. 但当线程数增加到一定数量后，线程开销主导了整体计算开销，线程切换和保证数据一致性的开销增大导致计算效率降低.

图  6  线程数对不同数据集的影响

Figure  6.  Effect of the number of threads on different datasets

下载: 全尺寸图片 幻灯片

5.   总　　结

本文研究了时态图上精确介数中心度计算问题，设计了一种高效的基于消息传播的2阶段迭代计算框架，提出了基于OpenMP框架的多线程并发算法FTBC，通过引理1理论证明了自底向上传播机制的正确性，并通过示例解释了FTBC算法的2阶段迭代计算过程. 基于8个真实的时态图数据集，与2种流行方法进行了时间效率对比，实验验证了FTBC算法的高效性与可扩展性. 通过理论与实验分析可以看出，精确计算时态介数中心度复杂性较高，设计高效的近似时态介数中心度算法是一个重要且值得研究的问题，后续工作计划对其展开深入研究.

作者贡献声明：张天明负责问题定义、方法设计、数据分析与论文撰写和修改工作；赵杰负责方法实现、实验验证、实验结果可视化；金露负责数据收集、实验测试和实验整理；陈璐负责实验分析与结果验证、论文写作指导；曹斌指导论文写作并提出修改建议；范菁指导论文写作并提出修改建议.