企业官网建设流程全解析

先从一句实话说起：我第一次看到“DAGs with NO TEARS”这个标题时，第一反应是“哪篇论文起名字还带谐音梗”。DAG是有向无环图（Directed Acyclic Graph）的缩写，因果推断领域的老熟人；NO TEARS字面是“没有眼泪”，细看又能拆成“NOT EARS”——一种“不需要靠暴力搜索、不用把图反复撕裂”的优化思路。等我把这篇论文的内容完整过了一遍，才意识到这个双关起得相当准确：它把“学习一个有向无环图”这个传统上被当成离散组合搜索的问题，改写成了一个光滑的连续优化问题，并给出了能稳定收敛的增广拉格朗日求解框架。

如果你需要在观测数据上还原变量之间的依赖关系和因果方向，或者在做贝叶斯网络、结构方程建模、因果发现相关的工作，这篇论文的思路绝对值得吃透。这篇文章我不打算复述原论文的公式推导，而是从“为什么以前的方案让人掉眼泪”“这个算法的关键等价关系是怎么想出来的”“我自己复现时踩过的坑”这几个角度展开，最后给出一份可以直接跑的简化实现。

1. 结构学习为什么难：组合爆炸与“环”的麻烦

1.1 结构学习到底在解什么题

结构学习（structure learning）的任务听起来很直白：给定一个观测矩阵 X，它有 n 行样本、d 列变量，我们希望学出一个 d×d 的加权邻接矩阵 W，其中 W_ij 非零表示变量 i 对变量 j 有直接影响。如果从因果角度解释，这个非零权重就是“i 是 j 的因”的数值化表达。

绝大多数实际场景里，我们拿到手的只是系统最终呈现出来的观测数据。比如一堆临床指标，我们想知道哪些指标是上游病因、哪些是下游症状；一座工厂里几十个传感器读数，我们想知道哪些参数驱动了哪些状态。只是拟合一个预测模型是不够的，因为预测模型可以任意地把相关的变量互相连起来，但因果结构要求的是一个“有向无环图”，也就是 DAG。

为什么要无环？因为真实的因果不能自相矛盾。如果 A 是 B 的原因，B 又是 A 的原因，那这个系统连数学定义都站不稳：到底谁先存在？谁决定谁？如果允许环存在，线性结构方程 X = XW + Z 中，X 的解可能不存在或者不唯一，整个推断就崩了。所以 DAG 不是人为偏好，而是因果模型本身的内在要求。

1.2 “无环”是所有麻烦的源头

如果抛开无环约束，学习 W 就是一个普通的最小二乘问题：X ≈ XW，用回归就能解，闭式解、梯度下降都行。麻烦全在“无环”这两个字上。

最朴素的结构学习器是这样工作的：定义一个打分函数，比如 BIC、BDeu 或者贝叶斯边际似然，然后遍历所有可能的 DAG，选出打分最高的那个。问题是 DAG 的数量增长率相当可怕。d 个节点的 DAG 数量满足一个经典的递归式：

a(n) = Σ_{k=1}^{n} (-1)^{k+1} C(n, k) · 2^{k(n-k)} · a(n-k)，其中 a(0)=1。

算一下就知道：d=10 的时候，DAG 数量已经达到千万量级；d=20 的时候，数量级在 10^20 以上。这种规模别说枚举，连设计一个聪明的剪枝策略都很吃力。所以早期方法要么用爬山法在相邻 DAG 之间跳来跳去，要么用条件独立性检验逐步判定边的有无和方向。这些方法不是不能用，而是稳定性差、结果对初始点敏感，而且在高维场景下很容易被噪声带跑。

1.3 为什么连续优化一直做不了

直觉上，把 W 当成实数矩阵，用梯度下降不是更自然吗？问题在于“无环”这个约束不是一个光滑条件，甚至不是一个普通的非凸约束。

你可以试着把“有环”当作一个惩罚项加到损失函数里，比如检测出环就加大惩罚。但环检测是一个全局性质：某条边的权重只需要改变一点点，就可能让整个图从无环变成有环，或者在某个局部修改中把环保留下来。梯度信息在这种“阶跃式”变化面前基本失效。更麻烦的是，所有有环矩阵构成的集合在矩阵空间里不是闭的；你很难定义一个稳定的投影，把任意矩阵投影到“最近的无环矩阵”上，因为最近的那个矩阵可能根本不存在。

所以这个问题的核心障碍不是损失函数，而是约束本身。要真正用连续优化去解它，就必须找到一个性质足够好的“无环指示器”：它的值等于 0 当且仅当图无环，而且它必须是 W 的光滑函数。这在当时看起来像一个很难满足的要求。

2. 核心技巧：把无环约束改写成光滑等式

2.1 一条边如何变成矩阵幂的计数

论文里最关键的一步，是找到一个令人拍案的矩阵等价关系。它利用了路径计数的基本事实：对于一个元素非负的方阵 M，M^k 的第 i 行第 j 列元素，正好等于从节点 i 出发、沿着有向路径走恰好 k 步到达节点 j 的路径权重乘积之和。

这里有个小细节：邻接矩阵 W 的权重可能有正有负，而我们要计数的是路径的“存在性”，不能被正负符号抵消。所以论文把 W 逐元素平方，得到 M = W ∘ W（Hadamard 乘积，即对应位置相乘）。平方之后所有元素非负，路径权重乘积也非负，就不会出现“一条正路径和一条负路径互相抵消”的假象。

然后考虑矩阵指数：

exp(M) = I + M + M²/2! + M³/3! + ...

它的对角元 exp(M)_ii 等于从 i 出发、走任意长度之后回到 i 的所有路径权重之和（按长度加权）。对于一个无环图，从任意节点出发，只要走至少一步，就不可能回到自己，所以所有 k ≥ 1 的项对角元都是 0。此时 exp(M) 的对角元只有 I 贡献的 1，整个矩阵的迹恰好等于节点数 d。反过来，如果图里有环，那么环上的节点在某个幂次之后必然出现非零对角元，矩阵的迹一定严格大于 d。

这就是那个著名的等价关系：

h(W) = tr(exp(W ∘ W)) - d，满足 h(W) = 0 当且仅当 W 对应的图是有向无环图。

用生活化的类比说：把每条有向边想象成一条单向通道，W∘W 的幂次就是在数“从某个路口出发，逛若干条路之后还能绕回这个路口”的带权路径条数。没有环的路网里，这种“回到原点”的路线根本不存在；有环的路网里，只要你愿意绕，回到原点的路线一定存在，而且权重为正。

2.2 迹等于节点数：一个干净的充要条件

这个等价关系为什么能成立，值得再仔细展开一下。它其实依赖两个方向都很干净。

一方面，如果 W 对应的图存在某个环，比如 v1→v2→…→vk→v1，那么这些边权重平方的乘积是一个正数。在 M^k 中，这个环会给对角元 (M^k)_v1v1 贡献一个正项。由于 k ≥ 环长，之后的所有更高次幂也都保留着这条回路的贡献，exp(M) 的对角元严格大于 1，迹严格大于 d。

另一方面，如果图是无环的，那么任意 k ≥ 1 都不存在长度为 k 的回路，所以 M^k 的所有对角元都是 0。矩阵指数展开后，对角元只来自 k=0 的单位阵，迹精确等于 d。

这里有一个特别好的性质：h(W) 是 W 的光滑函数。矩阵指数本身是解析函数，逐元素平方也是可微运算，所以 h(W) 对 W 的任意元素都有连续的导数。这意味着你可以正大光明地使用梯度信息，让优化器自己去找“无环区域”。

还有一个微妙但值得注意的点：h(W) 永远大于等于 0。因为 M 是非负矩阵，exp(M) 的对角元至少是 1，加上更高幂次的非负贡献，所以 h(W) ≥ 0 恒成立。这让优化过程的收敛判据变得很干净——你只需要看 h 是否下降到接近 0，不需要担心它“冲过头”变成负数再折返。

2.3 增广拉格朗日：比罚函数更聪明的做法

有了光滑的约束 h(W)=0，下一步就是怎么把约束融进优化目标。如果只是简单地在损失函数后面加一个 ρ·h(W)²，那就是传统罚函数法，问题在于 ρ 必须趋于无穷大才能让 h 精确为 0，而太大的 ρ 会让优化问题变得病态。

论文选择的是增广拉格朗日方法（ALM）。思路是这样：把目标函数写成

F(W) = L(W) + λ||W||₁ + α·h(W) + (ρ/2)·h(W)²

其中 L(W) = (1/(2n))||X - XW||²_F 是线性结构方程下的负对数似然，也就是最小二乘损失。优化时固定 α 和 ρ，先用梯度方法求解 W；然后更新乘子 α ← α + ρ·h(W)；如果 h 下降得不够快，就把 ρ 调大一些。如此重复，直到 h 收敛到 0。

为什么用 ALM 而不是普通拉格朗日？普通拉格朗日只对凸问题有比较好的收敛保证，而这里的目标函数非凸；罚函数又需要 ρ 无限增大。ALM 的巧妙之处在于，乘子 α 会自动累积“约束误差的信息”，即使 ρ 保持在一个有限值，也能在迭代过程中逐步逼出满足约束的解。这个套路在带约束的非凸优化里是相当实用的选择。

也顺带提一句，论文中讨论过其他等价形式，比如用 -log det(I - W∘W) 来刻画无环条件。它在谱半径小于 1 的区域与 exp 版本等价，数值上对某些病态场景更稳定；但 exp 版本的实现更直观，也是社区复现时用得最广的。我下文给出的代码也基于 exp 版本。

3. 自己上手复现：数据、代码与调参记录

3.1 先造一份带真实答案的数据

验证算法最好的方式，是先构造一个已知的 DAG，生成满足线性结构方程的数据，再让算法从纯观测数据里还原它。

我用这样的方式生成合成数据：

1. 固定节点数 d=20，随机生成一个拓扑顺序，只允许从序号小的节点指向序号大的节点，这样能天然保证无环。
2. 每条边以一定概率出现，比如 0.3，出现时权重在 [-2, -0.5] ∪ [0.5, 2] 里随机取值。
3. 按线性结构方程 X = XW_true + Z 生成 n=1000 个样本，其中 Z 是独立高斯噪声。

生成得到 X 之后，把每一列做标准化处理（减均值除以标准差）。这一步很多人会忽略，但它对后续的梯度优化影响很大：如果某个变量量纲特别大，平方损失会主导整个目标，导致算法只去拟合那个变量而忽略其他结构。

3.2 核心代码：几十行跑通 NOTEARS

下面这个版本是我为了教学和数据实验写的简化实现，使用 PyTorch 的自动微分来算梯度，逻辑上对应增广拉格朗日迭代。它不是一个高度工程化的库，但足够让你在 d ≤ 30 的小规模数据上跑通整个流程。

import numpy as np import torch def h_func(W): """无环约束：h(W)=0 当且仅当 W 对应的图是有向无环图""" M = W * W return torch.trace(torch.matrix_exp(M)) - W.shape[0] def notears_alm(X, lambda1=0.02, lr=0.01, rho=1.0, max_outer=50, max_inner=200, tol=1e-6): X = torch.as_tensor(X, dtype=torch.float64) n, d = X.shape W = torch.zeros(d, d, dtype=X.dtype, requires_grad=True) # 破对称：全零初始会让梯度为零，学不出任何边 with torch.no_grad(): W += torch.rand(d, d, dtype=X.dtype) * 1e-4 alpha = 0.0 h_prev = float("inf") for t in range(max_outer): # 内层：固定 alpha 和 rho，梯度下降优化 W for _ in range(max_inner): residual = X - X @ W loss = 0.5 / n * (residual ** 2).sum() loss = loss + lambda1 * W.abs().sum() h = h_func(W) loss = loss + alpha * h + 0.5 * rho * h * h loss.backward() with torch.no_grad(): W -= lr * W.grad W.grad.zero_() # 外层：更新 ALM 乘子与惩罚系数 with torch.no_grad(): h = h_func(W).item() alpha += rho * h if h > 0.1 * h_prev: rho = min(rho * 2, 1e6) h_prev = h print(f"outer {t}: h={h:.3e}, alpha={alpha:.1e}, rho={rho:.1e}") if abs(h) < tol: break return W.detach().numpy()

说明几点。内层我用的是纯梯度下降，每一步都重新计算一次 h 和它对应的梯度；这样写是为了逻辑清楚。实际跑的时候你会发现，内层迭代多了会有点慢，因为矩阵指数 expm 的计算复杂度是 O(d³)，d 小的时候无所谓，d 上百之后每一轮都肉眼可见地卡一下。工程上可以把内层换成 L-BFGS 或者 Adam，收敛会快不少。

还有一个细节：W 初始化为零矩阵会导致梯度为零，因为最小二乘损失和 L1 惩罚在零点的梯度都是零，h 的梯度也是零。所以我加了 1e-4 量级的随机扰动来破对称。这个坑很容易踩，很多人复现失败都是因为在全零初始处“原地不动”。

3.3 后处理与阈值化

算法输出的是一个稠密矩阵 W，绝大多数元素的绝对值都很小但不严格为 0。这是因为 L1 惩罚只会让权重压缩到接近 0，不会直接置零。所以要得到稀疏的图结构，还需要一步后处理：阈值化。

最简单的做法是设置一个阈值，比如 0.05，把绝对值小于阈值的元素直接置 0。更讲究的做法分两步：

1. 先用阈值或稀疏规则确定边的存在性，得到一个无环骨架。
2. 再对每个变量的“父节点集合”做一个不带无环约束的稀疏回归，重新精修权重。

第二步可以显著降低误检率，因为无环约束是一种全局软约束，它对权重的估计是“够用就好”，而不是为预测精度服务的。精修回归可以让保留下来的边权重更接近真实值。

文件里还应该有一个可复现的评估函数。如果你知道真实图 W_true，可以计算 SHD（Structural Hamming Distance），也就是“预测边集合与真实边集合的对称差数量”。这个指标在合成数据上很直观：SHD 越小，恢复得越好。

3.4 参数怎么定：λ、rho、lr

这三个参数是复现成败的关键，我之前的实测感受如下。

• λ 控制稀疏强度。太大会让所有权重都被压成 0，输出一个全零矩阵；太小则假边泛滥。建议在 [0.001, 0.01, 0.05, 0.1] 这样一个序列里扫描，观察 h 的收敛曲线和 SHD。数据标准化之后，0.01~0.05 通常是比较合理的起点。
• ρ 决定无环约束被满足的速度。初始值 1 是个稳妥起点。如果 h 下降很慢，可以在外层循环里按 2~10 倍速率放大。ρ 太小会导致需要很多轮外循环才能把 h 压下去；ρ 太大又会让优化器过于激进地满足无环约束，忽略拟合损失。
• lr 的取值依赖是否做了数据标准化。标准化之后，0.01 的梯度下降步长在 d=20 的合成数据上是可以工作的；如果数据没标准化，步长可能需要大幅缩小，否则损失直接发散。

这里还要强调一个容易出问题的地方：h(W) 是一个全局量，它的梯度会同时影响所有边。有些边可能只是为了压制某个环而存在，去掉它们会让 h 重新大于 0。所以阈值化之后不要再盯着原始的 h 值不放，应该重新计算一遍阈值化后的 h_func 来确认得到的图确实无环。

4. 实测中的坑与排查速查表

4.1 h 值死活降不下去怎么办

最常遇到的现象是：外循环跑了几十轮，h 始终停留在 1e-2 甚至更高的数量级，不再继续下降。这种情况我总结出三个排查方向。

第一，内层迭代次数不够。W 在固定的 α 和 ρ 下还没被充分优化，外层就急着更新乘子。可以加大 max_inner，或者换成 L-BFGS 这种二阶方法。第二，ρ 的更新太慢。如果 h 下降速度长期低于预期，应该更激进地放大 ρ，比如从 2 倍改成 10 倍。第三，学习率设置不当。梯度下降收敛慢的表现是损失曲线“磨蹭”，迟迟压不到最优。试着把 lr 放大到 0.05~0.1，或者换优化器。

还有一个容易被忽略的问题：数据噪声过大时，h 在收敛后期会有小幅振荡，绝对数值下不到 1e-8。这通常是正常的，不用死磕到机器精度。模型选择上，只要 h 降到比如 1e-4 及以下，图的拓扑结构已经基本稳定了。

4.2 输出全是零矩阵或者假边一堆

全零矩阵基本可以锁定在 λ 太大，或者 W 初始化为全零后梯度被卡住。前者调小 λ，后者加随机扰动。

假边太多则相反，说明 λ 太小，或者阈值取得太低。这种情况下，损失函数里拟合部分的权重太高，算法倾向于把能改善预测的边全部留下。一个实操技巧：先固定一个较大的阈值，比如 0.1，看保留下来的边集合里有多少是真实边；如果召回率不够，再逐步降低阈值。你会发现 SHD 往往在某个阈值附近出现一个明显的低值，这就说明数据质量没问题，只是你没找对后处理尺度。

4.3 阈值化之后反而出现环

这个事情我第一次遇到的时候很困惑：明明算法训练完 h≈0，阈值化之后却出现了一个小环。原因是无环约束是基于加权环的权重和来判定的，去掉一条弱边可能破坏原有环路上的“权重平衡”，让另一条反向路径成为新的环。

解决办法是阈值化之后重新计算 h_func，不要默认后处理会保留无环性。如果发现环，我的建议是：优先剔除环路上权重最小的那条边，而不是简单地把所有小权重全部清零。你可以写一个迭代剪枝函数，每次找一个环，删环上最弱边，直到图无环为止。这在稀疏图上的效果比单次阈值化稳定得多。

4.4 方向识别不了不等于代码错了

合成数据上还有一个特别容易打击信心的情况：SHD 不为 0，而且报告的边方向有几条标反了。这不是代码 bug。

因为 DAG 存在马尔可夫等价类：有些结构，比如 v1 → v2 → v3 和 v1 ← v2 ← v3，它们蕴含的条件独立关系完全相同，数据在统计上根本区分不了这两种方向。算法学出的可能是等价类中的任何一个方向，这属于理论不可识别问题，不是算法缺陷。评估的时候，建议同时报告骨架 SHD（只看边存在与否）和有向 SHD（还要看方向），两者对比能更准确地判断算法到底学到了什么程度。

4.5 数据不是线性 SEM 时怎么办

NOTEARS 的打分函数建立在线性结构方程假设之上。如果真实系统是非线性机制，直接套用线性损失会学到很多伪边。我在一个非线性机制的数据上试过，误差集中在噪声方差大的变量上，因为线性模型强行用线性项去逼近非线性效应，结果就是多加边、加大权重。

对付这种情况有两种思路：一是先检查残差，如果 X - XW 在某个变量上呈现明显的非线性模式，就说明线性假设不成立；二是转向非线性版本的结构学习方法，比如用神经网络建模变量之间的条件期望再做无环连续优化。论文组后续有一系列工作就是沿着这条路走的，用神经网络或加性模型替换线性得分部分。核心的无环约束 h(W)=0 并没有变，变的只是损失函数的表达方式。

结尾：一点复现体会

我自己复现这套算法时，最大的感受是“名字起得好不如数学选得巧”。NO TEARS 听起来像一句口号，但它真正厉害的地方在于把“环”这种离散的、全局的、不光滑的性质，翻译成了一个由矩阵指数定义的、光滑的、可以用梯度信号持续优化的等式约束。这个翻译一旦完成，后面的优化框架反而变得顺理成章。

如果只记住一个实操技巧，我会建议你一定要在动手写代码之前先把 h_func 验证明白：随便构造一个有环矩阵和一个无环矩阵，打印一下 h 值，你会立刻理解这个约束的几何意义。整个算法最核心的代码其实只占几十行，剩下的全是数据处理和参数调优。先在小规模 d=10~20 的合成数据上跑通，再上真实数据，你会少踩很多坑。