Graphlet Degree Vector方法的優化與并行

宋祥帥，楊伏長，謝 江*，張 武，2

（1.上海大學計算機工程與科學學院，上海200444；2.上海大學上海市應用數學與力學研究所，上海200444）

0 引言

比較生物網絡的相似和差異是當前計算生物學的一個主要問題［1］，生物網絡通常由圖來建立模型，圖中節點表示生物分子，如代謝物、蛋白質、基因等，而邊則表示各生物分子之間的相互作用［2］，研究生物網絡可以為疾病的發生機制和治療手段提供深刻的見解［3］。其中，一項很重要的研究就是尋找生物網絡中的自同構軌道。Graphlet Degree Vector（GDV）方法是Przulj在2003年提出的利用圖元及圖元向量來刻畫網絡中節點鄰域關系的方法，具體指在小連通非同構子圖中計算每個節點的自同構軌道，即每個節點所接觸的圖形數量［4］，這種方法基于網絡拓撲和鄰域定義了一系列非同構子圖和圖向量，用于識別網絡中結構相似的模塊［5］。人們利用這種方法進行了許多有意義的研究［6］，例如研究了生物網絡與隨機網絡的拓撲結構差異［7-8］，構建生物網絡的進化樹［9］，識別癌癥相關基因［4］，計算差異網絡的聚類系數［9］，生物網絡進行最優比對［10］和蛋白質功能分析［11］等。然而隨著小連通非同構子圖中節點數的增加，GDV 方法的計算時間復雜度會很高，它的擴展會受到很大的約束［3］。盡管Przulj［8］提出可以利用提高CPU的性能來提高擴展性，但是計算成本會變得越來越高，因此隨著生物網絡研究的規模以及小連通非同構子圖規模的不斷增大，參與枚舉的自同構軌道數量呈指數級別的增長，計算量越來越大，給圖元的擴展帶來了挑戰。

當前圖元方法仍以Przulj 于2003 年提出的GDV 方法為主流［12］，具體實現如Xie 等［5］于2017 年提出的基于2-4 nodes的枚舉方法，通過一個二維矩陣Net_Matrix 來存儲無向生物網絡，然后通過枚舉的方式找出2-4 nodes 連通非同構子圖的15 個自同構軌道，該算法通過枚舉的方式實現了軌道的查找，有效地完成了自同構軌道查找的任務。Ho?evar 等［13］在2014 年提出了一種新的計算網絡節點圖形和軌道特征的組合方法，取得了比較顯著的效果。此外，由于復雜度的原因，Ahmed 等［14］只研究了節點為3 和4 的圖元，利用4 個節點和3個節點的圖元在結構上相似性，減少判斷包含4 個節點的圖元向量的計算開銷。

目前，已有的GDV 方法在計算規模上都存在瓶頸［15］。隨著生物網絡數據獲取的渠道越來越多，生物網絡規模越來越大，對計算效率的要求也會越來越高［15］，因此，實現高效的并行化GDV 方法很有必要。本文從文獻［5］實現的串行的GDV方法著手，將該串行方法以消息傳遞接口（Message Passing Interface，MPI）為基礎實現并行化，并結合去除原來算法的重復運算部分和負載均衡策略改進并行算法，最后，通過仿真數據和真實數據進行了分析和討論。

1 GDV方法的主要思想

GDV 方法的主要思想是計算一個生物網絡中每個節點的自同構軌道數量，即每個節點所接觸的圖形的數量［4］。這在研究生物網絡的過程中發揮著重要的作用。

圖1 展示了包含2、3、4 個節點的非同構圖元。為了刻畫節點的拓撲等價性，Przulj把圖元中具有相同拓撲位置的節點標記為相同的記號，然后對其中具有不同拓撲位置的節點唯一標號。圖1 中包含了2-節點、3-節點、4-節點這三種圖元的15 個不同的拓撲位置，稱這些拓撲位置為自同構軌道，它們出現的頻率記錄為圖元向量［16］。

由于在大規模生物網絡中，非同構圖元的網絡結構差異各種各樣，下面將會以一個簡單無向網絡為例來對自同構軌道的查找進行說明。

圖1 圖元及圖元向量Fig.1 Graphlets and graphlet orbits

圖2 展示了網絡GA 所形成的無向網絡圖。在圖2 中，以節點1 為例，發現一共有4 個0 軌道向量，即（1，2）、（1，3）、（1，4）、（1，5），那么節點1 的0 軌道向量的數目就是4；同樣地，當以節點2 為例時，0 軌道向量的數目是4，即（2，1）、（2，3）、（2，4）、（2，5）。其他軌道向量數目的計算過程依此類推。

表1展示了圖2實例中每個節點的軌道向量數目，其中行代表了軌道向量編號，3 個圖元中軌道向量的總數目為14。列則代表了每個節點的編號。

圖2 無向生物網絡實例Fig.2 Example of undirected biological network

表1 圖2中5個節點在GA網絡中圖元向量的數量Tab.1 Number of graphlet orbits of the five nodes in the GA network of Fig.2

2 GDV方法的實現

GDV 方法是一種在連通生物網絡中枚舉各節點自同構軌道數量的方法，可以大致分為網絡初始化、自同構軌道查找和統計自同構軌道數量三個步驟。其中網絡的初始化是該方法的準備工作，需要將邊集形式轉換為一個無向生物網絡，之后再將網絡轉換成一個n×n（n 指的是無向網絡的節點數目）的鄰接矩陣用以存儲每個節點以及節點所對應的邊；GDV 方法所提出的自同構軌道查找保證了所枚舉圖元的唯一性和查找自同構軌道數量的準確性；最后是統計每個節點所對應的自同構軌道的數目，將其存儲在一個n×15 的矩陣中。查找每個節點的自同構軌道是整個GDV 方法的核心部分，因此，下文在介紹GDV 方法的前提下，將著重介紹每個節點自同構軌道的查找過程。

GDV方法的具體實現步驟（主要步驟如圖3）如下所示：

步驟1 網絡的初始化。GDV 方法在構建網絡時輸入為邊集的形式，節點的編號以連續的數值表示，邊由節點對來確定是否進行生成，若兩個點的節點值在同一個節點對中，那么就生成連接這兩個節點的邊。如圖2 就是由GA 的邊集所構建的無向網絡。然后將生成的網絡轉換成一個n×n 的鄰接矩陣，其中n 表示網絡的節點數。例如將圖2 的GA 無向圖轉化為鄰接矩陣A：

步驟2 查找每個節點的圖元向量的數量。對每個節點進行遍歷，查找其與相鄰節點之間的關系，進而來判定屬于哪一種自同構軌道。通過分析可以得到，在圖1中，G0圖元可以通過一個二維循環來對0 軌道進行查找，從而計算2-節點圖元中0 軌道的數量，但是當計算3-節點或者4-節點的圖元向量時，二維循環難以解決復雜的計算過程。本文采用了Xie等［5］于2017 年實現的基于2-4 nodes 的方法，該方法對不同的圖元向量進行了分類枚舉，巧妙地避開了在二維循環中計算過程復雜的問題，同時也確保了整個查找過程的嚴謹性，做到了精確查找。在該方法中對3-節點的圖元向量進行了三維循環操作，針對4-節點的圖元進行了四維循環操作，有效地解決了復雜的查找過程，但該方法的時間復雜度非常高。

步驟3 將步驟2 查找的結果存儲到一個n×15 的矩陣中，用以記錄每個節點的圖元向量的數量。

圖3 GDV方法的主要步驟Fig.3 Main steps of GDV method

3 GDV方法的優化及其并行化的實現

枚舉每個節點的非同構軌道數量的操作是整個GDV 方法的核心部分，同時也是整個方法中最耗時的部分，因此本文從查找每個節點的非同構軌道數量這一步驟上進行突破，完成了兩方面的工作：1）將串行GDV 方法實現并行化。2）分為兩步對GDV方法進行了優化：①改進GDV串行方法解決鄰接矩陣中重復計算的問題，同時進行并行化；②將改進后的并行化GDV 方法進行優化以解決各進程的負載不均衡問題，實現負載均衡。

3.1 GDV串行方法的并行實現

GDV 方法的主要任務就是尋找一個網絡中每個節點的自同構軌道的數量，由GDV方法步驟2的分析可知，在尋找每個節點的自同構軌道時，根據圖元向量的不同類別進行不同層次的循環查找就可以計算出不同節點的自同構軌道，所以可以將這類問題轉化為矩陣的運算來進行，針對矩陣的運算，本文實現的是按照行數來進行進程間任務的分配，最后再將子任務的結果規約至0號進程。

3.2 GDV串行方法的重復計算問題改善及其并行化實現

盡管進行了GDV 方法的并行化實現，但是該并行方法仍然耗時甚多，為了盡可能地提高該方法的運行效率，首先對GDV 串行方法進行了解決重復計算問題的改進，然后將改進后的方法實現了并行化。

自同構軌道數量的計算中，需要針對無向網絡中每個節點進行遍歷，查找該節點與鄰居節點的關系，進而確定以該節點為目標節點的自同構軌道的數量，查找的過程是在無向網絡轉換為鄰接矩陣后進行的。眾所周知，無向網絡轉換為鄰接矩陣后往往表現為是一個上（下）三角矩陣，以網絡GA 為例，很顯然它的鄰接矩陣A是一個上（下）三角矩陣，因此在計算的過程中只需要針對上（下）三角矩陣進行查找即可，然后再根據對稱關系，在相應的自同構軌道上記錄，最后得到結果。該過程的偽代碼如下所示：

偽代碼中，在進行第二次循環遍歷時，僅需要從第一個位置的下一個元素進行遍歷即可，不需要再從頭進行遍歷，這樣大大地縮短了對比所需要的時間；不過，在遍歷的同時，還需要對鄰接矩陣其對應位置的自同構軌道的數量進行記錄，這樣才能在記錄時避免出現遺漏的現象。由前面的討論可知，GDV 方法最耗時的部分是對每個節點進行枚舉自同構軌道的數量，因此進行并行化處理時需要針對這一問題展開分析，將矩陣按照進程數來進行行分，用以實現并行化。

3.3 改進后GDV方法的并行優化

傳統的矩陣行分較為規則化，但是在本文中由于改進后的GDV 方法中提供的是一個上（下）三角矩陣，使用傳統的方法進行行分時，就會出現每個進程負載極其不均衡的情況，因此在對矩陣進行行分時，需要改進策略，以解決負載不均衡的情況。改進策略屬于一個動態規劃的問題，程序需要根據進程數來進行合理行分。本文所采取的策略如下：

步驟1 根據矩陣的行數和進程的規模數進行劃分，具體劃分規則為：size*=n/(numprocs*2)。

步驟2 將得到的size*按照進程編號的順序分發給各個進程，此時所有進程運算的矩陣的行數為總體行數的一半。

步驟3 將得到的size*按照進程編號的逆序再次分發給各個進程，此時所有進程運算的矩陣的行數為總體行數的另外一半。

步驟4 根據主進程分發的規模，各進程開始進行計算。

步驟5 各進程將所得到的計算結果歸約求和發送給主進程，并行結束。

4 實驗與結果分析

4.1 實驗環境

本次研究使用的實驗平臺是上海大學高性能計算集群“自強4000”。實驗使用4 個內存節點，每個內存節點配置信息如下：2 顆Intel E5-2690 CPU（2.9 GHz/8-core），內存大小為64 GB。集群節點間使用標準的CLOS 二層Infiniband 網絡架構，MPI 庫版本為IntelMPI，實驗運行操作系統為CentOS 6.3，編程語言為C++。

4.2 實驗數據

實驗同時使用了模擬數據和真實生物網絡數據進行性能分析，其中模擬數據使用NetworkX［17］的python 包模擬了三類不同的網絡模型，分別是無標度網絡模型［18］、小世界網絡模型［19］和規則網絡模型［20］。

為了分析改進后的GDV 方法在不同網絡模型中的可拓展性以及泛化能力，實驗中使用了邊數相同（均為4 000）但節點數不同五種網絡模型進行實驗，具體情況如表2所示。

表2 五個模擬網絡數據集Tab.2 Five simulated network datasets

為了分析網絡的邊數對改進后方法的影響，真實生物網絡選取了酵母菌代謝網絡（Yeast Protein Interaction Network，YPIN）和人類基因調控網絡（Human Genetic Regulatory Network，HGRN）兩個生物網絡數據集，其中HGRN 數據來源于STRING（Search Tool for Recurring Instances of Neighbouring Genes）在線數據庫［21］，YPIN數據集來源于Uri Alon實驗室［22］。這兩個網絡的特點是節點數大致相同，但邊數不同：YPIN 的節點數為689，邊數為1 078；HGRN 的節點數為709，邊數為5 560。

圖4（a）為并行的GDV方法在5種網絡中所使用的時間比對。比較小世界模型、隨機模型和無標度模型，三種模型的節點數均為1 000，邊數為4 000，在圖4（a）中可以看出這三種網絡在相同核數下所花費的時間相差不大，因此可以認為在并行的GDV 方法中自同構軌道的查找與網絡的種類是不相關的。通過觀察圖4（a）可以看出，盡管兩種真實網絡的邊數相差很多，但它們的程序運行時間卻相差不多，因此可以認為并行的GDV方法與網絡的邊數并不相關。

圖4 GDV方法的并行性能Fig.4 Parallel performance of GDV method

圖4 （b）是并行的GDV 方法在5 種網絡中的加速比對比曲線。這5 種網絡雖然規模不相同，但它們的加速比曲線幾乎重合，加速比數值幾乎相同，說明了并行的GDV 方法應用范圍廣，在不同的模型中均有好的作用。

在GDV 方法中，查找自同構軌道的計算開銷比較大［12］，而且隨著網絡節點規模的增大，其運行時間消耗得也越來越多，以表3 的兩種生物網絡和表2 中編號3、4、5 的1 000 個節點的網絡為例。從實驗結果圖4（a）中可以看出，當單核運行程序查找自同構軌道數量時，YPIN 和HGRN 網絡查找時間將近1 h，而表2的1 000節點的網絡的運行時間長達近4 h，兩類網絡節點之差僅300 個節點左右。由分析可知，整個查找自同構軌道的運算過程時間復雜度達到了O（n4），因此其運行時間也會隨著網絡規模的增大，呈現出冪函數4 次方級別的增大，因此對于GDV方法的并行化計算是十分有必要的。

4.2.1 解決重復計算問題后的結果

為了驗證解決重復計算后GDV 方法的有效性，在表2 中編號3、4、5 的網絡和YPIN、HGRN 兩個真實網絡下進行多次測試，各種條件下的測試結果都很相似，因此，選取了生成網絡中的一個測試結果進行描述，生成網絡中選取的是編號為3 的無標度網絡，結果如圖5 所示。從圖5 可以明顯地看到，在同一個網絡中，在不同的核數并行情況下，改進后所消耗的時間遠比改進前所消耗的時間少，這在其他的網絡中也有所體現。

圖5 無標度網絡在解決重復計算前后的時間消耗Fig.5 Time consumption before and after solving double counting in scale-free network

但是，解決了重復計算的問題后還面臨著各進程資源分配極其不均勻的問題。為了驗證資源分配不均勻這一問題，選取了表2 中編號3、4、5 的網絡以及YPIN、HGRN 兩個真實網絡模型作為數據集，它們在不同并行核數下程序時間消耗以及加速比如圖6所示。

圖6 解決重復計算后的并行性能Fig.6 Parallel performance after solving double counting

由圖6（a）可以看出，5 個網絡模型在使用一個核和兩個核進行運算時，所消耗的時間相差無幾，而且由圖6（b）可以看出，該并行性能的加速比比較低，原因就是進程之間的資源分配不均勻。因為在多進程的情況下，某一進程所分得到的資源要遠比其他兄弟進程多，從而導致在計算時，該進程所花費的時間遠遠多于兄弟進程，所以造成了這種加速比極低的情況。

從穩定性方面進行分析，通過觀察圖6（b），可以得出五個網絡模型加速比一致，其穩定性良好。

4.2.2 采取負載均衡策略的實驗結果

本節實驗對4.2.1 節提到的各個進程之間資源分配不均衡做出了改進，通過對表2 中編號3、4、5 的網絡和YPIN、HGRN 兩個真實網絡進行了多次測試，與采取負載均衡策略前的結果進行了對比。實驗在5個網絡中分別進行了10次測試，選取了計算所需時間的平均值，并計算出改進前后的加速比，對比結果如圖7所示。

圖7 采取負載均衡策略后的實驗對比Fig.7 Experimental comparison after adopting load balancing strategy

圖7 展示了采取負載均衡策略后加速比的提升，上面的5條曲線展示的是前文提及的5 個網絡模型采取負載均衡策略后的加速比，而下面5條曲線則表示的是5個網絡模型未采取負載均衡策略的加速比，可以看出加速比提升較高，并且隨著并行核數的增大，加速比的提升空間也越來越大。

4.2.3 GDV方法改進后的并行性能

本節主要研究對GDV 方法進行了兩步優化策略后的并行性能。為了說明改進后的GDV 方法與網絡的節點和邊的關系以及該方法是否有良好的拓展性，本次實驗選取了表2中編號1、2、3 的網絡進行測試，該網絡選取的是隨機生成的無標度網絡，結果如圖8所示。

圖8 無標度網絡的并行性能Fig.8 Parallel performance of scale-free networks

在圖8 所顯示的實驗結果中，選取的網絡為無標度網絡，三個網絡的節點分別為500、800、1 000，其對應的網絡的邊數都為4 000。從圖8（a）來看：隨著核數的增加，運行時間變得越來越少，有較好的并行結果；在節點數與核數一定的情況下，從所耗時間角度來看，節點數越多，耗時越久。從圖8（b）來看：隨著核數的增加，加速比也在上升，但是加速比增加的效果并不是很理想，隨著核數的增加，加速比呈現出了非線性上升的狀態；縱向來看，盡管節點數不相同，但是它們的加速比曲線相互疊加，因此可以看出改進后的GDV 方法擁有良好的擴展性。

5 結語

本文實現了GDV 方法的并行計算，同時還提出了一種改進策略應用于GDV 方法：針對原有串行算法在計算自同構軌道時耗時較長的問題，提出了解決重復計算策略和負載均衡策略，大大地節省了程序運行時間并且實現了并行計算的負載均衡。本文使用了多種模擬網絡數據和真實生物網絡數據進行測試，測試結果表明，GDV 的并行方法及改進后的GDV并行方法在多個數據集上都能得到較好的加速比，有效解決了自同構軌道查找效率低的問題。