Learning Sparsest Linear Causal DAGs with Latent Confounders via Higher-Order Cumulants

# Learning Sparsest Linear Causal DAGs with Latent Confounders via Higher-Order Cumulants
Source: [https://arxiv.org/html/2607.05984](https://arxiv.org/html/2607.05984)
Kyoto UniversityKyotoJapan\\NameHisayuki Hara\\Emailhara\.hisayuki\.8k@kyoto\-u\.ac\.jp \\addrInstitute for Liberal Arts and SciencesKyoto UniversityKyotoJapan

###### Abstract

Recovering the exact directed acyclic graph \(DAG\) in linear non\-Gaussian acyclic models with latent confounders \(LvLiNGAM\) remains a challenging problem\. Although LvLiNGAM is identifiable only up to an observational equivalence class, each equivalence class is characterized by a unique sparsest DAG\. Recovering the sparsest DAG from finite samples, however, remains difficult\. Although existing methods are asymptotically consistent, they do not provide an explicit finite\-sample procedure for recovering the unique sparsest DAG, nor do they handle models with an arbitrary number of latent confounders\.

In this paper, we propose a finite\-sample method for recovering the sparsest DAG without imposing any restriction on the number of latent confounders\. Simulation studies and real\-data analyses demonstrate that the proposed method achieves superior finite\-sample performance compared with existing approaches\.

###### keywords:

Causal discovery; DAG; LiNGAM; Latent confounder; Cumulant;

## 1Introduction

Linear non\-Gaussian acyclic models \(LiNGAMs\) provide a powerful framework for causal discovery\(Shimizu2006;Shimizu2011\)\. In the absence of latent variables, LiNGAM enables complete identification of causal DAGs\. In many practical applications, however, latent confounders are unavoidable\.hoyer2008estimationintroduced LiNGAM with latent variables \(LvLiNGAM\) and demonstrated that any LvLiNGAM can be transformed into a canonical model in which all latent variables are mutually independent and causally precede the observed ones\. They also estimated the mixing matrix using overcomplete independent component analysis \(OICA; e\.g\.,eriksson2004identifiability\), assuming that the number of latent variables is known a priori\. However, because it relies on OICA, this approach is prone to converge to local optima\(shimizu14aBayesianestimation\)\.

To avoid relying on OICA, several methods have been proposed to estimate canonical LvLiNGAMs via residual independence tests, such as Pairwise LvLiNGAM\(Entner2011\), ParceLiNGAM\(tashiro2014parcelingam\), Repetitive Causal Discovery \(RCD\)\(Maeda2020;Maeda2022;maeda2022rcd\), and BANG\(wang2023BANG\)\. However, none of these methods can fully identify ancestral relationships or parent\-child relationships between observed variables that form the bow structures\(wang2023BANG\)\.

When bow structures are present,chen2024identificationuse cumulants of observed variables to identify their ancestral relationships in the bivariate setting with a single latent confounder\. Building on this,chen2025identificationextend the approach to multiple latent confounders in the observed bivariate case\.

schkoda2024causalproposed ReLVLiNGAM, a recursive cumulant\-based method that accommodates multiple observed variables and latent confounders\. Without relying on OICA or requiring prior knowledge of the number of latent variables, ReLVLiNGAM recovers the observational equivalence class of a canonical LvLiNGAM\.

Under the genericity assumption described below, the sparsest DAG within an observational equivalence class is uniquely determined\. The sparsest DAG is generic within its model class, whereas any denser DAG in the same observational equivalence class requires non\-generic parameter values\. This provides a natural justification for treating the sparsest DAG as the canonical representative of the observational equivalence class\. Although ReLVLiNGAM consistently estimates the mixing matrix, it does not provide a procedure for estimating the sparsest DAG from finite samples\.

Moreover, the original ReLVLiNGAM requires that, in each iteration, an observed variable with no observed parents \(an observed source hereafter\) has fewer latent parents than observed siblings; otherwise, the cumulant updates and mixing\-matrix estimation may fail \(Figure[1](https://arxiv.org/html/2607.05984#S2.F1); see Appendix[C](https://arxiv.org/html/2607.05984#A3)\)\. We refer to this issue as ReLVLiNGAM’s local restriction\.

In this paper, we develop a finite\-sample algorithm for recovering the sparsest DAG in the observational equivalence class of canonical LvLiNGAMs without the local restriction\. The proposed method builds on the top\-down framework of ReLVLiNGAM\. It successively infers parent\-child relationships among observed variables, starting from observed sources, to recover the sparsest DAG\. Within this framework, we introduce two key innovations\.

First, we introduce an update rule that directly residualizes the observed variables, rather than recursively updating higher\-order cumulants as in ReLVLiNGAM\. Updating the observed variables rather than their cumulants mitigates the recursive propagation of errors in higher\-order cumulant estimates, thereby improving finite\-sample performance\. Moreover, this update rule eliminates the need for the local restriction in ReLVLiNGAM, extending the applicability of the method\.

Second, we introduce a sequential procedure for identifying exact parent–child relationships between each observed source and its descendants from the estimated ancestral structure\. This procedure enables direct recovery of the sparsest DAG from finite samples\.

Our main contributions are as follows: \(1\) we propose a finite\-sample algorithm for recovering the sparsest DAG in the observational equivalence class of canonical LvLiNGAMs; \(2\) we introduce an update rule that directly residualizes the observed variables instead of recursively updating higher\-order cumulants, thereby removing the local restriction and improving finite\-sample performance; \(3\) we propose a new parent–child criterion that enables direct recovery of the sparsest DAG from finite samples; \(4\) experiments on synthetic and real data demonstrate the effectiveness of the proposed method, particularly when ReLVLiNGAM’s local restriction is violated\.

## 2Preliminaries

### 2\.1Canonical LvLiNGAM

Let𝑿=\(X1,…,Xp\)⊤\\bm\{X\}=\(X\_\{1\},\\ldots,X\_\{p\}\)^\{\\top\}be the observed variables and𝑳=\(L1,…,Lq\)⊤\\bm\{L\}=\(L\_\{1\},\\ldots,L\_\{q\}\)^\{\\top\}be the latent confounders\. Denote a causal DAG by𝒢=\(𝑽,E\)\\mathcal\{G\}=\(\\bm\{V\},E\), where𝑽=𝑿∪𝑳\\bm\{V\}=\\bm\{X\}\\cup\\bm\{L\}andE⊂𝑽×𝑽E\\subset\\bm\{V\}\\times\\bm\{V\}is the set of directed edges\. DefineEO=\(𝑿×𝑿\)∩EE^\{O\}=\(\\bm\{X\}\\times\\bm\{X\}\)\\cap EandEO​L=\(𝑳×𝑿\)∩EE^\{OL\}=\(\\bm\{L\}\\times\\bm\{X\}\)\\cap E\. We call𝒢O=\(𝑿,EO\)\\mathcal\{G\}^\{O\}=\(\\bm\{X\},E^\{O\}\)the observed DAG of𝒢\\mathcal\{G\}, and let𝒢O​L=\(𝑳,𝑿,EO​L\)\\mathcal\{G\}^\{OL\}=\(\\bm\{L\},\\bm\{X\},E^\{OL\}\)be the latent\-to\-observed bipartite graph of𝒢\\mathcal\{G\}\. ForXi∈𝑿X\_\{i\}\\in\\bm\{X\}, letAnc​\(Xi\)\\mathrm\{Anc\}\(X\_\{i\}\),Des​\(Xi\)\\mathrm\{Des\}\(X\_\{i\}\),Pa​\(Xi\)\\mathrm\{Pa\}\(X\_\{i\}\), andCh​\(Xi\)\\mathrm\{Ch\}\(X\_\{i\}\)denote the sets of ancestors, descendants, parents, and children ofXiX\_\{i\}, respectively\. For two variablesV,V′∈𝑽V,V^\{\\prime\}\\in\\bm\{V\}, we denote a directed edge fromVVtoV′V^\{\\prime\}byV→V′V\\to V^\{\\prime\}, and let𝒫​\(V,V′\)\\mathcal\{P\}\(V,V^\{\\prime\}\)be the set of directed paths fromVVtoV′V^\{\\prime\}\. For two observed variablesXiX\_\{i\}andXjX\_\{j\}, we define their \(possibly latent\) confounders as

Conf​\(Xi,Xj\)=\{V:∃π∈𝒫​\(V,Xi\),∃π′∈𝒫​\(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\\\}\.
In this paper, we employ the canonical LvLiNGAM\(hoyer2008estimation\), where latent variables are mutually independent, and each has at least two observed children\. The model defined by𝒢\\mathcal\{G\}is expressed as

𝑿=𝚲​𝑳\+𝑩​𝑿\+𝒆,\\displaystyle\\bm\{X\}=\\bm\{\\Lambda\}\\bm\{L\}\+\\bm\{B\}\\bm\{X\}\+\\bm\{e\},\(1\)where𝒆=\(e1,…,ep\)⊤\\bm\{e\}=\(e\_\{1\},\\ldots,e\_\{p\}\)^\{\\top\}, or, equivalently,

𝑿\\displaystyle\\bm\{X\}=\[\(𝑰−𝑩\)−1​𝚲,\(𝑰−𝑩\)−1\]​\[𝑳⊤,𝒆⊤\]⊤\.\\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\)The components of𝒖=\(𝑳⊤,𝒆⊤\)⊤\\bm\{u\}=\\left\(\\bm\{L\}^\{\\top\},\\bm\{e\}^\{\\top\}\\right\)^\{\\top\}are mutually independent, and each follows a continuous non\-Gaussian distribution with nonzero higher\-order cumulants\.𝚲=\{λj​i\}\\bm\{\\Lambda\}=\\\{\\lambda\_\{ji\}\\\}and𝑩=\{bj​i\}\\bm\{B\}=\\\{b\_\{ji\}\\\}collect the direct causal coefficients forLi→Xj∈EL\_\{i\}\\to X\_\{j\}\\in EandXi→Xj∈EX\_\{i\}\\to X\_\{j\}\\in E, respectively, and𝑴=\[\(𝑰−𝑩\)−1​𝚲,\(𝑰−𝑩\)−1\]\\bm\{M\}=\\left\[\\left\(\\bm\{I\}\-\\bm\{B\}\\right\)^\{\-1\}\\bm\{\\Lambda\},\\;\\;\\left\(\\bm\{I\}\-\\bm\{B\}\\right\)^\{\-1\}\\right\]is the mixing matrix whose entries are total effects\. Letmj​iO​Lm^\{OL\}\_\{ji\}andmj​iOm^\{O\}\_\{ji\}denote the total effects ofLiL\_\{i\}andXiX\_\{i\}onXjX\_\{j\}, respectively\. Since latent scales are arbitrary, without loss of generality, we fixλj​i=1\\lambda\_\{ji\}=1forXj∈Ch​\(Li\)X\_\{j\}\\in\\mathrm\{Ch\}\(L\_\{i\}\)with the highest causal order\. The same normalization is also used inschkoda2024causal\.

Throughout this paper, we assume that the coefficients in \([1](https://arxiv.org/html/2607.05984#S2.E1)\) and the higher\-order cumulants of𝒖\\bm\{u\}are generic\. We call this assumption the genericity assumption\. In particular, the results in this paper hold except for a set of coefficients and cumulants of Lebesgue measure zero\.

We denote byP​\(𝑽\)P\(\\bm\{V\}\)the joint distribution of𝑽\\bm\{V\}\. Suppose that, in \([2](https://arxiv.org/html/2607.05984#S2.E2)\), there existXiX\_\{i\}andLjL\_\{j\}such thatDes​\(Lj\)=Des​\(Xi\)∪\{Xi\}\\mathrm\{Des\}\(L\_\{j\}\)=\\mathrm\{Des\}\(X\_\{i\}\)\\cup\\\{X\_\{i\}\\\}\.salehkaleybar2020learningshowed that, under this condition, exchanging the columns of the mixing matrix𝑴\\bm\{M\}corresponding toLjL\_\{j\}andeie\_\{i\}yields another valid LvLiNGAM representation while preserving bothP​\(𝑽\)P\(\\bm\{V\}\)and the ancestral relationships of𝑿\\bm\{X\}\. Depending on the graph structure, the resulting representation may correspond to a different causal DAG\. Following their terminology, we refer to this column\-exchange operation as a swap betweenLjL\_\{j\}andeie\_\{i\}\.

Figure[1](https://arxiv.org/html/2607.05984#S2.F1)illustrates such an example\. DAG \(a\) is the original causal DAG𝒢\\mathcal\{G\}, whereas DAG \(b\) is obtained by a swap betweenL3L\_\{3\}ande2e\_\{2\}\. Although the two DAGs induce the sameP​\(𝑽\)P\(\\bm\{V\}\)and ancestral relationships, the swap introduces an additional directed edgeX1→X3X\_\{1\}\\to X\_\{3\}in DAG \(b\)\. This example shows that multiple causal DAGs may belong to the same observational equivalence pattern\.

Nevertheless, the sparsest DAG within an observational equivalence class is of particular interest\. In the example above, DAG \(a\) is sparser than DAG \(b\)\. Moreover, under the genericity assumption, a generic parameterization of DAG \(a\) corresponds to a measure\-zero subset of the parameter space of DAG \(b\)\. Thus, the sparsest DAG can be regarded as a canonical representative of the observational equivalence class\.

schkoda2024causalshowed that the mixing matrix of a canonical LvLiNGAM is identifiable under the genericity assumption, which in turn identifies the corresponding observational equivalence class\. Hence, the sparsest DAG within the class is also identifiable\. However, ReLVLiNGAM focuses on consistent estimation of the mixing matrix, rather than direct recovery of the sparsest DAG from finite samples\. Moreover, it requires that, at each iteration, every observed source has fewer latent parents than each of its observed siblings, a condition that we refer to as the local restriction\. The model shown in Figure[1](https://arxiv.org/html/2607.05984#S2.F1)violates this restriction, and therefore ReLVLiNGAM cannot be applied directly\.

\\subfigure

\[\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x1.png)\\subfigure\[\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x2.png)

Figure 1:An example of observationally equivalent models with different DAGs\.
### 2\.2Cumulants

In this section, we review several results established byschkoda2024causalconcerning the relationship between higher\-order cumulants of LvLiNGAM variables and total effects, which also play a central role in this paper\.

###### Definition 2\.1\(Cumulants\(Brillinger\)\)\.

Letℐ=\[p\]\\mathcal\{I\}\\\!=\\\!\[p\]be the set of indices\. For app\-dimensional variable𝐗=\(X1,…,Xp\)⊤\\bm\{X\}\\\!=\\\!\(X\_\{1\},\\\!\\dots,\\\!X\_\{p\}\)^\{\\top\}, thekk\-th order cumulantci1,…,ik\(k\)c^\{\(k\)\}\_\{i\_\{1\},\\dots,i\_\{k\}\}is defined by

ci1,…,ik\(k\)=∑Di∈\{D1,…,Dh\}\(−1\)h−1​\(h−1\)\!​∏di∈Di𝔼​\[∏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\],where\{i1,…,ik\}∈ℐk\\\{i\_\{1\},\\dots,i\_\{k\}\\\}\\in\\mathcal\{I\}^\{k\}and\{D1,…,Dh\}\\\{D\_\{1\},\\dots,D\_\{h\}\\\}is the set of all partitions of\{i1,i2,…,ik\}\\\{i\_\{1\},i\_\{2\},\\dots,i\_\{k\}\\\}\.

Wheni1=i2=⋯=ik=ii\_\{1\}=i\_\{2\}=\\cdots=i\_\{k\}=i,κ\(k\)​\(Xi\)\\kappa^\{\(k\)\}\(X\_\{i\}\)denotes thekk\-th order cumulant ofXiX\_\{i\}\. From \([2](https://arxiv.org/html/2607.05984#S2.E2)\), thekk\-th order cumulants of observed variables can be rewritten as

ci1,…,ik\(k\)=∑h=1qmi1​hO​L​…​mik​hO​L​κ\(k\)​\(Lh\)\+∑h=1pmi1​hO​…​mik​hO​κ\(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\)Fix two observed variablesXiX\_\{i\}andXjX\_\{j\}, and treat all variables in𝑽∖\{Xi,Xj\}\\bm\{V\}\\setminus\\\{X\_\{i\},X\_\{j\}\\\}as latent\. Applying Algorithm A inhoyer2008estimation, we obtain a canonical LvLiNGAM in which the latent confoundersConf​\(Xi,Xj\)=\{L1′,L2′,⋯,Lℓ′\}\\mathrm\{Conf\}\(X\_\{i\},X\_\{j\}\)=\\\{L^\{\\prime\}\_\{1\},L^\{\\prime\}\_\{2\},\\cdots,L^\{\\prime\}\_\{\\ell\}\\\}are mutually independent\. Without loss of generality, assumeXj∉Anc​\(Xi\)X\_\{j\}\\notin\\mathrm\{Anc\}\(X\_\{i\}\)\. Then, the canonical LvLiNGAM is expressed as

Xi\\displaystyle X\_\{i\}=∑h=1ℓmi​hO​L′​Lh′\+vi,Xj=∑h=1ℓmj​hO​L′​Lh′\+mj​iO​vi\+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\)whereviv\_\{i\}andvjv\_\{j\}are disturbances, andmi​hO​L′​and​mj​hO​L′\{m^\{OL^\{\\prime\}\}\_\{ih\}\}\\text\{ and \}\{m^\{OL^\{\\prime\}\}\_\{jh\}\}are total effects fromLh′L^\{\\prime\}\_\{h\}toXiX\_\{i\}andXjX\_\{j\}, respectively, in the canonical model overXiX\_\{i\}andXjX\_\{j\}\.ℓ\\ellis the number of confounders betweenXiX\_\{i\}andXjX\_\{j\}in the original model\.

schkoda2024causalestimate the canonical LvLiNGAM under the genericity assumption using higher\-order cumulants of observed variables\. They define the matrixAj,i\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{\{j\},\{i\}\}in \([12](https://arxiv.org/html/2607.05984#S2.E12)\) withk1<k2k\_\{1\}<k\_\{2\}\.

Aj,i\(k1,k2\)=\[ci,i,…,i\(k1\)ci,i,…,j\(k1\)…ci,j,…,j\(k1\)ci,i,i,…,i\(k1\+1\)ci,i,i,…,j\(k1\+1\)…ci,i,j,…,j\(k1\+1\)cj,i,i,…,i\(k1\+1\)cj,i,i,…,j\(k1\+1\)…cj,i,j,…,j\(k1\+1\)⋮⋮⋱⋮ci,…,i,i,i,…,i,i\(k2\)ci,…,i,i,i,…,i,j\(k2\)…ci,…,i,i,j,…,j,j\(k2\)⋮⋮⋱⋮cj,…,j,i,i,…,i,i\(k2\)cj,…,j,i,i,…,i,j\(k2\)…cj,…,j,i,j,…,j,j\(k2\)\]\.\\displaystyle A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{\{j\},\{i\}\}=\\left\[\\begin\{array\}\[\]\{cccc\}c^\{\(k\_\{1\}\)\}\_\{i,i,\\dots,i\}&c^\{\(k\_\{1\}\)\}\_\{i,i,\\dots,j\}&\\dots&c^\{\(k\_\{1\}\)\}\_\{i,j,\\dots,j\}\\\\ c^\{\(k\_\{1\}\+1\)\}\_\{i,i,i,\\dots,i\}&c^\{\(k\_\{1\}\+1\)\}\_\{i,i,i,\\dots,j\}&\\dots&c^\{\(k\_\{1\}\+1\)\}\_\{i,i,j,\\dots,j\}\\\\ c^\{\(k\_\{1\}\+1\)\}\_\{j,i,i,\\dots,i\}&c^\{\(k\_\{1\}\+1\)\}\_\{j,i,i,\\dots,j\}&\\dots&c^\{\(k\_\{1\}\+1\)\}\_\{j,i,j,\\dots,j\}\\\\ \\vdots&\\vdots&\\ddots&\\vdots\\\\ c^\{\(k\_\{2\}\)\}\_\{i,\\dots,i,i,i,\\dots,i,i\}&c^\{\(k\_\{2\}\)\}\_\{i,\\dots,i,i,i,\\dots,i,j\}&\\dots&c^\{\(k\_\{2\}\)\}\_\{i,\\dots,i,i,j,\\dots,j,j\}\\\\ \\vdots&\\vdots&\\ddots&\\vdots\\\\ c^\{\(k\_\{2\}\)\}\_\{j,\\dots,j,i,i,\\dots,i,i\}&c^\{\(k\_\{2\}\)\}\_\{j,\\dots,j,i,i,\\dots,i,j\}&\\dots&c^\{\(k\_\{2\}\)\}\_\{j,\\dots,j,i,j,\\dots,j,j\}\\end\{array\}\\right\]\.\(12\)DefineAi,j\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{i,j\}analogously by swappingiiandjj\. Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)allows us to identifyℓ\\ellin the model \([4](https://arxiv.org/html/2607.05984#S2.E4)\) and the causal order betweenXiX\_\{i\}andXjX\_\{j\}\.

###### Proposition 2\.2\(schkoda2024causal\)\.

Assume thatXiX\_\{i\}andXjX\_\{j\}are two observed variables whereXj∉Anc​\(Xi\)X\_\{j\}\\notin\\mathrm\{Anc\}\(X\_\{i\}\)\. Letd:=min⁡\(∑i=1k2−k1\+1i,k1\)d:=\\min\(\\sum^\{k\_\{2\}\-k\_\{1\}\+1\}\_\{i=1\}i,k\_\{1\}\)\. Then,

1\.Aj,i\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{\{j\},\{i\}\}generically has rankmin⁡\(ℓ\+1,d\)\\min\(\\ell\+1,d\)\.

2\. Ifmj​iO≠0m^\{O\}\_\{ji\}\\\!\\neq\\\!0,Ai,j\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{\{i\},\{j\}\}generically has rankmin⁡\(ℓ\+2,d\)\\min\(\\ell\+2,\\\!d\)\.

3\. Ifmj​iO=0m^\{O\}\_\{ji\}\\\!=\\\!0,Ai,j\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{\{i\},\{j\}\}generically has rankmin⁡\(ℓ\+1,d\)\\min\(\\ell\+1,\\\!d\)\.

According toschkoda2024causal, the smallest possible choice of\(k1,k2\)\(k\_\{1\},k\_\{2\}\)is\(ℓ\+2,\(ℓ\+2\)\+⌈\(−3\+8​ℓ\+17\)/2⌉\)\(\\ell\+2,\(\\ell\+2\)\+\\lceil\(\-3\+\\sqrt\{8\\ell\+17\}\)/2\\rceil\)\. DefineAj,i\(ℓ\)\{A\}^\{\(\\ell\)\}\_\{j,i\}asAj,i\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{j,i\}, where\(k1,k2\)\(k\_\{1\},k\_\{2\}\)is this choice\.

Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)provides a practical criterion for determiningℓ\\ell\. If the true number of confounders isℓ\\ell, then Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)implies that the matrixAj,i\(ℓ\)A\_\{j,i\}^\{\(\\ell\)\}generically has rankℓ\+1\\ell\+1\. Consequently, for anyr<ℓr<\\ell,rank⁡\(Aj,i\(r\)\)\>r\+1\\operatorname\{rank\}\\left\(A\_\{j,i\}^\{\(r\)\}\\right\)\>r\+1\. Therefore,ℓ\\ellis the smallest value ofrrsatisfyingrank⁡\(Aj,i\(r\)\)=r\+1\\operatorname\{rank\}\\left\(A\_\{j,i\}^\{\(r\)\}\\right\)=r\+1\. Moreover, Items 2 and 3 of Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)determine whether an ancestral relationship exists betweenXiX\_\{i\}andXjX\_\{j\}, and, if so, identify its direction\.

LetA~j,i\(ℓ\)\\tilde\{A\}^\{\(\\ell\)\}\_\{\{j\},\{i\}\}be a matrix obtained by adding\(1,m,…,mℓ\+1\)\(1,\\\!m,\\\!\\dots,\\\!m^\{\\ell\+1\}\)as the first row ofAj,i\(ℓ\)A^\{\(\\ell\)\}\_\{\{j\},\{i\}\}\.

###### Proposition 2\.3\(schkoda2024causal\)\.

Consider the determinant of an\(ℓ\+2\)×\(ℓ\+2\)\(\\ell\+2\)\\times\(\\ell\+2\)minor ofA~j,i\(ℓ\)\\tilde\{A\}^\{\(\\ell\)\}\_\{\{j\},\{i\}\}that contains the first row and treat it as a polynomial inmm\. Then, the roots of this polynomial aremj​iO,mj​1O​L′,⋯,mj​ℓO​L′m^\{O\}\_\{ji\},m^\{OL^\{\\prime\}\}\_\{j1\},\\cdots,m^\{OL^\{\\prime\}\}\_\{j\\ell\}\.

Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)identifies total effectsmj​iO,mj​1O​L′,⋯,mj​ℓO​L′m^\{O\}\_\{ji\},m^\{OL^\{\\prime\}\}\_\{j1\},\\cdots,m^\{OL^\{\\prime\}\}\_\{j\\ell\}up to permutation\.

###### Proposition 2\.4\(schkoda2024causal\)\.

Define𝐦j​iO:=\(1,mj​iO,\(mj​iO\)2,…,\(mj​iO\)k−1\)⊤\\bm\{m\}\_\{ji\}^\{O\}:=\\left\(1,m\_\{ji\}^\{O\},\(m\_\{ji\}^\{O\}\)^\{2\},\\ldots,\(m\_\{ji\}^\{O\}\)^\{k\-1\}\\right\)^\{\\top\}\. Similarly, define𝐦j​1O​L′,…,𝐦j​ℓO​L′\\bm\{m\}\_\{j1\}^\{OL^\{\\prime\}\},\\ldots,\\bm\{m\}\_\{j\\ell\}^\{OL^\{\\prime\}\}\. Then, the system of equations

\[𝒎j​iO,𝒎j​1O​L,…,𝒎j​ℓO​L\]​\[κ\(k\)​\(vi\),κ\(k\)​\(L1′\),…,κ\(k\)​\(Lℓ′\)\]⊤\\displaystyle\\left\[\\bm\{m\}\_\{ji\}^\{O\},\\bm\{m\}\_\{j1\}^\{OL\},\\dots,\\bm\{m\}\_\{j\\ell\}^\{OL\}\\right\]\\left\[\\kappa^\{\(k\)\}\(\{v\}\_\{i\}\),\\kappa^\{\(k\)\}\(L^\{\\prime\}\_\{1\}\),\\dots,\\kappa^\{\(k\)\}\(L^\{\\prime\}\_\{\\ell\}\)\\right\]^\{\\top\}=\[ci,i,…,i\(k\),ci,i,…,j\(k\),…,ci,j,…,j\(k\)\]⊤,\\displaystyle\\qquad=\\\!\\left\[c^\{\(k\)\}\_\{i,i,\\dots,i\},c^\{\(k\)\}\_\{i,i,\\dots,j\},\\dots,c^\{\(k\)\}\_\{i,j,\\dots,j\}\\right\]^\{\\top\}\\\!,\(13\)is generically uniquely solvable ifk≥ℓ\+1k\\geq\\ell\+1\.

When \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) is uniquely solvable, each root returned by Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)is associated with a uniquekk\-th order cumulant of a disturbance or latent variable,

\(mj​iO,κ\(k\)​\(vi\)\),\(mj​iO​L,κ\(k\)​\(L1′\)\),…,\(mj​ℓO​L,κ\(k\)​\(Lℓ′\)\)\.\\left\(m\_\{ji\}^\{O\},\\kappa^\{\(k\)\}\(v\_\{i\}\)\\right\),\\left\(m\_\{ji\}^\{OL\},\\kappa^\{\(k\)\}\(L\_\{1\}^\{\\prime\}\)\\right\),\\ldots,\\left\(m\_\{j\\ell\}^\{OL\},\\kappa^\{\(k\)\}\(L\_\{\\ell\}^\{\\prime\}\)\\right\)\.Although the roots and cumulants themselves are identified only up to permutation, the correspondence between a root and its associated cumulant is uniquely determined\. Hereafter, without any additional explanation, we focus on orderskkfor which \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) is solvable\.

## 3Proposed Method

We propose a method to recover the sparsest DAG under the genericity assumption\. In the following, let𝒢\\mathcal\{G\}denote the sparsest DAG in the observational equivalence class\. The proposed method first applies Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)to infer all ancestral relationships among the observed variables and identify those with no observed ancestors as observed sources\. The same proposition also estimates the number of latent confounders for each pair of observed variables\. We then use the total effects obtained from Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3), together with their correspondences to thekk\-th order cumulants of the disturbances and latent confounders established by Proposition[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4), to determine the parent–child relationships between each observed source and its descendants\. We then identify the sparsest latent\-to\-observed structure by selecting, among the candidate total effects obtained above, those that minimize the numbers of latent confounders between pairs of observed variables\. The identified sources are then residualized from the remaining observed variables, after which the above estimation procedure is repeated on the updated variables\. Recursively repeating this process recovers the sparsest DAG\. Proofs of all theorems in this section are given in Appendix[B](https://arxiv.org/html/2607.05984#A2)\.

Based on the pairwise ancestral relationships estimated using Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2), variables without observed ancestors are treated as observed sources\. LetXsX\_\{s\}be one such source\. LetCh~​\(Xs\)⊂Ch​\(Xs\)\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)\\subset\\mathrm\{Ch\}\(X\_\{s\}\)be the set of observed variables whose observed ancestors consist only of observed sources and includeXsX\_\{s\}\. To identify the remaining children ofXsX\_\{s\}, we recursively examine its descendants while maintaining two sets: the closed set𝑿closed\\bm\{X\}\_\{\\mathrm\{closed\}\}, containing variables already identified as children ofXsX\_\{s\}, and the open set𝑿open\\bm\{X\}\_\{\\mathrm\{open\}\}, containing descendants whose parent\-child relationships withXsX\_\{s\}have yet to be determined\. Initially,

𝑿closed=Ch~​\(Xs\),\\displaystyle\\bm\{X\}\_\{\\mathrm\{\\mathrm\{closed\}\}\}=\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\),𝑿open=\{Xj∈Des​\(Xs\)∖Xclosed:∀Xi∈𝑿closed,Anc​\(Xj\)∩Des​\(Xi\)⊆𝑿closed\}\.\\displaystyle\{\\bm\{X\}\}\_\{\\mathrm\{\\mathrm\{open\}\}\}=\\Big\\\{X\_\{j\}\\in\\mathrm\{Des\}\(X\_\{s\}\)\\setminus X\_\{\\mathrm\{\\mathrm\{closed\}\}\}:\{\\forall X\_\{i\}\\in\\bm\{X\}\_\{\{\\mathrm\{\\mathrm\{closed\}\}\}\},\\penalty 10000\\ \\mathrm\{Anc\}\(X\_\{j\}\)\\cap\\mathrm\{Des\}\(X\_\{i\}\)\\subseteq\\bm\{X\}\_\{\\mathrm\{\\mathrm\{closed\}\}\}\}\\Big\\\}\.\(14\)By definition of𝑿open\\bm\{X\}\_\{\\mathrm\{\\mathrm\{open\}\}\}, no observed variable in𝑿∖𝑿closed\\bm\{X\}\\setminus\\bm\{X\}\_\{\\mathrm\{\\mathrm\{closed\}\}\}lies on a directed path between𝑿closed\\bm\{X\}\_\{\\mathrm\{\\mathrm\{closed\}\}\}and𝑿open\\bm\{X\}\_\{\\mathrm\{\\mathrm\{open\}\}\}\.

###### Lemma 3\.1\.

For anyXj∈𝐗openX\_\{j\}\\in\\bm\{X\}\_\{\\mathrm\{open\}\},bj​sb\_\{js\}can be written as

bj​s=mj​sO−∑i:Xi∈Ch~​\(Xs\)bi​s​mj​iO\.\\displaystyle b\_\{js\}=m^\{O\}\_\{js\}\-\\sum\_\{i:X\_\{i\}\\in\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)\}b\_\{is\}\\,m^\{O\}\_\{ji\}\.\(15\)

According to Lemma[3\.1](https://arxiv.org/html/2607.05984#S3.Thmtheorem1), the parent\-child relationship betweenXsX\_\{s\}andXjX\_\{j\}can be determined by testing whetherbj​s=0b\_\{js\}=0\. However, Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)identifies the total effects appearing in \([15](https://arxiv.org/html/2607.05984#S3.E15)\) only up to permutation, sobj​sb\_\{js\}cannot be computed directly\.

Denote byℐ\\mathcal\{I\}and𝒥\\mathcal\{J\}the index sets of𝑿closed\\bm\{X\}\_\{\\mathrm\{closed\}\}and𝑿open\\bm\{X\}\_\{\\mathrm\{open\}\}, respectively\. Forj∈𝒥j\\in\\mathcal\{J\}andi∈ℐi\\in\\mathcal\{I\}, let𝒎j​sO\\bm\{m\}^\{O\}\_\{js\}and𝒃i​s\\bm\{b\}\_\{is\}denote the sets of candidate values ofmj​sOm^\{O\}\_\{js\}andbi​sb\_\{is\}, respectively, returned by Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)\. Initially,bi​s=mi​sOb\_\{is\}=m^\{O\}\_\{is\}\.

As discussed in the previous section, Proposition[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)associates each candidate total effect with the correspondingkk\-th order cumulant of a disturbance or latent confounder\. For anyi∈ℐi\\in\\mathcal\{I\}andj∈𝒥j\\in\\mathcal\{J\}, the systems \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) for the pairs\(Xs,Xi\)\(X\_\{s\},X\_\{i\}\)and\(Xs,Xj\)\(X\_\{s\},X\_\{j\}\)always share thekk\-th order cumulant ofese\_\{s\}\. Under the genericity assumption, distinct disturbances and latent confounders have distinctkk\-th order cumulants\. Therefore, we can chooseαj​s∈𝒎j​sO\\alpha\_\{js\}\\in\\bm\{m\}^\{O\}\_\{js\}andβi​s∈𝒃i​s\\beta\_\{is\}\\in\\bm\{b\}\_\{is\}so that they correspond to the samekk\-th order cumulant\. Such a choice is not necessarily unique, since the two systems may also share thekk\-th order cumulants of latent confounders common to the pairs\(Xs,Xi\)\(X\_\{s\},X\_\{i\}\)and\(Xs,Xj\)\(X\_\{s\},X\_\{j\}\)\.

Define the residualized variables

X~i=Xi−βi​s​Xs,X~j=Xj−αj​s​Xs\.\\tilde\{X\}\_\{i\}=X\_\{i\}\-\\beta\_\{is\}X\_\{s\},\\qquad\\tilde\{X\}\_\{j\}=X\_\{j\}\-\\alpha\_\{js\}X\_\{s\}\.By Lemma[A\.1](https://arxiv.org/html/2607.05984#A1.Thmtheorem1)in Appendix[A](https://arxiv.org/html/2607.05984#A1),X~i\\tilde\{X\}\_\{i\}andX~j\\tilde\{X\}\_\{j\}can be regarded as observed variables in an induced canonical LvLiNGAM\. Applying Propositions[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)and[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)to this induced model yields a candidate set for the total effect fromX~i\\tilde\{X\}\_\{i\}toX~j\\tilde\{X\}\_\{j\}\. We denote this set by𝒎j​i\.sO\\bm\{m\}^\{O\}\_\{ji\.s\}\.

Theorem[3\.2](https://arxiv.org/html/2607.05984#S3.Thmtheorem2)provides a criterion for determining whetherXs∈Pa​\(Xj\)X\_\{s\}\\in\\mathrm\{Pa\}\(X\_\{j\}\)for eachj∈𝒥j\\in\\mathcal\{J\}\.

###### Theorem 3\.2\.

Under the genericity assumption,Xs∉Pa​\(Xj\)X\_\{s\}\\notin\\mathrm\{Pa\}\(X\_\{j\}\)if and only if there exist

αj​s∈𝒎j​sO,𝜷ℐ,s∈∏i∈ℐ𝒃i​s,𝜸j,ℐ\.s∈∏i∈ℐ𝒎j​i\.sO\\displaystyle\\alpha\_\{js\}\\in\\bm\{m\}^\{O\}\_\{js\},\\quad\\bm\{\\beta\}\_\{\\mathcal\{I\},s\}\\in\\prod\_\{i\\in\\mathcal\{I\}\}\\bm\{b\}\_\{is\},\\quad\\bm\{\\gamma\}\_\{j,\\mathcal\{I\}\.s\}\\in\\prod\_\{i\\in\\mathcal\{I\}\}\\bm\{m\}^\{O\}\_\{ji\.s\}s\.t\.αj​s−𝜷ℐ,s⊤​𝜸j,ℐ\.s=0,\\displaystyle\\qquad\\text\{s\.t\.\}\\qquad\\alpha\_\{js\}\-\\bm\{\\beta\}\_\{\\mathcal\{I\},s\}^\{\\top\}\\bm\{\\gamma\}\_\{j,\\mathcal\{I\}\.s\}=0,\(16\)whereαj​s\\alpha\_\{js\}andβi​s\\beta\_\{is\},i∈ℐi\\in\\mathcal\{I\}are chosen so that their associatedkk\-th order cumulants are equal\.

Once the parent\-child relationship betweenXsX\_\{s\}andXjX\_\{j\}is determined using Theorem[3\.2](https://arxiv.org/html/2607.05984#S3.Thmtheorem2), we update𝑿closed\\bm\{X\}\_\{\\mathrm\{closed\}\},𝑿open\\bm\{X\}\_\{\\mathrm\{open\}\}, andCh~​\(Xs\)\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)as follows:

𝑿closed←𝑿closed∪\{Xj\},𝑿open←𝑿open∖\{Xj\},Ch~​\(Xs\)←\{Ch~​\(Xs\),Xs∉Pa​\(Xj\),Ch~​\(Xs\)∪\{Xj\},Xs∈Pa​\(Xj\)\.\\displaystyle\\bm\{X\}\_\{\\mathrm\{closed\}\}\\leftarrow\\bm\{X\}\_\{\\mathrm\{closed\}\}\\cup\\\{X\_\{j\}\\\},\\ \\bm\{X\}\_\{\\mathrm\{open\}\}\\leftarrow\\bm\{X\}\_\{\\mathrm\{open\}\}\\setminus\\\{X\_\{j\}\\\},\\ \\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)\\leftarrow\\begin\{cases\}\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\),&X\_\{s\}\\not\\in\\mathrm\{Pa\}\(X\_\{j\}\),\\\\ \\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)\\cup\\\{X\_\{j\}\\\},&X\_\{s\}\\in\\mathrm\{Pa\}\(X\_\{j\}\)\.\\end\{cases\}When𝑿open\\bm\{X\}\_\{\\mathrm\{open\}\}becomes empty, it is reinitialized using \([3](https://arxiv.org/html/2607.05984#S3.Ex5)\)\. As shown in the proof of Lemma[3\.1](https://arxiv.org/html/2607.05984#S3.Thmtheorem1)in Appendix[B](https://arxiv.org/html/2607.05984#A2), Lemma[3\.1](https://arxiv.org/html/2607.05984#S3.Thmtheorem1)still holds after updatingCh~​\(Xs\)\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)\. The coefficients of the directed edges fromXsX\_\{s\}to the variables inCh~​\(Xs\)\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)can also be computed by Lemma[3\.1](https://arxiv.org/html/2607.05984#S3.Thmtheorem1)\. Therefore, the parent\-child relationships ofXsX\_\{s\}and variables in𝑿open\\bm\{X\}\_\{\\mathrm\{open\}\}can still be identified by Theorem[3\.2](https://arxiv.org/html/2607.05984#S3.Thmtheorem2)\.

After determining all parent\-child relationships betweenXsX\_\{s\}and its descendants, we residualize each descendantXj∈Des​\(Xs\)X\_\{j\}\\in\\mathrm\{Des\}\(X\_\{s\}\)with respect toXsX\_\{s\}before proceeding to the next iteration:

X~j=Xj−αj​s​Xs\.\\displaystyle\\tilde\{X\}\_\{j\}=X\_\{j\}\-\\alpha\_\{js\}X\_\{s\}\.\(17\)For eachj∈𝒥j\\in\\mathcal\{J\}, chooseαj​s\\alpha\_\{js\}from triples\(αj​s,𝜷ℐ,s,𝜸j,ℐ,s\)\(\\alpha\_\{js\},\\bm\{\\beta\}\_\{\\mathcal\{I\},s\},\\bm\{\\gamma\}\_\{j,\\mathcal\{I\},s\}\)satisfying \([3\.2](https://arxiv.org/html/2607.05984#S3.Ex7)\) whenever available; and otherwise, choose it randomly from𝒎j​sO\\bm\{m\}^\{O\}\_\{js\}, so that all selected values ofαj​s\\alpha\_\{js\}correspond to the samekk\-th order cumulant\.

Letℛ=\[p\]∖\{s\}\\mathcal\{R\}=\[p\]\\setminus\\\{s\\\}\. For eachj∈ℛj\\in\\mathcal\{R\}, letαj​s\\alpha\_\{js\}be defined as above whenXj∈Des​\(Xs\)X\_\{j\}\\in\\mathrm\{Des\}\(X\_\{s\}\), and setαj​s=0\\alpha\_\{js\}=0whenXj∉Des​\(Xs\)X\_\{j\}\\notin\\mathrm\{Des\}\(X\_\{s\}\)\. Let

𝑿ℛ=\(Xj\)j∈ℛ⊤,𝜶ℛ,s=\(αj​s\)j∈ℛ⊤\.\\displaystyle\\bm\{X\}\_\{\\mathcal\{R\}\}=\(X\_\{j\}\)\_\{j\\in\\mathcal\{R\}\}^\{\\top\},\\qquad\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}=\(\\alpha\_\{js\}\)^\{\\top\}\_\{j\\in\\mathcal\{R\}\}\.Then, the residualized variables𝑿~ℛ\\tilde\{\\bm\{X\}\}\_\{\\mathcal\{R\}\}after removingXsX\_\{s\}are given by

𝑿~ℛ=\(X~j\)j∈ℛ⊤:=𝑿ℛ−𝜶ℛ,s​Xs\.\\tilde\{\\bm\{X\}\}\_\{\\mathcal\{R\}\}=\(\\tilde\{X\}\_\{j\}\)\_\{j\\in\\mathcal\{R\}\}^\{\\top\}:=\\bm\{X\}\_\{\\mathcal\{R\}\}\-\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}X\_\{s\}\.\(18\)By Lemma[A\.1](https://arxiv.org/html/2607.05984#A1.Thmtheorem1)and Corollary[A\.3](https://arxiv.org/html/2607.05984#A1.Thmtheorem3)in Appendix[A](https://arxiv.org/html/2607.05984#A1), we obtain the following theorem\.

###### Theorem 3\.3\.

Under the genericity assumption,𝐗~ℛ\\tilde\{\\bm\{X\}\}\_\{\\mathcal\{R\}\}admits an LvLiNGAM representation whose observed DAG is the induced subgraph of𝒢O\\mathcal\{G\}^\{O\}on𝐗∖\{Xs\}\\bm\{X\}\\setminus\\\{X\_\{s\}\\\}\.

By Theorem[3\.3](https://arxiv.org/html/2607.05984#S3.Thmtheorem3), after the update \([17](https://arxiv.org/html/2607.05984#S3.E17)\),𝑿~ℛ\\tilde\{\\bm\{X\}\}\_\{\\mathcal\{R\}\}admits an LvLiNGAM representation whose observed DAG coincides with the induced subgraph of𝒢O\\mathcal\{G\}^\{O\}on𝑿∖\{Xs\}\\bm\{X\}\\setminus\\\{X\_\{s\}\\\}\. Although the above selection procedure may yield different vectors𝜶ℛ,s\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}corresponding to differentkk\-th order cumulants, all such choices induce the same observed DAG\.

Let𝒜s\\mathcal\{A\}\_\{s\}denote the set of all vectors𝜶ℛ,s\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}obtained by selecting different commonkk\-th order cumulants in the above procedure\. Although every vector in𝒜s\\mathcal\{A\}\_\{s\}induces the same observed DAG over the remaining variables, different choices may yield different latent\-to\-observed structures\. Since𝒢\\mathcal\{G\}is assumed to be the sparsest DAG in its observational equivalence class, the remaining task is to identify the vector in𝒜s\\mathcal\{A\}\_\{s\}that yields the sparsest latent\-to\-observed structure\. Specifically, for each𝜶ℛ,s∈𝒜s\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}\\in\\mathcal\{A\}\_\{s\}, define

L​\(𝜶ℛ,s\)=∑i,j∈ℛℓi​j​\(𝜶ℛ,s\),L\(\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}\)=\\sum\_\{i,j\\in\\mathcal\{R\}\}\\ell\_\{ij\}\(\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}\),whereℓi​j​\(𝜶ℛ,s\)\\ell\_\{ij\}\(\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}\)denotes the number of latent confounders between the residualized variablesX~i\\tilde\{X\}\_\{i\}andX~j\\tilde\{X\}\_\{j\}identified by Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)\. Since every vector in𝒜s\\mathcal\{A\}\_\{s\}induces the same observed DAG and the same number of latent variables, minimizingL​\(𝜶ℛ,s\)L\(\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}\)is equivalent to selecting the sparsest latent\-to\-observed structure, and hence the sparsest DAG\. Accordingly, we select any minimizer

𝜶ℛ,s∗∈arg​min𝜶ℛ,s∈𝒜s⁡L​\(𝜶ℛ,s\)\.\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}^\{\*\}\\in\\operatorname\*\{arg\\,min\}\_\{\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}\\in\\mathcal\{A\}\_\{s\}\}L\(\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}\)\.\(19\)
###### Theorem 3\.4\.

Under the genericity assumption, the sparsest latent\-to\-observed structure over𝐗~ℛ\\tilde\{\\bm\{X\}\}\_\{\\mathcal\{R\}\}obtained from𝛂ℛ,s∗\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}^\{\*\}coincides with the structure obtained from the induced subgraph of𝒢O​L\\mathcal\{G\}^\{OL\}on𝐕∖\{Xs\}\\bm\{V\}\\setminus\\\{X\_\{s\}\\\}by absorbing every latent variable having only one observed child into the disturbance of that child\.

Thus, the update \([17](https://arxiv.org/html/2607.05984#S3.E17)\) becomes

X~j=Xj−αj​s∗​Xs,∀Xj∈𝑿∖\{Xs\},\\displaystyle\\tilde\{X\}\_\{j\}=X\_\{j\}\-\\alpha^\{\*\}\_\{js\}X\_\{s\},\\quad\\forall X\_\{j\}\\in\\bm\{X\}\\setminus\\\{X\_\{s\}\\\},\(20\)where𝜶ℛ,s∗=\(αj​s∗\)j∈ℛ\\bm\{\\alpha\}^\{\*\}\_\{\\mathcal\{R\},s\}=\(\\alpha^\{\*\}\_\{js\}\)\_\{j\\in\\mathcal\{R\}\}is selected by \([19](https://arxiv.org/html/2607.05984#S3.E19)\)\. If multiple such choices of𝜶ℛ,s∗\\bm\{\\alpha\}^\{\*\}\_\{\\mathcal\{R\},s\}exist, one is selected at random\.

At this iteration, the selected aligned group gives the total\-effect column associated with the removed source in the selected sparsest representation\. The remaining aligned groups of candidate total effects, whose entries correspond to the same cumulant, are recorded as total\-effect columns of latent variables in the mixing matrix\.

By Theorems[3\.3](https://arxiv.org/html/2607.05984#S3.Thmtheorem3)and[3\.4](https://arxiv.org/html/2607.05984#S3.Thmtheorem4), updating the variables using𝜶ℛ,s∗\\bm\{\\alpha\}\_\{\\mathcal\{R\},s\}^\{\*\}yields an LvLiNGAM whose DAG coincides with the induced subgraph of𝒢\\mathcal\{G\}induced by removingXsX\_\{s\}, where every latent variable having only one observed child is absorbed into the disturbance of that child\. Applying this procedure recursively therefore recovers the entire sparsest DAG𝒢\\mathcal\{G\}\. This result is formalized in the following theorem\.

###### Theorem 3\.5\.

Under the genericity assumption, the proposed top\-down procedure identifies all total\-effect columns of the mixing matrix up to permutation consistent with the sparsest DAG over𝐕\\bm\{V\}\.

According to the selected combination of total effects in the update \([20](https://arxiv.org/html/2607.05984#S3.E20)\), the proposed method obtains the mixing matrix among the observed variables, namely\(𝑰−𝑩\)−1\(\\bm\{I\}\-\\bm\{B\}\)^\{\-1\}, corresponding to the sparsest DAG over𝑿\\bm\{X\}\. Hence,𝑩\\bm\{B\}can be estimated from\(𝑰−𝑩\)−1\(\\bm\{I\}\-\\bm\{B\}\)^\{\-1\}\. In finite samples,𝑩\\bm\{B\}can be pruned by enforcing consistency with the estimated parent\-child relationships\. Moreover, the remaining total\-effect columns corresponding to latent sources form\(𝑰−𝑩\)−1​𝚲\(\\bm\{I\}\-\\bm\{B\}\)^\{\-1\}\\bm\{\\Lambda\}, which ensures the sparsest latent\-to\-observed bipartite graph between𝑳\\bm\{L\}and𝑿\\bm\{X\}\. Multiplying\(𝑰−𝑩\)\(\\bm\{I\}\-\\bm\{B\}\)by\(𝑰−𝑩\)−1​𝚲\(\\bm\{I\}\-\\bm\{B\}\)^\{\-1\}\\bm\{\\Lambda\}yields an estimate of𝚲\\bm\{\\Lambda\}, and thus identifies the directed edges from𝑳\\bm\{L\}to𝑿\\bm\{X\}\. In finite samples,𝚲\\bm\{\\Lambda\}can be pruned by setting entries whose absolute values are below a predefined threshold to zero\.

The update \([20](https://arxiv.org/html/2607.05984#S3.E20)\) residualizes the descendants by removing the effects of the observed source\. Higher\-order cumulants are then recomputed from the residualized variables, and the procedure based on Propositions[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)–[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)is recursively applied to the residualized variables\.

In practice, higher\-order cumulants are estimated from finite samples, and their estimation accuracy generally deteriorates as the order increases\. ReLVLiNGAM updates the cumulants of descendant variables after subtracting the contributions of the disturbance and latent confounders associated with the observed source\. Since this update explicitly relies on the estimated higher\-order cumulants of these latent variables and disturbances, errors in those estimates may directly affect the updated cumulants used in subsequent total\-effect estimation\.

The proposed method also uses higher\-order cumulants, but only to match candidate total effects across different variable pairs\. Once the matching is completed, the update is performed by residualizing the observed variables themselves rather than by updating cumulants\. Consequently, the update does not explicitly rely on the estimated higher\-order cumulants of individual disturbances or latent confounders, which is expected to reduce the impact of their estimation errors on downstream inference\.

In addition, unlike ReLVLiNGAM, the proposed method does not rely on low\-order cumulants to recursively estimate disturbance cumulants\. Instead, it only requires an orderkkfor which \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) is uniquely solvable\. Specifically, the value ofkkcan be determined by increasing it fromk=2k=2and choosing the smallest one for which the corresponding linear system in \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) has full column rank\. Hence, the proposed method does not suffer from the local restriction imposed by ReLVLiNGAM\.

The proposed procedure is summarized in Algorithm[1](https://arxiv.org/html/2607.05984#algorithm1)\. We also provide a detailed example illustrating the proposed algorithm on the models in Figure[1](https://arxiv.org/html/2607.05984#S2.F1)in Appendix[D](https://arxiv.org/html/2607.05984#A4)\.

1

2

Input :Observed data matrix

𝑿∈ℝn×p\\bm\{X\}\\in\\mathbb\{R\}^\{n\\times p\}
Output :Estimated causal graph

𝒢^\\widehat\{\\mathcal\{G\}\}
3

4Initialization

5

𝒢^O←\(𝑿,∅\)\\widehat\{\\mathcal\{G\}\}^\{O\}\\leftarrow\(\\bm\{X\},\\emptyset\),

6

𝑴^←𝑰p×p\\widehat\{\\bm\{M\}\}\\leftarrow\\bm\{I\}\_\{p\\times p\},

7

𝑴^latent←∅\\widehat\{\\bm\{M\}\}\_\{\\mathrm\{latent\}\}\\leftarrow\\emptyset
8

9Identify ancestral relationships among

𝑿\\bm\{X\}and observed sources

𝑿s\\bm\{X\}\_\{s\}by Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)

10

11while*𝐗s≠∅\\bm\{X\}\_\{s\}\\neq\\emptyset*do

12

𝑿s,next←∅\\bm\{X\}\_\{s,\\mathrm\{next\}\}\\leftarrow\\emptyset
13

14foreach*Xs∈𝐗sX\_\{s\}\\in\\bm\{X\}\_\{s\}*do

15Identify

Ch~​\(Xs\)\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\),

𝑿open\\bm\{X\}\_\{\{\\mathrm\{open\}\}\}, and

𝑿closed\\bm\{X\}\_\{\{\\mathrm\{closed\}\}\}
16

17Compute all possible total effects from

XsX\_\{s\}and the corresponding disturbance cumulants by Propositions[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)and[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)

18

19while*𝐗open≠∅\\bm\{X\}\_\{\{\\mathrm\{open\}\}\}\\neq\\emptyset*do

20Compute total effects from

Ch~​\(Xs\)\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)to

𝑿open\\bm\{X\}\_\{\{\\mathrm\{open\}\}\}after regressing out

XsX\_\{s\}under each possible total effect of

XsX\_\{s\}by Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)

21

22Determine the parent–child relationships between

XsX\_\{s\}and the nodes in

𝑿open\\bm\{X\}\_\{\\mathrm\{\\mathrm\{open\}\}\}by Theorem[3\.2](https://arxiv.org/html/2607.05984#S3.Thmtheorem2)

23

24Update

Ch~​\(Xs\)\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\),

𝑿open\\bm\{X\}\_\{\{\\mathrm\{open\}\}\}, and

𝑿closed\\bm\{X\}\_\{\{\\mathrm\{closed\}\}\}
25

26end while

27

28Add all identified edges

Xs→XjX\_\{s\}\\rightarrow X\_\{j\}for

Xj∈Ch~​\(Xs\)X\_\{j\}\\in\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)into

𝒢^O\\widehat\{\\mathcal\{G\}\}^\{O\}
29

30Remove the effects of

XsX\_\{s\}from

Des​\(Xs\)\\mathrm\{Des\}\(X\_\{s\}\)using the update \([20](https://arxiv.org/html/2607.05984#S3.E20)\) with the selected total effects

31

32Replace the corresponding entries in

𝑴^\\widehat\{\\bm\{M\}\}with all selected total effects in \([20](https://arxiv.org/html/2607.05984#S3.E20)\)

33

34Append other possible total effects not used in Equation \([20](https://arxiv.org/html/2607.05984#S3.E20)\) to

𝑴^latent\\widehat\{\\bm\{M\}\}\_\{\\mathrm\{latent\}\}
35

36Identify newly emerging observed sources according to the ancestral relationships by Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)and add them to

𝑿s,next\\bm\{X\}\_\{s,\\mathrm\{next\}\}
37

38end foreach

39

40

𝑿s←𝑿s,next\\bm\{X\}\_\{s\}\\leftarrow\\bm\{X\}\_\{s,\\mathrm\{next\}\}
41

42end while

43

44

𝑩^=𝑰−𝑴^−1\\widehat\{\\bm\{B\}\}=\\bm\{I\}\-\\widehat\{\\bm\{M\}\}^\{\-1\},

𝚲^=𝑴^−1​𝑴^latent\\widehat\{\\bm\{\\Lambda\}\}=\\widehat\{\\bm\{M\}\}^\{\-1\}\\widehat\{\\bm\{M\}\}\_\{\\mathrm\{latent\}\}, and construct

𝒢^O​L\\widehat\{\\mathcal\{G\}\}^\{OL\}from

Λ^\\widehat\{\\Lambda\}
45

46return*𝒢^=𝒢^O∪𝒢^O​L\\widehat\{\\mathcal\{G\}\}=\\widehat\{\\mathcal\{G\}\}^\{O\}\\cup\\widehat\{\\mathcal\{G\}\}^\{OL\}*

47

481ex

Algorithm 1Proposed Method
## 4Simulations

\\subfigure

\[case I\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x3.png)\\subfigure\[case II\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x4.png)\\subfigure\[case III\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x5.png)

Figure 2:Three models used for simulations\.This section reports simulation results111Code is available at[https://anonymous\.4open\.science/r/Test\_20260701](https://anonymous.4open.science/r/Test_20260701)on the three causal DAGs in Figure[2](https://arxiv.org/html/2607.05984#S4.F2), which are not handled by the original ReLVLiNGAM\. We compare the proposed method with ReLVLiNGAM\(schkoda2024causal\)to examine whether it overcomes ReLVLiNGAM’s local restriction\. To evaluate the proposed method under an oracle setting, we also report results for a variant of the proposed method that is provided with the true ancestral relationships and the true number of latent confounders\. The local restriction in the original ReLVLiNGAM arises because cumulant updates rely on low\-order cumulants\. Whenk<ℓ\+1k<\\ell\+1, the system \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) becomes underdetermined and is no longer uniquely solvable\. The local restriction is imposed to avoid this situation\. This restriction can be removed by simply using sufficiently high\-order cumulants when solving \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\)\. For completeness, Appendix[C](https://arxiv.org/html/2607.05984#A3)describes a simple modification of ReLVLiNGAM that replaces low\-order cumulants with sufficiently high\-order ones when solving \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\)\. We also include this modified version in the experimental comparison\.

### 4\.1Settings

All disturbances and latent variables are sampled fromLognormal​\(−1\.1,0\.8\)\\mathrm\{Lognormal\}\(\-1\.1,0\.8\), and then centered to have zero mean\. To avoid identical distributions among the disturbances and latent variables, consistent with the genericity assumption, each variable is further multiplied by an independent scale sampled fromUniform​\(0\.9,1\.1\)\\mathrm\{Uniform\}\(0\.9,1\.1\)\. For each latent variableLiL\_\{i\}, the coefficient fromLiL\_\{i\}to its observed child with the highest causal order is fixed to one\. All other coefficients in𝚲\\bm\{\\Lambda\}and𝑩\\bm\{B\}are independently drawn fromUniform​\(0\.5,0\.8\)\\mathrm\{Uniform\}\(0\.5,0\.8\)\. Since all causal coefficients are positive, the effects along different paths cannot cancel each other out\. Moreover, the selected distributions have nonzero higher\-order cumulants\. The sample sizeNNis set to 1K, 10K, 100K, and 1M, and each experiment is repeated5050times\. We evaluate the performance of the methods using the following metrics:

- •Ntol\\mathrm\{N\}\_\{\\mathrm\{tol\}\}andNobs\\mathrm\{N\}\_\{\\mathrm\{obs\}\}: the number of runs in which the DAGs of𝒢\\mathcal\{G\}and𝒢O\\mathcal\{G\}^\{O\}are correctly recovered, respectively \(see Figure[3](https://arxiv.org/html/2607.05984#S4.F3)\);
- •PRE\\mathrm\{PRE\},REC\\mathrm\{REC\}, andF1\\mathrm\{F1\}: the average precision, recall, and F1\-score of the estimated edges of𝒢\\mathcal\{G\}\(see Figure[4](https://arxiv.org/html/2607.05984#S4.F4)\)\.

In the proposed method, we consider \([3\.2](https://arxiv.org/html/2607.05984#S3.Ex7)\) in Theorem[3\.2](https://arxiv.org/html/2607.05984#S3.Thmtheorem2)to hold in finite\-sample settings if\|αj​s−𝜷ℐ,s⊤​𝜸j,ℐ\.s\|<τ0∈\{0\.2,0\.15,0\.125,0\.1\}\|\\alpha\_\{js\}\-\\bm\{\\beta\}\_\{\\mathcal\{I\},s\}^\{\\top\}\\bm\{\\gamma\}\_\{j,\\mathcal\{I\}\.s\}\|<\\tau\_\{0\}\\in\\\{0\.2,0\.15,0\.125,0\.1\\\}for increasing sample sizes\. We further prune an edge from𝑳\\bm\{L\}to𝑿\\bm\{X\}if its absolute coefficient is belowτ0,L=0\.3\\tau\_\{0,L\}=0\.3\. For the other methods, we also applyτ0\\tau\_\{0\}andτ0,L\\tau\_\{0,L\}to the estimated coefficient matrices to prune edges\. For the original and modified ReLVLiNGAM methods, we enumerate all candidate DAGs from the estimated mixing matrices, prune edges in the estimated coefficient matrices usingτ0\\tau\_\{0\}andτ0,L\\tau\_\{0,L\}, and select the sparsest DAG\.

Although both the proposed method and ReLVLiNGAM employ Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2), they differ in how the rank ofAj,i\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{\{j\},\{i\}\}is determined\. Letσr\\sigma\_\{r\}be therr\-th largest singular value ofAj,i\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{\{j\},\{i\}\}and letτs\\tau\_\{s\}andτc​s\\tau\_\{cs\}be two predefined thresholds\. ReLVLiNGAM treatsσr\\sigma\_\{r\}as zero ifσr/σ1≤τs\\sigma\_\{r\}/\\sigma\_\{1\}\\leq\\tau\_\{s\}\. In contrast, the proposed method setsσr=0\\sigma\_\{r\}=0if1−∑i∈\[r\]σi/∑i∈\[d\]σi≤τc​s1\-\\sum\_\{i\\in\[r\]\}\\sigma\_\{i\}/\\sum\_\{i\\in\[d\]\}\\sigma\_\{i\}\\leq\\tau\_\{cs\}, whereddis the number of singular values ofAj,i\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{\{j\},\{i\}\}and is defined as in Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)\. Followingschkoda2024causal, we also impose an upper boundℓhighest\\ell\_\{\\mathrm\{highest\}\}on the number of latent variables, withℓhighest=2\\ell\_\{\\mathrm\{highest\}\}=2in cases I and III, andℓhighest=4\\ell\_\{\\mathrm\{highest\}\}=4in case II\. We setτs=0\.008​\(i−1\)/N0\.125\\tau\_\{s\}=0\.008\(i\-1\)/N^\{0\.125\}for the original ReLVLiNGAM andτc​s=0\.002\+0\.0005​\(i−1\)\\tau\_\{cs\}=0\.002\+0\.0005\(i\-1\)for our method, whereiidenotes the depth from the observed source to reflect increasing estimation error with depth\. The settings of the modified ReLVLiNGAM follow those of the proposed method\. To estimate the cumulants of the disturbances and latent variables, both the proposed method and the modified ReLVLiNGAM increasekkfromk=2k=2until the system in \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) is uniquely solvable, and then use the resulting value ofkkin Proposition[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)\.

![Refer to caption](https://arxiv.org/html/2607.05984v1/x6.png)Figure 3:The performances inNtol\\mathrm\{N\}\_\{\\mathrm\{tol\}\}andNobs\\mathrm\{N\}\_\{\\mathrm\{obs\}\}\.![Refer to caption](https://arxiv.org/html/2607.05984v1/x7.png)Figure 4:The performances of PRE, REC, and F1\-score\.
### 4\.2Discussion

Figure[3](https://arxiv.org/html/2607.05984#S4.F3)shows that the proposed method outperforms the original ReLVLiNGAM on all evaluation metrics across all settings\. Since the DAGs in Figure[2](https://arxiv.org/html/2607.05984#S4.F2)do not satisfy the local restriction, the original ReLVLiNGAM fails to recover the correct DAG even as the sample size increases\. In contrast, the estimation accuracy of the proposed method improves steadily with increasing sample size\.

For reference, we also report the performance of the oracle version of the proposed method\. Compared with the standard version, the oracle version achieves substantially higher accuracy when the DAG is sparse or the sample size is small\. This result suggests that the accuracy of estimating ancestral relationships and the number of pairwise confounders has a non\-negligible impact on the accuracy of DAG recovery\.

The boxplots of PRE, REC, and F1 score in Figure[4](https://arxiv.org/html/2607.05984#S4.F4)show that, for the proposed method, both the mean \(triangles\) and the median \(horizontal bars\) of these metrics increase asNNincreases, indicating increasingly accurate recovery of parent–child relationships over𝑿\\bm\{X\}\.

The modified ReLVLiNGAM achieves performance comparable to that of the proposed method when the sample size is large in case II, but performs worse in sparse settings \(cases I and III\), likely because it lacks an effective edge\-pruning strategy from finite samples\. It is generally inferior to the proposed method in recovering the𝑳→𝑿\\bm\{L\}\\\!\\to\\\!\\bm\{X\}edges, except in case II with 100K samples\. This may be because errors in estimating higher\-order disturbance cumulants can propagate to subsequent cumulant computations\. In contrast, the proposed method uses estimated disturbance cumulants only to match candidate total effects with the same associated cumulant\. Since the cumulants are not recursively updated and propagated to subsequent iterations, the effect of cumulant estimation errors is expected to be less severe\.

## 5Real Data

We further evaluate the practical usefulness of the proposed method by applying it, together with ParceLiNGAM\(tashiro2014parcelingam\), RCD\(Maeda2020\), and the original and modified versions of ReLVLiNGAM\(schkoda2024causal\), to the Sachs protein dataset fromSachs2005\.

\\subfigure

\[Reference DAG\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x8.png)\\subfigure\[Canonical model\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x9.png)\\subfigure\[ParceLiNGAM\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x10.png)\\subfigure\[RCD\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x11.png)\\subfigure\[ReLL\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x12.png)\\subfigure\[modified ReLL\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x13.png)\\subfigure\[Proposed method\]![Refer to caption](https://arxiv.org/html/2607.05984v1/x14.png)

Figure 5:The results of different methods applied to the Sachs dataset\.The Sachs dataset\(Sachs2005\)records expression levels of phosphorylated proteins and phospholipids in human immune cells and contains 11 variables and 7,467 samples\. The dataset is accompanied by a reference signaling network constructed from biological knowledge\. Based on this reference network, we treat PKC and PKA as latent confounders and use Raf, Mek, Erk, and Akt as observed variables\. The DAG for the model is shown in Figure[5](https://arxiv.org/html/2607.05984#S5.F5)\(a\)\. Although the DAG in Figure[5](https://arxiv.org/html/2607.05984#S5.F5)\(a\) is not a canonical LvLiNGAM, it can be transformed into its canonical form in Figure[5](https://arxiv.org/html/2607.05984#S5.F5)\(b\) without changing the observed structure by applying Algorithm A ofhoyer2008estimation\.

The thresholdsτs\\tau\_\{s\}andτc​s\\tau\_\{cs\}for the proposed method and ReLVLiNGAMs are set in the same manner as in the numerical experiments in Section[4](https://arxiv.org/html/2607.05984#S4)\. The upper boundℓhighest\\ell\_\{\\mathrm\{highest\}\}is set to22for both methods\. Here, instead of employing a hard threshold, we verify \([3\.2](https://arxiv.org/html/2607.05984#S3.Ex7)\) in Theorem[3\.2](https://arxiv.org/html/2607.05984#S3.Thmtheorem2)using 99% bootstrap confidence intervals constructed from 400 bootstrap resamples\. ParceLiNGAM and RCD use the Hilbert–Schmidt independence criterion \(HSIC;Gretton2007\) to infer causal directions among observed variables\. RCD sets the HSIC significance level to 0\.01\. ParceLiNGAM additionally applies Fisher’s method to combine HSIC p\-values and uses a significance level of 0\.1 for Fisher’s test\. RCD also uses the Pearson test and the Shapiro\-Wilk test, both at the 0\.01 significance level\. For the original ReLVLiNGAM, since it outputs a very dense𝑳→𝑿\\bm\{L\}\\\!\\to\\\!\\bm\{X\}structure, we prune latent\-to\-observed edges whose absolute coefficients are below0\.10\.1\.

Figures[5](https://arxiv.org/html/2607.05984#S5.F5)\(c\)–\(g\) show the DAGs estimated by each method\. In \(e\) and \(f\), “ReLL” denotes ReLVLiNGAM\. As can be seen from these figures, the proposed method recovers an observed DAG that is closer to that in Figure[5](https://arxiv.org/html/2607.05984#S5.F5)\(a\) with only one extra edge,Mek→Akt\\mathrm\{Mek\}\\to\\mathrm\{Akt\}, although it estimates a denser𝑳→𝑿\\bm\{L\}\\\!\\to\\\!\\bm\{X\}structure\. ParceLiNGAM fails to identify the existence of latent confounders, and outputs redundant edgesRaf→Erk,Raf→Akt,Mek→Akt\{\\text\{Raf\}\\to\\text\{Erk\},\\text\{Raf\}\\to\\text\{Akt\},\\text\{Mek\}\\to\\text\{Akt\}\}\. RCD correctly concludes that there is no edge betweenRafandErk, but does not identify any directed edges among the observed variables\. Even with the pruning, the original ReLVLiNGAM yields an incorrect causal order over𝑿\\bm\{X\}and a dense𝑳→𝑿\\bm\{L\}\\\!\\to\\\!\\bm\{X\}structure\. The modified ReLVLiNGAM yields an𝑳→𝑿\\bm\{L\}\\\!\\to\\\!\\bm\{X\}structure closer to the reference graph than the original ReLVLiNGAM, but also fails to recover the true ancestral relationships\.

Overall, the proposed method recovers the observed DAG with only one extra edge,Mek→Akt\\mathrm\{Mek\}\\to\\mathrm\{Akt\}\. Although the original ReLVLiNGAM, the modified ReLVLiNGAM, and the proposed method all infer multiple latent variables, the proposed method and the modified ReLVLiNGAM produce the sparsest𝑳→𝑿\\bm\{L\}\\\!\\to\\\!\\bm\{X\}structure and most closely match the structure from latent variables to observed variables in Figure[5](https://arxiv.org/html/2607.05984#S5.F5)\(a\)\.

## 6Conclusion

In this paper, we proposed a method for recovering the sparsest causal DAG within the observational equivalence class from finite samples\. Unlike the original ReLVLiNGAM, the proposed method does not require the local restriction and is applicable to general canonical LvLiNGAMs\. Although ReLVLiNGAM consistently estimates the mixing matrix, recovering the sparsest DAG asymptotically still requires an appropriate permutation of its columns\. In contrast, the proposed method consistently estimates the sparsest DAG\.

The simulation results and the application to the Sachs protein data demonstrate the superiority of the proposed method over ReLVLiNGAM\. Although the modified ReLVLiNGAM can also recover the sparsest causal DAG from finite samples without requiring the local restriction, the proposed method exhibits better finite\-sample performance\.

Several limitations remain\. Since the proposed method relies on the estimation of higher\-order cumulants, its performance can deteriorate when the sample size is small or the data are noisy, which may degrade the accuracy of DAG estimation\. In addition, the computational cost can still be high when multiple candidate total effects must be examined, although this issue might not be severe when the DAG is sparse\.

Improving the accuracy and computational efficiency of the method is an important direction for future work\.

\\acks

This work was supported by JST SPRING under Grant Number JPMJSP2110 and JSPS KAKENHI under Grant Numbers 25K15017\.

## References

## Appendix APreservation of the LvLiNGAM Structure under Source Removal

Assume thatXsX\_\{s\}is an observed source of𝒢\\mathcal\{G\}\. Letℛ=\[p\]∖\{s\}\\mathcal\{R\}=\[p\]\\setminus\\\{s\\\}and fixus∈\{es\}∪Pa​\(Xs\)u\_\{s\}\\in\\\{e\_\{s\}\\\}\\cup\\mathrm\{Pa\}\(X\_\{s\}\)\. As discussed in Section[2\.2](https://arxiv.org/html/2607.05984#S2.SS2), the corresponding DAG might change when swappingusu\_\{s\}andese\_\{s\}\. Let𝒢\(us\)\\mathcal\{G\}^\{\(u\_\{s\}\)\}denote the DAG in the observational equivalence class of𝒢\\mathcal\{G\}that is obtained by swappingese\_\{s\}andusu\_\{s\}\. For eachj∈ℛj\\in\\mathcal\{R\}, choose a rootαj​s\(us\)\\alpha^\{\(u\_\{s\}\)\}\_\{js\}returned by Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)whose associatedkk\-th order cumulant isκ\(k\)​\(us\)\\kappa^\{\(k\)\}\(u\_\{s\}\)\. Then the update rule \([18](https://arxiv.org/html/2607.05984#S3.E18)\) is expressed as

𝑿~ℛ\(us\)=𝑿ℛ−𝜶ℛ,s\(us\)​Xs,𝜶ℛ,s\(us\)=\(αj​s\(us\)\)j∈ℛ\.\\displaystyle\\tilde\{\\bm\{X\}\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}=\\bm\{X\}\_\{\\mathcal\{R\}\}\-\\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}X\_\{s\},\\qquad\\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}=\\left\(\\alpha^\{\(u\_\{s\}\)\}\_\{js\}\\right\)\_\{j\\in\\mathcal\{R\}\}\.
###### Lemma A\.1\.

Under the genericity assumption,𝐗~ℛ\(us\)\\tilde\{\\bm\{X\}\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}admits an LvLiNGAM representation whose DAG is the induced subgraph of𝒢\(us\)\\mathcal\{G\}^\{\(u\_\{s\}\)\}obtained by deletingXsX\_\{s\}and all edges linked toXsX\_\{s\}\.

###### Proof A\.2\.

Denote the LvLiNGAM representation associated with𝒢\(us\)\\mathcal\{G\}^\{\(u\_\{s\}\)\}by

𝑿=𝚲\(us\)​𝑳\(us\)\+𝑩\(us\)​𝑿\+𝒆\(us\)\.\\bm\{X\}=\\bm\{\\Lambda\}^\{\(u\_\{s\}\)\}\\bm\{L\}^\{\(u\_\{s\}\)\}\+\\bm\{B\}^\{\(u\_\{s\}\)\}\\bm\{X\}\+\\bm\{e\}^\{\(u\_\{s\}\)\}\.\(21\)αj​s\(us\)\\alpha^\{\(u\_\{s\}\)\}\_\{js\}is the total effect fromusu\_\{s\}toXjX\_\{j\}\. Let𝐁ℛ,ℛ\(us\)\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}denote the submatrix of𝐁\(us\)\\bm\{B\}^\{\(u\_\{s\}\)\}with rows and columns indexed byℛ\\mathcal\{R\}\. Let𝐛ℛ,s\\bm\{b\}\_\{\\mathcal\{R\},s\}denote the column of𝐁ℛ,ℛ\(us\)\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}corresponding toXsX\_\{s\}\. By the definition of total effects,

𝜶ℛ,s\(us\)=𝒃ℛ,s\(us\)\+𝑩ℛ,ℛ\(us\)​𝜶ℛ,s\(us\)\.\\displaystyle\\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}=\\bm\{b\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}\+\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}\.\(22\)Let𝚲ℛ,:\(us\)\\bm\{\\Lambda\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},:\}denote the submatrix of𝚲\(us\)\\bm\{\\Lambda\}^\{\(u\_\{s\}\)\}consisting of the rows indexed byℛ\\mathcal\{R\}and let𝐞ℛ\\bm\{e\}\_\{\\mathcal\{R\}\}denote the subvector of𝐞\(us\)\\bm\{e\}^\{\(u\_\{s\}\)\}consisting of the entries indexed byℛ\\mathcal\{R\}\. Restricting the structural equations toℛ\\mathcal\{R\}, we have

𝑿ℛ=𝒃ℛ,s\(us\)​Xs\+𝑩ℛ,ℛ\(us\)​𝑿ℛ\+𝚲ℛ,:\(us\)​𝑳\(us\)\+𝒆ℛ\(us\)\.\\bm\{X\}\_\{\\mathcal\{R\}\}=\\bm\{b\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}X\_\{s\}\+\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\\bm\{X\}\_\{\\mathcal\{R\}\}\+\\bm\{\\Lambda\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},:\}\\bm\{L\}^\{\(u\_\{s\}\)\}\+\\bm\{e\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}\.Since𝐗ℛ=𝐗~ℛ\(us\)\+𝛂ℛ,s\(us\)​Xs\\bm\{X\}\_\{\\mathcal\{R\}\}=\\tilde\{\\bm\{X\}\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}\+\\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}X\_\{s\}, we obtain

𝑿~ℛ\(us\)\\displaystyle\\tilde\{\\bm\{X\}\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}=−𝜶ℛ,s\(us\)​Xs\+𝒃ℛ,s\(us\)​Xs\+𝑩ℛ,ℛ\(us\)​\(𝑿~ℛ\(us\)\+𝜶ℛ,s\(us\)​Xs\)\+𝚲ℛ,:\(us\)​𝑳\(us\)\+𝒆ℛ\(us\)\\displaystyle=\-\\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}X\_\{s\}\+\\bm\{b\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}X\_\{s\}\+\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\(\\tilde\{\\bm\{X\}\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}\+\\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}X\_\{s\}\)\+\\bm\{\\Lambda\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},:\}\\bm\{L\}^\{\(u\_\{s\}\)\}\+\\bm\{e\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}=𝑩ℛ,ℛ\(us\)​𝑿~ℛ\(us\)\+\(𝒃ℛ,s\(us\)\+𝑩ℛ,ℛ\(us\)​𝜶ℛ,s\(us\)−𝜶ℛ,s\(us\)\)​Xs\+𝚲ℛ,:\(us\)​𝑳\(us\)\+𝒆ℛ\(us\)\.\\displaystyle=\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\\tilde\{\\bm\{X\}\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}\+\\left\(\\bm\{b\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}\+\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}\-\\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}\\right\)X\_\{s\}\+\\bm\{\\Lambda\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},:\}\\bm\{L\}^\{\(u\_\{s\}\)\}\+\\bm\{e\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}\.By \([22](https://arxiv.org/html/2607.05984#A1.E22)\), the coefficient ofXsX\_\{s\}is zero\. Thus,

𝑿~ℛ\(us\)=𝑩ℛ,ℛ\(us\)​𝑿~ℛ\(us\)\+𝚲ℛ,:\(us\)​𝑳\(us\)\+𝒆ℛ\(us\)\.\\tilde\{\\bm\{X\}\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}=\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\\tilde\{\\bm\{X\}\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}\+\\bm\{\\Lambda\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},:\}\\bm\{L\}^\{\(u\_\{s\}\)\}\+\\bm\{e\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\}\}\.\(23\)This is an LvLiNGAM representation over𝐗ℛ\\bm\{X\}\_\{\\mathcal\{R\}\}, where its coefficient matrix among observed variables is exactly𝐁ℛ,ℛ\(us\)\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}, and its coefficient matrix from𝐋\\bm\{L\}is𝚲ℛ,:\(us\)\\bm\{\\Lambda\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},:\}\. Hence, the corresponding DAG is obtained from𝒢\(us\)\\mathcal\{G\}^\{\(u\_\{s\}\)\}by deletingXsX\_\{s\}and all edges linked toXsX\_\{s\}, namely the induced subgraph of𝒢\(us\)\\mathcal\{G\}^\{\(u\_\{s\}\)\}after deletingXsX\_\{s\}\.

Based on Lemma[A\.1](https://arxiv.org/html/2607.05984#A1.Thmtheorem1), we have the following Corollary\.

###### Corollary A\.3\.

Let𝐁ℛ,ℛ\\bm\{B\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}be the coefficient matrix among𝐗∖\{Xs\}\\bm\{X\}\\setminus\\\{X\_\{s\}\\\}in𝒢\(es\)\\mathcal\{G\}^\{\(e\_\{s\}\)\}\. Then, for eachus∈Pa​\(Xs\)∪\{es\}u\_\{s\}\\in\\mathrm\{Pa\}\(X\_\{s\}\)\\cup\\\{e\_\{s\}\\\},

𝑩ℛ,ℛ\(us\)=𝑩ℛ,ℛ\.\\displaystyle\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}=\\bm\{B\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\.That is, the observed parent\-child relationships among𝐗ℛ\\bm\{X\}\_\{\\mathcal\{R\}\}do not depend onusu\_\{s\}\.

###### Proof A\.4\.

The mixing matrix for𝐗\\bm\{X\}is expressed as

\(𝑰−𝑩\)−1=\[1𝟎𝒎ℛ,sO\(I−𝑩ℛ,ℛ\)−1\],\\displaystyle\(\\bm\{I\}\-\\bm\{B\}\)^\{\-1\}=\\left\[\\begin\{array\}\[\]\{cc\}1&\\bm\{0\}\\\\ \\bm\{m\}^\{O\}\_\{\\mathcal\{R\},s\}&\(I\-\\bm\{B\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\)^\{\-1\}\\end\{array\}\\right\],where𝐦ℛ,sO\\bm\{m\}^\{O\}\_\{\\mathcal\{R\},s\}is the vector of total effects fromXsX\_\{s\}to𝐗∖\{Xs\}\\bm\{X\}\\setminus\\\{X\_\{s\}\\\}\. For anyus∈Pa​\(Xs\)∪\{es\}u\_\{s\}\\in\\mathrm\{Pa\}\(X\_\{s\}\)\\cup\\\{e\_\{s\}\\\}, the corresponding swap affects only the total effects associated withXsX\_\{s\}\. Thus,\(𝐈−𝐁\(us\)\)−1\(\\bm\{I\}\-\\bm\{B\}^\{\(u\_\{s\}\)\}\)^\{\-1\}is obtained from\(𝐈−𝐁\)−1\(\\bm\{I\}\-\\bm\{B\}\)^\{\-1\}by replacing the column corresponding toXsX\_\{s\}with\(1,αℛ,s\(us\)\)⊤\(1,\\alpha^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}\)^\{\\top\}

\(𝑰−𝑩\(us\)\)−1=\[1𝟎𝜶ℛ,s\(us\)\(𝑰−𝑩ℛ,ℛ\)−1\],\\displaystyle\(\\bm\{I\}\-\\bm\{B\}^\{\(u\_\{s\}\)\}\)^\{\-1\}=\\begin\{bmatrix\}1&\\bm\{0\}\\\\ \\bm\{\\alpha\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}&\(\\bm\{I\}\-\\bm\{B\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\)^\{\-1\}\\end\{bmatrix\},where this replacement is trivial whenus=esu\_\{s\}=e\_\{s\}\.

Computing the inverse gives

𝑰−𝑩\(us\)=\[1𝟎−\(𝑰−𝑩ℛ,ℛ\)​αℛ,s\(us\)𝑰−𝑩ℛ,ℛ\]\.\\displaystyle\\bm\{I\}\-\\bm\{B\}^\{\(u\_\{s\}\)\}=\\begin\{bmatrix\}1&\\bm\{0\}\\\\ \-\(\\bm\{I\}\-\\bm\{B\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\)\\alpha^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}&\\bm\{I\}\-\\bm\{B\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\\end\{bmatrix\}\.Thus, we can obtain𝐁\(us\)\\bm\{B\}^\{\(u\_\{s\}\)\}

𝑩\(us\)=\[0𝟎\(𝑰−𝑩ℛ,ℛ\)​αℛ,s\(us\)𝑩ℛ,ℛ\],\\displaystyle\\bm\{B\}^\{\(u\_\{s\}\)\}=\\begin\{bmatrix\}0&\\bm\{0\}\\\\ \(\\bm\{I\}\-\\bm\{B\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\)\\alpha^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},s\}&\\bm\{B\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}\\end\{bmatrix\},implying𝐁ℛ,ℛ\(us\)=𝐁ℛ,ℛ\\bm\{B\}^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}=\\bm\{B\}\_\{\\mathcal\{R\},\\mathcal\{R\}\}from the lower\-right block\.

## Appendix BProofs of Theorems in Section[3](https://arxiv.org/html/2607.05984#S3)

Proof of Lemma[3\.1](https://arxiv.org/html/2607.05984#S3.Thmtheorem1)

Followingdrton2011global, the total effect from an observed sourceXsX\_\{s\}toXjX\_\{j\}is written as

mj​sO=bj​s\+∑π∈𝒫​\(Xs,Xj\)∖\{\(Xs,Xj\)\}∏l,h:\(Xl,Xh\)∈πbh​l\.\\displaystyle m^\{O\}\_\{js\}=b\_\{js\}\+\\sum\_\{\\pi\\in\\mathcal\{P\}\(X\_\{s\},X\_\{j\}\)\\setminus\\\{\(X\_\{s\},X\_\{j\}\)\\\}\}\\penalty 10000\\ \\prod\_\{l,h:\(X\_\{l\},X\_\{h\}\)\\in\\pi\}b\_\{hl\}\.Since every variable inCh~​\(Xs\)\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)is a descendant ofXsX\_\{s\}, and every node in𝑿open\\bm\{X\}\_\{\\mathrm\{open\}\}has ancestors that are either contained in\{Xs\}∪Ch~​\(Xs\)\\\{X\_\{s\}\\\}\\cup\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)or not influenced byXsX\_\{s\}, every directed path fromXsX\_\{s\}toXjX\_\{j\}other than the direct edgeXs→XjX\_\{s\}\\to X\_\{j\}must first pass through some nodeXi∈Ch~​\(Xs\)X\_\{i\}\\in\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)\. Thus,

mj​sO=bj​s\+∑i:Xi∈Ch~​\(Xs\)bi​s​∑π∈𝒫​\(Xi,Xj\)∏l,h:\(Xl,Xh\)∈πbh​l=bj​s\+∑i:Xi∈Ch~​\(Xs\)bi​s​mj​iO,\\displaystyle m^\{O\}\_\{js\}=b\_\{js\}\+\\sum\_\{i:X\_\{i\}\\in\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)\}b\_\{is\}\\sum\_\{\\pi\\in\\mathcal\{P\}\(X\_\{i\},X\_\{j\}\)\}\\penalty 10000\\ \\prod\_\{l,h:\(X\_\{l\},X\_\{h\}\)\\in\\pi\}b\_\{hl\}=b\_\{js\}\+\\sum\_\{i:X\_\{i\}\\in\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{s\}\)\}b\_\{is\}\\penalty 10000\\ m^\{O\}\_\{ji\},which establishes \([15](https://arxiv.org/html/2607.05984#S3.E15)\)\.■\\hfill\\blacksquare

Proof of Theorem[3\.2](https://arxiv.org/html/2607.05984#S3.Thmtheorem2)

By Lemma[3\.1](https://arxiv.org/html/2607.05984#S3.Thmtheorem1), the true direct effect satisfies

bj​s=mj​sO−∑i∈ℐbi​s​mj​i\.sO\.\\displaystyle b\_\{js\}=m^\{O\}\_\{js\}\-\\sum\_\{i\\in\\mathcal\{I\}\}b\_\{is\}m^\{O\}\_\{ji\.s\}\.Suppose first thatXs∉Pa​\(Xj\)X\_\{s\}\\notin\\mathrm\{Pa\}\(X\_\{j\}\)\. Thenbj​s=0b\_\{js\}=0\. Since the true valuesmj​sOm^\{O\}\_\{js\},bi​sb\_\{is\}, andmj​i\.sOm^\{O\}\_\{ji\.s\}are contained in the candidate sets𝒎j​sO\\bm\{m\}^\{O\}\_\{js\},𝒃i​s\\bm\{b\}\_\{is\}, and𝒎j​i\.sO\\bm\{m\}^\{O\}\_\{ji\.s\}, respectively, there exists a choice

αj​s∈𝒎j​sO,βi​s∈𝒃i​s,γj​i\.s∈𝒎j​i\.sO,\\displaystyle\\alpha\_\{js\}\\in\\bm\{m\}^\{O\}\_\{js\},\\qquad\\beta\_\{is\}\\in\\bm\{b\}\_\{is\},\\qquad\\gamma\_\{ji\.s\}\\in\\bm\{m\}^\{O\}\_\{ji\.s\},such thatαj​s\\alpha\_\{js\}andβi​s\\beta\_\{is\},i∈ℐi\\in\\mathcal\{I\}, correspond to the samekk\-th order cumulant and

αj​s−𝜷ℐ,s⊤​𝜸j,ℐ\.s=mj​sO−∑i∈ℐbi​s​mj​i\.sO=0\.\\displaystyle\\alpha\_\{js\}\-\\bm\{\\beta\}\_\{\\mathcal\{I\},s\}^\{\\top\}\\bm\{\\gamma\}\_\{j,\\mathcal\{I\}\.s\}=m^\{O\}\_\{js\}\-\\sum\_\{i\\in\\mathcal\{I\}\}b\_\{is\}m^\{O\}\_\{ji\.s\}=0\.Thus, \([3\.2](https://arxiv.org/html/2607.05984#S3.Ex7)\) holds\.

Conversely, suppose thatXs∈Pa​\(Xj\)X\_\{s\}\\in\\mathrm\{Pa\}\(X\_\{j\}\)\. Thenbj​s≠0b\_\{js\}\\neq 0\. Under the genericity assumption, no choice of candidates from𝒎j​sO\\bm\{m\}^\{O\}\_\{js\},𝒃i​s\\bm\{b\}\_\{is\}, and𝒎j​i\.sO\\bm\{m\}^\{O\}\_\{ji\.s\}, withαj​s\\alpha\_\{js\}andβi​s\\beta\_\{is\}corresponding to the samekk\-th order cumulant, can satisfy

αj​s−𝜷ℐ,s⊤​𝜸j,ℐ\.s=0\.\\displaystyle\\alpha\_\{js\}\-\\bm\{\\beta\}\_\{\\mathcal\{I\},s\}^\{\\top\}\\bm\{\\gamma\}\_\{j,\\mathcal\{I\}\.s\}=0\.Hence, \([3\.2](https://arxiv.org/html/2607.05984#S3.Ex7)\) fails generically\.

Therefore, \([3\.2](https://arxiv.org/html/2607.05984#S3.Ex7)\) holds if and only ifXs∉Pa​\(Xj\)X\_\{s\}\\notin\\mathrm\{Pa\}\(X\_\{j\}\)\.■\\hfill\\blacksquare

Proof of Theorem[3\.3](https://arxiv.org/html/2607.05984#S3.Thmtheorem3)

Immediate from Lemma[A\.1](https://arxiv.org/html/2607.05984#A1.Thmtheorem1)and Corollary[A\.3](https://arxiv.org/html/2607.05984#A1.Thmtheorem3)\.■\\hfill\\blacksquare

Proof of Theorem[3\.4](https://arxiv.org/html/2607.05984#S3.Thmtheorem4)

Letκ\(k\)​\(us\)\\kappa^\{\(k\)\}\(u\_\{s\}\)denote the cumulant corresponding to𝜶ℛ,s∗\\bm\{\\alpha\}^\{\*\}\_\{\\mathcal\{R\},s\}\. By \([21](https://arxiv.org/html/2607.05984#A1.E21)\), the corresponding latent\-to\-observed coefficient matrix isΛ\(us\)\\Lambda^\{\(u\_\{s\}\)\}\. Since𝜶ℛ,s∗\\bm\{\\alpha\}^\{\*\}\_\{\\mathcal\{R\},s\}is chosen to minimize the number of latent parents, the support ofΛ\(us\)\\Lambda^\{\(u\_\{s\}\)\}is the sparsest among all latent\-to\-observed structures in the observational equivalence class\. Therefore, it coincides with the structure of𝒢O​L\\mathcal\{G\}^\{OL\}\.

By \([23](https://arxiv.org/html/2607.05984#A1.E23)\), the latent\-to\-observed coefficient matrix for the model of𝑿~ℛ\\tilde\{\\bm\{X\}\}\_\{\\mathcal\{R\}\}isΛℛ,:\(us\)\\Lambda^\{\(u\_\{s\}\)\}\_\{\\mathcal\{R\},:\}\. Hence, its support coincides with that of the induced subgraph of𝒢O​L\\mathcal\{G\}^\{OL\}on𝑽∖\{Xs\}\\bm\{V\}\\setminus\\\{X\_\{s\}\\\}obtained by absorbing every latent variable having only one observed child into the disturbance of that child\.■\\hfill\\blacksquare

Proof of Theorem[3\.5](https://arxiv.org/html/2607.05984#S3.Thmtheorem5)

Immediate from Theorem[3\.3](https://arxiv.org/html/2607.05984#S3.Thmtheorem3)and Theorem[3\.4](https://arxiv.org/html/2607.05984#S3.Thmtheorem4)\.■\\hfill\\blacksquare

## Appendix CA Modification of the Original ReLVLiNGAM

As mentioned earlier, the origianl ReLVLiNGAM, proposed byschkoda2024causal, cannot be applied when the local restriction is violated\. In this section, we modify the original ReLVLiNGAM so that it can be applied even when the local restriction is violated\.

Like the proposed method, ReLVLiNGAM is a top\-down algorithm that recursively estimates total effects\. Here, letX1X\_\{1\}be an observed source and letXi∈Des​\(X1\)X\_\{i\}\\in\\mathrm\{Des\}\(X\_\{1\}\)\. The update rule forXiX\_\{i\}in ReLVLiNGAM is given by

Xi←Xi−mi​1O​e1−∑h:Lh∈Pa​\(X1\)mi​hO​L​Lh\.\\displaystyle X\_\{i\}\\leftarrow X\_\{i\}\-m\_\{i1\}^\{O\}e\_\{1\}\-\\sum\_\{h:L\_\{h\}\\in\\mathrm\{Pa\}\(X\_\{1\}\)\}m\_\{ih\}^\{OL\}L\_\{h\}\.\(24\)Although this update cannot be computed directly because neither the disturbances nor the latent variables are observed, the higher\-order cumulants of the updated variables satisfy

ci1,…,ik\(k\)←ci1,…,ik\(k\)−mi1​1O​⋯​mik​1O​κ\(k\)​\(e1\)−∑h:Lh∈Pa​\(X1\)mi1​hO​L​⋯​mik​hO​L​κ\(k\)​\(Lh\)\.\\displaystyle c^\{\(k\)\}\_\{i\_\{1\},\\ldots,i\_\{k\}\}\\leftarrow c^\{\(k\)\}\_\{i\_\{1\},\\ldots,i\_\{k\}\}\-m^\{O\}\_\{i\_\{1\}1\}\\cdots m^\{O\}\_\{i\_\{k\}1\}\\kappa^\{\(k\)\}\(e\_\{1\}\)\-\\sum\_\{h:L\_\{h\}\\in\\mathrm\{Pa\}\(X\_\{1\}\)\}m^\{OL\}\_\{i\_\{1\}h\}\\cdots m^\{OL\}\_\{i\_\{k\}h\}\\kappa^\{\(k\)\}\(L\_\{h\}\)\.Hence, once the higher\-order cumulants of the disturbances and latent variables have been estimated, the higher\-order cumulants after the update can be computed without explicitly updating the observed variables\.

For a source nodeX1X\_\{1\},schkoda2024causalproposed estimating

\(κ\(k\)​\(e1\),κ\(k\)​\(L1\),…,κ\(k\)​\(Lq\)\)\(\\kappa^\{\(k\)\}\(e\_\{1\}\),\\kappa^\{\(k\)\}\(L\_\{1\}\),\\ldots,\\kappa^\{\(k\)\}\(L\_\{q\}\)\)by solving the linear system

\[11⋯1m21Om21O​L⋯m2​qO​L⋮⋮⋱⋮mp​1Omp​1O​L⋯mp​qO​L\]​\[κ\(k\)​\(e1\)κ\(k\)​\(L1\)⋮κ\(k\)​\(Lq\)\]=\[c1​⋯​11\(k\)c1​⋯​12\(k\)⋮c1​⋯​1​p\(k\)\]\.\\begin\{bmatrix\}1&1&\\cdots&1\\\\ m\_\{21\}^\{O\}&m\_\{21\}^\{OL\}&\\cdots&m\_\{2q\}^\{OL\}\\\\ \\vdots&\\vdots&\\ddots&\\vdots\\\\ m\_\{p1\}^\{O\}&m\_\{p1\}^\{OL\}&\\cdots&m\_\{pq\}^\{OL\}\\end\{bmatrix\}\\begin\{bmatrix\}\\kappa^\{\(k\)\}\(e\_\{1\}\)\\\\ \\kappa^\{\(k\)\}\(L\_\{1\}\)\\\\ \\vdots\\\\ \\kappa^\{\(k\)\}\(L\_\{q\}\)\\end\{bmatrix\}=\\begin\{bmatrix\}c^\{\(k\)\}\_\{1\\cdots 11\}\\\\ c^\{\(k\)\}\_\{1\\cdots 12\}\\\\ \\vdots\\\\ c^\{\(k\)\}\_\{1\\cdots 1p\}\\end\{bmatrix\}\.However, when the local restriction is violated, this linear system becomes underdetermined and therefore cannot be solved\. By contrast, Proposition[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)enables the estimation of the higher\-order cumulants of the disturbance and the latent confounders even without the local restriction\. Inschkoda2024causal, Proposition[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)is used solely to establish the correspondence between total effects and higher\-order cumulants\. Once this correspondence has been identified, the higher\-order cumulants can in turn be estimated, making it possible to compute the cumulant update above\. Consequently, the top\-down update procedure remains applicable even when the local restriction is violated\.

The overall procedure of the modified ReLVLiNGAM is as follows\. At each iteration, it first constructsAj,i\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{j,i\}using sufficiently high\-order cumulants and applies Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)to estimate the number of latent confounders and the ancestral relationships among the observed variables\. Based on the estimated ancestral relationships, the observed source nodes are identified\. For each observed source, the candidate total effects from the source to the remaining observed variables are estimated by extending Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)to the present setting\. Next, Proposition[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)is used to estimate the higher\-order cumulants of the disturbances and latent variables whenever they are identifiable\. These estimates are then used to update the cumulants of the residualized variables\. The same procedure is repeated recursively until all columns of the mixing matrix corresponding to the total effects have been estimated\.

To avoid using unnecessarily high\-order cumulants, we choosek1k\_\{1\}andk2k\_\{2\}to be the smallest values that satisfy the requirements of Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)\. For each descendantXiX\_\{i\}of an observed sourceX1X\_\{1\}, letki,1k\_\{i,1\}denote the smallest order for which the linear system in Proposition[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)is generically solvable\.ki,1k\_\{i,1\}can be estimated by increasingkkfrom22and checking whether the corresponding linear system \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) is generically solvable\. Applying the same procedure to every descendant ofX1X\_\{1\}yields the corresponding valueski,1k\_\{i,1\}\.

Now suppose thatXiX\_\{i\}is the next observed source after removingX1X\_\{1\}\. For another descendantXjX\_\{j\}ofX1X\_\{1\}, let

k1=max⁡\{ℓ\+2,ki,1,kj,1\},k2=k1\+⌈−3\+8​ℓ\+172⌉,\\displaystyle k\_\{1\}=\\max\\\{\\ell\+2,k\_\{i,1\},k\_\{j,1\}\\\},\\qquad k\_\{2\}=k\_\{1\}\+\\left\\lceil\\frac\{\-3\+\\sqrt\{8\\ell\+17\}\}\{2\}\\right\\rceil,\(25\)and constructAi,j\(k1,k2\)A^\{\(k\_\{1\},k\_\{2\}\)\}\_\{i,j\}as in \([12](https://arxiv.org/html/2607.05984#S2.E12)\)\. Starting fromℓ=0\\ell=0,ℓ\\ellis estimated by increasingℓ\\elluntil the rank condition in Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)is satisfied\.

Similarly, defineA~j,i\(k1,k2\)\{\\tilde\{A\}\}^\{\(k\_\{1\},k\_\{2\}\)\}\_\{j,i\}by appending the row\(1,m,…,mk1−1\)\(1,m,\\ldots,m^\{k\_\{1\}\-1\}\)on top ofAj,i\(k1,k2\)\{A\}^\{\(k\_\{1\},k\_\{2\}\)\}\_\{j,i\}and then retaining the firstℓ\+2\\ell\+2columns\. In this case, an analogous result to Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3)holds for\(k1,k2\)\(k\_\{1\},k\_\{2\}\)defined in \([25](https://arxiv.org/html/2607.05984#A3.E25)\)\.

###### Theorem C\.1\.

Consider the determinant of an\(ℓ\+2\)×\(ℓ\+2\)\(\\ell\+2\)\\times\(\\ell\+2\)minor ofA~j,i\(k1,k2\)\{\\tilde\{A\}\}^\{\(k\_\{1\},k\_\{2\}\)\}\_\{j,i\}that contains the first row and treat it as a polynomial inmm\. Then, the roots of this polynomial aremj​iO,mj​1O​L′,⋯,mj​ℓO​L′m^\{O\}\_\{ji\},m^\{OL^\{\\prime\}\}\_\{j1\},\\cdots,m^\{OL^\{\\prime\}\}\_\{j\\ell\}\.

###### Proof C\.2\.

Similarly to the proof of Theorem 4 inschkoda2024causal, we can show that the determinant of a\(ℓ\+2\)×\(ℓ\+2\)\(\\ell\+2\)\\times\(\\ell\+2\)minor ofA~j,i\(k1,k2\)\\widetilde\{A\}^\{\(k\_\{1\},k\_\{2\}\)\}\_\{j,i\}containing the first row is generically not the zero polynomial inmm\. However, whenmmin the first row is set to any value in\{mj​iO,mj​1O​L′,…,mj​ℓO​L′\}\\\{m^\{O\}\_\{ji\},m^\{OL^\{\\prime\}\}\_\{j1\},\\ldots,m^\{OL^\{\\prime\}\}\_\{j\\ell\}\\\}, the determinant of this minor vanishes\. Since this determinant is a nonzero polynomial of degree at mostℓ\+1\\ell\+1inmm, and theseℓ\+1\\ell\+1values are generically distinct, they are exactly the roots of the polynomial, which completes the proof\.

After obtaining\{mj​iO,mj​1O​L′,⋯,mj​ℓO​L′\}\\\{m^\{O\}\_\{ji\},m^\{OL^\{\\prime\}\}\_\{j1\},\\cdots,m^\{OL^\{\\prime\}\}\_\{j\\ell\}\\\}, we choose the lowest orderkj,ik\_\{j,i\}for which \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) in Proposition[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)is solvable\. Then, thekj,ik\_\{j,i\}\-th and higher\-order cumulants of the residualized variables can be updated, and the procedure can proceed to the next iteration\.

Unlike the proposed method, the modified ReLVLiNGAM still relies on updates based on estimated disturbance cumulants\. Therefore, estimation errors in higher\-order cumulants may still propagate to subsequent iterations\.

## Appendix DIllustration of the Proposed Algorithm on the DAG in Figure[1](https://arxiv.org/html/2607.05984#S2.F1)\(a\)

The LvLiNGAM for the model in Figure[1](https://arxiv.org/html/2607.05984#S2.F1)\(a\) is expressed as

X1\\displaystyle X\_\{1\}=L1\+L2\+e1,\\displaystyle=L\_\{1\}\+L\_\{2\}\+e\_\{1\},X2\\displaystyle X\_\{2\}=λ21​L1\+λ22​L2\+L3\+b21​X1\+e2\\displaystyle=\\lambda\_\{21\}L\_\{1\}\+\\lambda\_\{22\}L\_\{2\}\+L\_\{3\}\+b\_\{21\}X\_\{1\}\+e\_\{2\}=\(b21\+λ21\)​L1\+\(b21\+λ22\)​L2\+L3\+b21​e1\+e2,\\displaystyle=\(b\_\{21\}\+\\lambda\_\{21\}\)L\_\{1\}\+\(b\_\{21\}\+\\lambda\_\{22\}\)L\_\{2\}\+L\_\{3\}\+b\_\{21\}e\_\{1\}\+e\_\{2\},X3\\displaystyle X\_\{3\}=λ33​L3\+b32​X2\+e3\\displaystyle=\\lambda\_\{33\}L\_\{3\}\+b\_\{32\}X\_\{2\}\+e\_\{3\}=b32​\(b21\+λ21\)​L1\+b32​\(b21\+λ22\)​L2\+\(b32\+λ33\)​L3\+b21​b32​e1\+b32​e2\+e3\.\\displaystyle=b\_\{32\}\(b\_\{21\}\+\\lambda\_\{21\}\)L\_\{1\}\+b\_\{32\}\(b\_\{21\}\+\\lambda\_\{22\}\)L\_\{2\}\+\(b\_\{32\}\+\\lambda\_\{33\}\)L\_\{3\}\+b\_\{21\}b\_\{32\}e\_\{1\}\+b\_\{32\}e\_\{2\}\+e\_\{3\}\.
By Propositions[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)–[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4), the proposed method identifies the observed source, the relevant total effects, and the associated disturbance cumulants\. Here,X1X\_\{1\}is the unique observed source, andCh~​\(X1\)=\{X2\}\\widetilde\{\\mathrm\{Ch\}\}\(X\_\{1\}\)=\\\{X\_\{2\}\\\}\. Hence,

𝒃21\\displaystyle\\bm\{b\}\_\{21\}=\{b21,\(b21\+λ21\),\(b21\+λ22\)\},\\displaystyle=\\\{b\_\{21\},\\ \(b\_\{21\}\+\\lambda\_\{21\}\),\\ \(b\_\{21\}\+\\lambda\_\{22\}\)\\\},𝒎31O\\displaystyle\\bm\{m\}^\{O\}\_\{31\}=\{b21​b32,b32​\(b21\+λ21\),b32​\(b21\+λ22\)\},\\displaystyle=\\\{b\_\{21\}b\_\{32\},\\ b\_\{32\}\(b\_\{21\}\+\\lambda\_\{21\}\),\\ b\_\{32\}\(b\_\{21\}\+\\lambda\_\{22\}\)\\\},together with the correspondences between cumulants and total effects

κ\(k\)​\(e1\)\\displaystyle\\kappa^\{\(k\)\}\(e\_\{1\}\)\\:\{b21,b21​b32\},\\displaystyle:\\ \\big\\\{b\_\{21\},\\ b\_\{21\}b\_\{32\}\\big\\\},κ\(k\)​\(L1\)\\displaystyle\\kappa^\{\(k\)\}\(L\_\{1\}\)\\:\{\(b21\+λ21\),b32​\(b21\+λ21\)\},\\displaystyle:\\ \\big\\\{\(b\_\{21\}\+\\lambda\_\{21\}\),\\ b\_\{32\}\(b\_\{21\}\+\\lambda\_\{21\}\)\\big\\\},\(26\)κ\(k\)​\(L2\)\\displaystyle\\kappa^\{\(k\)\}\(L\_\{2\}\)\\:\{\(b21\+λ22\),b32​\(b21\+λ22\)\},\\displaystyle:\\ \\big\\\{\(b\_\{21\}\+\\lambda\_\{22\}\),\\ b\_\{32\}\(b\_\{21\}\+\\lambda\_\{22\}\)\\big\\\},where \([2\.4](https://arxiv.org/html/2607.05984#S2.Ex3)\) is solvable at the orderkk\.

Among the candidates satisfying \([3\.2](https://arxiv.org/html/2607.05984#S3.Ex7)\), we choose

α31=b32​\(b21\+λ21\),β21=b21\+λ21,\\alpha\_\{31\}=b\_\{32\}\(b\_\{21\}\+\\lambda\_\{21\}\),\\qquad\\beta\_\{21\}=b\_\{21\}\+\\lambda\_\{21\},which correspond toκ\(k\)​\(L1\)\\kappa^\{\(k\)\}\(L\_\{1\}\)\. Then

X~2=X2−β21​X1\\displaystyle\\tilde\{X\}\_\{2\}=X\_\{2\}\-\\beta\_\{21\}X\_\{1\}=\(\(λ22−λ21\)​L2−λ21​e1\)\+L3\+e2,\\displaystyle=\\big\(\(\\lambda\_\{22\}\-\\lambda\_\{21\}\)L\_\{2\}\-\\lambda\_\{21\}e\_\{1\}\\big\)\+L\_\{3\}\+e\_\{2\},X~3=X3−α31​X1\\displaystyle\\tilde\{X\}\_\{3\}=X\_\{3\}\-\\alpha\_\{31\}X\_\{1\}=b32​\(\(λ22−λ21\)​L2−λ21​e1\)\+\(b32\+λ33\)​L3\+b32​e2\+e3\\displaystyle=b\_\{32\}\\big\(\(\\lambda\_\{22\}\-\\lambda\_\{21\}\)L\_\{2\}\-\\lambda\_\{21\}e\_\{1\}\\big\)\+\(b\_\{32\}\+\\lambda\_\{33\}\)L\_\{3\}\+b\_\{32\}e\_\{2\}\+e\_\{3\}=b32​X~2\+λ33​L3\+e3\\displaystyle=b\_\{32\}\\tilde\{X\}\_\{2\}\+\\lambda\_\{33\}L\_\{3\}\+e\_\{3\}By Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2), there is only one confounder betweenX~2\\tilde\{X\}\_\{2\}andX~3\\tilde\{X\}\_\{3\}\. Furthermore, by Proposition[2\.3](https://arxiv.org/html/2607.05984#S2.Thmtheorem3), we have

𝒎32\.1O=\{b32,\(b32\+λ33\)\}\.\\bm\{m\}^\{O\}\_\{32\.1\}=\\\{b\_\{32\},\\ \(b\_\{32\}\+\\lambda\_\{33\}\)\\\}\.Choosingγ32\.1=b32\\gamma\_\{32\.1\}=b\_\{32\}givesα31−β21​γ32\.1=0\\alpha\_\{31\}\-\\beta\_\{21\}\\gamma\_\{32\.1\}=0, and Theorem[3\.2](https://arxiv.org/html/2607.05984#S3.Thmtheorem2)concludesX1∉Pa​\(X3\)X\_\{1\}\\notin\\mathrm\{Pa\}\(X\_\{3\}\)\.

In fact, every aligned pair in \([D](https://arxiv.org/html/2607.05984#A4.Ex38)\) can satisfy \([3\.2](https://arxiv.org/html/2607.05984#S3.Ex7)\) and the number of latent confounders between remaining variables is always one\. For instance, suppose that the method selects

α21∗=b21\+λ22,α31∗=b32​\(b21\+λ22\)\\alpha^\{\*\}\_\{21\}=b\_\{21\}\+\\lambda\_\{22\},\\qquad\\alpha^\{\*\}\_\{31\}=b\_\{32\}\(b\_\{21\}\+\\lambda\_\{22\}\)to updateX2X\_\{2\}andX3X\_\{3\}\. The corresponding total effect fromX2X\_\{2\}toX3X\_\{3\}is alsob32b\_\{32\}, and \([3\.2](https://arxiv.org/html/2607.05984#S3.Ex7)\) still holds\. Since there are two unselected aligned groups of total effects, they correspond to two distinct latent variables and hence to two distinct total\-effect columns in the mixing matrix\. Then, the mixing matrices𝑴^\\widehat\{\\bm\{M\}\}and𝑴^latent\\widehat\{\\bm\{M\}\}\_\{\\mathrm\{latent\}\}in Algorithm[1](https://arxiv.org/html/2607.05984#algorithm1)are updated to

𝑴^←\[100\(b21\+λ22\)10b32​\(b21\+λ22\)b321\],𝑴^latent←\[11b21\(b21\+λ21\)b21​b32b32​\(b21\+λ21\)\]\.\\displaystyle\\widehat\{\\bm\{M\}\}\\leftarrow\\left\[\\begin\{array\}\[\]\{ccc\}1&0&0\\\\ \(b\_\{21\}\+\\lambda\_\{22\}\)&1&0\\\\ b\_\{32\}\(b\_\{21\}\+\\lambda\_\{22\}\)&b\_\{32\}&1\\end\{array\}\\right\],\\quad\\widehat\{\\bm\{M\}\}\_\{\\mathrm\{latent\}\}\\leftarrow\\left\[\\begin\{array\}\[\]\{cc\}1&1\\\\ b\_\{21\}&\(b\_\{21\}\+\\lambda\_\{21\}\)\\\\ b\_\{21\}b\_\{32\}&b\_\{32\}\(b\_\{21\}\+\\lambda\_\{21\}\)\\end\{array\}\\right\]\.In this case,X~2\\tilde\{X\}\_\{2\}andX~3\\tilde\{X\}\_\{3\}are

X~2\\displaystyle\\tilde\{X\}\_\{2\}=\(\(λ21−λ22\)​L1−λ22​e1\+e2\)\+L3,\\displaystyle=\\big\(\(\\lambda\_\{21\}\-\\lambda\_\{22\}\)L\_\{1\}\-\\lambda\_\{22\}e\_\{1\}\+e\_\{2\}\\big\)\+L\_\{3\},X~3\\displaystyle\\tilde\{X\}\_\{3\}=b32​\(\(λ21−λ22\)​L1−λ22​e1\+e2\)\+\(b32\+λ33\)​L3\+e3\.\\displaystyle=b\_\{32\}\\big\(\(\\lambda\_\{21\}\-\\lambda\_\{22\}\)L\_\{1\}\-\\lambda\_\{22\}e\_\{1\}\+e\_\{2\}\\big\)\+\(b\_\{32\}\+\\lambda\_\{33\}\)L\_\{3\}\+e\_\{3\}\.
Applying Proposition[2\.2](https://arxiv.org/html/2607.05984#S2.Thmtheorem2)again yieldsX2X\_\{2\}as the current observed source andX3X\_\{3\}as its unique descendant, henceX2∈Pa​\(X3\)X\_\{2\}\\in\\mathrm\{Pa\}\(X\_\{3\}\)\. Applying Proposition[2\.4](https://arxiv.org/html/2607.05984#S2.Thmtheorem4)to\(X2,X3\)\(X\_\{2\},X\_\{3\}\)gives candidatesb32b\_\{32\}and\(b32\+λ33\)\(b\_\{32\}\+\\lambda\_\{33\}\)\. Since the total effect fromX2X\_\{2\}toX3X\_\{3\}is already determined asb32b\_\{32\}, the remaining value\(b32\+λ33\)\(b\_\{32\}\+\\lambda\_\{33\}\)corresponds to a confounder betweenX2X\_\{2\}andX3X\_\{3\}\. Thus,

𝑴^latent←\[110b21\(b21\+λ21\)1b21​b32b32​\(b21\+λ21\)\(b32\+λ33\)\]\.\\displaystyle\\widehat\{\\bm\{M\}\}\_\{\\mathrm\{latent\}\}\\leftarrow\\left\[\\begin\{array\}\[\]\{ccc\}1&1&0\\\\ b\_\{21\}&\(b\_\{21\}\+\\lambda\_\{21\}\)&1\\\\ b\_\{21\}b\_\{32\}&b\_\{32\}\(b\_\{21\}\+\\lambda\_\{21\}\)&\(b\_\{32\}\+\\lambda\_\{33\}\)\\end\{array\}\\right\]\.Consequently,

𝑩^\\displaystyle\\widehat\{\\bm\{B\}\}=𝑰−𝑴^−1=\[000b21\+λ22000b320\],\\displaystyle=\\bm\{I\}\-\\widehat\{\\bm\{M\}\}^\{\-1\}=\\left\[\\begin\{array\}\[\]\{ccc\}0&0&0\\\\ b\_\{21\}\+\\lambda\_\{22\}&0&0\\\\ 0&b\_\{32\}&0\\end\{array\}\\right\],𝚲^\\displaystyle\\widehat\{\\bm\{\\Lambda\}\}=𝑴^−1​𝑴^latent=\[110−λ22λ21−λ22100λ33\]\.\\displaystyle=\\widehat\{\\bm\{M\}\}^\{\-1\}\\widehat\{\\bm\{M\}\}\_\{\\mathrm\{latent\}\}=\\left\[\\begin\{array\}\[\]\{ccc\}1&1&0\\\\ \-\\lambda\_\{22\}&\\lambda\_\{21\}\-\\lambda\_\{22\}&1\\\\ 0&0&\\lambda\_\{33\}\\end\{array\}\\right\]\.Let the columns of𝚲^\\widehat\{\\bm\{\\Lambda\}\}correspond to\(L1,L2,L3\)\(L\_\{1\},L\_\{2\},L\_\{3\}\)\. Then the recovered parent sets are

Pa​\(X1\)=\{L1,L2\},Pa​\(X2\)=\{L1,L2,L3,X1\},Pa​\(X3\)=\{L3,X2\}\.\\displaystyle\\mathrm\{Pa\}\(X\_\{1\}\)=\\\{L\_\{1\},L\_\{2\}\\\},\\quad\\mathrm\{Pa\}\(X\_\{2\}\)=\\\{L\_\{1\},L\_\{2\},L\_\{3\},X\_\{1\}\\\},\\quad\\mathrm\{Pa\}\(X\_\{3\}\)=\\\{L\_\{3\},X\_\{2\}\\\}\.
In finite samples, entries that are theoretically zero in𝑩^\\hat\{\\bm\{B\}\}and𝚲^\\hat\{\\bm\{\\Lambda\}\}may be estimated as nonzero\. For𝑩^\\hat\{\\bm\{B\}\}, whenever no parent\-child relationship is estimated between two observed variables, the corresponding entry of𝑩^\\hat\{\\bm\{B\}\}is set to zero\. For𝚲^\\hat\{\\bm\{\\Lambda\}\}, we prune the matrix by setting entries whose absolute values are below predefined thresholds to zero\.