基于SURE调优岭回归的等方差线性高斯DAG因果发现

# 通过 SURE 调谐岭回归在等方差线性高斯有向无环图中进行因果发现
来源：https://arxiv.org/html/2608.17132
## 通过 SURE 调谐岭回归在等方差线性高斯有向无环图中进行因果发现本文工作得到了以下一项或所有基金的资助：ARO W911NF1910269, ARO W911NF2410094, ONR N00014-22-1-2363, NSF CIF-2311653, NSF CIF-2148313, NSF RINGS-2148313, NSF DBI-2412522，并且还部分得到了 RINGS 项目中指定的联邦机构和行业合作伙伴资金的支持。

Sambit Mishra & Urbashi MitraAffiliation:Department of Electrical and Computer Engineering University of Southern California \{sambitmi, ubli\}@usc\.edu

###### 摘要

从观测数据中恢复结构方程模型（SEM）的有向无环图（DAG）是因果发现中的一个核心问题。连续优化方法的迭代梯度下降和逐问题超参数调优不太适合两个实践中重要的场景：样本受限场景（其中样本数量与 DAG 中节点数量相当或更小）和计算受限场景。本文提出 SURE-Ridge，这是一种针对等方差线性高斯 SEM 的非迭代、闭式估计器。该方法执行并行的逐节点回归，其正则化参数由 Stein 无偏风险估计（SURE）自适应选择，并对所得的软邻接矩阵应用自适应阈值过程以提取 DAG。数值结果表明，与 NOTEARS、DAGMA 和 GBNSL 基线相比，SURE-Ridge 在小样本场景下实现了最低的结构汉明距离，并在所有测试样本量下实现了最低的运行时间。

###### 索引术语：

因果发现，结构方程模型，有向无环图，Stein 无偏风险估计，岭回归

## I引言

线性结构方程模型（SEMs）为表示随机变量之间的因果关系提供了一个灵活的框架，而从观测数据中恢复底层有向无环图（DAG）是因果发现中的一个核心问题，其应用涵盖流行病学[12 (https://arxiv.org/html/2608.17132#bib.bib7)]、无线网络[9 (https://arxiv.org/html/2608.17132#bib.bib8)]和金融欺诈检测[16 (https://arxiv.org/html/2608.17132#bib.bib9)]。对于所有节点噪声方差相等的线性高斯 SEM，真实 DAG 仅从联合分布即可识别[11 (https://arxiv.org/html/2608.17132#bib.bib1)]。我们在此考虑两个关键场景：(a) 样本受限场景（如[13 (https://arxiv.org/html/2608.17132#bib.bib3), 8 (https://arxiv.org/html/2608.17132#bib.bib4)]中所述），其中样本数量通常与节点数量相当甚至更小，如[10 (https://arxiv.org/html/2608.17132#bib.bib13), 19 (https://arxiv.org/html/2608.17132#bib.bib14)]；(b) 计算受限场景，其中必须像在线和流式应用[17 (https://arxiv.org/html/2608.17132#bib.bib15)]中那样持续估计 DAG。

经典方法如 PC[14 (https://arxiv.org/html/2608.17132#bib.bib10)] 和 GES[5 (https://arxiv.org/html/2608.17132#bib.bib11)] 的计算成本随 DAG 中节点数量的增加而急剧增长。NOTEARS[20 (https://arxiv.org/html/2608.17132#bib.bib6)] 和 DAGMA[2 (https://arxiv.org/html/2608.17132#bib.bib12)] 引入的连续优化表述，通过可微的无环性泛函将结构学习重新表述为平滑的约束优化问题，但它们需要数千步梯度下降和逐问题超参数调优，因此不太适合上述考虑的两个场景。

直接利用等方差假设的工作主要考虑两种原则。第一种是从逆协方差矩阵进行迭代剥离：GBNSL[6 (https://arxiv.org/html/2608.17132#bib.bib16)] 逐个识别终端节点，而 BUILD[1 (https://arxiv.org/html/2608.17132#bib.bib17)] 自底向上修剪叶节点。第二种方法是分阶段回归，将排序和父节点选择分为顺序阶段：首先按条件方差对节点排序，然后在恢复的顺序上拟合父节点[4 (https://arxiv.org/html/2608.17132#bib.bib18)]。与 SURE-Ridge 不同，这两种策略都将排序恢复视为一个明确的中间步骤，并且当样本量与 DAG 中的节点数量相当时，会在各阶段累积有限样本误差。

我们提出 SURE-Ridge，这是一种针对等方差线性高斯 DAG 恢复的非迭代、闭式估计器。该方法执行并行的逐节点岭回归，其 L2 正则化参数通过最小化 Stein 无偏风险估计（SURE）[15 (https://arxiv.org/html/2608.17132#bib.bib2)] 来选择。SURE 提供了一种可处理的预测误差估计，产生了一种既不需要保留数据也不需要迭代搜索的调优规则，因此在样本稀缺时相对于交叉验证，以及在需要估计许多 DAG 时相对于基于梯度的调优，都具有决定性优势。将 SURE 与岭回归配对在等方差假设下是自然的，其中每个节点的结构方程本身就是其父节点上的线性高斯回归，并且子问题解耦，从而能够并行发现。本文的贡献如下：
1. 一个通过 SURE 调谐的逐节点回归恢复等方差线性高斯 DAG 的闭式表述，每个 DAG 的总计算成本为 O(d^4 + d^3 n + d M(d-1) + d^3 log d)；
2. 一个自适应阈值过程，仅需要两个尺度不变常数，并在典型情况下产生有效的 DAG；
3. 数值结果表明，与 NOTEARS、DAGMA 和 GBNSL 基线相比，SURE-Ridge 在小样本场景下实现了最低的结构汉明距离（SHD），并在所有测试样本量下实现了最低的运行时间。

本文组织如下。第 II 节 (https://arxiv.org/html/2608.17132#S2) 描述系统模型和符号。第 III 节 (https://arxiv.org/html/2608.17132#S3) 推导 SURE 调谐岭估计器。第 IV 节 (https://arxiv.org/html/2608.17132#S4) 介绍自适应阈值过程。第 V 节 (https://arxiv.org/html/2608.17132#S5) 报告数值结果，第 VI 节 (https://arxiv.org/html/2608.17132#S6) 总结工作。

## II系统模型

我们考虑一个 d 节点因果 DAG G = (V, E)，其中 V = {1, ..., d} 表示节点集，E = {(i,j): 存在边 i→j} 表示边集。我们使用加权邻接矩阵 W ∈ R^(d×d) 表示图，遵循约定
W_{i,j} ≠ 0 ⇔ (j,i) ∈ E。(1)
对于每个节点 i ∈ V，我们定义因果系数行向量 w_i ∈ R^{1×(d-1)}（W 的第 i 行，移除了 W_{i,i}），
w_i ≜ W_{i,-i}。(2)
令 X_i 为与节点 i 相关的随机变量，x_i ∈ R^{1×n} 为 X_i 的 n 个独立同分布观测的行向量。观测的连接记为 X ∈ R^{d×n}，其中 X_{i,:} = x_i。在线性高斯 SEM 中，噪声方差相等，有
x_i = w_i X_{-i,:} + n_i,(3)
其中 X_{-i,:} 表示移除了第 i 行的数据矩阵 X，且 n_i ∼ N(0, σ² I_n)，并且对于 i 独立同分布。根据[11 (https://arxiv.org/html/2608.17132#bib.bib1)]，这种情况是可识别的。

我们看到 (3) 可以解释为一个回归，以从由因果系数向量 w_i 控制的带噪观测中确定 X_i。因此，我们可以使用岭回归来估计 w_i，其中包含可调正则化参数 λ_i。令对应于节点 i 的格拉姆矩阵为 G_i = X_{-i,:} X_{-i,:}^T，则 w_i 的基于岭回归的估计器由下式给出
\hat{w}_i(λ_i) = (G_i + λ_i I_{d-1})^{-1} X_{-i,:} x_i^T。(4)
可调参数 λ_i 的一个合理选择至关重要，因为它直接控制所得估计器 \hat{w}_i(λ_i) 的偏差-方差权衡。我们采用 SURE 作为选择 λ_i 的数据驱动准则。以下命题形式化了逐节点的 SURE 目标。

###### 命题 1（逐节点 SURE 目标）。

令噪声向量满足 n_i^T ∼ N(0, σ² I_n)，令 \hat{w}_i(λ_i) 为 (4) 中定义的岭估计器。那么，期望平方误差风险 E[||x_i^T - X_{-i,:}^T \hat{w}_i^T(λ_i)||_2^2] - nσ² 的一个无偏估计由下式给出
SURE_i(λ_i) = -nσ²(5)
+ ||x_i^T - X_{-i,:}^T (G_i + λ_i I_{d-1})^{-1} X_{-i,:} x_i^T||_2^2
+ 2σ² tr[(G_i + λ_i I_{d-1})^{-1} G_i]。
该命题源自将[15 (https://arxiv.org/html/2608.17132#bib.bib2)]中的 SURE 构造应用于 (3) 中的回归模型。因此，(5) 中的目标函数为我们选择每个节点 i ∈ V 的最优参数 λ_i 提供了一个合理的风险度量，我们通过求解以下优化问题来获得它
λ_i^* = arg min_{λ_i > 0} SURE_i(λ_i)。(6)
我们进一步通过特征分解 G_i = U_i Γ_i U_i^T 对目标函数进行处理，其中 U_i 是 G_i 的特征向量正交矩阵，Γ_i = diag(γ_{i,1}, ..., γ_{i,d-1}) 收集特征值。定义投影 z_i = U_i^T X_{-i,:} x_i^T，将 SURE 目标转化为 λ_i 的标量函数之和，每个网格点的求值成本为 O(d-1)[7 (https://arxiv.org/html/2608.17132#bib.bib19)]。(5) 中的 SURE 目标具有以下闭式表达式
SURE_i(λ_i)(7) = -nσ² + ||x_i||_2^2 - 2 ∑_{k=1}^{d-1} \frac{z_{i,k}^2}{γ_{i,k} + λ_i}
+ ∑_{k=1}^{d-1} \frac{γ_{i,k} z_{i,k}^2}{(γ_{i,k} + λ_i)^2} + 2σ² ∑_{k=1}^{d-1} \frac{γ_{i,k}}{γ_{i,k} + λ_i}。
给定最优参数 λ_i^*，相应的闭式岭估计为
\hat{w}_i^* = z_i^T diag(\frac{1}{γ_{i,k} + λ_i^*})_{k=1}^{d-1} U_i^T。(8)
最后，将所有的 \hat{w}_i^* 按行连接，并设置额外的对角元素为 0，便得到了邻接矩阵的最终软估计 \tilde{W}_{est}。

SURE 调谐的岭估计器捕获了总体层面的回归结构，确保了在噪声等方差假设下的可识别性[11 (https://arxiv.org/html/2608.17132#bib.bib1)]，而第 IV 节 (https://arxiv.org/html/2608.17132#S4) 的自适应阈值阶段则从所得的软邻接矩阵中提取 DAG 估计。

## IV自适应阈值

软邻接矩阵 \tilde{W}_{est} 通常不是一个 DAG。为了提取有效的 DAG 结构，我们提出一个三阶段过程。我们从一个与信号幅度成比例的绝对下限阈值开始，然后对高于此下限的幅度执行离散二分搜索，最后在二分搜索变得过于激进时，考虑一个保留信息性边的备用方案。

### IV-A无环性测试

我们使用[18 (https://arxiv.org/html/2608.17132#bib.bib5)]中的多项式无环性泛函，
h(W) = tr[(I_d + \frac{W ⊙ W}{d})^d] - d,(9)
它是最初在[20 (https://arxiv.org/html/2608.17132#bib.bib6)]中提出的矩阵指数形式的数值稳定变体。可以验证，h(W)=0 当且仅当 W 表示一个 DAG。在候选阈值图上计算 h 的计算复杂度为 O(d^3)，并且可以轻松地跨一批候选矩阵进行向量化计算。

### IV-B信号缩放下限滤波器

简单的零阈值会保留所有岭收缩噪声，而从零开始的二分搜索在密集场景下可能会将阈值推高到真实边幅度之上，从而返回空图。因此，我们引入一个尺度自适应下限
η_{min} = max(η_0, β max_{i≠j} ||(\tilde{W}_{est})_{i,j}||)，(10)
其中 η_0 和 β 是可调超参数。