动态正则格兰杰因果学习方法

Dynamic Regularized Granger Causality Learning Method

摘要: 在医学和金融学等实际领域中，了解动态系统中的底层结构关系对于调节系统中的变量和预测系统未来状态至关重要。系统的动态变化会生成时间序列数据，通过观察时间序列数据可以分析系统的底层结构。格兰杰因果关系分析方法可以应用于一维或多维时间序列系统，现有的方法以组件式的建模方式分析每个系统变量特定的因果关系，受限于时间方向的强假设性和组件模型的单一性，其无法准确地挖掘出时间序列中的因果关系结构。本文提出了一种基于动态稀疏正则化的格兰杰因果发现方法DRGC (Dynamic Regularity Granger Causlity)。DRGC模型从卷积网络的输入权重中周期性地发掘变量在时间维度上的依赖信息，并以此为据向网络施加稀疏惩罚，以获得精确的格兰杰因果关系；同时，使用采样输入的循环网络提取数据中的长程依赖关系，同步优化卷积网络的权重，增强了模型发现因果关系的精确性和稳定性。在模拟数据集和真实系统生成的数据集上的实验表明，DRGC优于最先进的基线方法。

Abstract: In practical fields such as medicine and finance, understanding the underlying structural relationships in dynamic systems is crucial for regulating system variables and predicting the system’s future state. The dynamic changes of a system generate time series data, and by observing these time series, the underlying structure of the system can be analyzed. Granger causality analysis methods can be applied to univariate or multivariate time series systems. Existing methods analyze the specific causal relationships of each system variable using a modular modeling approach. However, these methods are limited by strong assumptions regarding the time direction and the simplicity of the modular models, which prevents them from accurately uncovering the causal relationship structure in multivariate systems. This paper proposes a Granger causality discovery method based on dynamic sparse regularization, called DRGC (Dynamic Regularized Granger Causality). The DRGC model periodically uncovers the temporal dependencies of variables from the input weights of a convolutional network and applies sparse penalties to the network accordingly, to obtain precise Granger causal relationships. Additionally, a cyclic network is used to extract long-range dependencies from the sampled input data, and the convolutional network’s weights are optimized simultaneously, which enhances the precision and stability of the model in discovering causal relationships. Experiments conducted on simulated datasets and real-world system-generated datasets show that DRGC outperforms state-of-the-art baseline methods.

1. 引言

从多维时间序列数据中挖掘变量间的因果关系、学习动态系统的内部结构，对时间序列的精准预测和干预分析等都有极大帮助。这种内部结构揭示了关于变量内部和变量之间的同期和滞后关系的信息。例如，在神经科学领域中，大脑活动通过各脑区域传播，不同区域的检查指标会随活动的传播而产生波动[1]，而定义这些指标的内部结构和活动的传播逻辑就是产生真实数据的原因。过去几十年，研究人员致力于从观测时间序列数据中发现因果关系[2]，取得了很大进展[3]-[6]。在诸多因果发现的方法中，格兰杰因果关系分析是最常用的一种框架[7] [8]，能够量化一个时间序列的过去值是否有助于预测另一个时间序列的未来演变趋势。

格兰杰因果关系取决于所研究的整个时间序列系统活动，故更适合分析高维的复杂序列数据。在多维的情况下，格兰杰因果分析可以捕获变量与其他变量过去值之间的因果关系并计算其大小。格兰杰因果关系分析有基于模型的方法和无模型的方法，大多数经典的基于模型方法都有“线性时间动力学”的假设，最具代表性的就是向量自回归模型[9]。在这种假设下，过去值对未来值的影响关系是线性的，线性系数的大小就是格兰杰因果效应的量化体现。使用向量自回归模型(vector autoregressive model, VAR)估计格兰杰因果关系需要设定一个最大考虑的时间滞后值，这种设定存在的问题是指定滞后值太大或太小都会影响格兰杰因果关系的评估效果。解决方法通常是通过添加Lasso或截断惩罚等策略来实现自动的滞后选择[4]，同时实现格兰杰因果关系的稀疏网络诱导，使得高维系统中每个变量只选择出少量与自己有因果关系的变量。无模型的方法不再有基于模型方法的强线性相关性假设[10]，可以在变量之间发现非线性因果关系[1] [11]，而缺陷在于估计的结果方差较大，且当序列数量增加时，会受到维度灾难影响[2]，在较高维度的场景下难以适用。

早期适应时间序列的神经网络变体在多变量时间序列预测方面表现出令人满意的效果[12] [13]。长期以来，神经网络在多维时间序列系统中的应用还存在两个问题：第一，在变量较多的情况下进行联合建模会导致巨量的网络参数，大大增加了模型训练时的开销与复杂度；第二，联合建模的方法等价于将多维时间序列系统视为矢量时间序列，其本质仍然是一种黑箱算法，内部的结构对所预测系统的结构没有相当的解释作用。以上两个问题的存在使得联合建模多维时间序列不适用于格兰杰因果关系分析。

神经格兰杰因果方法[4]的提出为解决上述问题提供了优秀的思路。它为每个变量构建单独的神经网络，使用MLPs和RNNs实现了“组件级”的因果分析。这种建模方式在大大减少网络参数数量、加快模型训练速度的同时，通过解析网络参数使得模型有了可解释性。然而，得益于循环体系结构模拟长程依赖，组件级循环网络完全回避了时滞因果选择的问题，只能获取变量间整体上的因果强度，而无法挖掘出变量之间的因果关系在时间跨度上的差异。同时，采用group Lasso或分层式的稀疏惩罚方式给数据在时间方向上的因果关系大小施加了假设，即距离当前时刻较近的时滞值的因果影响更大，忽略了实际系统中可能因为采样率差异而导致的序列值影响不一定随着时间推移而减弱的问题。

为了解决上述问题，本文提出了动态正则格兰杰因果学习模型(Dynamic Regularized Granger Causality, DRGC)，通过采样因果图的方式打通循环网络与线性网络输入层之间的通道，使得循环网络发挥其捕捉长程依赖的能力来协助时滞选择，确保在线性网络难以拟合短时间序列时提供可靠的拟合性能。同时，本文设计了一种动态的网络输入权重分级惩罚策略，以增强线性网络在时滞上选择因果关系的正确性。本文的主要贡献如下：

1) 构建单维度建模的线性因果发现网络，并对提取出的因果关系进行采样，使用采样处理的时间序列数据输入循环网络进行训练，增强了线性网络在时间维度上挖掘因果关系的合理性，并且充分利用循环网络长程依赖的特性来协助因果选择。

2) 设计基于分级群组Lasso惩罚的动态稀疏惩罚策略，模型在训练过程中自动判断不同时滞之间的数值关系，协调不同原因变量各个时滞上因果关系的大小。

3) 在具有代表性的模拟数据集和拟真数据系统上验证提出模型的有效性，与当下先进的格兰杰因果发现方法对比均取得了性能的提升。

2. 相关工作

Granger提出的格兰杰因果方法[7]，通过测试时间序列对预测另一个时间序列的帮助，被广泛应用于分析时间序列领域的因果关系。格兰杰因果分析最初假设线性模型和因果结构可以通过拟合向量自回归(Vector Autoregression, VAR)模型来发现。后来，格兰杰因果关系的概念被扩展到非线性情况[8]。由于与新兴的深度神经网络有着高度兼容性，格兰杰因果关系的研究可以被扩展到更多复杂数据的情况，分析更深层次或包括混杂在内的时间序列因果关系。

格兰杰因果分析可以分为基于模型的方法和无模型的方法，大多数基于模型的方法会假设线性相关，使用自回归模型[14]。它假设序列过去值对目标序列(被预测序列)的未来值具有线性效应，且非零系数可以量化格兰杰因果效应的大小。通过添加Lasso [15]或Group Lasso这类稀疏诱导正则化方法有助于将自回归模型中的线性格兰杰因果扩展到高维环境下[14] [16]。线性相关的假设可能会导致对实际中非线性关系的误解，并可能因为简化程度过高而产生对内部结构的不一致估计。无模型的方法可以克服在线性假设下实际观测之间的非线性依赖关系，它对潜在关系的假设最少，包括转移熵[1]或有向信息[11]。无模型方法估计的结果可能由于高自由度而具有高度不确定性，且对数据量的需求较大，不适用于维数大幅增多的情况[17]。

深度神经网络兴起并发展以来，研究人员尝试将神经网络应用于因果发现的方法中。DYNOTEARS方法[6]利用深度学习模型扩展了基于分数的方法来学习SVAR模型，也称为动态贝叶斯网络(Dynamic Bayesian Network, DBN)。随着节点数量增加，网络扩展效果非常好，但DYNOTEARS框架仍基于线性VAR模型。更多的研究人员利用神经网络来推断非线性格兰杰因果关系，克服了许多传统时间序列因果发现框架的缺点。Wu等人[18]的研究介绍了一种新的最小预测信息正则化方法，从时间序列推断因果关系，允许深度学习模型发现非线性因果关系。Xu等人[19]提出了一种基于深度神经网络的可扩展因果发现算法。Singh等人[20]进行了个体成对的格兰杰因果检验研究。为了将非线性相互作用纳入格兰杰因果关系检测，Tank等人[4]提出了一类非线性结构，包含多层感知器(MLP)和递归神经网络(RNN)。NTiCD算法[21]利用多种神经网络的能力组合，通过连续优化技术来捕获时间序列内部和之间不随时间变化的因果关系。CUTS算法[22]使用神经网络设计了潜在因果发现和数据填补两个交替阶段，用于将非线性格兰杰因果分析应用到不规则时间序列数据上，它的后续版本CUTS+ [23]通过引入C2FD技术和图神经网络提高了方法在不规则数据上的可扩展性。

在可解释性方面，神经格兰杰因果发现[4]使用多层感知机和递归神经网络，通过解析神经网络的参数来获得可解释的非线性格兰杰因果关系。获得可解释因果关系的关键在于单维度输出序列的独立模型。Horvath等人[24]的研究提出了学习型核函数LeKVAR和一种解耦时滞和个体时间序列的机制，实现了有延迟的时滞选择和因果解释，提供了更好的伸缩性。

总体来说，当下神经格兰杰因果发现方法实现了减轻模型参数量和时滞选择，同时大大增强了模型的可解释性。然而目前的方法仍然存在模型在短时间序列下性能不佳和时滞发掘不准确等问题。

3. 时间序列及其因果关系模型

3.1. 多维时间序列与因果图

多维时间序列数据可以是由一个实际的系统中多个变量随时间进行采样观测得到，其在形状上可以理解为一个二维矩阵结构。从序列内部的相互关系角度出发，可以理解为多个单维时间序列的拼接，也可以理解为一个矢量随时间变化而产生的一系列值。多维时间序列的数学模型由式(1)给出

$x_t=\left(x_{t1},x_t{}_2,\cdots ,x_{tp}\right).$ (1)

其中p是系统中的变量个数，$t=1,2,\cdots ,T$ 代表时间值。在数学模型中，由于没有考虑实际系统中的采样率因素，时间值上相邻的值仅代表数值的先后关系，而不包含时间维度上的距离信息。

在实际系统中，由于变量之间存在依赖，系统中各个变量会构成因果关系。因此可以为时间序列数据建立一个因果体系以指示变量之间相互影响的关系与程度大小，称为因果图。在数值上，因果图可以表示为一个$p\times p$ 大小的方阵，$G\in {ℝ}^{p\times p}$ 。G中的每一项$c_{ij}$ 代表变量j对变量i的影响程度，如图1(a)所示。从时间维度上看，原因变量对其结果变量的影响可以细化到不同的时刻，通过在时间维度上对全局因果图进行扩展，可以得到时滞因果图，“时滞”指时间滞后值，即当前时刻之前某时刻的值。时滞因果图是一个三维图形，可以对不同结果变量进行划分，得到单一变量的时滞因果图，如图1(b)所示。

3.2. 格兰杰因果关系模型

假设变量$x_t\in {ℝ}^p$ 是一个p维的平稳时间序列，并假设数据是在一定时间$1\le t\le T,t\in Z$ 内观测到的。时间序列的线性格兰杰因果关系一般使用VAR模型来进行研究。在VAR模型中，时间序列x在时刻t处的值$x_t$ 被认为是该序列在过去的K个时间滞后值的线性组合，如式(2)所示。

$x_t={\displaystyle \sum _{k=1}^K{A}_k{x}_{t-k}}+e_t$ (2)

其中K是一个正整数，可以被解释为假设的最大时间滞后阶数，即当$k>K$ 时，$x_{t-k}$ 在任何维度上不会影响$x_t$ 。$A_k\in {ℝ}^{p\times p}$ 是对应k滞后值的因果矩阵。$e_t\in {ℝ}^p$ 是当前时刻的噪声误差，可以服从特定的零均值噪声分布。通过对每一阶滞后的值$x_{t-k}$ 经过$A_k$ 变换得到这一阶滞后对当前时刻值的潜在影响。

Figure 1. (a) Global causal graph; (b) Lag causal graph of univariate time series

图1. (a) 时间序列数据的全局因果图；(b) 单变量时间序列时滞因果图

在类似金融系统等有高维度时间序列的领域中，数据之间往往不满足前面提到的线性关系，某变量过去的时间值可能会通过非线性函数作用于另一变量当前时刻的值。广义格兰杰因果关系见定义1。

定义1. 在多维时间序列x中，若对所有的$\left(x_{<t1},\cdots ,x_{<tp}\right)$ 和${x^{\prime }}_{<tj}\ne x_{<tj}$ ，有

$g_i\left(x_{<t1},\cdots ,x_{<tj},\cdots ,x_{<tp}\right)=g_i\left(x_{<t1},\cdots ,{x^{\prime }}_{<tj},\cdots ,x_{<tp}\right)$

则称变量j与变量i之间不构成格兰杰因果关系，即$g_i$ 对$x_{<tj}$ 是不变的。若某一对i和j不满足上述条件，则称变量j与变量i之间构成格兰杰因果关系。其中j为原因变量，i为结果变量。$x_{<tj}$ 代表序列j在t时刻之前的序列值。

4. 动态正则格兰杰因果学习模型

本文提出的模型总体结构如图2所示，由两个主要部分组成：线性因果发现网络和采样循环网络。线性网络以多层感知器的形式被构建，每个网络用于单独拟合时间序列数据中所有变量的K阶滞后值到某变量当前时刻值的映射函数；循环网络以长短期记忆网络的形式被构建，每个变量对应的循环网络都需要进行独立的因果图采样并处理其输入的原始序列数据。

两部分模型的训练是交替进行的，在每一个训练轮中，使用线性网络作用于序列数据，得到线性输出预测误差项；将线性网络的输入权重按照滞后阶数进行提取并归一化，使用动态正则策略对权重进行稀疏约束，得到正则项；从线性网络的输入权重中提取出目标变量的单维度因果图，对原始序列数据进行覆盖处理，使用循环网络作用于处理后的序列数据上，得到循环输出预测误差项。通过各个误差项求和得到损失函数并进行优化，使模型逐步收敛至可以精确模拟因果作用机制。

4.1. 线性因果发现网络

本文在线性部分为每个变量单独建立网络模型，如式(3)所示

$x_{ti}=g_i\left(x_{<t1},\cdots ,x_{<tp}\right)+e_{ti}$ (3)

其中，$x_{ti}$ 代表时间序列x_i在t时刻时的值，e_ti代表对应单维序列在t时刻的随机噪声值，g_i是一个函数，指定如何将多维时间序列x的t时刻前的“过去值”映射到序列i上。对多维序列中的每一个变量建立一个MLP，然后将所有的网络并行组合在一起得到输出，如图3(a)所示。

Figure 2. Overall model structure

图2. 总体模型结构

Figure 3. (a) Overall structure of component-level MLP; (b) Single-dimensional MLP structure

图3. (a) 组件级MLP的总体结构；(b) 单维度MLP结构

在单维度线性网络的第一层，有H个大小为$p\times K$ 的一维卷积核，其中H为可设置的隐藏单元维度，p为输入的多维时间序列的维度，K为预设定的最大可能时间滞后值。经过第一层卷积网络后得到维度为H的隐藏层序列，将隐藏序列输入到后续的隐藏线性层中，最后经过一个卷积核大小为$H\times 1$ 的一维卷积层，得到网络的最终输出${x^{\prime }}_i$ ，如图3(b)所示。

本文以上述结构来模拟式(3)中的$g_i(\cdot )$ 函数，在此结构中，网络的参数由每一层的权重W和偏置b确定，其中$W=\left\{W^1,\cdots ,W^L\right\}$ ，$b=\left\{b^1,\cdots ,b^L\right\}$ ，L代表包括第一层线性层、隐藏线性层和最后一层线性层在内的总层数。cMLP的可解释性主要体现在其第一层的网络权重中。第一层的权重结构是一个大小为$H\times p\times K$ 的三维张量，即$W^1\in {ℝ}^{H\times p\times K}$ 。对第一层权重进行分解，$W^1=\left\{W^{11},\cdots ,W^{1K}\right\}$ ，$W^{1k}\in {ℝ}^{p\times H}$ ，$k=1,\cdots ,K$ 。网络的其他参数维度为：$W^l\in {ℝ}^{H\times H}\left(l=2,\cdots ,L-1\right)$ ，$W^L\in {ℝ}^H$ ，$b^l\in {ℝ}^H\left(l=1,\cdots ,L-1\right)$ ，$b^L\in ℝ$ 。基于上述权重表达，输入数据在时间t处经过输入层后的隐藏向量由式(4)给出

$h_t^1=\sigma \left({\displaystyle \sum _{k=1}^K{W}^{1k}{x}_{t-k}}+b^1\right)$ (4)

其中，$\sigma$ 是激活函数，可以是logistic函数或relu函数，后续隐藏层中的隐藏向量使用$h_t^l$ 表示，使用相同的激活函数$\sigma$ ，$h_t^l=\sigma \left(W^l{h}_t^{l-1}+b^l\right)$ 。通过$L-1$ 个隐藏层后，输出的单维度序列$x_{ti}$ 由最后一个隐藏层所有单元的线性组合表示，$x_{ti}=W^L{h}_t^{L-1}+b^L+e_{ti}$ ，其中误差项$e_{ti}$ 由零均值高斯噪声分布建模。

4.2. 采样循环因果发现网络

循环网络能更好地提取时间维度上数据的变化趋势与不同时刻值之间的关系，因此当时间序列数据点有限时，循环网络呈现出比线性网络更好的因果图发现结果。

本文对MLP网络提取出的单维度因果图进行伯努利采样操作，用采样出的因果图对时间序列数据进行覆盖处理，使得输入LSTM网络的数据尽可能只保留对目标变量有因果影响的数据，消除其他变量的影响，如图4(a)所示。得到处理过的数据之后，将其输入对应维度的LSTM网络进行计算，并用输出的结果进行损失函数计算，与MLP网络的结果各自进行反向传播。通过LSTM网络类似“监督”的作用，组合网络会从静态样本数据和与时序相关数据的两个角度对时间序列维度间的因果关系进行学习。

Figure 4. (a) Sampling input LSTM process (b) Weighted hierarchical group lasso penalty

图4. (a) 采样输入LSTM流程；(b) 加权分级group lasso惩罚

单维度因果图的伯努利采样

在MLP网络中，每一个单维度组件模型提取出的因果图大小为$G_i\in {ℝ}^p$ ，其中$i=1,\cdots ,p$ ，多个单维度因果图拼接得到完整因果图G。因果图$G_i$ 中的项$c_{ij}$ 指示序列j对序列i的因果影响程度。使用Gumbel-softmax采样对$G_i$ 进行优化[25]，即

$s_{ij}=\frac{{\text{e}}^{\left(\log \left(c_{ij}\right)+g\right)/\tau }}{{\text{e}}^{\left(\log \left(c_{ij}\right)+g\right)/\tau }+{\text{e}}^{\left(\log \left(1-c_{ij}\right)+g\right)/\tau }}$ (5)

其中，$g=-\log \left(-\log \left(u\right)\right)$ ，$u~\text{Uniform}\left(0,1\right)$ 。$\tau$ 是一个“温度系数”，$\tau$ 值越小，采样就越接近伯努利采样。选择软采样的方式，在进行了softmax求期望操作之后，采样接近1的维度会保留其原始的序列数据，其对目标维度的影响不变；采样接近0的维度序列会被很大程度压缩，其对目标维度的影响在循环网络中也会相应大幅度减少。

使用Gumbel-softmax采样方式实现了采样操作可以反向传递，确保LSTM网络的损失信息可以回传至MLP的输入层权重，对其进行修正。

4.3. 网络损失函数与动态稀疏正则化

如图2所示，本文提出的网络模型在训练过程中的损失函数由三部分组成，分别是MLP网络部分的预测误差、LSTM网络部分的预测误差和MLP网络中输入权重的动态稀疏正则项。其中MLP网络的预测误差用来优化因果图提取网络的权重，使其输入权重快速收敛至正确的值，确保可以在数值上提取到小误差的因果图。LSTM网络的预测误差用于在处理过的输入数据下将循环网络的权重优化至正确的值，使LSTM网络可以通过采样的操作来“修正”全局因果图，更快排除无关维度对目标维度的因果影响。添加MLP网络输入权重的动态稀疏正则项是为了对输入权重的绝对值进行约束，使更多的权重值可以收敛到精确的0值，进而使得真正有因果关系的变量可以获得它们正确的因果关系值。

针对目标变量i，通过最小化式(6)给出的全局损失函数来优化模型

$\underset{W}{\min }{\displaystyle \sum _{t=K}^T{\left(x_{it}-g_i\left(x_{\left(t-K\right):\left(t-1\right)}\right)\right)}^2}+{\displaystyle \sum _{t=K}^T\left(x_{it}-l_i\left(x_{\left(t-K\right):\left(t-1\right)}\right)\right)}+\lambda \rho \left(W^1\right).$ (6)

其中，$g_i(\cdot )$ 代表目标变量i的MLP网络部分的近似函数，$l_i(\cdot )$ 代表目标变量i的LSTM网络部分的近似函数。$\lambda \rho (\cdot )$ 代表施加在MLP网络输入权重上的带系数正则项。在实际的神经网络中，对MLP网络的更高层权重和LSTM网络权重施加2-范数正则化，确保将权重约束在一定范围内。

遵照4.2.1节中的表述，MLP网络的输入权重结构是一个大小为$H\times p\times K$ 的三维张量，$W^1\in {ℝ}^{H\times p\times K}$ 。通过对该张量在0和1维度上求F-范数并进行归一化操作，得到维度平均时滞依赖向量 ${\lambda }_i=\left({\lambda }_{i1},{\lambda }_{i2},\cdots ,{\lambda }_{ik},\cdots ,{\lambda }_{iK}\right)$ ，其中$i=1,\cdots ,p$ 且$k=1,2,\cdots ,K$ 。${\lambda }_{ik}$ 代表第i个变量对过去第k个时间点上的值的平均依赖程度，平均的操作在所有变量的维度上执行。对${\lambda }_i$ 进行分级调整，得到维度平均分级时滞惩罚向量。计算过程由式(7)给出

${{\lambda }^{\prime }}_i=\frac{1}{{\lambda }_i},{{\lambda }^{″}}_i=\frac{{{\lambda }^{\prime }}_i}{{{\lambda }^{\prime }}_{i1}},r_{i1}={{\lambda }^{″}}_{i1},r_{ik}={{\lambda }^{″}}_{ik}-{{\lambda }^{″}}_{i\left(k-1\right)},r_i=\left(r_{i1},\cdots ,r_{iK}\right),k=2,\cdots ,K$ (7)

其中${\lambda }_i$ 是由网络权重求范数得到的维度平均时滞依赖向量，${{\lambda }^{\prime }}_i$ 和${{\lambda }^{″}}_i$ 是中间变量，$r_i$ 是变量i的维度平均分级时滞惩罚向量，用$r_i$ 对MLP网络输入权重值的正则操作进行校正。

动态分级惩罚策略下的$\rho \left(W^1\right)$ 规划由式(8)给出

$\rho \left(W^1\right)={\displaystyle \sum _{j=1}^p{\displaystyle \sum _{k=1}^K{r}_{ik}{\Vert W_j^{1k},\cdots ,W_j^{1K}\Vert }_F}}.$ (8)

其中$W_j^{1k}$ 代表输入变量j对应的输入权重在k时滞的列，$r_{ik}$ 是$r_i$ 的第k项，${\Vert \cdot \Vert }_F$ 是求F-范数操作，如图4(b)所示。

4.4. 优化带正则项的损失函数

本文采用迭代软阈值收缩算法[26]对网络的损失函数进行优化，迭代软阈值收缩算法是近端梯度下降算法的一种特殊形式，其损失函数的形式为预测均方误差加上一个带系数的稀疏正则项。该方法驱使参数矩阵的某些行或列陷入0阈值范围内，即获得精确的0解，这符合本文在3.2节中定义的非线性格兰杰因果关系所追求的目标，即发现并排除不存在格兰杰因果关系的两个变量。

在使用收缩算法进行目标优化时，采用线搜索的方式，保证损失函数可以收敛至局部极小值。网络权重由标准正态分布进行随机初始化为$W^{\left(0\right)}$ ，算法按照式(9)的方式从$W^{\left(0\right)}$ 开始迭代更新网络权重

$W^{\left(n+1\right)}={\text{prox}}_{d^{\left(m\right)}\lambda \rho }\left(W^{\left(n\right)}-d^{\left(n\right)}\nabla ℒ\left(W^{\left(n\right)}\right)\right),$ (9)

其中$W^{\left(n\right)}$ 代表第n个迭代步的网络权重值。$d^{\left(n\right)}$ 代表第n个迭代步的步长值。$ℒ\left(W\right)$ 是网络的预测误差，由线性和循环两部分构成。${\text{prox}}_{\lambda \rho }(\cdot )$ 是关于稀疏惩罚函数$\rho (\cdot )$ 和正则系数$\lambda$ 的近似算子。通过线搜索的方式，算法在每一步中驱动预测损失下降以得到“中间值”，再通过修正中间值得到下一步的数值。

在动态正则化策略中，动态稀疏约束只施加于MLP网络部分的输入权重上，MLP网络的其余层和LSTM网络的所有层权重的迭代步长都是单一的定值，施加在MLP网络输入权重上的动态群组Lasso惩罚的近似步长通过对输入权重进行分层加权软阈值运算得到，运算过程由式(10)给出

${\text{prox}}_{d^{\left(m\right)}\lambda \rho }\left(W_{:k}^1\right)={\left(1-\frac{r_k\lambda d^{\left(m\right)}}{{\Vert W_{:k}^1\Vert }_F}\right)}_{+}{W}_{:k}^1$ (10)

其中，$W_{:k}^1$ 是输入权重中代表对应前k阶滞后值的部分，$r_k$ 是对应每个滞后阶数的平均分级时滞惩罚值。${\left(\theta \right)}_{+}=\max \left(0,\theta \right)$ ，由于对权重取F范数后数值为非负值，所以只需要考虑权重由正值落入阈值范围内的情况。输入权重的近似步长是通过对惩罚函数中每个不同滞后范围的群组迭代应用群组软阈值操作来计算的。

5. 实验

本文针对提出的模型在模拟数据集和真实场景启发的数据集上进行验证，将因果图的学习结果与几种主流的基准方法的实验结果进行对比，以验证本文提出方法的有效性。其次，为模型设置不同的正则化参数来检验其对模型训练效果的影响。

5.1. 数据集

5.1.1. 向量自回归模型VAR数据

VAR数据[9]是通过设定好的因果图、时滞和预热数据模拟生成的可变长度时间序列数据，我们在VAR数据上验证模型在数据因果关系为线性关系时的效果。模拟序列维度为$p=10$ 的线性VAR(2)和VAR(3)数据。需要指出的是，VAR(n)数据集指明了在数据生成时当前时刻数据与过去的n个时刻值相关。如果变量i依赖于变量j，则设置因果图矩阵中对应项$A_{ij}^k=0.1$ ，其中$k=1,2,3$ 。因果图矩阵中的其他项均为零，标识对应变量之间无依赖关系。为了验证DRGC对于不同时间长度序列的有效性以及LSTM网络对时间信息捕捉的能力，分别生成长度为200、500和1000个时间点的序列数据。

5.1.2. Lorenz-96数据

Lorenz-96数据集是一个在混沌理论中广泛使用的合成数据集，用于研究非线性动力系统的复杂行为[27]。该数据集通过一个简化的三维流体动力学模型生成，尽管其维度较低，但却能够展现出高度复杂的混沌动态。p维连续的lorenz-96模型的变化方式由式(11)给出

$\frac{\text{d}{x}_{ti}}{\text{d}t}=\left(x_{t\left(i+1\right)}-x_{t\left(i-2\right)}\right)x_{t\left(i-1\right)}-x_{ti}+F$ (11)

其中$i=1,\cdots ,p$ ，为了使得式(11)对所有的i值都有意义，有式(12)给出的定义。

$x_{t\left(-1\right)}=x_{t\left(p-1\right)},x_{t0}=x_{tp},x_{t\left(p+1\right)}=x_{t1}$ (12)

通过式(12)的定义，变量形成一个循环链。F是一个“强迫常数”，决定了序列非线性的程度和“混乱”的程度。本文模拟$p=10$ 的Lorenz-96数据并采用两种不同的强迫常数，采样率设置为$\Delta t=0.05$ ，得到具有稀疏格兰杰因果关系的多变量非线性时间序列。

5.1.3. DREAM-3数据

我们在真实场景启发的时间基因表达数据集DREAM-3 [28]上验证提出的模型效果。该数据集在格兰杰因果关系检验上是一个较困难的非线性数据集。数据使用连续的基因表达和调控动力学来模拟，背后有多个隐藏因素是未被观察到的。数据共包含五个有着不同真实因果图的模拟数据集。每个数据集都包含10节点数据的4段序列，每段序列均有21个采样时间点，将不同序列前后拼接构成时间长度为84的总序列。

5.2. 基准方法

本文选择了4种主流的基准方法进行比较：1) 神经格兰杰因果发现方法NGC [4]，它利用多层感知机和循环神经网络两种组件模型结合群组惩罚来推断格兰杰因果关系。在本文实验中针对VAR数据集使用MLP模型，针对Lorenz-96数据使用RNN模型。2) PCMCI [29]，一种使用条件独立测试来检测非线性格兰杰因果关系的方法。3) economy-SRU [5]，一种使用基于统计回归单位(Statistical Regression Units, SRU)的分量时间序列预测模型进行非线性建模的方法，通过设计少量的可训练参数，提高了模型对预测数据的抗过拟合性能。4) CUTS [22]，一种通过交替进行因果发现和潜在数据预测两个阶段来学习数据中的因果关系的方法，这种方法可以对不规则数据(存在缺失值的时间序列)同时进行因果发现和数据填充。它同样适用于不存在缺失值的时间序列数据集上，并且反映出较好的学习结果。

5.3. 模拟数据实验结果

在定量评估方面，为了验证因果图复原数据底层结构的准确性，实验中以ROC曲线下面积指标(Area Under the Receiver Operating Characteristic Curve, AUROC)作为标准。具体来说，对本文提出的方法和NGC方法，通过将一定范围内变化的$\lambda$ 值作为动态正则项的系数进行模型训练，得到模型对数据集上无阈值因果图学习的ROC曲线，并计算AUROC值；对PCMCI和eSRU方法则是使用不同的阈值来获取AUROC值。

5.3.1. 因果图学习的AUROC比较

本部分实验生成VAR(1)、VAR(2)、VAR(3)数据和使用两种不同“强迫常数”的Lorenz-96数据，每个数据集都有200、500和1000三个不同时间长度的版本，以观察在数据集长度变化的情况下各个方法的学习效果，如表1所示。表中的DRGC-s方法代表在本文提出的模型基础上除去循环网络模块的消融模型。表1中所呈现的数值(如99.86 ± 0.12)均为AUROC的数值表示，其完整数值为对应的0.9986 ± 0.0012 (若换算为小数形式)或者98.66% ± 0.12% (若换算为百分比形式)。

Table 1. AUROC for causal graph using different methods on 4 time series datasets

表1. 不同方法在4个时间序列数据上因果图学习的AUROC对比

Model	VAR(2)			VAR(3)
T	200	500	1000	200	500	1000
DRGC	99.86 ± 0.12	100.00 ± 0.00	100.00 ± 0.00	98.86 ± 1.42	100.00 ± 0.00	100.00 ± 0.00
DRGC-s	99.17 ± 0.98	100.00 ± 0.00	100.00 ± 0.00	97.56 ± 2.05	100.00 ± 0.00	100.00 ± 0.00
NGC	96.01 ± 2.14	91.98 ± 3.87	98.52 ± 1.31	90.91 ± 4.15	97.53 ± 2.82	98.74 ± 0.98
PCMCI	71.40 ± 4.93	72.30 ± 5.73	72.19 ± 3.25	65.89 ± 5.23	71.30 ± 4.85	71.92 ± 3.61
eSRU	87.74 ± 6.71	89.97 ± 5.72	91.06 ± 4.86	81.25 ± 7.05	86.55 ± 5.48	90.38 ± 2.29
CUTS	98.92 ± 0.84	99.05 ± 1.03	100.00 ± 0.00	98.05 ± 1.64	99.91 ± 0.03	99.99 ± 0.02
Model	Lorenz-96 (F= 10)			Lorenz-96 (F = 20)
T	200	500	1000	200	500	1000
DRGC	94.91 ± 0.96	100.00 ± 0.00	100.00 ± 0.00	92.00 ± 0.97	95.49 ± 2.30	98.80 ± 0.45
DRGC-s	93.63 ± 0.55	98.64 ± 0.86	100.00 ± 0.00	90.51 ± 1.48	95.66 ± 1.62	98.01 ± 1.57
NGC	93.88 ± 1.64	98.56 ± 0.51	99.15 ± 0.33	84.35 ± 3.31	92.22 ± 3.72	92.82 ± 2.60
PCMCI	85.95 ± 4.55	91.65 ± 3.54	95.72 ± 2.23	80.87 ± 5.21	86.56 ± 3.18	86.39 ± 2.71
eSRU	87.42 ± 3.98	96.45 ± 3.06	97.51 ± 1.16	90.30 ± 3.63	97.57 ± 1.43	98.43 ± 1.74
CUTS	95.12 ± 1.13	99.64 ± 0.51	100.00 ± 0.00	89.20 ± 0.73	93.36 ± 0.82	97.15 ± 2.94

实验结果表明，DRGC在不同长度的VAR数据和Lorenz-96数据上基本取得了最优的结果。当时间序列的长度越长时，提供给模型的训练样本越多，所以基本所有方法在长度为500和1000数据上的结果都优于长度为200的数据。此外，由于VAR数据结构较为简单，DRGC可以完全准确地识别出系统中每个变量的真实原因变量，并在混乱度较低的Lorenz数据上保持优秀的识别率。相比于NGC方法，DRGC方法在长序列上大致有1%至5%的指标提升，这是动态正则化的作用。滞后值的正确选择使得模型能够准确量化目标变量对原因变量的依赖程度，进而更为精确地排除无关变量的微小影响，事实上，NGC方法有时选择出的错误因果关系相比于正确的那些在任何滞后值上都只有很小的权重，而传统的正则方法无法很好地排除它们的影响。

5.3.2. 不同依赖系数的VAR(3)数据上的AUROC比较

VAR数据的依赖系数p代表了在这个模拟系统中每个变量受几个变量的驱动，$p=0.2$ 代表驱动目标变量变化的“原因”变量个数占系统中总变量个数的20%，同时每个目标变量的“原因”变量一定包含其本身。p值越大，系统中相互依赖的变量数量就越多，对应的因果图边的数量也越多，系统变得复杂，模型发现正确的因果图则更加困难。本部分实验生成依赖系数分别为0.2、0.3和0.4的三种VAR(3)数据，且同样有三种不同的时间长度，以验证模型在面对不同复杂程度的数据集时发现因果图的能力。学习结果如表2所示，表中数据同样是百分比数值。

Table 2. AUROC for causal graph on VAR(3) with three dependency coefficients p

表2. 不同方法在3种依赖系数p下的VAR(3)数据上因果图学习的AUROC对比

Model	VAR(3)
	p = 0.2			p = 0.3			p = 0.4
T	200	500	1000	200	500	1000	200	500	1000
DRGC	98.86 ± 1.42	100.0 ± 0.00	100.0 ± 0.00	86.57 ± 6.07	90.63 ± 5.08	94.14 ± 5.34	78.67 ± 7.13	85.82 ± 2.27	88.21 ± 6.77
DRGC-s	97.56 ± 2.05	100.0 ± 0.00	100.0 ± 0.00	85.64 ± 4.09	99.43 ± 0.75	100.0 ± 0.00	76.06 ± 5.96	91.99 ± 1.38	95.72 ± 0.51
NGC	90.91 ± 4.15	97.53 ± 2.82	98.74 ± 0.98	81.18 ± 3.16	90.38 ± 3.91	93.84 ± 3.27	75.85 ± 5.02	82.84 ± 4.55	87.82 ± 3.67
PCMCI	65.89 ± 5.23	71.30 ± 4.85	71.92 ± 3.61	57.52 ± 4.83	54.95 ± 6.50	54.91 ± 5.81	56.49 ± 7.42	51.95 ± 7.75	53.49 ± 7.78
eSRU	81.23 ± 7.16	86.22 ± 5.49	90.10 ± 2.66	70.75 ± 7.33	75.68 ± 6.73	82.26 ± 4.11	63.95 ± 7.48	68.60 ± 5.13	73.99 ± 8.05
CUTS	98.97 ± 1.64	99.91 ± 0.03	99.99 ± 0.02	82.58 ± 4.79	95.33 ± 0.99	94.45 ± 2.28	75.91 ± 6.73	83.43 ± 1.65	92.59 ± 3.31

实验结果表明，随着依赖系数增大，时间序列系统各变量间的因果关系更为复杂，各个方法对因果关系学习的准确程度有所下降。当系统中每个变量仅有两个依赖变量时，DRGC方法可以准确预测出正确的因果图。对于较复杂的系统因果图，DRGC方法仍能领先其他方法表达出较好的效果。序列长度为500和1000时，得益于足够多的训练样本，消融模型学习到了最好的结果；而带有LSTM网络的完整模型则能够在时间序列长度较短时更好地捕捉时间维度上的依赖关系，表现出更好的性能。

5.3.3. 在VAR(3)数据上与NGC的时滞选择比较

VAR(3)数据在模拟生成时，每一时刻某变量的数据值仅受其原因变量过去最多三个时刻的值的影响，即$x_t$ 与$x_{t-4}$ 无因果关系。DRGC方法中的线性网络动态正则策略的最大优势在于可以定量还原不同滞后值之间的相对因果关系大小，同时不影响对无关变量的排除。在生成VAR数据时，本文设置当前值对过去3个滞后值的依赖程度相同，体现在因果图上则是具有相同的权重值。为了检验模型排除高阶滞后值的能力，在线性网络一维卷积时使用的卷积核大小为5。选取10维中的第2个和第5个变量为例，观察其时滞因果图的学习结果，如图5所示。第2变量的原因变量是第2和第5变量，第5变量的原因变量是第5和第7变量，两种方法学习到的单变量对单变量的时滞差异平均值的对比如表3所示。

Table 3. Mean value of the difference in the delay of a single variable on VAR(3)

表3. 在VAR(3)上单一变量时滞差异平均值对比

Model	变量2-变量2	变量2-变量5	变量5-变量5	变量5-变量7
cMLP	0.47255	0.30196	0.43922	0.46863
DRGC	0.09412	0.06564	0.11765	0.19216

Figure 5. Lagged causal graph on VAR(3) of DRGC and cMLP

图5. VAR(3)上的滞后因果图学习结果

由图5可得，cMLP的学习结果由于受到硬性的分级递增惩罚策略的限制，在选取的两个变量上的因果关系大小均随着滞后值的增大而有非常明显的减弱。同时当预设的最大滞后阶数增加时，惩罚的变化较为缓慢，不能较好地排除无关滞后值的影响。DRGC方法学习到了目标变量在不同滞后值上相对差异较小的因果关系大小，这符合VAR(3)数据在模拟生成时“过去3个时刻值具有相同大小的因果关系”的正确机制。

从具体数值的差异方面来看，DRGC方法对各个变量时滞权重之间的差异相比cMLP缩减了20%至40%。表3中所示的差异平均值计算方法是对两个单变量之间的三阶时滞值作差，然后对两个差值求平均，所有的时滞值在同一个目标变量下进行了归一化处理。DRGC方法将时滞因果强度之间的差异缩小到了10%，而cMLP方法的平均差异处于40%至50%之间，相当于最大时滞因果强度的一半。由于cMLP采用按时滞递增的惩罚方法，相邻时滞之间的惩罚差异相同，会导致学习到的时滞在一定程度上呈现出均匀下降的趋势。DRGC方法实现了最后一阶时滞与无因果时滞之间惩罚值的跳变，可以更加合理地拟合时滞的下降趋势。

5.3.4. 不同正则系数下VAR(3)数据的时滞选择

DRGC方法在线性层的输入权重上施加了动态的分级稀疏惩罚。对于对应各组滞后的输入权重，除了需要加入其维度平均分级时滞惩罚的值外，还需要为所有的惩罚值设置一个统一的系数$\lambda$ ，以实现神经网络对各个变量和各个滞后值对应权重的整体调控。当$\lambda$ 的值变化时，模型对数据的因果图提取结果也会变化。

本部分实验使用与5.3.3节中相同设置的VAR数据集，采用多种不同的动态正则系数$\lambda$ 进行时滞选择实验并与真实的滞后因果图进行对比，我们选择了其中一个对比较为明显的变量进行展示。

图6展示了一系列动态正则系数下模型对目标变量包含滞后因果图的学习结果和其真实的因果图，我们训练一个具有分级Lasso惩罚和最大滞后阶数为5的cMLP模型作为对比。

Figure 6. Lagged causal graph of DRGC and cMLP with different dynamic regularization coefficients

图6. 不同动态正则系数下DRGC和cMLP的时滞因果图学习结果

由图6可得，当$\lambda$ 取值适当时，NGC和DRGC方法都可以很好地发现目标变量正确的因果图，且DRGC相较于NGC可以更合理地分配不同滞后值的因果大小。对于较小的$\lambda$ 值，NGC方法会过多估计原因变量的滞后阶数；而DRGC方法则可以更多地保持住不同滞后值之间的相对平稳，在无关的变量上的估计值相对很小，可以通过设置相应的阈值来排除影响。对于较大的$\lambda$ 值，NGC方法会过多削弱更靠后的滞后值的大小，而DRGC方法仍可以较好地保持各滞后值的相对大小关系，同时更多地排除无关变量的影响。当$\lambda$ 值设置得过高时，由于输入权重惩罚值与预测误差之间的平衡，模型将学习不到目标变量的因果图，权重值被惩罚至0。实验中可以通过交叉验证来选择合适的$\lambda$ 值。

5.4. DREAM-3实验

DREAM-3数据集是用于严格比较格兰杰因果关系检测方法的较困难的非线性数据集，包含三个Yeast数据和两个E. Coli数据。我们分别在这5个数据集上使用DRGC方法进行学习，同时选用NGC方法的线性模型和循环模型分别作为对比。由于DREAM-3数据的序列长度较短，且数据中只有真实的全局因果图，我们将模型中的最大滞后阶数设置为2，隐藏层单元个数设置为10，以减轻模型的体量并加快训练过程，并对NGC方法采用相同的设置。在指标方面，仍使用AUROC来评估因果图的学习结果，在5个时间序列数据集上的评估结果如表4所示，对应的ROC曲线如图7所示。

Table 4. ROC for DRGC, cMLP, and cLSTM models on DREAM-3 datasets

表4. DRGC，cMLP和cLSTM在5个Dream数据集上的AUROC对比

Model	Yeast1	Yeast2	Yeast3	Ecoli1	Ecoli2
cMLP	0.5722	0.5629	0.5553	0.5741	0.5855
cLSTM	0.5933	0.58	0.5525	0.5132	0.5647
DRGC	0.6568	0.6485	0.6067	0.6236	0.6070

Figure 7. ROC for DRGC, cMLP, and cLSTM models on DREAM-3 datasets

图7. DRGC，cMLP和cLSTM在Dream数据集上的ROC曲线

由表4可得，DRGC方法在所有五个数据集上的结果优于cMLP和cLSTM两个模型，由于10节点的Dream数据集只包含4个长度为21的时间序列，用于训练的样本数量较少，模型学习到正确因果图的难度较高，各个模型的AUROC值均较低，但DRGC方法达到了60%的指标。DRGC方法融合了线性网络与循环网络的优点，可以使权重值向对应真实原因变量的单元集中。

6. 结论

本文提出了一种采用动态分级稀疏惩罚策略的线性与循环网络组合多维时间序列非线性格兰杰因果挖掘方法。为提高模型的可解释性，我们为系统中的每个变量单独建立线性网络和采样输入循环网络，依靠循环网络提取时间方向依赖关系的能力，对线性网络挖掘的整体因果关系进行监督和修正。为了提高模型在不同时滞上选择因果关系的准确性，通过从线性网络输入权重中提取变量对不同时滞的平均依赖程度，对模型施加动态变化的稀疏惩罚。在模拟数据集和真实基因调控子网络生成的数据集上进行的实验结果表明，DRGC方法在不同长度和复杂程度的数据集上都表现出领先于其他现有方法的性能。

未来的工作方向有两点，1) 探究模型如何在脱离数据集的真实全局因果关系和时间因果关系信息的情况下，配合因果图的阈值设置来进行更合理的参数调整与选取；2) 在更多的真实世界数据集上进行实验，探寻合适的底层系统结构来协助进行因果关系验证，以提高方法的可推广性。

基金项目

国家自然科学基金(62262016)；

中国澳门基金会2024年学术资助计划“建设横琴数据跨境传输安全管理试点(数据飞地)的多视角可行政策研究”(G01156-2309-262)；

中央高校基本科研业务费专项资金(2023JBZY035)资助。

NOTES

^*通讯作者。