時態(tài)圖頂點介數中心度計算方法

張?zhí)烀?趙 杰 金 露 陳 璐 曹 斌 范 菁

1（浙江工業(yè)大學計算機科學與技術學院 杭州 310023）

2（浙江大學計算機科學與技術學院 杭州 310013）

介數中心度（betweenness centrality,BC）通過計算經過頂點的最短路徑數來確定頂點在圖結構中的重要程度，是圖頂點重要性計算的一種流行方式，最早由Freeman[1]提出.頂點的介數中心度越大，則頂點對圖中其他頂點的控制能力越強，介數中心度算法被廣泛應用于社會網絡分析[2-3]、蛋白質交互網絡影響預測[4]、社區(qū)發(fā)現[5]等.

現有的介數中心度算法主要聚焦于普通圖，然而現實場景中建模的圖常為時態(tài)圖，邊上帶有時態(tài)信息.例如電子郵件網絡[6]邊上附帶有郵件發(fā)送/接收時間；社交網絡[7]邊上附帶有人與人的接觸時間.圖1 是一個電子郵件網絡示例，頂點表示用戶，邊表示用戶之間發(fā)送/接收郵件的關系，邊上時間戳表示郵件發(fā)送/接收或轉發(fā)/接收的時間.具體地，用戶1在3 個不同時間點向用戶2 共發(fā)送了3 份郵件，而后用戶2 將郵件1 和郵件2 轉發(fā)給用戶3 和用戶4，并將郵件3 轉發(fā)給用戶5.電子郵件網絡上介數中心度計算有助于精確推斷網絡結構[8]、用戶活躍程度等.圖2是一個農貿市場接觸網絡示例，頂點表示供貨商、商戶、消費者；邊表示他們之間的接觸時刻.具體地，消費者1 和消費者2 在商戶1 消費；消費者3 在商戶1和商戶2 消費.蔬菜供貨商為商戶1 送貨，水果供貨商為商戶1 和商戶2 送貨.在傳染病疫情暴發(fā)時，接觸網絡上介數中心度計算有助于識別超級傳播者[9]，控制傳染病的傳播.

Fig.1 An example of email network圖1 電子郵件網絡示例

Fig.2 An example of agricultural market contact network圖2 農貿市場接觸網絡示例

普通圖上介數中心度計算忽略了重要的時態(tài)信息，而時態(tài)信息對于信息流傳播擴散[10]具有重要作用.鑒于此，本文研究的時態(tài)圖頂點介數中心度計算方法與已有的普通圖介數中心度計算方法相比，時態(tài)圖包含時態(tài)信息，且1 對頂點之間存在多條被不同時間戳標記的邊，因此時態(tài)圖上介數中心度計算難度更大.具體原因可分為2 方面：

1）時態(tài)圖頂點介數中心度需要根據時態(tài)最短路徑計算，而時態(tài)最短路徑的定義方法多樣，且計算時需要考慮時態(tài)邊之間的時序依賴關系.例如圖3（a）給出的時態(tài)圖示例中，邊上的值為時間戳信息.如果不考慮時態(tài)信息，則a到d的最短路徑為a→c→d，但實際上a經c到d是不可達的，因為a到c的時間為4（＞3）.所以計算時態(tài)最短路徑時需要考慮時態(tài)邊與邊之間的時序依賴關系，記a到d的時態(tài)最短路徑為a→b→c→d.

Fig.3 Examples of temporal graph and general graph圖3 時態(tài)圖和普通圖示例

2）針對時態(tài)圖，普通圖上介數中心度的計算理論與方法已不再適用，需要設計全新的理論與方法.這是因為普通圖中的介數中心度計算方法主要根據Brandes 算法[11]設計，Brandes 算法有效的關鍵理論是最短路徑的子路徑依然是最短路徑，即最優(yōu)子結構特性.然而時態(tài)最短路徑并不滿足此特性.例如圖3（b）給出的普通圖中，a→c→d是1 條最短路徑，則子路徑a→c一定是最短路徑.而在時態(tài)圖中（如圖3（a）所示），時態(tài)最短路徑a→b→c→d的子路徑a→b→c很明顯不是時態(tài)最短路徑，a到c的時態(tài)最短路徑為a→c.

為了解決這2 個難點，本文根據時態(tài)邊之間的時序依賴關系，定義了嚴格（時態(tài)遞增）和非嚴格（時態(tài)非遞減）2 種時態(tài)路徑類型，提出了基于消息傳播的2 階段迭代計算框架以高效計算時態(tài)圖頂點介數中心度.其中，第1 階段采用自頂向下的廣度優(yōu)先遍歷方式計算時態(tài)最短路徑；第2 階段采用自底向上的方式計算頂點的后繼節(jié)點和孩子節(jié)點對其介數中心度的貢獻值，并設計了基于消息傳播機制的迭代累積計算方法.為了提高效率和可擴展性，實現了基于OpenMP（open multiprocessing）框架的多線程并發(fā)算法FTBC（fast temporal betweenness centrality）.概括而言，本文的主要貢獻有3 點：

1）提出了自頂向下和自底向上結合的2 階段計算框架，并設計了基于消息傳播機制的迭代累積計算方法以高效計算時態(tài)圖頂點介數中心度.

2）提出了基于OpenMP 框架的多線程并發(fā)算法FTBC，以提高時態(tài)圖頂點介數中心度計算的效率和可擴展性，并理論分析了算法的復雜度.

3）基于8 個真實的時態(tài)圖數據集進行了大量的實驗評估，驗證了多線程FTBC 算法相比目前流行的方法計算性能更優(yōu)，可擴展性更強.

1 相關工作

本節(jié)分別概述已有的普通圖和時態(tài)圖介數中心度計算方法.

1.1 普通圖介數中心度計算方法

針對普通圖，Brandes[11]推導了經典的成對依賴、迭代計算理論，并基于此提出了精確介數中心度計算方法.算法空間復雜度為O(n+m)，時間復雜度為O(nm)（無權圖）或O（nm+n2log2n）（有權圖），其中n和m分別表示頂點數和邊數.Erd?s 等人[12]提出了Brandes++算法，該算法采用分治策略，利用圖基礎結構加速計算.Sariyüce 等人[13]提出了BADIOS 算法，該算法針對無向無環(huán)圖，先將一些特殊頂點壓縮，并利用頂點、橋邊將圖劃分為多個子圖，而后在子圖上計算頂點的介數中心度.Baglioni 等人[14]提出利用圖的拓撲特征加速計算.但文獻[12-14]所提算法的最壞時間復雜度仍與Brandes 算法[11]相同.

為了提高介數中心度計算效率，研究人員提出了近似算法.為了避免計算所有頂點對之間的最短路徑，近似算法的總體思想是基于采樣方法計算部分頂點之間的最短路徑，并基于此估算所有頂點的介數中心度.為了保證結果質量，近似算法通常理論推導采樣數或迭代更新介數中心度值，直到滿足預置的終止條件.Brandes 等人[15]提出了基于頂點采樣的方法，其利用霍夫丁不等式[16]估計誤差概率.Riondato 等人[17]提出了基于最短路徑采樣的方法RK，其利用頂點直徑（vertex diameter,VD）和VC（vapnikchervonenkis）維[18]估算達到要求精度所需的最少樣本數.RK 首先采樣一對頂點，然后再基于采樣頂點對的一條最短路徑而不是所有最短路徑來近似計算頂點的介數中心度.此外，由于計算精確的頂點直徑需要所有頂點對之間的最短路徑，此操作非常耗時，因此RK 隨機采樣一個頂點，并根據該頂點到圖中其他頂點的單源最短路徑距離估計頂點直徑.基于RK，Borassi 等人[19]提出了KADABRA 算法，其采用雙向廣度優(yōu)先搜索方式來減少最短路徑的采樣時間，在介數中心度估計時允許為每個頂點設置不同的概率置信度.Riondato 等人[20]提出了ABRA 算法，其利用漸進式隨機抽樣，基于拉德馬赫平均值和偽維估計采樣數.Cousins 等人[21]提出了Bavarain 算法，其采用蒙特卡洛經驗拉德馬赫平均值[22]估計更加嚴格的樣本數上界，進而保證結果精度.Pellegrina 等人[23]提出了SILVAN 算法，與Bavarain 算法相同的是，SILVAN同樣采用蒙特卡洛經驗拉德馬赫平均值估計樣本數上界；與Bavarain 不同的是，Bavarain 僅適用于均勻邊界，而SILVAN 可以用于均勻和非均勻邊界.

然而，無論是普通圖介數中心度精確計算方法還是近似計算方法，方法有效的關鍵理論是最優(yōu)子結構性質，而時態(tài)圖不滿足最優(yōu)子結構特性，因此無法直接擴展到時態(tài)圖.

1.2 時態(tài)圖介數中心度計算方法

一些工作將時態(tài)圖看作是一系列圖的鏡像，即動態(tài)圖，而后研究動態(tài)圖介數中心度計算.Lee 等人[24]將動態(tài)圖分解為連通分量后再進行BC 值更新計算以減少搜索空間.Green 等人[25]提出最短路徑樹結構加速動態(tài)圖BC 值的更新.Kourtellis 等人[26]將最短路徑樹結構壓縮存儲以減少內存占用.Kas 等人[27]拓展了動態(tài)APSP（all pairs of shortest path）算法[28]以實現BC 值動態(tài)更新.Bergamini 等人[29]提出了半動態(tài)的近似計算方法以支持頂點和邊的插入操作.Hayashi 等人[30]提出了全動態(tài)的近似算法以支持邊和頂點的插入和刪除.然而文獻[25-30]算法均將時態(tài)圖視為圖鏡像集合，沒有考慮時態(tài)邊之間的時序依賴關系，本質上還是基于普通圖的最優(yōu)子結構性質設計的方法.

與本文研究問題最相似的工作為Bu?等人[31]提出的時態(tài)圖頂點介數中心度精確計算算法，其首先構建頂點的前置圖，使得前置圖上滿足最優(yōu)子結構性質，而后將Brandes 算法理論擴展，運用在前置圖上進行計算.另一篇較相似的工作為Tsalouchidou 等人[32]提出的算法，其將路徑長度與持續(xù)時間結合起來作為時態(tài)最短路徑定義標準，并在時態(tài)圖定長靜態(tài)窗口上計算介數中心度精確值.本文實驗測試了Bu?等人[31]提出的算法，發(fā)現當時態(tài)圖邊數為3 萬時，其在內存16GB 的機器上已無法運行.Tsalouchidou 等人[32]提出的算法對于時間離散程度較大的時態(tài)圖而言，需要設置較大的靜態(tài)窗口值，導致復制靜態(tài)窗口時耗費了大量的時間.可見，文獻[31-32]工作存在計算效率較低、可擴展性差的問題.

2 問題定義

本節(jié)主要介紹時態(tài)圖、時態(tài)路徑、時態(tài)最短路徑、時態(tài)介數中心度的相關概念，并給出問題定義.

定義1.時態(tài)圖.本文定義的時態(tài)圖既可以是有向的，又可以是無向的，表示為G=(V,E,T).其中，V表示頂點集合；E表示時態(tài)邊集合；T為時間戳集合.時態(tài)圖中2 點之間可以有多條時態(tài)邊.具體地，時態(tài)邊ei=（ui,vi,,ti）表示頂點ui與vi之間的事件發(fā)生時間為ti，ti∈T.

定義2.嚴格與非嚴格時態(tài)路徑.從頂點u到vk的時態(tài)路徑表示為其中，p的第1 條邊e1=（u,v1,t1），第i條邊ei=（vi-1,vi,ti）（1 ＜i＜k），第k條邊ek=（vk-1,vk,tk）使得對于任意的1 ≤i＜k，ti≤ti+1.特別地，當k=1，p也是1 條時態(tài)路徑.進一步地，如果對于任意的1 ≤i＜k，ti＜ti+1，則p稱為嚴格時態(tài)路徑；否則p稱為非嚴格時態(tài)路徑.令|p|=k表示時態(tài)路徑的長度.

定義3.時態(tài)最短路徑.給定從頂點u到v的時態(tài)路徑ps，如果不存在從u到v的其他時態(tài)路徑p滿足路徑長度|p|＜|ps|，則ps是時態(tài)最短路徑.特別地，當|ps|=1 時，ps也是一條時態(tài)最短路徑.

定義4.時態(tài)介數中心度.給定時態(tài)圖G=(V,E,T)，令σsf表示頂點s到頂點f的時態(tài)最短路徑數目；σsf（v）表示由頂點s經頂點v到頂點f的時態(tài)最短路徑數目，則對于所有頂點v∈V，v的時態(tài)介數中心度TBC（v）定義為TBC（v）=.

3 算法介紹

本節(jié)詳細闡述基于消息傳播的2 階段迭代計算框架以及基于OpenMP 框架的多線程并發(fā)算法FTBC.

3.1 2 階段迭代計算框架

為了清晰地說明框架的整體思路和2 階段計算過程，首先定義了分裂點集合.

定義5.分裂點集合.給定時態(tài)圖G=(V,E,T)，對于任意的頂點v∈V，v的分裂點集合表示為S（v）={（v,tm）|1≤m≤h}，其中tm是v入邊中到達v的時間實例，h表示不同的到達時間實例數.

圖4（a）給出了一個時態(tài)圖G示例.G由6 個頂點和15 條時態(tài)邊組成.以頂點b為例：b的分裂點集合S(b)={(b,0),(b,2),(b,4)}.

2 階段迭代計算框架中，第1 階段采用自頂向下的廣度優(yōu)先遍歷方式計算時態(tài)最短路徑.具體地，將每個頂點u作為源點，自頂向下計算源點u到其他頂點和分裂點的最短路徑數.這個階段需要保存的主要數據結構為：

1）σuv和σu（v,t）分別記錄源點u到v和源點u到分裂點(v,t)的時態(tài)最短路徑數.

2）Duv和Du（v,t）分別記錄源點u到v和源點u到分裂點（v,t）的時態(tài)最短路徑長度.

3）flag（v,t）標記（v,t）是否是源點u到v的時態(tài)最短路徑的終點，如flag（v,t）=1，則表示（v,t）是源點u到v的時態(tài)最短路徑的終點；否則flag（v,t）=0.

4）P（v,t）記錄（v,t）的前驅分裂點集合.對于時態(tài)最短路徑而言，（u,tk）為（v,tk+1）的一個前驅分裂點.

第2 階段基于消息傳播的機制，采用自底向上的方式迭代計算每個頂點的所有分裂點的時態(tài)介數中心度.這個階段需要保存的主要數據結構為：

TBC（v）和δu（v,t）分別記錄v的時態(tài)介數中心度和源點u經過分裂點（v,t）的時態(tài)最短路徑中，（v,t）的所有后續(xù)節(jié)點對（v,t）的貢獻值大小，其中δu(v,t)=

由于時態(tài)最短路徑不滿足子結構特性，因此需要分別計算分裂點的后繼節(jié)點和孩子節(jié)點對其時態(tài)介數中心度的貢獻值.具體地，對于時態(tài)最短路徑而言，（v,tk+1）稱為分裂點（u,tk）的孩子節(jié)點；（w,tk+2）…（z,tn）稱為（u,tk）的后繼節(jié)點.則對于分裂點（u,tk），需要計算2 部分貢獻值：第1 部分為其后繼節(jié)點（w,tk+2）…（z,tn）對其時態(tài)介數中心度的貢獻值，需要迭代計算得到；第2 部分為其孩子節(jié)點（v,tk+1）對其時態(tài)介數中心度的貢獻值.2部分貢獻值均通過消息傳播給分裂點（u,tk）.以圖4（a）為例：當根據時態(tài)最短路徑更新分裂點（b,0）的時態(tài)介數中心度時，需要考慮其孩子節(jié)點（c,1）和其后繼節(jié)點（e,2）的貢獻值.由于flag（c,1）=0，即（c,1）不是a到c的時態(tài)最短路徑的終點，因此（c,1）對（b,0）的貢獻值為0；只需計算第1 部分貢獻值，即后繼節(jié)點（e,2）對（b,0）的貢獻值.

3.2 FTBC 算法

基于2 階段迭代計算框架，本節(jié)提出了算法1.

算法1.時態(tài)圖介數中心度計算算法FTBC.

輸入：時態(tài)圖G=(V,E,T)，線程數#threadnum；

輸出：所有頂點的介數中心度{TBC(u),u∈V}.

FTBC 算法的輸入為時態(tài)圖G和自定義的線程數，輸出為G中所有頂點的時態(tài)介數中心度.首先，FTBC 創(chuàng)建線程，初始化TBC數組（行①）.然后，對于時態(tài)圖中的每一個頂點u，FTBC 委派空閑線程執(zhí)行Compute方法計算頂點的TBC值（行②～④）.最后，如果所有線程終止，FTBC 返回所有頂點的TBC值（行⑤～⑦）.

Compute方法是一個2 階段迭代計算的過程.階段1 通過廣度優(yōu)先搜索的方式遍歷時態(tài)圖以完成距離、前驅分裂點和最短路徑數的計算，并確定分裂點的flag值（行⑧～?）.具體地，Compute首先初始化棧S（行⑧）.然后從當前源點u出發(fā)遍歷嚴格或非嚴格時態(tài)路徑，將首次遍歷到的分裂點（w,t′）加入S中，并計算源點u到分裂點（w,t′）的最短路徑數σu（w,t'）、源點u到分裂點（w,t′）的距離Du（w,t'）、源點u到頂點w的最短路徑數σuw、源點u到w的距離Duw以及分裂點（w,t′）的前驅分裂點集合P（w,t′）（行⑨～⑩）.如果分裂點（w,t′）確定是u到w的時態(tài)最短路徑的終點，則令flag（w,t′）=1（行?～?）.階段2 根據引理1自底向上迭代計算分裂點的時態(tài)介數中心度，進而累加得到最終頂點的時態(tài)介數中心度（行?～?）.具體地，當棧S不為空時，從S中彈出分裂點（w,t′），如果flag（w,t′）=1，則計算（w,t′）對其前驅分裂點（v,t）的貢獻值（行?～?）；接著，累加計算（v,t）經（w,t′）到達的所有后繼節(jié)點對（v,t）的貢獻值（行?～?）.迭代計算這2 部分貢獻值直至S為空.最后，Compute累加計算頂點的時態(tài)介數中心度（行?～?）.

4 實驗分析

本節(jié)在真實的數據集中對FTBC 算法進行實驗測試，并與2 種流行算法進行對比，以驗證FTBC 的效率.

4.1 實驗環(huán)境和數據集介紹

本文采用了8 個真實的數據集進行實驗測試.email[8]數據集是一家中型制造企業(yè)員工之間的內部電子郵件通信網絡，hypertext[33]是參會者面對面的接觸網絡；以hs 為前綴的3 個數據集[34-36]是由高中生與朋友構成的聯絡網絡；hospital[37]，school[36]，infectious[33]分別為患者和護工、老師和學生、參展人之間構成的接觸網絡.其中email 數據集來自KONECT[38]，其他7 個數據集均來自ScocialPattern[39].表1 給出了數據集的統(tǒng)計信息，其中|V|表示頂點數，|E|表示邊數，|D|表示時態(tài)區(qū)間.時態(tài)區(qū)間為整個時態(tài)圖中最大時間戳和最小時間戳的差值.

Table 1 Dataset Information表1 數據集信息

本文將FTBC 與2 個流行算法SBT 算法[31]和SWTBC 算法[32]進行對比.相關實驗代碼分別從文獻[40-41]中獲取.在實驗測試過程中，為取得較好的計算性能，除infectious 外的所有數據集默認線程數均設置為8，infectious 數據集由于內存限制，線程數設置為1.本文所有實驗程序均使用C++語言編寫，實驗測試環(huán)境為一臺配置為英特爾至強CPU 處理器E5-2640 v4 2.40 GHz，128 GB 內存，Linux 系統(tǒng)版本為CentOS 7.9 的服務器.

4.2 實驗結果分析

4.2.1 時態(tài)介數中心度分布實驗

Fig.5 TBC distribution of different datasets圖5 不同數據集的頂點介數中心度分布

4.2.2 時間效率對比實驗

表2 給出了FTBC，SBT，SWTBC 算法的時態(tài)介數中心度計算時間.其中，FTBC 和SBT 算法名稱前加上前綴“N”為基于非嚴格時態(tài)路徑計算的時態(tài)介數中心度，加上前綴“S”為基于嚴格時態(tài)路徑計算的時態(tài)介數中心度.SWTBC 算法僅支持非嚴格的時態(tài)路徑.首先可以看出無論基于嚴格的還是非嚴格的時態(tài)最短路徑計算方式，FTBC 算法的計算效率最高.具體地，FTBC 計算時間比SBT 快0.7～3 倍，比SWTBC 算法最多快4 個數量級.這是因為，FTBC 算法采用2 階段迭代計算框架，基于引理1，運用并發(fā)機制高效迭代計算介數中心度.SBT 算法為了降低內存使用，需要不斷花費時間清空數據結構中的值；SWTBC 算法由于數據集中時間離散程度大，需要設置靜態(tài)窗口值較大，在復制靜態(tài)窗口中的圖數據方面花費了大量的時間，導致其效率最低.

Table 2 TBC Computation Time表2 頂點介數中心度計算時間 s

4.2.3 線程數對效率的影響實驗

實驗驗證了線程數對時態(tài)介數中心度計算時間的影響.圖6 給出了6 個數據集分別在線程數為1，8，16，24，32，40 進行實驗的結果.從圖6 可以看出時態(tài)介數中心度的計算時間隨著線程數的增加呈先減后增的趨勢.這是因為隨著線程數的增加，頂點并發(fā)計算數增多，計算效率提高.但當線程數增加到一定數量后，線程開銷主導了整體計算開銷，線程切換和保證數據一致性的開銷增大導致計算效率降低.

Fig.6 Effect of the number of threads on different datasets圖6 線程數對不同數據集的影響

5 總結

本文研究了時態(tài)圖上精確介數中心度計算問題，設計了一種高效的基于消息傳播的2 階段迭代計算框架，提出了基于OpenMP 框架的多線程并發(fā)算法FTBC，通過引理1 理論證明了自底向上傳播機制的正確性，并通過示例解釋了FTBC 算法的2 階段迭代計算過程.基于8 個真實的時態(tài)圖數據集，與2 種流行方法進行了時間效率對比，實驗驗證了FTBC 算法的高效性與可擴展性.通過理論與實驗分析可以看出，精確計算時態(tài)介數中心度復雜性較高，設計高效的近似時態(tài)介數中心度算法是一個重要且值得研究的問題，后續(xù)工作計劃對其展開深入研究.

作者貢獻聲明：張?zhí)烀髫撠焼栴}定義、方法設計、數據分析與論文撰寫和修改工作；趙杰負責方法實現、實驗驗證、實驗結果可視化；金露負責數據收集、實驗測試和實驗整理；陳璐負責實驗分析與結果驗證、論文寫作指導；曹斌指導論文寫作并提出修改建議；范菁指導論文寫作并提出修改建議.