Learning Sparsest Linear Causal DAGs with Latent Confounders via Higher-Order Cumulants Ming Cai*
[email protected] Graduate School of Informatics, Kyoto University, Kyoto, Japan
Hisayuki Hara
[email protected] Institute for Liberal Arts and Sciences, Kyoto University, Kyoto, Japan
arXiv:2607.05984v1 [cs.LG] 7 Jul 2026
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 finitesample performance compared with existing approaches. Keywords: Causal discovery; DAG; LiNGAM; Latent confounder; Cumulant;
1. Introduction Linear non-Gaussian acyclic models (LiNGAMs) provide a powerful framework for causal discovery (Shimizu et al., 2006, 2011). In the absence of latent variables, LiNGAM enables complete identification of causal DAGs. In many practical applications, however, latent confounders are unavoidable. Hoyer et al. (2008) introduced 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., Eriksson and Koivunen, 2004), 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 (Shimizu and Bollen, 2014). To avoid relying on OICA, several methods have been proposed to estimate canonical LvLiNGAMs via residual independence tests, such as Pairwise LvLiNGAM (Entner and Hoyer, 2011), ParceLiNGAM (Tashiro et al., 2014), Repetitive Causal Discovery (RCD) (Maeda and Shimizu, 2020, 2022; Maeda, 2022), and BANG (Wang and Drton, 2023). However, none of these methods can fully identify ancestral relationships or parent-child relationships between observed variables that form the bow structures (Wang and Drton, 2023). When bow structures are present, Chen et al. (2024) use cumulants of observed variables to identify their ancestral relationships in the bivariate setting with a single latent © M. Cai & H. Hara.
Cai Hara
confounder. Building on this, Chen et al. (2025) extend the approach to multiple latent confounders in the observed bivariate case. Schkoda et al. (2024) proposed 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; see Appendix C). 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.
2. Preliminaries 2.1. Canonical LvLiNGAM Let X = (X1 , . . . , Xp )⊤ be the observed variables and L = (L1 , . . . , Lq )⊤ be the latent confounders. Denote a causal DAG by G = (V , E), where V = X ∪ L and E ⊂ V × V
is the set of directed edges. Define E O = (X × X) ∩ E and E OL = (L × X) ∩ E. We call G O = (X, E O ) the observed DAG of G, and let G OL = (L, X, E OL ) be the latent-toobserved bipartite graph of G. For Xi ∈ X, let Anc(Xi ), Des(Xi ), Pa(Xi ), and Ch(Xi ) denote the sets of ancestors, descendants, parents, and children of Xi , respectively. For two variables V, V ′ ∈ V , we denote a directed edge from V to V ′ by V → V ′ , and let P(V, V ′ ) be the set of directed paths from V to V ′ . For two observed variables Xi and Xj , we define their (possibly latent) confounders as n o Conf(Xi , Xj ) = V : ∃ π ∈ P(V, Xi ), ∃ π ′ ∈ P(V, Xj ) s.t. π ∩ π ′ \ {V } = ∅ . In this paper, we employ the canonical LvLiNGAM (Hoyer et al., 2008), where latent variables are mutually independent, and each has at least two observed children. The model defined by G is expressed as X = ΛL + BX + e,
(1)
where e = (e1 , . . . , ep )⊤ , or, equivalently, h ih i⊤ X = (I − B)−1 Λ, (I − B)−1 L⊤ , e⊤ .
(2)
⊤ The components of u = L⊤ , e⊤ are mutually independent, and each follows a continuous non-Gaussian distribution with nonzero higher-order cumulants. Λ = {λji } and B = {bji } collecththe direct causal coefficients i for Li → Xj ∈ E and Xi → Xj ∈ E, respectively, and
M = (I − B)−1 Λ, (I − B)−1 is the mixing matrix whose entries are total effects. Let O mOL ji and mji denote the total effects of Li and Xi on Xj , respectively. Since latent scales are arbitrary, without loss of generality, we fix λji = 1 for Xj ∈ Ch(Li ) with the highest causal order. The same normalization is also used in Schkoda et al. (2024). Throughout this paper, we assume that the coefficients in (1) and the higher-order cumulants of 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 by P (V ) the joint distribution of V . Suppose that, in (2), there exist Xi and Lj such that Des(Lj ) = Des(Xi ) ∪ {Xi }. Salehkaleybar et al. (2020) showed that, under this condition, exchanging the columns of the mixing matrix M corresponding to Lj and ei yields another valid LvLiNGAM representation while preserving both P (V ) and the ancestral relationships of 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 between Lj and ei . Figure 1 illustrates such an example. DAG (a) is the original causal DAG G, whereas DAG (b) is obtained by a swap between L3 and e2 . Although the two DAGs induce the same P (V ) and ancestral relationships, the swap introduces an additional directed edge X1 → X3 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
Cai Hara
L1
L2
L3
L1
L2
e2
X1
X2
X3
X1
X2
X3
e1
e2
e3
e1
L3
e3
(a)
(b)
Figure 1: An example of observationally equivalent models with different DAGs. genericity assumption, a generic parameterization of DAG (a) corresponds to a measurezero subset of the parameter space of DAG (b). Thus, the sparsest DAG can be regarded as a canonical representative of the observational equivalence class. Schkoda et al. (2024) showed 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 violates this restriction, and therefore ReLVLiNGAM cannot be applied directly. 2.2. Cumulants In this section, we review several results established by Schkoda et al. (2024) concerning the relationship between higher-order cumulants of LvLiNGAM variables and total effects, which also play a central role in this paper. Definition 1 (Cumulants (Brillinger, 2001)) Let I = [p] be the set of indices. For a (k) p-dimensional variable X = (X1 ,. . . ,Xp )⊤ , the k-th order cumulant ci1 ,...,ik is defined by (k)
X
ci1 ,...,ik =
(−1)h−1 (h − 1)!
Y di ∈Di
Di ∈{D1 ,...,Dh }
E
Y
Xj ,
j∈di
where {i1 , . . . , ik } ∈ I k and {D1 , . . . , Dh } is the set of all partitions of {i1 , i2 , . . . , ik }. When i1 = i2 = · · · = ik = i, κ(k) (Xi ) denotes the k-th order cumulant of Xi . From (2), the k-th order cumulants of observed variables can be rewritten as (k) ci1 ,...,ik =
q X h=1
OL (k) mOL i1 h . . . mik h κ (Lh ) +
p X
O (k) mO i1 h . . . mik h κ (eh ).
(3)
h=1
Fix two observed variables Xi and Xj , and treat all variables in V \ {Xi , Xj } as latent. Applying Algorithm A in Hoyer et al. (2008), we obtain a canonical LvLiNGAM in which
the latent confounders Conf(Xi , Xj ) = {L′1 , L′2 , · · · , L′ℓ } are mutually independent. Without loss of generality, assume Xj ∈ / Anc(Xi ). Then, the canonical LvLiNGAM is expressed as Xi =
ℓ X
′
′ mOL ih Lh +vi ,
Xj =
h=1
ℓ X
′
′ O mOL jh Lh +mji vi +vj ,
(4)
h=1 ′
′
OL ′ where vi and vj are disturbances, and mOL ih and mjh are total effects from Lh to Xi and Xj , respectively, in the canonical model over Xi and Xj . ℓ is the number of confounders between Xi and Xj in the original model. Schkoda et al. (2024) estimate the canonical LvLiNGAM under the genericity assump(k ,k ) tion using higher-order cumulants of observed variables. They define the matrix Aj,i1 2 in (5) with k1 < k2 . (k1 ) (k1 ) (k1 ) ci,i,...,j ... ci,j,...,j ci,i,...,i (k1 +1) (k1 +1) 1 +1) c(k ci,i,i,...,j ... ci,i,j,...,j i,i,i,...,i (k1 +1) (k1 +1) c(k1 +1) c . . . c j,i,i,...,i j,i,i,...,j j,i,j,...,j .. .. .. (k1 ,k2 ) .. . Aj,i = (5) . . . . (k2 ) (k2 ) (k2 ) ci,...,i,i,i,...,i,i ci,...,i,i,i,...,i,j . . . ci,...,i,i,j,...,j,j .. .. .. .. . . . . (k2 ) (k2 ) (k2 ) cj,...,j,i,i,...,i,i cj,...,j,i,i,...,i,j . . . cj,...,j,i,j,...,j,j (k ,k )
Define Ai,j1 2 analogously by swapping i and j. Proposition 2 allows us to identify ℓ in the model (4) and the causal order between Xi and Xj . Proposition 2 (Schkoda et al. (2024)) Assume that Xi and Xj are two observed variP 2 −k1 +1 i, k1 ). Then, ables where Xj ∈ / Anc(Xi ). Let d := min( ki=1 (k1 ,k2 ) 1. Aj,i generically has rank min(ℓ + 1, d). (k ,k )
1 2 2. If mO generically has rank min(ℓ + 2,d). ji ̸= 0, Ai,j
(k ,k )
1 2 3. If mO generically has rank min(ℓ + 1,d). ji = 0, Ai,j
According to Schkoda et al. (2024), the smallest possible choice of (k1 , k2 ) is (ℓ + 2, (ℓ + 2) + √ (ℓ) (k ,k ) ⌈(−3 + 8ℓ + 17)/2⌉). Define Aj,i as Aj,i1 2 , where (k1 , k2 ) is this choice. Proposition 2 provides a practical criterion for determining ℓ. If the true number of (ℓ) confounders is ℓ, then Proposition 2 implies that the matrix Aj,i generically has rank ℓ + 1. (r) Consequently, for any r < ℓ, rank Aj,i > r + 1. Therefore, ℓ is the smallest value of r (r) satisfying rank Aj,i = r + 1. Moreover, Items 2 and 3 of Proposition 2 determine whether an ancestral relationship exists between Xi and Xj , and, if so, identify its direction. (ℓ) (ℓ) Let Ãj,i be a matrix obtained by adding (1,m,. . . ,mℓ+1 ) as the first row of Aj,i . Proposition 3 (Schkoda et al. (2024)) Consider the determinant of an (ℓ+2)×(ℓ+2) (ℓ) minor of Ãj,i that contains the first row and treat it as a polynomial in m. Then, the roots OL′ OL′ of this polynomial are mO ji , mj1 , · · · , mjℓ .
Cai Hara
′
′
OL OL Proposition 3 identifies total effects mO ji , mj1 , · · · , mjℓ up to permutation. ⊤ O O 2 O k−1 Proposition 4 (Schkoda et al. (2024)) Define mO . ji := 1, mji , (mji ) , . . . , (mji ) ′
′
OL Similarly, define mOL j1 , . . . , mjℓ . Then, the system of equations i⊤ O h (k) OL (k) ′ (k) ′ mji , mOL , . . . , m κ (v ), κ (L ), . . . , κ (L ) i j1 1 jℓ ℓ h i⊤ (k) (k) (k) = ci,i,...,i , ci,i,...,j , . . . , ci,j,...,j ,
(6)
is generically uniquely solvable if k ≥ ℓ + 1. When (6) is uniquely solvable, each root returned by Proposition 3 is associated with a unique k-th order cumulant of a disturbance or latent variable, (k) OL (k) ′ OL (k) ′ mO ji , κ (vi ) , mji , κ (L1 ) , . . . , mjℓ , κ (Lℓ ) . 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 orders k for which (6) is solvable.
3. Proposed Method We propose a method to recover the sparsest DAG under the genericity assumption. In the following, let G denote the sparsest DAG in the observational equivalence class. The proposed method first applies Proposition 2 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 3, together with their correspondences to the k-th order cumulants of the disturbances and latent confounders established by Proposition 4, 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. Based on the pairwise ancestral relationships estimated using Proposition 2, variables without observed ancestors are treated as observed sources. Let Xs be one such source. f s ) ⊂ Ch(Xs ) be the set of observed variables whose observed ancestors consist Let Ch(X only of observed sources and include Xs . To identify the remaining children of Xs , we recursively examine its descendants while maintaining two sets: the closed set Xclosed , containing variables already identified as children of Xs , and the open set Xopen , containing descendants whose parent-child relationships with Xs have yet to be determined. Initially, f s ), Xclosed = Ch(X n o Xopen = Xj ∈ Des(Xs ) \ Xclosed : ∀Xi ∈ Xclosed , Anc(Xj ) ∩ Des(Xi ) ⊆ Xclosed .
(7)
By definition of Xopen , no observed variable in X \ Xclosed lies on a directed path between Xclosed and Xopen . Lemma 5 For any Xj ∈ Xopen , bjs can be written as X bis mO bjs = mO ji . js −
(8)
f s) i:Xi ∈Ch(X
According to Lemma 5, the parent-child relationship between Xs and Xj can be determined by testing whether bjs = 0. However, Proposition 3 identifies the total effects appearing in (8) only up to permutation, so bjs cannot be computed directly. Denote by I and J the index sets of Xclosed and Xopen , respectively. For j ∈ J and O i ∈ I, let mO js and bis denote the sets of candidate values of mjs and bis , respectively, O returned by Proposition 3. Initially, bis = mis . As discussed in the previous section, Proposition 4 associates each candidate total effect with the corresponding k-th order cumulant of a disturbance or latent confounder. For any i ∈ I and j ∈ J , the systems (6) for the pairs (Xs , Xi ) and (Xs , Xj ) always share the k-th order cumulant of es . Under the genericity assumption, distinct disturbances and latent confounders have distinct k-th order cumulants. Therefore, we can choose αjs ∈ mO js and βis ∈ bis so that they correspond to the same k-th order cumulant. Such a choice is not necessarily unique, since the two systems may also share the k-th order cumulants of latent confounders common to the pairs (Xs , Xi ) and (Xs , Xj ). Define the residualized variables X̃i = Xi − βis Xs ,
X̃j = Xj − αjs Xs .
By Lemma 10 in Appendix A, X̃i and X̃j can be regarded as observed variables in an induced canonical LvLiNGAM. Applying Propositions 3 and 4 to this induced model yields a candidate set for the total effect from X̃i to X̃j . We denote this set by mO ji.s . Theorem 6 provides a criterion for determining whether Xs ∈ Pa(Xj ) for each j ∈ J . Theorem 6 Under the genericity assumption, Xs ∈ / Pa(Xj ) if and only if there exist Y Y αjs ∈ mO βI,s ∈ bis , γj,I.s ∈ mO js , ji.s s.t.
i∈I ⊤ αjs − βI,s γj,I.s = 0,
i∈I
(9)
where αjs and βis , i ∈ I are chosen so that their associated k-th order cumulants are equal. Once the parent-child relationship between Xs and Xj is determined using Theorem 6, we f s ) as follows: update Xclosed , Xopen , and Ch(X ( f s) ← Xclosed ← Xclosed ∪ {Xj }, Xopen ← Xopen \ {Xj }, Ch(X
f s ), Ch(X f s ) ∪ {Xj }, Ch(X
Xs ̸∈ Pa(Xj ), Xs ∈ Pa(Xj ).
When Xopen becomes empty, it is reinitialized using (7). As shown in the proof of Lemma 5 f s ). The coefficients of the directed in Appendix B, Lemma 5 still holds after updating Ch(X
Cai Hara
f s ) can also be computed by Lemma 5. Therefore, the edges from Xs to the variables in Ch(X parent-child relationships of Xs and variables in Xopen can still be identified by Theorem 6. After determining all parent-child relationships between Xs and its descendants, we residualize each descendant Xj ∈ Des(Xs ) with respect to Xs before proceeding to the next iteration: X̃j = Xj − αjs Xs .
(10)
For each j ∈ J , choose αjs from triples (αjs , βI,s , γj,I,s ) satisfying (9) whenever available; and otherwise, choose it randomly from mO js , so that all selected values of αjs correspond to the same k-th order cumulant. Let R = [p] \ {s}. For each j ∈ R, let αjs be defined as above when Xj ∈ Des(Xs ), and set αjs = 0 when Xj ∈ / Des(Xs ). Let XR = (Xj )⊤ j∈R ,
αR,s = (αjs )⊤ j∈R .
Then, the residualized variables X̃R after removing Xs are given by X̃R = (X̃j )⊤ j∈R := XR − αR,s Xs .
(11)
By Lemma 10 and Corollary 11 in Appendix A, we obtain the following theorem. Theorem 7 Under the genericity assumption, X̃R admits an LvLiNGAM representation whose observed DAG is the induced subgraph of G O on X \ {Xs }. By Theorem 7, after the update (10), X̃R admits an LvLiNGAM representation whose observed DAG coincides with the induced subgraph of G O on X \{Xs }. Although the above selection procedure may yield different vectors αR,s corresponding to different k-th order cumulants, all such choices induce the same observed DAG. Let As denote the set of all vectors αR,s obtained by selecting different common k-th order cumulants in the above procedure. Although every vector in As induces the same observed DAG over the remaining variables, different choices may yield different latentto-observed structures. Since G is assumed to be the sparsest DAG in its observational equivalence class, the remaining task is to identify the vector in As that yields the sparsest latent-to-observed structure. Specifically, for each αR,s ∈ As , define X L(αR,s ) = ℓij (αR,s ), i,j∈R
where ℓij (αR,s ) denotes the number of latent confounders between the residualized variables X̃i and X̃j identified by Proposition 2. Since every vector in As induces the same observed DAG and the same number of latent variables, minimizing L(αR,s ) is equivalent to selecting the sparsest latent-to-observed structure, and hence the sparsest DAG. Accordingly, we select any minimizer α∗R,s ∈ arg min L(αR,s ). (12) αR,s ∈As
Theorem 8 Under the genericity assumption, the sparsest latent-to-observed structure over X̃R obtained from α∗R,s coincides with the structure obtained from the induced subgraph of G OL on V \ {Xs } by absorbing every latent variable having only one observed child into the disturbance of that child.
Thus, the update (10) becomes ∗ X̃j = Xj − αjs Xs ,
∀Xj ∈ X \ {Xs },
(13)
∗ ) ∗ where α∗R,s = (αjs j∈R is selected by (12). If multiple such choices of α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 7 and 8, updating the variables using α∗R,s yields an LvLiNGAM whose DAG coincides with the induced subgraph of G induced by removing Xs , 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 G. This result is formalized in the following theorem.
Theorem 9 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 V . According to the selected combination of total effects in the update (13), the proposed method obtains the mixing matrix among the observed variables, namely (I − B)−1 , corresponding to the sparsest DAG over X. Hence, B can be estimated from (I − B)−1 . In finite samples, B can be pruned by enforcing consistency with the estimated parent-child relationships. Moreover, the remaining total-effect columns corresponding to latent sources form (I − B)−1 Λ, which ensures the sparsest latent-to-observed bipartite graph between L and X. Multiplying (I − B) by (I − B)−1 Λ yields an estimate of Λ, and thus identifies the directed edges from L to X. In finite samples, Λ can be pruned by setting entries whose absolute values are below a predefined threshold to zero. The update (13) 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–3 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 order
Cai Hara
Algorithm 1: Proposed Method Input : Observed data matrix X ∈ Rn×p Output : Estimated causal graph Gb c ← Ip×p , M clatent ← ∅ bO ← (X, ∅), M 1 Initialization G 2 Identify ancestral relationships among X and observed sources Xs by Proposition 2 3 while Xs ̸= ∅ do 4 Xs,next ← ∅ 5 foreach Xs ∈ Xs do f s ), Xopen , and Xclosed 6 Identify Ch(X 7 Compute all possible total effects from Xs and the corresponding disturbance cumulants by Propositions 3 and 4 8 while Xopen ̸= ∅ do f s ) to Xopen after regressing out Xs under each possible total effect 9 Compute total effects from Ch(X of Xs by Proposition 3 10 Determine the parent–child relationships between Xs and the nodes in Xopen by Theorem 6 f s ), Xopen , and Xclosed 11 Update Ch(X 12 end f s ) into GbO 13 Add all identified edges Xs → Xj for Xj ∈ Ch(X 14 Remove the effects of Xs from Des(Xs ) using the update (13) with the selected total effects c with all selected total effects in (13) 15 Replace the corresponding entries in M clatent 16 Append other possible total effects not used in Equation (13) to M 17 Identify newly emerging observed sources according to the ancestral relationships by Proposition 2 and add them to Xs,next 18 end 19 Xs ← Xs,next 20 end clatent , and construct GbOL from Λ b =I −M c−1 , Λ b =M c−1 M b 21 B 22
return Gb = GbO ∪ GbOL
k for which (6) is uniquely solvable. Specifically, the value of k can be determined by increasing it from k = 2 and choosing the smallest one for which the corresponding linear system in (6) 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. We also provide a detailed example illustrating the proposed algorithm on the models in Figure 1 in Appendix D.
4. Simulations L1
L2
L3
L1
L2
L3
L1
L2
L3
X1
X2
X3
X1
X2
X3
X1
X2
X3
(a) case I
(b) case II
(c) case III
Figure 2: Three models used for simulations.
X4
This section reports simulation results1 on the three causal DAGs in Figure 2, which are not handled by the original ReLVLiNGAM. We compare the proposed method with ReLVLiNGAM (Schkoda et al., 2024) 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. When k < ℓ + 1, the system (6) 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 (6). For completeness, Appendix C describes a simple modification of ReLVLiNGAM that replaces low-order cumulants with sufficiently high-order ones when solving (6). We also include this modified version in the experimental comparison. 4.1. Settings All disturbances and latent variables are sampled from 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 from Uniform(0.9, 1.1). For each latent variable Li , the coefficient from Li to its observed child with the highest causal order is fixed to one. All other coefficients in Λ and B are independently drawn from 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 size N is set to 1K, 10K, 100K, and 1M, and each experiment is repeated 50 times. We evaluate the performance of the methods using the following metrics: • Ntol and Nobs : the number of runs in which the DAGs of G and G O are correctly recovered, respectively (see Figure 3); • PRE, REC, and F1: the average precision, recall, and F1-score of the estimated edges of G (see Figure 4). In the proposed method, we consider (9) in Theorem 6 to hold in finite-sample settings ⊤ γ if |αjs − βI,s j,I.s | < τ0 ∈ {0.2, 0.15, 0.125, 0.1} for increasing sample sizes. We further prune an edge from L to X if its absolute coefficient is below τ0,L = 0.3. For the other methods, we also apply τ0 and τ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 and τ0,L , and select the sparsest DAG. Although both the proposed method and ReLVLiNGAM employ Proposition 2, they (k ,k ) differ in how the rank of Aj,i1 2 is determined. Let σr be the r-th largest singular value of (k ,k )
Aj,i1 2 and let τs and τcs be two predefined thresholds. ReLVLiNGAM treats P P σr as zero if σr /σ1 ≤ τs . In contrast, the proposed method sets σr = 0 if 1 − i∈[r] σi / i∈[d] σi ≤ τcs , (k ,k )
where d is the number of singular values of Aj,i1 2 and is defined as in Proposition 2. Following Schkoda et al. (2024), we also impose an upper bound ℓhighest on the number 1. Code is available at https://anonymous.4open.science/r/Test_20260701
Cai Hara
Case I 43 39
Number
30
29
0
6
29
25
29
30
0
100K
Sample Size
0
1M
Ntol of Proposed Nobs of Proposed
27 24 2221
2323
20 10
10K
35
13 2
0
50
40
14
5
00
1K
35
19
17 4
36 26
23
19
20 10
33
46 43
Case III
50
0
6
26 23
5
3
0
1K
Ntol of Proposed(True) Nobs of Proposed(True)
7
2
10K
25
19 14 15
16 9
19 13 6
0
4
0
100K
Sample Size
40
3535 30
Number
46
40
Case II
48 44
Number
50
0
1M
Ntol of ReLVLiNGAM(Revised) Nobs of ReLVLiNGAM(Revised)
3536
30
25
20 10 0
16 12
13 77
3
7
27
27 24
23 19
16
2
6
6
00
10
1K
10K
100K
00
Sample Size
10
5 0
1M
Ntol of ReLVLiNGAM Nobs of ReLVLiNGAM
Figure 3: The performances in Ntol and Nobs . of latent variables, with ℓhighest = 2 in cases I and III, and ℓhighest = 4 in case II. We set τs = 0.008(i − 1)/N 0.125 for the original ReLVLiNGAM and τcs = 0.002 + 0.0005(i − 1) for our method, where i denotes 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 increase k from k = 2 until the system in (6) is uniquely solvable, and then use the resulting value of k in Proposition 4. 4.2. Discussion Figure 3 shows that the proposed method outperforms the original ReLVLiNGAM on all evaluation metrics across all settings. Since the DAGs in Figure 2 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 show that, for the proposed method, both the mean (triangles) and the median (horizontal bars) of these metrics increase as N increases, indicating increasingly accurate recovery of parent–child relationships over 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 L → 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
REC of Case I
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4 0.2
1K
10K
100K
PRE of Case II
0.0
1M
0.4 0.2
1K
10K
100K
REC of Case II
0.0
1M
1.0
1.0
0.8
0.8
0.8
0.6
0.6
0.6
0.4 0.2 0.0
Score
1.0
Score
0.4 0.2
1K
10K
100K
PRE of Case III
0.0
1M
1K
10K
100K
REC of Case III
0.0
1M
1.0
0.8
0.8
0.6
0.6
0.6
Score
1.0
0.8
0.2 0.0
0.4 0.2
1K
10K
100K
1M
Proposed
0.0
10K
100K
1M
10K
100K
1M
10K
100K
1M
F1 of Case II
0.4
1.0
0.4
1K
0.2
Score
Score
0.0
Score
1.0
0.2
Score
F1 of Case I
1.0
Score
Score
PRE of Case I 1.0
1K
F1 of Case III
0.4 0.2
1K
Proposed(True)
10K
100K
1M
ReLVLiNGAM(Revised)
0.0
1K
ReLVLiNGAM
Figure 4: The performances of PRE, REC, and F1-score. to subsequent iterations, the effect of cumulant estimation errors is expected to be less severe.
5. Real Data We further evaluate the practical usefulness of the proposed method by applying it, together with ParceLiNGAM (Tashiro et al., 2014), RCD (Maeda and Shimizu, 2020), and the original and modified versions of ReLVLiNGAM (Schkoda et al., 2024), to the Sachs protein dataset from Sachs et al. (2005). The Sachs dataset (Sachs et al., 2005) 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(a). Although the DAG in Figure 5(a) is not a canonical LvLiNGAM, it can be transformed into its canonical form in Figure 5(b) without changing the observed structure by applying Algorithm A of Hoyer et al. (2008). The thresholds τs and τcs for the proposed method and ReLVLiNGAMs are set in the same manner as in the numerical experiments in Section 4. The upper bound ℓhighest is set
Cai Hara
PKC
PKA
PKC
PKA Raf
Raf
Mek
Erk
Akt
Raf
(a) Reference DAG L1
L2
Raf
Mek
Mek
Erk
L4
L5
Erk
Akt
(e) ReLL
L1
Raf
Erk
Akt
Raf
Mek
Erk
Akt
Akt
(b) Canonical model L3
Mek
(c) ParceLiNGAM
L2
Mek
L3
Erk
L1
Akt
(f) modified ReLL
Raf
(d) RCD L2
Mek
L3
Erk
Akt
(g) Proposed method
Figure 5: The results of different methods applied to the Sachs dataset.
to 2 for both methods. Here, instead of employing a hard threshold, we verify (9) in Theorem 6 using 99% bootstrap confidence intervals constructed from 400 bootstrap resamples. ParceLiNGAM and RCD use the Hilbert–Schmidt independence criterion (HSIC; Gretton et al., 2007) 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 L → X structure, we prune latent-to-observed edges whose absolute coefficients are below 0.1. Figures 5(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(a) with only one extra edge, Mek → Akt, although it estimates a denser L → X structure. ParceLiNGAM fails to identify the existence of latent confounders, and outputs redundant edges Raf → Erk, Raf → Akt, Mek → Akt. RCD correctly concludes that there is no edge between Raf and Erk, but does not identify any directed edges among the observed variables. Even with the pruning, the original ReLVLiNGAM yields an incorrect causal order over X and a dense L → X structure. The modified ReLVLiNGAM yields an L → 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. 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 L → X structure and most closely match the structure from latent variables to observed variables in Figure 5(a).
6. Conclusion 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 higherorder 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.
Acknowledgments This work was supported by JST SPRING under Grant Number JPMJSP2110 and JSPS KAKENHI under Grant Numbers 25K15017.
References David Ross Brillinger. Time series: data analysis and theory. Society for Industrial and Applied Mathematics, 2001. Wei Chen, Zhiyi Huang, Ruichu Cai, Zhifeng Hao, and Kun Zhang. Identification of causal structure with latent variables based on higher order cumulants. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38(18), pages 20353–20361, 2024. Wei Chen, Linjun Peng, Zhiyi Huang, Haoyue Dai, Zhifeng Hao, Ruichu Cai, and Kun Zhang. Identification of causal direction under an arbitrary number of latent confounders, 2025. Mathias Drton, Rina Foygel, and Seth Sullivant. Global identifiability of linear structural equation models. The Annals of Statistics, 39(2), 2011. Doris Entner and Patrik O. Hoyer. Discovering unconfounded causal relationships using linear non-Gaussian models. In Takashi Onada, Daisuke Bekki, and Elin McCready, editors, New Frontiers in Artificial Intelligence, pages 181–195, Berlin, Heidelberg, 2011. Springer Berlin Heidelberg. ISBN 978-3-642-25655-4. Jan Eriksson and Visa Koivunen. Identifiability, separability, and uniqueness of linear ICA models. IEEE signal processing letters, 11(7):601–604, 2004.
Cai Hara
Arthur Gretton, Kenji Fukumizu, Choon Teo, Le Song, Bernhard Schölkopf, and Alex Smola. A kernel statistical test of independence. In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems, volume 20, pages 1–8. Curran Associates, Inc., 2007. Patrik O. Hoyer, Shohei Shimizu, Antti J. Kerminen, and Markus Palviainen. Estimation of causal effects using linear non-Gaussian causal models with hidden variables. International Journal of Approximate Reasoning, 49(2):362–378, 2008. ISSN 0888-613X. Special Section on Probabilistic Rough Sets and Special Section on PGM’06. Takashi Nicholas Maeda. I-RCD: an improved algorithm of repetitive causal discovery from data with latent confounders. Behaviormetrika, 49(2):329–341, 2022. Takashi Nicholas Maeda and Shohei Shimizu. RCD: Repetitive causal discovery of linear non-Gaussian acyclic models with latent confounders. In International Conference on Artificial Intelligence and Statistics, pages 735–745. PMLR, 2020. Takashi Nicholas Maeda and Shohei Shimizu. Repetitive causal discovery of linear nonGaussian acyclic models in the presence of latent confounders. International Journal of Data Science and Analytics, pages 1–13, 2022. Karen Sachs, Omar Perez, Dana Pe’er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(5721):523–529, 2005. Saber Salehkaleybar, AmirEmad Ghassami, Negar Kiyavash, and Kun Zhang. Learning linear non-Gaussian causal models in the presence of latent variables. Journal of Machine Learning Research, 21(39):1–24, 2020. Daniela Schkoda, Elina Robeva, and Mathias Drton. Causal discovery of linear nonGaussian causal models with unobserved confounding. arXiv preprint arXiv:2408.04907, 2024. Shohei Shimizu and Kenneth Bollen. Bayesian estimation of causal direction in acyclic structural equation models with individual-specific confounder variables and non-gaussian distributions. Journal of Machine Learning Research, 15(76):2629–2652, 2014. Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. A linear non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7:2003–2030, 2006. Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvärinen, Yoshinobu Kawahara, Takashi Washio, Patrik O. Hoyer, and Kenneth Bollen. Directlingam: A direct method for learning a linear non-Gaussian structural equation model. Journal of Machine Learning Research, 12:1225–1248, 2011. Tatsuya Tashiro, Shohei Shimizu, Aapo Hyvärinen, and Takashi Washio. ParceLiNGAM: A causal ordering method robust against latent confounders. Neural Computation, 26(1): 57–83, 2014.
Y. Samuel Wang and Mathias Drton. Causal discovery with unobserved confounding and non-gaussian data. Journal of Machine Learning Research, 24(271):1–61, 2023.
Appendix A. Preservation of the LvLiNGAM Structure under Source Removal Assume that Xs is an observed source of G. Let R = [p] \ {s} and fix us ∈ {es } ∪ Pa(Xs ). As discussed in Section 2.2, the corresponding DAG might change when swapping us and es . Let G (us ) denote the DAG in the observational equivalence class of G that is obtained by (u ) swapping es and us . For each j ∈ R, choose a root αjs s returned by Proposition 3 whose associated k-th order cumulant is κ(k) (us ). Then the update rule (11) is expressed as (u ) (u ) (u ) (u ) . X̃R s = XR − αR,ss Xs , αR,ss = αjs s j∈R
(u )
Lemma 10 Under the genericity assumption, X̃R s admits an LvLiNGAM representation whose DAG is the induced subgraph of G (us ) obtained by deleting Xs and all edges linked to Xs . Proof Denote the LvLiNGAM representation associated with G (us ) by X = Λ(us ) L(us ) + B (us ) X + e(us ) . (u )
(14)
(u )
s αjs s is the total effect from us to Xj . Let BR,R denote the submatrix of B (us ) with rows
(u )
s and columns indexed by R. Let bR,s denote the column of BR,R corresponding to Xs . By the definition of total effects,
(u )
(u )
(u )
(u )
s αR,ss = bR,ss + BR,R αR,ss .
(15)
(u )
Let ΛR,:s denote the submatrix of Λ(us ) consisting of the rows indexed by R and let eR denote the subvector of e(us ) consisting of the entries indexed by R. Restricting the structural equations to R, we have (u )
(u )
(u )
(u )
s XR = bR,ss Xs + BR,R XR + ΛR,:s L(us ) + eR s .
(u )
(u )
Since XR = X̃R s + αR,ss Xs , we obtain (u )
(u )
(u )
(u )
(u )
(u )
(u )
(u )
s X̃R s = −αR,ss Xs + bR,ss Xs + BR,R (X̃R s + αR,ss Xs ) + ΛR,:s L(us ) + eR s (us ) (u ) (u ) (us ) (us ) (u ) (u ) (u ) = BR,R X̃R s + bR,ss + BR,R αR,s − αR,ss Xs + ΛR,:s L(us ) + eR s .
By (15), the coefficient of Xs is zero. Thus, (u )
(u )
(u )
(u )
(u )
s X̃R s = BR,R X̃R s + ΛR,:s L(us ) + eR s .
(16)
Cai Hara
This is an LvLiNGAM representation over XR , where its coefficient matrix among observed (us ) (u ) variables is exactly BR,R , and its coefficient matrix from L is ΛR,:s . Hence, the corresponding DAG is obtained from G (us ) by deleting Xs and all edges linked to Xs , namely the induced subgraph of G (us ) after deleting Xs . Based on Lemma 10, we have the following Corollary. Corollary 11 Let BR,R be the coefficient matrix among X \ {Xs } in G (es ) . Then, for each us ∈ Pa(Xs ) ∪ {es }, (u )
s BR,R = BR,R .
That is, the observed parent-child relationships among XR do not depend on us . Proof The mixing matrix for X is expressed as
−1
(I − B)
=
1 0 −1 mO (I − B R,R ) R,s
,
where mO R,s is the vector of total effects from Xs to X \ {Xs }. For any us ∈ Pa(Xs ) ∪ {es }, the corresponding swap affects only the total effects associated with Xs . Thus, (I −B (us ) )−1 (u ) is obtained from (I − B)−1 by replacing the column corresponding to Xs with (1, αR,ss )⊤ " (I − B (us ) )−1 =
1 (u )
αR,ss
# 0 , (I − BR,R )−1
where this replacement is trivial when us = es . Computing the inverse gives "
# 1 0 . I − B (us ) = (u ) −(I − BR,R )αR,ss I − BR,R Thus, we can obtain B (us ) " B
(us )
=
0
0 (u )
(I − BR,R )αR,ss
BR,R
# ,
(u )
s implying BR,R = BR,R from the lower-right block.
Remark 12 The updated model (16) may contain latent variables in Pa(Xs ) having only one observed child in R. Such latent variables can be absorbed into the corresponding disturbance term, as in the canonicalization described by Hoyer et al. (2008).
Appendix B. Proofs of Theorems in Section 3 Proof of Lemma 5 Following Drton et al. (2011), the total effect from an observed source Xs to Xj is written as X Y bhl . mO = b + js js π∈P(Xs ,Xj )\{(Xs ,Xj )} l,h:(Xl ,Xh )∈π
f s ) is a descendant of Xs , and every node in Xopen has ancestors Since every variable in Ch(X f s ) or not influenced by Xs , every directed path that are either contained in {Xs } ∪ Ch(X from Xs to Xj other than the direct edge Xs → Xj must first pass through some node f s ). Thus, Xi ∈ Ch(X X X Y X b bis mO mO = b + b = b + is js js hl js ji , f s) i:Xi ∈Ch(X
π∈P(Xi ,Xj ) l,h:(Xl ,Xh )∈π
f s) i:Xi ∈Ch(X
which establishes (8). Proof of Theorem 6 By Lemma 5, the true direct effect satisfies X bjs = mO bis mO js − ji.s .
■
i∈I O Suppose first that Xs ∈ / Pa(Xj ). Then bjs = 0. Since the true values mO js , bis , and mji.s are O contained in the candidate sets mO js , bis , and mji.s , respectively, there exists a choice
αjs ∈ mO js ,
γji.s ∈ mO ji.s ,
βis ∈ bis ,
such that αjs and βis , i ∈ I, correspond to the same k-th order cumulant and X ⊤ αjs − βI,s γj,I.s = mO − bis mO js ji.s = 0. i∈I
Thus, (9) holds. Conversely, suppose that Xs ∈ Pa(Xj ). Then bjs ̸= 0. Under the genericity assumption, O no choice of candidates from mO js , bis , and mji.s , with αjs and βis corresponding to the same k-th order cumulant, can satisfy ⊤ αjs − βI,s γj,I.s = 0.
Hence, (9) fails generically. Therefore, (9) holds if and only if Xs ∈ / Pa(Xj ). ■ Proof of Theorem 7 Immediate from Lemma 10 and Corollary 11. ■ Proof of Theorem 8 Let κ(k) (us ) denote the cumulant corresponding to α∗R,s . By (14), the corresponding latent-to-observed coefficient matrix is Λ(us ) . Since α∗R,s is chosen to minimize the number
Cai Hara
of latent parents, the support of Λ(us ) is the sparsest among all latent-to-observed structures in the observational equivalence class. Therefore, it coincides with the structure of G OL . (u ) By (16), the latent-to-observed coefficient matrix for the model of X̃R is ΛR,:s . Hence, its support coincides with that of the induced subgraph of G OL on V \ {Xs } obtained by absorbing every latent variable having only one observed child into the disturbance of that child. ■ Proof of Theorem 9 Immediate from Theorem 7 and Theorem 8. ■
Appendix C. A Modification of the Original ReLVLiNGAM As mentioned earlier, the origianl ReLVLiNGAM, proposed by Schkoda et al. (2024), 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, let X1 be an observed source and let Xi ∈ Des(X1 ). The update rule for Xi in ReLVLiNGAM is given by X Xi ← Xi − mO e − mOL (17) 1 i1 ih Lh . h:Lh ∈Pa(X1 )
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 X (k) (k) O (k) OL (k) ci1 ,...,ik ← ci1 ,...,ik − mO · · · m κ (e ) − mOL 1 i1 1 ik 1 i1 h · · · mik h κ (Lh ). h:Lh ∈Pa(X1 )
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 node X1 , Schkoda et al. (2024) proposed estimating (κ(k) (e1 ), κ(k) (L1 ), . . . , κ(k) (Lq )) by solving the linear system (k) (k) c1···11 1 1 ··· 1 κ (e1 ) (k) mO mOL · · · mOL κ(k) (L1 ) c1···12 21 2q 21 . .. = . .. .. .. .. . .. . . . . OL · · · mOL (k) mO κ(k) (Lq ) pq p1 mp1 c
1···1p
However, when the local restriction is violated, this linear system becomes underdetermined and therefore cannot be solved. By contrast, Proposition 4 enables the estimation of the higher-order cumulants of the disturbance and the latent confounders even without the local restriction. In Schkoda et al. (2024), Proposition 4 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, (k ,k ) it first constructs Aj,i1 2 using sufficiently high-order cumulants and applies Proposition 2 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 3 to the present setting. Next, Proposition 4 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 choose k1 and k2 to be the smallest values that satisfy the requirements of Proposition 2. For each descendant Xi of an observed source X1 , let ki,1 denote the smallest order for which the linear system in Proposition 4 is generically solvable. ki,1 can be estimated by increasing k from 2 and checking whether the corresponding linear system (6) is generically solvable. Applying the same procedure to every descendant of X1 yields the corresponding values ki,1 . Now suppose that Xi is the next observed source after removing X1 . For another descendant Xj of X1 , let √ −3 + 8ℓ + 17 k1 = max{ℓ + 2, ki,1 , kj,1 }, k2 = k1 + , (18) 2 (k ,k )
and construct Ai,j1 2 as in (5). Starting from ℓ = 0, ℓ is estimated by increasing ℓ until the rank condition in Proposition 2 is satisfied. (k ,k ) (k ,k ) Similarly, define Ãj,i1 2 by appending the row (1, m, . . . , mk1 −1 ) on top of Aj,i1 2 and then retaining the first ℓ + 2 columns. In this case, an analogous result to Proposition 3 holds for (k1 , k2 ) defined in (18). (k ,k )
Theorem 13 Consider the determinant of an (ℓ + 2) × (ℓ + 2) minor of Ãj,i1 2 that contains the first row and treat it as a polynomial in m. Then, the roots of this polynomial OL′ OL′ are mO ji , mj1 , · · · , mjℓ . Proof Similarly to the proof of Theorem 4 in Schkoda et al. (2024), we can show that the e(k1 ,k2 ) containing the first row is generically determinant of a (ℓ + 2) × (ℓ + 2) minor of A j,i not the zero polynomial in m. However, when m in the first row is set to any value in OL′ OL′ {mO ji , mj1 , . . . , mjℓ }, the determinant of this minor vanishes. Since this determinant is a nonzero polynomial of degree at most ℓ + 1 in m, and these ℓ + 1 values are generically distinct, they are exactly the roots of the polynomial, which completes the proof. ′
′
OL OL After obtaining {mO ji , mj1 , · · · , mjℓ }, we choose the lowest order kj,i for which (6) in Proposition 4 is solvable. Then, the kj,i -th and higher-order cumulants of the residualized variables can be updated, and the procedure can proceed to the next iteration.
Cai Hara
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 D. Illustration of the Proposed Algorithm on the DAG in Figure 1(a) The LvLiNGAM for the model in Figure 1(a) is expressed as X1 = L1 + L2 + e1 , X2 = λ21 L1 + λ22 L2 + L3 + b21 X1 + e2 = (b21 + λ21 )L1 + (b21 + λ22 )L2 + L3 + b21 e1 + e2 , X3 = λ33 L3 + b32 X2 + e3 = b32 (b21 + λ21 )L1 + b32 (b21 + λ22 )L2 + (b32 + λ33 )L3 + b21 b32 e1 + b32 e2 + e3 . By Propositions 2–4, the proposed method identifies the observed source, the relevant total effects, and the associated disturbance cumulants. Here, X1 is the unique observed f 1 ) = {X2 }. Hence, source, and Ch(X b21 = {b21 , (b21 + λ21 ), (b21 + λ22 )}, mO 31 = {b21 b32 , b32 (b21 + λ21 ), b32 (b21 + λ22 )}, together with the correspondences between cumulants and total effects κ(k) (e1 ) :
b21 , b21 b32 ,
(k)
(L1 ) :
(b21 + λ21 ), b32 (b21 + λ21 ) ,
κ(k) (L2 ) :
(b21 + λ22 ), b32 (b21 + λ22 ) ,
κ
(19)
where (6) is solvable at the order k. Among the candidates satisfying (9), we choose α31 = b32 (b21 + λ21 ),
β21 = b21 + λ21 ,
which correspond to κ(k) (L1 ). Then X̃2 = X2 − β21 X1 = (λ22 − λ21 )L2 − λ21 e1 + L3 + e2 , X̃3 = X3 − α31 X1 = b32 (λ22 − λ21 )L2 − λ21 e1 + (b32 + λ33 )L3 + b32 e2 + e3 = b32 X̃2 + λ33 L3 + e3 By Proposition 2, there is only one confounder between X̃2 and X̃3 . Furthermore, by Proposition 3, we have mO 32.1 = {b32 , (b32 + λ33 )}. Choosing γ32.1 = b32 gives α31 − β21 γ32.1 = 0, and Theorem 6 concludes X1 ∈ / Pa(X3 ).
In fact, every aligned pair in (19) can satisfy (9) 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 )
to update X2 and X3 . The corresponding total effect from X2 to X3 is also b32 , and (9) 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 c and M clatent in Algorithm 1 are updated to matrix. Then, the mixing matrices M 1 0 0 1 1 c ← (b21 + λ22 ) clatent ← b21 1 0 , M (b21 + λ21 ) . M b32 (b21 + λ22 ) b32 1 b21 b32 b32 (b21 + λ21 ) In this case, X̃2 and X̃3 are X̃2 = (λ21 − λ22 )L1 − λ22 e1 + e2 + L3 , X̃3 = b32 (λ21 − λ22 )L1 − λ22 e1 + e2 + (b32 + λ33 )L3 + e3 . Applying Proposition 2 again yields X2 as the current observed source and X3 as its unique descendant, hence X2 ∈ Pa(X3 ). Applying Proposition 4 to (X2 , X3 ) gives candidates b32 and (b32 + λ33 ). Since the total effect from X2 to X3 is already determined as b32 , the remaining value (b32 + λ33 ) corresponds to a confounder between X2 and X3 . Thus, 1 1 0 clatent ← b21 . (b21 + λ21 ) 1 M b21 b32 b32 (b21 + λ21 ) (b32 + λ33 ) Consequently,
0 0 0 b =I −M c −1 = b21 + λ22 0 0 , B 0 b32 0 1 1 0 b =M c −1 M clatent = −λ22 λ21 − λ22 1 . Λ 0 0 λ33 b correspond to (L1 , L2 , L3 ). Then the recovered parent sets are Let the columns of Λ Pa(X1 ) = {L1 , L2 },
Pa(X2 ) = {L1 , L2 , L3 , X1 },
Pa(X3 ) = {L3 , X2 }.
In finite samples, entries that are theoretically zero in B̂ and Λ̂ may be estimated as nonzero. For B̂, whenever no parent-child relationship is estimated between two observed variables, the corresponding entry of B̂ is set to zero. For Λ̂, we prune the matrix by setting entries whose absolute values are below predefined thresholds to zero.