Causal Graph Learning via Distributional Invariance of Cause-Effect Relationship
Authors: Nang Hung Nguyen
Phi Le Nguyen
Thao Nguyen Truong
Trong Nghia Hoang
Masashi Sugiyama
Venue: Accepted by TMLR ()

Abstract: This paper introduces a new framework for recovering causal graphs from observational data, leveraging the fact that the distribution of an effect, conditioned on its causes, remains invariant to changes in the prior distribution of those causes. This insight enables a direct test for potential causal relationships by checking the variance of their corresponding effect-cause conditional distributions across multiple downsampled subsets of the data. These subsets are selected to reflect different prior cause distributions while preserving the effect-cause conditional relationships. Using this invariance test and exploiting an (empirical) sparsity of most causal graphs, we develop an algorithm that efficiently uncovers causal relationships with quadratic complexity in the number of observational variables, reducing the processing time by up to $25\times$ compared to state-of-the-art methods. Our empirical experiments on a varied benchmark of large-scale datasets show superior or equivalent performance compared to existing works, while achieving enhanced scalability.

Published in Transactions on Machine Learning Research (01/2026)
Causal Graph Learning via Distributional Invariance of Cause-Effect Relationship
Nang Hung Nguyen
The University of Tokyo, Japan
Phi Le Nguyen
lerpp@soict.hut..commm
Institute for AI Innovation and Societal Impact (AIİLIFE), Vietnam Honoi University Science ond Technology, Vietnam
Thao Nguyen Truong
ngguyen.truonggaist..gg
National Institute of Advanced Industrial Science and Technology (AIST), Japan
Trong Nghia Hoang
trongnghiaa hoong Cunsu.uuu....m edu
School of Electrical Engineering and Computer Science, Voiland College of Engineering and Architecture, Washington State University, Pullman, Washington, US
Masashi Sugiyama
RIKEN, Japan
The University of Tokyo, Japan
Reviewed on OpenReview: https: openreveew. net/ forum? id- By380882Xe
Abstract
This paper introduces new framework for recovering causal graphs from observational data, leveraging the observation that the distribution of effect, conditioned its causes, remains invariant to changes in the prior distribution of those causes. This insight enables direct test for potential causal relationships by checking the variance of their corresponding effect-cause conditional distributions across multiple downsampled subsets of the data. These subsets are selected to reflect different prior cause distributions, while preserving the effect-cause conditional relationships. Using this invariance test and exploiting an (empirica(") sparsity of most causal graphs, we develop an algorithm that efficiently uncovers causal relationships with quadratic complexity in the number of observational variables, reducing the processing time by up 25x compared to state-of-the-art methods. Our empirical experiments varied benchmark of large-scale datasets show superior equivalent performance compared to existing works, while achieving enhanced scalability.
1 Introduction
A core challenge in causal learning is finding a directed acyclic graph (DAG) that captures the cause-effect relationships between variables in an observational dataset (Zanga et al. 2022 A direct approach to uncover these causal links is through intervention which conducts more experiments confirm whether changes set of variables will consequently affect the outcome distribution of another variable (Peters al. 2016 Guo et al. 2024h. However, such interventional experiments be prohibitively expensive. Furthermore, recovering the entire causal graph requires running intervention on an exponential number of candidate
'mpirical sparsity refers to the observation that most benchmark datascts used for causal discovery evaluation, the augmented hidirectional graph (see Theorem 7 indicating whether nodes each other's Markov blanket p-degenerate defined Definition 4 with small p.

---

Published in Transactions on Machine Learning Research (01/2026)
cause-effect relationships among subsets of variables, thus rendering this approach computationally costly and impractical (Pearl 2009).
To sidestep such expensive interventions, numerous approximation approaches have been developed to instead find an equivalence class of DAGs (Pearl et al. 2016 Peters et al.) 2017]. all of which are compatible with a set of statistical evidence or constraints derived from the observational data? Most of these approaches are formulated as graph searches that proceed either by finding statistical evidence to eliminate incompatible
graph candidates (Spirtes & Glymour 1991 Spirtes 2001 or by optimizing for heuristic score defined on graphs (Hauser & Bithlmann 2012 Rolland et al. 2022 Montagna et al. 2023 Such approaches however are either (a) not scalable due to the expensive computation cost of running statistical tests on an exponentially large number of cause-effect relationship candidates (i.e., subsets of variables); or (b) less accurate due to the heuristic nature of the score function defined over the graph space.
For examples, Spirtes & Glymour (1991) and [Spirtss (2001) required performing local conditional independence lests between eflect variables and their corresponding sets of causes across all candidate causal relationships, incurring a process that scales exponentially with the number of features/ /variaeev In contrast, Hauser & Bühhmann (2012), Rolland et al. (2022). and Montagna et al. (2023) proposed direct optimization of a global hcuristic score on graphs, avoiding the need to solve such constraint satisfaction problems involving an exponentially large set of local constraints. There are also other approaches in this direction which further impose simplified assumptions on the functional structure of the causal relationship (e.g., linear relationship between effects and causes perturbed with Gaussian noises). Based on those assumptions, the problem of causal learning can be formulated as a continuous optimization task and can be solved by more effective solution techniques (Kalainathan et al., 2022 Lachapelle et al., 2019 Ng et al. 2022 2019h. However, these approaches often perform less robustly when the heuristic scores modeling assumptions do not fit well with the nature of the observational data (see Section 5).
To mitigate the aforementioned limitations, we propose a. new solution perspective based on a novel causal test. Such a method avoids imposing additional assumptions on the effect-cause data generation process while also achieving better performance with affordable computation costs. Our approach leverages the invariance of the effect-cause conditional distribution P(effect cause) to changes in the prior cause distribution P(cause). This inspires a principled test for causal relationships via estimating the variance of P(X Z) across syuthetic data augmentations that reflect different cause distributions P(Z). That is. (Z + X) potentially represents a causal relationship if the ieta-induced P(X Z) does not change muct across different choices of P(Z).
This test can then be integrated with a systematic search that identifies candidates for the causal parents of all variables with quadratic complexity, achieving improved scalability and performance over previous methods. Accordingly, our proposal method relies on the contrapositive statement: invariance does not occur, causality does not exist. Furthermore, with mild assumptions, we can recover Pa[X] among the set of plausible candidates by Theorem 1 Our technical contributions that substantiate the above include:
1. An invariance test that reliably determines - PP[] for each effect variable X and candidate set of causes Z via checking the variance of P(X Z) against changes in P(Z): In particular, the invariance test identifies and constructs the most informative data augmentations to reliably approximate the variance of P(X Z) across potential changes to PZz). This is achieved via a downsampling scheme of the observational data that modifies P(causc) without changing the conditional P(effect | cause) for true (cause, (ffusc. tuples (Section 4.2h.
2. A practical findit-iindig algorithm that (i) adopts a previous approach (Edera et al. 2014) to find the Markov blankets of all variables using obscrvational data; and (ii) uses this information to construct an augmented bidirectional graph for each variable whose maximal cliques correspond to its most plausible candidate parent sets: Due to the sparsity of such augmented graphs, the number of maximal cliques (correspondingly, parent candidates Section (.3) is quadratic in the number of variables and can be enumerated by an effective depth-first first search (DFS) algorithm. The developed invariance test (Section 4.2
Finding the true causal graph is not possible without rumning intervention since different graphs can entail the set of statistical constraints induced by observational data. Thee variance of the true effect-cause conditional distribution P(effect | cause) with respect to changes in P(cause) is theoretically zero but in practice, P(effect cause) has be estimated using observational data which causes (small) additional variance due to (slight) variations in its est imation across augmented datasets.
2

---

Published in Transactions on Machine Learning Research (01/2026)
can then be used to find the true parent set among the plausible candidates for each variable, thus recovering the true causal graph.
3. An extensive empirical evaluation of our proposed causal discovery framework on a variety of diverse benchmark datasets including both synthetic and large-scale real-world datasets: The reported results consistently show that our framework performs better than or comparable to prior work in terms of causal discovery performance while achieving generally better scalability, with an average of up to 73.3% reduction in spurious rate and (up to) 25x reduction in processing time (Section 5
2 Related Works
Existing causal discovery methods that aim to recover causal graphs from observational data without using interventional data can be categorized in three main groups (Glymour et al. 2019
First, constraint- based methods aim to recover an equivalence class of causal graphs via deriving statistical evidence from the observational data to eliminate incompatible candidates as much as possible. Their main goal is to minimize the chance of mistaking correlation for causation and hence, maximizing the reliability of the out put graph. For example, Peter-Clark (PC) (Spirtes & Glymour 1991) and Fast Causal Inference (FCI) (Spirtes 2001] Spirtes et al. 2013) run all possible conditional independence tests between all potential (effect, cause) tuples to find the most reliable relationship candidates that pass all tests. In practice, while such approaches often produce reliable results, they do not scale well to high-dimensional datasets since the number of (effect,, cause) candidate tuples often grows exponentially in the number of variables aariables/featurrs /eatuures (Spirtes et al. 2001).
Alternatively, score-based methods instead use heuristic scores defined on graphs to reformulate causal learning as an optimization task which associate true causal graphs with those that maximize the score (Heckerman) let al. 1995 Chickering 2002; Teyssier & Koller 2012 Solus et al.) 2021). Thus, unlike constraint- methods which cast causal learning set of local constraints, soored-based methods recast it as global optimization task, which often admits more scalable solutions. As a result, existing score-based methods (Hauser & Büihlmann 2012 Rolland eet al. 2022 Montagna et al.] (022) often scale better to larger datasets. However, the heuristic design of the score function imposes implicit assumptions on the causal structure of data which are often violated in practice. Consequently, the reliability (i.c., not mistaking correlation for causation) of methods this group is relatively lower than those of constraint-basod methods (Section 5.1 Figures 2 and 3
constraint satisfaction task that involves an exponentially large
To reconcile the above conflicting goals of reliability and scalability, model-based methods adopt additional assumptions on the effect-cause generation model of the observed data. (e.g., linear relationship perturbed with Gaussian noise (Zheng et al. 2018)). This often allows a provable reformulation of the true causal structure as an optimal solution to a cont inuous optimization task that can be solved efficiently with numerous modern and scalable machine learning algorithms (Scanagatta et al. 2015 Zheng et al., 2018 2020). However, their performance often not stable in scenarios where such assumptionn ions on the data generation process do not hold, e.g., non-linear effect-cause relationships perturbed with non- non-Gaussian noises.
Existing Limitations. Overall, existing methods are either limited by expensive processing costs (as seen in constraint-basee approaches), unreliable performance due to the use of heuristic score function and greedy nature of the optimization algorithms (as with score- score-based methods), or inapplicable to scenarios involving unknown data-generation models (as is the case with model-basod approaches). To mitigate these limitations, we investigate an alternative approach to causal learning by exploiting the invariance of the effect-cause conditional distribution across different data augmentations that induce changes to the prior distribution of causes. This inspires highly scalable and effective causal discovery procedure, as detailed in Section 4
3 Problem Formulation and Background
Let D denote a dataset comprising n observations (1) ..... (") of set X (X,, Xx) of d random variables. Supposing these observations are drawn independently from an unknown distribution
3

---

Published in Transactions on Machine Learning Research (01/2026)
Effect-Cause Distributional Invariance Test
Section 4.2.2 Resampling Data
D1
D2
D3
Dm
Causal Graph
PB)
Section 4.2.3 Choosing
P(B)
Z
Source Priors (B)
(B)
(2)
Section 4.2.1 Finding Source Variables
Z
Original Data
Z₁
Part I: Section 4.2
Part Il: Section 4.3
Invariance Test Preparation Plausible Parent Sets
Finding Plausible Parent Sets
Figure 1: Overall workflow of our proposed GLIDE framework which comprises two main sleps: (a) an algorithmic configuration - the key effect-cause distributional invariance test that helps test potential parent-child relationships (Section 4.2 and (b) a graph scarch algorithm exploiting prior knowledge of cach node's Markov blanket and an (empirically verified) sparsily of causal graphs to provably reduce the number of tests to recover the true causal graph (Section 4.3
P(X)= 1(X,, X2...., Xa), want to learn from D DAG G (X, E) over (X1, X2,..., ... Xa) which is causal model of P(X) following Definition
Definition 1 (Causal Model (Peters al. 2017 ) A direct acylic graph G G= (X, E) is causal model of P(X) if every conditional independence derived from P(X) can be derived from G, and vice versa.
A conditional independence X dLp Y Z derived from P(X) means P(X Y, 2, = P(x Z) where P(X YY and P(X Z) are marginal likelihoods derived from P(X). On the other hand, a conditional independence X l Y Z derived from G means given Z, X and Y are d-separated following the below definition d-separation.
Definition 2 (D-separation Pearl 2009)). ) Given two variables X,Y E X and a set of variables Z S X (X, Y), X and Y are d-separated given Z (i.e., X 1 Y Z) if any path between them contains either: (a) a fork A t B + with B E Z, (b) a chain A +
C such that Z; and (c) a collider A B C such that B or any of its descendants Desc[ [B] is not in Z.
Definition limpiiss G is a causal model of P(X) when X _LP Y Z + X Le Y Z. Thus, supposing G a
causal model of P (X), X L (X) Desc[X]] Pa[X] since X LL. (X \ Desc[[]] Pa[X] due to d-separation
see Definition
2 where Pa[X] denotes the set of parents of X in G.
The above means that each node X is independent of its non-descendants \ Desc[X X] given its parents Pa[X] according to P(X), resulting in the following factorization: (I) P(X)= = I P(X; Pa (X:) which implies (IT) X 1 (x MM(X)) Markov blanket of variable X that consists of its immediate parents, its children, and its children's other parents. (II) can also be verified via checking d-separation on G.
Local Markov Conditions.
M(X), where M(X) denotes the
Core Idea. Supposing that Markov blanket condition (II) is known, an augmented bidirectional graph for each variable can be constructed such that cach of its maximal cliques corresponds to plausible candidate Z of the true causal parents Pa[X]. We can leverage this representation to develop an effective DES procedure that systematically enumerates through each candidate with quadratic time complexity in the number of variables (Section 4.3). Each candidate can be tested using the developed causal test via synthetic data augmentation (Section
For this, we adopt an existing algorithm that can recover the Markov blankets for all X e G from obser- vational data D with O(d) complexity (Edera al. (0044). This is guaranteed under the causal sufficient. assumption (Pearl al., (0166) that there no unobserved common causes (i.e., confounders) that affect the causal mechanism generating the observed data.
4

---

Published in Transactions on Machine Learning Research (01/2026)
4 Invariance of P(Effect Cause)
First and foremost, we state our premise regarding the invariance of the causal relationships, which is well aligned with principles of causality discussed in (Peters et al. 2016) and is commonly accepted as an intrinsic property in Bayesian structures:
Observation 1 (Causal Invariance). Let X E X be any variable and B C X be the set of source variable D The causal conditional distribution P(X Pa(X]) remains invariant with respect changes of the joint distribution P(B) over the source variables B.
Remark 1. A direct implication of Observation is that the conditional distribution of any variable conditioned on its parents remains invariant with respect to the change in the joint distribution P(X) if the change is originated from P(B) since P(X)= = P(B) IIxØ P(X Pa(X]) (local Markov Conditions).
Relying on this logic, we design our proposed framework accordingly. To commence with, our overview is presented in Scction 4.1 where we devisc an efficient test to identify the causal parents of any variable. However, utilizing such method requires that we can simulate conditional distributions across different joint distributions P(X). As we will show in Section 4.2 this can be done by sampling environments of distinct. joint distributions of the source variables P(B). For this, we design an approximation search method to locate the source variables, then we devise a sampling strategy that is both theoretically sound and efficient. Finally, we present plausible parent selecting strategy in Section 4.3 The bird's eye view of our workflow is depicted in Figure
4.1 An Overview
The core idea of our causal discovery framework is based on the invariance of the (marginal) effect-cause conditional distribution P(X Pa(X]) induced from P(X) across different choices of the prior P(B) over the source variables B - {B) E X, Pa[B] = 0} C X of the true causal graph G (X, E) see Theorem 1 below (a detailed proof is provided in Appendix A.I).
Theorem 1 (Effect-Cause Distributional Invariance). Let PI P2(B),.... Pm (B) denote set of different priors over B. Let P.(X)= ((B)) P(X P(X B B) denote the corresponding angmentation of the true data distribution P(X) when we replace its marginal prior P(B) with P,( For each variable X E X B, its true causal parents Pa[X], its spouses Sp[X], and plausible candidate set Z = [Z Zz.... ., assuming that Pa[X] n Sp[X] = Sp[X] n B = U, we have:
(B)
3Zj E 2 VP,()) P
[ (2(
0 <> Z, =a=Pa[X]
(1)
(x.)). Here, VP,(x)~ (P: (X] 2,)1 denotes the and P.(X) is drawn from P & (P) PXX)., statistical voriance of the induced Z,) over the random choice of the joint distribution (X)
Intuitively, we argue that the parent set of a variable X is one of those when conditioned on, the probability of observing X becomes unchanged, hence zero statistical variance with respect to distributional changes of the source variables. With mild assumptions (see Appendix A.I], we can assure the uniqueness of such Zj in Equation 1 This result reveals principled lest of whether a subset Z,ε 2 C X is the parent set of X E X according to the true causal graph G. The intuition is that if we have access to D2. Dm such that D: ~ P.))), we can test whether Z, Pa[X] via checking whether the sample variance of the empirical conditional P Zj) is distinguishably small. As we show in Theorem 5 later, the aforementioned {D; y" can be simulated by sampling from the observational data D. Further details on the sampling method is discussed in Section 4.2.2 Given D2...... Dm above, we perform the invariance test as formalized below:
A. Difect-Causs Invariance Test. Given m augmented datasets D1, Dz...., Dm which are resampled from D P(X) such that D; P.(X), then Zj 2 is the causal parent set of X, i.c., Z₁ Pa(X], when
Σ
m
P. x 2,,
P x 2, )} 0, where P(X12))
1 Σ m m
P
x1,)
(2)
4A source variable/node (a.ka. root node, exogeneous node) is node that has no parent nades
5

---

Published in Transactions on Machine Learning Research (01/2026)
Empirically, we consider = 0 equivalent to where 10-3, This test will become more accurate with more source priors. When m is sufficiently large and (P:(B) |= are sufficiently representative of the entire space of source priors, the test might not be perfect but remains highly accurate (see Section 5 In the case multiple 2 E 2 satisfy Equation
2 the one with minimum variance is chosen.
B. High-Level Framework. Supposing that we know how to generate D1. D2...., Pm (see Section 4.2 for which D; ~ (.) as required in Theorem [ solving Equation 2 can be achieved via exhaustively checking all subsets Z as candidates for Pa[X]. Repeating this for all X allows to recover the causal graph. This process is, however, impractical since its complexity is exponential in d.
Fortunately, this can be avoided using the local Markov condition (II) (see Section 3) which implies that (i) Z C M)) if Z = Pa[[ and (ii) there are at most. p = ((d) candidates for Pa[X] which can be provably identified via finding maximum cliques on augmented bidirectional graphs. Under an empirically verified assumption on causal graph sparsity, this can be achieved with a customized DFS procedure with (dd)
complexity (see Theorem 7 Section 4.3). This results in an O(d2) total complexity for our causal discovery framework as detailed below.
C. Time Complexity. As there are d observational variables, our framework needs to make d calls to the parent-finding routine in part (B). As each routine will consider at most p (d)) plausible parent candidates, the per-call complexity is 0(md). This amounts to a total cost. of ((md²) to recover the full causal graph. This can be enabled with (a) a practical (dd) overhead to find the Markov blankets for all variables following the prior work of Edera et al. (2014); and (b) another O + m|D| + m B) overhead construct D1, D2 ..... Pm for Equation (Section 4.2.4 The overall complexity (including the overhead) is therefore (md2 + D mll))
D. Identifiability. Generally, our proposed method does not guarantee a unique solution given solely the observational dataset. However, assuming that the basis set found by Theorem 2 matches the set of source variables, GLIDE can recover the correct causal parent set for all nodes that satisfy Thcorem D assumptions. In practice, we can only estimate the basis set (Theorem 3 and there exist. variables who do not satisfy assumptions in Theorem | Regardless, as shown in Section our performance remains superior over stato-of-tho-art (SOTA) basclines.
To substantiate the above high-level framework, we need to (i) choose representative variants of the source priors to improve the test reliability (Section 4.2 and (ii) find all plausible sets of candidate for Pa[X] for all X E X which are guaranteed to have sizes at most (d) (Section 4.3)
4.2 Augmenting Source Prior
To enable practical use of the parent-finding routine in Equation
2 we need to (i) find the set B of sources, Pm B) so that the invariance test is reliable;
(ii) choose a representative set of source priors Pi (B), and (iii) re-sample D; from D so that D2 ~ (x( = (B) P(X B B).
The above (iii) is essential because we do not have direct access P(X \ B B). This means that we cannot compute (X) and induce Px Z) to assess its variance directly. However, as P(X B B) can be simulated using the original data D, we re-sample D, from D so that D, P,(X), which then allows to use D; to estimate Pa( Z) and consequently, its variance.
For this, we will develop an algorithm to find the source variables B (Section 4.2.1 We will then derive a re-sampling procedure to obtain a downsampled dataset D; from D) such that D, ~ (X) (Section 4.2.2 Last, we will detail a criterion choose informative P.( (for the invariance test) basod on the above re-sampling procedure and how determine a set of representative source priors (Section 4.2.
4.2.1 Finding Source Variables
To find the set of source variables (i.e., those with no parent in the true causal graph), we introduce below the concept of a basis of a DAG.
Definition 3. A basis B C X of DAG G (X, E) is set of mutually d-separated (independent) variables
(see Definition 2 such that for each X & B, there exists X' E B such that X' J X.
6

---

Published in Transactions on Machine Learning Research (01/2026)
As finding the true sources is not possible without intervention, will use basis set as surrogate. This is because the basis has similar properties to the set of sources. First, similar to source variables, basis variables are mutually independent. Second, each source variable X is either in the basis B or shares the same dependence set ((X)= = $(X') - as defined in Theorem 3 with a basis variable X' E B (sce Appendix B.2p. Furthermore, Theorem 2 confirms that the maximum size of a basis set is equal to the number of sources the true causal graph. This means that changing the prior over basis variables will have similar effect changing priors over source variables. Hence, we can use the priors over basis variables instead of priors over sources to test the effect-cause invariance.
Theorem 2. The marimum size of a basis set of a DAG is upper-bounded by the number of its sources.
The proof of Theorem 2 is deferred to Appendix
A.2 A direct corollary of Theorem is that the maximum
sized basis has the same cardinality as the source set of the graph. Theorem
3 provably provides a systematic
method to identify such maximum sized basis set with O(d2) computational complexity. See Appendix AA.3
Theorem 3. Let V = X and (X) A (YEV : Y M X3. A set B is constructed via iteratively (i) inserting into B variable X E V with the lowest (+(X)); ; (t) setting V + V V ((X); and (tii) stopping when ==0. Then, B forms mazimun-sized basis.
Note that X ((X) by default since X H X, VX. The intuition of Theorem is that given a causal graph, the set of dependence nodes for a source does not include other sources while most non-source nodes are dependent on all their upstream (including all source nodes) and downstream nodes. Thus, the set of dependence nodes for source node will not be larger than that of non-source node. Computing (XX requires checking Y A X which can be practically achieved using existing reliable pairwise independence tests provided by the causal-learn open library (Zheng et al. 2024). The independence checking for (Y 1 X) is done by the ((n) Fishor test (Fishor et al. 1921). As there are ((d2) such pairs, the total cost to compute all (X) is Onnd)). Furthermore, (XX does not. need recomputing after step (ii) in Theorem 3 Hence, computing & is a one-time overhead that is not expensive. It is worth noting that the con-uniqueness of the basis set has negligible empirical impact on the overall performance of our proposal. In fact, across multiple random seeds, data generating models, and topologies, the highest recorded variance in term of SHD is less than 6% of the mean value across a variety of experimental scenarios.
4.2.2 Re-sampling Observational Data
Given a target (BB,, we want to find a resampled dataset D; from the original observation data D such that Di ~ P.( To guarantee this, D; must be a downsampled version of D via sampling with no replacement to avoid introducing duplicates and hence, false causal biases into D,. Intuitively, because this operation only removes information rather than generating new samples, it cannot introduce any spurious patterns into the data. To elaborate, downsampling can be viewed as introducing an auxiliary selection variable for each observational snapshot of the underlying DAG. Spurious dependence can arise only if this selection variable conditions on a collider of the form U + X + V where U and V originate from different source ancestrics, in which case the corresponding trail changes from inactive to active. In our construction, however, downsampling uses B = b as the sole selection criterion, so the selection variable depends only on basis variables. Each basis variable is either a source or node with single-source ancestry (see Lemma 2 In the latter case, such basis node can participate in collider U + X +- V only when U and V share the same source ancestry: otherwise, X would depend on multiple sources and, by Lemma 2 could not be selected into the basis. Since U and V are already dependent via their common source, act ivating this collider does not introduce any new dependence. On the other hand, among valid downsamples D, ~ (X(),, we want to choose the one that has the minimal downsampled rate |D] to preserve (as much as possible) obscrvations of the effect-cause conditionals in P(X B B). Interestingly, it can be shown that given the target P. PBB), the minimum downsampling rate can be computed and achieved via. the following results.
Theorem 4. Suppose that D2 is the downsampled dataset with minimum downsampling rate that satisfies the condition D2 ~ P.(X). Then, il follows that.
DD/|Daa = 7, where 7t- = min ==mmin (P(B=
(3)
b
7

---

Published in Transactions on Machine Learning Research (01/2026)
Given the (computable) optimal downsampling rate in Theorem 4 the corresponding downsampling procedure that achieves it can be derived via Theorem below.
Theorem 5. Let % be defined in Theorem Suppose that D; is created via sampling without replacement
P.(B- b) DD- -% "i points from D where B b. Then, DI/ID 1/7 and D; (.()).
The proofs of Theorems and are deferred to Appendices A.4 and A.5 To elaborate more the sampling
of Dr from D, we determine the size of D; via Theorem This allows us to compute (for each candidate value b of the source variables B) how many samples (where B = b) we need to draw (without replacement) from D. Those samples (across all observed candidate values B bb constitute D;. Then, Theorem 5 guarantees that such D, is statistically equivalent to a direct sample from P.. In a practical implementation, we associate P(B b) with its empirical estimate P(B b) R |D = b]l/|D) For continuous data, we categorize the input data into fixed-size bins.
4.2.3 Choosing Source Priors
We can now leverage the insight of Theorem to choose the most informative source priors to enhance the reliability of the invariance test in Equation 2 Intuitively, we want D, & D so % should not be too large. Otherwise, as 7/8 + 1, D; + D, and P,( + P(B) which camnot be used to test the effect-cause invariance. On the other hand, if y, is too small, D2 might drop too much information from which might obscure some effect-cause relationship.
As such, we want to choose P. (B) such that its inverse downsampling rale 74 (see Theorem is above a certain threshold %o where % E (0, is adjustable parameter. Our ablation studies in Table 3 investigate the impact of % on the overall causal discovery performance. The question is now how to sample representative (B( from the subspace of source priors whose (inverse) optimal downsampling rale 7/ 2 7o- To answer this question, we establish Theorem
i below which characterizes its convex hull (see Appendix A.
Theorem 6. Let r A [Dom(B)] denote the number of (categorical) candidate values of B. The subspace of (B( that satisfies 7i 2 "o is a conver subspace Cr(r) of the r-dimensional simplex Dr over Dom(B) which cuts D, at r points (B), p(2) (()) (B)
.t)...,(p))(. (B) representing its convex hull:
(B) = as P(B) + (1 o..) (B) where O% 10-77" )/11) -qx,
(4)
Gh = PPP - b()),, and (B) is a point mass function that assigns 1 when = (b))) and otherwise. Here, (kk) is the k-th candidate value in Dom(B).
PP-BBBbb))
Each source prior PB) with % 2 %o belongs to this convex and can be represented as a linear combination of the above points:
P.(B) = ag (B).
(5)
k=l
where El= Og 1 and ag 0. Hence, P.( can be sampled via drawing (a1, (2,...,,,, u,) from the simplex Dr and using Equation 5
4.2.4 Practical Implementation
The parent-finding procedure involves three main steps: finding basis variables, sampling source priors, and re-sampling observational data. First, finding the basis variables has a time complexity of Odd). Next, we sample a from a Dirichlet distribution to compute P,(B), with 10a samples clustered using K K-means. The K = m centroids are selected as source priors, with complexity of O(mr) O (Bl) where e is the maximum number of candidate values of single variable. To avoid exponential computation costs in [B,, we resample for each basis variable individually, reducing the complexity to O(mc Bl). Finally, generating the augmented datasct for cach (BB incurs a cost of O(m[Dl) via Theorem 5 This amounts to the ( + me| B| + m DI) total time complexity. As for hyper- parameters, we recommend choosing m as high as the incurred additional runtime is acceptable. As for Y, we find 7 = 0.5 is an empirically reasonable selection. As a rule of thumb, 7 should be high enough so that estimations induced from downsampled data are adequately accurate, but also should be small enough to promote diversity between downsampled datasets.
8

---

Published in Transactions on Machine Learning Research (01/2026)
4.3 Finding Plausible Parent Sets
To ensure that the parent-finding routine is scalable, we customize DFS algorithm which leverages prior work (Edera et al. 2014) in Markov blanket identification to provably find O(d) candidate parent sets for each effect variable X. The true causal graph can then be recovered via solving Equation L for each X with respect to the discrete set. of O(d) plausible parents found above.
As previous work in Markov blanket identification incurs an O(d2) cost (Edera et al. 2014). the total cost of our parent finding phase also O(d)) This helps avoid the bruto-force search through all subsets of variables candidates for Pa[X] whose complexity is otherwise exponential in d (Peters et al. (016). To achieve this, we leverage the following result which relates candidate sets of causal parents to maximal cliques on augmented bidircctional graph.
Theorem 7 (Plausible Parent Sets). Let M(X) denole the Markow blankel of X and (X) - M SS[[X]
denote the set of non-spouse variables within MX).. For each variable X, let GX)) (V,E") denote bidirectional graph where V (X) and (U.V) E E iff V E MI) and U E MVV. Then, Pa[X] S (X)
belongs to at least clique in (X)
Here, a clique is defined to be set of nodes in the bidirectional graph G' such that there is an edge between any two nodes. See Appendix A.7 for a detailed proof. Notice that M' (X) can be derived directly from M(X by using conditional independence tests. Indeed, since AY E Sp Y IL X Pa[X] and Pa[X] C M(X), we are guaranteed to detect and remove all such Y E M.)),, hence M (X). Finding MI (X) for any X E X costs a computational complexity of ((2M) where M < is is the maximum number of elements in a Markov blanket. Leveraging Theorem 7 we can find all plausible parent sets of X as subsets M (X), each corresponding to a clique in G'(X). Restricting the set of plausible parent sets that of maximal cliques in G', we can reduce the task of finding plausible parent sets to finding maximal cliques in a bidirectional graph. For a sparse graph (*) with a low degeneracy constant (scc Definition 4) this can be achieved effectively via a customized version of the DFS-based Bron- Kerboshh (Bron & Kerbosch 1973) algorithm (see Appendix B.3 which finds all maximal eliques with (dp3p/3) time complexity.
Definition 4 (Degeneracy [Kirousis & Thilikos 1996)). A bidirectional graph is p-degenerate if every subgraph has at least one node with degree less than or equal to p. The graph's degeneracy is the smallest p for which it is p-degenerate.
Considering p as small constant compared to d, the cost of finding all plausible parent sets for each variable X is effectively (d). Hence, the total cost of finding all plausible parent sets for all variables is (dd). Furthermore, it is also established in Bron & Kerbosch (1973) that the (worst-case) number of maximum cliques is O(dd - p) 32/3). Our empirical studies in Appendix C.4.3 show that S 13 across all benchmark datascts even for the largest graphs, confirming that the set of maximal cliques in G'(X) is practically Odd)
with a constant factor smaller than 313/3 < 117.
5 Experiment
This section evaluates and compares the performance of our proposed method, Causal Graph Learning via Distributional Invariance of Cause-Effect Relationship (GLIDE), against existing state-of-the-art (SOTA) baselines in causal graph learning.
Basclines. Our empirical evaluations are conducted with respect to a variety of data models and diverse
suite of both classical and recent baselines (Appendix C. Table [], including PC (Spirtes & Glymour
1991h, GIES (Hauser & Bühlmann 2012), FCI (Spirtes 2001). NOTEARS (Zheng et al. 2018). MLP.
NOTEARS (Zheng et al. 2020), DAS (Montagna et al. 2023) and SCORE (Rolland et al. 2022
-Each
baseline is configured with its best hyper-parametery (see Appendix C.I.
Datasets. Our experiments are based on both synthetic and real-world datasets. For synthetic experiments (Sections 5.1 and Appendix C.2), we follow the commonly used protocols in previous work to generate
observationals data based on the Erdos Erdos-Renyi (Zheng et al. 2018 Heinze-Dema et al., 22018), bipartite and
scale-free (Zheng et al. 2020) classes of causal graphs. Due to limited space, we only present the empirical
9

---

Published in Transactions on Machine Learning Research (01/2026)
Notears
MLPNotears
DAS
SCORE
GIES
FCI
GLIDE
L-G
nL-nG
L-G
nL-nG
20
10
SHD (x10')
SHD (x10p)
6
Lx103)
15
SHD (x10a)
8
2
2
SHD
0
100 200 300 400 500
100 200 300 400 500
600 700 800 900 1000
600 700 800 9001000
Number nodes
Number of nodes
Number of edges
Number of edges
L-G
Spurious Rate (%6
nL-nG
Spurlous Rate (%)
L-G
Spurious Rate (%)
nL-nG
75
Spurious Rate
30
50
20:
20
25
25
10
100 200 300 400 500
100 200 300 400 500
600 700 800 9001000
600 700 800 900
Number nodes
Number of nodes
Number of edges
Number of edges
Runtime (x10² mins)
L-G
Runtime (x100² mins)
nL-nG
Runtime Lx102 mins)
L-G
Runtime (x102 mins)
nL-nG
& S
20
15
w
2
2
100
1-
08
100 200 300 400 500
100 200 300 400 500 800 700 800 9001000 Number of edges
600 700 800 9001000
Number nodes
Number of nodes
Number of edges
Figure 2: Baseline performance in normal setting (continuous data). Lower metrics are better.
Figure 3: Baseline performance in extreme setting
(continuous data). Lower metrics are better.
results on synthetic data generated from the Brdos-Renyi class of graphs. Other experiments with the bipartite and scale-free classes of graphs are deferred to Appendix C.44 We also provide evaluation on the SACHS
dataset (Section 5.2 provided by [aahs et al. (005).. This dataset has been widely used in many existing literatures on causal discovery Zheng al. 2018 2020 The SACHS datasel (11 variables, 17 edges) contains records of continuous measurements of expression levels of proteins and phospholipids in human immune system cells and is widely accepted by the biology community. Additionally, we also conduct experiments real-world graphs (see Appendix C.3 using the datasets provided in the bulearn package (Scutari 2009 which includes 7 small to large real-world graphs. Notably, the Munin dataset features a large-scale graph comprising 1041 variables.
Evaluation Metrics. We use the structural Hamming distance (SHD), spurious rate (percentage), and running time (minutes) to measure both the effectiveness and scalability our proposed method. The SHD and rumning time metrics are commonly used in the evaluation of existing causal discovery methods (Zheng et al. 2018 Montagna et al. 2023). The spurious rate measures the ratio of the number of false causal re relationships among all causal relationships discovered by each algorithm. This helps compare the reliability of causal prediction across baselines. All performance metrics are reported with their mean and confidence interval 95% averaged over 10 independent runs.
5.1 Synthetic Continuous Data
We consider two experiment settings. First, in the normal setting, the graph's number of edges is the same its number of nodes, which varies from 100 to 500. Second, in the extreme setting, fix the number of nodes at 500 and increase the number of edges from 600 to 1000. In each setting, we further consider two dala-generation models for the elfect-cause relationship: (1) linear relat relationship perturbed with Gaussian noise (annotated as L-G); and (2) non-liness and non- non-Gaussian relationship (annotated as nL-nG). Such data generation mochanism implemented using the public code (Zheng et al. 2018 2020). In each exporiment, we first generate random Erdos- Erdos-Renyi DAGs and then simulate (synthetic) observational data from it using the above mechanisms. Our reported results and observations in each setting detailed below.
10

---

Published in Transactions on Machine Learning Research (01/2026)
Normal Setting. Figure 2 reports the performance of our method GLIDE in comparison those of other baselines. In the L-G cases, GLIDE, NOTEARS, and MIP-NOTEARS outperform the rest of the baselines significantly in terms of both SHD and spurious rate. Notably, GLIDE achieves remarkable 64.17% reduction rate on SHD over FCI. Furthermore, GLIDE incurs much less computational cost than both NOTEARS and MLP-NOTEARS (see the runtime plots) while achieving comparable or better performance. For example, GLIDE runs 96.54x and 15.66x faster than both NOTEARS and MLP-NOTEARS in 100- and 500-node graphs while being second to MLP-NOTEARS in terms of SHD (with small gap of 16.2%) and achieving best spurious rate in the largest graph setting with 500 nodes (4.2% vs the second best of 11.75% of NOTEARS and third best of 13. 19% of MLP-NOTEARS) In the nI-nG cases, we also have a similar observations where GLIDE is again the second best in SHD and best in spurious rate (in the largest graph settings) while being much more scalable than the best and second best baselines in SHD and spurious rate, respectively. Our reported results also confirm our intuition earlier that score- based methods, despite being more scalable than constraint- basod approaches, tend to perform much less robustly in large graph settings where the score function becomes less accurate. For example, both SCORE and DAS reach over 30% of spurious rate in 500-node graphs.
Extreme Setting. Figure setting. In the L-G cases, it is observed that GLIDE achieves the best performance in both SHD and spurious rate while also achieving the fastest running time in most graph settings. In terms of SHD, the performance gaps between GLIDE and the second best (NOTEARS) and worst baselines (SCORE) are 11.74% and 108.5%, respectively. In terms of spurious rate, GLIDE also improves over the basclines with substantial performance gaps ranging from 5.46% and 48.61% % against the second best and worst baselines, rospectively. In the nI-nG cases, GLIDE and MLP-NOTEARS perform comparably as best basclines in SHD but GLIDE outperforms it by gap of 8.23% in term of spurious rate. Furthermore, against the most high-performing baselines in this case, NOTEARS and MLP-NOTEARS, GLIDE achieve 9.6x and 25.52x faster processing time, respectively.
3 reports the performance of GLIDE in comparison to other basclines in this
5.2 Real-World Data
The results show that GLIDE achieves the lowest SHD (8.7 = 0.48) and spurious rate (0%), significantly outperforms the second-best baseline, SCORE, whose SHD is 13.2 + 0.66 and spurious rate of 22 + 4%. It is worth noting that GLIDE only takes 10.46 + 0.71 seconds to accomplish this task, which is three times faster than SCORE (at 30.2 + 3.98 secouds). FCI and DAS respectively stand at the third and forth place with the average SHD being 17 and 22.27 + 1.65, and the corresponding spurions rates are 44% and 55% + 4%.
6 Conclusion
This paper presents a new perspective of causal learning via a new invariance test for causality that inspires a reliable and scalable algorithm for recovering causal graphs from observational data. Our approach explores a key insight that the effect-cause conditional distribution remains invariant under changes in the prior cause distribution, leading to a parent-finding procedure for each variable via syuthetic data-auggentttion. This procedure is further coupled with an effective search algorithm that exploits prior knowledge of each effect variable's s Markov blanket and the sparsity of the causal graphs to significantly reduces the overall complexity. This, in turn, results in significant reduction in complexity and marked improvements both in speed and accuracy over existing SOTA methods on various benchmark datasets.
Impact Statement
Our work focuses on improving the scalability and reliability of the causal discovery procedure. The findings in this work offer a promising direction for further research and practical applications in causal inference, especially in scenarios involving large and complex datascts. As such, could have significant. broader impact on future research due to the ever-growing ubiquity of larger datasets. The source code and the datasets used for evaluation are publicly accessible, thus raising ethical concerns.
11

---

Published in Transactions on Machine Learning Research (01/2026)
Acknowledgments
MS was supported by JST ASPIRE Grant Number JPMJAP2405. PIN was funded by Hanoi University of Science and Technology (HUST) under grant number T2024-TD-002. In addition, this work utilized GPU compute resource at SDSC and ACES through allocation CIS230391 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services and Support (ACCESS) program (Boerner et al. 2023 which is supported by U.S. National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. The authors would also like to thank anonymous reviewers for their careful reviews and insightful comments.
References
Timothy J. Boerner, Stephen Deems, Thomas R. Furlani, Shelley L. Knuth, and John Towns. Access: Advancing innovation: Nsf's advanced cyberinfrastructure coordination ecosystem: Services & support. In Practice and Experience in Advanced Research Compuling 2023: Computing for the Common Good, PEARC '23, pp. 173-176, New York, NY, USA, 2023. Association for Computing Machinery. ISBN 9781450399852. doi: 10.1145/35500 3556955133597538
URL https:///oi.oo. 11110.114//366955/ 3597559
Coen Bron and Joep Kerbosch. Algorithm 457: finding all cliques of undirected graph. Communications of the ACM, 16(9):575 577, 1973.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3(Nov):507 554, 2002.
Alejandro Edera, Yanela Strappa, and Facundo Bromberg. The grow-shrink strategy for learning markov network structures constrained by context-spocific independenoes. In Advances in Artificial Intelligence- IBERAMIA 2014: 14th AbericAmeriaan Conference on AI, Santiago de Chile, Chile, November 24-27, 2014, Proceedings 14, pp. 283- 13-2994.2 Springer, 2014.
Ronald Aylmer Fisher et al. 014: On the' probable error' of a coefficient of correlation deduced from a small sample. 1921.
Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery methods based on graphical models. Frontiers in genetics, 10:524, 2019.
Siyuan Guo, Viktor Tóth, Bernhard Schölkopf, and Ferene Huszár. Causal de finetti: On the identification of invariant causal structure in exchangeable data. Advances in Neural Information Processing Systems, 36, 2024.
Alain Hauser and Peter Bühlmann. Characterization and greedy learning of interventional markov equivalence classes of directed acyelic graphs. The Journal of Machine Learning Research, 13(1) (((1))2409-24444 2012.
David Heckerman, Dan Geiger, and David M Chickering. Learning bayesian networks: The combination of knowledge and statistical data. Machine learning, 20:197-243, 1995.
Christina Heinze- Dema, Jonas Peters, and Nicolai Meinshausen. Invariant causal prediction for nonlinear models. Journal of Causal Inference, 6(2):20170016, 2018.
Diviyan Kalainathan, Olivier Goudet, and Ritik Dutta. Causal discovery toolbox: Uncovering causal relationships in python. Journal of Machine Learning Research, (1()))11-5, 2020.
Diviyan Kalainathan, Olivier Goudet, Isabelle Guyon, David Lopez-Paz, and Michèle Sebag. Structural agnostic modeling: Adversarial learning of causal graphs. Journal of Machine Learning Research, 23(219): 1 1-62, 2022.
Lefteris M Kirousis and Dimitris M Thilikos. The linkage of graph. SIAM Journal on Computing, 25(3): 626-66 647, 1996.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-basse neural dag learning. arxio preprint arNiv: 1906 1900.9000000 006.02226, 2019.
12

---

Published in Transactions on Machine Learning Research (01/2026)
Francesco Montagna, Nicoletta Noceti, Lorenzo Rosasco, Kun Zhang, and Francesco Locatello. Scalable causal discovery with score matching. In Conference on Causal Learning and Reasoning, pp- 752-771. The Proceedings of Machine Learning Research, 2023.
Ignavier Ng, Shengyu Zhu, Zhitang Chen, and Zhuangyan Fang. A graph autoencoder approach to causal structure learning. arXiv preprint arXiv:1911.074.30 2019.
Ignavier Ng, Shengyu Zhu, Zhuangyan Fang, Haoyang Li, Zhitang Chen, and Jun Wang. Masked gradient- based causal structure learning. In Proceedings of the 2022 SIAM International Conference on Data Mining (SDM), pp. 424-432. SIAM, 2022.
Judea Pearl. Causality. Cambridge university press, 2009.
Judea Pearl, Madelyn Glymour, and Nicholas P Jewell. Causal inference in statistics: A primer. John Wiley & Sons, 2016.
Jonas Peters, Peter Bühlmann, and Nicolai Meinshauson. Causal inference by using invariant prediction: identification and confidence intervals. Journal of the Royal Stulistical Sociely Series B: Statislical Methodology, 78(5):94 1012, 2016.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
Alexander Reisach, Christof Seiler, and Sebastian Weichwald. Beware the simulated dag! causal discovery benchmarks may be easy to game. Advances in Neural Information Processing Systems, 34:27772-27784, 2021.
Alexander Reisach, Myriam Tami, Christof Seiler, Antoine Chambaz, and Sebastian Weichwald. A scale- invariant sorting criterion to find a causal order in additive noise models. Advances in Neural Information Processing Systems, 96:785-807, 2023.
Paul Rolland, Volkan Cevher, Matthäus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score matching enables causal discovery of nonlinear additive noise models. In International Conference Machine Learning, pp. Research, 2022.
The Proceedings of Machine Learning
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 (08(5)2)):523 529, 2005. doi: 00.112///cience.110550... URL https:/ http://wwwwiiicccce/// www.a/iinncc..gg///////////////a http:/10..111///ccccccc 1105809
Mauro Scanagatta, Cassio P de Campos, Giorgio Corani, and Marco Zaffalon. Learning bayesian networks with thousands of variables. Advances in neural information processing systems, 28, 2015.
Marco Scutari. Learning bayesian networks with the bulearn package. arXiv preprint arXiv:090s.3817, 2009.
Liam Solus, Yuhao Wang, and Caroline Uhler. Consistency guarantees for greedy permutation-based causal inference algorithms. Biometrika, 108(4):795-814 -14, 2021.
Peter Spirtes. An anytime algorithm for causal inference. International Workshop on Artificial Intelligence and Statistics, pp. 278-285. The Proccedings of Machine Learning Rescarch, 2001.
Peter Spirtes and Clark Glymour. An algorithm for fast recovery of sparse causal graphs. Social science compuler review, 9(1):62 72, 1991.
Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, prediction, and search. MIT press, 2001.
Peter Spirtes, Christopher Meek, and Thomas S Richardson. Causal inference in the presence of latent variables and selection bias. arXiv preprint arXio: 132:::02.4400000 2013.
13

---

Published in Transactions on Machine Learning Research (01/2026)
Mare Teyssier and Daphne Koller. Ordering-based search: A simple and effective algorithm for learning bayesian networks, Xrviv preprinl ar Xrx:::1077 1429, 2012.
Alessio Zanga, Elif Ozkirimli, and Fabio Stella. A survey on causal discovery: Theory and practice. International Journal of Approzimate Reasoning, 151:101-129, 2022.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Fric P Xing. Dags with no tears: Continuous optimization for structure learning. Advances in neural information processing systems, 31, 2018
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning sparse nonparametric dags. In International Conference on Artificial Intelligence and Statistics, 3414 3425. The Proceedings of Machine Learning Research, 2020.
Yujia Zheng, Biwei Huang, Wei Chen, Joseph Ramsey, Mingming Gong, Ruichu Cai, Shohei Shimizu, Peter Spirtes, and Kun Zhang, Causal-learn: Causal discovery in python. Journal of Machine Learning Research, 25(60):1-8 8, 2024.
14

---

Published in Transactions on Machine Learning Research (01/2026)
A PROOFS
A.1 Proof of Theorem 1
First, we establish the following auxiliary result which will be leveraged to establish our main result.
Lemma 1. VZ C X B: (Z =) B). Intuitively, this means the distribution of any set of non-source variables Z conditioned on source variables B remains invariant, i.e,, P. B) P(Z B), when we perturb the source marginals, i.e., replacing P(B) with P.(B). As direct consequence, we have P.( f (Z (B [ P(z ()),(B)dB.
Proof of Lemma + We have P(X)= = P(B) B B) where P(X B B = =xxxx P(X Pa(X)) is constant. Therefore: P.(Z B) (BU B B)dB [X((BUZ) P(X B) = P(Z B). ]
Key result: For each node X, its candidate parent set is therefore a power set of Pa(X) U Ch(X) (see Theorem 7) Our key result then establishes that over the random choice (B as prior over sources B: (1) the variance over the induced conditional P.(X Pa(X)) zero, (P) Pa(X))] O; and (2) the variance over Pa(X Z) where Z S Pa(X) Ch(X) and Z + Pa(X) is not This means
use
(P:() Z)] where Z S U Ch(X) to determine the true causal parents of X. The proof is following:
Since we only perform our invariance test on non-source variables X, we consider the following tuple (X, Z)
where Z1 U 22, Z Z, = 0, Z B 0, and Z2 C B. Then we have:
P,(X, Z1, Σ2) iZ Z2) JB P(Z B)P Z2)
JA.Z, P(x, Z1 B)) \ Z2)
( Z) P 1 22)
(6)
Here, the last equality follows from Lemma This mcans if P.) [X 2)] = 0, then 2. P.) must
be a constant function aXX, , Z1, Z2) with respect to P. Given its expansion above, we must have:
P. Z1 B) = 2) B) (X B) B) (X Z1 Zz, B Z2) (X, Z1, Z2),
(7)
where a(X, 11 Z2) a scalar function independent of the choice P(B for source prior. Equivalently, this mcans knowing Zz does not change the condition of X on (Z1, Z2). Following the definition of conditional independence, we have X [ Z2) (Z) Z2). Hence,
VP(B) (P Z)] 0 X T B 22) (Z1, Zz) On the other hand, if X 1 (B) Z2) Z,, Zz) holds, P(X (B) 22), Z1, Z2) P(X Z1, Z2). Plugging
(8)
this into Equation 0 results in:
P((x Z)
P.(X Z1, Zz)
JBYZ, P(X, Zi B)).(B)d \ 22)
Apply Lemma 1 for P(X, Zi B)
JBYZ. P(Z, B)). Z2)
SMYZ, P(X Zi, 22, B (2) (BP B Z2)
SmZ, P(Z, B)P,(B)d Z2)
JB, P(X Z1, (Z, ()) Zz)
Conditional independence holds
JBZZ P(Z1 (B)d( Z2)
SBIZ: P(Z, (B) dd B Zz)
P(X Z1. Zz)
JANZ, P(Z, B) Zz)
P(X Zi, Zz)
which does not depend P2 (B). This means:
X B Z2) (Z,, 22) VP,(B) 2)) =
(9)
15

---

Published in Transactions on Machine Learning Research (01/2026)
Combining Equations 8 and 9
we obtain bidirectional equivalence:
P ( Z1, =2] X (B) Z2) 11, Zz.
(10)
Assuming that Pa(X) n Spouse X) 0 and Spouse( X) n B = 0, the above bidirectional equivalence in Equation 10 has these direct consequences:
21 Z2 Pa(X) + X (B) Z2) (Z) 22) = (P =
=0
b. Z 22 (a(X) U (Y) where Y Ch(X), then when conditioned Y, there exists active
(backdoor) path B
+ Spouse ()) Y + X. Therefore, X 1 (B \ Z2) (Z) U Z2)
does not hold. Hence, VP, 2)] > 0 because otherwise, (P)( 2)] 0 implies
X [B] 22) (2) U Z2) which contradicts the above backdoor consequence.
c. Pa(X) 2 Z J 22 also implies the existence of an active backdoor connecting B and X. Using the same argument as in (b), we know this case also implies (P.(B Z)] > O.
Overall, (b) and (c) colloctively ensure that long as Z is not Pa(X), (R Z)] > 0. Otherwise, (a) ensures that PP,, [P:(X Z)] = O.
We also note that the structural assumptions Pa(X) n Spouse = 0 and Spouse(XX) n B = hold with relatively high probability. For example, regarding real-world graphs (provided by bnlearn (Scutari 2009) we presented, the percentage of nodes satisfying these assumptions is over 75%.
A.2 Proof of Theorem 2
Let S(G) C X and ( S X denote respectively the set of sources and arbitrary basis set of the DAG
G = (X, E). We will prove the following inductive statement:
P(n): "Given a DAG G of n (n 2 1) nodes, we have [BG) < [S()))'
To see this, note that the base case n = 1, G has single node which is obviously both a source and basis node. Hence, [BG [S(G) and P11) is true.
Now, suppose P(n) is true, we will complete induction by showing that P(n + 1) is also true. To show this, let X be terminal node in G' that has no children. Removing X from G thus results in another DAG G' with n nodes.
Let B(G") denote an arbitrary basis set of G.. By definition, nodes in B(G') are mutually independent and for any nodes X' & X (B(G)) U (X)), there exists a node in B(G") such that X' y Z following the definition of d-separation.
Choosing X' E Pa[X], there must exist Z E B(G)) such that Z is connected to X' via d-separation. Since X' is parent of X, Z is also connected to X via d-separation. This means each node in G' is either part of B(G") or connected to another node in ((G") via d-separation. Hence, B(G') is also a basis set of G. In this case, [B(G)] = |S(G)) < [S(G)) since P(n) is true.
Otherwise, suppose (G) a basis set of that is not in G'. In this case, B(G) must contain X and by definition of basis (see Definition 3) and its choice, X must be an isolated node with no parents and children. This also means B(G) \ {X} is basis set of G². Otherwise, there exists node in G that is not connected to any nodes in the basis set BGG,, resulting in contradiction. As such, we have [BG] [BGG) {X] +1 [S(G') + SSG)) where the first inequality is true due to P(n) and the last equality is true since X is an isolated node which is also a souree node.
As such, B(G) < [S(G) meaning P(n + 1) is true P(n) is true. Induction completes.
It worth noting that for every graph G, there exists basis sot B such that [B(G) = |S(G)). The trivial solution = since the set of sources satisfies Definition 3 of a basis set (i.c., source variables are mutually independent and every other variable is dependent on least one source).
16

---

Published in Transactions on Machine Learning Research (01/2026)
A.3 Proof of Theorem 3
To prove that the procedure in Theorem produces basis set with maximum size, it suffices to prove that each step in this procedure removes exactly one source from the graph (see Lemma 2) In this case, the returned basis set has the same size as the source set, which is also the maximum basis set according to Theorem 2
Lemma 2. If X has the smallest number of dependent nodes, i.e,, X arg (+(Y), then X is dependent eractly one source in G = (V, E).
Proof of Lemma 2 Let statement (P) be "X E V has the smallest number of dependent nodes in G", and (Q) be "X is dependent exactly one source in G To prove (P) = (Q), we will prove that its negation (P) A -Q) is false. Specifically, (P) A (Q) reads "X has the smallest number of dependent nodes and X is dependont on two or more sources in G
To see this, when -QQ) is true, there must exist two source variables S₁ and S2 which are dependent on X (note that S and S2 are independent by definition). This means there must exist an active path from Si to X, and another from S2 to X. Consequently, VY ((Si), there exists an active path from Si to Y and hence, X and Y are connected by active path via common cause trail at S1 which implies Y E ((X). Likewise, Y E ((S2) also implies Y E ((X). This means ((X)] 2 (*())) U (S2)). Finally, note that (1) S2 & ((S)) as both are sources; and (2) S2 E ((22) by definition since S2 M. S2. Given this, (()) U (S2)) 2 ((()) + 1 [s,)) Combining this with the earlier established fact that [*(X)] 2 ((S)) (S2) ((2)) leads to (X)) [S:)] which is a contradiction to (P). Thus (P) A -(Q) false.
A.4 Proof of Theorem 4
Let D; be any downsampled dataset of D. Consider candidate value b of shuree variable B that occurs D. Let n(b) and (b(b) denote the corresponding numbers of data points in D and D, that have B = b. Since D, is downsampled from D, nb( I n(b). Hence, let _ DI and n, - (ril
P(B b)/PB) b) z (m(b)/nn/(/(())///// (b)b)),(bb))))).. (n//n) 2 n/nn DD/I/DD
(11)
which follows from the facts that P(B = b) (bb/nn P3 = b) = n. (b)/n; and n(b) 2 n;(b) due to the downsampled nature of Dj. Taking the minimum over all candidate values b of B,
min (P(B= 2 D
(12)
As a result, for all downsampled dataset D; of D, the downsampling rate DD/ID % where % defined in terms of the original and target marginal over the source variable B stated in Eq. 12). Hence, if D, is the downsampled data with minimum downsampling rate |D/|Da 2 7-
A.5 Proof of Theorem 5
Since D, is created via sampling with no replacement P,(B = b)- ID] % data points from D where B = b,
we know that n;(b)= = P.(B = b) [ |DI P: = b) nn 7.. Thus,
Em(b) = Σ Pi(B= =b) IDI b
n1
b
(B = b) n
n.n Σ P3 = b) (B-(B=)) b
n n
(13)
As result, |D.I/ID ni/n- = 7 and n. (b)/n, (P. = b) (n) P.((). Thus, D, P.( and the downsampling rate DD achieves minimum expected.
17

---

Published in Transactions on Machine Learning Research (01/2026)
A.6 Proof of Theorem 6
A. Convexity. We will first prove that the subspace of P.( that salislies vi 2 is convex. To see this, consider (B) and (B) whose (inverse) downsampling rates (see Theorem 4 7 and 7i are both larger than %/o This means both (B) and (B) belong to the aforementioned subspace.
Now, let a E 10 and P = a. (B) + (1- PB()). The (inverse) downsampling rate of PP is defined
P(B b) P3 (B]=))] b))
)
P(B = b)
71
min
min
a. b) + (1 (1- BBBB
(14)
P(B= =8=b)
1
min
%o
(15)
0.-70 P(B -0) PP =)
70
where the inequality follows from the facts that (1) P(B P) b) 2 2 %o: and P(B bB
%o which follows from the (inverse) downsampling rate's definition in Theorem 4 This means ( b) %o- P(B = b) and PB b) < Tot P(B = b), which can plugged into Eq. 14] to arrive at Eq. (1). This in turn implies (B) belongs to the subspace of source priors with the induced (inverse) downsampling rate above %o- Thus, following the definition of convexity, we know that this subspace is convex.
b) > 1 2
B. Convex Hull. Note that any source prior PB) is a point in an r-dimensional simplex where is the number of candidate values of the source variables B. The convex hull of the aforementioned convex set then be determined by finding its intersection with the edges of the simplex. This comprises points of the following form:
p(k) = ax P(BB +
(1 - as) (B)
(16)
where og is selected such that following equation:
7/ (k) mins p(k) [ = b)) = 7o- Solving for as results in the
P P B bb())
P B b(K)
%o = min
(17)
P B b(() p(B=b(k)
+ (1 (1-ax)
O6 P PB= b(k) + (1-08)
where (8) denote the k-th candidate value of B. The second equality in the above holds since aj E (0, which means 0% P(B b) S 0% P(B b) + 0% equivalently, 1/ou 2 P(B b)/(ax P(B b) ag).
This allows us to get rid of the min operator in Eq. [17] and consequently, compute ag in closed-form:
O = (1 (1-70- P(B = B(x)
(1-(BB-BB b²)) --))bb))
(18)
This completes our derivation for the convex hull in Theorem b
A.7 Proof of Theorem 7
Suppose U, V Pa[X], must follow that U E M(V) and V M(U since U and V are spouses. This means (U,V) E E. Hence, there is bidirectional edge in G'(X) between any parents of X which means there exists at least clique in (X) containing Pa[X].
J
B ALGORITHMS & ANALYSIS
B.1 Pseudo-code for GLIDE
B.2 Find the Basis of graph
This section provides visualization and pseudocode of the basis finding procedure in Theorem
18

---

Published in Transactions on Machine Learning Research (01/2026)
Algorithm 1 GLIDE
Input: Data D, the set of all variables X, down-sampling ratio threshold 7, number of down-sampling sub-datascts m.
Output: DAG G.
1: Find the basis set B using Algorithm
Compute m priors (B) by Theorem
Sample m sub-datasets (D, P.(B) from D (Section 4.2.2
Use the GSMB algorithm to recover the Markov blanket for each node M(X) AA X. Find the plausible parent set 2 for each variable X using Algorithm 3
5:
6: Find causal parent Pa[X] for each X E X by Theorem using {D. its plausible parent set Z.
7: Build G using (X, Pa(]) for each X E X.
8: Return G.
Algorithm Maximun-Sized Basis Search
Input: Data D and the set of all variable indices V. Maximum-Sizza Basis B.
Output:
Compute { ;j II(X; + Xj) using existing statistical independence test.
while V 0 do
Choose = arg minjev (
4: Update B = B U (X,
5: Updale V + V { 1 1 1}
6: end while
7: Return
A. Pseudocode. First, the pseudocode of the basis finding algorithm is detailed below.
This algorithm is guaranteed to find the maximum-sized of basis variables since each of its iteration will remove exactly one source node from V as proved in Appendix A.3 Thus, the number of nodes added to the basis is also the number of source nodes. Such basis is guaranteed to has maximum size following Theorem 2 A visualization of this algorithm provided in Figure
B. Time Complexity. Algorithm 2 consists of (i) computing the dependence network matrix + and (ii) selecting the node with minimum number of dependents from V per iteration. The complexity for step (i) is O(d²) WHILE the complexity for step (ii) is O(d) per iteration. As there are at most d iterations, the total cost of Algorithm 2 is therefore (dd))
B.3 Find Plausible Parent Sets
Theorem implies that if two nodes in the Markov blanket of X do not belong to each other's Markov blanket, then they can not both be parents of X. According to this logic, we wish to find all sets nodes whose Markov blanket contains one another. To this, we develop a tree-based recursive algorithm that works as follows. Given the Markov blanket of all nodes, we perform Algorithm 3 which modification of the Bron-Kerbosch algorithm (Bron & Kerbosch, 1973). Instead of running Bron-Kerbose for each node X and its Markov blanket M(XX), we add virtual node that directly connects to all other nodes in the graph and start building the set of plausible parent sets for all nodes recursively from this virtual node.
Note that Algorithm 3 does not create a sub-tree (branch) if the search space of the root of that sub-tree is guaranteed to re-produce the content of a previously visited leaf and generate no new information see lines 4-7. Otherwise, plausible parent of is the of variables the path from the search node containing X leaf node, which is computed following the recursion of Algorithm g step-by-step visualization of this algorithm is illustrated in Figure 5
19

---

Published in Transactions on Machine Learning Research (01/2026)
X1
X3
X1
X3
X1
X3
5
5
6
5
6
X2
X4
X5
X2
X4
X5
X2
X4
X5
X6
X8
X6
X8
X6
X8
8
8
8
8
X7
X7
X7
8
8
8
B X1
B E
X1, X3
X1, X4
X1, X5
X2 B E
X2, X3
X2, X4
X2, X5
Figure A step-by-step visualization of Algorithm From left to right: each step finds the node with minimum number of dependents and remove along with its dependent set from the graph. Each node is accompanied by number in red that indicates the number of variables being dependent that node, The removed nodes are colored gray WHILE active nodes are blue. The example graph is taken from the ASIA dataset (Scutari 2009).
Algorithm Finding Plausible Parent Sets
Inputs: root node тo containing the virtual variable X., Markov blanket M(X) of all variable X. Outputs: a plausible-pareen-sed tree, the path from each intermediate search node of variable X each leaf node represents plausible candidate for Pa[X]
1: Initialize an empty list of leaves L < 0;
2: search(r.) + the set all variables X;
3: Procedure Recursive-build(7) for do
if path(r, + r) search(r) path(ro + E) then Return end if
8: end for
9: S + sort X; E search(r) in the decreasing order of (search(r) n M(X))) #initialize the set of visited nodes
10: V +-0
11: [S] > then
12: while S do
13:
X + S.pop)
14:
q no)e(XX)
15:
scarch(g) + (M(X n scarch(r)) V
16:
path(q) - path(r)) {X}
17:
Recursive-Build(g)
18: V [u
19: end while
20: else
21: L L-
22: end if
23: End Procedure
X0
X2
X4
X5
X6
X2
X4
X5
X6
X7
X8
X6
X8
X2
X4
X5
X7
X8
X4
X6
X7
X7
X4
X8
X6
X8

---

Published in Transactions on Machine Learning Research (01/2026)
Table 1: Baselines and Data Models. The table indicates the applicability of each baseline to each data models. Applicability is determined based the baseline's empirical effectiveness on recovering the causal graph from observational data under the corresponding data model. Note that except for our versatile method, most other baselines are not applicable to all data models (i.e., producing very poor performance that is not meaningful for comparison).
Data Models
PC
GIES
FCI
Notears
MLP-Nootar
DAS SCORE
GLIDE(Our)
Linear, Gaussian (L-G)
Non-linear, non-Gaussian (nI-nG)
/
/
V
Categorical, Synthetic
V
Categorical, Real-world
V
+ PC
GIES
GLIDE
12
[E| = 3|V|
Spurious (%)
|E| = V
Runtime (x10 mins)
9
EE| = 3
|E| =
4|Vl
Spurious Rate (%)
440
|E| 4
Runtime (x102 mins)
EE| = 4
6
SHD (x10+)
SHD (x100²)
8
30
8
20
20
2
O
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
M
M
M
M
M
M4
Figure 6: Reported performance of the basclines across two cases:
(1)
E
3V and (2)
4V
C ADDITIONAL EXPERIMENTS
C.1 Experiment settings
Baselines & Hyper-parameters for baselines. In this section, present the experimental settings for each baseline used in the empirical report in the Section V in the manuscript and Appendix sections C.2 C.3
Firstly, the code for PC is from the Bnlearn package (Scutari 2009), provided in Python. We reported results of the PC-stable, which an upgraded version of PC. On the other hand, the implementation of GIES and FCI are provided from the edt package (Kalainathan el al. 2020 and the causal-learn package (Zheng et al. 2024). respectively. These algorithms require no tuned parameters as input and therefore can be used as recommended from the package from which they are provided. Notcars and MLP-Notcars are open-source Github, and are provided by Zheng et al. (2018). For experiments, reuse the parameters recommendation in the original papers: Notears (X1 = 0.1, Wis = 0.3) and MTPP-oteers (A1 X2 0.01, Wen 0.3). Tt is worth noting that we do not modify the original architecture of the network used in MLP-Notears. Lastly, DAS and SCORE are also open-source and are available from (Montagna et al. 2023 Both of these baselines share the same operational parameters as follows: TJG = 7H 0.001, K 10, pns 10, threshold = 0.05 and camcutoff 0.001.
Syper-parameters for our proposed method GLIDE. Our proposal has the following paper-pprramttrs:: (1) The number of prior distributions used for the invariance test; and (2) 7 factor which controls the ratio of the data generated from weak interventions. Both of these parameters are concerned with the data sampling procedure. Other than these two parameters, there are two more thresholds which control the level of tolerable variance we consider as invariant and the confidence interval for the Conditional Independence Tests (CITs). These thresholds are set 10-3 and 0.05, respectively, across all experiments reported this study. We recommend tuning the former threshold depends on the nature of the data whereas the latter should remain unchanged.
21

---

Published in Transactions on Machine Learning Research (01/2026)
Table 2: Performance on real-world datasets. N/A denotes no results within 48 hours.
Sachs (11)
Insurance (27)
Water (32)
Alarm (37)
Barley (48)
Pathfinder (186)
Munin
Metrics
Method
(1041)
GIES
14.0 = 0.0
36.0 + 0.0
49.0 + 0.0
44.0 + 0.0
65.0 + 0.0
1156.0 = 0.0 N/A
1235.0 + 0.0
SHD
PC
10.3 # 0.7
33.7 + 0.7
54.3 + 0.7
35.7 + 2.4
58.3 + 1.7
N/A
GLIDE
5.2 1 0.4
18.0 = 2.8
41.6 = 1.8
27.8 1 2.0
45.8 1 2.4
59.1 +: 1.9
883.2 1 21,8
Spurious rate (%)
GIES
34.8 +
31.9 0.0
35.8
44.3 0.0
35.8 0.0
87.0 + 0.0
42.4 0.0
=
0.0 = 0.0
17.1 = 2.1
27.7 3.6
17.8 1.9
N/A
N/A
GLIDE
0.0 1 0.0
3.6 16:1:27 2.7
23.0 1 0.9
13.1 1 0.9
8.2 1 1.1
1.9 0.5
1.8 1 0.2
GIES
1.5 + 0.4
2.7 = 0.6
3.1 = 0.6
3.3
0.4
3.4 + 0.8
19.3 0.4
61.5 + 15.2
Runtime (seconds)
PC
52.0 1 188
450.3 1 19.4
573.8 16.1
503.2 1 21.2
761.1 9.5
N/A
N/A
GLIDE
23.0
34.2 0.9
49.1 3.3
43.3 0.4
61.7 2.0
197.0 21.2
6200.3 + 88.7
C.2
Main result: Synthetic categorical data
Categorical data are generated in the manner continuous data in Section 5.1 The only difference is that instead of using Gaussian models, we first randomize conditional probabilities for each node and then use Gibbs sampling to simulate the data. The number of categories are randomly chosen from 2 to 5 for each variable. We compare our proposal with two competitive baselines, GIES and PC, while excluding the rest since they fail to produce meaningful results for comparison (i.c., very poor performance). To focus on settings where the PC baseline can produce results within the time limit of 48 hours, we restrict our evaluation to two classes of graphs where the number of edges is (i) 3x, and (ii) 4x the number of nodes, which ranges between 100 and 500.
The results are reported in Figure 6 which shows that GLIDE consistently outperforms the baselines in terms of SHD with substantial gaps of 10.49% and 27.95% in the (E 3V) case; and 11.56% and 13.3% in the (E 4V) case. In both cases, GLIDE the second best in terms of spurious rale, increasing the false causal detection (spurious) rate of the best baseline (PC) by a mere margin 4%. In exchange, GLIDE achieves a 30x faster processing time than PC. In contrast, GIES achieves the fastest rumning time but incurs 40% spurious rate (i.e., low reliability). Overall, GLIDE has the best trade-off between scalability (processing time) and performance (SHD, spurious rate).
C.3 Main result: Real-world graphs
Table 2 reported the performance of PC, GIES, and GLIDE on 7 real-world graphs provided by the Bulearn open library acutari (2009). Note that GIES is deterministic method and has zero deviation across different. seeds on the same test case. It is observed that GLIDE achieves the best SHD performance across 7/7 datasets and also achieves the best spurious rate in 5/7 datasets, especially on large datasets such as Barley, Pathfinder, and Munin. On the largest dataset (Munin), GLIDE achieves a spurions rate 1.8% which is remarkably better than GIES' (42.36%). Again, GLIDE has the best. balance between performance and scalability.
C.4
Ablation studies
C.4.1 Studies on the influence of and 7
The number of prior distributions m corresponds to the number of periionmennt/ssu-daataas that proposal generates. The higher this figure, the better the performance of invariance test. However, this comes with a trade-off on time consumption: the invariance test loops through all generated sub-datasets: therefore, it consumes increasing time linearly with the increase K. On the other hand, factor virtually cont rols the error in estimations produced by sub-datasets because y bounds the volume the downsampled datasets.
In this study, we conduct experiments on the 100-variable Erdos-Renyl categorical data graph and continuous data graph. In each data model, examine different values of (0.2, 0.4, 0.6, 0.8), along with increasing
22

---

Published in Transactions on Machine Learning Research (01/2026)
Table Influence of hyper-paramntters. Mean value and confidence interval 95% are reported.
Spurious rate
Runtime
7
SHD
(%)
(minutes)
0.2
122.2 + 9.788
5.91 1.678
1.95 = 0.037
0.4
127.4 + 10.31
5.41 0.951
2.02 + 0.127
10
0.6
113.6 + 2.525
4.64 1.01
2.39 = 0.113
0.8
121.0 + 7.717
4.19 + 1.552
2.45 = 0.114
0.2
116.8 + 9.139
5.53 0.841
4.54 + 0.221
0.4
109.8 = 8.796
4.30 = 0.792
4.18 = 0.086
20
0.6
123.0 + 2.772
5.96 + 0.957
3.88 = 0.098
0.8
112.8 + 5.958
3.27 + 1.448
3.93 + 0.095
0.2
117.0 + 4.679
5.19 + 0.505
5.67 = 0.068
0.4
112.8 + 7.874
4.37 0.613
5.66 = 0.012
30
0.6
102.0 + 6.263
3.85 1.603
5.58 + 0.218
0.8
106.0 + 9.461
4.15 0.981
5.75 = 0.194
0.2
114.0 + 5.714
4.79 0.753
7.44 = 0.032
0.4
112.8 + 5.088
3.94 0.529
7.50 + 0.028
40
0.6
108.6 + 8.072
3.97 + 1.536
7.30 = 0.051
0.8
108.8 + 4.444
4.03 + 1.279
7.36 = 0.048
0.2
100.8 - 2.432
5.03 + 0.518
8.97 = 0.024
0.4
106.6 + 5.834
4.64 + 0.958
9.35 = 0.061
50
0.6
99.01
= 6.261
3.68 1.589
9.21 - 0.021
0.8
107.6 + 7.606
4.61 0.955
9.22 + 0.039
values of m E {10, 20, 20,30, 40, 50). The input data volume 10000 and we examine 10 runs to report mean values and 95% confidence intervals as Table3 The results show a siguificant decrease (averagely, over 18%) in SHD when we increase the number of environments m from 10 50. However, such performance comes the expense of Runtime. As we discussed, the time complexity is lincar with m, thus, see that the runtime is about 5 times higher when is 5 times higher. As for y. The increase of y does not guarantee better performance, as can be seen in settings m = 10, m = 30, and m 50 where the best is neither the smallest. nor the biggest. Despite the effect of y is secondary to that of m, a good can boost the performance roughly 8% in term of SHD without incurring additional time complexity.
C.4.2 Studies on different topologies
In this section, we test the performance of baselines and our proposal on Scale-free graphs and Bipartite graphs (Zheng et al., 2018). The code used to generate the datasets for experiments is introduced by Zheng et al. 2018). We also evaluate our proposal and baselines in normal and extreme cases the former has the number nodes increasing from 100 to 500 and the number of edges equals the number of nodes, the latter has the number of nodes fixed at 500 while the number of edges increases from 600 to 1000. In each case, we also generate linear Gaussian data and non-linear non- non-Gaussian data scenarios. The results are depicted in Figures and g for Bipartite graphs and Scale-free graphs, respectively.
Bipartite graphs. Regarding Bipartite graphs, our proposal shows superior performance in almost all cases and data generation scenarios. As can be seen in Figure 7(a) upper, proposal achieves relatively similar to the strongest bascline (MLP-Notears) in terms of SHD and spurious rate, but costs over an order of magnitude smaller in runtime, about 15.68 times.
The performance gap is noticeable when comes to the non-linear non- non-Gaussian data scenario - Figure 7(a) lower. Our proposal now has clear advantage terms of SHD (roughly 10% performance gap) and spurious rate (17.21% performance gap) against MLP-Notears, while the runtime difference as significant as the previous scenario. It is worth noticing that, none of the baselines can achieve both fast runtime and relatively
23

---

Published in Transactions on Machine Learning Research (01/2026)
Notears
MLPNotears
DAS
SCORE
GIES
FCI
GLIDE
Runtime (minutes) Runtime (minutes)
Spurious Rate (%) Spurious Rate
400
SHD
200
100
20
100
300
500
100
400
500
200
300
500
800
100
300
600
SHD
60-
400
40
100
200
20
100
Number of Nodes
Number of Nodes
Number of Nodes
(a) Evaluation linear Gaussian (upper) and non-linear non-Caussian (lower) data
models. The number of edges cquals the number of nodes.
600
Spurlous Rate (%)
Runtime (minutes) Runtime (minutes)
SHD
2.5
100
200
14
700
900
1000
Spurious Rate (%)
15
800
200
10-
SHD
600
400
100
200
700
900
1000
700
800
900
600
700
Number of Edges
Number of Edges
Number of Edges
(b) Evaluation on linear Gaussian (upper) and non-linear non- non-Gaussian (lower) data models at extremes. The number of nodes is fixed at 600. Best baselines are selected to compare with our proposal.
Figure 7: Performance continuous data models in normal and extreme cases on Bipartite graphs. Apply for all metrics: Lower is better.
high performance. For example, GIES and DAS have the same runtime proposal, but their performance in terms of SHD and spurious rate are significantly worse (at an aaverage of 61 61.55% and 46.22%, respectively).
As for the extreme cases in Bipartite graphs, we isolate best baselines, Notears and MLP-Notears, to compare with our proposal. The overall result in Figure 7(b] is that our proposal outperforms both baselines in both scenarios and in all metries, albeit with an exception. Figure upper shows the result in linear Gaussian data scenario. We can see that our proposal has a noticcable gap the other basclines both SHD (12.81% and 27.18% lower than MLP-Notears and Notears, respectively) and runtime (16.62 and 14.18 times less than Mot-Notears and Notears, respectively). However, Notears marginally outperforms proposal in term of spurious rate, about 1%.
Figure 7(b) lower shows that our proposal returns relatively stable performance even in the non-linear non- Gaussian data scenario with an average of 0.69% spurious rate. Furthermore, our proposal consistently outperforms both Notcars and MIP-Notcars in all cases term of SHD. In contrast, Notcars and MLP- Notcars seem to be heavily affected by non- cineaiiyyy in term of spurious rate. MLP-Notears has an average of 10.06% spurious rate about. 5 times higher than linear Gaussian data scenario. Notears also suffers
24

---

Published in Transactions on Machine Learning Research (01/2026)
Notears
MLPNotears
DAS
SCORE
GIES
FCI
GLIDE
600
Spurious Rate (%) Spurious Rate (%)
100
Runtime (minutes) Runtime (minutes)
1000
80
400
750
SHD
60
500
40
250
20:
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
800
100
300
600
80
SHD
60
200
400
40
100
200
20
100
200
300
400
500
100
200
300
400
500
100
200
300
400
500
Number Nodes
Number of Nodes
Number Nodes
Evaluation on linear Gaussian (upper) and non-linear non- non-Gaussian (lower) data models on Scale-free graphs. The number of edges equals the number nodes.
Figure 8:
Degeneracy measure
20
Degeneracy measure
50
Degeneracy measure
200
40-
15
30
20
100
10
5
500
1000
1500
2000
500
1000
1500
2000
500
1000
1500
2000
Number of Edges
Number of Edges
Number of Edges
(a) Erdos Renyl graphs
(b) Bipartite graphs
(c) Scale-free graphs
Figure 9: Degeneracy measure different topologies. The number of nodes is fixed at 500.
average 3% higher spurious rate. Interestingly, both baselines seem to benefit from the increasing number of edges in the graph as the results show an downward trend in spurious rate. However, it is worth noticing that Notcars incurs a lincarly incrcasing in SHD as the mumber edges increases.
Scale-free graphs. Regarding Scale-free graphs, Figure X shows the results of algorithms normal cases. Overall, despite not being the best performer in all metrics, GLIDE maintains a good trade-off between all metrics. Regarding the L-G case (Figure uupper) shows different trends in performance of basclines compared to GLIDE in term of SHD. In details, the SHD of our proposal gradually increases until plateaus out in the 400 and 500-node graphs contrasting to other baselines (e.g., SCORE, and MLP-NOTEARS that show an accelerating rate in SHD towards 500-node graphs. Furthermore, these baselines are noticeably worse than GLIDE when it comes to spurious rate (GLIDE produces an average of 5.52% less spurious relationships than SCORE, 28.57% than MEP-NOTEARS while being comparable NOTEARS and DAS). As for the nL-nG case (Figure g lower), we have the same observation where our proposal consistently being comparable with other baselines term of SHD, having a low spurious rate and runtime, simultaneously.
As regarding the extreme cases the Scale-free graphs, when the number of edges exceeds 700, encounter the following problem: the average number of parents for each node becomes exponentially large due to the nature of the graph. Such that the performance of MLP-Notears is heavily impacted: MLP-Notears cannot. produce meaningful results within 48 hours. Our proposal GLIDE is also impaired by this setting. The reason that the data set is not sufficiently large to perform invariance test with adequate accuracy due the exponentially large number of parents. Therefore, we do not include the report on performance of GLIDE as well other baselines this particular setting.
25

---

Published in Transactions on Machine Learning Research (01/2026)
Table Degeneracy measures on real-world causal graphs.
Datasets
Sachs
Insurance
Water
Alarm
Barley
Pathfinder
Munin
Number of Nodes
11
27
32
37
48
186
1041
Degeneracy p
4
6
4
C.4.3 Studies on the Degeneracy measure
In this ablation study, we investigate the degeneracy measure p on different topologies and connect them the performance of GLIDE. As mentioned in Section 4.3 the time complexity that GLIDE requires to find the plausible parent. sets of all d variables is O(pd ((pd-3P//) As such, want to assert that, the degeneracy p of most graphs is indeed insignificant compared to the number of variables d of those graphs. To show this, we design the following ablation study: we generate DAG with fixed number of nodes (500) and an increasing number of edges (500 to 2000). For each setting of the number of edges, we randomly generate 10 graphs and record their degeneracy measure. Notice that most baselines presented in this research are incapable of rumning on graphs with 500 nodes and more than 1000 edges as we presented in Main text's Section 5 and Appendix C.4 However, for the sake of the argument of this section, we increases the number of edges 2000 to investigate the range value of the degeneracy measure p. Figure shows the mean and confidence interval 95% of the degeneracy measure different topologies (Erdos Renyi, Bipartite, and Scale-free graphs) as the nunber of edges increases.
Clearly, Bipartite graphs return the most stable degeneracy measure out of the three topologies. With very small variations, the degencracy measures on increasingly donse Bipartite graphs demonstrate an almost linear growth. We speculate that this stable behavior of the degeneracy measure of the Bipartite graph benefits the performance of GLIDE, we can see in Figure 7[b] where GLIDE consistently and noticeably outperforms the other two prominent baselines, especially in term of runtime. It is worth noting that at the (500-node, 1000-edge) setting, the degeneracy measure is roughly 15, which is much less than the number of nodes. This shows that, at cases where most baselines take hours to solve, GLIDE can still achieve almost. quadratic time complexity (Section 4.3) and thus only takes minutes to recover the causal graph.
Regarding the Erdos-Ro Mitos-RRnyiii graphs, we can see it random nature manifest in a noticeable yet not too significant variation of the degeneracy measure. Nonetheless, these values (even at maximum, e.g., p = 1000-edge or p 22 at 2000-edge graphs) are insignificant compared to the number of nodes. This partially explains the scalability of GLIDE the Erdos-Reny Renyi graphs as reported in Main text's Section 5 Contrasting to the previous two, the degeneracy measure Scale-free graphs shows an unstable behavior and wide range of fluctuation. The degeneracy measure on 500-edge graphs upto 700-edge graphs has an average of 52.4 + 14.27, which is highly unstable and generally much higher than that of the same set but on other topologies. As the number of edges gradually reaches 700, the degeneracy has a sudden leap upto an average of 112 27.43. This unusual (and perhaps, non-linear) behavior may stem from the implementation (Zheng ct al. 2018) or may need further studics address. Regardless, we can see that the range for degeneracy Scale-free graphs are significantly higher than that in Erdos-Renyi or Bipartite graphs, which potentially affect the runtime of Algorithm 3 However,
will see in Appendix C.4.4 the number of plausible sets in Scale-free graphs is still negligible compared to the number of nodes.
Finally, we also investigate the degeneracy of real-world causal graphs. These graphs are the same as the one presented in the Main text's Section C.3 As it turns out, real-world causal graphs are indeed sparse (see Table 4). In details, the degeneracy of the causal graph of these dataset are: (Sachs: 3, Insurance: 4, Water: 6, Alarm: 4, Barley: 5, Pathfinder: 5, Munin: 4). In which the Munin dataset has over 1000 variables. This shows that, in practice, indeed often encounter large causal graphs that sparse.
C.4.4 Studies on the number of plausible parent sets
The number of plausible parent sets plays an essential role the time complexity of our proposal. As per Algorithm 3 an enhanced Brrb-Shr algorithm (Bron & Kerbosch, 1973) each plausible set corresponds to maximal clique. Theoretical results from Bron & Kerbosch (173 give us the bound for the number
26

---

Published in Transactions on Machine Learning Research (01/2026)
Number Plausible Sets
Number of Plausible Sets
Number of Plausible Sets
2.0
1.8
15-
10
1.4
1.2
2-
1 6 7 9:
1
5
90
1
6 9
Number of Nodes (x10²)
Number of Nodes (x10²)
Number of Nodes (x10?)
(a) Erdos Erdos-Renyl graphs
(b) Bipartite graphs
(c) Scale-free graphs
Figure 10: Investigation on the number of plausible parent sets in different
topologics. The number of edges equals the number of nodes.
Number Plausible Sets
Number of Plausible Sets
Number of Plausible Sets
300
20-
100
200
100
5
1
6
90
2
5 8 90
6 8 10
Number Nodes (x10²)
Number of Nodes (x10²)
Number Nodes (x10²)
(a) Erdos- Renyl graphs
(b) Bipartite graphs
(c) Scale-free graphs
Figure Investigation on the number of plausible parent sets in different
topologies. The number of edges doubles the number of nodes.
Number of Plausible Sets
Number of Plausible Sets
6
100
20
Number Plausible
150
N0ggg
60
10-
000336
1000000
40
000000
00
50
500 600 700
9001000
500 700 800 90m 000
500 600 700 800 90m 000
Number of Edges
Number of Edges
Number of Edges
(a) Erdos-Rent Renyi graphs
(b) Bipartite graphs
(c)
Scale-free graphs
Figure 12: Investigation the number of plausible parent sets in different topologies. The number of nodes is fixed at 500.
of maximal clique, which O(dd-)) 33/3) 3p/33 where p is the degeneracy of the graph. As we have shown Appendix Scction C.4.3, p is indced insignificant compared to the common Erdos-Renyi and Bipartito class of graphs. In this section, empirically investigate the number of plausible parent sets in different. graph scenarios. For each scenarios, randomly generate 10 graphs using the source code provided by (Zheng) ct al. 2018) and perform Algorithm then record the number of plausible parent sets of cach node. We use the violin plot to show the distribution of the number of plausible sets at each setting. The maximum, minimum, mean, and median are reported.
Figure 10 shows the number of plausible sets on graphs whose number of nodes equals that of edges, which gradually increases from 100 1000. Notice that this setting mimics that of experiments that we reported in normal cases in Main text's Section 5.1 but the number of nodes grows from 100 500. As we see
27

---

Published in Transactions on Machine Learning Research (01/2026)
Table A case study on sortability, Comparison with Varsort and R2sort
SHD
Spurious rate (%)
Time (min)
Varsort- R2sort-
EE|
Varsort R2sort
Ours
Varsort R2sort Ours
Varsort R2sort Ours
ability
ability
L-G
500
0.98
0.76
92
814
249.12
12.74
55.28
1.98
11.01
8.56
27.25
L-G
1000
0.98
0.92
739
4600
539.6
40.83
81.43
5.78
10.76
10.21
15.19
nL-nG
500
0.88
0.32
587
1316
340.66
50.99
68.02
4.18
27.57
26.21
7.25
uL-nG
1000
0.87
0.23
1642
4157
710.66
59.35
78.41
7.82
22.28
21.99
5.43
in Figure 10 a); the number of plausible sets on Erdos-Renyi graphs peaks at 10, with mean being roughly 2. This supports the excellent runtime of GLIDE as we reported. When it comes to Bipartite graphs, the number of plansible sets is approximately double that in the same settings on the Erdos-Reng graphs. However, in all cases, this figure is much lower than the number nodes. Noticably, in the case of Scale- free graph, we see an almost similar distribution of the number of plausible sets across graphs with different nodes. In details, most nodes have only plausible parent set, and other have 2. This observation combined with nature of the Scalo-free graphs, suggest that each plausible set contains large number of nodes. Consequently, given that the observational data volume is limited at 10000, the accuracy of the invariance lest of GLIDE reduces significantly, which explains the performance that we reported in Appendix Section C.4
We also conduct experiments on graphs that are twice as dense (the number of edges is doubled that of nodes) as in the previous setting, and the results are reported in Figure 11 Interestingly, while the number of plausible sets on Erdos-Reny's roughly doubles (both maximum and medium) compared the previous setting, the that on Bipartite graphs increases significantly. In details, in most setting, the maximum number of plausible ranges from under 100 (on 100-node graphs) to over 300 (on 400-node graphs), which is innegligible compared to the number of nodes. Nontheless, these figures are still strictly less than the number of nodes. Tt. is worth noting that most nodes in the graphs have about 10 to 15 plausible parent sets. As for the Scale-free graphs (Figure [c)))) despite the median number of plausible sets is 10, the maximum increases lincarly with the number of nodes and then plateaus out at about 150 when the number of nodes reaches 900.
Lastly, we study the number of plausible sets extreme scenarios (as in Main text's Section 5.1 we fix the number of nodes 500, and increase the number edges from 500 to 1000. The results are depicted by box plot in Figure 12] The number of plausible sets the Erdos-Reny" graphs is strictly less than 25 even in the hardest sotting (1000-edge graphs) while the mean value is roughly 3.8 and 75%-percentile is at 5. In contrast, the number of plausible sets on the Bipartite graphs grows non-linearly, reaching a maximum (outlier) 180 on 1000-edge graphs. However, since most nodes do not have more than 5 plausible parent sets, we see that the outliers have little effect on the overall performance of GLIDE, as can be seen in Figure 7( On the other hand, the number of plausible sets on Scalo-free graphs is very unstable, and (intcrestingly) has the same trend as the degeneracy (reported in Figure 9 c))- In details, when the number of edges is below 700, we have the number of plausible sets ranges from 1 to 3 with the mean 1.4. But when the number of edges exceeds 700, despite the mean of 5 and 75%-percentile of 10, the maximum number of plausible sets may reach to over 100.
C.4.5 Independence of varsortability
In this section, we include baselines Varsort (Reisach et al. 2021) and R2sort (Reisach et al., 2023) into experiments, to show that the performance of GLIDE is independent of sortability measures proposed in these works. Sortability is the score defined in the interval [0, 1] where 0 implies extreme difficulty to find the topological order over the variables, and implics extreme casiness. Following the same practice in the main texl, we run Varsort and R2sort on Erdos-Renyi graphs, with 500 nodes and (500, 1000} edges. The results are shown in Table 5 Generally, despite high var-sortability, Varsort is outperformed by GLIDE in 3/4 cases. Specifically, GLIDE is evidently superior when it comes to extremely large graphs, as can be seen in a stark reduction in term of both SHD (up to 2.2x) and Spurious rate (up to 8x). R2-sortability, on the other hand, plummets we move from Linear Gaussian models to non-Linear Gau-Saas models, which results
28

---

Published in Transactions on Machine Learning Research (01/2026)
Table Varsort, R2sort, and GLIDE on real-world graphs.
Spurious
Runtime
Datasets
Baselines
Sortability
SHD
rate (%)
(min)
Varsort
0.64
105
62.01
2.54
Insurance
R2sort
0.18
154
69.93
2.71
GLIDE
18
3.6
34.2
Varsort
0.69
93
51.18
4.03
Water
R2sort
0.41
145
62.96
3.39
GLIDE
41.6
23
49.1
Varsort
0.55
108
65.87
5.29
Alarm
R2sort
0.59
118
68.57
5.38
GLIDE
27.8
13.1
43.3
Varsort
0.76
127
54.33
7.78
Barley
R2sort
0.36
181
61.05
GLIDE
45.8
8.2
61.7
Varsort
0.79
1496
89.29
43.91
Pathfinder
R2sort
0.05
2433
92.94
44.94
GLIDE
59.1
1.9
197
in its poor performance. Overall, GLIDE is more preferable compared both of these baselines in synthetic scenarios.
We also conduct the experiments on real-world scenarios (graphs provided in the bulearn package), and the results are shown in Table 6 It is clear that real-world graphs are much harder to solve as evidently shown by the low sortability of both Varsort and R2sort. Regardless, GLIDE outputs a relatively good result with respect to such difficulties. Especially, rate of 1.9% whereas that of Varsort and R2sort are 89.29% and 92.94%, respectively. Furthermore, we also record significantly lower in SHD, almost 42x lower than R2sort and 25x lower than Varsort.
the Pathfinder dataset, GLIDE achieves a remarked low spurious
In conclusion, the performance of GLIDE is relatively independent of sortability (Reisach et al. 2021 2023 empirical evidences have shown. Furthurmore, these experiments also show that our experimental design is extensive and not biased.
D SUPPORT MATERIALS
6(B b*)
8(B = b*)
(
p
(.6r) P(B)
(.6%)
P(B)
P(B)
(
p()
p(u
0g
p
8(B b))
5(B = b²)
8(B b¹)
8(B b²)
P(
(a) Illustration the space (r(r)) original prior P(B), and point mass distributions 8(B).
(b) Step 1: Computing the boundary points (B)(B via Eq. 17, which forms
(c) Step 2: Performing weighted aver- age of the houndary points to generate midpoints P.(B).
a convex hull.
Figure 13: The procedure generating prior distributions.
29

---

Published in Transactions on Machine Learning Research (01/2026)
0.6
0.6
0.4
0.3
0.2
0.1
0.1
0.05
0.3
0.1
0.15
2,
96
0.2
0.55
0.30
0.9
0.8
0.9
0.8
(a) Before applying K-means.
(b) After applying K-means.
Figure 14: A mock-up experiment: Representative distributions selection via K-means.
D.1 On the generation of prior distributions
The fact that we can sample an infinite number of prior distributions P((B) in C+(70) is not helpful from the computational point of view. Because larger number of prior distributions incurs corresponding increase in the processing time of the invariance lest. Furthermore, when the number of sampled prior distributions grows too large, it is inevitable that there exist similar prior distributions, which in fact does not improve the quality of the invariance test. On the other hand, an insufficient number of prior distributions reduces the reliability of the invariance test. Therefore, we wish to select only an adoquate number of representative prior dist risbttitionns within (70)).
Inspired by the aforementioned logic, we perform the following:
1. Sampling PB( where i goes up to 10 (Figure 14(a) to overwhelmingly fill in the convex hull C,(7)) We use the Dirichlet distribution to sample weighted vectors (
2. Using the K-means algorithm with parameter K = m to cluster them into m partitions, each corrosponds to a representative controid. These centroids are the output prior distributions of the procedure (Figure 14(b)).
To this end, have sampled m prior distributions that as variants of the original P(B) and can be used to fuel the invariance test by applying Theorem 6
D.2 Further Discussion
Parallelism Prospect. Our proposed framework GLIDE is readily applicable to scenarios where the data is distributed across multiple, private local devices such as the popular federated learning setting. In such scenarios, GLIDE will benefit directly from the fact that the distributed datasets, which are presumed be governed by the same causal model, are already admitting identical underlying conditional distributions. As such, GLIDE can use those local datasets synthetic augmented datasets which are essential to the effect-cause distributional invariance test. This save time and help avoid unnecessary sampling errors.
Limitations. The proposed framework GLIDE has the following limitations:
1. Our performance partially relies on how well Algorithm 2 can find the source nodes for the basis. As previously mentioned, not guaranteed that Algorithm 2 will find all the basis variables. However, our extensive experiments have shown that when the number of data augmentations increases, the chance that non-source nodes are included in the basis set is lessen.
30

---

Published in Transactions on Machine Learning Research (01/2026)
2. GLIDE's performance might be hampered on scenarios where the observational data is insufficient to discern the true complexity of the underlying causal graph can be seen with scale-free graphs.
E REPRODUCIBILITY
Software. Our implementation is in Python. The requirements include the installation of the causal-learn package (Zheng et al. 2024) for the use of CITs, scikit-learn Python package, along with pandas and numpy, which are standard libraries commonly used in Python.
Support Module for Continuous Data. For the use of our proposal for continuous data models, we use a Discretizer module provided by the Python standard scikit-learn library. The number of discretizing bins is fixed at 4 and the width bins is equal to capture the marginal distribution of each variable the observational data. As for the categorical data models, this module is deactivated.
Hardware. All experiments are run machine with 64-core Intel(R) Xeon(R) Gold 6242 CPU @ 2.80GHz. The running of GPU-based baselines is conducted on an NVIDIA GeForce RTX 4090.
31


Reviews:
[{"review_id": "XxYLTPg6yh", "reviewer": "Reviewer_Gmtf", "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": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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": "This paper introduces GLIDE, a new framework for causal discovery from observational data1. Its core idea is that the conditional distribution $P(\\text{effect} | \\text{causes})$ is invariant to changes in the prior $P(\\text{causes})$. GLIDE proposes a novel \"invariance test\" that uses strategic downsampling to simulate different priors and check if the variance of $P(X|Z)$ is near-zero. By combining this test with Markov blanket identification and graph sparsity assumptions, it efficiently finds parent sets, claiming an $O(d^2)$ complexity. Experiments show it matches or beats existing methods on accuracy while being up to 25x faster.\n\n**Key Strengths:**\n\n- Novel application of distributional invariance principles to causal discovery\n- Interesting downsampling approach\n- Comprehensive empirical evaluation across multiple graph types and datasets\n- Clear computational complexity analysis with claimed significant improvements\n\n**Key Weaknesses:**\n\n- Gaps in theoretical proofs, particularly Theorem 1 which forms the foundation\n- Limited principled guidance for hyperparameter selection (m, γ, ε parameters)\n- Unclear positioning relative to other invariance-based causal discovery methods", "claims_and_evidence_yn": "No", "claims_explanation": "The theoretical claims are not adequately supported due to gaps in the mathematical proofs: \n\n**Theorem1:**\n- Point 2 is insufficiently justified: The proof claims that if Z contains a child Y of X, then P(X|Z) has variance > 0 across different source priors. The reasoning provided that P(Y|Z\\Y) \"lacks X in its conditioning set\"doesn't rigorously establish why this creates dependence on source variable distributions. The connection between d-separation properties and distributional variance across resampled datasets needs much stronger mathematical foundation.\n- Point 4's spouse detection claim is unsubstantiated: The proof asserts that the PrunedX operator can reliably detect and remove spouse variables. The reasoning for why \"X ⊥⊥ Y | Z'\" would hold for spouses in the specific cases outlined is not convincingly established.\n- Point 5 makes an unjustified leap: Why is Pa[X] necessarily the *unique* maximum-sized set among PrunedX(Z₀)? The proof doesn't rule out multiple sets of equal maximum size\n  - The proof fails to rule out the following critical possibilities:\n    - That another pruned, invariant set exists with the same size as Pa[X]\n    - That a set Z containing Pa[X] and a spouse Y happens to pass the invariance test (𝕍 ≈ 0), and Y is precisely the \"pruning failure\" case shown in Figure 5. In this scenario, this set Z would be larger than Pa[X], causing the algorithm to fail.\n\n**Theorem 2:**\n- The induction only shows |B(G)| ≤ |S(G)|, but proving the maximum basis size equals the number of sources requires showing equality, not just inequality\n\n**Theorem 3**\n- Lemma 1's claim about nodes with minimum dependents being sources lacks rigorous justification\n- The argument about Type-1 colliders and why there \"must exist\" a node Z with smaller |Φ(Z)| is not mathematically sound\n\n**Practical Implementation:**\n- The method shows substantial sensitivity to hyperparameters but provides limited guidance for practical parameter setting", "audience_yn": "Yes", "audience_explanation": "Causal discovery is a fundamental problem in machine learning with broad applications across scientific domains. The paper tackles scalability issues that limit practical deployment of causal discovery methods, which is highly relevant to the TMLR audience. The distributional invariance perspective, while building on existing work, offers a novel algorithmic approach that could inspire further research. The claimed computational improvements, if validated, would be significant for practitioners working with large datasets. Even with the theoretical and experimental limitations, the core ideas are sufficiently interesting to warrant attention from researchers working on causal inference, graphical models, and related areas.", "requested_changes": "**Critical Changes:**\n\n1. **Fix theoretical foundations**: Rigorously prove Theorem 1, particularly Point 2's claim about variance behavior and Point 4's spouse detection mechanism. \n\n2. **Provide principled hyperparameter selection**: Develop theoretical or empirical guidelines for setting m, γ, and ε parameters. The current ad-hoc approach limits practical applicability.\n\n3. **Strengthen basis-finding algorithm**: Provide rigorous proof for Theorem 3, particularly Lemma 1's claim about minimum dependency nodes being sources.\n\n**Changes that would strengthen the work:**\n\n1. **Improve positioning relative to related work**: Better compare against other invariance-based causal discovery methods and explain the specific advantages of the proposed approach.\n\n2. **Add failure case analysis**: Analyze scenarios where the method fails and provide guidance on detecting such cases in practice.", "broader_impact_concerns": "This work presents a methodological improvement to causal discovery algorithms, focusing on computational efficiency and scalability rather than introducing new application domains or fundamentally changing how causal discovery is used. The potential risks (misinterpretation of causal relationships, over-confidence in results) are inherent to all causal discovery methods and are not specifically amplified by this approach.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "badcPVUpQA", "reviewer": "Reviewer_pMr2", "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": "", "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": "This work introduces a framework for causal graph discovery that exploits the invariance of the conditional distribution of effects given their causes to changes in the prior distribution of the causes, as implied by the local Markov property. The method employs a sampling strategy to generate data with altered cause distributions relative to the original sample, followed by an invariance test that evaluates whether the conditional distribution of effects given causes remains unchanged across these generated samples. Experimental results demonstrate that the proposed framework achieves performance that is superior or comparable to existing benchmarks, while maintaining quadratic time complexity.", "claims_and_evidence_yn": "Yes", "claims_explanation": "I found this paper engaging and nice to read. The main idea of testing invariance across generated samples seems a clever and scalable approach to causal discovery. The theoretical developments are carefully made and presented and appears to have broader applications beyond causal discovery, such as the proposed sampling procedure. While several theoretical compromises are made for practical implementation, for instance, substituting basis variables for source variables and assuming sparsity in the graph, these choices are well justified, clearly explained, and applicable to a broad class of graphs.", "audience_yn": "Yes", "audience_explanation": "The experimental results show favorable performance of the proposed framework in uncovering causal graphs from observational data. Although it does not outperform MLP-NOTEARS in some settings, this is a reasonable trade-off given the significant gains in computational efficiency.", "requested_changes": "My suggestions primarily concern the presentation. Some definitions and notations are either inconsistent or not clearly specified. In addition, the exposition does not always follow the order of implementation steps, which makes the paper somewhat less straightforward to read:\n\n1. It is not entirely clear whether Assumption 1 is truly an assumption or rather a property implied by the local Markov property. It seems more like the latter. By contrast, the statement “Accordingly, our proposed method relies on the contrapositive: if invariance does not occur, causality does not exist” functions as an assumption, but it is not explicitly emphasized or formalized outside of the introduction.\n2. The notation is inconsistent in several places, which makes the context difficult to follow. For instance, the paper appears to use $\\boldsymbol{Z}_i \\in \\mathcal{Z}_0$ and $\\boldsymbol{Z} \\in \\mathcal{Z}_0$ interchangeably, which can be confusing. In Equation (1), $\\prune_X(\\boldsymbol{Z}_i)$ is defined as the operator that removes all $Y \\in \\boldsymbol{Z}$ such that $\\exists \\boldsymbol{Z}' \\in \\boldsymbol{Z}_i : X \\perp Y \\mid \\boldsymbol{Z}_i$. It is unclear what the $\\boldsymbol{Z}$ in $Y \\in \\boldsymbol{Z}$ refers to, is this a typo, or does it denote a different set from $\\boldsymbol{Z}_i$? Also, does $\\boldsymbol{Z}'$ include the empty set? In addition, the distinction between Equation (1) and Equation (3) is not clear, as they appear to be repetitive. Moreover, Equation (1) uses $\\boldsymbol{Z} \\in \\mathcal{Z}_0$, whereas Equation (3) uses $\\boldsymbol{Z}_i \\in \\mathcal{Z}_0$, which adds to the confusion.\n3. The meaning of $\\mathbb{V}_{P_{+}(\\boldsymbol{X})} \\sim \\mathcal{P}\\left[P_{+}\\left(X \\mid \\boldsymbol{Z}_i\\right)\\right]=0$ in Equation~(1) in not stated or explained.\n4. What is the intuition behind initializing $\\boldsymbol{B}$ with $X \\in \\boldsymbol{V}$ having the smallest $|\\Phi(X)|$? It seems that the selected $X$ is the variable least connected with others in the graph, but it is not clear why this makes $X$ a favorable candidate for inclusion in $\\boldsymbol{B}$.\n5. Could you expand on this ``To guarantee this, $D_i$ must be a downsampled version of $D$ via sampling with no replacement to avoid introducing duplicates and hence, false causal biases into $D_i$.''? In particular, why is sampling without replacement necessary in this context?\n6. What does $|D|/|D_i|$ mean? Does it denote the ratio of unique observations, or something else?\n7. Below Equation~(3), what does this line mean? ``Furthermore, it is worth noting that $\\pa[X]\\in\\Z_0$ is always guaranteed while $\\boldsymbol{Z}_i\\notin \\pa[X]$ only exists in $\\Z_0$ in some specific scenarios.''\n8. Clarify what is meant by ``(empirical) sparsity'' in the statement ``Using this invariance test and exploiting an (empirical) sparsity of most causal graphs''? An expanded description or discussion would be helpful, particularly in relation to the application scenarios of the proposed framework.\n\nTherefore, in general I request changes:\n1. Providing an algorithm summarizing the proposed framework for causal discovery. Make clear the applicable scenario, and mark where theoretical compromises are made.\n2. Use consistent notation throughout and to clearly define each notation when first introduced, so that the presentation is fully self-contained.\n3. Expand on key definitions.", "broader_impact_concerns": "None", "additional_comments": "None", "extra_scores": {}, "extra_text": {}}, {"review_id": "K1DeWXQWBG", "reviewer": "Reviewer_gpEM", "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": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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": "**Summary:**\n\nThis work studies the problem of recovering a causal DAG from observational data. Relying on observational data alone limits one to only being able to recover the DAG up to Markov equivalence class, but it does allow one to probe causal relationships even in the absence of interventional experiments.\n\nExisting approaches to this problem typically fall into one of three categories: constraint-based, score-based, and model-based approaches. All have drawbacks: constraint-based approaches can have exponential computational cost, score-based approaches rely on heuristic objectives, and model-based approaches can be unstable when the model assumptions are violated. The authors propose a new model-based method, which they claim mitigates all three above drawbacks.\n\nThe authors’ core insight is that if one assumes effect-cause conditional distributions, e.g., P(effect | cause), are invariant to changes in the prior causal distribution P(cause), then one can efficiently leverage this relationship. For example, under this assumption, if one can construct sub-sampled datasets that reflect different causal priors P(cause) then one can rule out cause-effect pairs by concluding that no causal relationship (cause -> X) exists for any variable X such that P(X | cause) is not invariant across the different data subsets.\n\nSurprisingly, the authors argue that under the causal invariance assumption, not only can they rule out cause-effect pairs in the above manner, but actually that with high probability, one can provably recover the full set of cause-effect pairs from such a procedure.\n\nThe procedure relies on an existing method for recovering the Markov blanket of every variable from observational data alone. In order for the recovery guarantees for this pre-existing method to apply, the authors must assume that no latent confounders impact the causal relationships reflected in the observed variables, which is a reasonable assumption within causal discovery literature.\n\nThe resultant algorithm recovers the underlying causal graph with computational complexity quadratic in the number of observed variables, which is very efficient by the standards of causal discovery.\n\nThe core computational speedup (and the bulk of the computational cost of the proposed procedure) comes from the execution of a quadratic-time algorithm, which finds the Markov blankets (subsets, one per variable, which contain all parent nodes, child nodes, and child nodes’ other parents).\n\nThe work includes empirical results showing that on the synthetic and real-world datasets considered, the proposed novel method provides both speedup and a reduction in false-positive causal relationships.  \n\n**Weaknesses:**\n\nMy primary concern with the manuscript is that I am confused about some of the core theoretical claims. Here are some of my most fundamental confusions.\n- Identifying the parent set through invariance. I am confused about the claim made in Theorem 1, namely that the parents of a node can be identified as the largest set of variables for which “when conditioned on, the probability of observing X becomes unchanged.” Isn’t it possible that X would be invariant when conditioned on other ancestors, not necessarily only parents? For example, consider some path graph A-> B -> X such that B is a copy of A. Then conditioning on A should have the same impact as conditioning on B. Are the Markov blankets somehow being leveraged to prevent considering non-parent ancestors? Are they implicitly constraining the candidate set $\\mathcal{Z}$?\n- Identifying parent sets despite only targeting recovery up to Markov equivalence. As acknowledged by the authors, “Finding the true causal graph is not possible without running intervention[s] since different graphs can entail the same set of statistical constraints induced by observational data” (footnote, page 1). Thus when using only observational data, one’s goal must be restricted to identifying causal DAGs up to Markov equivalence. How then are the parent sets being recovered in Theorem 1? Is there some notion of parent set up-to-equivalence that the authors are employing?\n- Reliance on correctly identifying source nodes. Does the proof of Theorem 4 rely on the fact that B is assumed to be a source variable? The proof references the fact that B is a source variable, but as acknowledged by the authors, “finding the true sources is not possible without intervention.”", "claims_and_evidence_yn": "No", "claims_explanation": "As discussed in the Weaknesses section of “Summary of Contributions,” I am confused about the correctness of some of the core theoretical claims.\n\nIf the authors can correct my confusions about the theoretical claims, I am open to changing this answer.\n\nThe empirical results provide evidence for the practicality of the method.", "audience_yn": "Yes", "audience_explanation": "If the paper’s theoretical claims are correct, then I think that readers would be excited by a novel approach to causal discovery.", "requested_changes": "Some concrete changes that would improve readability would include:\n- Clarifying that causal invariance is necessarily an assumption. When first introduced in the text, the authors’ presentation may mislead readers into thinking that causal invariance is a fact rather than an assumption. See for example, in the abstract “leveraging the fact that the distribution of an effect, conditioned on its causes, remains invariant to changes in the prior distribution…”, and the text on page 2, which states “Our approach leverages the invariance of the effect-cause conditional distribution … to changes in the prior cause distribution.”\n- Minor: re-organize proofs so that each result only relies on previously-proven results. E.g. the conclusion of the proof of Theorem 4 seems to rely on Theorem 5. \n\nTypos:\n- In Theorem 1, the line beginning “For each variable…” does not seem to parse grammatically to me. Should the authors remove the word “then” following the comma, or is there some quantifier or modifier missing from the statement?", "broader_impact_concerns": "I am not concerned about the ethical implications of the proposed method.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)