利用高阶累积量学习含潜在混杂变量的最稀疏线性因果DAG

# 通过高阶累积量学习存在潜在混淆变量的最稀疏线性因果DAG 来源: https://arxiv.org/html/2607.05984 京都大学，日本\\NameHisayuki Hara\\Emailhara\.hisayuki\.8k@kyoto\-u\.ac\.jp \\addrInstitute for Liberal Arts and Sciences，京都大学，日本 ###### 摘要 在带有潜在混淆变量的线性非高斯无环模型（LvLiNGAM）中恢复精确的有向无环图（DAG）仍然是一个具有挑战性的问题。尽管LvLiNGAM仅能识别到观测等价类，但每个等价类都由唯一的最稀疏DAG刻画。然而，从有限样本中恢复最稀疏DAG仍然困难。现有方法虽然渐近一致，但它们没有提供恢复唯一最稀疏DAG的显式有限样本程序，也无法处理带有任意数量潜在混淆变量的模型。在本文中，我们提出了一种有限样本方法，无需对潜在混淆变量的数量施加任何限制即可恢复最稀疏DAG。仿真研究和真实数据分析表明，与现有方法相比，所提方法在有限样本下性能更优。 ###### 关键词: 因果发现; DAG; LiNGAM; 潜在混淆变量; 累积量; ## 1 引言 线性非高斯无环模型（LiNGAM）为因果发现提供了强大的框架（Shimizu2006; Shimizu2011）。在没有潜在变量的情况下，LiNGAM能够完全识别因果DAG。然而，在许多实际应用中，潜在混淆变量是不可避免的。hoyer2008estimation引入了带有潜在变量的LiNGAM（LvLiNGAM），并证明了任何LvLiNGAM都可以转化为一个标准模型，其中所有潜在变量相互独立且在因果顺序上先于观测变量。他们还假设潜在变量的数量先验已知，通过超完备独立成分分析（OICA；例如，eriksson2004identifiability）估计混合矩阵。然而，由于依赖于OICA，这种方法容易收敛到局部最优（shimizu14aBayesianestimation）。为了避免依赖OICA，已经提出了几种通过残差独立性测试来估计标准LvLiNGAM的方法，例如 Pairwise LvLiNGAM（Entner2011）、ParceLiNGAM（tashiro2014parcelingam）、重复因果发现（RCD）（Maeda2020; Maeda2022; maeda2022rcd）和BANG（wang2023BANG）。然而，这些方法都无法完全识别形成弓形结构（wang2023BANG）的观测变量之间的祖先关系或父子关系。当存在弓形结构时，chen2024identification使用观测变量的累积量来识别双变量设置中单个潜在混淆变量的祖先关系。在此基础上，chen2025identification将方法扩展到观测双变量情况下多个潜在混淆变量的情形。schkoda2024causal提出了ReLVLiNGAM，一种基于累积量的递归方法，适用于多个观测变量和潜在混淆变量。无需依赖OICA或预先知道潜在变量的数量，ReLVLiNGAM能够恢复标准LvLiNGAM的观测等价类。在下面描述的一般性假设下，观测等价类内的最稀疏DAG被唯一确定。最稀疏DAG在其模型类内是一般性的，而同一观测等价类中任何更稠密的DAG则需要非一般的参数值。这为将最稀疏DAG视为观测等价类的标准代表提供了自然的合理性。尽管ReLVLiNGAM能够一致地估计混合矩阵，但它没有提供从有限样本中估计最稀疏DAG的程序。此外，原始的ReLVLiNGAM要求，在每次迭代中，一个没有观测父节点的观测变量（以下称为观测源）具有比其观测兄弟节点更少的潜在父节点；否则，累积量更新和混合矩阵估计可能会失败（图1（https://arxiv.org/html/2607.05984#S2.F1）；参见附录C（https://arxiv.org/html/2607.05984#A3））。我们将此问题称为ReLVLiNGAM的局部限制。在本文中，我们开发了一种有限样本算法，用于在标准LvLiNGAM的观测等价类中恢复最稀疏DAG，且不受局部限制限制。所提方法建立在ReLVLiNGAM的自上而下框架之上。它从观测源开始，逐步推断观测变量之间的父子关系，以恢复最稀疏DAG。在此框架内，我们引入了两个关键创新。首先，我们引入了一种直接对观测变量进行残差化的更新规则，而不是像ReLVLiNGAM那样递归更新高阶累积量。更新观测变量而不是其累积量，减轻了高阶累积量估计中误差的递归传播，从而提高了有限样本性能。此外，这种更新规则消除了ReLVLiNGAM中对局部限制的需要，扩展了方法的适用性。其次，我们引入了一个顺序程序，用于从估计的祖先结构中识别每个观测源与其后代之间精确的父子关系。该程序能够直接从有限样本中恢复最稀疏DAG。我们的主要贡献如下：（1）提出了一种有限样本算法，用于恢复标准LvLiNGAM观测等价类中的最稀疏DAG；（2）引入了一种直接对观测变量进行残差化而非递归更新高阶累积量的更新规则，从而消除了局部限制并提高了有限样本性能；（3）提出了一种新的父子准则，能够直接从有限样本中恢复最稀疏DAG；（4）在合成数据和真实数据上的实验证明了所提方法的有效性，特别是在违反ReLVLiNGAM局部限制时。 ## 2 预备知识 ### 2.1 标准LvLiNGAM 令X=\(X1,...,Xp\)⊤\\bm\{X\}=\(X_{1},\ldots,X_{p}\)^{\top}为观测变量，L=\(L1,...,Lq\)⊤\\bm\{L\}=\(L_{1},\ldots,L_{q}\)^{\top}为潜在混淆变量。将因果DAG表示为G=\(V,E\)\\mathcal\{G\}=\(\\bm\{V\},E\)，其中V=X∪L\\bm\{V\}=\(\\bm\{X\}\\cup\\bm\{L\}\)，E⊂V×V\\subset\\bm\{V\}\\times\\bm\{V\}为有向边集合。定义EO=\(X×X\)∩E^\{O\}=\(\\bm\{X\}\\times\\bm\{X\}\)\\cap E和EOL=\(L×X\)∩E^\{OL\}=\(\\bm\{L\}\\times\\bm\{X\}\)\\cap E。我们称GO=\(X,EO\)\mathcal\{G\}^\{O\}=\(\\bm\{X\},E^\{O\}\)为G\mathcal\{G\}的观测DAG，并令GOL=\(L,X,EOL\)\mathcal\{G\}^\{OL\}=\(\\bm\{L\},\\bm\{X\},E^\{OL\}\)为G\mathcal\{G\}的潜在到观测二分图。对于Xi∈XX_{i}\\in\\bm\{X\}，令Anc\(Xi\)\\mathrm\{Anc\}\(X_{i}\)、Des\(Xi\)\\mathrm\{Des\}\(X_{i}\)、Pa\(Xi\)\\mathrm\{Pa\}\(X_{i}\)和Ch\(Xi\)\\mathrm\{Ch\}\(X_{i}\)分别表示其祖先、后代、父节点和子节点集合。对于两个变量V,V′∈VV,V^{\\prime}\\in\\bm\{V\}\)，我们用V→V′\\to V^{\\prime}表示从V到V′的有向边，并令P\(V,V′\)\mathcal\{P\}\(V,V^{\\prime}\)为从V到V′的有向路径集合。对于两个观测变量Xi和Xj，我们定义它们的（可能是潜在的）混淆变量为Conf\(Xi,Xj\)=\{V:∃π∈P\(V,Xi\),∃π′∈P\(V,Xj\)s\.t\.\(π∩π′\)∖\{V\}=∅\}\。\\displaystyle\\mathrm\{Conf\}\(X_{i},X_{j}\)=\\Big\\\{V:\\exists\\,\\pi\\in\\mathcal\{P\}\(V,X_{i}\),\\exists\\,\\pi^{\\prime}\\in\\mathcal\{P\}\(V,X_{j}\)\\text\{ s\.t\. \}\\big\(\\pi\\cap\\pi^{\\prime}\\big\)\\setminus\\\{V\\\}=\\emptyset\\Big\\\}。 在本文中，我们采用标准LvLiNGAM（hoyer2008estimation），其中潜在变量相互独立，且每个潜在变量至少有两个观测子节点。由G\mathcal\{G\}定义的模型表示为 X=ΛL\+BX\+e，\\displaystyle\\bm\{X\}=\(\\bm\{\\Lambda\}\\bm\{L\}+\\bm\{B\}\\bm\{X\}+\\bm\{e\}\)，\(1\)其中e=\(e1,...,ep\)⊤\\bm\{e\}=\(e_{1},\ldots,e_{p}\)^{\top}\)，或者等价地， X\\displaystyle\\bm\{X\}=\[(I−B)−1Λ,(I−B)−1\][L⊤,e⊤]⊤。\\displaystyle=\\left\[\\left\(\\bm\{I\}-\\bm\{B\}\\right\)^{-1}\(\\bm\{\\Lambda\},\\;\\;\\left\(\\bm\{I\}-\\bm\{B\}\\right\)^{-1}\\right\]\\left\[\\bm\{L\}^{\top},\\bm\{e\}^{\top}\\right\]^{\top}。\(2\)u=\(L⊤,e⊤\)⊤\\bm\{u\}=\(\\bm\{L\}^{\top},\\bm\{e\}^{\top}\)^{\top}的分量相互独立，且各自服从具有非零高阶累积量的连续非高斯分布。Λ=\{λji\}\\bm\{\\Lambda\}=\{\\lambda_{ji}\}和B=\{bji\}\\bm\{B\}=\{b_{ji}\}分别收集Li→Xj∈E和Xi→Xj∈E的直接因果系数，M=\[(I−B)−1Λ,(I−B)−1\]\\bm\{M\}=\\left\[\\left\(\\bm\{I\}-\\bm\{B\}\\right\)^{-1}\(\\bm\{\\Lambda\},\\;\\;\\left\(\\bm\{I\}-\\bm\{B\}\\right\)^{-1}\\right\]是混合矩阵，其元素为总效应。令mjiOL^\{OL\}_{ji}和mjiO^\{O\}_{ji}分别表示Li和Xi对Xj的总效应。由于潜在尺度是任意的，不失一般性，我们对具有最高因果顺序的Xj∈Ch\(Li\)固定λji=1。相同的归一化也用于schkoda2024causal。 在本文中，我们假设(1 (https://arxiv.org/html/2607.05984#S2.E1))中的系数和u\\bm\{u\}的高阶累积量是一般性的。我们称此假设为一般性假设。特别地，本文的结果在系数和累积量集合的勒贝格测度为零的集除外成立。 我们记P\(V\)P\(\\bm\{V\}\)为V\\bm\{V\}的联合分布。假设在(2 (https://arxiv.org/html/2607.05984#S2.E2))中，存在Xi和Lj使得Des\(Lj\)=Des\(Xi\)∪\{Xi\}\\mathrm\{Des\}\(L_{j}\)=\\mathrm\{Des\}\(X_{i}\)\\cup\\\{X_{i}\}。salehkaleybar2020learning表明，在此条件下，交换混合矩阵M\\bm\{M\}中对应Lj和ei的列，会得到另一个有效的LvLiNGAM表示，同时保留P\(V\)P\(\\bm\{V\}\)和X\\bm\{X\}的祖先关系。根据他们的术语，我们将此列交换操作称为Lj和ei之间的交换。图1 (https://arxiv.org/html/2607.05984#S2.F1)展示了这样一个例子。DAG (a)是原始因果DAGG\\mathcal\{G\}，而DAG (b)是通过交换L3和e2得到的。尽管两个DAG诱导相同的P\(V\)P\(\\bm\{V\}\)和祖先关系，但交换在DAG (b)中引入了一条额外的有向边X1→X3。这个例子表明，多个因果DAG可能属于同一观测等价模式。然而，观测等价类内最稀疏的DAG特别令人感兴趣。在上面的例子中，DAG (a)比DAG (b)更稀疏。此外，在一般性假设下，DAG (a)的一般参数化对应于DAG (b)参数空间的一个测度零子集。因此，最稀疏DAG可以被视为观测等价类的标准代表。schkoda2024causal表明，在一般性假设下，标准LvLiNGAM的混合矩阵是可识别的，这进而识别了对应的观测等价类。因此，类内的最稀疏DAG也是可识别的。然而，ReLVLiNGAM专注于混合矩阵的一致估计，而不是直接从有限样本中恢复最稀疏DAG。此外，它要求每次迭代中，每个观测源的潜在父节点数量少于其每个观测兄弟节点的潜在父节点数量，我们称此条件为局部限制。图1 (https://arxiv.org/html/2607.05984#S2.F1)所示的模型违反了此限制，因此ReLVLiNGAM无法直接应用。 \\subfigure \[\]Refer to caption\\subfigure\[\]Refer to caption 图1：具有不同DAG的观测等价模型示例。 ### 2.2 累积量 在本节中，我们回顾schkoda2024causal建立的关于LvLiNGAM变量的高阶累积量与总效应之间关系的几个结果，这些结果在本文中也起着核心作用。 ###### 定义2.1（累积量（Brillinger））。令I=\[p\]\mathcal\{I\}=\[p\]为索引集合。对于pp维变量X=\(X1,...,Xp\)⊤\\bm\{X\}=\(X_{1},\dots,X_{p}\)^{\top}\)，kk阶累积量ci1,...,ik\(k\)c^{\(k\)}_{i_{1},\dots,i_{k}}定义为 ci1,...,ik\(k\)=∑Di∈\{D1,...,Dh\}\(−1\)h−1\(h−1\)\!∏di∈DiE\[∏j∈diXj\]，\\displaystyle c^{\(k\)}_{i_{1},\dots,i_{k}}=\\sum_{D_{i}\\in\\{D_{1},\dots,D_{h}\}}\(-1\)^{h-1}\(h-1\)!\\prod_{d_{i}\\in D_{i}}\\mathbb{E}\\left\[\\prod_{j\\in d_{i}}X_{j}\\right\]，其中\{i1,...,ik\}∈Ik\\{i_{1},\dots,i_{k}\\}\\in\\mathcal\{I\}^{k}且\{D1,...,Dh\}\\{D_{1},\dots,D_{h}\}是\{i1,i2,...,ik\}\\{i_{1},i_{2},\dots,i_{k}\}的所有划分的集合。当i1=i2=⋯=ik=i时，κ\(k\)\(Xi\)\\kappa^{\(k\)}\(X_{i}\)表示Xi的kk阶累积量。 从(2 (https://arxiv.org/html/2607.05984#S2.E2))出发，观测变量的kk阶累积量可以重写为 ci1,...,ik\(k\)=∑h=1qmi1hOL...mikhOLκ\(k\)\(Lh\)\+∑h=1pmi1hO...mikhOκ\(k\)\(eh\)。\\displaystyle c^{\(k\)}_{i_{1},\dots,i_{k}}=\\sum_{h=1}^{q}m^{OL}_{i_{1}h}\\dots m^{OL}_{i_{k}h}\\kappa^{\(k\)}\(L_{h}\)+\\sum_{h=1}^{p}m^{O}_{i_{1}h}\\dots m^{O}_{i_{k}h}\\kappa^{\(k\)}\(e_{h}\)。\(3\) 固定两个观测变量Xi和Xj，并将V∖\{Xi,Xj\}\\bm\{V\}\\setminus\\{X_{i},X_{j}\}中的所有变量视为潜在变量。应用hoyer2008estimation中的算法A，我们得到一个标准LvLiNGAM，其中潜在混淆变量Conf\(Xi,Xj\)=\{L1′,L2′,⋯,Ll′\}\\mathrm\{Conf\}\(X_{i},X_{j}\)=\\{L^{\\prime}_{1},L^{\\prime}_{2},\\cdots,L^{\\prime}_{\\ell}\\}相互独立。不失一般性，假设Xj∉Anc\(Xi\)X_{j}\\notin\\mathrm\{Anc\}\(X_{i}\)。那么，标准LvLiNGAM表示为 Xi\\displaystyle X_{i}=∑h=1lmihOL′Lh′\+vi,Xj=∑h=1lmjhOL′Lh′\+mjiOvi\+vj，\\displaystyle\\!=\\!\\sum_{h=1}^{\\ell}\{m^{OL^{\\prime}}_{ih}\}L^{\\prime}_{h}\\!+\\!v_{i},\\quad X_{j}\\!=\\!\\sum^{\\ell}_{h=1}\{m^{OL^{\\prime}}_{jh}\}L^{\\prime}_{h}\\!+\\!m^{O}_{ji}v_{i}\\!+\\!v_{j}\)，\(4\)其中vi和vj是扰动项，mihOL′和mjhOL′\{m^{OL^{\\prime}}_{ih}\}\\text\{ and \}\{m^{OL^{\\prime}}_{jh}\}分别是原模型中Xi和Xj上的潜在混淆变量Lh′到Xi和Xj的总效应。l\\ell是Xi和Xj之间混淆变量的数量。 schkoda2024causal在一般性假设下使用观测变量的高阶累积量估计标准LvLiNGAM。他们定义了矩阵Aj,i\(k1,k2\)A^{\(k_{1},k_{2}\)}_{\{j\},\{i\}}在(12 (https://arxiv.org/html/2607.05984#S2.E12))中，满足rank⁡\(Aj,i\(r\)\)\>r\+1\\operatorname{rank}\\left\(A_{j,i}^{\(r\)}\\right\)\\>r\+1。因此，l\\ell是满足rank⁡\(Aj,i\(r\)\)=r\+1\\operatorname{rank}\\left\(A_{j,i}^{\(r\)}\\right\)=r\+1的最小r值。