第 22卷第 10期2023年 10月
Vol.22 No.10
Oct.2023软件导刊
Software Guide
基于GRU网络的格兰杰因果网络重构
杨官学，王家栋
（江苏大学电气信息工程学院，江苏镇江 212013）
摘要：传统格兰杰因果依赖线性动力学，无法适应非线性应用场景的需求，因此提出一种基于GRU网络的格兰杰
因果网络重构方法。该方法将整个网络重构划分为每个目标节点的邻居节点选择问题，针对每个目标节点构建基于
GRU网络的格兰杰因果模型，在循环神经网络中引入简单的门控机制控制信息的更新方式，并对网络输入权重施加
组稀疏惩罚以提取节点间的格兰杰因果关系。然后集成每一个子网络，获得最终完整的因果网络结构，并在GRU网
络建模训练过程中考虑采用正则化的优化方法。通过线性矢量自回归、非线性矢量自回归、非均匀嵌入时滞矢量自
回归、Lorenz-96模型及DREAM3竞赛数据集的实验表明，所提网络鲁棒性较强、有效性较高，在网络重构性能上具有
明显的优越性。
关键词：网络重构；因果推断；循环神经网络；格兰杰因果；门控循环单元
DOI：10.11907/rjdk.231360开放科学（资源服务）标识码（OSID）：
中图分类号：TP183  文献标识码：A文章编号：1672-7800（2023）010-0049-09
Network Reconstruction via Granger Causality Based on GRU Network
YANG Guanxue， WANG Jiadong
（School of Electrical and Information Engineering， Jiangsu University， Zhenjiang 212013， China）
Abstract：Reconstruction method of Granger causality network based on GRU network is proposed to address the traditional Granger causality that relies on linear dynamics and cannot meet the needs of nonlinear application scenarios. This method divides the entire network reconstruc⁃tion into neighbor node selection problems for each target
node， constructs a Granger causality model based on GRU network for each target node， introduces a simple gating mechanism to control the update of information in the recurrent neural network， and applies a sparse penalty to the network input weight to extract the Granger causality between nodes. Then integrate each sub network to obtain the final complete causal network structure， and consider using regularization optimization methods during the GRU network modeling and training process. The experi⁃ments on linear vector autoregressive， nonlinear vector autoregressive， non-uniformly embedded time-delay vector autoregressive， Lorenz-96 model， and DREAM3 competition dataset show that the proposed network has strong robustness， high effectiveness， and obvious superiority in network reconstruction performance..
Key Words：network reconstruction； causal inference； recurrent neural network； Granger causality； gated recurrent unit
0 引言
现实生活中，许多复杂系统均可在网络角度被抽象表达，其中网络节点代表系统变量，连边代表各变量间的相互作用关系。基于关联的属性和特征，网络可被划分为不同种类，例如根据边是否有权重或方向，网络可被分为有权网络和无权网络、有向网络和无向网络。典型的网络包含电力网络、脑网络、生物网络、交通网、人际网［1-2］等。然而，实际的网络结构往往是未知或难以直接观察。在当今大数据时代，复杂网络涵盖的信息量越来越多，如何基于隐藏在节点背后
的数据挖掘节点间的作用关系，就显得至关重要，这也使得复杂网络重构成为当前研究的热点方向之一［3-4］。复杂网络重构不仅能使人们更好地了解系统的动力学行为和演化机制，也是后续网络结构研究和分析的基础。例如，针对脑神经网络研究大脑神经元间如何相互激活［5］；针对基因调控网络从基因表达水平的时间序列推断基因调控网络［6］等。
收稿日期：2023-04-06
基金项目：国家自然科学基金项目（61903161）
作者简介：杨官学（1984-），男，博士，江苏大学电气信息工程学院讲师，研究方向为复杂网络重构与分析、复杂系统建模与智能控制、机器学习等；王家栋（1996-），男，江苏大学电气信息工程学院硕士研究生，研究方向为深度学习、神经网络、网络重构。
2023 年软件导刊
当下，复杂网络重构算法包括相关系数［7-8］、互信息［9-11］、布尔网络［12-13］、贝叶斯网络［14］、格兰杰因果［15］、压缩感知［16］等。尤其在有向网络的推断方面，格兰杰因果（Granger Causality， GC）及其拓展形式得到了广泛应用。传统格兰杰因果是一种基于线性矢量自回归模型的判定方法，能判断两个时间序列变量间的因果关系，但当网络中存在两个以上变量时可能会导致虚假连边。为此，提出条件格兰杰因果来解决多个变量间的间接影响［17］。
此外，针对样本量不足的问题，在考虑实际网络稀疏性的前提下，很多学者将稀疏特征引入格兰杰因果中，提出LASSO格兰杰因果［18］、组稀疏格兰杰因果［19］、基于贪婪算法的格兰杰因果［20］等方法，这些方法虽然结构形式简单、运算方便，但无法匹配非线性因果关系的处理情形。因此，很多学者将核函数引入传统格兰杰因果中提出核格兰杰因果方法［21-22］，但本质上还只是针对某种特定的非线性情形，适用范围存在一定限制。
考虑到神经网络在非线性建模方面能非常擅长地处理输入和输出间复杂的非线性关系。在因果预测方面结合格兰杰因果模型、多层感知机（Multi-Layer Perceptron，
MLP）、循环神经网络（Recurrent Neural Network， RNN）及长短期记忆元（Long Short-Term Memory， LSTM）已得到了应用。Chivukula等［23］使用RNN分析雅虎财经获取的真实股市数据，结果表明因果特征显著改善了现有的深度学习回归模型。Pathod等［24］使用RNN和LSTM处理多元脑连接性检测问题提出了RNN-GC模型，能对非线性和变长时延信息传输进行建模，并在方向性脑连接性估计方面十分有效。Tank等［25］使用MLP和LSTM在时序数据上进行建模，并融合LASSO提取节点间的因果关系，两种方法均能自动探索最大滞后阶数，且在复杂非线性DREAM3数据上取得了较好的结果。
然而，MLP、RNN存在收敛慢、容易过拟合、可解释性差、计算复杂度高等缺点。其中，RNN的循环结构使其可能无法很好地处理长期依赖关系，易出现梯度消失或梯度爆炸现象；LSTM在RNN上进行改进，不易出现梯度问题，但在处理噪声较大的数据时性能较差。相较于LSTM网络，GRU网络结构仅有更新门和重置门，更容易训练。因此，考虑到变量间影
响关系的非线性和因果性，本文提出一种基于门控循环单元（Gated Recurrent Unit， GRU）网络的格兰杰因果网络重构方法（GRUGC）。首先围绕每个目标节点，构建基于GRU网络的格兰杰因果模型；然后对网络输入权重进行组稀疏惩罚约束；最后基于Adam的梯度下降网络训练法，获取节点之间的格兰杰因果关系。
在仿真验证方面，首先基于模型生成的数据集进行相关仿真研究，例如线性矢量自回归、非线性矢量自回归、非均匀嵌入时滞矢量自回归、Lorenz-96模型。然后，采用经典的DREAM3竞赛数据集（Ecoli数据集和Yeast数据集）分析模型性能。1 格兰杰因果模型
在格兰杰因果分析中，传统分析方法通常采用线性矢量自回归模型（Vector-Autoregression， VAR），设模型的最大时滞阶数为P。
X t=∑p=1P A()p X t-p+e t（1）式中：X t为模型中所有变量在t时刻的样本矩阵。可表示为：
X t=
é
ë
ê
ê
ê
êê
ê
ê
êù
û
ú
ú
ú
úú
ú
ú
ú
x1，t
x2，t
⋮
x N，t
=
é
ë
ê
ê
ê
êê
ê
ê
êù
û
ú
ú
ú
úú
ú
ú
ú
x11，t x21，t⋯x M1，t
x12，t x22，t⋯x M2，t
⋮⋮⋱⋮
x1N，t x2N，t⋯x M N，t
∈R N×M（2）
式中：N为模型变量个数；M为样本数目；x i，t∈R1×M 为第i个变量在t时刻的样本向量；x m i，t为第i个变量在t时刻的第m个样本值；i=1，2，…，N；m=1，2，…，M；p为模型阶数；A
()p∈R N×N为t-p时刻X t-p对应的系数矩阵。可表示为：
A()p=
é
ë
ê
ê
ê
êê
ê
ê
êù
û
ú
ú
ú
úú
ú
ú
ú
A(p)1，1A(p)1，2⋯A(p)1，N
A(p)2，1A(p)2，2⋯A(p)2，N
⋮⋮⋱⋮
A(p)N，1A(p)N，2⋯A(p)N，N
∈R N×N（3）式中：p=1，2，…，P；A
()p
i，j
表示在t-p时刻第i个变量和第j个变量间的影响关系；X t-p∈R N×M为模型变量在t-p时刻的样本矩阵，形式如X t；e t∈R N×M为服从标准正态分布的噪声。
在VAR中，当且仅当对所有阶数p有A
()p
i，j=0，即可获得变量j不是影响变量i的格兰杰原因；反之，变量j是影响变量i的格兰杰原因，即存在一条从变量i指向变量j的有向边。在高维情况下考虑到稀疏性，格兰杰因果模型可被认为是一个带有组稀疏正则化的回归问题。
min
A，A，…，A
X t-∑p=1P A(p)X t-p22+λ∑i，j
A(1)i，j，…，A(P)i，j
2
（4）式中：  •2为ℓ2正则化，目的是防止模型发生过拟合现象；λ为惩罚系数。
2 神经网络格兰杰因果模型
为了解决传统格兰杰因果模型无法有效处理非线性关系的问题，提出了一种基于神经网络的格兰杰因果模型。首先构建一个通用的非线性矢量自回归模型，假设该模型中变量的个数为N，最大时滞阶数为P，则模型t时刻的输出y t的通用表达式如式（5）所示。
y t=f(x-t)+e t（5）式中：f(•)为输入输出间的非线性映射函数关系；x-t为
··50
第 10 期杨官学，王家栋：基于GRU 网络的格兰杰因果网络重构在t 之前时刻模型各输入时滞分量的集合，即x -t =
{x t -1，x t -2，…，x t -P }，x t -p ={x 1，t -p ，x 2，t -p ，…，x N ，t -p }，1≤p ≤P ；e t 为服从标准正态分布的高斯白噪声。
目前，通常使用神经网络对非线性函数f (•)
进行预测，常见方法包括MLP 、RNN 、LSTM 等。由于神经网络为黑盒模型，且模型中各输入共享隐藏层，难以直接进行因果推断的分析研究。因此，在式（5）中并未直接采用常规的多输入多输出映射来获取时序变量间的关系。
为了进一步详细阐述模型的原理，将整个网络重构任务分解为每个目标节点的邻居节点选择问题。针对每个目标节点i 分别使用单独的模型f i (•)
，以清晰地提取与邻居节点间的格兰杰因果关系。节点i 在t 时刻的表达值x i ，t 由式（6）所示。
x i ，t =f i (x
-t
)+e
i ，t
（6）
式中：f i (•)
采用的是GRU 网络。
同理，当且仅当所有阶数p (
1≤p ≤P )
在有关x i ，t 的预测时不依赖节点变量j 的时滞分量，即在预测模型中加入
{x
j ，t -1
，x j ，t -2，…，x j ，t -P }
后并不能提升对x i ，t 的预测精度，
则认为节点j 不是影响节点i 的原因；反之，节点j 是影响节点i 的原因，
即存在j →i 。2.1　基于GRU 的循环神经网络
RNN 属于一种具有短期记忆能力的神经网络，尤其适
用于处理和预测时间序列数据，作为一种常用的深度学习环路网络结构模型，RNN 由一个或多个循环层组成，每个循环层包含多个神经元，不仅能接收其他神经元信息，还能接收自身信息。理论上，RNN 可逼近任意的复杂非线性
动力系统。在RNN 网络参数训练过程中，由于链式法则，随着错误信息的反向传播，当输入序列时间步长较长时，梯度值可能会趋近于0（梯度消失）或非常大（梯度爆炸）。这两种情况均会导致训练效果不佳。为了解决RNN 的长程依赖问题，在本文提出神经网络格兰杰因果模型的框架中，采用了基于门控循环单元（Gated Recurrent Unit ， GRU ）的循环神经网络结构模型，如图1所示。
GRU 内部结构主要由两个门构成，分别为重置门
r t ∈[0，1]D
（reset gate ）和更新门z t ∈[0，1]D
（update gate ），
通过调节这两个门的激活值来控制信息流动，解决梯度消
失问题，以能更好地捕捉序列中的长期依赖关系。GRU 单元的内部计算式如式（7）—式（10）所示：
r t =σ(W r x t +U r h t -1+b r
)
（7）z t
=σ(W z x t
+U z
h t -1
+b z
)（8）h 't
=tanh (W h
x t
+U h
(r t
⊙h t -1
)+b h
)（9）h t
=(1-z t
)
⊙h 't
+z t
⊙h
t -1
（10）
式中：x t ∈R N ×1为t 时刻网络输入；h t ∈R D ×1为t 时刻
隐含层的状态；W *∈R D ×N 、U *∈R D ×D 、b *∈R D ×1为待学习的神经网络参数；*属于集合{r ，z ，h }中的元素；
⊙为Ha⁃damard Product ，表示矩阵对应位置元素的乘积。
重置门r t 和更新门z t 状态的迭代更新来源于上一时刻
的隐藏层h t -1和当前时刻的输入节点x t ，待获得h t -1、x t 后分别通过式（7）、式（8）计算重置门r t 和更新门z t 的值。式
（9）首先使用r t ⊙h t -1更新数据得到h 't -1；然后与各自权重相乘后将W h x t 与U h h 't -1相加完成特征融合；
最后使用tanh 函数将特征值映射到［-1，1］，得到h 't 、h 't 中具有上一时刻的隐含层状态h t -1和当前输入状态x t 的信息。
式（10）为GRU 最重要的一步，在这个阶段遗忘和记忆同时进行，既忘记h t -1中部分信息，又记住h 't 中节点输入的部分信息。GRU 网络直接用一个门来控制输入和遗忘间的平衡关系，当z t 越接近1表示记住的数据越多，反之遗忘的越多。
针对时间序列输入x t ，将x t 在不同时刻输入到循环神经网络中得到相应的隐含层状态h t ，循环神经网络的输出y t 最终由h t 的加权线性求和表示，
如式（11）所示。y t =f i (x -t
)正则匹配关键词
+e t =W T o h t +e t （11）
式中：W o ∈R D ×1为输出层权重；W o 为待学习的网络
参数。
2.2　基于GRU 网络的格兰杰因果
在因果网络重构问题中，根据式（6）、式（11）针对每个目标节点i 构建子任务优化模型，分别得出相应的邻居节点。具体表达形式如下：
min W
Loss total +λ∑j =1
N
M ：，j
2
（12）
式中：W ={W r ，W z ，W h ，U r ，U z ，U h ，W o }为待求解的参
数集合；Loss total =∑t =2
T (
x i ，t -f i (x -t
))
2
为总训练误差；T 为时
间序列的总长度。
针对输入节点变量x t ，将输入层至隐含层的权重矩阵
W r 、W z 、W h 进行堆叠拼接以便于后续优化，由M =[W T r W T z W T h ]T ∈R 3D ×N 表示总输入权重矩阵。为了捕
捉驱动节点j 对目标节点i 的因果影响，当且仅当总输入权重矩阵M 的第j 列均为0，即M ：，j =0得出驱动节点j 对目标节点i 不存在因果关系；反之驱动节点j 是影响目标节点i 的原因，即存在j →i 的连边。此外，合理加入正则项能提升神经网络模型的泛化能力，ℓ2正则化项的惩罚因子λ
同
Fig. 1　Gated loop unit structure of GRU networks
图1　GRU 网络的门控循环单元结构
·
·51
2023 年
软件导刊
时也是控制网络结构稀疏性的参数［26］。
2.3　基于Adam的梯度下降法
当下，最流行的神经网络学习优化方法大都始于随机
梯度下降法（Stochastic Gradient Descent， SGD）［27］，但SGD
每次只用一个样本更新梯度，使得SGD并非每次迭代均向
最优方向更新参数，容易造成准确度下降且无法线性收敛
的情况。同时，由于只用一个样本更新梯度并不能代表全
部样本的趋势，因此易陷入局部最小值。
为此，本文采用基于以往梯度信息和动量的Adam［28］
优化方法。该方法是RMSProp和动量法的结合，主要纠正
两项偏差和平均梯度滑动的方法，具体计算过程如下：
步骤1：初始化学习率lr、平滑常数β1、β2（分别用于平
滑m t、v t），可学习参数θ0=0、m0=0、v0=0、t=0。
步骤2：当未停止训练时，更新训练次数t=t+1，计
算梯度g t（所有的可学习参数都有各自梯度，因此g t指全
部梯度的集合）。
m t=β1m t-1+(1-β1)g t
v t=β2v t-1+(1-β2)(g t)2
m t=m t
1-()β1t
v t=v t1-()β2t
θt=θt-1-m t
v t+ϵlr
（13）
本文神经网络模型的训练方法采用了以Adam为基础的优化方法。为了研究正则项的影响，训练采用了无正则项的AdamU和有正则项的Adam。
3 仿真实验
为了验证本文所提基于GRU网络的格兰杰因果网络重构算法的有效性，分别基于线性VAR、非线性VAR、
Lorenz-96、非均匀嵌入时滞VAR模型及DREAM竞赛数据集进行仿真实验研究。网络重构性能的具体量化性能指标采用ROC曲线下的面积（Area Under Receiver-Operat⁃ing-Characteristic Curve， AUROC）和PR曲线下的面积（Ar⁃ea Und
er Precision-Recall Curve， AUPR）。
3.1　线性VAR
首先基于N=10的VAR网络模型，利用式（1）随机生成的稀疏转移矩阵生成仿真时间步长T=1 000的时间序列数据。噪声服从高斯分布N(0，σ2)，σ=0.1，基准因果网络如图2所示。由此可见，若节点与节点间存在GC关系，则对应区间的数值大小为1；若不存在GC关系，则对应区间的数值大小为0。例如，图2中（1，1）位置和（1，9）位置的数值大小均为1，分别表示存在一条由节点1指向节点1的GC边和存在一条由节点1指向节点9的GC边。
具体参数设置如下：针对每一个目标节点i，创建自我依赖关系（自环），并在其他N-1个节点中随机选择一个驱动节点构建因果关系，即在j→i的情况下设置A
()p
ij= 0.1，p=1，2，3，同时令其他A()p ij=0。GRU网络模型的隐藏
层单元数D设定为100，分别比较Adam、AdamU两种训练优化方法，最大滞后阶数P=5，正则系数λ=0.002，学习率lr =0.05，训练步长为20 000，AdamU无需设置正则系数λ。
VAR模型的两种训练方法仿真结果如图3所示，图中区域数值大小表示经过GRUGC推理后存在GC边的概率值。与图2的基准网络比较可见，考虑正则化的Adam相较于考虑正则化的AdamU效果更好。
为了量化每种训练方法的性能，计算图3不同方法的AUROC和AUPR，具体数值如表1所示。由此可知，Adam 的AUROC和AUPR均达到1，实现了完美的因果网络重构，而AdamU的AUROC、AUPR仅为0.225、0.125。
3.2　Lorenz-96
设置网络节点个数N=10，F=10，添加服从高斯分布的噪声N(0，0.012)，获得序列长度T=1 000的多变量时间序列矩阵。Lorenz-96模型的微分方程表达式（14）所示，可得N=10的Lorenz-96模型的基准网络如图4所示。
dx i
dt=(x i+1-x i-2)x i-1-x i+F（14）式中：i=1，2，…，N；F为惩罚力度，值越大表示时间序列的非线性越强。
设置GRU网络模型的隐藏层单元数D=100，比较Ad⁃am和AdamU两种训练优化方法，最大滞后阶数P=10，正则系数λ=0.2，学习率lr=0.001，训练步长为10 000。Lorenz-96模型的两种训练方法结果及相关性能指标如图5、表2所示。由此可见，经AdamU训练获取的网络结构较差，Adam获取的网络结构几乎全部预测正确；表2中Ad⁃amU方法的AUROC和AUPR 分别为0.825、0.744，Adam的AUROC和AUPR均等于1，达到了完美重构效果。
通过线性VAR模型和Lorenz-96模型的仿真发现，有无正则项对模型最终的训练结果影响较大，因此在仿真后Fig. 2　Benchmark network of VAR model
图2　VAR模型基准网络
··52
第 10 期杨官学，王家栋：基于GRU 网络的格兰杰因果网络重构续部分将仅基于Adam 方法进行训练。
此外，为了考察不同非线性程度F 和隐藏层单元数D
对模型的影响，基于Adam 方法对不同F 和D 数值进行消融实验，比较结果如表3、表4所示。由此可知，在不同D 的情况下随着F 增加，AUROC 、AUPR 均会下降；在F =
10
Fig. 5　Results of two training methods for Lorenz-96 model
图5　Lorenz-96
模型的两种训练方法结果
Fig. 4　Benchmark network of Lorenz-96 model
图4　Lorenz-96模型的基准网络
Table 4　AUPR for different F and D 表4　不同F 与D 情况下的AUPR
F\\D 1020304050.9850.7980.7240.708100.9980.7930.7690.751
250.9990.8140.8090.796
500.9930.8610.8220.7961000.9910.8690.8480.802
Table 3　AUROC for different F and D 表3　不同F 与D 情况下的AUROC F\\D
10203040
50.9910.8890.8280.783
100.9980.8960.8710.821
250.9990.8970.8950.874
500.9960.9120.8970.881
1000.9940.9150.9080.879
Table 2　Performance comparison of two training methods for
Lorenz-96 model
表2　Lorenz-96模型的两种训练方法性能比较
Method Adam
AdamU AUROC 1.0000.825
AUPR 1.000
0.744
Fig. 3　Results of two training methods for VAR model
图3　VAR 模型的两种训练方法结果
Table 1　Performance comparison of two training methods for VAR
model
表1　VAR 模型的两种训练方法性能比较
Method Adam
AdamU AUROC 1.0000.225
AUPR 1.0000.125
·
·53

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。