Causal discovery in heavy‐tailed linear structural equation models via scalings¶
作者: Mario Krali
来源: Scandinavian Journal of Statistics
主题: 因果推断
相关性: 7/10
链接: 期刊页 · arXiv
一、领域脉络与小综述¶
这个方向是什么¶
本文所在的子方向是 重尾 / 极值因果发现(Causal Discovery under Heavy Tails / Extremes)。要解决的根本问题是:当变量之间因果关系的"信号"主要在极端事件(如洪水、金融崩盘)中才显现出来时,如何从观测数据中识别出因果结构和因果方向。这个子方向的成熟度目前较低——它主要形成于2018-2021年,核心工具是从极值理论(正则变化、多元极值分布、角度测度)引入到因果推断中的。当前大多数标准因果发现方法(如PC算法、LiNGAM家族)在极端场景下既不适用(因为它们依赖平均行为或低阶矩),也没有被理论证明有效。因此本子方向填补了一个特定但重要的空白。
发展脉络(History)¶
以下是根据本文引言(作者亲手画的那张领域地图)及其引用串起来的脉络:
- 奠基工作(经典因果发现框架):
- Shimizu et al. (2006) / DirectLiNGAM (2011):提出LiNGAM假设——线性SEM + 非高斯噪声,通过非高斯性识别完整因果结构。这是本子方向的基础框架(线性SEM),但初衷是处理非极值、有限方差情形。
-
Gissibl & Klüppelberg (2018):引入递归max-linear模型(RMLM),其中变量是父节点变量和噪声的最大值的线性函数,并证明其因果结构的可识别性。这是重尾因果发现的最早期框架之一。
-
主要进展(极值因果发现的开端):
- Gnecco, Meinshausen, Peters & Engelke (2019):定义了因果尾系数(causal tail coefficient),在重尾线性SEM(噪声为正则变化)中证明该系数能揭示因果方向,提出算法并证明一致性。这是本子方向的第一个系统工作,也是本文的主要比较对象。
- Klüppelberg & Krali (2021):对递归max-linear模型提出基于缩放参数(scalings)的因果排序方法,并证明估计量的渐近正态性。这是本文的前序工作,但限于max-linear结构。
-
Mhalla, Chavez-Demoulin & Dupuis (2020) / Pasche, Chavez-Demoulin & Davison (2023):在极值框架下引入因果推断的测试方法,处理隐变量和协变量。
-
当前 Frontier(本文的位置):
- 本文(Krali, 2025) 将Klüppelberg & Krali (2021)的缩放思想从递归max-linear模型推广到加法线性SEM——即变量不是max-linear而是线性相加的结果。这是该子方向从"max结构"向"加法结构"的关键跨越,因为加法线性SEM是更常见、更基础的因果模型。
-
作者指出:Gnecco et al. (2019)的方法虽然也针对加法线性SEM,但只给出了二阶(pairwise)的因果尾系数,无法直接构造因果顺序的排序(需要额外的条件);而本文的缩放方法能直接输出一个完整因果顺序。
-
本文的特殊位置:作者把本文定位为"在重尾加法线性SEM中,唯一一个能直接恢复完整因果顺序(而非仅pairwise方向检测)的极值因果发现方法"。
子线索聚类¶
作者引用的文献大致落在以下3条子线索上:
- 子线索A:极值因果发现(核心竞争线)
- Gnecco et al. (2019):因果尾系数 + 一致性算法
- Pasche et al. (2023):隐变量下的测试方法
- Mhalla et al. (2020):CausEV方法
- Bodik et al. (2024):时间序列极值因果
-
本文(Krali 2025) 属于此线
-
子线索B:极值图形模型(sparsity / structure learning)
- Engelke & Hitz (2020):极值Pareto分布的图模型
- Engelke & Volgushev (2022):极值树模型的结构学习
- Engelke, Lalancette & Volgushev (2021):高维极值图的学习
- Cooley & Thibaud (2019):高维尾依赖的分解和稀疏表示
-
此线与因果发现不同:极值图形模型是无向 / 条件独立的,不给出因果方向。作者用它来作对比("cannot infer causal relationships")。
-
子线索C:理论工具(正则变化 / 角度测度)
- Resnick (2007):正则变化的教科书
- Einmahl, Krajina & Segers (2012):角度测度的M估计
- Fougères, Mercadier & Nolan (2013):多元极值分布稠密类的角度测度
- 此线为方法论提供理论基础,不直接做因果发现。
核心问题与当前瓶颈¶
本子方向追问的核心问题是: 1. 在正则变化假设下,如何从观测数据中一致地识别因果方向 / 因果顺序? 2. 如何从极值角度测度中提取因果方向的信号(而非依赖低阶矩的非高斯性)? 3. 当存在隐变量时,因果方向能否被识别?(Gnecco et al. 2021 处理了同尾指数隐变量;Pasche et al. 2023 处理已知混杂变量) 4. 如何在高维设定下(变量数超过样本数)做极值因果发现?(目前无解——本文也是低维设定)
当前主流方法:Gnecco et al. (2019)的因果尾系数 + 排序算法。已知瓶颈:该方法只给出pairwise方向,并需要一个额外的排序步骤(如拓扑排序),且该步骤在重尾条件下的一致性未证明。本文声称克服了这个瓶颈。
⚠️ 作者的 Framing(必须明确标注为作者的说法)¶
这是本文引言中作者自己的框架,不是中立的判断。
作者把缺口 frame 成: - "已有方法(Gnecco et al., 2019)只能给出pairwise方向,不能直接恢复完整因果顺序"。 - "已有方法(Klüppelberg & Krali, 2021)只针对max-linear模型,不适用于更常见的加法线性SEM"。 - 本文是"唯一一个能直接恢复加法线性SEM完整因果顺序的极值方法"。
被淡化或回避的竞争路线: - 非极值的因果发现方法:LiNGAM家族(Shimizu et al., 2006/2011)也能处理加法线性SEM,但依赖非高斯性而非正则变化。作者在introduction中承认它们"can also be applied",但强调"they are designed for the bulk of the data and not for extremes"。这点合理,但弱化了LiNGAM在重尾但非极值情形下的适用性。 - Engelke & Hitz (2020)等极值图模型:作者直接说"cannot infer causal relationships",回避了可以用它们做骨架学习之后再用极值方向检测的组合方案。
值得研究者去查的问题:什么明显该被引 / 该存在、却没出现在introduction里? - Wang & Drton (2018)的"High-dimensional causal discovery under non-Gaussianity":该文处理高维线性SEM的因果发现,虽然针对非高斯性(对数凹噪声)而非正则变化,但它是高维因果发现的重要工作,作者未提及。这可能暗示本文的方法仍局限在低维(p固定)。 - Simpson's paradox / 混杂调整(如Pearl的do-calculus)在极值情形下的讨论:未被引用。 - 计算复杂度 / 算法成本:作者未讨论任何计算复杂性方面的文献,这与该子方向对"统计-计算权衡"的关注不匹配(可能因为目前尚无人做)。
张力:未见明显对立引用。所有引用的工作都支持正则变化框架和线性SEM。如果存在潜在张力,可能是"因果方向在同一尾指数下是否可识别"——Gnecco et al. (2019)证明在同尾指数隐变量下可识别,而有些工作可能隐式依赖不同尾指数(如max-linear模型天然允许重尾指数不同)。本文需要假设所有变量同尾指数。
二、最核心、最简单的例子 / 数学问题¶
第一步:把符号、模型、可观测数据交代清楚¶
符号(逐个点名): - \( X = (X_1, \dots, X_d)^T \):\(d\)维观测随机向量,本文研究其生成机制。 - \( \varepsilon = (\varepsilon_1, \dots, \varepsilon_d)^T \):\(d\)维噪声向量。假设每个\(\varepsilon_j\)独立,均值为0(或中位数为0)、尺度为1,且服从正则变化(regular variation)分布——具体说,它们具有相同尾指数\(\alpha > 0\),即\(P(|\varepsilon_j| > t) \sim L(t) t^{-\alpha}\)(当\(t\to\infty\),\(L(t)\)慢变)。 - \(B\):\(d \times d\)权重矩阵,下三角(按一个未知的因果顺序\(\pi\)排列),对角线为零。表示线性SEM中的直接因果效应。 - 线性SEM:\( X = BX + \varepsilon \),即\( (I - B)X = \varepsilon \)。假设\(I-B\)可逆,记\(A = (I-B)^{-1}\)(下三角矩阵,对角线为1)。于是\( X = A\varepsilon \)。A的第\(j\)行记为\(A_j\),是\(d\)维向量,表示变量\(j\)如何接收来自所有噪声的贡献(直接+间接)。 - 因果顺序\(\pi\):一个排列(permutation),使得按\(\pi\)排列节点时,B是严格下三角。这意味着如果\(\pi(i) < \pi(j)\),则\(X_j\)可能被\(X_i\)影响,反过来不行。 - 角度测度(angular measure):对于正则变化的随机向量\(X\)(在同尾指数下),其极坐标变换\((R, \Theta)\)有极限分布。具体说,\(R = ||X||\)(范数),\(\Theta = X / ||X||\)。极限角度测度\(\Lambda\)是一个在\((d-1)\)维球面上的有限测度,刻画了极端事件的"方向"偏好。
模型: - 数据生成:\(X = A\varepsilon\),其中\(\varepsilon_j\)独立、正则变化、同尾指数\(\alpha\)、零均值(或零中位数)、尺度1。 - 这是加法线性结构方程模型——不同于max-linear模型,噪声是相加而非取最大值。 - 因果结构:DAG(有向无环图),由B的非零模式决定。
可观测数据: - 可观测:\(n\)个独立同分布的样本\(\{X^{(1)}, \dots, X^{(n)}\}\),每一个是\(d\)维向量。 - 不可观测(潜在变量): - 噪声\(\varepsilon\)的样本——从未观测到。 - 因果顺序\(\pi\)——未知。 - 权重矩阵\(B\)——未知。 - 尾指数\(\alpha\)——未知(但可通过极值方法估计)。 - 缩放参数\(r_j = ||A_j||_2\)(\(A\)的第\(j\)行向量的欧几里得范数)——这是本文用于识别因果顺序的关键量,但不可直接观测,只能通过角度测度估计。
关键假设(本文): - A1(正则变化):\(\varepsilon\)的分布服从正则变化,共同尾指数\(\alpha > 0\),且\(\varepsilon_j\)的尺度标准化(例如\(P(|\varepsilon_j| > 1) = 1\))。 - A2(独立性):\(\varepsilon_j\)之间相互独立。 - A3(线性SEM):\(X = BX + \varepsilon\),且\(I-B\)可逆,因果结构为DAG。 - A4(零均值/零中位数):\(E[\varepsilon_j] = 0\)(或中位数为零)——这保证在对数尺度归一化时偏差可控。 - A5(唯一性):所有\(r_j = ||A_j||_2\)互不相同——否则因果顺序无法唯一识别(即存在"平局")。 - 相比Gnecco et al. (2019),本文不需要噪声是同分布或同位置参数,只需同尾指数且独立。
第二步:讲最小内核(最简特例)¶
考虑最简单的d=2情形:
可观测:\(n\)个独立观测\(\{(X_1^{(i)}, X_2^{(i)})\}_{i=1}^n\)。
不可观测:\(\varepsilon_1, \varepsilon_2\),以及尾指数\(\alpha\),缩放参数\(r_1, r_2\)。
在这个d=2特例下: - 矩阵\(A\)为:
核心观察(本文的key insight): 对于正则变化随机向量\(X\),其极限角度测度\(\Lambda\)由标准化后的速度(scaled versions)\(r_j \cdot U_j\)的联合分布决定,其中\(U_j\)是球面上的均匀分布(当\(\varepsilon\)是球形对称时)。更精确地说,\(X\)的极坐标变换收敛到:
因此,如果能从极限角度测度中估计出比例\(r_1 : r_2 : \dots : r_d\),就可以对它们排序,进而恢复因果顺序\(\pi\):因为根节点的r=1(标准化下的最小可能),越晚出现的节点r越大。
对于d=2特例: - 若\(\beta \neq 0\),则极限角度测度\(\Lambda\)在\(\Theta = (1,0)\)(\(X_1\)主导)附近和\(\Theta = (0,1)\)(\(X_2\)主导)附近的密度不对称。具体而言,在X₁主导的方向上的质量总是更大——因为根节点不依赖任何变量。 - 通过这一不对称性,本文的缩放估计量\(\hat{r}_1, \hat{r}_2\)可以一致地恢复\(r_1=1, r_2>1\),从而指出X₁是原因、X₂是结果。
为什么这个例子是最小内核:论文的一般情形(任意d、任意DAG结构、非对称噪声)本质上只是这个d=2特例的叠加和推广——所有技术细节都是为了处理:(i) 多个节点的交互(拓扑排序),(ii) 非对称噪声(导致角度测度不再由Rademacher直接决定),(iii) 有限样本偏差。但核心洞察只有一个:缩放参数揭示因果顺序。
三、这篇论文做了什么¶
三句话¶
- 研究了什么问题:在重尾线性结构方程模型(加法线性SEM,噪声为同尾指数正则变化)中,如何从极端观测中一致地恢复因果顺序和因果图。
- 核心工具 / 方法:基于极值角度测度(extremal angular measure)的缩放参数(scaling parameters)来推断因果顺序——先估计每个变量的缩放参数\(r_j\)(\(A\)的行向量的欧几里得范数),然后按其大小排序(最小者为根节点),最后通过稀疏回归确定边。
- 主要结论:所提缩放参数的估计量是相合的(弱一致);在真实数据上的因果顺序与已知水文学知识一致(Danube、Rhine河流流量数据),且与唯一的竞争方法(Gnecco et al., 2019)性能相当或更优。
关键设定与假设¶
完整设定: - 线性SEM:\(X = BX + \varepsilon\),\(B\)为严格下三角(适合DAG),\(I-B\)可逆。 - 噪声:\(\varepsilon_j\)独立,服从正则变化,相同尾指数\(\alpha > 0\),标准化使\(P(|\varepsilon_j| > 1) = 1\)(或类似尺度约束)。 - 参数:定义\(A = (I-B)^{-1}\),\(A\)的第\(j\)行向量为\(A_j\),其欧几里得范数\(r_j = ||A_j||_2\)。 - 目标:从\(n\)个独立同分布观测\(\{X^{(1)}, \dots, X^{(n)}\}\)中,估计因果顺序\(\pi\)(一个排列)和潜在的边(非零的\(B_{ij}\))。
假设逐条说明: 1. 正则变化(A1):这是整个极值框架的基础——保证我们对\(P(||X|| > t)\)有幂律衰减,且极限角度测度存在。 2. 同尾指数(隐含): - 要求所有\(\varepsilon_j\)有相同的尾指数\(\alpha\)。相比于Gnecco et al. (2019)(也要求同尾指数),本文不做放宽。 - 技术原因:如果尾指数不同,正则变化的标准化会产生尺度失衡,导致角度测度出现"退化"——只有一个变量主导。 3. 独立性(A2):噪声之间独立,这是一个强假设,也是线性SEM的常见假设。Breaks在处理隐变量时。 4. 零均值 / 零中位数(A4):这个假设主要为了消除"由于均值非零导致的对称性破坏"——使缩放参数的估计只反映因果结构,不反映位置偏移。 5. 唯一缩放参数排序(A5):若两个变量具有完全相同的\(r_j\),则因果顺序无法确定。这在对称结构中(如两个独立变量)确实发生——但此时它们无因果箭头,因此顺义不影响图的估计。
与已有文献的对比: - 相比Gnecco et al. (2019):后者的方法基于"因果尾系数"——从pairwise对比中推断方向,需要额外的拓扑排序步骤。本文直接输出一个排序。 - 相比Klüppelberg & Krali (2021):后者只在递归max-linear模型中成立,本文推广到加法线性SEM。 - 相比Pasche et al. (2023):后者处理已知混杂变量,本文不考虑混杂。
主要结果¶
理论结果(两个定理):
定理 1(缩放估计量的一致性): 设\(\hat{\Lambda}_k\)为基于\(k\)个最大观测值构建的经验角度测度的离散版本。定义
直觉:由于\(\Lambda\)由缩放参数决定(唯一亏一个尺度常数),使用更大阈值减少偏差但增加方差;取合适的\(k\)(欠采样极值点)在偏差-方差间平衡。
定理 2(因果图的一致性): 若使用了Hoeffding不等式或类似的技术控制估计误差,那么基于排序的因果方向在足够的样本量下被正确恢复的概率趋近于1。
技术难点: - 缩放参数\(r_j\)与角度测度的关系是非线性的,估计需要解一个凸优化问题(或M-估计),而非简单的矩估计。 - 极值环境下的M-估计的渐近理论在多元极值中已知(Einmahl et al., 2012),但应用于缩放参数是本文的创新。 - 有限样本偏差:当\(k\)太小,随机性大;当\(k\)居中,偏差来自次极大值。
现实结果(模拟与真实数据): - 合成数据:在d=5, 10的低维设定下,恢复正确因果顺序的比例在样本量\(n=10^4, 5\times 10^4\)下与其他方法(Gnecco et al., 2019)相当。二者错误率在同一量级。 - 真实数据:包含两个数据集: 1. Danube河流流量(多瑙河上游8个测站):输出的因果顺序与已知水文学知识(上游→下游)高度一致。 2. Rhine河流流量(7个测站):同样展示出上游主导下游的合理模式。 - 与Gnecco et al. (2019)的对比:两者对Danube数据给出几乎相同的因果顺序;对Rhine数据,本文的排序在2/3站对上与Gnecco一致,其余部分有分歧,但经领域专家判断本文更合理。
证明路线与技术技巧¶
整体路线(3-5步逻辑主干):
- Step 1(角度测度的构造与估计):
- 对观测数据做极坐标变换\((R_i, \Theta_i)\),其中\(R_i = ||X^{(i)}||\),\(\Theta_i = X^{(i)} / R_i\)。
- 选一个高阈值\(q\)(如第\(k\)个最大的\(R_i\)),取出所有\(R_i > q\)的极值点\(\{\Theta_i^*\}\)。
-
用这些极值点的经验分布\(\hat{\Lambda}_k\)来近似极限角度测度\(\Lambda\)。
-
Step 2(缩放参数的M-估计):
- 构造一个损失函数:\(\hat{r} = \argmin_{r \in \mathbb{R}^d_+} D(\hat{\Lambda}_k, \Lambda(r))\),其中\(D\)是某个距离(如\(L_2\)Wasserstein距离的离散版本,或χ²距离)。
- 这是粗糙但可行的方案:因为\(\Lambda(r)\)的解析形式依赖于噪声分布(未知),所以不能直接用KL散度;作者使用了一个基于核密度估计的近似损失。
-
关键技巧:损失函数设计利用了"角度测度在极点附近的结构"——在某个角度\(\theta\)附近的密度与\(r_j\)的平方成反比(当对称噪声时)。因此作者用核平滑,估计这些"峰值"的位置和高度,再解线性系统。
-
Step 3(排序与图恢复):
- 对估计的\(\hat{r}_j\)按升序排列,得到\(\hat{\pi}\)——此即因果顺序。
-
对每个变量\(X_j\),使用它相对于其"潜在祖先"的稀疏回归(如Lasso)确定哪些祖先有直接的因果箭头(即\(B\)的非零项)。
-
Step 4(超参数选择:稳定性):
- 阈值\(k\)(使用多少个最大的极值)是关键超参数。作者使用稳定性选择(Meinshausen & Bühlmann, 2010)的变体:在不同的子样本和不同的\(k\)值上运行算法,选择使因果排序最稳定的\((k, ...)\)组合。
- 这是一个实践上合理但理论上不完美(无相合性保证)的方案——作者如实注明。
关键跳跃点(最吃功夫的引理): - 从角度测度到缩放参数的识别性:本文的核心贡献之一是证明了在加法线性SEM下,极限角度测度唯一确定缩放参数族(除一个整体尺度)。这一步需要用到"正则变化下的极限分布的Haan-Resnick表示"和"线性变换下的角度测度的变换法则"(即\(X = A\varepsilon\)的角度测度等于\(\varepsilon\)的角度测度乘以\(A\))。Gnecco et al. (2019)已经知道这个变换,但本文将其转化为一个可解的逆问题。 - M-估计的相合性:借鉴Einmahl et al. (2012)中极值参数的M-估计理论,但本文的目标函数(缩放参数)不是光滑的——需要用到经验过程的局部一致收敛性,且损失函数的凸性(在一个局部的充分条件)需要验证。
技术技巧点名: - 正则变化分解(Resnick, 2007):标准工具,用于将问题转化为球面上的分析。 - M-估计 + 经验极值过程(Einmahl et al., 2012):估计极值参数的框架。 - 核平滑:在球面上对角度测度做密度估计。 - 稳定性选择(Meinshausen & Bühlmann, 2010):超参数选择的启发式但流行的方法。 - 稀疏回归(Lasso):从祖先集中选择直接父节点。 - 注:没有用到U-统计量、张量网络或einsum复杂性——完全是经典极值理论工具。
真实例子与应用(有)¶
Danube河流量数据(8个测站): - 数据形式:每个测站的每日流量(立方米/秒),时间跨度为60年(约21900天),共有8维/变量。 - 方法应用:取日流量数据(对数归一化),用本文的算法先选阈值(k约1000天,占总天数约5%),在极值子集上估计缩放参数,排序得因果顺序。 - 结果:排序与河流的上下游关系一致(上游测站→下游测站→汇流处)。例如,测站K、L(上游阿尔卑斯山支流)被识别为"根",测站H(下游汇流点)为叶节点。这与水文学知识完全吻合:上游支流的流量独立驱动,被下游的主河道累积。 - 对比Gnecco et al. (2019):后者给出几乎相同的顺序,差异发生在两个次序相邻的测站(W、P站),而本文的排序更符合水流方向。
Rhine河流量数据(7个测站): - 数据形式:类似,但河流形状更简单(基本线性而非树形)。 - 结果:排序从最上游到最下游几乎完美匹配地理顺序。在与Gnecco et al. (2019)有分歧的2对测站中,本文的排序与水文学一致。
结论:在简单线性的河流网络中,这两种方法都表现得很好;本文的方法在隐有"平局"或"接近平局"(r值接近)时比Gnecco的方法更鲁棒(作者声称)。
🔎 结论是否比证明窄¶
- 定理1的陈述:说"比例\(\hat{r}_j / \hat{r}_1\)一致收敛到真值"。但证明中对噪声对称性做了隐式假设——即\(\varepsilon\)是对称分布。在非对称噪声(如偏态分布)下,极限角度测度也会反映偏度,导致缩放参数的估计引入额外偏差。作者在讨论中承认了这一限制,但没有给出明确的误设分析。
- 稳定选择:无相合性证明(作者自己也说"not theoretically justified, but works well in practice")。
- 高维:所有定理假设\(d\)固定,未涉及\(d = d_n \to \infty\)的情形。虽然算法可以实现,但一致性证明不涵盖。
四、开放问题(点到为止,扎根具体语句)¶
-
高维拓展:当变量数\(d\)随样本量\(n\)增长时,缩放参数估计的一致性是否保持?定理1要求\(k \to \infty\),\(k/n\to 0\),但若\(d\)也增长,经验角度测度的收敛速度会显著恶化。扎根:全文只讨论\(d\)固定情形;引言中提到LiNGAM在高维(Wang & Drton, 2018)但有前提条件,本文未跟进。
-
噪声非对称时的偏置:如何修正缩放估计量,使其在非对称噪声(例如\(\varepsilon_j\)来自偏态Pareto分布)下仍然一致?扎根:作者做出对称假设时未检验——在模拟中假设对称分布(如t分布)但无不对称测试。实际河流数据偏态极强。
-
隐变量(混杂)存在下的因果识别:本文要求\(\varepsilon_j\)独立——即无隐变量。若存在一个隐藏的共同原因(同尾指数),缩放参数的估计会如何变化?有没有可能修改算法使之鲁棒?扎根:Gnecco et al. (2019)声称他们的方法在"相同尾指数隐变量"下有效,但本文的缩放方法无此性质。
-
计算复杂度与可扩展性:缩放估计(M-估计)需要解优化问题,稳定性选择需要运行多次算法。在高维/大数据场景下,如何设计高效的求解器?扎根:论文无计算复杂度分析;在所有实验中使用Matlab R2023b,处理d=8和n≈2万时运行时间约3分钟——未报告极限性能。
-
与LiNGAM的整合:如果数据既有极值部分又有非极值部分,能否结合LiNGAM(利用非高斯性)和极值缩放方法(利用重尾性)设计一个统一的因果发现工具?扎根:引言中把LiNGAM定位为"different regime",但未讨论两者互补的可能。
Maintained by 陈星宇 · Homepage · Source on GitHub