在前面的章节中，我们学习了如何在已知因果结构的情况下进行因果推断。但在许多实际场景中，因果结构本身是未知的。我们需要从观察数据中自动发现变量之间的因果关系——这就是因果发现（Causal Discovery）的任务。本章将介绍主流的因果发现算法，包括基于约束的方法、基于分数的方法和连续优化方法，并探讨如何评估发现的因果结构的质量。
因果发现的目标是从观察数据 $\mathcal{D} = {X^{(1)}, …, X^{(n)}}$ 中学习变量集合 $\mathbf{X} = {X_1, …, X_p}$ 之间的因果图结构 $G$。与传统的关联分析不同，因果发现试图揭示数据生成过程的真实机制，而不仅仅是统计相关性。
这个任务面临几个根本性挑战：
马尔可夫等价类：多个不同的DAG可能对应相同的条件独立性关系。例如，链结构 $X \rightarrow Y \rightarrow Z$ 和 $X \leftarrow Y \leftarrow Z$ 在观察数据下无法区分，除非我们有额外的假设或干预数据。
隐变量问题：未观测的混杂因素可能导致虚假的因果关系。当两个变量 $X$ 和 $Y$ 都被未观测变量 $U$ 影响时，即使 $X$ 和 $Y$ 之间没有直接因果关系，它们也会表现出相关性。
因果充分性假设：我们通常需要假设观测到了所有相关变量。违反这个假设会导致错误的因果结论。然而，在实际应用中，完全满足因果充分性几乎是不可能的。
忠实性假设：DAG中的条件独立性关系需要在数据分布中得到体现。参数的特殊取值可能导致额外的独立性，使得因果结构无法唯一确定。
样本复杂度：随着变量数量增加，可能的DAG数量呈超指数增长。对于 $p$ 个变量，可能的DAG数量为： \(|DAG_p| = \sum_{k=0}^{p(p-1)/2} (-1)^k \binom{p(p-1)/2}{k} 2^{p(p-1)/2-k} \cdot p!\)
例如，3个变量有25个可能的DAG，10个变量就有超过 $4.2 \times 10^{18}$ 个可能的DAG。
统计功效问题：有限样本下的假设检验可能缺乏统计功效，特别是在高维情况下需要进行大量的条件独立性测试。
因果发现的一个核心概念是马尔可夫等价类（Markov Equivalence Class, MEC）。理解MEC对于正确解释因果发现算法的输出至关重要。
定义：两个DAG $G_1$ 和 $G_2$ 属于同一个马尔可夫等价类，当且仅当它们蕴含相同的条件独立性关系集合。
等价性判定定理（Verma and Pearl, 1990）： 两个DAG属于同一个马尔可夫等价类，当且仅当它们具有：
例子1：链结构的等价类
(a) A → B → C (b) A ← B ← C (c) A ← B → C
这三个DAG具有相同的骨架A-B-C，且都没有v-结构
它们属于同一个马尔可夫等价类
条件独立性：A ⊥ C | B
例子2：v-结构的识别
(d) A → C ← B (e) A → C → B (f) A ← C → B
↓
D
(d)包含v-结构 A → C ← B（A和B不相连）
(e)和(f)没有v-结构
这三个属于不同的马尔可夫等价类
CPDAG（完成部分有向无环图）： 马尔可夫等价类可以用CPDAG表示，其中：
本质图（Essential Graph）： CPDAG也称为本质图，它编码了从观察数据可以确定的所有因果信息。计算本质图的算法：
成功的因果发现依赖于几个关键假设。理解这些假设的含义和局限性对于正确应用因果发现算法至关重要。
假设1：因果马尔可夫条件（Causal Markov Condition） 每个变量在给定其直接原因（父节点）的条件下，与其非效应（非后代节点）条件独立。
形式化表述：对于DAG $G = (V, E)$ 和联合分布 $P$， \(X_i \perp\!\!\!\perp \text{NonDesc}(X_i) | \text{Pa}(X_i)\)
这个假设意味着DAG捕获了所有相关的因果关系，没有遗漏重要的直接因果联系。
假设2：忠实性假设（Faithfulness Assumption） 数据分布 $P$ 对DAG $G$ 是忠实的，如果 $P$ 中的所有条件独立性关系都对应于 $G$ 中的d-分离关系，反之亦然。
| 形式化：$X \perp!!!\perp Y | Z \text{ in } P \Leftrightarrow X \perp_d Y | Z \text{ in } G$ |
忠实性假设排除了参数的”病态”取值。例如，在路径 $X \rightarrow Z \leftarrow Y$ 上，如果两条路径的效应恰好相消，可能导致 $X \perp Y$，违反忠实性。
假设3：因果充分性（Causal Sufficiency） 不存在未观测的共同原因（混杂因素）同时影响两个或多个观测变量。
违反因果充分性的后果：
假设4：无选择偏差（No Selection Bias） 数据采样过程不依赖于变量之间的因果关系。如果数据收集受到某些变量值的影响（如只观察幸存者），会引入选择偏差。
假设5：无测量误差 观测变量准确反映了真实的潜在变量。测量误差可能掩盖或扭曲真实的因果关系。
基于约束的方法（Constraint-based methods）通过系统地测试条件独立性关系来构建因果图。这类方法的核心思想是：因果结构蕴含特定的条件独立性模式，通过识别这些模式可以推断因果结构。
PC算法是最经典的基于约束的因果发现算法，由Peter Spirtes和Clark Glymour于1991年提出。它通过渐进式地测试条件独立性，从数据中学习因果图的等价类。
算法概览： PC算法分为两个主要阶段：骨架学习（发现哪些变量之间有边）和方向确定（确定边的方向）。
阶段1：骨架学习
骨架学习采用自顶向下的策略，从完全图开始逐步删除边：
算法：PC-骨架学习
输入：变量集V，显著性水平α
输出：无向图骨架G，分离集Sep
1. 初始化：G = 完全无向图
2. l = 0 // 条件集大小
3. while 存在节点的邻居数 ≥ l do
4. for 每对相邻节点 (Xi, Xj) in G do
5. for 所有大小为l的条件集 S ⊆ Adj(Xi)\{Xj} ∪ Adj(Xj)\{Xi} do
6. if 条件独立性测试 Xi ⊥ Xj | S (p-value > α) then
7. 删除边 Xi - Xj
8. Sep(Xi, Xj) = S // 记录分离集
9. break
10. l = l + 1
11. return G, Sep
优化技巧：
阶段2：方向确定
方向确定分为两步：识别v-结构和应用传播规则。
算法：PC-方向确定
输入：骨架G，分离集Sep
输出：CPDAG
1. // 识别v-结构
2. for 每个三元组 (Xi, Xj, Xk) where Xi-Xj-Xk in G and Xi与Xk不相连 do
3. if Xj ∉ Sep(Xi, Xk) then
4. 定向为 Xi → Xj ← Xk // v-结构
5.
6. // 应用Meek规则
7. repeat
8. 应用规则R1-R4（见下文）
9. until 没有新的边被定向
10. return CPDAG
Meek规则（完整版）：
时间复杂度分析：
FCI算法是PC算法的扩展，由Spirtes等人于1999年提出。它放松了因果充分性假设，允许存在隐变量和选择偏差，是处理不完全观测系统的重要工具。
PAG（部分祖先图）表示： FCI输出PAG而非DAG，使用四种端点类型：
边的类型及含义：
FCI算法流程：
算法：FCI主流程
输入：数据D，显著性水平α
输出：PAG
1. // 阶段1：初始骨架发现
2. 运行PC算法的骨架学习阶段，得到初始骨架G和分离集Sep
3.
4. // 阶段2：识别v-结构
5. 对所有不相邻的Xi, Xk，如果存在Xj使得Xi-Xj-Xk：
6. if Xj ∉ Sep(Xi, Xk) then
7. 标记为v-结构：Xi *→ Xj ←* Xk
8.
9. // 阶段3：可能的d-分离（Possible d-sep）
10. repeat
11. for 每对不相邻的节点 (Xi, Xj) do
12. 找到所有可能的d-分离路径PDS(Xi, Xj)
13. for S ⊆ PDS(Xi, Xj) do
14. if Xi ⊥ Xj | S then
15. 更新Sep(Xi, Xj) = S
16. 如果Xi-Xk-Xj且Xk ∉ S，标记新的v-结构
17. until 没有新的独立性被发现
18.
19. // 阶段4：最终方向确定
20. 应用完整的FCI定向规则R0-R10
21. return PAG
FCI定向规则（部分）：
PC-Stable和FCI-Stable变体： 原始PC/FCI算法的骨架可能依赖于变量顺序。Stable版本通过以下修改确保顺序无关：
条件独立性测试是基于约束方法的核心组件。测试的质量直接影响因果发现的准确性。
测试框架： 零假设 $H_0: X \perp!!!\perp Y | Z$ 备择假设 $H_1: X \not\perp!!!\perp Y | Z$
连续变量的测试方法：
偏相关测试（线性高斯假设）： \(\rho_{XY|Z} = \frac{\rho_{XY} - \rho_{XZ}\rho_{YZ}}{\sqrt{(1-\rho_{XZ}^2)(1-\rho_{YZ}^2)}}\)
| Fisher’s Z变换：$Z = \frac{1}{2}\ln\frac{1+ | \rho | }{1- | \rho | }$ |
| 在零假设下：$\sqrt{n- | Z | -3} \cdot Z \sim \mathcal{N}(0,1)$ |
优点：计算快速，理论成熟 缺点：仅适用于线性关系
核条件独立性测试（KCIT）： 使用再生核希尔伯特空间（RKHS）中的统计量： \(\text{KCIT} = n \cdot \text{tr}(K_X K_Y K_Z^{-1})\)
其中 $K_X, K_Y, K_Z$ 是相应的核矩阵。
优点：可检测非线性关系 缺点：计算复杂度高 $O(n^3)$
条件互信息测试： \(I(X;Y|Z) = \mathbb{E}_{p(x,y,z)}\left[\log\frac{p(x,y|z)}{p(x|z)p(y|z)}\right]\)
使用k近邻或核密度估计来估算互信息。
优点：理论优美，无分布假设 缺点：高维时估计不准确
离散变量的测试方法：
条件 $\chi^2$ 测试： \(\chi^2 = \sum_{x,y,z} \frac{(O_{xyz} - E_{xyz})^2}{E_{xyz}}\)
| 自由度：$df = ( | X | -1)( | Y | -1) | Z | $ |
G-test（似然比测试）： \(G = 2\sum_{x,y,z} O_{xyz} \ln\frac{O_{xyz}}{E_{xyz}}\)
渐近分布：$G \sim \chi^2_{df}$
条件互信息测试： \(\hat{I}(X;Y|Z) = \sum_{x,y,z} \frac{n_{xyz}}{n} \log\frac{n \cdot n_{xyz}}{n_{xz} \cdot n_{yz}}\)
混合变量的处理：
多重假设检验校正： PC算法需要进行大量的独立性测试，需要控制总体错误率：
基于分数的方法（Score-based methods）将因果发现表述为优化问题：在所有可能的DAG空间中搜索使某个评分函数最大化的结构。这类方法的优势在于可以直接比较不同的因果假设，并自然地处理模型复杂度与拟合优度的权衡。
评分函数需要满足几个关键性质才能保证因果发现的一致性：
理想性质：
BIC分数（贝叶斯信息准则）：
BIC在似然和模型复杂度之间寻求平衡： \(\text{BIC}(G) = \log P(\mathcal{D}|G, \hat{\theta}_G) - \frac{|G|}{2}\log n\)
其中：
| $ | G | $ 是模型参数数量 |
对于线性高斯模型： \(|G| = p + |E| + p = 2p + |E|\) （$p$ 个截距，$|E|$ 个回归系数，$p$ 个噪声方差）
分解形式： \(\text{BIC}_i(Pa_i) = -\frac{n}{2}\log(RSS_i/n) - \frac{|Pa_i|+2}{2}\log n\)
其中 $RSS_i$ 是节点 $i$ 对其父节点回归的残差平方和。
BDeu分数（贝叶斯Dirichlet等价uniform）：
专门用于离散变量，基于Dirichlet先验的边际似然：
\[\text{BDeu}(G) = \log P(\mathcal{D}|G) = \sum_{i=1}^{p} \sum_{j=1}^{q_i} \left[\log\frac{\Gamma(\alpha_{ij})}{\Gamma(\alpha_{ij}+n_{ij})} + \sum_{k=1}^{r_i} \log\frac{\Gamma(\alpha_{ijk}+n_{ijk})}{\Gamma(\alpha_{ijk})}\right]\]其中：
等价样本大小（ESS）选择： \(\alpha_{ijk} = \frac{\text{ESS}}{r_i \cdot q_i}\)
常用 $\text{ESS} \in {1, 10, 50}$，ESS=1 对应无信息先验。
BGe分数（贝叶斯高斯等价）：
用于连续高斯变量，基于Normal-Wishart先验：
\[\text{BGe}(G) = \sum_{i=1}^p \text{BGe}_i(Pa_i)\]其中每个局部分数为： \(\text{BGe}_i(Pa_i) = \log\frac{\Gamma_p(\frac{\nu+n}{2})}{\Gamma_p(\frac{\nu}{2})} + \frac{\nu}{2}\log|T| - \frac{\nu+n}{2}\log|T+S_{Pa_i}|\)
参数说明：
混合变量的评分函数： 对于包含连续和离散变量的系统，可以使用：
GES算法在马尔可夫等价类空间中进行贪婪搜索：
前向阶段（Forward Phase）：
后向阶段（Backward Phase）：
GES算法的优势：
对于小规模问题（$p \leq 25$），可以使用动态规划找到全局最优解：
\[\text{score}^*(U) = \max_{X_i \in U} \left\{ \max_{Pa_i \subseteq U \setminus \{X_i\}} \text{score}_i(Pa_i) + \text{score}^*(U \setminus \{X_i\}) \right\}\]这种方法保证找到最优DAG，但计算复杂度为 $O(p \cdot 2^p)$。
近年来，将组合优化问题转化为连续优化问题成为因果发现的新趋势。
NOTEARS（NO TEARS - 无环约束的平滑表示）将DAG学习表述为连续优化问题：
关键创新：无环约束的代数表示 \(h(W) = \text{tr}(e^{W \circ W}) - p = 0\)
其中 $W$ 是邻接矩阵，$\circ$ 表示Hadamard积。当且仅当 $W$ 对应无环图时，$h(W) = 0$。
优化问题： \(\min_W \frac{1}{2n}\|X - XW\|_F^2 + \lambda \|W\|_1\) \(\text{s.t. } h(W) = 0\)
使用增广拉格朗日方法求解： \(\mathcal{L}_{\rho}(W, \alpha) = f(W) + \alpha h(W) + \frac{\rho}{2}h(W)^2\)
加性噪声模型（ANM）： \(X_j = f_j(Pa_j) + N_j\)
其中 $f_j$ 是非线性函数，$N_j$ 是独立噪声。
识别性结果：在温和条件下，如果 $f$ 是非线性的，则因果方向是可识别的（除了高斯噪声+线性函数的情况）。
RESIT算法（Regression with Subsequent Independence Test）：
因果VAE（Causal VAE）： 将因果结构编码到变分自编码器中：
编码器：$q_{\phi}(Z|X)$ 解码器：$p_{\theta}(X|Z)$ 因果层：$Z_i = f_i(Pa_i^Z) + \epsilon_i$
目标函数： \(\mathcal{L} = \mathbb{E}_{q_{\phi}(Z|X)}[\log p_{\theta}(X|Z)] - \text{KL}(q_{\phi}(Z|X)||p(Z)) + \lambda \cdot \text{DAG\_penalty}\)
评估因果发现算法的性能是一个挑战，因为真实的因果结构通常是未知的。
结构汉明距离（SHD）： 编辑操作（添加、删除、翻转边）的最小数量： \(\text{SHD} = \text{FP} + \text{FN} + \text{边方向错误数}\)
精确率与召回率：
结构干预距离（SID）： 衡量两个因果图在干预分布上的差异： \(\text{SID}(G, G') = \sum_{i,j} \mathbb{I}[\text{Pa}_j^{do(X_i)}(G) \neq \text{Pa}_j^{do(X_i)}(G')]\)
当可以进行干预实验时，可以更准确地评估因果结构：
合成数据：
半合成数据：
真实数据：
华为在部署5G网络时面临复杂的故障诊断挑战。5G网络包含数万个相互依赖的组件，当出现性能下降或服务中断时，快速准确地定位根本原因至关重要。传统的基于规则的故障诊断系统难以处理复杂的依赖关系和新型故障模式。
监控指标：
挑战：
阶段1：离线因果图构建
使用历史数据构建基础因果图：
数据预处理：
1. 时间序列对齐和重采样（5分钟粒度）
2. 异常值检测和处理
3. 特征工程（滑动窗口统计量）
算法选择：
- PC-Stable算法：处理高维稀疏图
- 时间序列扩展：考虑时滞效应
- 领域知识约束：禁止某些边（如用户指标→硬件指标）
阶段2：在线根因定位
当检测到异常时，使用因果图进行根因分析：
输入：异常KPI集合 A = {KPI_1, ..., KPI_k}
输出：根因节点集合 R
算法：
1. 在因果图中标记所有异常节点
2. 计算每个节点的责任分数：
RS(v) = P(A|do(v=abnormal)) - P(A)
3. 使用反向追溯找到最上游的异常节点
4. 验证：模拟修复该节点，预测其他节点恢复情况
特征工程：
因果强度评估： 使用部分相关系数衡量直接因果强度： \(\rho_{ij|S} = \frac{\text{Cov}(X_i, X_j | S)}{\sqrt{\text{Var}(X_i|S) \cdot \text{Var}(X_j|S)}}\)
定量结果：
业务价值：