3-因果发现代码结构

概述与定位

一句话定位: 因果发现算法的代码实现遵循标准架构模式，理解核心函数有助于快速实现和调试算法。

** prerequisites **

因果推断基础（PC算法、GES算法）
概率图模型基础（DAG、贝叶斯网络）
线性代数基础（矩阵运算、优化）
Python编程基础（numpy, scipy）

核心直觉

因果发现算法分为三大类：基于约束的方法（PC算法）、基于分数的方法（GES算法）、基于功能因果的方法（LiNGAM）。虽然实现细节各异，但核心代码架构相似——都涉及条件独立性测试、图结构优化和DAG约束。

代码模块化的核心直觉：因果发现算法可以分解为独立的模块——CI测试模块负责图骨架学习，分数计算模块负责结构优化，DAG约束模块负责确保结果是有向无环图。这种模块化设计便于调试和比较不同算法。

关键函数的核心直觉：每个算法都有几个核心函数决定其性能和适用范围。理解这些函数的作用和实现，是深度掌握因果发现算法的关键。

数学推导

结构方程模型

线性非高斯无环模型（LiNGAM）：

其中是带约束的系数矩阵，是潜在独立成分。

NOTEARS优化目标

其中：

是加权邻接矩阵估计
是L1正则项
是DAG约束

DAG约束的数学形式

使用矩阵指数确保非奇异：

GES评分函数

贝叶斯信息准则（BIC）评分：

对于高斯线性模型：

代码架构

核心模块结构

# causal_discovery/
# ├── core/
# │   ├── ci_test.py          # 条件独立性测试
# │   ├── score_func.py       # 评分函数
# │   ├── dag_constraint.py  # DAG约束
# │   └── structure_search.py # 结构搜索
# ├── methods/
# │   ├── pc.py              # PC算法实现
# │   ├── ges.py             # GES算法实现
# │   ├── noteals.py        # NOTEARS实现
# │   └── lingam.py          # LiNGAM实现
# └── utils/
#     ├── graph_utils.py     # 图工具函数
#     └── eval_metrics.py    # 评估指标

条件独立性测试模块

# ci_test.py
def ci_test(X, Y, S, method='partial_corr', alpha=0.05):
    """
    Conditional Independence Test
    
    Parameters:
    -----------
    X : array, variable X
    Y : array, variable Y
    S : array, conditioning set
    method : str, 'partial_corr' or 'gp'
    
    Returns:
    --------
    p_value : float, p-value of CI test
    """
    if len(S) == 0:
        # Marginal independence (Pearson correlation)
        corr, p_value = pearsonr(X, Y)
    else:
        # Partial correlation
        corr, p_value = partial_corr(X, Y, S)
    return p_value

DAG约束模块

# dag_constraint.py
import numpy as np
 
def dag_constraint(W, n):
    """
    Check if W corresponds to a DAG
    h(W) = tr(e^{W◦W}) - n = 0
    """
    M = np.linalg.matrix_power(np.eye(n) + W * W, n)
    h = np.trace(M) - n
    return h
 
def notears_loss(W, X, lambda1):
    """
    NOTEARS objective with DAG constraint
    """
    N, n = X.shape
    loss = 0.5 / N * np.sum((X - X @ W) ** 2)
    penalty = lambda1 * np.sum(np.abs(W))
    constraint = dag_constraint(W, n)
    return loss + penalty, constraint

核心函数实现

PC算法骨架学习

def pc_skeleton(data, alpha=0.05):
    """
    Learn skeleton of causal graph using PC algorithm
    
    Returns:
    --------
    skeleton : dict, {pair: separator_set}
    """
    n_vars = data.shape[1]
    skeleton = {(i, j): set() for i in j in range(n_vars)}
    
    # Iterate over separation set sizes
    for sep_set_size in range(n_vars):
        for (i, j) in skeleton.keys():
            if i == j:
                continue
            for S in combinations(range(n_vars) - {i, j}, sep_set_size):
                p_value = ci_test(data[:, i], data[:, j], S)
                if p_value > alpha:
                    # i and j are conditionally independent
                    del skeleton[(i, j)]
                    del skeleton[(j, i)]
                    skeleton[(i, j)] = set(S)
                    break
    return skeleton

GES搜索

def ges_score(data, parent_set):
    """
    Compute BIC score for a given parent set
    
    BIC = -N/2 * log|S_G| - d/2 * log N
    """
    N, n = data.shape
    d = len(parent_set)
    
    if d == 0:
        cov = np.cov(data[:, parent_set + []])
        score = -N / 2 * np.log(np.linalg.det(cov + 1e-6))
    else:
        # Fit linear model and compute residual covariance
        residual = fit_linear_model(data, parent_set)
        score = -N / 2 * np.log(np.linalg.det(np.cov(residual.T) + 1e-6))
    
    score -= d / 2 * np.log(N)
    return score

NOTEARS优化

def notears_optimize(X, lambda1=0.1, max_iter=100):
    """
    NOTEARS algorithm for DAG learning
    """
    N, n = X.shape
    W = np.zeros((n, n))
    
    for iteration in range(max_iter):
        # Compute loss and gradient
        loss, grad = notears_loss(W, X, lambda1)
        
        # Add DAG constraint gradient
        M = np.linalg.matrix_power(np.eye(n) + W * W, n)
        h_grad = M.T @ (2 * W)
        
        # Update W
        W -= 0.01 * (grad + h_grad)
        
        # Check convergence
        h = dag_constraint(W, n)
        if h < 1e-6:
            break
    
    return W

依赖库

核心库

库名	用途	安装命令
numpy	矩阵运算	`pip install numpy`
scipy	统计测试	`pip install scipy`
networkx	图操作	`pip install networkx`

因果发现专用库

库名	用途	安装命令
gcastle	主流因果发现库	`pip install castle`
bnlearn	贝叶斯网络	`pip install bnlearn`
causallearn	因果学习工具	`pip install causallearn`

复现注意事项

关键实现点

CI测试实现: PC算法的核心是高效的CI测试
DAG约束数值稳定性: 矩阵指数计算需要数值稳定技巧
优化收敛判定: NOTEARS需要准确判断DAG约束满足
方向确定规则: V-结构确定边方向需要正确实现

常见问题

内存爆炸: 高维数据CI测试组合爆炸
数值精度: 矩阵求逆在高维时数值不稳定
优化局部最优: GES可能陷入局部最优

超参数参考

参数	典型值	注意事项
alpha (CI阈值)	0.01-0.1	太小遗漏边，太大多假边
lambda1 (NOTEARS)	0.1-1.0	控制稀疏性
max_iter	100-1000	视数据规模而定

本章要点总结

因果发现代码分模块：CI测试、评分、DAG约束、搜索
PC算法核心是条件独立性测试
GES使用BIC评分引导结构搜索
NOTEARS使用连续约束转化为优化问题
DAG约束 $◦$ 确保无环
gcastle/bnlearn/causallearn是主流实现库
CI测试效率是PC算法瓶颈
数值稳定性在高维实现中关键
GES需要正确实现插入-删除操作
LiNGAM基于独立成分分析识别因果方向

关键词

因果发现代码结构 PC算法 GES算法 NOTEARS LiNGAM CI测试 DAG约束结构学习代码实现

Huajishe Notes

Explorer