企业官网建设流程全解析

1. 从"组合爆炸"谈起：为什么结构学习曾经这么痛苦

做因果推断或概率图模型的朋友，对结构学习（Structure Learning）应该都不陌生。我们手里有一堆观测数据，希望推断出变量之间的依赖关系，甚至是有向的因果关系——也就是一张有向无环图（DAG）。这个问题的应用范围非常广，从基因调控网络重构、金融风险传导分析，再到推荐系统里的因果纠偏，本质上都绕不开"先学一张图"这一步。

但传统DAG学习最大的麻烦在于：搜索空间是离散的、组合爆炸式的。一个包含d个节点的DAG数量近似为：

[ f(d)=\sum_{k=1}^{d}(-1)^{k+1}\frac{d!}{(d-k)!}2^{d(d-k)} ]

d=10时，这个数字大约是 (4.8 \times 10^{18})，到了d=20就已经是天文数字。如果你用经典的贪心搜索或动态规划去枚举图结构，只能应对个位数的节点；用PC算法这类基于条件独立性检验的方法，在高维、非线性、样本量小的情况下又极不稳定。所以相当长一段时间里，大家觉得结构学习是个"可以定义清楚、但很难规模化"的问题。

"DAGs with NO TEARS"这篇工作，正是从根上改写了这个问题的"打开方式"。它的全称我习惯叫NOTEARS（Non-combinatorial Optimization via Trace Exponential and Augmented lagRangian for Structure learning），核心思路用一句话概括就是：把"图是否是无环的"这个离散组合条件，等价转化成一条光滑、可微的连续等式约束，然后把整个结构学习变成标准的连续约束优化问题。这样一来，梯度下降、增广拉格朗日这类经典连续优化工具箱都可以直接搬过来用，不再需要跟离散DAG空间硬碰硬。

我第一次读到这个思路时确实挺惊讶的——这个"将离散约束连续化"的切入点看着很简单，但之前就是没人想到。这篇博文想把它的核心机制、完整复现步骤以及我实际跑实验时踩过的坑都拆开讲一讲，适合正在做因果发现、图结构学习，或者想给自己的模型引入无环约束的朋友参考。

2. 核心思路拆解：NO TEARS 究竟是怎么把"环路"写进数学公式的

2.1 三个关键成分的化学反应

要理解NOTEARS，得先明确它针对的标准设定：假设变量 (x = (x_1, ..., x_d)) 之间满足线性结构方程模型（Linear SEM）：

[ x_i = w_i^T x + b_i + z_i, \quad i=1,...,d ]

其中 (z_i) 是噪声项，(w_i) 是第i个变量对其他变量的回归系数向量。把所有系数拼成一个加权邻接矩阵 (W \in \mathbb{R}^{d\times d})，那么结构学习就变成：从数据X中学习 (W)，同时也要求 (W) 对应的是一个DAG。

于是问题可以写成：

[ \min_{W} ; \ell(W; X) \quad \text{s.t.} \quad G(W) \in \text{DAGs} ]

这里的 (\ell(W;X)) 是损失函数，比如最小二乘损失 (\frac{1}{2n}||X - XW||^2_F)。难点全在约束条件上："G(W)是DAG"不是一个可以用梯度直接处理的约束。

NO TEARS做了一件非常聪明的事——构造了一个标量函数来刻画"无环程度"：

[ h(W) = \mathrm{tr}(e^{W \circ W}) - d ]

看到这个公式先别慌，我们拆开看。(W \circ W) 是逐元素平方，得到一个非负矩阵 (A = W \circ W)；(e^A) 是矩阵指数（注意：是矩阵指数，不是逐元素指数，这是整个构造的核心）；(\mathrm{tr}(e^A)) 是矩阵指数对角线元素之和；d 是节点数。

这个函数的奇妙性质在于：** (h(W)=0) 当且仅当 (W) 对应的图是无环的。** 换句话说，DAG约束被等价地转换成了一个光滑、可微的等式约束，可以直接塞进连续优化框架里。我刚开始也没想通它为什么能捕捉环路，后来用一个简单例子给自己验证了一遍。

2.2 用两节点三节点小例子验证：环路是怎么被"抓"出来的

考虑只有两个节点的情况。如果 (W = \begin{bmatrix} 0 & a \ 0 & 0 \end{bmatrix})，这是一个从1指向2的DAG，矩阵指数的性质此时很好算：

[ A = W \circ W = \begin{bmatrix} 0 & a^2 \ 0 & 0 \end{bmatrix} ]

这是个严格上三角矩阵，它的矩阵指数等于 (e^A = I + A)（高次项为0），所以 (\mathrm{tr}(e^{W \circ W}) = 2)，正好等于d，于是 (h(W)=0)。

但如果出现一条回路：(W = \begin{bmatrix} 0 & a \ b & 0 \end{bmatrix})，那么：

[ A = \begin{bmatrix} 0 & a^2 \ b^2 & 0 \end{bmatrix} ]

计算矩阵指数时，注意到：

[ A^2 = \begin{bmatrix} a^2 b^2 & 0 \ 0 & a^2 b^2 \end{bmatrix} ]

所以 (\mathrm{tr}(e^A) = 2 + a^2 b^2 + \dots)，显然大于2。只要 (a\neq 0) 且 (b\neq 0)，这个环就一定会被捕捉到，(h(W)>0)。更一般的结论是：矩阵指数展开式中，对角线元素会累积"环的强度"，非零环路不可能同时让对角线上所有累计项为零。

2.3 梯度计算：从 h(W) 到优化器能用的 ∇h(W)

有了约束函数，还得有梯度，优化器才能跑起来。NO TEARS给出的梯度公式非常优雅：

[ \nabla h(W) = (e^{W \circ W})^T \circ (2W) ]

这里的 (\circ) 是哈达玛积（逐元素乘积），(e^{W\circ W}) 依然是矩阵指数。这个梯度的直觉是：对于每一条边，它的梯度会结合"经过该边的环路总强度"和"边自身的权重"，如果某条边处于大量环路中，优化器会给它更大的压力来打破环路。实际工程实现时需要注意：矩阵指数这一步千万别手写泰勒展开，直接用NumPy的 scipy.linalg.expm 或 PyTorch的 torch.matrix_exp，数值稳定性远胜手工实现。

到这里，整个问题的形态已经完全变了：

[ \min_{W} ; \ell(W; X) \quad \text{s.t.} \quad h(W)=0 ]

一个标准的带等式约束的光滑优化问题。

3. 完整实操过程：从损失函数到增广拉格朗日的落地实现

如果你只是理解了上面的数学，距离真正跑通NOTEARS还有不少距离。我先说结论：整个可复现的流程包含三个模块：损失函数构造、增广拉格朗日外层循环、结果阈值化处理。下面逐个拆开讲。

3.1 损失函数怎么选：最小二乘是默认，但不是唯一

原文默认使用线性SEM下的最小二乘损失：

[ \ell(W; X) = \frac{1}{2n} |X - XW|_F^2 ]

这里 (X) 是 (n\times d) 的数据矩阵。最小二乘的好处是：在高斯加性噪声假设下，它对应最大似然估计，与线性SEM的假设完全一致。我自己的实验中也发现，最小二乘损失的收敛相对平滑，是新手最不容易翻车的选择。

但需要注意：如果数据明显是非高斯的，最小二乘依然能用，但会丢掉一些分布信息。后续有一些工作改成使用似然损失甚至秩损失，不过初次复现NOTEARS不需要在这上面过度纠结，先用最小二乘把框架跑通，后面再按需更换。

实际经验：当变量维度 d 超过50时，我建议把损失改成带L1正则的版本 (\ell + \lambda_1 |W|_1)，因为纯最小二乘很可能学到全连接稠密图，阈值化后虽然也能用，但SHD（结构汉明距离）指标通常会变差。

3.2 增广拉格朗日外层循环：这是整个算法的核心发动机

带等式约束的优化问题，最经典的处理方法是增广拉格朗日法（Augmented Lagrangian Method）。NOTEARS并没有上来就暴力地把 (h(W)=0) 罚进损失里，而是构造了一个每个迭代步都在变化的子问题：

[ W^{(k+1)} = \arg\min_W ; \ell(W) + \lambda_k h(W) + \frac{\rho_k}{2} h(W)^2 ]

其中 (\lambda) 是拉格朗日乘子，(\rho) 是二次惩罚系数。整个外层循环套路如下：

1. 初始化：(W=0)，设置初始 (\lambda=0)、(\rho=1)，以及一个增长因子比如 (c=10)，还有收敛阈值比如 (\epsilon=10^{-8})。
2. 重复：
   • 用梯度优化器（L-BFGS最稳，Adam也行）求解上面的无约束子问题，直到收敛；收敛判据是梯度范数小于某个阈值。
   • 更新对偶变量：(\lambda \leftarrow \lambda + \rho \cdot h(W))
   • 增大惩罚：(\rho \leftarrow \rho \cdot c)
   • 直到 (h(W)) 的绝对值小于设定的阈值（比如 (10^{-8})）。
3. 最后一轮迭代结束后，对 (W) 做阈值化：小于阈值的权重置为0，得到最终的稀疏DAG。

这套流程有个很关键的内在逻辑：一开始 (\rho) 很小，优化器主要在降低损失；随着 (\rho) 增大，约束逐渐变硬，最后强制把图推向无环。整个过程中 (\lambda) 也在不断累积，防止 (\rho) 单方面过大导致问题病态。我模拟过这个动态过程，前几个外层循环损失下降很快，后面几个循环的损失甚至略有回升，但只要h(W)不断逼近0，最后结果大概率是好的——这是正常的，不用慌张。

3.3 给一份可直接抄作业的伪代码

我用Python风格写一个核心骨架，自己实操时可以参考这个版本：

import numpy as np from scipy.linalg import expm from scipy.optimize import minimize def h_function(W): # W: d x d 加权邻接矩阵 d = W.shape[0] M = expm(W * W) return np.trace(M) - d def grad_h(W): M = expm(W * W) return 2.0 * W * M.T def loss(W, X): n, d = X.shape W = W.reshape(d, d) XW = X @ W diff = X - XW return 0.5 * np.sum(diff**2) / n def grad_loss(W, X): n, d = X.shape W = W.reshape(d, d) return (X.T @ (X @ W - X)) / n def obj_func(W_flat, X, lam, rho): W = W_flat.reshape(-1, d) l = loss(W, X) h = h_function(W) return l + lam * h + 0.5 * rho * h * h def grad_obj(W_flat, X, lam, rho): W = W_flat.reshape(-1, d) g = grad_loss(W, X).flatten() h = h_function(W) gh = grad_h(W).flatten() return g + lam * gh + rho * h * gh # 主流程 d = X.shape[1] W_flat = np.zeros(d * d) lam = 0.0 rho = 1.0 c = 10.0 h_thresh = 1e-8 while True: res = minimize(obj_func, W_flat, args=(X, lam, rho), method='L-BFGS-B', jac=grad_obj, options={'maxiter': 50, 'ftol': 1e-7}) W_flat = res.x h_val = h_function(W_flat.reshape(d, d)) lam = lam + rho * h_val rho = rho * c if abs(h_val) < h_thresh: break W_est = W_flat.reshape(d, d) W_est[np.abs(W_est) < 0.3] = 0.0 # 阈值化

关于这段代码，有几个亲自验证过的要点：

• 内层优化器选L-BFGS-B：二次惩罚项让目标函数局部曲率变化比较剧烈，L-BFGS这类拟牛顿法对这类问题的适应性远超普通SGD。如果你坚持用Adam，学习率要调到非常小（比如1e-4以下），否则外层迭代很难收敛到无环解。
• maxiter不要太长：子问题有时并不需要充分优化，充分优化的性价比不高。我把maxiter设为50左右，整体速度明显快于原文默认的100，且最终SHD基本无差别。
• 阈值0.3不是固定值：这个超参数需要根据数据噪声水平调整。噪声大时阈值要适当提高，否则一小堆非零噪声权重会被误判为边。

3.4 阈值化算不算作弊？谈谈为什么这一步是必要的

很多人会问：你不是把DAG约束直接写进优化里了吗，为什么最后还要手动阈值化？答案是：h(W)=0约束保证的是"图无环"，但并没有保证"图稀疏"。一个完全稠密的无环图（比如所有边按拓扑序从上往下排列）也是满足约束的，但它显然不是我们想要的因果结构。

阈值化处理正是为了在"无环"的基础上再施加"稀疏先验"。实际操作中，我见过有人跳过阈值化直接用绝对权重排序来选边，效果也差不多。只要不是把阈值设得太离谱，这一步对整个算法是锦上添花而非埋雷。

4. 实验对比与调参避坑：我在真实数据上遇过的那些事

理论说得头头是道，真上手跑的时候坑不少。我说几个自己实际踩过的坑，都是原文读不出来、必须跑了才会知道的细节。

4.1 关于收敛阈值：别把h(W)压得太狠

刚跑NOTEARS时，我心想既然约束可以趋近0，那就把收敛阈值设成1e-14好了，结果非但没让结果更好，反而训练慢了几倍。后来总结出经验：h(W)在1e-6到1e-8之间就足够了。因为h(W)的尺度跟矩阵指数有关，当它小于1e-8时，图中残存的"有效环路强度"在浮点精度下已经可以忽略不计。压到1e-14除了让优化器多做几十轮，不会带来任何精度收益。

4.2 噪声水平与阈值选择之间强耦合

我做过一组对比实验：数据由线性SEM生成，噪声标准差分别取0.1、1.0和5.0，固定不变地使用0.3阈值。结果在噪声标准差为5.0的数据上，阈值0.3会把许多真实系数（比如0.4的边）直接置零，导致召回率断崖下跌；而噪声0.1时，阈值0.1以下才能保留真实弱边。所以一个靠谱的做法是：先在无监督情况下确定噪声尺度，再设定阈值；或者直接对W估计值做非参数聚类，把"明显大于0"和"趋近0"的权重划开。

4.3 常见问题速查表

我把遇到的问题整理成了一张表，方便大家直接对照排查：

4.4 初始化方式对结果的影响

原文是使用全零矩阵初始化，实测下来全零初始化确实是最稳的——它的梯度无偏、不会预先偏向某个方向。但我试过用随机初始化，结果发现如果随机初始化矩阵的谱半径偏大，矩阵指数会迅速爆炸，h(W)初始值飙升到几千，之后优化过程要花很多轮次才能把约束拉回来。建议保持全零初始化，如果要加扰动，扰动幅度要控制在1e-3以内。

5. 评估指标怎么用：怎样科学地"证明"你学到了DAG

学习出来的W不会自己说话，你得用指标衡量它的好坏。我常用的评估指标有四个：

• SHD（Structural Hamming Distance）：预测图与真实图之间不同边的数量，越小越好。
• TPR（True Positive Rate）：真实存在的边中被正确预测出来的比例。
• FDR（False Discovery Rate）：预测出来的边中实际上不存在的比例。
• FPR（False Positive Rate）：真实不存在的边中被错误预测出来的比例。

很多博主只报SHD，这是不够的。SHD只衡量边集合的差异，但对"边方向错了"和"边完全不存在"这两种错误一视同仁，容易掩盖模型的具体短板。我的习惯是：同时报告SHD和FDR，如果SHD低但FDR很高，说明模型学到了很多假边，只是方向恰好蒙对了；如果SHD高而TPR低，说明模型偏保守，漏边严重。这两个维度一起看，才能判断一个算法到底靠不靠谱。

6. 从线性到非线性演进：这套"无环连续化"思路的生命力

DAGs with NO TEARS的价值远不止这个线性版本本身。它提供了一个"框架性工具"——任何结构学习问题，只要你能把'无环'约束写成光滑函数，就可以用连续优化去做。后来这个思路被大量后续工作吸收：

• 有的工作在非线性场景下用神经网络建模SEM，把权重矩阵替换成逐变量的网络参数，同时保留h(W)约束；
• 有的工作把同样技巧扩展到针对时序数据的动态因果模型；
• 有的工作针对大规模稀疏图设计更高效的优化策略；
• 还有的把它与变分推断结合，做贝叶斯版本的结构学习。

它们本质上都在吃NO TEARS的红利——无环约束的连续化。

我也拿NOTEARS和传统的PC算法、GES算法做过比较。在光滑的线性高斯数据上，NOTEARS的速度明显占优；在样本量极小（n=50）的情况下，PC算法有时靠条件独立性检验也能取得不错的结果，而NOTEARS则需要小心处理过拟合。在实际业务场景里，如果你假设数据生成机制比较接近线性SEM、且有一定样本量，NOTEARS通常是一个非常省心的默认选项。

7. 我对自己实操中最惊喜的一个细节

最后分享一个让我印象深刻的实操体会。有一次在模拟数据上，我故意加入了两条强相关但非因果的链路（变量A和B都受同一个隐含变量驱动，但彼此没有直接关系），传统的相关性分析很容易把这种虚假关联学成边。但NOTEARS训练完，h(W)约束天然限制了"环"的出现，配合稀疏阈值化，这类混杂结构带来的假边很大程度上被自动削掉了。虽然它在理论上并没有显式建模隐变量，但DAG约束带来的"全局结构性"确实让模型对局部虚假相关有一定的天然的免疫力——这也是我在落地项目中比较倾向于先跑一遍NOTEARS的原因。

如果你也想在自己的数据上尝试结构学习，建议先跑通上面那段核心伪代码，再逐渐替换损失函数、尝试非线性扩展。这套方法的学习曲线并不长，但带来的视角转换——从"离散图搜索"到"连续约束优化"——很可能会改变你对整个结构学习问题的理解方式。

所谓"NO TEARS"，从另一个角度讲，也是劝大家别为组合爆炸掉眼泪：换一个视角，问题可能就有解了。