Differentiable Causal Discovery of Linear Non-Gaussian Acyclic Models Under Unmeasured Confounding
Authors: Yoshimitsu Morinishi
Shohei Shimizu
Venue: Accepted by TMLR ()

Abstract: We propose a score-based method that extends the framework of the linear non- Gaussian acyclic model (LiNGAM) to address the problem of causal structure estimation in the presence of unmeasured variables. Building on the method pro- posed by Bhattacharya et al. (2021), we develop a method called ABIC LiNGAM, which assumes that error terms follow a multivariate generalized normal distribu- tion and employs continuous optimization techniques to recover acyclic directed mixed graphs (ADMGs). We demonstrate that the proposed method can esti- mate causal structures, including the possibility of identifying their orientations, rather than only Markov equivalence classes, under the assumption that the data are linear and follow a multivariate generalized normal distribution. Additionally, we provide proofs of the identifiability of the parameters in ADMGs when the er- ror terms follow a multivariate generalized normal distribution. The effectiveness of the proposed method is validated through simulations and experiments using real-world data.

Published in Transactions Machine Learning Research (8/2025)
Differentiable Causal Discovery of Linear Non-Gaussian Acyclic Models Under Unmeasured Confounding
s7024104 shigg u.o... ac.jj Graduate School of Data Science, Shiga University; Strategy and Consulting, Data and AI, Accenture Co.
Yoshimitsu Morinishi
Ltd.
Shohei Shimizu
SANKEN; University of Osaka; Shiga Unéversity
Reviewed on OpenReview: https: // openreveew. net/forunt id- 4056
Abstract
We propose score-based method that extends the framework of the lincar non- Gaussian acyclic model (LiNGAM) to address the problem of causal structure estimation in the presence of unmeasured variables. Building on the method pro- posed by Bhattacharya et al. (2021), we develop method called ABIC LINGAM, which assumes that error terms follow a multivariate generalized normal distribu- tion and employs continuous optimization ion techniques to recover acyelic directed mixed graphs (ADMGs). We demonstrate that the proposed method can esti- mate causal structures, including the possibility of identifying their orientations, rather than only Markov equivalence classes, under the assumption that the data are linear and follow multivariate generalized normal distribution. Additionally, provide proofs of the identifiability of the parameters in ADMGs when the er- ror terms follow multivariate generalized normal distribution. The effectiveness of the proposed method is validated through simulations and experiments using weal-world data.
1 Introduction
1.1 Background and Motivation
Uncovering causal relationships from observational data- referred to as causal discovery is critical objective in numerous disciplines, including life sciences, social sciences, and marketing. While randomized controlled trials (RCTs) are the gold standard for identifying causal effects, they are frequently infeasible in practice due to prohibitive costs or ethical constraints. Therefore, it is essential to develop robust methods for inferring causal structures from observational data.
Many cxisting algorithms assume the absence of unmeasured variables, implying that all relevant factors are fully observed (Spirtes et al., 2000; Chickering, 2002; Shimizu et al., 2006). For example, the linear non- non-Gaussian acyclic model (LiNGAM) (Shimizu al., 2006) demonstrates that causal directions can be uniquely identified when errors are non-Gaussian and no hidden confounders

---

Published in Transactions Machine Learning Research (8/2025)
exist. However, this assumption of no unmeasured variables is often unrealistic. In practical applications such consumer purchasing behavior analysis, latent factors like income or occupation can significantly confound observed relationships.
To address this issue, acyclic directed mixed graphs (ADMGs) (Richardson & Spirtes, 2002) have been introduced. By including both directed and bidirected edges, ADMGs can represent the more intricate dependencies induced by latent variables, capturing confounding effects that traditional directed acyclic graphs (DAGs) cannot. Additionally, ADMGs support nonparametric equality constraints, such as the Verma constraints (Verma & Pearl, 1990), providing flexible framework for causal discovery in the presence of hidden confounders.
Key Point (1): ADMGs allow us to capture both direct causal effects and latent confounding within single graph, which is crucial in many real-world scenarios where some variables are unobserved.
1.2 Related Work
Methods for estimating ADMG structures can be broadly classified into three categories:
Constraint-based methods. These use repeated conditional independence tests to build partial ancestral graphs (PAGs) (Spirtes et al., 2000). Prominent examples include FCT (Spirtes et al., 2000), RFCI (Colombo et al., 2012), and GFCI (Ogarrio et al., 2016), extended to hign-dimensional data in approaches like IFCI (Chen et al., 2021). While effective under sufficient sample sizes, these methods can be sensitive to the multiplicity and order of tests.
Score-based methods. By optimizing a global objective (e.g., log-likelihood or BIC), score- based approaches can provide a more coherent framework for model selection (Nowzohour et al., 2017; Bernstein et al., 2020; Chen et al., 2021; Claassen & Bucur, 2022; Ng et al., 2024). However, many existing score-based techniques assume Gaussian error terms, thus limiting their ability to recover strict causal orientations when latent variables induce equivalence classes.
Hybrid methods.
These combine constraint-based and score-based strategies, aiming to balance local statistical decisions with global model fit. GFCI (Ogarrio et al., 2016) exemplifies this line of research.
Recent developments emphasize strict orientation recovery under latent confounding by leveraging non-Gaussianity or nonlinearity (Shimizu et al., 2006; Tashiro et al., 2014; Wang Drton, 2024). Notably, the BANG method (Wang & Drton, 2024) exploits higher-order moments a constraint- based procedure to identify bow-free ADMGs beyond Markov equivalence classes. Further refining these ideas, Wang & Drton (2024) (sce also references thercin) incorporate non-Gaussian assump- tions into framework that can recover directions more strictly, though specialized moment condi- tions are still crucial for direction recovery.
In parallol, Bhattacharya et al. (2021) introduced ABIC, a continuous score-based method for learn- ing ancestral, arid, or bow-free ADMGs under the linear Gaussian setting. Although ABIC unilies different structural constraints within a differentiable optimization framework, strict orientation recovery remains clusive under Gaussianity. Subsequent work, such as Ng ct al. (2024), has fo-
2

---

Published in Transactions Machine Learning Research (8/2025)
cused on scalable continuous relaxations for large ADMGs, but again, the challenge of uniquely determining causal directions persists without non-Gaussian assumptions.
Key Point (2): In Gaussian settings, we can typically only recover up to a Markov equivalence class. Non-Gaussian error assumptions are needed to fully determine directions.
1.3 Our Goal and Contributions
Motivated by these developments, propose ABIC LiNGAM, novel approach for strictly identifying bow-free ADMGs under linear non-Gaussian error distributions. Specifically, we assume error terms follow a multivariate generalized normal distribution (MGGD), thereby generalizing ABIC (Bhattacharya et al., 2021) beyond Gaussianity. Our primary contributions are:
1. Orientation Recovery Beyond Gaussian Equívalence. We show that non- Gaussianity (i.e., * in MGGD) enable more precise identification of causal di- rections in bow-free ADMGs, going beyond the limitations of Gaussian-based methods that recover only up to Markov equivalence classes.
2. Unified Framework. When the shape parameter 3 = 1 (Gaussian case), our method reduces to the original ABIC scheme, ensuring compatibility with purely Gaussian models while offering broader scope for non-Gaussian scenarios.
3. Identifiability Proof. Building on Brito & Pearl (2002) and Wang & Drton (2024), we rigorously prove identifiability of bow-free ADMGs under MGGD assumptions, solidifying the theoretical foundation of our approach.
4. Empirical Validation and Scalability. We conduct extensive simulations (up to 50 variables) and apply our method to real-world social survey data, illustrating higher accu- racy than constraint-based (c.g., FCI) and non- Con-Gaussian methods (c.g.. BANG). Morcover, the fully differentiable formulation scales effectively to larger graphs.
Key Point (3): In summary, we extend a continuous score-based approach (ABIC) from Gaussian to non- Non-Ganssian error assumptions, potentially enabling more thorough orientation recovery for bow-free ADMGs.
the first continuous score-based method that achieves full identification of bow-free ADMGs under linear non-Gaussian assumptions. The remainder of this paper is organized ADMGs, and the MGGD. Section 3 introduces our differentiable constraints and likelihood-based objective. Section 4 presents the ABIC LiNGAM algorithm and theoretical properties. Sections 5- 6 report experimental results on synthetic and real datasets. We conclude in Section 7 with a discussion future directions and broader implications. Unlike a simple relaxation of assumptions, the key is the introduction of non- Gaussianily which provides crucial higher-order information. Hence, our framework can strictly identify causal directions in linear models with latent confounders, rather than settling for Markov equivalence classes.
To our knowledge, ABIC LINGAM
follows: Section 2 reviews the linear structural equation model, bow-frec
3

---

Published in Transactions Machine Learning Research (8/2025)
2 Problem Definition
2.1 Representation by Linear SEM
In this section, review linear SEMs and their graphical representations. We use uppercase letters (e.g., .X) to denote variables or nodes in the graph and indexed uppercase letters (c.g., X;) to denote specific variables or nodes. We also use the following standard matrix notation: Ag denotes the element the ith row and jth column of matrix A, A-i,- denotes the submatrix of A obtained by removing the ith row and jth column, and A.; denotes the ith column of A.
Additionally, for each vertex belonging set V, let (pa(t) | i E V} and {sp(i) | E V} be two families of index sets. The vertex set of G is the index set V, and G contains the edge j i if and only if j E Pa(i) and the edge ++ i if and only if j E Sp(i) (or equivalently, i Sp(j)). Furthermore, {sp(i) holds if and only if E sp(j). These two families of sets {pa(i) V and (sp(1) | VV define the system of structural equations.
}} satisfies the following symmetry condition: for any j E EVV, E sp(i)
2.1.1 Linear SEM
We consider linear SEMs for d variables, parameterized by weight matrix 0 E Rdddd For each variable X E X, the structural equation is
X,
Σ jEpa(s)
Oo + ci. ieV
Here, the noise terms are mutually independent. In this case, sp(1) = for all i, since no unmeasured variables exist. The graph G and corresponding binary adjacency matrix D E (0, 1]ddd are defined as follows: An edge X; > X; exists G if and only 0f + 0, which case Di, = The graph is acyclic if and only if f can be created as an upper triangular matrix by the permutation of vertex labeling.
Key Point (4): Even in linear SEM, once unneasured variables are discovered, we may need an ADMG (not just DAG) to handle latent confounding via bidirected edges.
2.1.2 Linear SEM with unmeasured variables
An observed set of variables causally insufficient if there exist unmeasured variables that are the ancestors of two more observed variables in the system. In linear structural equation model (SEM), these unmeasured variables often manifest as dependencies among the error terms Pearl (2009). Consider d-dimensional random vector X = (XI...,,A, represented by real-valued matrices 6, f. E dd For each X,, the structural equation takes the following form:
x;= jEPa())
Σ X + iEV.
(2)
Here, c c (a......ca) is a vector of error terms with zero mean without loss generality and not necessarily Gaussian. Allowing non-Gaussian noise terms accommodates a wider class of underlying distributions and may improve identifiability via higher-order moments nonsymmetric distributional features Shimizu et al. (2006); Wang & Drton (2024).
4

---

Published in Transactions Machine Learning Research (8/2025)
In the special case where a given variable X, has no unmeasured variables, its error term e, may be independent of all the others. However, if unmeasured variables inffuence multiple observed variables, their corresponding error terms become dependent on each other. These dependencies are captured by the matrix Ω = EBeec"] ", which does not need to be diagonal. For the Gaussian noise, the marginalized distribution of X is a zero-mean multivariate normal with a covariance matrix as follows:
Σ= (1- ())
And the same covariance structure be considered for non-Gaussian errors, least at the level of second moments. In non-Gaussian setting, higher-order moments and distributional asymmetries can be exploited to identify causal directions and latent structures.
The induced graph G is an ADMG that includes both directed (-) and bidirected (++) edges. The graph G and associated binary adjacency matrices E {0, 1] 1]axd and B E {0, are defined as follows: a directed edge X;- + X; exists in G if and only if d.j # O, in which case Di 1. A bidirected edge Xj ++ X; exists in G if and only if G4 # 0 (symmetry ensures Ωz + 0), in which case Bij B,s = 1. In the special case where there are unmeasured variables, the ADMG reduces to DAG, and the B matrix is zero matrix.
In summary, this framework does not restrict the noise terms to be Gaussian, allowing a broader class of SEMs that can represent latent variable-induced dependencies through non-Gaussian dis- tributions. By leveraging non-Gaussinnity, one can potentially achieve stronger identifiability and more robust causal inferences than would be possible under Gaussian assumptions alone.
Key Point (5): Gon-Gaussian error terms often break Markov equivalences that would hold under purely Gaussian assumptions, thus allowing clearer causal direction recovery.
2.2 Motivation Example
As discussed in Section 2.1.1, when there are no unmeasured variables, the observed variables can be represented as a DAG. Thus, the problem reduces to estimating the structure of DAG, However, when unmeasured variables are present, a DAG cannot adequately represent the rela- tionships between variables while accounting for such unmeasured variables. Therefore, we use an ADMG, which can represent latent variuble-induced covariation and confounding through directed and bidirected edges. Consequently, in the presence of unmeasured variables, the problem reduces to estimating the structure of an ADMG. This section builds on the work of Bhattacharya et al. (2021).
Figure I(a) depicts a DAG, which represents the relationships among the variables in the absence of unmeasured variables. Figures (b), (c), and (d) illustrate examples of ADMGs that we aim to estimate in this study. Figure (b) shows ancestral ADMG, where directed path X, + and bidirected edges X; < Xj exist simultaneously in G for any pair of vertices X,, X, E X. Figure (c) shows an arid ADMG that does not contain any c-trees. A c-tree is subgraph of G where its directed edges form a directed tree, and its bidirected edges form a single bidirectional connected component within the subgraph. For details on c-trees, see Shpitser & Pearl (2006). Figure (d) shows a bow-free ADMGs, where no pair of vertices X; + X, and X; <> X, both exist in G. These three types of ADMGs exhibit an inclusion relationship, with bow-free ADMGs being the most general type of ADMGs.
Xj
5

---

Published in Transactions Machine Learning Research (8/2025)
B
B
A
A
D
D
c
D
(a)
(b)
(c)
Figure 1: (a) DAG without unmeasured variables. (b) Ancestral ADMGs. (c)
Arid ADMGs. (d)
bow-free ADMGs.
Ancestral Arid Bow-free
In Bhattacharya et al. (2021), these ADMGs are expressed as differentiable constraints, allowing the selection of the appropriate ADMGs type to be estimated based on the data. This study adopts the differentiable constraints proposed by Bhattacharya al. (2021), enabling the selection of suitable ADMGs types according to the data.
2.3 Identifiability in the Model
This section discusses the identifiability of the parameters in bow-free ADMGs where the error terms follow a multivariate generalized normal distribution. Brito & Pearl (2002) proved that given a bow- free ADMG model, the parameters are almost everywhere identifiable from the observed covariance matrix. Since this fact is often utilized under the assumption that the error terms are Gaussian, we show this study that also applics bow-free ADMG models when the crror terms follow a multivariate generalized normal distribution. Furthermore, we draw on Wang & Drton (2024), who demonstrated that the model identify causal directions, and not only Markov equivalence classes, using the non-Gaussia of error terms. This study also provides evidence that causal directions can be estimated.
2.3.1 Definition and Key Terms
bow-free ADMG. A bow-free ADMG is type of directed mixed graph (ADMG) that may include latent variables. For any two nodes X, and X2, it disallows having both
x and X, X
simultaneously (i.e,, no 'bow" shape). This ensures that any bidirected edge (representing covari- ance structures induced by latent variables) does not overlap with a directed edge between the same pair of nodes.
Markov equivalence class. Markov equivalence class of DAGs (or ADMGs) is a of graphs that encode the same set of conditional independence relationships. When the error terms
6

---

Published in Transactions Machine Learning Research (8/2025)
assumed be Gaussian, one often can only identify the Markov equivalence class, meaning that the directions of certain edges cannot be distinguished.
Non- Non-Gaussianity. Con-Gaussinnity refers to error terms that do not follow a Gaussian distri- bution (e.g., distributions with skewness or high kurtosis). Such non-Gaussian characteristics can allow one to distinguish causal directions that remain indistinguishable under purely Gaussian assumptions, by exploiting higher-order moments or skewness.
2.3.2
On the Identifiability of Parameters in bow-Free ADMGs with Multivariate Generalized Normal Distributions
Brito & Pearl (2002) demonstrated that a bow-free ADMG model almost always identifiable from the observed covariance matrix. As the argument in Brito & Pearl (2002) primarily assumes that the error terms are Gaussian, we follow Brito & Pearl (2002) to demonstrate that when the error terms have multivariate generalized normal distribution (MGGD), bow-free ADMG model is almost always identifiable from the observed covariance matrix. (See the Appendix for the proof.)
Theorem
Let be a bow-free ADMG with error terms following a multivariate generalized normal distri- bution, and let the set of parameters of be {8.2}. Then, for almost all 0, the following holds:
E()) = I(0')
implics
0 = 6'.
In other words, if two parameter sets 9 and 0' yicld the same covariance matrix E, then A and #' must be identical, except possibly when 0 belongs to set of Lebesgue measure zero.
This result follows from extending the argument of Brito & Pearl (2002) (which relies Wright s path analysis and the inductive construction of linear SEMs) to the MGGD sctting. Because the MGGD includes the Gaussian distribution (as special case with 3 = 1) and is closed under linear transformations ions, imposing the structural constraints of a bow-free ADMG yields a unique parameter set {6, 2t..
Moreover, if the error terms are non-Gaussian (i.e., 3 + 1). then not only are the parameters identifiable, but there is also potential to determine the directions of causal edges, as suggested by Wang & Drton (2024) and others. In this paper, we leverage this non Non-Gaussianity to propose method for estimating both the structure and orientations in bow-free ADMGs.
Intuition. mntuution.Because Because the generalized normal distribution includes the Gaussian as a special case, each valid (6, ()) pair generally induces a unique covariance structure, making parameters identifiable in the bow-free setting.
2.3.3 Identifiability of the structure in bow-Free ADMGs Using Non- Non-Gaussianity
Wang Drton (2024) demonstrated that when the data are non-Gaussian and correspond to a bow-free ADMG, the bow-free ADMG can be consistently recovered, including both the Markov cquivalence class and causal directions. The multivariate generalized normal distribution assumed
7

---

Published in Transactions Machine Learning Research (8/2025)
for the observed data in this study is non- Gaussan except in the special case. Therefore, un- like Bhattacharya et al. (2021), who were limited to recovering the Markov oquivalence class, it is expected that the causal directions can also be estimated.
3 Decomposition of the Log-Likelihood Function of the Multivariate Generalized Normal Distribution
3.1 Probability Density Function of the Multivariate Generalized Normal Distribution
As defined by Gómez et al. (1998), the probability density function of the multivariate generalized Gaussian distribution (MGGD) given by:
T(2) π 2
B
f(X n1 8, =
exp
-2(((X - - (x-
(4)
2
|>|
where X is a p-dimensional random vector (p 2 1) that follows a power-exponentiil distribution with parameters μ, Σ, and B. Specifically, j E IR, Σ is a (pxp) positive-defiitte symmetric matrix, and E (0, r()) denotes the gamma. function. Notably, the MGGD reduces multivariate normal distribution when 1. In this distribution, if Σ diagonal matrix, then the correlation coefficients wetween the components become zero. However, because the multivariate generalized normal distribution belongs to the elliptical family, a correlation does not imply independence. Nonetheless, if the components are assumed to be generated independently, they can be considered truly independent rather than merely uncorrelated.
Gómez et al. (1998) show that the MGGD is invariant under affine transformations. More precisely, if f(X 1x,Exxx 3) is the probability density function, then for the affine transformation
==CX² CX+b,
where C' is nonsingular matrix, a vector in R2, and the transformed variable Y follows J(Y Cux + b, Cxx CT, B). This indicates that the family of distributions remains within the same class under any nonsingular linear transformation and translation.
This affine transformation property is particularly useful in linear structural equation models (cs- pecially (2)). In such models, the relationships among the observed variables and the covariance structure of the error terms are modeled, leading to the covariance matrix of the observed vari- ables in the form
Σ == -111)))10(110))
Furthermore, discussed in Section 3.2.1, one estimate the parameters and Ω from the observed covariance matrix. In other words, if the data follow an MGGD, then in principle and can be estimated from an observed covariance matrix of the form Σ = (I - 2 5)-
Key Point (6): The elliptical nature and affine-invariance of the MGGD are crucial for applying linear SEMs. They ensure that under linear transformations (e.g. confounding), the distributional form remains within the same family, aiding our later decomposition.
8

---

Published in Transactions Machine Learning Research (8/2025)
3.2 Log-Likelihood Function
Assuming that in the ADMG graph G (V, of the linear model (2), N observations are drawn, where all the components of y zero (the mcan vector is zero). The reason for setting all the components of to zero is to prevent notational clutter without loss of generality. In this case, the ligliike function is given by (4) as follows:
(11 B|x)) N log r
2
P
+ N log
3
P N log
P
P
N
N log r
N log 2
log |
23
28
(6)
N
Σ E - 8
(1)
2
3.3 Decomposition of the Log-Likelihood Function
The main component of the proposed algorithm is the decomposition of the log-likelihood function of the multivariate generalized normal distribution, inspired by Dorton et al. (2009), who decomposed the log-likelihood function of the multivariate normal distribution.
Let X; R denote the ith of the observation matrix X and X-: Xv/s be the (V\(i) x N submatrix of X. We adopt the abbreviated notation Xo to represent the C x N submatrix of the D x N matrix X, where C < D.
Theorem 3
Let i X be a variable node in the ADMG graph G (V, of the linear model (2). Let ||>|| Tr x and define Qu.-s as the conditional variance of E, given E-; as follows:
au. = Ω2 - -0,
(7)
Here, D4 is the row vector obtained by removing the ith element from the ith row, 2_ a, is the column vector obtained by removing the ith clement from the ith column, and 2- denotes the submatrix obtained by removing the ith row and ith column from Ω. Additionally, let 2 (0-1,-4)-2 Then, the log-likelihood function ((B, Q, 3) of the graph G (V, E)) can be decomposed
N log Pú.-i- log det (ot
N
(1,,,11X))
1 (xx!") Si.pa(1) x()) x(I)
pa(i) Qe (0_)
sp()))
(1) s 121 _001 -
(8)
Proof. By rearranging the log- likelloodd function described in equation (6), considernng the constant parts unrelated to the cocfficient matrices and Ω noting that the determinant det(7 - 6) 1

---

Published in Transactions Machine Learning Research (8/2025)
since 8 is acyelie, and that (1 6) e, and considering that the covariance matrix can be
expressed using the adjacency matrix B and P as E Var(X) = (I - 1 5) T obtain the following:
N
(p,
N
log |حح
1
x
(("TT-1X1)
1=1
N
log [I - 6)0((I - 5) 1
1 N
Σ
X (I- (X 8) (- ))
(9)
=
2
2
1=1
N
log
|S2|
1 N
((0)
100TT--0 (20))
2
2
l=1
We can partition Q as block matrix:
Uu Q-i Qii
D4,-i
Ω
(10)
Based on equations (7) and (10), we can rearrange log [[]] as shown in equation (11):
log [2]] = log (Sz - -4 0-1,)) + log 12-
(11)
log Qúi-i + log (
Next, we rearrange the term e14T in Dm (2 We partition Q-1 as a block matrix:
P -2 2-
Ω --
Qui P
Pi--
(12)
Ω --1
2 0-)--
-i 11-1522,--
Considering that 12-1 is block matrix, rearrange e as follows:
(U) (1T00-1(00 = =(e!") ( (- (-2:) -- 0 0
-
Q -_, 0
= - -s (5, It -2.- 2y (U)
(1)
"0 (1)
(13)
From this rearrangement, the log- ligliikllhood function becomes equation (14):
N
(14, E,8 = 2
log Qú.-i
log (et
2-1 1 ((e - -
(T UT
(2)) (14)
+
1=1
10

---

Published in Transactions Machine Learning Research (8/2025)
By definition, the error term & x(1) 81,pa(1) x-1) pa(ts) bow-free ADMGs, have Qr - ! sisp(()) (0 (=) This yields the claimed
Morcover, since we are dealing with
decomposition.
The decomposition of the log- lilkiiiihoodd function is based on decomposing the joint distribution of into the marginal distribution of e-; and conditional distribution (c, E=_)- In particular, shown in (14), the squared term (c.") -,, represents the deviation e from its conditional expectation given E-4, which plays a key role in deriving the likelihood decomposition. This idea leads to an approach similar to that of Dorton et al. (2009), who decomposed the lig-likelihood function of the multivariate normal distribution and proposed an iterative algorithm. The steps of this algorithm distribution. To fix the marginal distribution of e--,, must fix the submatrix of f excluding the ith row and ith column and the submatrix of 8_1,v excluding the ith row. This is because e- is determined depending on Q_ and 8_iv-
based fixing the marginal distribution of E-; and estimating the conditional
(Intuition.) By isolating each node's contribution and conditioning out the rest, the log-likelihood can be decomposed into local pieces that depend on the parents and siblings (in terms of bidirected edges). Such an approach is crucial for efficiently optimizing over ADMG structures under the MGGD.
Before presenting the full decomposition, we first define the pseudo-variable Z-; as follows:
z_, 9 -_,_,
(15)
where Q-ii-i is the submatrix of obtained by removing the i-th row and column, and E-i the error vector excluding the i-th component. We will use Z_, in the subsequent derivations to simplify the notation of the conditional distribution.
Note that - E-r in (equation 15) captures the same quantity that appears in the decomposition of E; given E-i. In later equations, if we write Ω Q4 1 e-i, it should be understood that the vector part [_] E-v is precisely Z_:.
N
log Qz. -
1
N
Qu.
(x(0)
-
6 -bipa(i)
(10) Qisp(() 1 (2_1,,
2
2
sp(()
1=1
N
log Pii--i
x
Σ
x,
Σ Qi
20
Hi,-i
1=1
jEpa(i)
kEep(1)
(16)
Assuming 8 2 1, we rearrange equation (16) using Hölder's inequality (for details on the application of Hölder's incquality, see the APPENDIX).
11

---

Published in Transactions Machine Learning Research (8/2025)
N Σ
log Q4
x
Σ
Σ
23 2(0) 2
2
2n
11.
1=1
jEPa(S)
kEep(i)
21
N
log Aii-i
N
x
Σ
Σ
20 (3.
N8-1
jEpa(i)
kEsp(()
2
N
log Sti
1
x
Σ
Σ
li 2,11
2n
N8-1
1i.-
jEpa(1)
kEsp(i)
29
N log Pii-i
1
X;
Σ (E
8
E kEep(()
Di,Ks
20 23.
N8-1
(17)
Maximizing equation (17) to estimate 8 and Ω is equivalent performing regression of X; (as the target variable) on X, (the parent variable of X;) and the pseudo-variale variables Ze, considering the shape parameter B.
Utilizing these observations, in the next section, we propose a method for causal structure estimation using continuous optimization, considering the presence of unmeasured variables and assuming that the error variables follow multivariate generalized normal distribution, using the decomposition results of the log-likelihood function organized in this section to estimate 8, Ω.
4 Proposed Method
4.1 Causal Discovery Based on Differentiable Scores
Score-based methods aim to estimate causal structures by maximizing a graph's score (e.g., log- likelihood) given the data. Learning DAGs from data is an NP-hard problem because it is chal- lenging to efficiently enforce combinatorial acyclicity constraints ((Chickoring, 1996)).
Zheng et al. (2018) proposed a new approach for score-based DAG learning by converting the traditional combinatorial optimization problem (18) into a continuous optimization problem (19).
min F(0) dee
subject to G(D) E DAGs
(18)
min F(0) subject to h(0) 0,
(19)
850
where G(O) is a d-node graph induced by the weight matrix 0 1adxd and F : Rddd - R score function. h Rdxd + R is smooth function over real matrices, and the constraint h(0) 0
12

---

Published in Transactions Machine Learning Research (8/2025)
can precisely characterize the acyclicity of the graph. Causal structure estimation via continuous optimization climinates the need for specialized algorithms to explore the combinatorial space of DAGs and instead allows the use of standard numerical algorithms for constrained problems, making implementation particularly straightforward, as mentioned in Zheng et al. (2018).
The acyclicity constraint is defined as follows:
h(0) trace (c -
(20)
Here, 0 o denotes the Hadamard product (element-wise multiplication). trace (c"o0) is the trace (sum of the diagonal elements) of the matrix exponential, and d is the number of variables. This constraint ensures that the matrix forms a DAG.
Zheng al. (2018) used the augmented Lagrangian method as continuous optimization technique. This method solves constrained optimization problems using objective function that includes penalty terms and is formulated as
min L(Ø) + ||ll|| + ah(O) + h())² subject to h(0) 0.
(21)
Here, À is the weight of the regularization term, is the Lagrange multiplier, and the penalty coefficient.
Key Point (7): Score-based methods combined with differentiable constraints allow us to leverage standard optimization routines (e.g. gradient-based) instead combinatorial searches, thus scaling more casily to larger graphs.
4.2 ABIC
Although the differentiable score-based causal discovery mothod using continuous optimization pro- posed by Zheng et al. (2018) has been successful in estimating causal structures DAGs, it cannot be directly applied to ADMGs. This because ADMGs require two adjacency matrices, and B, represent the directed and bidirected edges, respectively. To extend the differentiable algebraic characterization to ADMGs, Bhattacharya et al. (2021) proposed differentiable constraints for each causal graph. In Bhattacharya et. al. (2021), three differentiable constraints (i.e., ancestral, arid, and bow-free)are proposed for ADMGs as shown in Table 1.
ADMGs
Algebraic Constraint
Ancestral
trace(eP) d sum(e" BB
Arid Bow-free
(race(e d Greenery (D, D, B) (race(e d + sum D B) 0
Table 1: Differentiable constraints for each causal graph in ADMGs (Bhattacharya el al., 2021).
For example, estimate the causal structure of bow-free ADMGs, the constraint equation becomes (22).
13

---

Published in Transactions Machine Learning Research (8/2025)
h(0) = -r(ep)) d sum(D B)
(22)
Here, tr(A) denotes the trace of the square matrix A, i.c., the sum of diagonal entries.
Here, sum(-) denotes the sum of all the elements in matrix. It has been proven that when (()) 0, the estimated graph corresponds to an ADMG ((Bhattacharya et al., 2021)). Essentially, tr(eP) -d signifies the standard acyclicity constraint for directed edges, and the latter term sum(Do ensures that bidirected edges are not introduced when directed edge exists (i.e., it enforces the bow-free ADMGs property).
Bhattacharya et al. (2021) used the augmented Lagrangian method, similar to Zheng et al. (2018), to convert the problem into an optimization problem with quadratic penalty term, and proposed ABIC, which solves the following primal equation at each iteration:
min ABICA(X;))) + (1(0))²++ + ah((),
(23)
where p the weight of the penalty term, and a is the Lagrange multiplier. Then, the Lagrange multiplier is updated as + a ph(0). Intuitively, by optimizing the primal equation with a large P. force h(0) to be very close to zero, thus satisfying the equality constraint.
(Intuition.) By formulating ADMG constraints as smooth penalty terms, we can directly incor- porate them into standard optimization frameworks, rather than enumerating all possible graphs. This approach generalizes well beyond DAGs to bow-free or ancestral ADMGs.
4.3 ABIC LiNGAM
In this study, we extend the method proposed by Bhattacharya et al. (2021) and present algorithm2 that can estimate causal structures when the error terms follow a multivariate generalized normal distribution. The basic framework is the same as that in Bhattacharya et al. (2021), but we consider that the error terms follow multivariate generalized normal distribution and incorporate the shape parameter into the loss function. The shape parameter must be estimated from the observed data in advance. We use the estimated shape parameter B. In addition, since the multivariate normal distribution corresponds to = the proposed method can handle the normal distribution case, thereby generalizing the method of Bhattacharya et al. (2021). Depending on the data, it is possible to switch between non- Gaussianity and Gaussianity.
Key Point (8): By gencralizing from ABIC (Gaussian) to ABIC LiNGAM (non-Gaussian), we can exploit higher-order statistios to potentially identify the orientations of edges in ADMGs. This is especially powerful when hidden confounders introduce bidirected edges.
Comparison of ABIC vs. ABIC LINGAM:
1. Residual norm:
ABIC uses || 11 (Ganssian assumption).
ABIC LINGAM uses || ||20 for non-Gaussian errors if 1.
14

---

Published in Transactions Machine Learning Research (8/2025)
Algorithm 1: Algorithm 1: ABIC
Input X E Rexd Observed data
Ω E Raxd Initial error covariance
tol > 0: Convergence threshold
maa iters: Max iterations
h(()= Differentiable constraint (e.g. bow-free)
P. a Augmented Lagrangian coefficients
A: Regularization weight
Output: 8,, Ω Final estimates
Initialization:
1. 80, .2" + Ia) or random.
2. c < In(n) (used in lanh(c(0)).
3. Define
LS(0) 2n
Σ -
Xo., z S2.
(Here, Z(0) are pscudo-voriables updated per iteration.) for I 0 to mar ilers - 1 do
(A) Pseudo-variables: for == d do
- X,, - +
- 0, - + 5-. (0(,,,))))
(B) Parameter optimization: Solve
(8 15+1, + arg min LS(()) 2 2
+ah(() + A Σ tanh (c10,1) }.
for i = to d do
6 - X - X 81+1
2 var(
(C) Convergence check:
if | 5'1 + || 211 < tol then
break
return
8¹, Ω
Tips:
M0) enforce DAG ADMG constraint (e.g. now-free).
Increase p gradually ((0) stays large.
If 80 0, f20 Ta. residuals start simply.
2. Shape parameter 3:
ABIC LINGAM requires 8 >
If 3 1, it becomes exactly the Gaussian case (Algorithm 1).
15

---

Published in Transactions Machine Learning Research (8/2025)
Algorithm Algorithm 2: ABIC LiNGAM
Input X E Rexd Observed data
N E Raxd Tnitial error covariance
tol > 0: Convergence threshold
maa iters: Max iterations
h((0) Differentiable constraint
P. a Augmented Lagrangian coeff.
A Regularization weight.
B > 0: Shape parameter (if + 1, non-Gaussian)
Output: 8
Initialization:
1. 80, Ω <- (0, Ia) or random.
2. c c- -m(n).
3. Define
d Σ
20
LS(0) 2n
X, X X6., - 7(0) 22
==1
(Note: 3 ==1 recovers the Gaussian case.)
for = 0 to mar iters - do
(A) Pscudo-variables: for = to d do
- x-x,,-- -
- 0, + E_t +
(B) Parameter optimization:
Solve
+ arg LS(O) + 51000 +a0)) + Σ
tanh (0,1) }
for == to
6 + X., X 8
2 + var
(C) Convergence check:
if ||+11 8'1| || 22l <
tol then
break
return
8
Tips:
Tf # 1, non- non-Gaussian stats help identify directionality strictly.
= 1, same as Algorithm
Use gradient-based methods or coordinate descent for the inner optimization.
3. Identifiability: Orientation Recovery Beyond Gaussian Equivalence:
16

---

Published in Transactions Machine Learning Research (8/2025)
Under non- Non-Gaussian assumptions (3 + 1), ABIC LINGAM can potentially identify causal directions.
Gaussian ABIC recovers structure up to Markov equivalence.
5 Experiments
5.1
Simulation
5.1.1
Simulation Setup and Evaluation
We aim to evaluate whether our proposed method (ABIC LINGAM) can accurately recover bow- free ADMG structures under various conditions, including non-Gaussian error terms. Specifically, we focus on:
How well the method identifies both the skeleton (presence absence of edges) and the direction (arrowhead/tail) (/ail) of cach edge.
The effect of different shape parameters 3 of the multivariate generalized normal distribu- tion (MGGD), which control the deviation from Gaussianity.
Following and extending the setup Bhattacharya et al. (2021), we generate synthetic data from a bowefree ADMG as follows:
1. Node pairs and edge assignment. For each pair of nodes (i.j) with i j, draw uniform random value in [0, 1).
If this value below a predefined threshold for directed edges, X₃ Xj and sample the coefficient dij uniformly in [-2.0, -0.5] U (0.5, 2.0].
If this value is within the threshold for bidirected edges, set X, + Xy and assign Pay = Szz uniformly from [-0.7, -0.4] J 0.74.
Otherwise, no edge is placed between (i.j).
2. Diagonal entries of Q. Each f4 is sampled from an interval =10.4, 0.7]. To ensure Ω positiveddelinte, we add an adjustment term proportional to E( plus an offset in 0.51.
3. Shape parameter and error terms. Let B (1,3,5) be the shape parameter of the MGGD. For each node i, generate error terms from multivariate generalized normal distribution (mean 0, covariance Ω shape B). Note that 8 = recovers the Gaussian case.
4. Generate observed data. Construct X, via
X;
diz X + ei i l. jspa(s)
Σ
Repeating this for samples yields the data matrix X Rexd
17

---

Published in Transactions Machine Learning Research (8/2025)
This probabilistic framework ensures that the randomly generated ô and Q conform bow-free ADMG while allowing flexible edge structures.
We vary:
Sample size (100, 500, 1000);
Number of variables d E (5, 10),
Shape parameter E (1,3,5) (MGGD).
Hence, we obtain 3 x 2 3 = 18 total conditions. For each condition, we run 50 trials. All
simulations use Python 3.8 with NumPy/Scilyy. We used the hyperparameters recommended by ABIC htttps//gitllab.com caw.aaa.ccombbaattaaaa rrbaabattttttooto as a reference when setting up our implementation.
Compared Methods.
We compare the proposed ABIC LiNGAM to:
bow-free ABIC (Bhattacharya et al., 2021), which assumes Gaussian errors (B 1) and recovers structures up to Markov equivalence.
FCI (Spirtes al., 2000), a constraint-based algorithm for ADMGs, outputting partial ancestral graph (PAG).
BANG (Wang & Drton, 2024), a constraint-based method that exploits higher-order mo- ments to identify bow-free ADMG directions (5% significance).
Because synthetic data are bow-free, evaluate ABIC LINGAM under the bow-free constraint. We assess:
Skeleton accuracy (presence absence of edges): Precision, Recall, F1-score.
Arrowhead accuracy (correctly oriented edges): Precision, Recall, Fi-score.
Tail accuracy (the tail end of directed odges): Precision, Recall, Fi-score.
For FCI, which outputs a partial ancestral graph (PAG), we interpret Ao+ B as A + B and a-o as A ++ B.
Key Point (9): By introducing non-Gaussian errors (when f 1), we can potentially distinguish causal directions beyond what is possible in purely Gaussian scenarios. Specifically, ABIC LINGAM leverages score-based continuous optimization under multivariate generalized normal assumption, providing the possibility of identifying orientations rather than merely recovering Markov equiva- lence classes.
5.1.2 Simulation Results and Discussion
In what follows, we report the performance of ABIC LINGAM terms of Recall, Precision, and F1-score three aspects of the learned graph: Skeleton, Arrowhead, and Tail. We compare
18

---

Published in Transactions Machine Learning Research (8/2025)
Category SKELTON
Category ARROWHHEAD
1. Category TAIL
ABICrbow-freet
ABIC LINGAM beta true
Sampll Sa 100
Sample Sue 500
Sample Size 1600
Sampe SizW
Sample Size 500 Sample 1000
Sample 2 300
Sample Sizw 500
Sample Size 1000
Cetegory
Category
Category
Sample
Sampl B 5, Category ARRDWHEAD
Category SKELTON
5, Category TAIL
Sample sm 300
Semple Siew 500
Sample Size 1000
Samp Size 100
Sample Siee 500
Sample Sa 1800
Sample 2: 300
Sample Size 50D
Sample Size 1000
Figure The precision results for each method with five variables.
6 1. Category SKELTON
1. Category ARRDWHEAD
1, Category TAIL
ABIC(bow free)
BANG
ABIC LINGAM bota true
beta
Sample San 100
Simple Sue 500
Sample Size 1000
Samp Size 100
Sample Size 500 Sample Saze 1800
Sample 5:m 300
Sample Siee 500 Sample Stee 1000
SKELTON
ARROWHEAD
6 Category SKELTON
B Category ARROWHEAD
5, Category TAIL
Samole Sze
Somole Size
Sample 1000
Sampl Size 100
Sample Siee 500
Sample Sze 1000
Semple zze 100
Semple Size 500
Sample Size 1000
Figure 3: The recall results for each method with five variables.
19

---

Published in Transactions Machine Learning Research (8/2025)
Category SKELTON
Category ARROWHEAD
1, Category TAIL
ABICrbow-freet
ABIC LINGAM beta true
Samplle 100
Sample Sue 500
Sample Size 1600
Sampe Size
Sample Size 500 Sample 1000
Sample 2 300
Samale Size 500 Sample Size 1000 Category
Cetegory
Category
Sample
Sampl B 5, Category ARRDWHEAD
Category SKELTON
5, Category TAIL
Sample 5m 300
Semple Siew 500
Sample Size 1000
Samp Size 100
Sample Size 500
Sample Sa 1800
Sample 2: 300
Sample Siee 50D
Sample Size 1000
Figure 4: The F1-score results for each method with five variables.
B 1. Category SKELTON
B 1. Category ARRDWHEAD
- 1, Category TAIL
ABIC(bow free)
BANG
ABIC LINGAM bota true
LINGAM beta
Semple San 100
Sample Sue 500
Sample Size 1000
Samp Size 100
Sample Size 500 Sample Sazn 1800
Sample 2: 300
Samale Sie 500 Sample Stee 1000
SKELTON
ARROWHEAD
6 Category SKELTON
== Category ARROWHEAD
5, Category TAIL
Samble Sze 100
Somble Size
Sample 1000
Sampl Size 100
Sample Siee 500
Sample Sze 1000
Sample zze 100
Semple Size 500 Sample Size 1000
Figure The precision results for each method with ten variables.
20

---

Published in Transactions Machine Learning Research (8/2025)
Category SKELTON
Category ARROWHEAD
1, Category TAIL
ABIC(bow-free)
H
ABIC LINGAM beta true
est
Samplle 100 Sample Sue 500 Sample Sizu 1600
Sampe Size 100
Sample Size 500 Sample 1000
Sample 2 300
Sample Size 500 Sample Size 1000
Cetegory
Category
Category
Sample
Sampl B 5, Category ARRDWHEAD
Category SKELTON
5, Category TAIL
Sample Sm 300
Simple Sie 500
Sample Size 1000
Samp Size 100
Sample Siee 500
Sample Sa 1800
Sample 2: 300
Sample Sie 500
Sample Siee 1000
Figure 6: The recall results for each method with ten variables.
6 1. Category SKELTON
B Category ARRDWHEAD
1. Category TAIL
ABiC(bow free)
BANG
ABIC LINGAM bota true
Sample San 100
Sample Sun 500
Sample Size 1000
Sampe Size 100
Sample Size 500
Sample Saze 1800
Sample 2m 300
Sample Siee 500 Sample Siee 1000
SKELTON
ARROWHEAD
6 Category SKELTON
B Category ARROWHEAD
5, Category TAIL
Samble Sze 100
Somble Size
Sample Size 1000
Sampl Size 100
Sample Size 500
Sample Sze 1000
Semple ze 100
Semple Size 500
Sample Size 1000
Figure 7: The Fl-score results for each method with ten variables.
21

---

Published in Transactions Machine Learning Research (8/2025)
method with BANG (Wang & Drton, 2024), as well ABIC (Bhattacharya al., 2021) and FC1 (Spirtes et al., 2000), under various conditions:
n E (100, 500, 1000); de (5,10, B (1,3,5).
Key numerical results are shown in Figures (2)-(7)-
(1) Gaussian case: B ==1.
Skeleton accuracy. ABIC LiNGAM achieves almost the same Skeleton Re- call Call/PrecisionF1 as ABIC. This agrees with the fact that, in the Gaussian regime, both methods can recover structures only up to the Markov equivalence class.
Directional accuracy (Arrowhead/Taii Since = implies normal errors, ABIC LiNGAM does not gain additional information to identify directions uniquely. For in- stance, at (n, d) (500, 5), both methods yield high Skeleton-F1 but yield only limited Arrowhead/ Tail correctness, consistent with standard LiNGAM theory.
(2)
case: B * 1.
Improved direction recovery. Once diverges from 1, ABIC LiNGAM can exploit non- Gaussianity to break Markov equivalences. For example, at (d, n, (5, 500, 3), ABIC LINGAM obtains Arrowhead Recall 0.760, Precision 0.733, F1 0.735, surpassing BANG's Arrowhead Recall of 0.651 (though BANG reaches a high Precision of 0.898).
Consistency with increasing B. We also observe stable performance as 3 grows further, e.g., (d, n, 3) 500,5) yields Arrowhead F1 0.707, close to the = case. This suggests that the proposed method diroction estimation remains robust across different levels of Ganssiasiantty.
(3) Effect of dimensionality (d 5 vs. d d= 10).
Stability of ABIC LINGAM higher d. Comparing d = 5 and 10 n 500, see little degradation in ABIC LiNGAM s Arrowhcad Recall and F1 (c.g., 0.726 0.760). Hence, the method appears well-suited for moderately larger graphs.
Comparison to BANG at d = 10. As d increases, the gap in directional accuracy be- tween ABIC LiNGAM and BANG tends to narrow. For instance, at (d,n, = (10,500, 3), ABIC LiNGAM' s Arrowhead Recall (0.726) exceeds BANG's s (0.590), though BANG still sometimes shows slightly better Tail precision certain sottings.
(4) Overall comparison to existing methods.
BANG generally exhibits high donaltional accuracy, especially for d 5. In particular, with Tail metrics at (d, n, B) (5, 500, 3), BANG shows Recall of 0.945, higher than ABIC LiNGAM 0.684. That said, the difference decreases as grows.
22

---

Published in Transactions Machine Learning Research (8/2025)
ABIC (Gaussian-only) matches ABIC LiNGAM under 1, but under non- non- Gaussianity, it struggles to recover directions (o.g., Tail Recall 0.055 at d 5, 500, 3 3).
FCI works well in identifying skeletons and partial ancestral structures, but reading off precise arrowheads from the PAG can be more involved.
(5) Summary of key findings.
In summary, ABIC LiNGAM can accurately estimate causal structures structures- inchuding orientations when errors follow a non-Gaussian distribution, even the presence of unmeasured variables. behaves similarly to ABIC for 8 = (Gaussian), and approaches slightly surpasses BANG in certain hon-Gaussian conditions, especially the dimensionality increases. Meanwhile, BANG still shows very high directional accuracy in lover-dimensional settings. Hence, our results suggest that ABIC LiNGAM is competitive with the state-of-the-art.
These observations are based on the aggregated numerical values (Figures 1-6). In particular, we highlight Arrowhead and Tail metrics because direction estimation is crucial for causal infer- ence. By exploiting score-based unified criterion rather than repeated independence lests, ABIC LiNGAM maintains stable and consistent accuracy across different sample sizes, dimensions, and non-Gaussian parameters.
Although larger samples theoretically lead to more accurate recovery of the true structure, we ob- served cases where increasing n improves recall but lowers precision (Figure 3). Further experiments and discussions on this phenomenon are provided in Appendix.
Key Point (10): In many conditions, non-Gaussion error terms (3 * 1) can help ABIC LiNGAM surpass purely Gaussian-based methods (e.g., ABIC) in direction identification, because they break the Markov equivalences that persist under Gaussian assumptions.
Key Point (11): Our experiments suggest that as dimensionality increases, the performance gap between ABIC LiNGAM and BANG often narrows, indicating potential scolability and robust performance of the proposed approach higher-direensional settings.
5.2 Performance Evaluation on weal-world Data
We evaluated our proposed method using a sociological data repository (https:/ htaps://gss.norc.rrg.. norc.org/), which has also been studied in the context of DirectLiNGAM (Shimizu al., 2011). The dataset contains 1380 samples of sociological variables, such parental pntraationoccupatiion and off- spring's outcomes. In this study, we focus the subset of variables and their presumed causal directions depicted in Figure 8, following the domain knowledge and temporal ordering from Duncan & Featherman (1972).
Compared Methods.
Similar to our simulation experiments, we compared:
bow-free ABIC LINGAM (proposed), which handles non-Gaussian errors under
bow-free ADMG constraint;
bow-free ABIC (Bhattacharya et al., 2021), the bascline Gaussian version;
23

---

Published in Transactions Machine Learning Research (8/2025)
FCI (Spirtes al., 2000), constraint-based algorithm for ADMGs, producing a partial ancestral graph (PAG);
BANG (Wang & Drton, 2024), which also identifies bow-free ADMG directions using Gaussianity at a 5% significance level.
These methods were applied in the same fashion as the simulation, locusing on how well they recover known directions and possible confounding relationships.
Estimation of the Shape Parameter B. To remain consistent with earlier assumption that the error terms follow a multivariate generalized normal distribution (MGGD), we needed to estimate the shape parameter 3 from the observed GSS data. Initially, we attempted to estimate a single global 3 using the entire dataset, but the resulting value was below 1 and produced unstable estimates (ō became all zeros, and turned into a near-diagonal matrix). Instead, we adopted a variable- variable
barriabbb cstimation approach:
1. For each variable X,, fit an MGGD- model separately to approximate its marginal distribution and infer an individual shape parameter 3r.
2. Choose the maximum among (3:) as a conservative estimate, B = max; &t.
3. Use this 3 the shape parameter for the error distribution in the bow-free ADMG model.
This procedure avoided the instability observed the global estimation method, ensuring all vari- ables were effectively represented in the final choice of B.
Implementation Details. With 3 fixed, ran bow-free ABIC LINGAM under the same aug- mented Lagrangian optimization framework described previously. For ABIC, BANG, and FCI, we applied their respective default or recommended settings, mirroring the simulation setup. Finally, we evaluated each learned graph against the domain-consistent directions in Figure 8
Remarks. By leveraging both domain knowledge (Duncan Featherman, 1972) and our non- Gaussian ADMG framework, we expected ABIC LiNGAM to capture latent confounders via bidi- rected edges and to estimate orientations more accurately than purely Gaussian methods. In the next section, we present the quantitative results and discuss how they compare to the simulation findings.
Table 2: Experimental Results Bollen Data
Method
SKELETON
ARROWHEAD
TAIL
Recall
Precision Fl-score
Recall Precision Fl-score
Recall
Precision
Fl-score
FCI
0.727
1.000
0.842
0.643
0.692
0.666
0.250
0.667
0.363
BANG
0.727
0.889
0.800
0.643
0.643
0.643
0.250
0.500
0.333
ABIC
0.455
1.000
0.625
0.429
0.600
0.500
0.000
0.000
0.000
ABIC LINGAM
0.818
1.000
0.900
0.714
0.800
0.740
0.500
0.800
0.615
bota est.
Table2 reports the performance metrics for bow-free ABIC LINGAM (beta est) compared to other methods. We focus three evaluation criteria: Skeleton, Arrowhead, and Tail.
24

---

Published in Transactions Machine Learning Research (8/2025)
Father's Education (x3)
Son's Education (x5)
Father's Occupation (x1)
Son's Income (x2)
Son's Occupation
Number of Siblings
(x4)
(x6)
Figure 8: Variables and causal relations in the General Social Survey dalaset used for the eva ivaluation.
Father's Education (x3)
Son's Education (x5)
Father's Occupation (x1)
Son's Income (x2)
Son's Occupation (x4)
Number of Siblings
(x6)
Figure Causal graph produced by ABIC LINGAM: The dashed lines represent predicted arrows that differ from the true arrows.
25

---

Published in Transactions Machine Learning Research (8/2025)
Skeleton Performance. Our proposed method achieves Recall of 0.818, Precision of 1.000, and an Fi-score of 0.900, surpassing all baseline methods in identifying the pres- ence absence of edges. In contrast, ABIC attains Recall of 0.455 with an F1 of 0.625. This highlights that ABIC LaNGAM (beta est) can more consistently recover the correct adjacency relationships than the purely Gaussian version of ABIC.
Arrowhead and Tail Performance (Direction). When examining directed-edge ori- entation, ABIC LINGAM (beta est) also shows promising results:
Arrowhead: Recall 0.714, Precision 0.800, F1 0.740. Notably, this exceeds BANG Arrowhead metrics (Recall 0.643, Precision 0.643, Fl 0.643).
Tail: Recall 0.500, Precision = 0.800, F 0.615. Again, this is higher than BANG' Recall (0.250), Precision (0.500), and F1 (0.333).
As such, even for orientation-specific measures, ABIC LINGAM (bela est) appears to out- perform both BANG and ABIC on this particular dataset. Figure ?? (formerly Figure 9) suggests that the true causal structure is largely recovered, including plausible directional- ity.
Overall, these results demonstrate that ABIC LINGAM (beta est) can estimate causal structures (including edge directions) with higher accuracy than competing methods the GSS datasct. In particular, its Skeleton F1-score (0.900) and the improvement in Arrowhead/ Tail measures over BANG and ABIC indicate the potential effectiveness of adopting a non- non-Gaussian error assumption in a score-based, bow-free ADMG framework. This is consistent with our simulation findings that exploiting non-Gausiayytty is crucial for precise direction inference. ABIC LINGAM (beta est) can flexibly incorporate prior knowledge. In this dataset as well, incorporating prior knowledge led to improved accuracy, See Appendix for further details.
Intuition: The shape parameter B can significantly influence direction recovery. In particular, when 3 > 1, the distribution exhibits heavier luils or skewness, thereby highlighting causal direc- tions that might remain indistinguishable under purely Gaussian assumptions (3 = 1).
6 Conclusion
In this study, we built on the score-based continuous optimization method proposed by Bhattacharya et al. (2021) and introduced ABIC LINGAM, an extension of LiNGAM designed for causal structure estimation in the presence of unmeasured variables. By assuming that error terms follow a multivariate generalized normal distribution, we demonstrated that we can poten- tially identify not only the Markov equivalence class but also the directions in the causal structure. Moreover, because ABIC LiNGAM accurately estimate the SKELETON even when the data follow Gaussian distribution, our approach can be viewed a more generalized version of the method in Bhattacharya el al. (2021).
We also proved that parameters in bow-free ADMGs are almost everywhere identifiable from the covariance matrix under the multivariate generalized normal distribution, thereby extending the results of Brito & Pearl (2002) to non-Gaussian settings. Through simulations and experiments with real-world data, we confirmed that the proposed method achieves accuracy comparable to existing methods (such as BANG or FCT) in recovering causal structures, including causal directions. These
26

---

Published in Transactions Machine Learning Research (8/2025)
findings suggest that our approach provides useful framework for causal discovery in practical settings where unmeasured variables may be present.
Concretely, our main contributions can be summarized as follows:
Direction Identification under Non-Gaussianity: Whereas many Gaussian-based methods only recover the Markov equivalence class, method leverages the mul- tivariate generalized normal distribution (B + 1) to rigorously identify causal directions.
Thcorctical Identifiability: We proved that parameters in bow-free ADMGs are al- most everywhere identifiable from the covariance matrix, thereby exlending the results of Brito & Pearl (2002) to non-Gaussian settings.
Empirical Validation: Using both simulations and real datasets, we showed that our method achieves accuracy comparable to that of existing methods (such as BANG or FCI) in recovering causal structures and directions.
In future work, we will explore ways to reduce estimation time, investigate extensions to nonlinear data structures, and apply our method mixed data that include discrete variables. Advancing these directions will help establish more efficient and versatile methodology for causal discovery across diverse real-world scenarios.
In this paper, we extended the differentiable score-based approach from the Gaussian to the non- Gaussian setting. This does not imply removing all assumptions; rather, the non-Gaussian error terms furnish additional identifiable information beyond second moments, thereby enabling us to determine causal directions more strictly. The results underscore that adopting nou-Caussiy not just a weaker constraint but critical tool enhance identifiability causal discovery under latent confounding.
7
Acknowledgments
This work was supported by the Japan Science and Technology Agency (JST) under CREST Grant Number JPMJCR22D2.
References
D. Bernstein, B. Saeed, C. Squires, and C. Uhler. Ordering-based causal structure learning the presence of latent variables. In Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS), pp. 4098-4108, 2020.
R. Bhattacharya, T. Nagarajan, D. Malinsky, and I. Shpitser. Differentiable causal discovery under unmeasured confounding. In International Conference Artificial Intelligence and Statistics, pp. 2314 2322, March 2021.
C. Brito and J. Pearl. new identification condition for recursive models with correlated errors. Structural Equation Modeling, 9(4):459 474, 2002.
W. Chen, R. Cai, K. Zhang, and Hao. Causal discovery in linear non-gaussian acyelic model with multiple latent confounders. IEEE Transactions on Neural Networks and Learning Systems, 33(7):2816-288 2827, 2021.
27

---

Published in Transactions Machine Learning Research (8/2025)
D. M. Chickering, Learning Bayesian networks is NP-complete, pp. 121 130. 1996.
D. M. Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3(No)):507 554, 2002.
T. Claassen and 1. G. Bucur. Greedy equivalence search in the presence of latent confounders. In Proceedings of the 38th Conference on Uncertainly in Artificial Intelligence (UAT), pp. 2022.
D. Colombo, M. H. Maathuis, M. Kalisch, and T. S. Richardson. Learning high digh-dimensional directed acyelic graphs with latent and selection variables. The Annals of Statistics, pp. 294 321, 2012.
M. Dorton, M. Eichler, and T. Richardson. Computing maximum likelihood estimates in recursive linear models with correlated errors. Journal of Machine Learning Research, 10, 2009.
O. D. Duncan and D. L. Featherman. Psychological and cultural factors in the process occupa- tional achievement. Social Science Research, 1(2):121-145, 1972.
Gomez-Villegas, E. Gómez, M. A. Gomez- Villegas, and J. M. Marín. A multivariate generalization of the power exponential family of distributions. Communications in Statistics- Seeiistics-Thorry and Methods, 27(3): 589-600, 1998.
1. Ng, B. Huang, and K. Zhang. Structure learning with continuous optimization: sober look and beyond. In Proceedings of the 3rd Conference on Causal Learning and Reasoning (CLeaR), pp. 1-35, 2024.
C. Nowzohour, M. H. Maathuis, R. J. Evans, and P. Bühlmann. Distributional equivalence and structure learning for bow-free acyclic path diagrams. Flectronic Journal of Statistics, 11(2): 5342 5374, 2017.
J. M. Ogarrio, P. Spirtes, and J. Ramsey. A hybrid causal scarch algorithm for latent variable models. In Proceedings of the 8th International Conference on Probabilistic Graphical Models (PGM), pp. 368 (88-379, 2016.
M. Okamoto. Distinctness of the eigenvalues of quadratic form in multivariate sample. The Annals of Statistics, pp. 763 765,765,11 1973.
Pearl. Causality. Cambridge University Press, 2009.
T. Richardson and P. Spirtes. Ancestral graph markov models. Anmals of Statistics, 30(4):962-103) 1030, 2002.
S. Shimizu, O. Hoyer, A. Hyvärinon, and A. Kerminen. A lincar non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Reseurch, 7(10), 2006.
S. Shimizu, T. Inazumi, Y. Sogawa, A. Hyvarinen, Y. Kawahara, T. Washio, and P. Hoyer. Di- rectlingam: A direct method for learning linear non-Gaussian structural equation model. Jour- 1111222-888,,--111111111111222
nal of Machine Learning Research, 12:1225 1248, 2011.
1. Shpitser and J. Pearl. Identification of joint interventional distributions in recursive semi- Markovian causal models. Proceedings of AAAI, pp. 1219-1226, July 2006.
28

---

Published in Transactions Machine Learning Research (8/2025)
Spirtes, C. N. Glymour, and R. Scheines. Causation, Prediction, and Search. MIT Press, 2000.
T. Tashiro, S. Shimizu, A. Hyvärinen, and T. Washio. Parcelingam: A causal ordering method robust against latent confounders. Neural Compulation, 26(1):57 83, 2014.
T. Verma and J. Pearl. Equivalence and synthesis of causal models. Proceedings of the 6th Annual Conference on Uncertainty in Artificial Intelligence, 1990.
Y. S. Wang and M. Drton. Causal discovery of linear Mon-Gaussian causal models with unobserved confounding. Journal of Machine Learning Research, 21:1329 1359, 2024.
S. Wright. Path coefficients and path regressions: alternative complementary concepts? Bio- metrics, 16(2):189 202, 1960.
X. Zheng, B. Aragam, P. K. Ravikumar, and E. P. Xing. DAGs with no lears: Continuous optimiza- tion for structure learning, Advances in Neural Information Processing Systems, 31:9472-990 199422--483, 2018.
A Appendix
A.1 Proof of the Identifiability of Parameters in Bow-Free ADMGs with Multivariate Generalized Normal Distributions
Theorem 1
Let G be a bow-free ADMGs with error terms following a multivariate generalized normal distri- bution, and let the set of parameters of be = {6,2). Then, for almost all 0, the following holds:
E()) E(0)
implies 0 = 0'.
In other words, if two parameter 0 and 0' give the same covariance matrix Σ, then 0 and 0' must be identical, except possibly when 0 belongs to set of Lebesgue measure zero. If the two lemmas described later be proven for the case where the error terms follow multivariate generalized normal distribution, this theorem can be demonstrated using the same proof as in Brito & Pearl (2002).
Definition
path in a graph is sequence of edges (directed or bidirectional), where each edge starts from the node where the previous edge ends. A directed path consists only of directed edges all pointing in the same direction. A node X is called to Y. A path is said to be blocked if there is node Z on the path such that there are consecutive edges pointing toward Z (e.g., ... + Z + ...). In this case, Z is called collider.
ancestor of a node Y if there is directed path from X
Definition
In DAG, the depth of node defined as the length (number of edges) of the longest path directed from its ancestors to that node.
29

---

Published in Transactions Machine Learning Research (8/2025)
Lemma 1
Let X and Y be the nodes in bow-free ADMG with (X) 2 depth(Y). Then, all paths between X and Y that include node Z satisfying depth 2 depth(X) are blocked by colliders. This lemma is based on graph theory and does not depend on the distribution of the error terms. It is quoted from Brito & Pearl (2002).
Definition
For each node Y, the set of edges directed to Y, denoted by I(Y,, is defined as the union of the following two sets: (a) the set of all directed edges pointing to Y. (b) the set of all bidirectional edges between X and Y, where depth(X) < depth(Y).
Lemma 2
Let be a variable at depth & in a bow-free ADMG. Assume that the parameters of all edges connocting variables depth less than are identifiable. Then, almost all cases, the parameters of each edge in the set I(Y) are identifiable.
Proof. In Brito & Pearl (2002), the identifiability of bow-free models was established under the assumption that the error terms follow a multivariate normal distribution. In this study, we extend this identifiabiiity result to the case where the error terms follow the aforementioned multivariate generalized normal distribution. The proof itself draws heavily from Brito & Pearl (2002).
Wright's 's method Wright (1960) relies on linear relationships between variables and covari- ance structures. Since the multivariate generalized normal distribution is closed under linear transfornations((Gmme et al., 1998)), Wright's method is applicable beyond the normal distri- bution as long as the necessary linear conditions are satisfied. Indeed, Wright (1960) also mentions that. Wright's 's method can be applied to distributions other than the normal distribution.
Let X X2,, ; } be the set of variables with depth less than k, and suppose that these variables are connected to Y by directed undirected edges. By the properties of bow-free ADMGs, a one-to-one correspondence exists between each variable in X and the edges in T(Y). Therefore, ((Y) can be expressed as
I(Y) = ((x,,) ().... (Xm.:..)..))))) Y)).
Applying Wright's mothod to each pair (X,, Y) yiclds the following equations:
OX,Y = Σ T(P), i= 1,...,m m
where ox,y represents the covariance between X, and Y, the sum over all paths pt botween X; and Y that have direct indirect effects associations, and T(p:) represents the product of parameters along the path p,.
For each i, let X; be the parameter corresponding to the edge (X, The equation can be rewritten as
OX. >; +
Σ Xyaij 1, 1,..., m
30

---

Published in Transactions Machine Learning Research (8/2025)
where the coefficients di are functions of identifiable parameters corresponding to edges connecting variables of depth less than k. These coefficicnts reflect contributions from direct or indirect effects associations involving known parameters, excluding the direct edge (X,
Under the assumption, by the induction hypothesis, that all parameters of edges connecting vari- ables of a depth less than k are identifiable, the coefficients aij are known quantities. Therefore, we obtain a system of m linear equations with m unknowns , Am; which can be written a matrix form as
AA
(24)
where
xxiY xx2Y
A1 A
d12
01m
Q21
1
a2m
A
A
xx...
0mml
0m2
To establish the identifiability of the paramcters A, it suffices to show that matrix A is invertible in almost all cases, that is, (A) jf except on a set measure zero, considering that the left-hand side is observable. The matrix A has all diagonal elements equal and off-diagonal elements depending the model parameters. The determinant can be expressed in terms of the diagonal and off- diagonal elements, shown in (25).
det ()) = 1
(25)
where T is either zero or a polynomial in the model parameters that do not contain any constant term.
According a well-known result in algebraic geometry Okamoto (1973), the set of parameter values where (A)(A) has Lebesgue measure zero in the parameter space. This is because det(A) 0 defines algebraic variety of lower dimension within the parameter space. Therefore, the matrix A is invertible in almost all cases, and the system of linear equations has unique solution.
Thus, under the given assumptions, each parameter X; is identifiable in almost all cases.
A.2 Hölder's Inequality
Hölder's inequality is a fundamental result in analysis that provides estimates for sequences (or more generally, measurable functions on a mcasure space (4) in terms of their L Specifically, for p. g 2 satisfying = 1, Hölder's inequality states that for any two sequences (ax) and (bx)::
1/p
1/4
80
a byl $
Σ 1=1
laml))
by
(26)
1=1
31

---

Published in Transactions Machine Learning Research (8/2025)
Furthermore, by taking 1, we obtain useful inequality for finite sums as follows:
Σ
Jax
S
aall..
(27)
This special case reflects how the 7--norm behaves finite setting and central understanding the interplay between norms and summation.
In specific problem, Hölder's inequality to handle the terms involving 3th powers of squared residuals. Considering the log-likelihood expression after rearrangement,
{|xx)
N
log Qzi
1
N
x,0
Σ
Σ
2
20
11.-
1=1
jEPa())
k(ep(1)
N
log Qúi-i
N8-1
Σ
x
Σ
Σ
Di 2(10 Z(
2
2n
N8-1
11.
jEpa(i)
kEap(i)
(28)
We introduce the factor the sequence as
No No-T to rewrite the sum in a form amenable to Hölder's inequality. Define
2/3
x
Σ
Σ
2(00
Q1
jEpa(1)
kEsP(1)
and let by = 1. we choose = B, and hence q = gives us
+ 1 1), Hölder's incquality B-1 (so that
$
larbel
#
by
at
N 11
Rearranging this inequality, obtain lower bound on Σ a₁ in terms of 8 and the L -orm of the residuals
ΣN IX! -
Σ (xx))...)))
N3-1
Substituting this bound back into the expression for E(u, E, X(X), have
32

---

Published in Transactions Machine Learning Research (8/2025)
28
N
((u, ((xx) 2 log fai-ii 2n2 14
X;
Σ jEpa(5)
K,
Σ
Di,k
(29)
N8-1
kEsp(1)
Thereby, Hölder's inequality is employed to provide a nontrivial lower bound on the log-likelihood by relating sums of 5th powers of squared terms to the 5th power their L -norm, scaled appropriately by N This facilitates a more tractable analysis of the growth behavior and bounding propertics of the likelihood function.
A.3 Additional Experiments with Larger Sample Sizes
In addition to the simulations deseribed in Section of the main text, we conducted fur- ther experiments to investigate how our score-based continuous oplimization approach (ABIC LiNGAM) behaves under larger sample sizes. Specifically, we generated datasets of size {100, 500, 1000, 2000, 5000) with five observed variables (dimension & = 5) and fixed non-Caussian parameter 5 = 5. We considered both cases: beta ture (i.e., 3 is known) and beta est (i.e., 8 is estimated from the data). Our primary goal was to examine whether increasing n leads to higher recall while potentially decreasing precision due to "over-detection" of minor dependencies.
Table 3 presents the results for three types of edges (SKELETON, ARROWHEAD, TAIL), evalu- ated by recall, precision, and F1-score. One would expect, from theoretical standpoint, that as the sample size grows, the estimated causal structure should approach the true structure. However, seen in Table 3, alt althou recall does lend to increase with larger n, precision can decline in some settings, leading to F' Fl-scores that do not strictly increase. This phenomenon arises because a larger sample size increases the statistical power our score-based mothod, causing it to detect even subtle or spurious edges as false positives. Such a trend is common flexible models when the penalty strength or thresholds are not tuned for very large sample sizes.
As Table 3 shows, recall generally improves with (i..,, fewer missod edges), but precision (< the proportion of correctly identified edges among those detected) can decrease, indicating rise in false positives. Consequently, the F1-scores (F1 = 2 x recall precision (reciion///rcalll (recall + precision)) do not necessarily increase monotonically. Moreover, because our causal graph is randomly generated for each run, the "difficulty" of each dataset can vary, causing some fluctuations (e.g., at n 500
vs. n = 1000).
These results help explain why large-sample outcomes can sometimes appear deviate from the theoretical expectation of "convergence to the true graph. In principle, if the regularization or threshold hyperparameters are adjusted for the sample size, false positives may be reduced, thus im- proving precision alongside recall. However, in our current implementation, these hyperparamcters are held constant, leading to potential "over-detection when n grows large.
In Section 5 (Figure 3), the last subgraph similarly shows case in which incrcasing sample size does not strictly improve Fi-scores, despite recall gains. As clarified above, our score-based approach tends pick up spurious edges with large n,, thus reducing precision and causing mismatch from purely theoretical expectations. We believe that this phenomenon arises from not tuning the penalty term to suppress minor spurious connections in very large samples.
Relation to Figure 3.
33

---

Published in Transactions Machine Learning Research (8/2025)
Table 3: Extended experiment results of ABIC LINGAM (8 = - 5) with sample sizes n E {100, 500, 1000, 2000, 5000] and dimension = Metrics are reported for three edge categories (SKELETON, A ARROWHEAD, TAIL), under both the "B-est" (top) and "3-true" (bottom) con- ditions. Numbers are averaged over multiple random simulations. Note that data generation relies on randomly assigned causal graphs and random noise, causing minor variability across runs.
SKELETON
ARROWHEAD
TAIL
Sample Size
Recall
Precision
FT
Recall Precision
F1
Recall
Precision
F1
(a) ABIC LiNGAM beta true
N = 100
0.903
0.776
0.826
0.722
0.682
0.691
0.669
0.443
0.502
n = 500
0.948
0.774
0.843
0.727
0.676
0.691
0.643
0.357
0.430
n 1000
0.965
0.695
0.800
0.739
0.584
0,641
0,647
0.315
0.404
n 2000
0.985
0.713
0.821
0.768
0.619
0.673
0.699
0.328
0.429
TI 5000
0.994
0.720
0.826
0.771
0.597
0.661
0.673
0.303
0.403
(b) ABIC LINGAM M. bela est
71 100
0.883
0.806
0.832
0.706
0.682
0.680
0.613
0.452
0.490
500
0.942
0.821
0.869
0.742
0.727
0.720
0.673
0.409
0.485
т 1000
0.928
0.727
0.808
0.739
0.634
0.672
0.603
0.359
0.431
π 2000
0.953
0.750
0.833
0.734
0.670
0.689
0.660
0.364
0.451
n 5000
0.992
0.731
0.832
0.744
0.616
0.660
0.711
0.324
0.424
Future work will focus on dynamically adjusting thresholds and penalties to maintain a balance between recall and precision, ensuring more consistent alignment with the theoretical asymptotic behavior.
A.4 Performance Evaluation on Real Data Using Prior Knowledge
In ABIC LINGAM, prior knowledge can be incorporated into the inference process following the code implementation of ABIC available at https://gitlab..oo http://aaa..ooobbaaataaaaaa/a//e Specifically, by defin- ing hicrarchical causal orders (tiers) and incorporating prior knowlodge that certain variables are unconfounded, we can restrict the parameter ranges (bounds) considered during the estimation. Consequently, edges that contradict the causal order, as well bidirectional edges among uncon- founded variables, are represented with predetermined constraints in the parameter space (e.g., fixed to zero) and are thus automatically excluded during the inference process. Thereby, explicitly reflecting prior knowledge in the model ensures that assumptions regarding causal directionality and the absence of confounding factors are maintained, enabling an efficient structural estimation.
For example, the General Social Survey dataset obtained from sociological data repository httpss://sssnoor..rgg includes variables pertaining to parents and children. Tf causal path exists from the parents to the children, causal path from the children to the parents is not possible. In the following section, consider the inferences with hierarchical causal structure. Specifically, we divided the variables into two groups with two-tiered structure-Fatte s Education and Number of Siblings, Son Education, Son Occupation, Son Income-atd and prohibited the existence directed edges from children' s variables to parents' variables.
Occupation, Father'
As demonstrated by the results shown in Figure 10 and Table incorporating prior information into ABIC LINGAM leads to higher ostimation accuracy compared to the approach without such
34

---

Published in Transactions Machine Learning Research (8/2025)
Father's Education (x3)
Son's Education (x5)
Father's Occupation (x1)
Son's Income (x2)
Son's Occupation (x4)
Number of Siblings
(x6)
Figure 10: Causal graph produced by ABIC LiNGAM incorporating prior knowledge: The dashed lines represent predicted arrows that differ from the true arrows.
Table 4: Experimental Results for Bollen Data incorporating prior knowledge
Method
SK SKELETON
ARROWHEAD
TAIL
Recall Precision F1 Fl-score
Recall Precision -l-ccore
Recall Precision F1-score
ABIC LiNGAM
0.818
1.000
0.900
0.714
0.800
0.740
0.500
0.800
0.615
beta est
ABIC LiNGAM
beta est. prior knowledge
0.818
1.000
0.900
0.692
0.900
0.782
0.857
0.750
0.800
The estimated graph aligns closely with the presumed true causal structure. This suggests that incorporating prior knowledge into ABIC LiNGAM feasible and that leveraging such information can yield improved accuracy in practical applications.
information.
35


Reviews:
[{"review_id": "qBB9NTlOR4", "reviewer": "Reviewer_KWDp", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The main weakness are twofold.\n\n1) The contribution itself is a trivial extension of existing prior work that derived continuous differentiable objective for bow-free ADMGS. There are also no new theoretical results that show that optimization of the continuous objective will derive the optimal ADGM under specific assumptions.\n\n2) The utility of multivariate generalized normal distributions is not motivated and dubious as far as I can tell. I'm not sure what advantage this confers over standard Gaussian distributions, and doesn't belie any standard practical situation (e.g., nonparametric conditions like sub-Gaussian)\n\nHence, I don't feel like there is a significant enough contribution to merit acceptance.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors provide a continuously differentiable objective for learning causal structures for multivariate generalized normal distribution, that is a direct extension of prior work that derived learning of bow-free graphs (a generalization of acyclic graphs to account for undirected edges).", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "I think a fundamentally new contribution has to be made for this paper to be worth accepting --- perhaps addressing one of:\n1. Practicality: what are the real data/simulations where relaxing the assumptions on the noise terms make a big difference? What are assumptions on the noise terms that actually match a practical situation?\n\n2. Theoretical contribution: is it possible to prove stronger theoretical results about the guarantees of existing algorithms in novel settings, or to improve their guarantees?", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "6CYjWDGs9R", "reviewer": "Reviewer_DUVv", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\n1. The manuscript addresses the important and challenging problem of causal structure estimation in the presence of unmeasured confounders. \n\n2. The authors propose a score-based method for estimating bow-free acyclic directed mixed graphs (ADMGs) under a linear non-Gaussian acyclic model. \n\n3. They also provide theoretical guarantees by proving identifiability results when the error terms follow a multivariate generalized normal distribution. \n\n4. The overall structure of the manuscript is clear and well-organized.\n\n**Weaknesses**\n\nA notable limitation of the study is the simplicity of the model setup. In particular, the assumption that the noise terms follow a multivariate generalized normal distribution may be too restrictive and could limit the generalizability of the proposed method.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose a novel score-based method, ABIC-LiNGAM, for identifying causal structures with latent variables under the linear non-Gaussian acyclic model (LiNGAM). By assuming that the error terms follow a multivariate generalized normal distribution, the method enables the estimation of not only the Markov equivalence class, but also the causal directions within the structure. The authors further provide theoretical guarantees by proving the identifiability of parameters in acyclic directed mixed graphs (ADMGs) under the assumed error distribution. To demonstrate the effectiveness of their approach, they conduct extensive simulation studies as well as experiments on real-world datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. The manuscript includes a limited discussion of relevant literature. It would benefit from citing more related work on score-based methods for causal discovery, particularly earlier contributions on identifying ADMGs. For the cited works discussed in more detail (e.g., Wang & Drton (2024); Bhattacharya et al. (2021)), the authors should explicitly clarify the differences from the present study and highlight the contributions and advantages of their method.\n\n2. There are several instances of repetitive descriptions that could be streamlined. For example, the statement “Brito & Pearl (2002) demonstrated that a bow-free ADMG model is almost always identifiable from the observed covariance matrix” appears in both Section 2.3 and Section 2.31. In addition, some key technical terms lack sufficient explanation and should be better clarified for accessibility.\n\n3. Some elements, such as Algorithm 2 (ABIC-LiNGAM), appear to be reused or closely resemble content from the article Differentiable Causal Discovery Under Unmeasured Confounding. The authors should clarify whether these are original contributions or adapted from prior work.\n\n4. In Equation (15), is the notation $Z_{-i}$ consistent with the $Z_{-i}$ used in the first line following the equation? Clarification would help prevent potential confusion.\n\n5. The trace operator notation, $tr()$, should be clearly defined when first introduced, such as in Equation (22).\n\n6. The algorithms presented would benefit from more detailed explanations or better formatting (e.g., spacing or step annotations) to improve readability. Specifically, the differences between Algorithm 1 and Algorithm 2 should be clearly articulated.\n\n7. Some descriptions in the manuscript lack structure. It is recommended to use \"Remarks\" to emphasize key references, results, or insights to enhance the clarity and readability of the manuscript.\n\n8. The simulation results do not appear to exhibit expected large-sample properties. For instance, the last subgraph in Figure 3 seems inconsistent with theoretical expectations. Can the authors provide further explanation?\n\n9. It would strengthen the theoretical contribution if the authors could include an additional theorem to show that ABIC-LiNGAM can consistently recover the skeleton of the causal graph.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "vR5mP4pL2V", "reviewer": "Reviewer_4ipP", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nThe paper is overall strong, although I will state upfront that causal discovery is slightly outside my primary area of expertise. I will let other reviewers comment on the technical details of proofs. That said, the authors make tightly scoped technical claims, provide what appears to be sufficient theoretical and empirical evidence, and provide a thorough discussion of how this work fits within the broader context of structure learning. From a somewhat outsiders perspective, it appears to be a natural extension of two existing and powerful ideas; casting structure learning as continuous optimization and exploiting non-Gaussian noise to disambiguate edge orientation.\n\nWeaknesses:\nWhile the paper is presented clearly, it's not clear how the empirical study substantiates the claims made. The differences between ABIC LINGAM and BANG appears very minimal in almost all of the empirical results. Is there another reason one should prefer ABIC LINGAM? The results do appear to be better in the \"real data\".\n\n\nAt various points it was a bit difficult to distinguish between what was background and what were novel contributions. Concretely, I would recommend moving some of the discussion of Bhattacharya et al. (2021) in Section 4.2 to an early section on Background, but this is really a minor suggestion.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper, the authors present ABIC LINGAM, a novel method that brings insights from recent continuous optimization-based approaches to graph structure learning to the non-Gaussian setting, extending to multivariate generalized normally distributed error terms. In doing so, the proposed methods go beyond just Markov equivalence classes, breaking what would otherwise be symetric in the likelihood under Gaussian errors. The authors prove identification, and present experiments validating their approach.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See Weaknesses above.\n\nI don't see a section in this review portal for clarifying questions, so I will ask here.\n\nCould the authors please comment on how ABIC LINGAM will perform when the parametric assumptions of multivariate generalized normally distributed error terms are violated? Specifically, what happens when errors are actually Gaussian? What about when structural functions are nonlinear? What about structural functions are linear, but errors are non-Gaussian in some other way?", "broader_impact_concerns": "I have no broader impacts concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)