楊官學(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)系?；陉P(guān)聯(lián)的屬性和特征，網(wǎng)絡(luò)可被劃分為不同種類(lèi)，例如根據(jù)邊是否有權(quán)重或方向，網(wǎng)絡(luò)可被分為有權(quán)網(wǎng)絡(luò)和無(wú)權(quán)網(wǎng)絡(luò)、有向網(wǎng)絡(luò)和無(wú)向網(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ù)時(shí)代，復(fù)雜網(wǎng)絡(luò)涵蓋的信息量越來(lái)越多，如何基于隱藏在節(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)的動(dòng)力學(xué)行為和演化機(jī)制，也是后續(xù)網(wǎng)絡(luò)結(jié)構(gòu)研究和分析的基礎(chǔ)。例如，針對(duì)腦神經(jīng)網(wǎng)絡(luò)研究大腦神經(jīng)元間如何相互激活［5］；針對(duì)基因調(diào)控網(wǎng)絡(luò)從基因表達(dá)水平的時(shí)間序列推斷基因調(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)格蘭杰因果是一種基于線(xiàn)性矢量自回歸模型的判定方法，能判斷兩個(gè)時(shí)間序列變量間的因果關(guān)系，但當(dāng)網(wǎng)絡(luò)中存在兩個(gè)以上變量時(shí)可能會(huì)導(dǎo)致虛假連邊。為此，提出條件格蘭杰因果來(lái)解決多個(gè)變量間的間接影響［17］。

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

考慮到神經(jīng)網(wǎng)絡(luò)在非線(xiàn)性建模方面能非常擅長(zhǎng)地處理輸入和輸出間復(fù)雜的非線(xiàn)性關(guān)系。在因果預(yù)測(cè)方面結(jié)合格蘭杰因果模型、多層感知機(jī)（Multi-Layer Perceptron，MLP）、循環(huán)神經(jīng)網(wǎng)絡(luò)（Recurrent Neural Network，RNN）及長(zhǎng)短期記憶元（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 處理多元腦連接性檢測(cè)問(wèn)題提出了RNN-GC 模型，能對(duì)非線(xiàn)性和變長(zhǎng)時(shí)延信息傳輸進(jìn)行建模，并在方向性腦連接性估計(jì)方面十分有效。Tank 等［25］使用MLP 和LSTM 在時(shí)序數(shù)據(jù)上進(jìn)行建模，并融合LASSO 提取節(jié)點(diǎn)間的因果關(guān)系，兩種方法均能自動(dòng)探索最大滯后階數(shù)，且在復(fù)雜非線(xiàn)性DREAM3 數(shù)據(jù)上取得了較好的結(jié)果。

然而，MLP、RNN 存在收斂慢、容易過(guò)擬合、可解釋性差、計(jì)算復(fù)雜度高等缺點(diǎn)。其中，RNN 的循環(huán)結(jié)構(gòu)使其可能無(wú)法很好地處理長(zhǎng)期依賴(lài)關(guān)系，易出現(xiàn)梯度消失或梯度爆炸現(xiàn)象；LSTM 在RNN 上進(jìn)行改進(jìn)，不易出現(xiàn)梯度問(wèn)題，但在處理噪聲較大的數(shù)據(jù)時(shí)性能較差。相較于LSTM 網(wǎng)絡(luò)，GRU 網(wǎng)絡(luò)結(jié)構(gòu)僅有更新門(mén)和重置門(mén)，更容易訓(xùn)練。因此，考慮到變量間影響關(guān)系的非線(xiàn)性和因果性，本文提出一種基于門(mén)控循環(huán)單元（Gated Recurrent Unit，GRU）網(wǎng)絡(luò)的格蘭杰因果網(wǎng)絡(luò)重構(gòu)方法（GRUGC）。首先圍繞每個(gè)目標(biāo)節(jié)點(diǎn)，構(gòu)建基于GRU 網(wǎng)絡(luò)的格蘭杰因果模型；然后對(duì)網(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)仿真研究，例如線(xiàn)性矢量自回歸、非線(xiàn)性矢量自回歸、非均勻嵌入時(shí)滯矢量自回歸、Lorenz-96 模型。然后，采用經(jīng)典的DREAM3 競(jìng)賽數(shù)據(jù)集（Ecoli 數(shù)據(jù)集和Yeast 數(shù)據(jù)集）分析模型性能。

1 格蘭杰因果模型

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

式中：Xt為模型中所有變量在t時(shí)刻的樣本矩陣?？杀硎緸椋?/p>

式中：N為模型變量個(gè)數(shù)；M為樣本數(shù)目；xi，t∈R1×M為第i個(gè)變量在t時(shí)刻的樣本向量為 第i個(gè) 變量 在t時(shí)刻的第m個(gè)樣本值；i=1，2，…，N；m=1，2，…，M；p為模型階數(shù)為t-p時(shí)刻Xt-p對(duì)應(yīng)的系數(shù)矩陣?？杀硎緸椋?/p>

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

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

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

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

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

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

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

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

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

同理，當(dāng)且僅當(dāng)所有階數(shù)p(1 ≤p≤P)在有關(guān)xi，t的預(yù)測(cè)時(shí)不依賴(lài)節(jié)點(diǎn)變量j的時(shí)滯分量，即在預(yù)測(cè)模型中加入{xj，t-1，xj，t-2，…，xj，t-P}后并不能提升對(duì)xi，t的預(yù)測(cè)精度，則認(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ù)測(cè)時(shí)間序列數(shù)據(jù)，作為一種常用的深度學(xué)習(xí)環(huán)路網(wǎng)絡(luò)結(jié)構(gòu)模型，RNN 由一個(gè)或多個(gè)循環(huán)層組成，每個(gè)循環(huán)層包含多個(gè)神經(jīng)元，不僅能接收其他神經(jīng)元信息，還能接收自身信息。理論上，RNN 可逼近任意的復(fù)雜非線(xiàn)性動(dòng)力系統(tǒng)。在RNN 網(wǎng)絡(luò)參數(shù)訓(xùn)練過(guò)程中，由于鏈?zhǔn)椒▌t，隨著錯(cuò)誤信息的反向傳播，當(dāng)輸入序列時(shí)間步長(zhǎng)較長(zhǎng)時(shí)，梯度值可能會(huì)趨近于0（梯度消失）或非常大（梯度爆炸）。這兩種情況均會(huì)導(dǎo)致訓(xùn)練效果不佳。為了解決RNN 的長(zhǎng)程依賴(lài)問(wèn)題，在本文提出神經(jīng)網(wǎng)絡(luò)格蘭杰因果模型的框架中，采用了基于門(mén)控循環(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ò)的門(mén)控循環(huán)單元結(jié)構(gòu)

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

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

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

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

針對(duì)時(shí)間序列輸入xt，將xt在不同時(shí)刻輸入到循環(huán)神經(jīng)網(wǎng)絡(luò)中得到相應(yīng)的隱含層狀態(tài)ht，循環(huán)神經(jīng)網(wǎng)絡(luò)的輸出yt最終由ht的加權(quán)線(xià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)問(wèn)題中，根據(jù)式（6）、式（11）針對(duì)每個(gè)目標(biāo)節(jié)點(diǎn)i構(gòu)建子任務(wù)優(yōu)化模型，分別得出相應(yīng)的鄰居節(jié)點(diǎn)。具體表達(dá)形式如下：

針對(duì)輸入節(jié)點(diǎn)變量xt，將輸入層至隱含層的權(quán)重矩陣Wr、Wz、Wh進(jìn)行堆疊拼接以便于后續(xù)優(yōu)化，由M=表示總輸入權(quán)重矩陣。為了捕捉驅(qū)動(dòng)節(jié)點(diǎn)j對(duì)目標(biāo)節(jié)點(diǎn)i的因果影響，當(dāng)且僅當(dāng)總輸入權(quán)重矩陣M的第j列均為0，即M：，j=0得出驅(qū)動(dòng)節(jié)點(diǎn)j對(duì)目標(biāo)節(jié)點(diǎn)i不存在因果關(guān)系；反之驅(qū)動(dòng)節(jié)點(diǎn)j是影響目標(biāo)節(jié)點(diǎn)i的原因，即存在j→i的連邊。此外，合理加入正則項(xiàng)能提升神經(jīng)網(wǎng)絡(luò)模型的泛化能力，?2正則化項(xiàng)的懲罰因子λ同時(shí)也是控制網(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每次只用一個(gè)樣本更新梯度，使得SGD 并非每次迭代均向最優(yōu)方向更新參數(shù)，容易造成準(zhǔn)確度下降且無(wú)法線(xiàn)性收斂的情況。同時(shí)，由于只用一個(gè)樣本更新梯度并不能代表全部樣本的趨勢(shì)，因此易陷入局部最小值。

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

步驟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)練時(shí)，更新訓(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)練采用了無(wú)正則項(xiàng)的AdamU 和有正則項(xiàng)的Adam。

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

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

3.1 線(xiàn)性VAR

首先基于N=10 的VAR 網(wǎng)絡(luò)模型，利用式（1）隨機(jī)生成的稀疏轉(zhuǎn)移矩陣生成仿真時(shí)間步長(zhǎng)T=1 000 的時(shí)間序列數(shù)據(jù)。噪聲服從高斯分布N(0，σ2)，σ=0.1，基準(zhǔn)因果網(wǎng)絡(luò)如圖2 所示。由此可見(jiàn)，若節(jié)點(diǎn)與節(jié)點(diǎn)間存在GC 關(guān)系，則對(duì)應(yīng)區(qū)間的數(shù)值大小為1；若不存在GC 關(guān)系，則對(duì)應(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è)置如下：針對(duì)每一個(gè)目標(biāo)節(jié)點(diǎn)i，創(chuàng)建自我依賴(lài)關(guān)系（自環(huán)），并在其他N-1 個(gè)節(jié)點(diǎn)中隨機(jī)選擇一個(gè)驅(qū)動(dòng)節(jié)點(diǎn)構(gòu)建因果關(guān)系，即在j→i的情況下設(shè)置0.1，p=1，2，3，同時(shí)令其他=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)練步長(zhǎng)為20 000，AdamU 無(wú)需設(shè)置正則系數(shù)λ。

VAR 模型的兩種訓(xùn)練方法仿真結(jié)果如圖3 所示，圖中區(qū)域數(shù)值大小表示經(jīng)過(guò)GRUGC 推理后存在GC 邊的概率值。與圖2 的基準(zhǔn)網(wǎng)絡(luò)比較可見(jiàn)，考慮正則化的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)個(gè)數(shù)N=10，F(xiàn)=10，添加服從高斯分布的噪聲N(0，0.012)，獲得序列長(zhǎng)度T=1 000 的多變量時(shí)間序列矩陣。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為懲罰力度，值越大表示時(shí)間序列的非線(xiàn)性越強(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)練步長(zhǎng)為10 000。Lorenz-96 模型的兩種訓(xùn)練方法結(jié)果及相關(guān)性能指標(biāo)如圖5、表2 所示。由此可見(jiàn)，經(jīng)AdamU 訓(xùn)練獲取的網(wǎng)絡(luò)結(jié)構(gòu)較差，Adam 獲取的網(wǎng)絡(luò)結(jié)構(gòu)幾乎全部預(yù)測(cè)正確；表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é)果

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

此外，為了考察不同非線(xiàn)性程度F和隱藏層單元數(shù)D對(duì)模型的影響，基于Adam 方法對(duì)不同F(xiàn)和D數(shù)值進(jìn)行消融實(shí)驗(yàn)，比較結(jié)果如表3、表4 所示。由此可知，在不同D的情況下隨著F增加，AUROC、AUPR 均會(huì)下降；在F=10時(shí)隨著D增加AUROC 和AUPR 變化不大；當(dāng)F={20，30，40}時(shí)，隨著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ù)的非線(xiàn)性增強(qiáng)，GRU 網(wǎng)絡(luò)的擬合能力變?nèi)?。在非線(xiàn)性程度較低時(shí)，增加隱藏層單元數(shù)無(wú)法大幅度提升模型精度；在非線(xiàn)性程度較高時(shí)，增加隱藏層單元數(shù)能明顯提升模型精度，尤其在D從5 增加到10、從10 增加到25 時(shí)，AUROC 和AUPR 提升最明顯，這也符合實(shí)際經(jīng)驗(yàn)，在數(shù)據(jù)較復(fù)雜時(shí)模型精度往往會(huì)降低，可通過(guò)采用更復(fù)雜的模型提升表達(dá)能力。并且，在神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)中增加網(wǎng)絡(luò)深度和寬度是提高精度的常見(jiàn)方法［29］。

3.3 非線(xiàn)性VAR

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

Fig.6 Benchmark network of nonlinear VAR model圖6 非線(xiàn)性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)練步長(zhǎng)為5 600。

訓(xùn)練完成后提取格蘭杰因果矩陣如圖7 所示，然后繪制ROC 曲線(xiàn)和PR 曲線(xiàn)如圖8 所示。由圖8 進(jìn)一步計(jì)算ROC 曲線(xiàn)和PR 曲線(xiàn)的面積，即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 非線(xiàn)性VAR模型仿真的PR、ROC曲線(xiàn)

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

3.4 非均勻嵌入時(shí)滯VAR

根據(jù)式（16）研究一種非均勻嵌入時(shí)滯VAR 模型，其與非線(xiàn)性VAR 模型的不同之處在于節(jié)點(diǎn)間存在多階滯后關(guān)系，生成T=1 000 的多變量時(shí)間序列矩陣，非均勻嵌入時(shí)滯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)練步長(zhǎng)為5 600。訓(xùn)練完成后提取格蘭杰因果矩陣如圖10 所示，然后繪制ROC 曲線(xiàn)和PR 曲線(xiàn)如圖11 所示，經(jīng)過(guò)計(jì)算AUROC、AUPR 分別為0.904和0.921。仿真實(shí)驗(yàn)表明，GRUGC可處理非均勻嵌入時(shí)滯的復(fù)雜非線(xiàn)性數(shù)據(jù)，且在整體上性能較好。

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

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

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

式中：xi，t為第i個(gè)節(jié)點(diǎn)在t時(shí)刻的值；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)賽中的兩個(gè)Ecoli 數(shù)據(jù)集和3 個(gè)Yeast 數(shù)據(jù)集（https：//doi.org/10.7303/syn2853594）進(jìn)行仿真研究。實(shí)驗(yàn)采取N=10 的小網(wǎng)絡(luò)模型，在5 個(gè)數(shù)據(jù)集中每個(gè)數(shù)據(jù)集分別包含10 個(gè)不同的時(shí)間序列。針對(duì)每個(gè)時(shí)間序列重復(fù)試驗(yàn)4次，每次采集21個(gè)時(shí)間點(diǎn)，總共為84個(gè)樣本點(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)練步長(zhǎng)為5 000。首先畫(huà)出ROC 曲線(xiàn)和PR 曲線(xiàn)，如圖12 所示；然后計(jì)算AUPR、AUROC，并分別與DREAM3挑戰(zhàn)賽的最終獲獎(jiǎng)榜單的第一名bteam、第二名Team291和第三名Team304 進(jìn)行比較；最后畫(huà)出AUPR 和AUROC 的柱狀圖，如圖13所示。

Fig.12 PR curve and ROC curve in the simulation of DREAM3圖12 DREAM3仿真的PR曲線(xiàn)和ROC曲線(xiàn)

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

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

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

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

4 結(jié)語(yǔ)

網(wǎng)絡(luò)重構(gòu)旨在基于測(cè)量得到的數(shù)據(jù)推斷網(wǎng)絡(luò)節(jié)點(diǎn)間的相互作用關(guān)系，是分析系統(tǒng)動(dò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)系的非線(xiàn)性和因果性，從建模和分析中發(fā)現(xiàn)正則項(xiàng)對(duì)最終結(jié)果影響較大，有正則項(xiàng)的優(yōu)化方法Adam 的性能相較于無(wú)正則項(xiàng)的AdamU 效果更好，而且當(dāng)輸入數(shù)據(jù)為長(zhǎng)時(shí)間序列時(shí)GRUGC 不會(huì)發(fā)生梯度消失或梯度爆炸現(xiàn)象。

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

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