深度因果发现：GraN-DAG

可视化演示

展示GraN-DAG将神经网络作为非线性结构方程模型的架构、连续化无环约束原理、以及非线性vs线性因果发现性能对比，交互切换三种模式探索深度因果发现。

因果发现困境

城市规划师拿到的那张“数据快照”，本质上就是一堆观测值——没有干预，没有实验，只有同时记录下来的数字。因果发现要做的，就是从这些数字里反推出“谁是谁的因”。

传统方法（如 PC 算法、FCI 算法）的典型思路，是先做条件独立性检验，再用一组规则去“拼接”出因果关系图。这套思路的问题在于两个“硬伤”：第一，它通常假设变量之间的关系是线性的——如果一个变量的平方、对数或指数会影响另一个变量，检验的效力就会大打折扣；第二，随着变量数量增加，可能的图结构数量呈现超指数增长（3 个变量就有 25 种 DAG，10 个变量则超过 \(4 \times 10^{18}\) 种），穷举式的搜索在 20 个变量的场景下就彻底失效了。

这个组合爆炸，实际上对城市规划师意味着：他根本不可能一条条尝试所有可能的道路网络方案，必须找到一种“自动收敛”的绘制方式——这正是可微因果发现登场的契机。

可微范式革命

如果把“画地图”这件事从“手工拼路”换成“机器自动调优”，关键突破在于：让图结构本身变成连续可微的参数。

2018 年，Zheng 等人提出 NOTEARS 方法，首次将 DAG 学习转化为连续优化问题。这里的核心魔法是：用一个加权邻接矩阵 \(W\) 来表示图结构，\(W_{ji}\) 表示从 \(X_i\) 到 \(X_j\) 的边权重。只要能够设计一个光滑函数 \(h(W)\)，使得 \(h(W)=0\) 当且仅当 \(W\) 对应一个有向无环图，那么原本“离散的搜索”就变成了“连续的约束优化”：

其中 \(\ell\) 是拟合损失（如最小二乘或负对数似然）。

GraN-DAG 正是沿着这条“可微革命”的道路前行，但把“线性结构方程”升级为“非线性神经网络建模”，让自动画笔不仅能画直线，还能画弯曲的河流与盘山公路。

有了这套连续化的“坐标系”，神经网络就可以无缝嵌入因果发现流程——这正是 GraN-DAG 的立足之本。

连续无环约束

城市规划师最担心的事，是画着画着，道路形成了一个“圆形回路”——从 A 到 B，从 B 到 C，再从 C 回到 A，这在因果图中是荒谬的（一个变量不可能直接或间接成为自己的因）。所以“禁止环路”是硬性规定。

问题在于：“有环 / 无环”本身是一个离散的组合概念，怎么把它写成能求导的公式？

答案在于邻接矩阵的迹函数。对于加权邻接矩阵 \(W\)（元素非负），\(W^k_{ii}\) 的值等于从节点 \(i\) 出发经过 \(k\) 条边回到节点 \(i\) 的“路径强度”之和。因此，矩阵的迹：

其中 \(d\) 是节点数。当且仅当 \(W\) 对应 DAG 时，\(h(W)=0\)。这是 NOTEARS 的经典结论，GraN-DAG 直接沿用了这个光滑约束。

理解了“连续无环约束”这枚指南针，我们再看 GraN-DAG 如何把神经网络装进这艘船。

Neq优化原理

GraN-DAG 的全称是 Gradient-based Neighborhood Optimization of DAGs，它的核心设计可以概括为一句话：用一组相互独立的神经网络子模块 + 一个全局邻接矩阵掩码，构建一个可微的非线性因果模型。

具体而言，对于 \(d\) 个变量 \(X_1, \dots, X_d\)，GraN-DAG 为每个变量 \(X_j\) 配备一个神经网络子模块 \(f_j\)，这个子模块的输入理论上是“所有其他变量的取值”，但通过一个逐元素的掩码（mask）操作来控制哪些输入真正影响输出：

其中 \(W_{\cdot j}\) 是邻接矩阵的第 \(j\) 列，\(\odot\) 是逐元素乘法。当 \(W_{ij}=0\) 时，变量 \(X_i\) 对 \(X_j\) 的预测没有任何贡献，等价于“没有这条因果边”。

这套设计相当于城市规划师给每个街区配了一个“独立规划小组”，每个小组负责设计自己街区的道路连接，但所有小组共用一张全局蓝图（邻接矩阵），并通过整体优化不断修订蓝图。

算法流程拆解

现在，我们把 GraN-DAG 的完整算法流程像“地图绘制工序”一样拆开来看。

第一步：初始化。 邻接矩阵 \(W\) 初始化为零或微小随机值，所有子模块的神经网络参数随机初始化。

第二步：联合优化。 GraN-DAG 使用增广拉格朗日方法（Augmented Lagrangian）来求解带约束的优化问题。每一步迭代，先固定无环惩罚系数，用 Adam 优化器更新网络参数和邻接矩阵，使预测损失下降；然后逐步增大惩罚系数，迫使 \(h(W)\) 趋向于 0：

第三步：权重裁剪与投影。 由于邻接矩阵在更新过程中可能出现负值，GraN-DAG 会对 \(W\) 施加非负约束（投影到非负象限），并对神经网络权重的范数进行裁剪，保证模型的表达能力受限但稳定，避免过拟合。

第四步：阈值化输出。 优化收敛后，将邻接矩阵中绝对值小于阈值的边视为不存在，得到最终的稀疏 DAG。

整个过程像一位耐心的绘图师：先粗线条勾勒全局，再反复修正曲线走向，最后擦除多余的辅助线，留下一张干净清晰的因果地图。

实践应用与局限

GraN-DAG 这样的“自动制图师”在实践中大显身手的典型场景包括：

基因调控网络推断。基因表达数据通常有数千个基因维度，且基因间的调控关系高度非线性。GraN-DAG 能够从基因表达谱中推断出调控网络，其非线性建模能力优于传统的线性假设方法。

金融风险传导分析。在金融网络中，各机构之间的风险传导关系错综复杂，GraN-DAG 可以通过观察市场数据推断出“谁的风险传染给谁”。

但自动绘图师也有自己的软肋：

• 计算开销大：每个变量一个神经网络，参数量随变量数线性增长，训练时间较长。

• 局部最优风险：增广拉格朗日是非凸优化，对初始值敏感，可能收敛到局部最优解。

• 因果充分性假设：GraN-DAG 假设不存在未观测的混杂变量，这在真实场景中往往难以满足。

城市规划师不会因为工具不够完美就放弃规划——他会理解工具的边界，并选择性地使用它。

---

AI 应用场景

Python 演示

import numpy as npimport torchimport torch.nn as nnimport torch.optim as optim torch.manual_seed(42)np.random.seed(42)  def simulate_linear_sem(n, d, edge_prob=0.3, noise_std=0.1):    """生成一个线性结构方程模型（SEM）的观测数据"""    W_true = np.zeros((d, d))    for i in range(d):        for j in range(i + 1, d):            if np.random.rand() < edge_prob:                W_true[j, i] = np.random.uniform(0.5, 1.5)    X = np.random.randn(n, d)    for j in range(1, d):        X[:, j] = X @ W_true[j, :] + np.random.randn(n) * noise_std    return torch.tensor(X, dtype=torch.float32), W_true  class GraNDAGSubmodel(nn.Module):    """每个变量对应一个神经网络子模块"""    def __init__(self, d, hidden=64):        super().__init__()        self.net = nn.Sequential(            nn.Linear(d, hidden), nn.ReLU(),            nn.Linear(hidden, hidden), nn.ReLU(),            nn.Linear(hidden, 1),        )        self.mask = nn.Parameter(torch.zeros(d), requires_grad=False)     def forward(self, x):        return self.net(x * self.mask).squeeze(-1)  def h_func(W):    """连续无环约束：h(W) = trace(exp(W)) - d"""    d = W.shape[0]    return torch.trace(torch.matrix_exp(W)) - d  def train_gra_ndag(X, n_iters=500, lr=1e-2, lambda_init=1.0, rho=1.0):    n, d = X.shape    models = [GraNDAGSubmodel(d) for _ in range(d)]    params = [p for m in models for p in m.parameters()]    mask_params = [m.mask for m in models]    for m, i in zip(models, range(d)):        m.mask.data = torch.rand(d) * 0.1    optimizer = optim.Adam(params + mask_params, lr=lr)    lam = lambda_init     for it in range(n_iters):        optimizer.zero_grad()        W = torch.stack([m.mask for m in models], dim=1)  # d x d        W = torch.relu(W)  # 非负投影        loss = 0.0        for j, m in enumerate(models):            pred = m(X)            loss += torch.mean((X[:, j] - pred) ** 2)        h_val = h_func(W)        total = loss + lam * h_val + 0.5 * rho * h_val ** 2        total.backward()        optimizer.step()        lam += rho * h_val.item()        if it % 100 == 0:            print(f"iter {it}: loss={loss.item():.4f}, h={h_val.item():.4f}")     W_final = torch.stack([m.mask for m in models], dim=1)    return torch.relu(W_final).detach().numpy()  if __name__ == "__main__":    X, W_true = simulate_linear_sem(n=500, d=6)    print("真实邻接矩阵:\n", W_true)    W_est = train_gra_ndag(X, n_iters=300)    print("学习到的邻接矩阵:\n", np.round(W_est, 2))

1. mask = nn.Parameter(torch.zeros(d), requires_grad=False)每个子模块持有一个掩码向量，该向量的值就是邻接矩阵的一列，通过梯度更新逐步逼近真实因果边。
2. self.net(x * self.mask)掩码与输入逐元素相乘，实现“软剪枝”：掩码接近 0 的输入被屏蔽，不参与预测。
3. torch.trace(torch.matrix_exp(W)) - d连续无环约束的计算，基于矩阵指数求迹。
4. W = torch.relu(W)非负投影，保证邻接矩阵权重非负（因果边没有负方向）。
5. lam += rho * h_val.item()增广拉格朗日更新：每轮迭代增大惩罚系数，迫使 h(W) 收敛到 0。

常见误区

面试题