基于GRU網(wǎng)絡(luò)的格蘭杰因果網(wǎng)絡(luò)重構(gòu)

楊官學(xué)，王家棟

（江蘇大學(xué) 電氣信息工程學(xué)院，江蘇 鎮(zhèn)江 212013）

0 引言

現(xiàn)實(shí)生活中，許多復(fù)雜系統(tǒng)均可在網(wǎng)絡(luò)角度被抽象表達(dá)，其中網(wǎng)絡(luò)節(jié)點(diǎn)代表系統(tǒng)變量，連邊代表各變量間的相互作用關(guān)系。基于關(guān)聯(lián)的屬性和特征，網(wǎng)絡(luò)可被劃分為不同種類，例如根據(jù)邊是否有權(quán)重或方向，網(wǎng)絡(luò)可被分為有權(quán)網(wǎng)絡(luò)和無權(quán)網(wǎng)絡(luò)、有向網(wǎng)絡(luò)和無向網(wǎng)絡(luò)。典型的網(wǎng)絡(luò)包含電力網(wǎng)絡(luò)、腦網(wǎng)絡(luò)、生物網(wǎng)絡(luò)、交通網(wǎng)、人際網(wǎng)［1-2］等。然而，實(shí)際的網(wǎng)絡(luò)結(jié)構(gòu)往往是未知或難以直接觀察。在當(dāng)今大數(shù)據(jù)時代，復(fù)雜網(wǎng)絡(luò)涵蓋的信息量越來越多，如何基于隱藏在節(jié)點(diǎn)背后的數(shù)據(jù)挖掘節(jié)點(diǎn)間的作用關(guān)系，就顯得至關(guān)重要，這也使得復(fù)雜網(wǎng)絡(luò)重構(gòu)成為當(dāng)前研究的熱點(diǎn)方向之一［3-4］。復(fù)雜網(wǎng)絡(luò)重構(gòu)不僅能使人們更好地了解系統(tǒng)的動力學(xué)行為和演化機(jī)制，也是后續(xù)網(wǎng)絡(luò)結(jié)構(gòu)研究和分析的基礎(chǔ)。例如，針對腦神經(jīng)網(wǎng)絡(luò)研究大腦神經(jīng)元間如何相互激活［5］；針對基因調(diào)控網(wǎng)絡(luò)從基因表達(dá)水平的時間序列推斷基因調(diào)控網(wǎng)絡(luò)［6］等。

當(dāng)下，復(fù)雜網(wǎng)絡(luò)重構(gòu)算法包括相關(guān)系數(shù)［7-8］、互信息［9-11］、布爾網(wǎng)絡(luò)［12-13］、貝葉斯網(wǎng)絡(luò)［14］、格蘭杰因果［15］、壓縮感知［16］等。尤其在有向網(wǎng)絡(luò)的推斷方面，格蘭杰因果（Granger Causality，GC）及其拓展形式得到了廣泛應(yīng)用。傳統(tǒng)格蘭杰因果是一種基于線性矢量自回歸模型的判定方法，能判斷兩個時間序列變量間的因果關(guān)系，但當(dāng)網(wǎng)絡(luò)中存在兩個以上變量時可能會導(dǎo)致虛假連邊。為此，提出條件格蘭杰因果來解決多個變量間的間接影響［17］。

此外，針對樣本量不足的問題，在考慮實(shí)際網(wǎng)絡(luò)稀疏性的前提下，很多學(xué)者將稀疏特征引入格蘭杰因果中，提出LASSO 格蘭杰因果［18］、組稀疏格蘭杰因果［19］、基于貪婪算法的格蘭杰因果［20］等方法，這些方法雖然結(jié)構(gòu)形式簡單、運(yùn)算方便，但無法匹配非線性因果關(guān)系的處理情形。因此，很多學(xué)者將核函數(shù)引入傳統(tǒng)格蘭杰因果中提出核格蘭杰因果方法［21-22］，但本質(zhì)上還只是針對某種特定的非線性情形，適用范圍存在一定限制。

考慮到神經(jīng)網(wǎng)絡(luò)在非線性建模方面能非常擅長地處理輸入和輸出間復(fù)雜的非線性關(guān)系。在因果預(yù)測方面結(jié)合格蘭杰因果模型、多層感知機(jī)（Multi-Layer Perceptron，MLP）、循環(huán)神經(jīng)網(wǎng)絡(luò)（Recurrent Neural Network，RNN）及長短期記憶元（Long Short-Term Memory，LSTM）已得到了應(yīng)用。Chivukula 等［23］使用RNN 分析雅虎財(cái)經(jīng)獲取的真實(shí)股市數(shù)據(jù)，結(jié)果表明因果特征顯著改善了現(xiàn)有的深度學(xué)習(xí)回歸模型。Pathod 等［24］使用RNN 和LSTM 處理多元腦連接性檢測問題提出了RNN-GC 模型，能對非線性和變長時延信息傳輸進(jìn)行建模，并在方向性腦連接性估計(jì)方面十分有效。Tank 等［25］使用MLP 和LSTM 在時序數(shù)據(jù)上進(jìn)行建模，并融合LASSO 提取節(jié)點(diǎn)間的因果關(guān)系，兩種方法均能自動探索最大滯后階數(shù)，且在復(fù)雜非線性DREAM3 數(shù)據(jù)上取得了較好的結(jié)果。

然而，MLP、RNN 存在收斂慢、容易過擬合、可解釋性差、計(jì)算復(fù)雜度高等缺點(diǎn)。其中，RNN 的循環(huán)結(jié)構(gòu)使其可能無法很好地處理長期依賴關(guān)系，易出現(xiàn)梯度消失或梯度爆炸現(xiàn)象；LSTM 在RNN 上進(jìn)行改進(jìn)，不易出現(xiàn)梯度問題，但在處理噪聲較大的數(shù)據(jù)時性能較差。相較于LSTM 網(wǎng)絡(luò)，GRU 網(wǎng)絡(luò)結(jié)構(gòu)僅有更新門和重置門，更容易訓(xùn)練。因此，考慮到變量間影響關(guān)系的非線性和因果性，本文提出一種基于門控循環(huán)單元（Gated Recurrent Unit，GRU）網(wǎng)絡(luò)的格蘭杰因果網(wǎng)絡(luò)重構(gòu)方法（GRUGC）。首先圍繞每個目標(biāo)節(jié)點(diǎn)，構(gòu)建基于GRU 網(wǎng)絡(luò)的格蘭杰因果模型；然后對網(wǎng)絡(luò)輸入權(quán)重進(jìn)行組稀疏懲罰約束；最后基于Adam 的梯度下降網(wǎng)絡(luò)訓(xùn)練法，獲取節(jié)點(diǎn)之間的格蘭杰因果關(guān)系。

在仿真驗(yàn)證方面，首先基于模型生成的數(shù)據(jù)集進(jìn)行相關(guān)仿真研究，例如線性矢量自回歸、非線性矢量自回歸、非均勻嵌入時滯矢量自回歸、Lorenz-96 模型。然后，采用經(jīng)典的DREAM3 競賽數(shù)據(jù)集（Ecoli 數(shù)據(jù)集和Yeast 數(shù)據(jù)集）分析模型性能。

1 格蘭杰因果模型

在格蘭杰因果分析中，傳統(tǒng)分析方法通常采用線性矢量自回歸模型（Vector-Autoregression，VAR），設(shè)模型的最大時滯階數(shù)為P。

式中：Xt為模型中所有變量在t時刻的樣本矩陣。可表示為：

式中：N為模型變量個數(shù)；M為樣本數(shù)目；xi，t∈R1×M為第i個變量在t時刻的樣本向量為 第i個 變量 在t時刻的第m個樣本值；i=1，2，…，N；m=1，2，…，M；p為模型階數(shù)為t-p時刻Xt-p對應(yīng)的系數(shù)矩陣。可表示為：

式中：p=1，2，…，P表示在t-p時 刻第i個 變量和第j個變量間的影響關(guān)系；Xt-p∈RN×M為模型變量在t-p時刻的樣本矩陣，形式如Xt；et∈RN×M為服從標(biāo)準(zhǔn)正態(tài)分布的噪聲。

在VAR 中，當(dāng)且僅當(dāng)對所有階數(shù)p有=0，即可獲得變量j不是影響變量i的格蘭杰原因；反之，變量j是影響變量i的格蘭杰原因，即存在一條從變量i指向變量j的有向邊。在高維情況下考慮到稀疏性，格蘭杰因果模型可被認(rèn)為是一個帶有組稀疏正則化的回歸問題。

式中：‖?‖2為?2正則化，目的是防止模型發(fā)生過擬合現(xiàn)象；λ為懲罰系數(shù)。

2 神經(jīng)網(wǎng)絡(luò)格蘭杰因果模型

為了解決傳統(tǒng)格蘭杰因果模型無法有效處理非線性關(guān)系的問題，提出了一種基于神經(jīng)網(wǎng)絡(luò)的格蘭杰因果模型。首先構(gòu)建一個通用的非線性矢量自回歸模型，假設(shè)該模型中變量的個數(shù)為N，最大時滯階數(shù)為P，則模型t時刻的輸出yt的通用表達(dá)式如式（5）所示。

式中：f(?)為輸入輸出間的非線性映射函數(shù)關(guān)系為在t之前時刻模型各輸入時滯分量的集合，即=1 ≤p≤P；et為服從標(biāo)準(zhǔn)正態(tài)分布的高斯白噪聲。

目前，通常使用神經(jīng)網(wǎng)絡(luò)對非線性函數(shù)f(?)進(jìn)行預(yù)測，常見方法包括MLP、RNN、LSTM 等。由于神經(jīng)網(wǎng)絡(luò)為黑盒模型，且模型中各輸入共享隱藏層，難以直接進(jìn)行因果推斷的分析研究。因此，在式（5）中并未直接采用常規(guī)的多輸入多輸出映射來獲取時序變量間的關(guān)系。

為了進(jìn)一步詳細(xì)闡述模型的原理，將整個網(wǎng)絡(luò)重構(gòu)任務(wù)分解為每個目標(biāo)節(jié)點(diǎn)的鄰居節(jié)點(diǎn)選擇問題。針對每個目標(biāo)節(jié)點(diǎn)i分別使用單獨(dú)的模型fi(?)，以清晰地提取與鄰居節(jié)點(diǎn)間的格蘭杰因果關(guān)系。節(jié)點(diǎn)i在t時刻的表達(dá)值xi，t由式（6）所示。

式中：fi(?)采用的是GRU 網(wǎng)絡(luò)。

同理，當(dāng)且僅當(dāng)所有階數(shù)p(1 ≤p≤P)在有關(guān)xi，t的預(yù)測時不依賴節(jié)點(diǎn)變量j的時滯分量，即在預(yù)測模型中加入{xj，t-1，xj，t-2，…，xj，t-P}后并不能提升對xi，t的預(yù)測精度，則認(rèn)為節(jié)點(diǎn)j不是影響節(jié)點(diǎn)i的原因；反之，節(jié)點(diǎn)j是影響節(jié)點(diǎn)i的原因，即存在j→i。

2.1 基于GRU的循環(huán)神經(jīng)網(wǎng)絡(luò)

RNN 屬于一種具有短期記憶能力的神經(jīng)網(wǎng)絡(luò)，尤其適用于處理和預(yù)測時間序列數(shù)據(jù)，作為一種常用的深度學(xué)習(xí)環(huán)路網(wǎng)絡(luò)結(jié)構(gòu)模型，RNN 由一個或多個循環(huán)層組成，每個循環(huán)層包含多個神經(jīng)元，不僅能接收其他神經(jīng)元信息，還能接收自身信息。理論上，RNN 可逼近任意的復(fù)雜非線性動力系統(tǒng)。在RNN 網(wǎng)絡(luò)參數(shù)訓(xùn)練過程中，由于鏈?zhǔn)椒▌t，隨著錯誤信息的反向傳播，當(dāng)輸入序列時間步長較長時，梯度值可能會趨近于0（梯度消失）或非常大（梯度爆炸）。這兩種情況均會導(dǎo)致訓(xùn)練效果不佳。為了解決RNN 的長程依賴問題，在本文提出神經(jīng)網(wǎng)絡(luò)格蘭杰因果模型的框架中，采用了基于門控循環(huán)單元（Gated Recurrent Unit，GRU）的循環(huán)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)模型，如圖1所示。

Fig.1 Gated loop unit structure of GRU networks圖1 GRU網(wǎng)絡(luò)的門控循環(huán)單元結(jié)構(gòu)

GRU 內(nèi)部結(jié)構(gòu)主要由兩個門構(gòu)成，分別為重置門rt∈[0，1]D（reset gate）和更新門zt∈[0，1]D（update gate），通過調(diào)節(jié)這兩個門的激活值來控制信息流動，解決梯度消失問題，以能更好地捕捉序列中的長期依賴關(guān)系。GRU 單元的內(nèi)部計(jì)算式如式（7）—式（10）所示：

式中：xt∈RN×1為t時刻網(wǎng)絡(luò)輸入；ht∈RD×1為t時刻隱含層的狀態(tài)；W*∈RD×N、U*∈RD×D、b*∈RD×1為待學(xué)習(xí)的神經(jīng)網(wǎng)絡(luò)參數(shù)；*屬于集合{r，z，h}中的元素；⊙為Hadamard Product，表示矩陣對應(yīng)位置元素的乘積。

重置門rt和更新門zt狀態(tài)的迭代更新來源于上一時刻的隱藏層ht-1和當(dāng)前時刻的輸入節(jié)點(diǎn)xt，待獲得ht-1、xt后分別通過式（7）、式（8）計(jì)算重置門rt和更新門zt的值。式（9）首先使用rt⊙ht-1更新數(shù)據(jù)得到；然后與各自權(quán)重相乘后將Whxt與相加完成特征融合；最后使用tanh函數(shù)將特征值映射到［-1，1］，得到中具有上一時刻的隱含層狀態(tài)ht-1和當(dāng)前輸入狀態(tài)xt的信息。

式（10）為GRU 最重要的一步，在這個階段遺忘和記憶同時進(jìn)行，既忘記ht-1中部分信息，又記住中節(jié)點(diǎn)輸入的部分信息。GRU 網(wǎng)絡(luò)直接用一個門來控制輸入和遺忘間的平衡關(guān)系，當(dāng)zt越接近1 表示記住的數(shù)據(jù)越多，反之遺忘的越多。

針對時間序列輸入xt，將xt在不同時刻輸入到循環(huán)神經(jīng)網(wǎng)絡(luò)中得到相應(yīng)的隱含層狀態(tài)ht，循環(huán)神經(jīng)網(wǎng)絡(luò)的輸出yt最終由ht的加權(quán)線性求和表示，如式（11）所示。

式中：Wo∈RD×1為輸出層權(quán)重；Wo為待學(xué)習(xí)的網(wǎng)絡(luò)參數(shù)。

2.2 基于GRU網(wǎng)絡(luò)的格蘭杰因果

在因果網(wǎng)絡(luò)重構(gòu)問題中，根據(jù)式（6）、式（11）針對每個目標(biāo)節(jié)點(diǎn)i構(gòu)建子任務(wù)優(yōu)化模型，分別得出相應(yīng)的鄰居節(jié)點(diǎn)。具體表達(dá)形式如下：

針對輸入節(jié)點(diǎn)變量xt，將輸入層至隱含層的權(quán)重矩陣Wr、Wz、Wh進(jìn)行堆疊拼接以便于后續(xù)優(yōu)化，由M=表示總輸入權(quán)重矩陣。為了捕捉驅(qū)動節(jié)點(diǎn)j對目標(biāo)節(jié)點(diǎn)i的因果影響，當(dāng)且僅當(dāng)總輸入權(quán)重矩陣M的第j列均為0，即M：，j=0得出驅(qū)動節(jié)點(diǎn)j對目標(biāo)節(jié)點(diǎn)i不存在因果關(guān)系；反之驅(qū)動節(jié)點(diǎn)j是影響目標(biāo)節(jié)點(diǎn)i的原因，即存在j→i的連邊。此外，合理加入正則項(xiàng)能提升神經(jīng)網(wǎng)絡(luò)模型的泛化能力，?2正則化項(xiàng)的懲罰因子λ同時也是控制網(wǎng)絡(luò)結(jié)構(gòu)稀疏性的參數(shù)［26］。

2.3 基于Adam的梯度下降法

當(dāng)下，最流行的神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)優(yōu)化方法大都始于隨機(jī)梯度下降法（Stochastic Gradient Descent，SGD）［27］，但SGD每次只用一個樣本更新梯度，使得SGD 并非每次迭代均向最優(yōu)方向更新參數(shù)，容易造成準(zhǔn)確度下降且無法線性收斂的情況。同時，由于只用一個樣本更新梯度并不能代表全部樣本的趨勢，因此易陷入局部最小值。

為此，本文采用基于以往梯度信息和動量的Adam［28］優(yōu)化方法。該方法是RMSProp 和動量法的結(jié)合，主要糾正兩項(xiàng)偏差和平均梯度滑動的方法，具體計(jì)算過程如下：

步驟1：初始化學(xué)習(xí)率lr、平滑常數(shù)β1、β2（分別用于平滑mt、vt），可學(xué)習(xí)參數(shù)θ0=0、m0=0、v0=0、t=0。

步驟2：當(dāng)未停止訓(xùn)練時，更新訓(xùn)練次數(shù)t=t+1，計(jì)算梯度gt（所有的可學(xué)習(xí)參數(shù)都有各自梯度，因此gt指全部梯度的集合）。

本文神經(jīng)網(wǎng)絡(luò)模型的訓(xùn)練方法采用了以Adam 為基礎(chǔ)的優(yōu)化方法。為了研究正則項(xiàng)的影響，訓(xùn)練采用了無正則項(xiàng)的AdamU 和有正則項(xiàng)的Adam。

3 仿真實(shí)驗(yàn)

為了驗(yàn)證本文所提基于GRU 網(wǎng)絡(luò)的格蘭杰因果網(wǎng)絡(luò)重構(gòu)算法的有效性，分別基于線性VAR、非線性VAR、Lorenz-96、非均勻嵌入時滯VAR 模型及DREAM 競賽數(shù)據(jù)集進(jìn)行仿真實(shí)驗(yàn)研究。網(wǎng)絡(luò)重構(gòu)性能的具體量化性能指標(biāo)采用ROC 曲線下的面積（Area Under Receiver-Operating-Characteristic Curve，AUROC）和PR 曲線下的面積（Area Under Precision-Recall Curve，AUPR）。

3.1 線性VAR

首先基于N=10 的VAR 網(wǎng)絡(luò)模型，利用式（1）隨機(jī)生成的稀疏轉(zhuǎn)移矩陣生成仿真時間步長T=1 000 的時間序列數(shù)據(jù)。噪聲服從高斯分布N(0，σ2)，σ=0.1，基準(zhǔn)因果網(wǎng)絡(luò)如圖2 所示。由此可見，若節(jié)點(diǎn)與節(jié)點(diǎn)間存在GC 關(guān)系，則對應(yīng)區(qū)間的數(shù)值大小為1；若不存在GC 關(guān)系，則對應(yīng)區(qū)間的數(shù)值大小為0。例如，圖2 中（1，1）位置和（1，9）位置的數(shù)值大小均為1，分別表示存在一條由節(jié)點(diǎn)1 指向節(jié)點(diǎn)1的GC 邊和存在一條由節(jié)點(diǎn)1指向節(jié)點(diǎn)9的GC 邊。

Fig.2 Benchmark network of VAR model圖2 VAR模型基準(zhǔn)網(wǎng)絡(luò)

具體參數(shù)設(shè)置如下：針對每一個目標(biāo)節(jié)點(diǎn)i，創(chuàng)建自我依賴關(guān)系（自環(huán)），并在其他N-1 個節(jié)點(diǎn)中隨機(jī)選擇一個驅(qū)動節(jié)點(diǎn)構(gòu)建因果關(guān)系，即在j→i的情況下設(shè)置0.1，p=1，2，3，同時令其他=0。GRU 網(wǎng)絡(luò)模型的隱藏層單元數(shù)D設(shè)定為100，分別比較Adam、AdamU 兩種訓(xùn)練優(yōu)化方法，最大滯后階數(shù)P=5，正則系數(shù)λ=0.002，學(xué)習(xí)率lr=0.05，訓(xùn)練步長為20 000，AdamU 無需設(shè)置正則系數(shù)λ。

VAR 模型的兩種訓(xùn)練方法仿真結(jié)果如圖3 所示，圖中區(qū)域數(shù)值大小表示經(jīng)過GRUGC 推理后存在GC 邊的概率值。與圖2 的基準(zhǔn)網(wǎng)絡(luò)比較可見，考慮正則化的Adam 相較于考慮正則化的AdamU 效果更好。

Fig.3 Results of two training methods for VAR model圖3 VAR模型的兩種訓(xùn)練方法結(jié)果

為了量化每種訓(xùn)練方法的性能，計(jì)算圖3 不同方法的AUROC 和AUPR，具體數(shù)值如表1 所示。由此可知，Adam的AUROC 和AUPR 均達(dá)到1，實(shí)現(xiàn)了完美的因果網(wǎng)絡(luò)重構(gòu)，而AdamU 的AUROC、AUPR 僅為0.225、0.125。

Table 1 Performance comparison of two training methods for VAR model表1 VAR模型的兩種訓(xùn)練方法性能比較

3.2 Lorenz-96

設(shè)置網(wǎng)絡(luò)節(jié)點(diǎn)個數(shù)N=10，F(xiàn)=10，添加服從高斯分布的噪聲N(0，0.012)，獲得序列長度T=1 000 的多變量時間序列矩陣。Lorenz-96 模型的微分方程表達(dá)式（14）所示，可得N=10的Lorenz-96模型的基準(zhǔn)網(wǎng)絡(luò)如圖4所示。

Fig.4 Benchmark network of Lorenz-96 model圖4 Lorenz-96模型的基準(zhǔn)網(wǎng)絡(luò)

式中：i=1，2，…，N；F為懲罰力度，值越大表示時間序列的非線性越強(qiáng)。

設(shè)置GRU 網(wǎng)絡(luò)模型的隱藏層單元數(shù)D=100，比較Adam 和AdamU 兩種訓(xùn)練優(yōu)化方法，最大滯后階數(shù)P=10，正則系數(shù)λ=0.2，學(xué)習(xí)率lr=0.001，訓(xùn)練步長為10 000。Lorenz-96 模型的兩種訓(xùn)練方法結(jié)果及相關(guān)性能指標(biāo)如圖5、表2 所示。由此可見，經(jīng)AdamU 訓(xùn)練獲取的網(wǎng)絡(luò)結(jié)構(gòu)較差，Adam 獲取的網(wǎng)絡(luò)結(jié)構(gòu)幾乎全部預(yù)測正確；表2 中AdamU 方法的AUROC 和AUPR 分別為0.825、0.744，Adam 的AUROC 和AUPR 均等于1，達(dá)到了完美重構(gòu)效果。

Table 2 Performance comparison of two training methods for Lorenz-96 model表2 Lorenz-96模型的兩種訓(xùn)練方法性能比較

Fig.5 Results of two training methods for Lorenz-96 model圖5 Lorenz-96模型的兩種訓(xùn)練方法結(jié)果

通過線性VAR 模型和Lorenz-96 模型的仿真發(fā)現(xiàn)，有無正則項(xiàng)對模型最終的訓(xùn)練結(jié)果影響較大，因此在仿真后續(xù)部分將僅基于Adam 方法進(jìn)行訓(xùn)練。

此外，為了考察不同非線性程度F和隱藏層單元數(shù)D對模型的影響，基于Adam 方法對不同F(xiàn)和D數(shù)值進(jìn)行消融實(shí)驗(yàn)，比較結(jié)果如表3、表4 所示。由此可知，在不同D的情況下隨著F增加，AUROC、AUPR 均會下降；在F=10時隨著D增加AUROC 和AUPR 變化不大；當(dāng)F={20，30，40}時，隨著D增加AUROC、AUPR 均得到了不同程度的提升。

Table 3 AUROC for different F and D表3 不同F(xiàn)與D情況下的AUROC

Table 4 AUPR for different F and D表4 不同F(xiàn)與D情況下的AUPR

總之，隨著數(shù)據(jù)的非線性增強(qiáng)，GRU 網(wǎng)絡(luò)的擬合能力變?nèi)酢Ｔ诜蔷€性程度較低時，增加隱藏層單元數(shù)無法大幅度提升模型精度；在非線性程度較高時，增加隱藏層單元數(shù)能明顯提升模型精度，尤其在D從5 增加到10、從10 增加到25 時，AUROC 和AUPR 提升最明顯，這也符合實(shí)際經(jīng)驗(yàn)，在數(shù)據(jù)較復(fù)雜時模型精度往往會降低，可通過采用更復(fù)雜的模型提升表達(dá)能力。并且，在神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)中增加網(wǎng)絡(luò)深度和寬度是提高精度的常見方法［29］。

3.3 非線性VAR

一種非線性VAR 模型如式（15）所示，生成T=1 000的多變量時間序列矩陣，非線性VAR 模型基準(zhǔn)網(wǎng)絡(luò)如圖6所示。

Fig.6 Benchmark network of nonlinear VAR model圖6 非線性VAR模型的基準(zhǔn)網(wǎng)絡(luò)

網(wǎng)絡(luò)訓(xùn)練采用Adam 優(yōu)化方法，GRU 網(wǎng)絡(luò)模型的隱藏層單元數(shù)=100，最大滯后階數(shù)P=10，正則系數(shù)λ=0.14，學(xué)習(xí)率lr=0.001，訓(xùn)練步長為5 600。

訓(xùn)練完成后提取格蘭杰因果矩陣如圖7 所示，然后繪制ROC 曲線和PR 曲線如圖8 所示。由圖8 進(jìn)一步計(jì)算ROC 曲線和PR 曲線的面積，即AUROC、AUPR 分別為0.915和0.876。

Fig.7 Granger causal matrix inferred by GRUGC圖7 GRUGC的格蘭杰因果矩陣

Fig.8 PR and ROC curve in the simulation of nonlinear VAR model圖8 非線性VAR模型仿真的PR、ROC曲線

式中：xi，t為第i個節(jié)點(diǎn)在t時刻的值；i=1，2，…，10；εi，t～N(0，0.012)。

3.4 非均勻嵌入時滯VAR

根據(jù)式（16）研究一種非均勻嵌入時滯VAR 模型，其與非線性VAR 模型的不同之處在于節(jié)點(diǎn)間存在多階滯后關(guān)系，生成T=1 000 的多變量時間序列矩陣，非均勻嵌入時滯VAR 模型基準(zhǔn)網(wǎng)絡(luò)如圖9 所示。網(wǎng)絡(luò)訓(xùn)練采用Adam優(yōu)化方法，GRU 網(wǎng)絡(luò)模型的隱藏層單元數(shù)D=100，最大滯后階數(shù)P=5，正則系數(shù)λ=0.14，學(xué)習(xí)率lr=0.001，訓(xùn)練步長為5 600。訓(xùn)練完成后提取格蘭杰因果矩陣如圖10 所示，然后繪制ROC 曲線和PR 曲線如圖11 所示，經(jīng)過計(jì)算AUROC、AUPR 分別為0.904和0.921。仿真實(shí)驗(yàn)表明，GRUGC可處理非均勻嵌入時滯的復(fù)雜非線性數(shù)據(jù)，且在整體上性能較好。

Fig.9 Benchmark network of non-uniformly embedded time-delay VAR model圖9 非均勻嵌入時滯VAR模型基準(zhǔn)網(wǎng)絡(luò)

Fig.10 Granger causality matrix for non-uniformly embedded timedelay VAR model圖10 非均勻嵌入時滯VAR模型的格蘭杰因果矩陣

Fig.11 PR，ROC curve in the simulation of non-uniformly embedded time-delay VAR model圖11 非均勻嵌入時滯VAR模型仿真的PR、ROC曲線

式中：xi，t為第i個節(jié)點(diǎn)在t時刻的值；i=1，2，…，5；εi，t～N(0，0.012)。

3.5 DREAM3數(shù)據(jù)實(shí)驗(yàn)

為了驗(yàn)證GRUGC 在實(shí)際網(wǎng)絡(luò)數(shù)據(jù)集上的性能，基于DREAM3 挑戰(zhàn)賽中的兩個Ecoli 數(shù)據(jù)集和3 個Yeast 數(shù)據(jù)集（https：//doi.org/10.7303/syn2853594）進(jìn)行仿真研究。實(shí)驗(yàn)采取N=10 的小網(wǎng)絡(luò)模型，在5 個數(shù)據(jù)集中每個數(shù)據(jù)集分別包含10 個不同的時間序列。針對每個時間序列重復(fù)試驗(yàn)4次，每次采集21個時間點(diǎn)，總共為84個樣本點(diǎn)。

神經(jīng)網(wǎng)絡(luò)訓(xùn)練方法采用Adam 優(yōu)化方法，GRU 網(wǎng)絡(luò)模型的隱藏層單元數(shù)D=5，正則系數(shù)λ=0.18，學(xué)習(xí)率lr=0.001，訓(xùn)練步長為5 000。首先畫出ROC 曲線和PR 曲線，如圖12 所示；然后計(jì)算AUPR、AUROC，并分別與DREAM3挑戰(zhàn)賽的最終獲獎榜單的第一名bteam、第二名Team291和第三名Team304 進(jìn)行比較；最后畫出AUPR 和AUROC 的柱狀圖，如圖13所示。

Fig.12 PR curve and ROC curve in the simulation of DREAM3圖12 DREAM3仿真的PR曲線和ROC曲線

Fig.13 Performance comparison among GRUGC and top three methods in DREAM3圖13 GRUGC與DREAM3前三名方法的性能比較

由圖13 可見，GRUGC 在Yeast2、Yeast3 數(shù)據(jù)集 上的AUROC 和AUPR 均超過前3 名，而在Yeast1 數(shù)據(jù)集上GRUGC 的AUROC 不如前3 名，但AUPR 超過了Team304。同時，在Ecoli1、Ecoli2 數(shù)據(jù)集上GRUGC 的AUPR 與bteam獲得的AUPR 值更接近，同時GRUGC 的AUROC 也能保持不錯的效果。

表5 為對5 個數(shù)據(jù)集的AUROC、AUPR 取平均值進(jìn)行綜合分析。由此可知，GRUGC 的整體性能超過了第二名Team291，僅次于第一名bteam，尤其是AUPR 的平均值與bteam 值非常接近。通過DREAM3 挑戰(zhàn)賽的研究分析發(fā)現(xiàn)，GRUGC 方法性能較好，具有一定的實(shí)際競爭力。

Table 5 Comparison analysis of mean AUROC and AUPR表5 AUROC和AUPR的均值對比分析

4 結(jié)語

網(wǎng)絡(luò)重構(gòu)旨在基于測量得到的數(shù)據(jù)推斷網(wǎng)絡(luò)節(jié)點(diǎn)間的相互作用關(guān)系，是分析系統(tǒng)動力學(xué)行為、結(jié)構(gòu)特性和影響機(jī)制的前提和基礎(chǔ)。本文將GRU 神經(jīng)網(wǎng)絡(luò)模型和格蘭杰因果理論相結(jié)合，提出基于GRU 網(wǎng)絡(luò)的格蘭杰因果網(wǎng)絡(luò)重構(gòu)方法（GRUGC）。該方法考慮了變量間影響關(guān)系的非線性和因果性，從建模和分析中發(fā)現(xiàn)正則項(xiàng)對最終結(jié)果影響較大，有正則項(xiàng)的優(yōu)化方法Adam 的性能相較于無正則項(xiàng)的AdamU 效果更好，而且當(dāng)輸入數(shù)據(jù)為長時間序列時GRUGC 不會發(fā)生梯度消失或梯度爆炸現(xiàn)象。

首先，通過對線性VAR、非線性VAR、Lorenz-96 和非均勻嵌入時滯VAR 模型的仿真研究證實(shí)了所提方法的有效性，然后在不同網(wǎng)絡(luò)參數(shù)下對Lorenz-96 模型進(jìn)行了消融實(shí)驗(yàn)，以進(jìn)一步驗(yàn)證本文方法的性能。最后，基于DREAM3 挑戰(zhàn)賽的Ecoli 和Yeast 數(shù)據(jù)集，與最終榜單的前3 名算法進(jìn)行比較分析來驗(yàn)證GRUGC 方法的優(yōu)越性和實(shí)用性。

在GRUGC 方法中使用的正則項(xiàng)對結(jié)果影響較大，因此需要更深入地研究正則項(xiàng)的影響機(jī)制，并提供更優(yōu)的正則化方法。在應(yīng)用GRUGC 進(jìn)行網(wǎng)絡(luò)重構(gòu)時，還需考慮數(shù)據(jù)的質(zhì)量以保證結(jié)果的準(zhǔn)確性和可靠性。下一步，將深入拓展GRUGC 的適用范圍，探索圖像、文本等不同數(shù)據(jù)類型的應(yīng)用。此外，對于存在多種因果關(guān)系的網(wǎng)絡(luò)，探索如何選擇最優(yōu)的因果關(guān)系，并提供相應(yīng)的可解釋方法也是未來研究的熱點(diǎn)之一。