GRACE：面向高维时间序列精准因果边发现的门控精化方法

# GRACE：高维时间序列中精确因果边发现的门控精炼方法
来源：https://arxiv.org/html/2606.23880

###### 摘要

从气候遥相关到基因调控，现代时间序列数据集包含数十或数百个交互变量，使得因果发现日益具有挑战性。基于约束的方法提供了统计学严谨性，但其非线性条件独立检验在大规模场景下不可行，而基于评分的方法避免了条件独立检验，但需要人为设定阈值来将连续的边评分二值化。我们提出GRACE（基于门控精炼的精确因果边发现方法），该方法利用带有L0正则化的Hard Concrete门控机制来精炼基于约束的发现结果——每个候选边拥有一个独立的门控，其值集中在0或1附近，产生清晰的二模分离，使得二值决策鲁棒——这与L1和基于注意力的方法所产生的狭窄、重叠的评分分布截然不同。一个快速的线性条件独立骨架提供了高召回率的候选边；然后，一个单一的门控模型通过学习哪些边能真正改善预测来剪除误报，其正则化自动适应问题维度和骨架密度。在涵盖多种图拓扑结构（无标度、Erdős–Rényi、小世界）和高达d=100维度的合成基准上的系统实验表明，GRACE在保持高精度的同时，显著提升了其基础条件独立方法的F1分数，并且优于基于注意力和基于评分的方法。GRACE以极低的成本（快75倍）达到或超过了昂贵的非线性条件独立检验。在一个真实世界的河流流量数据集上（其中降雨混杂因子、变量传播延迟和分布偏移违反了标准假设），GRACE的时间自举变体成功恢复了易北河上11条因果边中的9条，仅产生1个误报（F1=0.86，AUROC=0.99），将原始骨架的106个误报减少了99%。

## 1 引言

从观测时间序列中推断因果关系是气候科学、神经科学、金融学、基因调控网络重建、工业根因分析以及许多其他领域的核心挑战 (Peters et al., 2017)。给定一个包含T次观测的多变量时间序列xt = (xt^1, ..., xt^d)，目标是恢复*滞后因果图* G，其中 G[c, ℓ, e] = 1 表示变量 x^c 在滞后 ℓ 处是变量 x^e 在时间 t 的直接原因。

时间序列因果发现主要有两类方法，各有其优缺点。*基于评分的方法*——神经格兰杰因果关系 (Tank et al., 2021)、NAVAR (Bussmann et al., 2021)、DYNOTEARS (Pamfil et al., 2020)、CUTS+ (Cheng et al., 2024)——采用全局优化视角，搜索在某个评分准则下最拟合数据的图。它们可以捕捉非线性依赖关系，并能扩展到中等维度，但通过连续惩罚（L1或Gumbel-Softmax）来强制稀疏性，需要事后对权重进行阈值化才能得到二值图——而模型设定错误或评分假设不佳会悄然降低恢复图的质量。基于注意力的方法，如 TCDF (Nauta et al., 2019)，同样优化预测目标，但将注意力权重用作因果代理；然而，softmax 归一化在候选父节点间造成人为竞争，且注意力阈值需要针对数据集进行调整。*基于约束的方法*——PCMCI/PCMCI+ (Runge et al., 2019; Runge, 2020)、CDNOTS (Sadeghi et al., 2025)、SyPI+ (Fesanghary and Gopal, 2025)、LPCMCI (Gerhardus and Runge, 2020)——通过条件独立检验构建因果图，提供了可解释的、统计原则性的推断，无需全局评分函数。然而，这些检验对样本量和条件集维度敏感：非线性条件独立检验如 CMI (Runge, 2018)、KCIT (Zhang et al., 2011)（复杂度 O(n^3)）和 RCoT (Strobl et al., 2019)（对于 df 个傅里叶特征，复杂度 O(df^2 n)）在高维 d 时仍然成本高昂，迫使实践者依赖快速但会遗漏非线性依赖关系的线性检验。更广泛的综述可参见 Assaad et al. (2022) 和 Gong et al. (2023)。

在这项工作中，我们提出 GRACE（基于门控精炼的精确因果边发现），一个用于改进高维约束因果发现的框架。GRACE 分两阶段运行：首先，一个高召回率的约束方法（如 CDNOTS (Sadeghi et al., 2025) 或 PCMCI (Runge et al., 2019)）生成一个候选骨架；然后，一个门控神经模型通过为每个候选边分配一个独立的、经过 L0 正则化训练的 Hard Concrete 门控 (Louizos et al., 2018; Maddison et al., 2017) 来精炼此骨架，其门控值集中在 0 或 1 附近，从而在自然的 0.5 阈值下产生鲁棒的二值图。门控模型学习哪些骨架边能真正改善预测，在剪除误报的同时保留真实的因果联系——以极低的计算成本缩小了与昂贵非线性条件独立检验的差距。尽管我们的实验集中于滞后效应（τ ≥ 1），该框架原生支持同期（滞后 0）边：启用滞后 0 门控允许相同的 L0 机制从骨架提供的瞬时候选中进行选择。

我们的贡献如下：

1.  1. GRACE框架：一个门控精炼阶段，通过带有 L0 正则化的 Hard Concrete 门控将高召回率的骨架转换为高精度的因果图。门控机制独立于骨架来源：我们使用 CDNOTS 和 PCMCI 骨架进行了演示。
2.  2. 一个经验推导的正则化规则 (Eq. 12)，使 λ 适应问题维度（d, T）和骨架密度，并结合 L0 归一化，将正则化与批量大小解耦。
3.  3. 在易北河网络 (Stein et al., 2025) 上的真实世界评估，其中一种时间自举变体——将 GRACE 应用于随机时间窗口并保留一致出现的边——无需领域特定启发式方法即可处理非平稳性、隐藏混杂因子和可变滞后。

## 2 方法

### 2.1 问题定义

考虑一个 d 维平稳时间序列 {xt}_t=1^T，其中 xt ∈ ℝ^d。我们假设一个具有最大滞后 L 的结构因果模型：

x_t^e = f_e( { x_{t-ℓ}^c : G[c, ℓ, e] = 1 } ) + ε_t^e, ε_t^e ∼ N(0, σ_e^2), (1)

其中 G ∈ {0,1}^{d × (L+1) × d} 是滞后因果图，f_e 是效果变量 e 的（可能非线性的）因果机制。目标是从观测 {xt}_t=1^T 中恢复 G。

门控精炼阶段假设：(i) *因果充分性*：不存在未观测到的共同原因；(ii) *有限滞后*：所有因果效应发生在滞后 0, ..., L 内；(iii) *加性结构*：效应在父节点间加性组合，尽管每个个体效应 f_{c,ℓ → e} 可以是非线性的。其他假设（例如平稳性）继承自骨架发现算法；当使用 CDNOTS 时，非平稳性被原生处理。

基于约束方法（例如 CDNOTS, PCMCI）骨架 S 掩码 Stage 1 Stage 2 x_{t-1}^1 x_{t-1}^2 ⋮ x_{t-L}^d W_{1,1} W_{2,1} W_{d,L} ⋮ z_1 z_2 z_D × × × ⋮ Σ MLP_e μ_e, σ_e 输入 编码器 门控 解码器 共享 每个效果 e Hard Concrete + L0 惩罚
图 1: GRACE 两阶段流水线。Stage 1: 基于约束的方法（例如 CDNOTS 或 PCMCI）识别候选边。Stage 2: 门控模型精炼骨架——只有 S 中的边具有可学习的 Hard Concrete 门控；所有其他边被掩码为零。L0 惩罚驱使门控精确趋近于零或一，产生双峰分布，使 0.5 决策边界鲁棒。

### 2.2 门控因果模型

我们为每个效果变量 e 构建一个预测模型，其中候选因果输入被单独门控。该架构由三个组件组成（图 1）：共享编码器、每个效果的门控聚合、每个效果的解码器。

#### 共享编码器。

一个线性投影将每个滞后输入映射到一个所有效果变量共享的隐藏表示：

h_{c,ℓ} = W_{c,ℓ} x_{t-ℓ}^c + b_{c,ℓ} ∈ ℝ^H, (2)

其中 H 是编码器隐藏维度，W_{c,ℓ} ∈ ℝ^{H × 1}，b_{c,ℓ} ∈ ℝ^H。跨效果共享编码器将参数从 O(d^2 L H) 减少到 O(d L H)；效果特定性由门控提供。

#### 门控聚合。

对于每个效果 e，一个 Hard Concrete 门控 z_{c,ℓ,e} ∈ {0,1} 独立控制每个候选父节点：

\hat{x}_e = ∑_{c,ℓ} z_{c,ℓ,e} · h_{c,ℓ}. (3)

门控由可学习的 log-odds log α_{c,ℓ,e} 参数化（第 2.3 节）。

#### 每个效果的解码器。

每个效果变量都有一个独立的两层 MLP 解码器，带有 SiLU 激活函数，将聚合表示映射到高斯参数：

(μ_e, log σ_e) = MLP_e(\hat{x}_e), MLP_e: ℝ^H → ℝ^2. (4)

使用每个效果的解码器（而不是共享解码器）可以防止混淆不同的因果机制——例如，二次效应 x_{t-1}^c → x_t^e 和线性效应 x_{t-1}^c → x_t^{e'} 可以学习独立的响应函数。

### 2.3 Hard Concrete 门控

我们采用 Hard Concrete 分布 (Louizos et al., 2018)，它通过拉伸二元 concrete (Gumbel-Softmax) (Maddison et al., 2017; Jang et al., 2017) 分布来产生精确的零和一。

#### 训练（随机）。

在训练期间，通过重参数化技巧采样门控：

u ∼ Uniform(ε, 1-ε), (5)
s = sigmoid( (log u - log(1-u) + log α) / τ ), (6)
\bar{z} = s · (ζ - γ) + γ, (7)
z = min( max( \bar{z}, 0 ), 1 ), (8)

其中 τ = 2/3 是温度，γ = -0.1，ζ = 1.1 定义了拉伸区间。拉伸超过 [0,1] 然后进行裁剪，使得精确 0 和精确 1 上放置了正概率质量。

#### 评估（确定性）。

在测试时，我们使用 Louizos et al. (2018) 中的确定性 Hard Concrete 估计量：

z = min( max( sigmoid(log α) · (ζ - γ) + γ, 0 ), 1 ). (9)

至关重要的是，拉伸和裁剪机制在*精确*零上放置了正概率质量，因此训练后的门控收敛到 z=0（边不存在）或 z≈1（边存在）；附录 A 通过门控值直方图经验证实了这一点，显示真边和假边之间存在清晰的双模分离。最终因果图通过在 0.5 处对确定性门控值进行阈值化获得——这是一个自然边界，落在两个模式之间的间隙中，使得二值决策对确切的阈值选择鲁棒。与 L1 和基于注意力的方法所需的*事后权重阈值化*不同——这些方法中连续边得分聚集在一个狭窄范围内，结果对所选截止值敏感——双模门控分布使 0.5 边界成为一个原则性默认值，而非任意调谐参数。门控初始化（log α = -0.5）和 L0 惩罚强度（λ）仍然是设计选择，但在训练前一次性设置，而不是在检查学习到的权重后进行调整。

#### L0 惩罚。

门控非零的概率具有闭式表达式：

P(z ≠ 0) = sigmoid( log α - τ log(-γ/ζ) ). (10)

这实现了无需蒙特卡洛估计的分析 L0 正则化器。

#### 初始化。

我们将所有门控的 log α 初始化为 -0.5，产生大约 0.36 的确定性门控值（低于活跃阈值 0.5）。门控初始大多关闭，要求边展示预测价值才能打开。与在 log α ≈ 0 附近的对称初始化相比，这种非对称初始化减少了误报。

### 2.4 损失函数

训练目标结合了预测质量和稀疏性：

L = 预测项 [ 1/(Bd) Σ_{t=1}^{B} Σ_{e=1}^{d} NLL(x_t^e | μ_t^e, σ_t^e) ] + L0 稀疏项 [ λ · s Σ_{c,ℓ,e} P(z_{c,ℓ,e} ≠ 0) ]
     + 滞后集中项（可选） [ λ · λ_{lag} Σ_{c,e} ReLU( Σ_ℓ P(z_{c,ℓ,e} ≠ 0) - k ) ], (11)

其中 NLL(x | μ, σ) = (x-μ)^2/(2σ^2) + log σ + (1/2) log 2π 是高斯负对数似然，B 是批量大小，d 是变量数量。