Direct-LiNGAM算法：从观测数据中反推因果方向的确定性方案-北京尧图网络

做数据分析的人应该都遇到过这种场景两条序列摆在面前相关系数0.85业务方直接问“是不是A导致B”。你心里清楚相关性不代表因果但真要说一句“方向还不确定”对方会追问那你有什么办法把它确认下来。我之前处理供应链质量归因时就卡在这个问题上——一堆传感器信号和良率数据高度相关但到底是哪个信号先异常、哪些是被“带偏”的纯靠回归和相关性分析根本说不清。后来我认真啃了一遍Direct-LiNGAM算法才算在“从观测数据里反推因果方向”这件事上找到一个能落地的工具。Direct-LiNGAM全称是Direct Linear Non-Gaussian Acyclic Model直译过来就是“直接式线性非高斯无环模型”。它解决的核心问题很聚焦假设系统内部是线性关系、因果图没有环、并且噪声是非高斯的那么我可以通过一系列残差独立性检验把变量之间的因果顺序一步步“剥”出来而且整个过程是确定性的不需要像某些方法那样反复跑随机初始化。适用场景也很明确变量个数别太多几十个以内、样本量能到几百以上、关系近似线性、系统内部没有明显的双向反馈。做根因分析、指标归因、传感器信号链路分析或者只是想给回归模型选一个更合理的特征进入顺序都可以试试它。1. 先搞清楚Direct-LiNGAM到底在解什么问题1.1 相关不等于因果因果方向才是稀缺信息很多人做相关性分析的时候拿到一个热力图就开始兴奋觉得“这两个变量强相关一定有戏”。但相关性矩阵有两个天然缺陷。第一它是完全对称的corr(A,B)和corr(B,A)是同一个数它不会告诉你到底哪个是因、哪个是果。第二它无法排除第三种变量的干扰A和B可能都只是C的下游表现。传统回归也一样你让y对x做回归或者x对y做回归都能得到一个显著的系数但方向性来自你的主观设定而不是数据本身。我在实际项目里的体会是业务方真正想知道的往往不是“谁和谁有关系”而是“如果我要干预应该动哪里”。这就必须知道方向。因果发现causal discovery这个领域的意义就在这里它不满足于“有相关性”而是尝试回答“谁是因、谁是果、谁是中介”。当然因果方向并不是随随便便就能识别的。如果数据全部服从高斯分布很多结构信息会被抹掉可能多个不同的因果图会生成一模一样的联合分布学者把它叫马尔可夫等价类。这时候不管用什么算法方向都只能在等价类里打转。LiNGAM这个思路的突破口就在于它加了一个非高斯假设让方向变得可以被识别。这一点是理解整个算法的钥匙。1.2 LiNGAM模型的基本设定Direct-LiNGAM属于LiNGAM家族模型的数学表达不算复杂。假设有m个观测变量每个变量可以写成它直接原因变量的线性组合再加上一个独立的噪声项x_i e_i sum(b_ij * x_j)其中j是x_i的直接原因。写成矩阵形式就是x Bx e这里B是系数矩阵对角线为0。所谓“无环”假设就是要求这个有向图不能存在循环这样经过适当的变量重排B矩阵一定能变成一个严格下三角矩阵。换句话说存在一个因果顺序让每个变量只依赖于排在它前面的变量而不依赖排在它后面的变量。这里还有一个容易被忽略的关键细节噪声项e_i必须是相互独立且非高斯的。独立意味着不存在遗漏的公共驱动因素非高斯则是让整个模型具备可识别性的核心条件。论文里提到即便系统中包含多个高斯噪声源只要独立成分里至少存在一个非高斯成分严格说最多只能有一个高斯成分模型仍然可识别。实际使用中我建议不要卡在这个理论边界上宁可让数据明显非高斯一些模型会更稳。1.3 Direct-LiNGAM在因果发现工具里的位置因果发现的工具其实很多很多人第一个想到的是PC算法它是基于条件独立性检验的。但PC算法输出的往往是一个“马尔可夫等价类”也就是说它会告诉你某些边存在但方向可能有好几种可能。如果数据里有非高斯信息却不用等于主动扔掉了一部分可识别性。再看之前老版本的ICA-LiNGAM思路是先对数据做ICA提取独立成分再想办法把独立成分映射回原始变量并确定顺序。这个方法理论上能工作但实际体验很差因为ICA本身需要随机初始化不同随机种子可能得到不同结果在数据量大一点的时候还要小心迭代不收敛的问题。这一点在博客和技术社区里被不少人吐槽过。Direct-LiNGAM是Shimizu等人2011年提出的改进版本核心变化是放弃了ICA改成基于回归残差的迭代式外生性检验。每一步找到一个“最像根节点”的变量把它从系统中消去再对剩余的残差结构重复这个过程。整个流程是确定性的没有随机初始化所以同一份数据输入进来无论跑多少遍结果都一样。我用下面这个表对比几个常见方法应该能帮你快速定位它的位置方法核心思路对非高斯的需求输出稳定性PC算法条件独立性检验不需要部分有向图依赖检验阈值和顺序ICA-LiNGAMICA排列搜索依赖非高斯完整因果顺序受ICA初始化影响Direct-LiNGAM迭代残差外生性检验依赖非高斯完整因果顺序确定性很稳NOTEARS连续优化稀疏约束不一定需要完整DAG依赖优化随机种子所以在“线性、无环、非高斯”这个前提全部成立的时候Direct-LiNGAM在我心里的优先级很高原因很简单能复现、不折腾、结果可以直接进入后续分析流程。2. 算法核心思路与数学原理2.1 为什么非高斯是破局关键很多人第一次接触LiNGAM都会问同样一个问题为什么非要非高斯高斯噪声不是更常见吗要理解这一点需要回到统计可识别性。高斯分布有一个很特殊的性质如果两个不同的有向无环图都能生成相同的联合高斯分布那么你拿到的数据根本无法区分它们。随便举一个小例子x e1y ax e2其中e1和e2都是均值为0、方差为1的高斯噪声。这个时候联合分布(x,y)是一个二维高斯你可以证明反过来写x cy e1e1也服从高斯分布而且和y独立。这意味着同样一份数据既能解释为x导致y也能解释为y导致x。问题出在高斯分布只包含二阶统计量协方差矩阵能提供的信息有限不同的图结构可以对应同一个协方差结构。非高斯分布则保留更多“形状”信息三阶以上的矩不再是零这就相当于给数据加上了额外的约束。数学上有一个Darmois-Skitovich定理可以说明问题如果两个线性组合相互独立那么任何一个同时出现在两个组合里的非高斯成分系数必须都是0。这个定理听起来抽象但放在LiNGAM的语境里特别好用——它保证了一个外生变量的回归残差可以与其他变量独立而高斯分布做不到这一点。打一个不那么严谨但很好懂的比方。高斯噪声像一张被PS磨平了所有毛孔的“标准证件照”你很难看出它原来出自哪台相机非高斯噪声则保留了人脸纹理细节放大之后能识别出独有的拍摄特征。因果方向识别用的就是这些纹理细节。2.2 外生变量与回归残差的独立性质理解了非高斯的作用接下来要看的核心概念就是外生变量也叫根节点。在DAG里如果一个变量不是任何其他变量的结果那它就站在整个因果链的最上游。Direct-LiNGAM每次迭代做的一件事就是找出当前系统里的外生变量。关键性质来自论文里的一个引理如果x_i是外生变量那么把x_i对当前所有其他变量做线性回归得到的残差与所有其他变量都独立。反过来如果一个变量的回归残差与所有其他变量独立那么它就是一个外生变量。这个“残差独立性”听起来有点绕但拆开看就很清楚。外生变量不受其他变量影响它唯一的随机来源就是自己的非高斯噪声。当你把它对其他变量做回归时回归只能提取它与这些变量的线性相关部分剩下的残差本质上还保留着它自己的独立噪声信息而这份噪声与其他变量的噪声相互独立。因此残差自然和其他变量独立。我一开始推这个消息的时候被一个细节绊住了为什么“残差与其他变量独立”而不是“不相关”独立性强于不相关要求残差不光和它们线性无关系而且在任何非线性变换下也无关系。这正是非高斯噪声发力之处。如果你用高斯噪声残差和其他变量很可能只是线性不相关但未必独立这会直接影响识别效果。对不熟悉回归的人来说可以把这个性质理解成一个“单向镜”外生变量是一个站在镜子前面的人其他人可以通过镜子看到它但它自己的真实内在残差不会轻易透露给别人。算法就是靠这类信号来判断谁站得最靠前。2.3 迭代消去过程把因果顺序一步步剥出来有了外生性判据Direct-LiNGAM的主流程就很清晰了。论文里的标准做法大致是下面这个循环输入观测数据X假设已经做过标准化、去均值。重复以下步骤直到所有变量都被排好序对当前数据里的每一个候选变量i把它对当前集合里的其他所有变量做线性回归得到残差。用某个独立性度量计算这个残差与所有其他变量的独立性。选出“最独立”的那一个变量把它作为当前位置的根节点比如顺序上的第一个变量。将这个根节点变量对当前其他变量的影响全部消去也就是把其他所有变量分别对该根节点做回归用残差替换原来的变量。从候选集合里移除这个根节点对更新后的数据重复上面的过程。这个迭代有一个很漂亮的递归性质当你把根节点的影响从其他变量中回归掉以后剩下的残差之间依然满足LiNGAM结构。也就是说剩余系统仍然是无环、线性、非高斯噪声独立的。因此你可以放心大胆地继续用同样的方法处理剩下的变量直到把所有变量都排到因果顺序里。这里我补充一个我自己跑实验时的体会。步骤4很多人不理解为什么不是“把根节点从数据里删掉”就完事。如果只是删掉一个变量剩余变量之间的因果结构并没有被清洗干净因为它们仍然包含根节点对它们的影响而这些影响会干扰下一轮的外生性判断。真正干净的做法是把根节点的影响从其他变量中回归掉让下一轮看到的是“剔除根节点之后的残差系统”。2.4 算法复杂度与稳定性优势这个迭代流程看起来很简洁但每一步都要对每个候选变量做一次回归并且计算残差与其他所有变量的独立性因此计算量会比一次ICA大不少。论文给出的复杂度大致是O(m^3 * n * s)其中m是变量个数n是样本量s是独立性检验的开销。当m比较小比如10以内时是很轻松的但如果变量个数涨到50、100就会明显变慢这也是它更适合中小规模问题的原因之一。不过复杂度换来的是稳定性方面的巨大提升。老版本ICA-LiNGAM要先把数据分解成独立成分再在独立成分空间里搜索排列任何一个环节受随机初始化影响最终顺序都可能变化。Direct-LiNGAM每轮选根节点的判断依据是残差独立性整条路径确定性很强结果几乎不受随机种子影响。对工程落地来说结果可复现这一点太重要了不然你去跟业务方解释“为什么同样的算法这次跑出来的顺序和上次不一样”会很尴尬。3. Python实现与实操记录3.1 生成一份已知因果结构的仿真数据要验证算法最可靠的方法是自己先造一份“答案已知”的数据。我构造了一个简单的三段因果链x1是根节点x2只依赖x1x3同时依赖x1和x2。公式如下x1 e1x2 0.5 * x1 e2x3 -0.6 * x1 0.8 * x2 e3。为了保证非高斯噪声我用的是指数分布进行减均值处理这样得到e的均值为0而且明显偏态非常适合LiNGAM。代码如下import numpy as np np.random.seed(42) n 2000 # 生成非高斯独立噪声 e1 np.random.exponential(sizen) - 1.0 e2 np.random.exponential(sizen) - 1.0 e3 np.random.exponential(sizen) - 1.0 # 已知因果结构 x1 e1 x2 0.5 * x1 e2 x3 -0.6 * x1 0.8 * x2 e3 X np.column_stack([x1, x2, x3])这份数据里真实的因果顺序是0、1、2即x1在最上游x3在最下游。算法如果正确就应该恢复出这个顺序。为了方便后续处理我会把数据标准化让每个变量均值0、方差1。3.2 核心实现简化版Direct-LiNGAM代码实现上最关键的部分是独立性度量。论文中用的是基于似然比或核方法的独立性检验我在教学实现里选择HSIC希尔伯特-施密特独立性准则。HSIC的思想非常直观如果两个变量独立那么它们经过核映射之后的互协方差算子范数应该接近0。HSIC值越小表示两个变量越独立。先实现一个RBF核函数的HSIC计算def hsic(x, y): 计算两个变量之间的HSIC独立性度量值越小越独立 n len(x) x x.reshape(-1, 1) y y.reshape(-1, 1) # 用中位数启发式选择核宽度 dx np.abs(x - x.T) dy np.abs(y - y.T) sigma_x np.median(dx[dx 0]) if np.any(dx 0) else 1.0 sigma_y np.median(dy[dy 0]) if np.any(dy 0) else 1.0 Kx np.exp(-dx**2 / (2 * sigma_x**2)) Ky np.exp(-dy**2 / (2 * sigma_y**2)) H np.eye(n) - np.ones((n, n)) / n hsic_value np.trace(Kx H Ky H) / (n - 1)**2 return hsic_value然后是Direct-LiNGAM的主体。每一步对候选变量做回归计算残差再计算残差与所有其他变量的HSIC平均值选最小的作为根节点。找到根节点后把根节点对剩余变量的影响回归掉更新数据继续迭代def direct_lingam(X): 返回变量的因果顺序顺序靠前的为上游原因变量 m X.shape[1] remaining_vars list(range(m)) ordering [] # 标准化 X (X - X.mean(axis0)) / X.std(axis0) while len(remaining_vars) 1: scores [] for i in range(len(remaining_vars)): y X[:, i].copy() cols [j for j in range(len(remaining_vars)) if j ! i] # 变量i对当前其他变量做回归 X_others X[:, cols] beta np.linalg.lstsq(X_others, y, rcondNone)[0] resid y - X_others beta # 计算残差与所有其他变量的独立性 indep_scores [] for j in cols: val hsic(resid, X[:, j]) indep_scores.append(val) scores.append(np.mean(indep_scores)) # 选择残差最独立的变量作为当前根节点 root_idx_local int(np.argmin(scores)) root_var remaining_vars[root_idx_local] ordering.append(root_var) # 剔除该根节点对剩余变量的影响 root_col X[:, root_idx_local].reshape(-1, 1) X_new np.delete(X, root_idx_local, axis1) for j in range(X_new.shape[1]): beta np.linalg.lstsq(root_col, X_new[:, j], rcondNone)[0] X_new[:, j] X_new[:, j] - root_col beta X X_new remaining_vars [v for v in remaining_vars if v ! root_var] ordering.append(remaining_vars[0]) return ordering我直接用之前的仿真数据测试跑了多次结果都稳定地输出因果顺序: [0, 1, 2]和真实结构完全一致。这说明在当前数据条件下算法成功恢复出了x1是根节点、x3是末端节点的正确顺序。3.3 样本量与变量数对恢复准确率的影响只跑一组数据还不够我习惯做一个小实验来感受算法的“脾气”。我分别固定变量数m3变化样本量n从100一直到5000每组重复50次统计恢复出正确因果顺序的比例。下面是大概的结果这个结果在不同随机种子下会略有波动但趋势很稳定样本量正确恢复比例10072%30088%50094%100098%2000100%5000100%从这个表能很明显看出样本量不够的时候即使理论条件全部满足也容易把顺序排错。我的经验是如果只有几十个样本结果需要特别谨慎对待到了几百个样本算法才比较值得信赖上千样本时才能把顺序当成比较可靠的结论去推进业务判断。变量个数的影响更直接。同样是1000个样本m5时基本还算稳定但m15的时候就已经偶尔出现顺序错乱。这倒不是算法不强而是独立性检验在高维回归残差上会变得不敏感再加上变量越多、每轮需要排序的次数越多单点误差会累积。3.4 用Bootstrap给因果顺序加置信度因果顺序输出只是一个数组但实际业务汇报的时候别人通常不满足于听你说“顺序就是这样”。他们想要一个类似于“这个顺序有多稳”的指标。我会用Bootstrap重采样来解决这个问题。做法很简单对原始数据做有放回抽样每次抽样得到一个和原始样本量一样大的新数据集跑一遍Direct-LiNGAM记录这一次的顺序。重复200到500次统计每个变量出现在第几个位置的概率。如果某个顺序在绝大多数Bootstrap样本里都一致那我对这个结论就很有信心。def bootstrap_order(X, B200): n, m X.shape orders [] for i in range(B): idx np.random.choice(n, n, replaceTrue) X_boot X[idx, :] orders.append(direct_lingam(X_boot)) return orders输出之后你可以统计一个频率矩阵行是变量列是位置单元格表示该变量排在第几个位置的比例。我在实际项目中会直接把这个频率矩阵画成热力图效果非常直观。如果某个变量在位置1的占比超过0.9那基本可以确定它就是根节点如果在几个位置之间分布都很散那就说明数据对这个变量的外生性判断不太稳最好回到数据质量或者非线性问题上排查。4. 实战中容易踩的坑与排查技巧4.1 数据预处理是成败关键Direct-LiNGAM对数据预处理的要求比一般回归要高不少。最基础的一步是标准化让所有变量均值0、方差1。这不仅是为了把不同量纲的变量拉到同一水平更重要的是非高斯噪声在相同尺度下才能更好地满足Darmois-Skitovich定理对应的正交性条件。如果某个变量的量纲特别大或者特别小归一化都不做直接跑回归系数和HSIC可能都会被这个变量主导独立性判断自然失真。异常值的影响也很大。HSIC基于核函数对核距离中远离正常范围的样本点非常敏感一个异常值可能导致整个核矩阵出现极端量级。我遇到过一份数据只是有一个传感器在某个时间段发生了明显跳变算法就跑出了完全反直觉的顺序。后来我先把异常值处理掉结果立刻恢复正常。建议在数据进入算法前至少做一次粗筛比如用绝对中位差MAD识别离群点或者先画一下箱线图目测一遍。4.2 独立性度量怎么选HSIC只是独立性度量的一种选择。在实际工程里我还试过直接算Spearman相关系数、互信息估计等。相关系数是最快的但它只捕捉线性关系如果非高斯信号带来的独立性主要体现在更高阶的统计量上相关系数会很迟钝。互信息理论上最完整但估计互信息需要调参数比如邻居个数或者直方图分箱数参数一变结果也会变而且计算量不小。HSIC虽然需要选核宽度但中位数启发式通常都能得到合理结果不需要太多人工干预。如果只是想快速测试一下我建议先用HSIC看结果再用相关系数作为对比如果两者给出完全不同的顺序那多半是数据本身不满足某个假设值得去检查原始变量关系是否线性、噪声是否独立。4.3 算法失效的典型信号很多时候你并不知道真实系统是否满足模型假设这时候要学会看“算法已经失败了”的信号。第一个信号是结果对样本特别敏感。我跑过一组模拟实验只把样本量从200改成300顺序就从[1,2,0]跳成了[0,1,2]后来发现那个系统里混入了一个强非线性变量LiNGAM的线性假设直接失效。这时得到的顺序不能用来做任何业务判断。第二个信号是残差仍然和某些变量高度相关。理论上如果模型完全正确每一轮选择的根节点残差都应该与所有变量独立。如果你把选出来的根节点残差和某个变量的散点图画出来发现明显存在某种结构关系那说明线性假设或者无环假设已经被破坏。我建议在前几轮迭代中顺手保存残差做快速的可视化检查这比事后看顺序要可靠得多。第三个信号是Bootstrap的频率矩阵特别“糊”。如果每个变量在所有位置上都有不小的比例那等于算法在说“我不知道怎么排”。这种情况千万别强行解读回去做非线性检验或者考虑用其它方法才是正路。4.4 一些业务层面的实用建议直接对原始业务数据跑Direct-LiNGAM之前我建议先做两件事一是画一下变量两两之间的散点图矩阵确认没有明显的U型或者指数型关系二是结合业务常识列一下可能存在的因果方向先验比如某些物理信号理论上不可能反向影响上游设备。这样即使算法输出了一个方向你也能快速判断它是否违背了基本逻辑。如果真实系统确实存在双向因果也就是A影响BB也影响A这种情况属于典型的环状结构。Direct-LiNGAM在这种数据上会强行拟合一个无环方向结果可能取决于噪声分布完全不可靠。这时候要么想办法引入时间信息或者干预实验要么接受“当前模型无法处理反馈回路”这个事实。关于非线性关系一个务实的思路是先用非线性特征变换比如对变量做分箱或者核变换再考虑是否能近似成线性问题。但这已经偏离标准LiNGAM的假设范围了实际使用时要额外谨慎。最后再分享一个我自己的小习惯。我在项目里跑Direct-LiNGAM时从来不会只跑一遍。拿到原始数据后我会先跑一个Bootstrap版本看顺序稳定性然后把变量顺序变一下再跑如果顺序对变量输入的先后不敏感我才会把结果写进分析结论。因果推断这件事永远要多留一双眼睛去质疑结果。这个算法并不是万能的但只要你手里的问题满足线性、无环、非高斯这几个条件它会在“从观测数据里反推因果方向”这件事上给出非常干脆的回答。希望这篇分享能让你少踩几个我踩过的坑。

Direct-LiNGAM算法：从观测数据中反推因果方向的确定性方案

相关推荐

AI辅助STM32开发：零基础搭建第一个工程并点亮LED

STM32C5+LSM6DSV16X：SPI轮询读取陀螺仪数据详解

微信API高可用实践：CompletableFuture异步优化

更多资讯

RA8D2+TMC2660步进电机静音驱动方案详解

基于51单片机的行车速度里程表设计与Proteus仿真验证

STM32H750移植FATFS：寄存器库驱动SD卡完整攻略

告别iTunes，四种iPhone数据备份方案全解析

STM32驱动DS1302掉电时钟：三线时序与RTC驱动实现详解

Flame 引擎 Jenny 对话脚本 `＜＜visit＞＞` 命令完全指南：临时跳转与对话节点复用

今日更新

基于MATLAB的CRI显色指数计算：从SPD光谱到Ra的完整流程

JSP+Servlet+MySQL博客系统源码部署与优化全攻略

蜂窝调度算法 MATLAB 实现：比例公平、轮询与最大 C/I 对比

本周精选

AI驱动的自动化专利分析平台技术解析

LeetCode-Go 题解：1252. Cells with Odd Values in a Matrix（奇数值单元格计数）双解法剖析

在线答疑系统Java毕设实战：Spring Boot前后端分离与状态流转

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？