Differentiable Constraint-Based Causal Discovery
Authors: Jincheng Zhou
Mengbo Wang
Anqi He
Yumeng Zhou
Hessam Olya
Murat Kocaoglu
Bruno Ribeiro
Venue: NeurIPS 2025 poster ()
Keywords: Causal Discovery
Constraint-based
d-Separation
Continuous optimization
Structure learning
Primary area: probabilistic_methods
TL;DR: We propose a percolation-based differentiable d-separation framework that bridges constraint-based causal discovery methods with differentiable DAG structure learning.

Abstract: Causal discovery from observational data is a fundamental task in artificial intelligence, with far-reaching implications for decision-making, predictions, and interventions. Despite significant advances, existing methods can be broadly categorized as constraint-based or score-based approaches. Constraint-based methods offer rigorous causal discovery but are often hindered by small sample sizes, while score-based methods provide flexible optimization but typically forgo explicit conditional independence testing. This work explores a third avenue: developing differentiable $d$-separation scores, obtained through a percolation theory using soft logic. This enables the implementation of a new type of causal discovery method: gradient-based optimization of conditional independence constraints. Empirical evaluations demonstrate the robust performance of our approach in low-sample regimes, surpassing traditional constraint-based and score-based baselines on a real-world dataset. Code implementing the proposed method is publicly available at [https://github.com/PurdueMINDS/DAGPA](https://github.com/PurdueMINDS/DAGPA).

Differentiable Constraint-Based Causal Discovery
Jincheng Zhou'
Mengbo Wang'
Anqi He? Yumeng Zhou? Hessam Olya?
Murat Kocaoglu Bruno Ribeiro'
Purdue University
Ford Motor Company 3Johns Hopkins University
{zhou791, wa1g48an ribeiroi @purdue edu y moy ma @fod.. com
mkocaoglu@jhu. edu
Abstract
Causal discovery from observational data is fundamental task in artificial in- telligence, with far-reaching implications for decision-making predictions, and interventions. Despite significant advances, existing methods be broadly cate- gorized as constraint-basee or score-based approaches. Constraint- based methods offer rigorous causal discovery but are often hindered by small sample sizes, while score-based methods provide flexible optimization but typically forgo explicit conditional independence testing. This work explores a third avenue: developing differentiable d-separation scores, obtained through a percolation theory using soft logic. This enables the implementation of new type of causal discovery method: gradient-basse optimization of conditional independence constraints. Empirical evaluations demonstrate the robust performance of our approach low-sample regimes, surpassing traditional constraint based and score-based baselines on a real-world dataset. Code and data of the proposed method are publicly available bttps://giive titb.bb... wwtm/oommu/ueeee httmm/urduuMMNND/DDDDM
Introduction
Causal discovery-the task of inferring causal relationships from observational data-is funda- mental problem in machine learning and statistics with far-reaching implications across multiple disciplines, including biology, economics, social sciences, and medicine [19, 25, 28, 33, 45). Directed Acyelic Graphs (DAGs) provide powerful framework for representing causal relationships, enabling the estimation the effects interventions or actions. However, many complex systems, the underlying causal graphs remain unknown, and conducting randomized controlled trials (RCTs) to establish causal relationships can prohibitively expensive or impractical. This has significant impli- cations in various domains, including planning, explainability, and faimeess [26, 47, 49]. highlighting the need for reliable methods to infer causal structures from observational data alone.
The task of causal discovery is inherently challenging, specifically under small sample sizes. Consid- ering Markov, acyclicity, and faithfulness assumptions [25, 34), some recent works have focused on two popular avenues [13, 39, 44]: differentiable acore-based methods, which leverage differentiability to optimize objective functions, and constraint- based approaches, which rely on conditional indepen- dence tests to infer causal structures. Common methods for constraint-basel causal discovery, such as the PC algorithm [34], are vulnerable to errors in datasets with small sample sizes, due to uncertainty in conditional independence tests (CI tests). More recent constraint-bace methods, such LOCI and k-PC [18, 42] among others, have ameliorated these challenges with smaller conditioning sets.
In recent years, significant progress has been made addressing the challenges of small sample sizes in causal discovery via differentiable score-based methods, casting the combinatorial graph search problem as a continuous optimization problem on weighted adjacency matrices representing directed
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

graphs [4, 36, 51, 52]. These differentiable methods rely either on linear models or increasingly complex neural networks to model the underlying functional relationship between variables.
Contributions. promising avenue for advancing causal discovery, complementing existing ap- proaches, is the development of hybrid methods [13, 37, 39). Our work tackles key challenge in this domain by introducing differentiable functions of d-separation and d-connection, capable of scoring conditional independencies within probabilistic causal graph, bridging constraint-based and gradient-based methods. To achieve this, we propose a novel framework grounded in a percolation measure applied to a continuous relaxation of the causal graph structure. This new measure effec- tively captures the inherent dependencies among paths in causal graph, addressing a significant limitation of conventional diffusion- based methods, such as matrix powers, which would be unable to adequately model d-separation.
The distinction between diffusion and percolation lies the heart of deriving differentiable metric for d-separation, which essentially boils down to measuring reachability over random graphs [16]. To illustrate this concept, consider a typical setup in gradient-base graph optimization, where matrix A E [0. _| represents edge probabilities [41, 51], with each entry AxY encoding the likelihood of an edge from node X node Y. A straightforward approach assess connectivity between X and Y might involve computing A which represents the sum of all path probabilities of length from X Y. However, this method implicitly assumes independence between paths, akin diffusion process. For instance, in the paths X - Z - Y and X + Z + W + Y, the simultaneous absence of the sampled X + Z edge would invalidate both paths, introducing a probabilistic dependence that is not captured by matrix powers or diffusion in general.
In contrast, d-separation measure over A is type of percolation measure, which bounds the probability of sampling valid paths between X and Y, accounting for edge dependencies. Although characterizing percolation involves combinatorial complexity, we derive differentiable bound for d- separation using soft logic, enabling gradient-ban based optimization of causal structures. We instantiate this differentiable d-separation framework in an algorithm (DAGPA, for DAG Percolation Apartness), where we combine CI- Cl-based scores, state-of-the-as techniques in multi-task learning, and Bayesian sampling. While DAGPA instantiates this novel framework, we believe differentiable d-separation is of independent interest and can be integrated into existing gradient-bast based methods.
Empirically, DAGPA shows that differentiable d-separation yields strong results in the small-sample regime, demonstrating good robustness and accuracy compared to popular constraint-base and differentiable score-based baselines. Moreover, on the real-world Sachs dataset [28), DAGPA shows that differentiable d-separation offers accurate modeling of the independence patterns in the data, outperforming baselines in our metrics.
2 Notation and background
Notation. We denote the integer set {1...., dd as [al] and dataset of n samples and d variables as D. We represent a weighted directed graph with d nodes and edge weights in the range [0, [0 by a real square matrix W E [0, _] Given a node 2 E [d]. we denote W-3 the submatrix obtained by removing the z-th row and column of W, which is equivalent to removing the node and their connecting edges from the graph. When the square matrix is binary, e.g. A E {0, 1] dxdd we interpret it unweighted directed graph. Here, we write x A y to say "y is connected from x via directed edge in A, and x ~ A y/ to say "y is reachable from r via a directed path (a path where every edge has direction from of these statements. Finally, we use 11 D for (conditional) independence statements in the dataset D I and if A is acyclic, II. A for d-separation statements in the DAG A.
y) in A. "" Conversely, x A y and x A y stand for negations
Constraint-based causal discovery. Methods in this category aim to identify collection of DAGs from a dataset D such that each DAG A in it is an I-Map of D, meaning all d-separations in A imply conditional independencies in D with sufficiently large dataset. Namely, Yz, y E [d], C [d].
(x 11 A Z)
(x 1L by 2).
The collection of all DAGS satisfying this criterion constitutes the Markov Equivalence Class (MEC). Our method, as most methods for learning DAGs, relies on the Causal Faithfulness Assumption, which
Assuming infinite samples in D. IL y means P( z) P(x |x) in the true data distribution.

---

posits the unique existence of DAG A' such that the observed data D exhibits only the conditional independencies represented by the d-separations in A. This assumption equivalent the converse direction of Equation (1). However, challenges emerge when dealing with large conditioning sets 2, as statistical tests for conditional independence become increasingly unreliable. To mitigate this issue, recent approaches like LOCI [42] and k-PC [18] restrict their search to DAGs with low-order conditioning sets (typically limited to a cardinality of or less). As a result, their solution space is not the entire MEC, but rather the keequivalence class (or k-essential graphs) [18], which comprises DAGs whose low-order d-separation statements are reflected as low-order conditional independencies in D. In similar vein, our approach focuses conditioning sets with cardinality of 1 or less ([2) < 1), allowing for more reliable and efficient causal discovery.
Continuous DAG Acyelicity Constraint. Pioneered by NOTEARS [36], recent score- score-based meth- ods provide an interesting alternative in DAG search by reformulating as a continuous optimization problem over weighted adjacency matrices. This transformation relies the development of DAG acyclicity regularizers, which quantify the degree of cyclicity in weighted adjacency matrices. Our work follows this practice and uses the log-determinant DAG regularizer introduced by DAGMA [4), defined as:
CDAG(W: s) = log dst - W) + d log
where s is a hyperparameter that controls the size of the valid domain of this function, denoted as "" { E gdxd 2 (W)), with p(W) representing the spectral radius of W. As established by Bello et al. [4], for W E W*, the condition LDAG(W, s) 0 implies that W maps to DAG.
3 Differentiable d-Separation Framework
Our core contribution lies in the transformation of d-separation, which is discrete and binary computed via combinatorial algorithms [8. 19, 32), into a continuously differentiable function over W E [0, 1]dxd, This transformation satisfies three central properties: (1) it exactly recovers conventional d-separation statements when evaluated on discrete DAGs A E {0, 1) (2) it provides principled probabilistic interpretation when W is viewed as parameterization of a distribution over graphs, and (3) it is differentiable function admitting gradient- based optimization
W.
The key insight grounding the score in percolation based graph connectivity measure. Specifically, we quantify how well a probabilistic graph structure (parameterized by W) aligns with conditional independence (CI) patterns observed in data, such as those measured through statistical tests yielding p-values, which pave the way to comprehensive set of differentiable objective functions that measure violations of low-order I-Mapness, Faithfulness, and Acyclicity constraints.
We achieve this differentiable framework through three key steps: (1) recasting d-separation as first-order logic (FOL) formulae operating on graph reachability (Section 3.1), (2) applying soft logic operators to transform these discrete statements into continuously differentiable functions with valid probabilistic interpretations (Section 3.2), and (3) combining these soft d-separation measures from the graph with conditional independence measures from data to construct objective functions that enforce low-order I-Mapness and Faithfulness constraints while maintaining acyelicity (Section 4).
3.1 FOL Formulation of d-Separation
Conventionally, d-separations are checked by discrete, combinatorial algorithms using data structures such as hash sets and double-ended queues [8, 19, 32). The inherently sequential and combinatorial nature of these algorithms presents significant challenges in developing continuous relaxations and enabling gradient-
optimization. To address this limitation, we introduce a First-Order Logic (FOL) framework that provides an equivalent characterization of d-separation computations in discrete settings while maintaining differentiability properties. This reformulation bridges discrete reasoning and continuous optimization paradigms, it allows us to leverage established techniques from probabilistic soft logic and fuzzy logic systems for end-to-end differentiation of the computation.
Specifically, we observe that the conventional notion of (low-order) d-separation and d-connection on any given discrete DAG A E (0, 1)xd can be precisely captured by FOL. formulæ based solely on graph reachability information. That is, given the graph reachability logical predicate
3

---

a?
b
a
backdoor path
frontdoor path
X
(a) Oth-order d-connection pattern: checking exis- tence of common ancestor
(b) Ist-order d-connection pattern A: (c) 1st-order d-connection pattern B: checking existence of backdoor path that does not involve =
checking existence of frontdoor path that involves
Figure Illustration the graph connectivity patterns checked by the low-order d-separation/d- connection FOL formulæ (Definition 3.1) to determine whether the given query ((I, y) or (x, y =)) d-separated or d-connected. The double-dashed line- = refers to a Oth-order d-connecting path, and the arrowed line and could be the query node x or y themselves.
refers to a directed path. In either pattern, the variable nodes
RA [d] + 10 defined as
Ra(a,y) = and only if ~e y in A.
We can define the FOL formulæ expressing the Oth-order (i.e. unconditional) and 1st-order (i.e. conditioning on one variable) d-separation and d-connection statements as follows:
Definition 3.1. (Low-order d-separation and d-connection FOL formulæ) Given discrete directed
graph A E {0, 1]dxd the Oth-order d-separation formula 5 [d] - {0, 1} and the 1st-order d-separation formula ( :dd3 + {0, 1} are defined as:
S (x, v) = Va [d]. - (Ra(a, Ra(a,
((1) (I, == 1) (x, ^ Va E [d] \ {=), 2 (2, v V -Ra(a,2))
(3)
v v(v E [d] (z). S0 (2.b) RA (b,
Equivalently, the O-th and 1-st order d-connection statements C/2 and c22 are:
x [z, y) 3a [d] R) (a,
(4)
2 (x) 00 (5, 3a [d] (=), 00 (I, A Ra(a,z)
(5)
(3b E id (2), 00 (v.b) RA(.)
To understand why these formulæ correctly capture d-separation and d-comnection when A is a DAG, consider first the Oth-order formulae ( (x, y) and C00) (ar, y). These formulæ essentially determine whether any node (including possibly x or y themselves) serves as common ancestor to both æ and y (e.g. the node a? in Figure la). If such common ancestor exists, it creates either chain structure (when a = x or a = meaning x is the ancestor of y or vice versa) or a fork structure (when a & x and a # y, and third node is the common ancestor). In either case, x and y are unconditionally d-connected, yielding C00) (x., y) = as expressed in Equation (4).
The 1st-order formulæ ((1) (z, y. 2) and (-11) (Ix =) involve more intricate logic but follow an intuitive two-step process. Take C(1) ( 2) for instance. first examines, via the first term (.0)) (z, y). whether and remain d-connected in the subgraph A_, where conditioning node removed effectively checking for any backdoor path that bypasses (Figure 1b). If no such backdoor path exists, the formulæ then evaluate whether there exists a valid frontdoor path through z. For this frontdoor path to d-connect x and when conditioning on z, z must a collider or a descendant of collider. This translates to checking whether there is any Oth-order d-connecting path that does involve z and that connects & to any ancestor of z or z itself, and similarly whether there such path connecting any ancestor of or itself. As illustrated in Figure Ic, this is checking

---

if the Oth-order d- d-connecting path 2 = == a? and the directed path a? - .... - 2 both exist the same time. The following theorem formalizes the correctness of these logical characterizations.
Theorem 3.2. For any DAG A with d nodes and any three nodes I. y, E
[d] that are distinct,
I IL. y if and only if ( (x, y) = 1, and x 11 V 2 if and only if (1)) (z, y 2) = 1. Similarly,
A
H A if and only if C (x, y) = L, and x 1 A V if and only if 1 (z, y 2) 1.
The d-separation and d-connection formulæ presented in Theorem 3.2 rely fundamentally on the graph reachability predicate RĄ. Crucially, we observe that RA itself can be expressed through FOL formulæ that operate directly on the adjacency matrix A, by taking a recursive form analogous to the generalized Bellman-Ford algorithm [3].
Definition 3.3. (Graph Reachability FOL formulæ) Given a discrete directed graph A {0, 1]dxd we define the reachability formulæ A (I, [d]2 + {0,1} for path lengths no greater than ! recursively as follows:
(0) (I, 1x y),
(6)
(2,
(z. v)
R R(l-1) (x, A Aus)) v R
((1-1) (x,
(z,1).
KEId
In practice, we take RA(x, y) == 1 (z. y) where is the number of nodes, since the longest directed path in DAG with d nodes has length d - 1. The correctness of these formulæ is articulated in the following lemma:
Lemma 3.4. For any discrete graph A E {0, 1)dxd with maximum directed path length and for all puir of nodes E [d].
y if and only if (1) (x, v) =1 1.
By integrating the formulæ in Definition 3.1 with those in Definition 3.3, we establish a complete framework that can systematically determine d-separation and d-connection relationships in DAGs using purely FOL operations, laying the foundation for the differentiable relaxation that follows.
3.2 Continuous Relaxation via Soft Logic
Having established d-separation and d-connection as FOL formulæ, we can now leverage the rich literature on probabilistic soft logic (PSL) and fuzzy logic to transform these discrete statements into continuously differentiable functions. These soft logic frameworks systematically relax Boolean operations (conjunction, disjunction, existential and universal quantification) into continuous func- tions that operate over the interval [0.1] via various t-norms and contormss-cottuuuuus analogs of logical AND and OR operations [2, 5, 17, 38J. In doing so, the logical semantics necessary to defined reachability (a percolation metric) are preserved while enabling gradient-based optimization and sampling.
In this work, we adopt the LogLTN framework proposed by Badreddine et al. [2], which implements the product t-norm (for logical AND) and max t-conorm (for logical OR) in the logarithmic space. Given set of m values [i(]i each v; E [0, 1], and their logarithmic representations {2/ [mm]
each x log(x;), the standard product t-norm and max t-conorm Tmx Om [0, 10,1)T _ + [0, 1]. and
LogLTN's logarithmic versions Tm, Öm I + R_,, where R_ {x R S 0), are:
Tm (m(fx))ie(m))= = II Om(lashie(m)) = maxfzi)ie(m)
1=1
(8)
Ōm (m((4)femm)) = (c+ log
Ei E Eia (iaiet/o-)))
Tm Tm((fa)ieemm) =
x
where C' = max(xaaaaa.... am//a) and a E (0, is the temperature controlling the approximation accuracy (the closer is to 0 the more accurate). Notably, the logarithmic t-norm T.,. exactly
"The proof' of Theorem 3.2 and all following theorems and lemmas can be found in Appendix A.
5

---

represents its standard counterpart, i.e., (.)()) = lo)( whereas the logarithmic t-conorm On implements the Log- Log-Mean-Ean Expnenaxlneetial operation that provides a lower bound approximation of log(Om(() (-)) with bounded error [2]:
log(On, (x)) a log(m) S Ōm(x') - log(Om(x)).
Our choice of the LogLTN operators motivated by two key considerations. First, a comparative analysis by van Krieken et al. [38] shows that the product t-norm provides the most stable and informative gradients for differentiable learning. while Badreddine et al. [2] shows that the Log- McanSum approximation of the max t-conorm in the logarithmic space further enhances gradient stability. Second, as we will demonstrate next, this specific combination of max- product operators consistently yields principled lower bounds on the quantities being estimated and allows us derive continuously differentiable functions, which we name the differentiable ceaeparatiooddccooneeee that compute lower bounds expected d-separation/ddcomneett DAGs parameters by the weighted adjacency matrix W. Specifically, we now transform the FOL formulæ in Definition 3.1 and Definition 3.3 into their continuous counterparts that operate on the weighted adjacency matrices W E [0, I] 1jaxd, using the LogLTN logical operators.
statements over distribution of
Definition 3.5. (Differentiable Graph Reachability and Unreachability Percolations) Given a weighted adjacency matrix W E [0, 1jaxd, we define the relaxed reachability and unreachabil- ity functions, RW (x, v), Ü
(1) (I, v) : [d[² + R_, for path lengths no greater than / recursively as
follows:
u), 1o U ( 2 (z, (o),)))); (10) w (x.y) log( # F)) Üw( v) = Tati ((02(UW) (1)))))) (z, log w
RW (z. v) = log(1(( = v))
( v) = Oaii
(z, W. Wwyl)loutg U (Ū "x,"((,,)), wll)
where we take (0g = -cc
Definition 3.6. (Differentiable d-Separation and d-Connection) Given a weighted adjacency matrix W [0, _jdxd, we define the Oth-order differentiable d-separation and d-connection function ((0) (I,, y), CW (z,w) [d]e + R__ and the sst-order differentiable d-separation and d-connection function 5 SW (z, y yl2,, CW (a, v12) [d + R_ as follows:
S (z. - Ta O2 (oW) (a, G (a.w)) )eeg)) (a :- Ou ((Ta (a. (a. v)) ).eli)
(W = (4.) Oz Te-s (000500 (2. (..))) besvee)) Fer (10 (.0) (0. hessvae)
(a == On CW Ji Oa
(Ts (20 (z. MW (Q )esoorcc) Oc ((7: (GW (v. (.. 1010. )
Notably, for the d-separation functions relying on graph unreachability, we directly derive unreacha- bility from the negation of reachability FOL formule rather than computing log(1 exp( R (I, This is essential for preserving the correct probabilistic interpretation, as we now demonstrate.
Consider W as representing parametric distribution over discrete graphs A. Namely, to sample A, we sample each edge x + y independently via Azy ~ Bern Wzy) (i.e. Bernoulli distribution). We denote the entire graph's 's sampling as A Bern W). We first notice that the differentiable graph reachability and graph unreachability functions on W given in Definition 3.5 always yield a lower bound to the expected reachability and unreachability on A respectively, when A ~ Bern (W). Intuitively, this is because the max operator provides a lower bound on the probability of a union of events (logical OR), while the product operator correctly computes the probability of an intersection of events (logical AND) when said events are mutually independent or provides a lower bound when they are positively correlated, which is precisely the case for overlapping paths in graphs needed in reachability (a percolation metric). As these operators compose recursively in our reachability computation, they preserve the lower bound properties, leading to the following lemma.
Lemma 3.7 (Reachability Percolation Lower Bound). Given a weighted adjacency matrix W E [0, l]axd for any 0 S d, and for any pair of nodes x, y E [d], we have
R:00) (5,y)) S log E BeBemW(W) [ (1, w)] and ( (x. v) S log E A. [ (Ul((,),))

---

The lower bound property of both our reachability and unreachability functions is critical-non note that naively computing unreachability as log(1 exp( (I, u))) would yield an upper bound instead. This consistency in providing lower bounds is essential for extending our theoretical guarantees to the differentiable d-separation and d-connection functions. As the following lemma demonstrates, these functions also preserve the lower bound property relative to their expected values:
Theorem 3.8 (Lower Bound on Expected d-Separation Statements). Given a weighted adjacency x, for any three nodes x, y, 2 E [d].
matrix W E [0, 1]dxd,
(: WW (a. y) I log E A- (W)
1 (z) v)
sW (z. y 2) S log A.w W)
((()) (z, y 2)
CW (r.y)= S log E. (AwBee(W))
(1, w)
CW (z, y 2) S log BA~Bem(W)
(z. y 2)
These lemmas establish that our differentiable functions provide lower bounds on the logarithm of expected d-connection and d-separation statements when viewing W as parameterizing a distribution over discrete graphs. This mathematical guarantee opens a principled pathway for gradient- based optimization by maximizing these lower bounds.
4 DAGPA (DAG Percolation Apartness)
We now introduce DAGPA, practical instantiation to demonstrate the validity of our differentiable d-separation scores (Section 3.2). This implementation employs certain heuristics and simplifications to facilitate optimization. The primary goal is not to claim DAGPA's 's optimality, but rather to provide empirical validation of the framework's 's potential for causal discovery and open avenues for future research exploring more sophisticated realizations of the core principles.
In DAGPA, we parameterize the weighted adjacency matrix W with parameters 0 Rdxdd via the sigmoid transformation W = ( We develop objective functions that encourage alignment statements implied by the graph structure W and conditional independence patterns in data. Rather than using binary decisions based on thresholded p-values, DAGPA directly incorporates raw p-values from statistical independence tests as soft CI labels from data. This approach makes the model naturally robust to uncertainty in CI testing, as gradient-based optimization can prioritize clear cases (p- (p-values near 0 indicating dependence or near indicating independence) over borderline ones.
between the d-separationn/-cooneeee
Our objective functions consists of three components: Acyclicity, "True Positive," " and "True Negative" losses A true positive occurs when the model predicts a high probability of d-separation for variable pairs showing high independence p-values in the data. Conversely, a true negative occurs when the model predicts high probability of d-connection for pairs showing low p-values (strong dependence). canstraint-based CI Statement Losses of DAGPA). Given dataset D with n samples and d variables, and the p-values (-(-) of some statistical conditional independence tests on all the Oth-order and 1st-order statements. Let Bo
Definition 4.1 (Multi-task
{(x,y) v,y y [d], x y) and In = (1,, y, 2) x. y. 2 E [d], > y. X + 2, y + 2}. The CI Statement "true positive (TP)" losses LTP. and LTP. the "true negative (TN)" losses LTN-0 and LTN.I, well as the log-determinant DAG acyclicity loss LDAG from DAGMA [4], defined on the model parameters 0 E Raxd, are as follows:
Lypo(0, D) = Σ
(1, (z. s)pp(a v).
LTp, D) =
Σ
5())
mij (a, =)pp(E,y =))
(xa))l
L1N D) = - Σ
(40)
016) (z. (M - p)(x,))) ()) =
Σ
C (z.y 2) - pp(a.y =)),
(mgn))ET
Lnsus0, (0,s) = log det (f)) + dlogs.
where Mo maxi, ( y). M1 = maxi, PD (z, 2)) is the sigmoid function, and is the hyperparameter controlling the domain of the log-determinant DAG loss of DAGMA [4].
We note that, assuming the p-values accurately reflect the true CI patterns, minimizing the "TP" and "TN" losses maximizes S and C' on the correct queries. Furthermore, since Š and C' are lower bounds of the expected d-separation and d-connection statements over Bern (Theorem 3.8), we are maximizing lower bounds the correct optimization direction. The following lemmas show the consistency of the losses and the validity of taking sum aggregation over individual CI queries.
We use these terms loosely, p-vuuu in practice may not perfectly reflect true CI statements.
7

---

DADMA ELaee11
(k-1 MFC
NOTEARS (hosieAn))
Modet's Predicted DAG
oth- and 1al-orderr di d-Segaration Statements
Ground Truch DAG GP
Mand??
,
oF 8
A4.C
Ale0
1F D
X
C14
84 D A B. D
CD.D g PIA DIA
Number of
FN
CL-OC
0.8025
Lew wow-rreee Stateesien MCG (Cta6ic))
PP 6 TN -N=40
Beller
Figure 2: The Low-order CI Statement Matthews Corre- lation Coefficient (CI-MCC) metric measures alignment between the low-order d-separation statements predicted by the causal discovery model, and those implied from the ground-truth DAG (symmetric statements considered).
The empirical CDF CI- MCC on Sachs [28]. The closer the curves to the bottom-right (2), the better the performance. DAGPA achieves the best performance compared to all baseline methods.
Figure 3:
Lemma 4.2 (Consistency of Multi-Task CI Losses). Given a faithful data D as the sample size n + oo and the Eap temperature a + 0+ (Equation (8)). the optimal DAG A achieves the minimum value on each of Lre-o, CTP-I, LrN-O, LTN-1, and DDGG
Lemma 4.3 (Mutual Independence of Low-order CI Statements). All Oth- and Ist-order CI statements over the variable sets Iu = {(x,)) y) y. y E [d]. > yt and I = {(x, y, x. y. E [d], x > y. * 2.3 y # 2) mutually independent. In other words, none of these CI statements can be implied from any other of these CI statements via the graphoid axioms [24].
In Appendix C, we provide the complete algorithm for DAGPA along with detailed descriptions
of additional heuristics employed to address optimization challenges. including: PCGrad [43] for resolving conflicting gradients in multi-task learning. Discrete Langevin Proposal (DLP) [50] for efficient exploration of the weight space, and an optimization-time DAG selection score based on CI statement measures from the data.
Experiments
This section presents an empirical evaluation of DAGPA's ability to discover DAGs whose separation statements are consistent with those derived from the underlying causal structure that generated the data. Our objective is not to test superiority over baselines in terms of exact structural recovery (that is left to Appendix F), but rather to demonstrate that DAGPA successfully achieves its intended purpose: learning DAGs that entail the same low-order conditional independence statements as the ground truth. The experimental results indicate that our differentiable d-separation framework offers promising direction for causal discovery.
Evaluation. To evaluate causal discovery from CI statement alignment perspective, propose a metric that aligns with this objective: the Low-order CI Statement Matthews Correlation Coefficient (CI-MCC), illustrated in Figure 2. Given a model's predicted causal graph (which may be a DAG, PDAG, CPDAG, or k-essential graph), we derive the set of all Oth- and sst-order d-separation statements implied by the structure. We then compare these against the corresponding d-separation statements derived from the ground-truth DAG that generated the data. Treating this as a binary classification task, we compute the Matthews Correlation Coefficient, which ranges from [-1, 1] with higher values indicating better alignment between model- peeeel-predittd and ground- truth CI statements. We select MCC over alternatives like Fl score due its robustness to class imbalance critical consideration since the true ratio of independence to dependence relationships varies widely across different causal structures and is typically unknown in real- sal-world settings. For completeness, results using conventional causal discovery metrics, such as graph structure F1 and Structural Hamming Distance (SHD) scores, are provided in Appendix F.
Baselines. We compare DAGPA with representative causal discovery methods spanning different paradigms: constraint-based approaches including PC [34] and k-PC [18]; score-based methods such
8

---

INSERS BNGINA, (LANN]] INGMA (Ncelioor)
DAOPA (D
686
NOTRARS (Evwa))
S-11800
NOTEARS (Nieirenn
DAGIAA (LaWAR)
NOTEARS (hoelanne)
Ciesteent
Eayinca CDP
Destitt
Besiaer
mmportal
Lowerrdee C1 Shatement MCC (CL-MEC)
Low-ander CI Stabemin: MEC ICEAOCO)
Betber
Belber
ER d 10, 2, 100
(b) ER 10, n, 10000
Figure 4: The empirical cumulative distribution functions (CDFs) of the low-order CI Statement Matthews Correlation Coefficient (CI-MCC) on two synthetic binary settings. The closer the curves to the bottom-right ( (2)), the better the performance. Our method achieves on-par performance to the constraint-based baselines in high-sample regime (Figure 4b) while demonstrating much better robustness in low-sample regime (Figure 4a).
as GES [6]; and continuous optimization approaches including NOTEARS (linear and nonlinear) [36] and DAGMA (linear and nonlinear) [4]. This selection provides a comprehensive comparison land- scape across methodological families. Detailed baseline configurations are provided in Appendix D.5.
Experimental Settings. We evaluate methods on both synthetic and real-world datasets. For synthetic data, we generate binary data requiring no normalization. We follow a protocol similar to Kocaoglu [18], generating binary datasets from randomly constructed DAGs. We use both Erdós-Rééns Rényi (ER) and Scale-Fre Free (SF) graph models via the pyAgrum package [10], with varying complexity parameters: nodes d E {10, 50), edge-to-node ratio E {2. 4), and sample sizes {100, 1000, 10000, 100000]. For each configuration, we generate 10 datasets. All constraint- based methods employ Fisher-z tests [11] (with potential of other CI tests as alternative; details in appendíx), for computing conditional independence p-values, ensuring fair comparison. For real-world validation, we use the Sachs dataset [28], a benchmark protein signaling network with 11 variables and a known tround-truth causal structure derived from experimental interventions. This dataset represents a challenging real-world case where the underlying causal mechanisms are likely nonlinear and complex. We employ the standard discretized version preprocessed via the Hartemink discretization method, which converts continuous protein concentrations into 3-level categorical variables (low, average, high) while preserving dependence structure.
We note that our Cl-constraint- based approach offers inherent robustness to preprocessing artifacts that can affect score-based methods: standard CI tests (e.g. Chi-squared for discrete data, Fisher-Z. for continuous data) are invariant to scaling transformations, whereas likelihood- objectives can exploit variance differences as spurious causal signals [27]. This makes our framework particularly robust across different preprocessing pipelines.
Results on Sachs. To evaluate performance on real-world data with potentially complex underlying mechanisms, we test DAGPA on the Sachs dataset [28], a benchmark protein signaling network with 11 variables and 853 samples. Figure displays the empirical CDF of CI-MCC scores across methods. DAGPA demonstrates dramatically superior performance, achieving CI-MCC scores of 0.75-0.98, while all baseline methods cluster between 0.3-0.65. For our method, we select the 10 DAGs with highest DAG selection scores to reflect the uncertainty inherent causal discovery from real data. For baselines that depend on random seeds, we run each method with 10 different initializations to ensure fair comparison. The substantial performance gap highlights DAGPA's 's ability to effectively discover causal structures in real-world data where the underlying data- eeta-generaatig process is complex and potentially nonlinear. This confirms that our differentiable d-separation framework. combined with gradient-informels discrete sampling, provides robust causal discovery capabilities that generalize well beyond synthetic settings to real-world applications where ground truth is available but unknown during learning.
Results on synthetic binary data. Figure 4 presents the empirical cumulative distribution functions (CDFs) of CI-MCC scores across methods for different sample sizes. For each baseline, we include

---

one result per data instance in the batch (10 results total), while for DAGPA, we select the 5 discovered DAGs with the highest DAG selection score (Appendix C.3) per instance, yielding 50 total evaluations.
In low-sample settings (n=100, Figure 4a), DAGPA demonstrates superior performance compared to most baselines and performs on-par with DAGMA-Nonlinear [4], with CDF curves positioned toward the bottom-right indicating consistently higher CI-MCC scores. The advantage of DAGPA is particularly pronounced when compared to constraint-basee methods (PC, k-PC), which struggle in low-sample regimes due to their reliance on hard conditional independence statements derived from uncertain p- palluess. Most score- score-based methods, including DAGMA- -lnnear and both NOTEARS variants [36], also underperform in this setting, likely due overfitting to limited data.
As sample size increases to n=10000 (Figure 4b), most methods improve as expected, with DAGPA maintaining competitive performance alongside constraint-basee approaches. This pattern holds consistently across different graph structures and dimensions (see Appendix F), confirming that DAGPA effectively addresses RQI by exhibiting strong robustness low-sample settings while maintaining competitive performance when data is abundant.
6 Related work
Constraint-based Methods. These methods discover causal structures through conditional indepen- dence (CI) testing. The PC algorithm [33, 34] produces CPDAGs representing Markov equivalence classes but struggles with small samples due unreliable high-order CI tests. Recent advances like LOCI [42] and k-PC [18] demonstrate that low-order CI statements can sufficiently identify causal structures in many practical settings. Our work builds on these insights but uniquely transforms discrete CI decisions into a differentiable learning framework, further enhancing robustness in the small sample regime.
Score-based Methods. These methods optimize a score function over the DAG space. Traditional approaches like GES [6] use greedy search with cntemaation-thhorette criteria but face scalability challenges due to the discrete DAG space. NOTEARS [51] introduced a breakthrough with differ- entiable acyclicity characterization, enabling gradient-based optimization weighted adjacency matrices and spurring extensions like nonlinear variants [52] and improved acyclicity constraints in DAGMA [4]. Our framework builds upon these continuous optimization advances but redirects the objective toward maximizing agreement with CI constraints rather than fitting functional models. This creates a hybrid approach combining strengths from both paradigms and opening a new avenue for causal discovery.
A more comprehensive discussion related work can be found in Appendix E.
7 Conclusions
We introduced novel causal discovery framework that bridges constraint-based and gradient-based continuous-optimizatio methods through differentiable d-separation. By proposing percolation- based differentiable d-separation scores, we enabled gradient-based learning of causal DAG structures from conditional independence patterns in finite datasets. Our model instantiation, DAGPA. demon- strates competitive performance across synthetic and real-world settings, with strong robustness in low-sample regimes where traditional methods often struggle, while maintaining competitive performance as data becomes abundant.
Limitations and Future Work. Several limitations in current approach present opportunities for future research. The reliance of DAGPA on p-values to measure both conditional dependence and independence from data can be improved, as using constant minus the p-value to assess dependence strength is a simple heuristic. Fortunately, our framework is measurement-agnosti and supports alternatives, which recommend future work investigate. In addition, our method focuses on only low-order conditioning sets and assumes no confounding variables in the data. Future directions include extending the framework handle latent confounders, efficiently incorporating higher-order CI statements, and simultancously improving scalability. Finally, as suggested by reviewer, a particularly promising avenue integrating our differentiable di d-separation framework with acore-based methods like NOTEARS and DAGMA, enabling causal structure learning from both data likelihood signals and CI statement fitness signals (case studies in Appendix G. These advancements would further bridge constraint- and score- score-based paradigms in causal discovery.
10

---

Acknowledgment
The authors thank Ruqi Zhang and Patrick Pynadath for the insightful discussions and valuable feed- back Bayesian sampling methods. This work was funded in part by the Ford Motor Company under the Purdue-Ford Alliance, the National Science Foundation (NSF) awards CAREER IIS-194333 INS-2212160 and CAREER 2239375, IIS-2348717, Amazon Research Award, AnalytiXIN, Adobe Research, Intuit, the Wabash Heartland Innovation Network (WHIN). Walther Cancer Foundation and Purdue Institute for Cancer Research (PICR) (Grant P30CA023 168). Computing infrastructure supported in part by AMD under the AMD HPC Fund program and PICR. Any opinions, findings and conclusions recommendations expressed in this material those of the authors and do not necessarily reflect the views of the sponsors.
References
[1] Yashas Annadani, Nick Pawlowski, Joel Jennings, Stefan Bauer, Cheng Zhang, and Wenbo Gong. Bayesdag: Gradient-based posterior inference for causal discovery. Advances Neural Information Processing Systems, 16:117881766 2023.
[2] Samy Badreddine, Luciano Serafini, and Michael Spranger. logitn: Differentiable fuzzy logic in the logarithm space. ArXiv, abs/2306. 11e233066114,,, 14546, 2023.
[3] John Baras and George Theodorakopoulos. Path problems in networks. Synthesis Lectures on Communication Networks, 3:1-77, 2010.
[4] Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. Dagma: Learning dags via m-matrices and a log-determinant acyclicity characterization. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 8226-8239. Curran Associates, Inc., 2022.
[5] Matthias Bröcheler, Lilyana Mihalkova, and Lise Getoor. Probabilistic similarity logic. In Proceedings of the Twenty-Sixth Conference on Uncertainty in Artificial Intelligence, pages 73-82, 2010.
[6] David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 30Nov):507-554, 2002.
[7] Gabor Csardi and Tamas Nepusz. The igraph software package for complex network research. InterJournal, Complex Systems:1699 11 2005.
[8] Adnan Darwiche. Modeling and reasoning with Bayesian networks. Cambridge university press, 2009.
[9] Iddo Drori, Anant Kharkar, William R Sickinger, Brandon Kates, Qiang Ma, Suwen Ge, Eden Dolev, Brenda Dietrich, David P Williamson, and Madeleine Udell. Learning to solve combinatorial optimization problems on real-world graphs in linear time. In 2020 19th IEEE International Conference on Machine Learning and Applications (ICMLA), pages 19-24. IEEE, 2020.
[10] Gaspard Ducamp, Christophe Gonzales, and Pierre-Henri Wuillemin. aGrUM/pyAggrum Toolbox to Build Models and Algorithms for Probabilistic Graphical Models in Python. In 10th International Conference on Probahilistic Graphical Models, volume 138 of Proceedings of Machine Learning Research, pages 609-612, Skørping, Denmark, September 2020.
[11] Ronald Fisher. Frequency distribution of the values of the correlation coefficient in samples from an indefinitely large population. Biometrika, 10(4):507-521, 1915.
[12] HI. Frisch and JM Hammersley. Percolation processes and related topics. Journal of the society for industrial and applied mathematics, 11(4):894-918 1963.
[13] Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery methods based on graphical models. Frontiers in genetics, 10:524, 2019.

---

[14] Isabelle Guyon, Constantin Aliferis, Gregory Cooper, André Elisseeff, Jean Philippe Pellet, Peter Spirtes, and Alexander Statnikov. Causality workbench. In Causality in the sciences. Oxford University Press, 2011.
[15] W Keith Hastings. Monte carlo sampling methods using markov chains and their applications. Biometrika, 57(1):97-109, 1970.
[16] Barry D Hughes. Random walks and random environments. Oxford University Press, 1996.
[17] George Klir and Bo Yuan. Fuzzy sets and fuzzy logic, volume 4. Prentice hall New Jersey, 1995.
[18] Murat Kocaoglu. Characterization and learning of causal graphs with small conditioning sets. Advances Neural Information Processing Systems, 36, 2023.
[19] Daphne Koller and Nir Friedman. Probabilistic Graphical Models: Principles and Techniques. The MIT Press, 2009, ISBN 0262013193.
[20] Bo Liu, Xingchao Liu, Xiaojie Jin, Peter Stone, and Qiang Liu. Conflict-averse gradient descent for multi-task learning. Advances in Neural Information Processing Systems, 14:18888-18880, 2021.
[21] Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller. Equation of state calculations by fast computing machines. The journal of chemical physics, 21(6):1077-1099 1953.
[22] Aviv Navon, Aviv Shamsian, Idan Achituve, Haggai Maron, Kenji Kawaguchi, Gal Chechik, and Ethan Fetaya. Multi-task learning bargaining game. In International Conference Machine Learning, pages 16428-16446. PMLR, 2022.
[23] Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and dag constraints for learning linear dags. H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 17943-17954. Curran Associates, Inc., 2020.
[24] Judea Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988. ISBN 1558604790.
[25] Judea Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press, USA, 2nd edition, 2009. ISBN 052189560X.
[26] Judea Pearl and James M. Robins. Probabilistic evaluation of sequential plans from causal models with hidden variables. In Conference Uncertainty in Artificial Intelligence, 1995.
[27] Alexander Reisach, Christof Seiler, and Sehastian Weichwald. Beware of the simulated dag! causal discovery benchmarks may be easy game. Advances in Neural Information Processing Systems, 34:27772-27784, 2021.
(28) Karen Sachs, Omar Perez, Dana Pe'er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(5221): 523-529, 2005.
[29] Gideon Schwarz. Estimating the dimension of a model. The Annals of Statistics, 6(2):461-464. 1978.
[30] Marco Scutari and Jean-Baptiste Denis. Bayesian Networks with Examples in R. Chapman and Hall, Boca Raton, 2nd edition, 2021. ISBN 978-0367366513.
[31] Ozan Sener and Vladlen Koltun. Multi-task learning multi-objective optimization. Advances in neural information processing systems, 31, 2018.
[32] Ross D Shachter. Bayes-ball: Rational pastime (for determining irrelevance and requisite information in belief networks and influence diagrams). Proceedings of the Fourteenth conference on Uncertainty in artificial intelligence, pages 480-487, 1998.
12

---

[33] Peter Spirtes and Clark Glymour. An algorithm for fast recovery of sparse causal graphs. Social science computer review, (1):622-7, 1991.
[34] Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, prediction, and search. MIT press, 2000.
[35] Eric V Strobl, Kun Zhang, and Shyam Visweswaran. Approximate kernel-based conditional independence tests for fast non panametraee causal discovery. Journal of Causal Inference, 7(1): 201 20180017, 2019.
[36] Xiangyuan Sun, Oliver Schulte, Guiliang Liu, and Pascal Poupart. Nts-notears: Learning nonparametrie dbns with prior knowledge. In International Conference on Artificial Intelligence and Slutistics, 2021.
[37] Ioannis Tsamardinos, Laura E Brown, and Constantin F Aliferis. The max-min hill-climbing bayesian network structure learning algorithm. Machine learning, 65:31-78, 2006.
[38] Emile van Krieken, Erman Acar, and Frank van Harmelen. Analyzing differentiable fuzzy logic operators. Artificial Intelligence, 302:103602 103602, 2022.
[39] Matthew J Vowels, Necati Cihan Camgoz, and Richard Bowden. D'ya like dags? survey structure learning and causal discovery. ACM Computing Surveys, 55(4):1-36, 2022.
[40] Zirui Wang and Yulia Tsvetkov. Gradient vaccine:
Investigating and improving multi-task optimization massively multilingual models. In Proceedings of the International Conference on Learning Representations (ICLR), 2021.
[41] Dennis Wei, Tian Gao, and Yue Yu. Dags with fears: A closer look at continuous optimization for learning bayesian networks. Advances in Neural Information Processing Systems, 33:3895- 3906, 2020.
[42] Marcel Wienöbst and Macicj Liskiewiez. Recovering causal structures from low-order con- ditional independencies. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 10302-10309, 2020.
[43] Tianhe Yu, Saurabh Kumar, Abhishek Gupta, Sergey Levine, Karol Hausman, and Chelsea Finn. Gradient surgery for multi-task learning. Advances in Neural Information Processing Systems, 33:5824-5836, 2020.
[44] Alessio Zanga, Elif Ozkirimli, and Fabio Stella. survey on causal discovery: Theory and practice. International Journal of Approximate Reasoning, 151:101-129, 2022.
[45] Bin Zhang, Chris Gaiteri, Liviu-Gabriel Bodea, Zhi Wang, Joshua McElwee, Alexei A Podtelezh- nikov, Chunsheng Zhang, Tao Xie, Linh Tran, Radu Dobrin, et al. Integrated systems approach identifies genetic nodes and networks in late-onset alzheimer's disease. Cell, 153(3):707-720, 2013.
[46] Jiji Zhang. Causal inference and reasoning in causally insufficient systems. Philosophy of Science, 75(5):960-970, 2008
[47] Junzhe Zhang and Elias Bareinboim. Fairness in decision-making the causal explanation formula. In AAAI Conference on Artificial Intelligence, 2018.
[48] Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kernel-based conditional independence test and application in causal discovery. In Proceedings of the 27th Conference on Uncertainty in Artificial Intelligence (UAI), pages 804-813, 2012.
[49] Kun Zhang, Bernhard Schölkopf, Krikamol Muandet, and Zhikun Wang. Domain adaptation under target and conditional shift. In International conference on machine learning, pages 819-827 Pmlr, 2013.
[50] Ruqi Zhang, Xingchao Liu, and Qiang Liu. langevin-like sampler for discrete distributions. In International Conference on Machine Learning, pages (6375-26396 PMLR, 2022.
13

---

[51] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning. Advances in neural information processing systems, 2018.
[52] Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning sparse nonparametric dags. In International Conference on Artificial Intelligence and Statistics, pages 3414-3425 PMLR, 2020.
[53] 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, 2024.
14

---

A Proofs
Definition 3.1. (Low-order d-separation and d-connection FOL formulæ) Given a discrete directed graph A E {0,1) the Oth-order d-separation formula s(0) d-separation formula 5 [d] + {0, 1} are defined as:
A [d]2 + {0, 1 and the Ist-order
S (z, y) va [d]. - (R (a(aa,)))
(2)
(a,y =)= 2 (x, (a E [d] \ ( 22, V= Ra((a,)))
(3)
V [d| (=). (v.b) (RA ((,(
Equivalently, the O-th and 1tst order d-connection statements /¹4²/ and c are:
c (z. y) Fa [d]. Ra x) ya,
(4)
(Z,y =)==
00 A_, (2, v (33a E [d] \ {2), (0 A_, A Ra(a,z)
(5)
[db R (), ^ Ra (b. (
Theorem 3.2. For any DAG A with d nodes and any three nodes I, y. E [d] that are distinct, IL. A if and only if 1 (x, y) 1, and 11. A V and only if S (x, y 2) 1. Similarly, x V A if and only if C (x. y) = 1, and x V1 A if and only if (1) (x, 2) = 1.
Proof. First, we note that the d-separation formulæ ((0) (( and the d-connection formulæ C0 . 0. are the negation of each other and can be obtained from each other via De Morgan's 's law. Thus, the statement "ær IL y if and only if ( (3c, v) = 1 is equivalent to "a V. A if and only C (z, y) = 1", and vice versa, the statement "<x IL A 2 if and only if s(1)
A (z, y 2) = 1" equivalent "x 1/ A | if and only if c/ (z. y 2) 1." Thus, we proceed to prove the statements involving the d-connection formulæ A c(())
(0)
Part 1: x 1 A if and only if C100) (z. (2 y) 1
We show the direction VL AU -
(2, 1.
æ JL_ A y. i.c., x is d-connected to y with an empty conditioning set, then there exists a path between and that does not contain any colliders, because otherwise the path would be blocked. Denote the sequence of nodes in this path p = (p...... PPm) with p1 = x and pm = y. In other words, for any three consecutive nodes a, b, c p. have either the chain structure - b or a + b- c, or the fork structure a +- - c.
Because there is no collider structure in p. we will show that there must be most one fork structure in p. To see this by contradiction, consider p contains more than one fork, and let Pt-2 + Pt-1 - Pt and p. + Pr+1 + Pr-2 be two forks (1 < r) in p such that there is no other fork between Pt and pr. Then, either there are some colliders between pl and pr, or pt and p, are connected by some directed paths, i.e., either Pl Pr p, Pt. If pt we Pr, then p. is a collider, and vice versa if pr ~ Pt. then pi collider. Hence, in all scenarios, p will have at least one collider, contradicting the fact that p should not have any collider.
Now, we discuss case by case based on whether p has no fork, or it has one fork. Figure la illustrates the structure such path p.
Case 1: Suppose there is no fork on p: If there is no fork all, then p is directed path. Without loss of generality, assume = x - ... - y. Then, let a = x, we have a ~ x. (a node always
'The proofs of Theorem 3.2 and all following theorems and lemmas can be found in Appendix A.
15

---

naively reachable from itself) and a ~ y. Hence, (Ra(a,z) A Ra v)) = 1, and therefore the formulæ CAA (z. y) Ea E [d], RA 2) A Ra(a, y) is 1.
Case 2: Suppose there is one fork on When there is one and only one fork, it means the center node of the fork is the common ancestor of both x and y. In other words, let be this center node, then we have ~.. x and a ~- y- Thus, similarly, we have (Ra(a, x) A Ra (a, v)) = 1, and therefore the formulæ C (z. v) 3a [d], (a(a, x) A (a v) is 1.
Thus, we have shown that r 1 A = 1 (z, y) =1
Now we show the other direction (1) A (z, y) = = r x 1 AV.
and a an y. Let a be such a node. Thus, a* is a common ancestor of a and Let Pa+ = a + ... - x be the directed path that starts at and ends x, and similarly let Py = a + + y. Then, the path form by the sequence of nodes in PN and py: i... p (I.....,a, a....... y), is clearly path consisting of at most one fork structure with other wise chain structures. When a or a y. p is a directed path without fork, and otherwise has one fork. In either case, is a d-connecting path between and y. making r 1 A
Given C/ (x,
we have 3a [d] such that a
Part x V4 A if and only if C [) =
We show the direction 14 A - ( =1)=1
We are given x + A 2, i.c., x and are d-connecting when conditioning on z. Thus, there must exist a non-empty set of d-connecting paths P between and y that is not blocked given 2. We proceed with a case-by- casee discussion depending on whether all of these paths contain colliders.
1. First, consider the case where there is some path in P that does not contain any collider. Let P E P be such a path. Then, 2 must not be on p. because otherwise, since is not a collider and it is being conditioned, p would be blocked and therefore p P.
Since p does not include it remains unchanged in the subgraph A-2, in which node is removed from the graph A, and thus it is still a d-connecting path in A- Hence, because of p. we have x and y are d-connecting in the subgraph A_g. In other words, we have Ca(0) (x,
1. Since c(O) (xz, 1/) is the one of two terms from the outmost disjunction operator V in the formula for (.)) (Ir, y =)) we have (00 (xx y) = renders c/2 (az, y 2) 1.
In Figure 1b, this pattern is illustrated via the "backdoor path" structure.
2. Otherwise, all paths P contain colliders. We first show that this implies that we can find path p in P that has only one collider.
To see this, pick any path p E P. If p indeed has only one collider, then we are done. Hence, we consider the non-trivial case where p has more than one collider. We further divide the cases depending on whether node is included in path p.
(a) If is included in p. then z must itself be a collider, because otherwise p would be blocked by 2.
Now, let 2' be another collider in p. Then, must be descendant of s', because otherwise would be blocked by 2' even though 2 is in the conditioning set. In other words, is reachable from z' with directed path, i.e., z' 2.
Using this knowledge, we can construct the following alternative path par that still has z as a collider but converts to a non-collider. Specifically, without loss of generality, we assume the path p consists of the following sequence of nodes: p (z..... pur, p',Pr..... Pt, 2, Pr..... ..) for some indices l, r' with collider structures pv + z' <- Pr and pl responding directed path + qi + + qm + 2. We can obtain the new path as follows:
+- pr. We have shown that z' ~4 2, and denote the cor-
Ps' (z...., pr, z' a1, ....... Ims 2, Pr. ...r....,?7 y).
where we replace the middle section sub-path (z', .r...... Pt, 2) in p to (2', (z',q),... qq..... ......... Ans to form Ps'
16

---

Now, note that in Par, z' is no longer collider but forms a chain structure, because Pr + 2' + Q1- Nevertheless, is still collider, because 9m - + Pr- Furthermore, Pz' is still d-connecting path. This is because the sub-paths (Ir,... " x²) and (z, Pr.pr....:)) V) are d-connecting conditioning on because they are parts of the original p which is a
d-connecting path, and the "new" sub-path (z', f1.-... Im: 2) by definition is d-connecting conditioning on 2 because it is a directed path. Hence, taken together, P is d-connecting given so Par E P.
Thus, by iteratively constructing these alternative paths P for any colliders 'n p that is not 2,
arrive at path p" E P that has the only collider.
(b) Now, suppose is not included in p and has more than one collider. Let z' and =" be any two such colliders, and without loss of generality, let the path consists of the sequence of nodes (I. ..... .P' z', .Pr'...... ,ir z", -p, ..... y). Similarly, both 2' and 2" must be ancestors of z, because otherwise one of them would block the path. Hence, we know that z' ~* and z" ~n 2. Let the directed paths be a' + 91
+ 9m + and O2 7 respectively. We can construct the following new path:
z" + 01
P2' (z..... .r. x/, ........ 9m 2, ................................ .1, ........)..)..
Here, note that both z' and z" are no longer collider, because pv + z' q1 and On < z²" + Pr. Nevertheless, 2 is a collider since Qm + 2 + Oz. Furthermore, 2= is still a d-connecting path conditioned on z. As have established in the discussion of the previous case, that the sub-paths "c,,...P,,
and (=", ,rr:...,,.) y) are d-connecting. Since (z'. Sz.91...... 9m 2) and (z, Oz..... 01, 2²) are directed paths, we have that the sub- paths (I.... In. 2) and (2, ........). v)) are also d-connecting. Taken together, the whole path is d-connecting given 2, i.c., P=',z" E In other words, we have converted both colliders 2' and z" in the original p into non- colliders while introducing 2 as a new collider. This means that the resulting path p satisfies the initial condition of path discussed in the previous case. Hence, by further following the procedure described the previous case, we can eventually arrive at a path p" E P that has z as the only collider.
Now that we know P will necessarily include a path p with only one collider, denote this collider 2", and we know either 2" = or 2" is an ancestor of 2. In other words, we have 2. Hence, we can denote p" as a sequence of nodes: p' = (I....,Pi ,i, p,, .....) y) for some indices 1, r. where P? + 2 t p,,
Now, let a = x if there is no fork structure in the sub-path (z...., =*) of p", otherwise let a be the rightmost fork in (z...., (i.e. there is no other fork in the sub-path (a..... 2*) 2²) of p*). Then, we have ~ 2n because there is neither any fork nor any collider in the sub-path (a...., .....,,,), and we know has constituent edge pi + 2" with an edge direction pointing towards =². Moreover, because 2² ~4 we have ~ 2, and equivalently (A(a, 2).
On the other hand, the sub-path (x...., a) is a d-connecting path conditioning on because its super-path p* E P is d-connecting path conditioning on z, does not have any collider because =" is the only collider in p" . and does not include z. Hence, in the subgraph A_; where is removed from the graph A. the sub-path (I.... a) remains unchanged and still a d-connecting path, and in this case without conditioning on z. Therefore, x and a is d-connecting in A_g. i.c., x H A_, a. Hence, we have found a node a such that (rr K A_, a) (a ~y 2). or equivalently, we have the following:
(3a
[d] \ (2), UL A_. (),(,,) 1.
Similarly, let b = y if there is no fork structure in the sub-path (z, (z..... y) of p", otherwise let b* be the lefimost fork in (2,..., v). We then have also have the following:
3b E [d] \ (2), (y 4 A_, b) A (A(b, =.
Combining everything, we fulfill second term in the outermost disjunction operator in the FOL formula and thus have (02 (z, v 2) = 1.
In Figure 1c, this pattern is illustrated as the "frontdoor path" structure.
Now we show the other direction C/ (z,vl v 2) = x If A 2.
17

---

If T (2, = 1, then we have
(0) C (2.11) =1
we have
((3a e [d (=), C2 (x) A Ra(a,))) a (36 E [d] fz), 22, (v.b) Ra(0,,))) 1.
We discuss it case by case.
1. If it is the first case, C (I, 1, then there exists d-connecting path p between and in the subgraph A-2. 2- p then does not include 2, and must not include any colliders. Hence, back in the original graph A. is a "backdoor" path that does not involve that renders æ and y d-connecting. Thus, we have 1 A y and also x 1f A
2. We now consider the other case. Since the left-hand-side term is a conjunction (A), both term in the conjunction must be true, meaning there exists some nodes and b satisfying either term respectively. Taking a step further, this means that all of the four terms C A_ (z, a). RA (a, z).
C-0 A_. (v, b), and RA (b, 2) are true.
Using similar reasoning mentioned in the previous case, C (x), a) implies that we have d-connecting path between x and a and that x 4 a 2. Similarly, C(O) (v. b) implies that y M. b Furthermore, because (a =) and (b, =) = 1, that is, ~ 2 and b ~ the path between a and b, formed by the corresponding directed paths that renders b ~ respectively, is a d-connecting path conditioning z, serves the only collider.
2 and
Hence, we have a 1K A bb] Hence, by transitivity of d-connection statements, we have that x V. A finishing the entire proof.
Lemma 3.4. For any discrete graph A E {0, 12dxd with maximum directed path length and for all pair of nodes E [d].
if and only if R (x, =1
Proof. We prove this lemma by induction on the maximum directed path length (.
Basecase: Suppose the maximum directed path length in A is 0. This means there is path and A an empty graph. Thus, for all x [d], x ~ and for all x, y [d] with x * x + y.
Induction: Suppose the statement in the lemma holds true with a maximum directed path length of | 1. We now prove the statement for
We prove the direction x ~ MAY =
Given that A y. let p be such directed path starting from and ending y. Let be the second last node in p with an edge u - y, or equivalently, An,y = 1. Then, we also have x A
Since the path from x to u will have length smaller than or equal 1, using the induction hypothesis, we have 12(1-1) (x, u) 1.
Thus, have R11 (11-1) u) Au,y and consequently Vued R satisfying the first term in the outmost disjunction in the recursive formula for R!). rendering R(l) (z, 8) =
R((-1) (zx, u) A Au,y
Now we prove the other direction ( = = CMAA y.
Given R20 (x, y) 1, there possible scenarios: that Vued R((-1) ().1))( (xx, u) A Auy 1 or that
(1l-1) (I.z) 1. we discuss case by case.
18

---

A 1. Suppose V 12 (2r, (z, u) A Auy = 1, then there exists at least one node such that
AAA Rll-1 (I, A Au.. = L, meaning both - the induction hypothesis, x AA u with maximum path length / - 1. Furthermore, have u - is an directed edge. Thus, there exists directed path from to y, i... ~M with maximum path length I.
(x, u) and Au.y = 1. This means that, by
2. Suppose H(l-11(a
(I. y) = 1, then by the inductive hypothesis, x ~ y with a maximum path length / - I, which also makes the statement true under the setting of maximum path length l.
Thus, in either case, we have x ~ A y.
The proof of Lemma 3.7 and subsequent Theorem 3.8 requires the following lemma showing the lower-bound nature of LogLTN's [2] t-norm and (-conorm logical operators. We state this lemma and its proof as follows.
Lemma A.1. Given the t-norm and 1-conorm operators of LogLTN [2], which we restate as follows,
I'm ((aA)eelm)) > (mm((r()ss(m)m =a (c+ log
m
Ei= edlla-0
m
Let X1. 12..2..... Xm be Bernoulli random variables that are mutually independent or positively correlated. Denote their probabilities as P1, P2...... .2.....pm.. respectively. Then, we have the following:
Tm ((log(p))) log }}} log P (x1 A Xm)
(m(llo((p)).... Log(,............ }} S log (x1 V xm)
Proof. We first show the inequality involving the t-norm T We first observe that T,, is the logarithm of the product t-norm, meaning
exp(TI((flg(pn)) (flog(p1). log(pm))
II i=1
Pr-
Now, since x; are mutually independent or correlated, we naturally have II P: P (x1 A ... X.). This is because for any two random variables and Xj.
S
((x; A Ax) Exx; x,] =
+ Cov(x;, x,) PPP + Cov(xi, 2 PiPj:
where the last inequality sign holds because Cov(x,, x,) 2 0 for mutually positively correlated or independent random variables.
Thus, we have
(Llog(p)).... (mllll(()) log_,,,,l.)m))) < log P (xl A... xm)
< II Pi < P (x1 A... Xm)
For the t-conorm Öm. the inequality holds because Ôm is lower bound of the logarithm of max, and max is then also lower bound the union of random variables. Specifically, as shown in the main text, Öm as the LogMeanExp operation has the following property:
exp(ôm (flog(p)) log(p.)))) < max(pi..... .Pn)
On the other hand, max is in general the lower bound of the probability of the union of random variables.
max(log(PI).. Pm) < IP (x1 V xm).
19

---

Thus, taken together, we have
exp(Om ((log(p))...... .... logi))........mm))m))m maxxpp.....,Pii S Pxx V ... VXm) log (x1 V.... Xm).
On (f)og((pi))........ (ogpp)))...,,,(pmm))) log(pm)3)
Lemma 3.7 (Reachability Percolation Lower Bound), Given a weighted adjacency matrix W [0, 1jaxd for any 0 S < d, and for any pair of nodes y E [d], we have
(w, (2, S log Be-Be)(W) (c..u) and (.0) U (2, S log IE Bem()) [ (c,00)
(00) Proof. We first show that Rw(2; (z. y) s log A W)
R (z. v)]
Intuitively, this inequality holds because Bernoulli random variable parameterized by Wa.y. the paths as random variables are either mutually independent, when they consist of distinct edges, or mutually positively correlated, when they share some common edges. Thus, using the results of Lemma A. the t-norm of the path probabilities yields lower bound to their joint probability, and the t-conorm yields lower bound to their union R()) probability. Thus, the continuous reachability score ww consisting of disjunctions and conjunctions of the path and edge probabilities, also gives lower bounds.
probabilistic graph W, where each edge x + y is
More specifically, we can observe this fact via induction. In the base case, we obviously have
NW (x, v) log(1 [ (xr y)) log P A>~Bem(W ( (s,, u)) Pw y) log W Wa.y.
In addition, we can obviously that the random variable R2 (2r, u) is independent of the random variable An,y for any nodes u. E [d].
Then, assuming the inequality holds for graph with maximum path length of - 1, and that the random variable f((-1) (x, u) either independent of or positively correlated with the random variable A Au,y for any nodes u E [d], we now aim to show that these two statements also hold for I.
First, we see that, for any node E [d].
(11-1) (x, u), (ogg(WWW)))) log Aw Th (
(RQ-") (ar, u) A Au.u]
from results in Lemma A.I regarding the t-norm T3 and the induction hypothesis that R (1-11) cither independent or positively correlated with A
(z, u) is
Then, again using the result in Lemma A. regarding the t-conorm Ömz we have
RW (z. v) = Öott (F2
((1-") 11 (Ix u). log uy)))eejdd U (RW ( (z. (c.1)3)
V
(RM-") (z, A Aun) v (2.-1((2))) (2.1)
< log E A.-
wEid
R (I. (1.8))
log E AemBe
Regarding the relationship between R( (x, u) and An,,, we can see that the probability of R2 (ar, u) increases if the probability of An.,, also increases. This is because the formula, there is no negation
and A makes possible contribution to the value of R()) (xr,u.. u). Thus, Rlll (x, u) and Au, positively correlated.
Now, we show for the unreachability inequality, UW (a. v) I log E A-BB
(U(((, v))
Similarly, we show by induction. In the base case, we have
ū- (0) w (x, y) = (og(1(1 (x * v)) log I (e (z, v)) = log(1 Wa,y)
20

---

In addition, we can see that the random variable 1 (x, u) is independent of the random variable -R A (x,y). y)-
Then, assuming the inequality holds for graph with maximum path length of - 1, and that the random variable U=(""(x,, (ax, u) either independent of positively correlated with the random variable A - (x, v). We aim to show that these two properties still hold for I.
First, using the result in Lemma A.I regarding the t-conorm Öm, we have, for all nodes E [d].
W Óz(Uw " (z. log(1 Wwy)) Bew (W)
U (1--111 U
-R (u. v)
(ur, "I: independent or positively correlated. They are also independent of or positively correlated with U.A (=((-1) (xt y). This is similar to the mutual independence positive correlation among paths, since for these "non-paths" if they consist of distinct "non-edges", then they are independent. Otherwise, they positively correlated, since increasing the probability of the shared "non-edge" (or decreasing the probability of the shared edge) would simultaneously increase the probabilities said "non-pan "hon-paths." Thus, again using the result in Lemma A.1 regarding the t-norm Ta+1, we have
Then, we that the random variable (U7--1)(
V V (u, v)) for different are mutually
TW w (xx, y) = Tati (102(U (2, u), log(1 - Way))aaa U UlOW (x. (c.1)))
((A, (U"l"1(2, v -vvu)))) ( a ^UX (x,
I log E Am~Bem(W) (W)
(c.
= log A~ (W)
Theorem 3.8 (Lower Bound Expected d-Separation Statements). Given weighted adjacency matrix W E [0, 1j11d for any three nodes y, E [d].
(10) W (x, y) S log A A-Bem W)
S (z. vi
5(1) ww (z. S log E A~-Beei((W
(1)) (z, 2)
C (x, < log Ae-Be(WW)
(0) (Ir, (Ir,v)
(z, CW < log BA~Bem(W)
(1) (x, =)
2(1) Proof. First, we reiterate that, in the computation of differentiable d-separation 5(0) and W we use the continuous unreachability score ü:dd) directly, rather than taking the negation of the continuous reachability score.
Using similar reasoning as in the proof of the previous lemma, we first observe that, when treating S0 s(1) ((O) and C( as random variables when A ~ Ber(W ), all terms in the d-separation/d- connection formulæ (Definition 3.1) are either mutually independent or positively correlated. In particular, this holds true for S (I, a) and AR 2) for any a e [d]. as well as for C:(0) (y. 6) and RA (b, 2) for any b E [d]. Thus, using the results from Lemma A.1, we have the lower bounds as stated in the theorem.
Lemma 4.2 (Consistency of Multi-Task CI Losses). Given a faithful data as the sample size + co and the LogMeanExp temperature a - O+ (Equation (8)), the optimal DAG A achieves the minimum value on each of Lrr-o, CTP-1, Lixo, LiN-I, and LDAG-
Proof. Let 0* E Rdxd be the parameter that approximates the optimal binary DAG A of the given data D via o(0)) + A
When the sample size n - oc in a faithful data D, we have PD v) + 0 and PD (z. 2) + 0 when 1L, 3 and x 4 A- y 2, and conversely PD (z, y) + and po(2, y 2) + when IL A- and x IL A 2 for any v.Y, E [d], where A is the optimal causal graph of D. When the LogM temperature
0+, we have log max J (((lgg(())))))) + 0+.
21

---

In addition, we note that when given sequence of binary random variables Xi...., Xm, whose Bernuuli probabilities approach in the limit to either or the product t-norm Tm and max t- conorm Om over x; s prohabilities approach in the limit to the probability of [P(x1 A... A Xm) and P(x; V... V Xm) respectively. This is because the limit case, the random variables reduce to constant 0's and constant I's, in which case the product t-norm and t-conorm reduce to the vanilla Boolean operations of conjunction and disjunction. Thus, combining the fact that LogLTN's [2] t-norm T. is equivalent to Tm in the logarithm, and its t-conorm On approaches in limit to On in the logarithm, we have the following:
log PP (x1 A A... A xm) T' (flog F(x+) - 0+
log IP (x1 V V... V xm) O. (flog P(x;) | + -+
Using this property, we can see that for the continuous reachability and unreachability, for any given maximum path length I and any nodes x, y [d]. they satisfy
R (z, v) R2!) (000- (z, y) + 0+
U T (z, v) 1 a(o-)
0(0-)(221))) (I.
- 0+
Consequently, for the differentiable d-separation and d-connection scores, for any nodes x, V, E [d], they satisfy
100 (2x, v)
5:00 ((00) (z. y) + 0+
(z. (x,y
50) 0000 (c. 2) + +0
C 10 (x, (x,,)) - Col0 (x.y) - -0 C (x,yl - ((00)(( (I, y 2) - -
Furthermore, recall from Theorem 3.2, since A is a discrete DAG, ( s c and c reduces exactly to the ground-truth d-separation and d-connection statement values in A Thus, combining with the result above, we have that
((0) (20,) v)
{1 when 11. A y
5(1) (00)) ((0-)( (I, 22)
{1 when x 11 A y
(0 otherwise
0 otherwise
(:00) (000)(111)1 +
1 when x +1 A- otherwise
( (x,1 (0,0)(0 +
1 when 14 A- y
otherwise
Therefore, combining with the values of the p-values pp, have the following value matching between the model's predicted d-separation and d-connection score and the p- p-values. For all (x. v) e Lo (Iu [(x, v) x, y [d]. x >> v/) and for all (æ, y. 2) E I1 (I1 (1, v. x, U. 2 E [d]. x
2 # z # z}).
$(0) 000-,(2 (z.v) v) + 00 when (x) ++
I when (2, 3) +
50) (1.0
1 2)
[1 when pD (z, (r
2)
10 when pp [z, 2) +0
(40))
(2,-)(2,1)) (2,8) v) + o when Mo - p2( +0 0
1 when Mo -pD(,,)) +
0(0-)(2 (z,y1 v +
(o when Mi pD(2, y | +
when Mj pD(E,y =) +1
where Mo- maxi, PD (2. 3) 1, and similarly Mi maxı, pD(a, =) + 1. In other words, the differentiable deseparationddccomeciioo scores over the optimal parameter 0*) and the p-values from the data have zero mismatch.
Thus, the TP and TN losses have the following values:
Cir-o(0, D) - Σ
1/
11 A-
vl.
(Ir-10,,
--
Σ
1[x IL A y 2].
(a.g,z)el,
LiNo(H, (). +
Σ
14 vl.
(0,
Σ
1/x A-
2)
(w.wllo
(z.H.))TT
which are the lowest possible values that these loss functions can achieve, because any other configu- ration of the values for the differentiable dessparation/d-coomees scores would result larger loss value due to possible mismatches with the p-vauuss. Thus, the TP and TN losses are consistent.
22

---

Finally, since A' is DAG, the log-de det acyclicity loss DAG also achieves the minimum value, for any value of the hyperparameter s. This property is proved in Bello et al. [4]. Thus, all loss functions in Definition 4.1 are consistent.
I
Lemma 4.3 (Mutual Independence of oower CI Statements). All Oth- and Ist-order CI statements over the variable sets Iu = {(x,)) z. y E [d] > yt and I = {(x, y.2) x, x, 2 E [d], > y, # 2. # 2) mutually independent. In other words, none of these CI statements can be implied from any other of these CI statements via the graphoid axioms [24].
from
Proof. Let X' . S [d] be any non-empty node sets, and 2, W S [d] be node sets that could be empty. The graphoid axiom [24] states the following five rules describing the relationship and dependencies between different CI statements in a graphoid dependency model:
1. Symmetry: X L L 2 = y' 1L X
2. Decomposition: X' 11 yUW 2 (X 11 2) II W 2)
3. Weak Union: X 11 yUW| 2 (X 11 y ZUW) (X TL W ZU))
4. Contraction: (X' 11. y 2) A (X IL W 2U)) X 11 YUw| 2
5. Intersection: (x 11 y 2 U A (x n W ZU)) X II yUw! 2
We check rule by rule whether any of the low-order CI statements considered in ll and Il can be implied by any others.
1. Symmetry: No CI statement in Il and 1 can be inferred from others via the symmetry rule. This because Tg and Tn consists of only asymmetric statements where always x > y.
Decomposition: In the non-trive non-trivial case when | 2 and [W] 2 1, this rule does not apply because its left-hand side (LHS) statement involves nodes set Y W with cardinality of at least 2. This type of statement is not considered in Tn or Th.
Weak Union: Similarly, in the non-trivial case when |}]| 2 and [W] 2 1, this rule does not apply because its LHS statement is not included in In or It.
4. In the non-trivial case when | 2 and [W] 2 1, the LHS applies, but the right-hand-side (RHS) involve node set y U W with cardinality of at least 2, which is not included in mo or Il
Intersection: Similarly, in the non-trivial case when [}] and [W] 2 1, the RHS involve node set yu U W with cardinality of at least 2, which is not included in To or TT.
Thus, no CI statements in Ly and ll can be inferred from other statements in Hg and 1₁ via the graphoid axioms. Thus, the CI statements we considered in our approach mutually independent from the perspective of the graphoid dependency structure, justifying the use of the sum aggregation over these CI statements in the TP and TN losses.
B Further Discussions on Percolation versus Diffusion
The distinction between graph percolation and graph diffusion is fundamental to understanding the probabilistic interpretation of our differentiable d-separation framework. Percolation theory studies graph properties (such as connectivity or reachability) when the graph structure itself is randomized [12, 16]. In contrast, diffusion theory typically assumes a fixed graph structure, with randomness arising from particles making probabilistic traversal decisions across this deterministic environment.
In the context of our weighted adjacency matrix W E Rdxd these perspectives offer different interpretations. From a percolation standpoint, each entry W,.y represents the probability that edge + exists in the graph. The reachability question then becomes: what is the probability that directed path exists from one node another, considering all possible graph configurations? From a diffusion perspective, W would instead represent transition probabilities on a complete graph, where
23

---

Random Graph W
True Reachability Probability (Percolation):
P(x ~ y) = 0,4 (0.2 0.7 0.6 0.5 -0.2 0.7 0.6 0.5)
0.2
0.7
0.1592
X
Z
y
Max-Product Bellman- Ford Lower Bound: R (x. y) 0.4 max 0.7, 0.5)
0.4
0.6
= 0.12
0.5
Random Walk (Diffusion) is Upper Bound:
W
p(RW) (x y) 0.4 (0.2 0.7 0.6 0.5)
= 0.176
Figure 5: Example random graph to illustrate that the Max-Product Bellman- FFor reachability computes lower bound of the true percolation- based reachability probability, whereas diffusion- based random walk computes an upper bound. Marked numbers are the edge probabilities. For the random walk computation, we assume a model with self-loop probabilities, e.g. P(x - x) 0.6 so that the transition probabilities to other nodes and the self-loop probabilities sum to one.
at node v, the weights (W W..hhuedd dj (possibly normalized) determine the probability distribution of particle's s next position. The key distinction is that percolation randomizes the environment (the graph structure), while diffusion randomizes the trajectory through a fixed environment.
It is therefore most appropriate to frame our problem as a graph percolation problem. Given weighted adjacency matrix W. we interpret as parametrizing a distribution Bern from which random discrete graphs are sampled, A ~ Bern(W). Our goal is to estimate the expected reachability (with path lengths up to d), (e-BBe(WW) structure itself.
R(d) (x. v) where the randomness comes from the graph
This insight guides our estimation approach for the percolation-base reachability. We observe that generalized Dellman-FFord-basee differentiable reachability score (Definition 3.5) provides a lower bound on this expectation, while graph diffusion methods like Random Walk algorithms typically yield upper bounds. This distinction becomes clear when considering two potential paths P1 and P2 from node æ to node that share some edges. In percolation theory, reachability requires only that at least one path forms completely. Our Bellman-Ford approach with operators captures this intuition-when one path's 's formation probability approaches 1, additional paths contribute minimally to the overall reachability probability. The lower bound property stems from the approximation gap between product-max t-norms/t-conorin and the true intersectioo/union of events, combined with the approximation in LogLTN's s [2] LogMeanExp operator (Equation (8)). Formal proof of this lower bound property can be found in the proof for Lemma 3.7 in Appendix A. Consequently, our method will systematically underestimate true percolation baaeed reachability.
Conversely, in diffusion models, particles make independent decisions at each node. When two paths overlap, the probability of particle reaching the target via either path is treated additively rather than as a union probability. This independence assumption fails to account for the correlation between overlapping paths in the percolation setting, causing diffusion methods overestimate the true percolation- based reachability.
Figure provides concrete example illustrating the distinction between max max-product Bellman- Ford reachability and graph diffusion scores. The example shows a weighted adjacency matrix W with 5 nodes, where edge weights are indicated in the figure (with zeros elsewhere). For the diffusion- based method, we choose the typical Random Walk algorithm, where we assume self-loops each node so that the of self-loop probability and outgoing edge probabilities equals (e.g.. node x has a self-loop probability of 0.6).
Consider two possible paths from to y. path Pi -
y and path P2 = w - y- P1 and pa are two random variables under the graph distribution Bern W). The true percolation-basee reachability probability, denoted by P( ~4 y). is computed the union of the
24

---

events of these two paths forming: P(P1 V P2). Since these paths share the edge + but diverge afterward, we must carefully apply probability theory on union of events. For the divergent segments. the probability is:
P((2 - u - 8) V (z + w B) = (= - v) + ((= w v) P((= - u n) w v))
Because these two sub-paths share no edges, they are independent, giving P((= + y) (z + w > y) P(z + v + y) IP(= + w - y)- Calculating the full expression yields y( 0.1992.
Our max-product Bellman-Ford approach approximates this union using the max operator: for the - y), P(= -> y)) Since the max operator always provides lower bound for the union of events, our approach yields the final value RW (z. y) = 0.12, which underestimates the true reachability probability. In contrast, the iittusion-haeedd Random Walk method, denoted by p(RW) (æ: ~e y). adds probabilities additively across possible paths: p(ww) (Iz ~* F(PI) + (P22). This leads to an overestimation of the true reachability probability, giving p(RW) (I ~ y) 0.176. This example clearly demonstrates why our max-product approach provides principled lower bound on percolation-base reachability, while diffusion- based methods typically yield upper bounds.
divergent sub-paths from to y, we compute nW (2. max (ax(P(=
Our max-product Bellman-Ford approach thus partieularly well-suited for our proposed model instantiation, DAGPA, as it provides consistent lower bounds throughout computation, ensuring a coherent probabílistic interpretation. While we focus on lower-bound estimation, an alternative approach could theoretically estimate upper bounds on expected d poparaationnd-coonettoon statements. For such method, diffusion-basee based approaches like Random Walk would be more appropriate for estimating graph reachability. However, this would require identifying a different pair of t-norm/t- conorm differentiable logical operators that yield upper bounds instead of lower bounds. Based on van Krieken et al. [38]'s comprehensive analysis of t-norm and t-conorm suitability for differentiable learning, we have not identified a configuration that simultaneously provides upper bounds and maintains gradient stability comparable to our max-product operators. This represents interesting direction for future research, potentially expanding the theoretical foundations of differentiable d-separation frameworks, and we encourage future work to explore this avenue.
C Full Algorithm and Implementation Details of DAGPA
DAGPA, our proposed instantiation of the differentiable d-separation framework, employs multiple techniques to address optimization challenges and result validation. This section provides comprehen- sive implementation details for all techniques used in our approach. We first discuss the multi-task gradient conflict resolution technique (Section C.I) and the Bayesian sampling method for enhanced weight space exploration (Section C.2). Since DAGPA uses Bayesian sampling, each optimization step generates new candidate causal graph. Consequently, we require a principled procedure to select the best DAGs from those sampled. using only information available from the input data. We present this DAG selection heuristic in Section C.3. Finally, we provide the complete algorithm pseudocode in Section C.4.
C.1 Multitask optimizer to address conflicting gradients
The five distinct loss functions- 0thh and 1st- order TP/TN losses plus the DAG acyclicity constraint (Definition r.nder our opti- mization a multi-objective problem. Recent research in multi-task learning has demon- strated that naively summing these loss func- tions leads to theoretically suboptimal con- vergence, as contradictory gradient direc- tions among tasks can impede optimization
Algorithm 1 PCGrad [43] Gradient Projection
Require: Model parameters 0, number of tasks K, loss functions {Cx) 1K
1: gu + VoLu (0) Vk E [K]
2: g + 9k Vk E [K]
for 1: K do
for j if PC
uniformly
[K] \ in random order do
5:
gj < 0 then
progress [20, 22, 31, 40, 43]. This challenge particularly pronounced in our framework, where the TP losses encouraging d-separation
6:
b Subtract projection of g. onto gj
7:
gt + g. -
g.
8: return update (((PC) = = 1i=1 g.e
Σk
25

---

naturally conflict with the TN losses that require d-connection. To address these inherent gradient conflicts, we adopt PCGrad [43]. a gradient projection technique that resolves conflicting directions by projecting one gradient onto the normal plane of the other when they conflict. We reproduce the PCGrad algorithm [43] adopted to our setting in Algorithm
C.2 Gradient-informed discrete Bayesian sampling of W
Rather than using conventional stochastic gradient descent (SGD), we adopt Discrete Langevin Proposal (DLP) [50]. mundient-informee discrete Bayesian sampling technique, to sample parameters 0 that yield low-loss configurations'. Our choice of DLP is motivated by several key factors. First, our empirical observations show that SGD often becomes trapped in local minima, possibly due to the highly non-convex loss landscape, in which case Bayesian sampling can naturally escape local continuous
minima and explores the weight space more effectively. Second, we favor discrete sampling because the most likely DAG A from weighted adjacency matrix W is obtained via thresholding 0.5. This threshold divides the parameter space into discrete regions, and continuous optimization or sampling approaches spend significant computational effort exploring within a single region, repeatedly yielding the same causal graph structure. In contrast, DLP's discrete support enables efficient exploration across different regions while leveraging gradient information from the continuous space guide proposals.
To adopt DLP for DAGPA, we first restrict the parameter space from the real space Rddd to a space of discrete supports, Daxx where IT> (D1,.... D,, D; Dj. Vi j} is finite set of values. Furthermore, we require D1 < 0 and Dm > 0, so that for any pair of nodes (I; y). if 0g. = D1, then W x..y (00..) 0.5, which allows us to interpret æ, y as more likely not having an edge. Similarly, if Oz.g = Dm. then Way = ol2..)) > 0.5, which allows us to interpret I, y as more likely to have an edge x + y.
In reality through empirical experiments, we found that D = { 2.0, 0.0, 2.03 works quite well. We hypothesize this could be due to 2.0 and 2.0 as sigmoid logits are not too extreme that they enable smooth and meaningful gradients, while also sufficiently far apart for the differentiable d-separation formulæ. Furthermore, the "middle point" 0.0 allows DAGPA model uncertain edges, whose final values can only to be decided after other edges are settled.
The next step to adopt DLP [50] is to de- fine the proposal sampling distribution and the acceptance-rejection step [15, 21J. This is the discrete Metropolis- adjusted Langevin algorithm (DMALA) vari- ant of DLP algorithm [50]. One of DLP's novel innovations is the parallel sampling of the parameters 0 when its proposal distribu- tion can be factorized along the dimensions, meaning at every step, multiple 0g for differ- ent i's may have their values updated. This technique significantly enhances DLP's effi- ciency compared to other discrete sampling methods. DAGPA fits this requirement as we can formulate factorizable proposal distribu-
Algorithm DLP (DMALA) [50] Sampling Step
Require: Model parameters 8, step size B8, gradient AO(PC), current energy (())
1: // Proposal Step
2: for 1 do
> Can be done in parallel 3: construct q () as in Equation (9)
4: sample 02 ~ qi(- 0)
5: // MH Acceptance Step
compute (
compute (PC) via PCGrad (Algorithm 1)
8: compute (('
(0' 0)
compute q(0 0') = -I, (00 0')
10: set 0 4- 0' with probability in Equation (10)
tion. Specifically, we adopt Equation (2) from Zhang et al. [50], but change the gradients those obtained via PCGrad [43] from Algo-
11: return new sample 0'
rithm 1. That is, let q0' 0) be the proposal distribution. can be factorized along the d dimensions,
q(0' ()) I (0)) 0). where (0%) () is a categorical distribution of the form:
Softmax (
(0: - 0,)
qa (0( 0) == Categorical
(9)
2
where AO(C) the projected gradients given by the PCGrad algorithm (Algorithm 1), and (3 is hyperparameter controlling the DLP step size.
'Note that this differs from the discrete DAG sampling A Bern W) in Section 3.2, which serves to theoretically establish the probabilistic interpretation of our continuous d-separation and d-connection values.
26

---

After a proposal 6' is sampled from @ according to Equation (9), we perform an MH acceptance- rejection step. Intuitively, this step is to ensure that we are not taking a step too far and landing into a "bad region" in the parameter space. Mathematically, this step is to ensure the Markov chain is reversible. Here, we adopt Equation (3) from Zhang et al. [50], but for the (negative) energy function, we simply adopt the heuristic of summing the multi-task losses (Definition 4.1). Thus, for DAGPA, we accepts the proposal 0' with probability
min (I. (1exp(U()) 9(0 07) )
(10)
8(0' 10)
where ((6) Lir-o(0, D) + (I,)(0, LINOO(0,)) D) + LI,-(0, D) + Losa(O,s) with D the dataset.
We reproduce the DLP sampling algorithm adopted to DAGPA in Algorithm 2. We note, however, that the combination of PCGrad[433 projected gradients (P(CC) and the energy function (()) as the sum of multi-task losses may not form a mathematically well-defined and reversible Markov chain for the sampling, because the gradients are not directly derived from (+) but have been post-hoc modified via gradient projections. Nevertheless, we argue both theoretically and from empirical observations that the gradient modification step with either PCGrad or some alternative multi-task learning method essential, as it addresses gradient conflicts and navigates the gradient landscape much more efficiently. Furthermore, we empirically found that the MH acceptance step using the energy function U (0) obtained via summing the multi-task losses is sufficiently performative and, more importantly, provides the acceptance rate as important indicator for adjusting the DLP step size B. which is one of the most important hyperparameters in DAGPA. We detail the hyperparameter choice in Appendix D.4. We invite future research to tackle this challenging of integrating multi-task gradients with MH acceptance step in more principled way.
C.3 Training-time DAG Selection
Finally, DAGPA employs a heuristic score to evaluate the quality of causal graphs obtained from sampled model parameters. This evaluation step is essential because DAGPA adopts DLP [50], and each DLP sampling iteration generates new model parameters and corresponding causal graphs, making necessary to identify the highest-quality DAGs for final output. Importantly, this score not a traditional "validation score" since it requires neither a separate validation dataset nor knowledge of the ground-truth causal structure. Instead, it operates solely the input dataset D that is already used during optimization and sampling. This approach offers a practical advantage over model- methods like NOTEARS [51] and DAGMA [4], which typically require splitting the dataset to create separate validation sets, thereby reducing the amount of data available for model training. Our score makes full use of all available data while providing a principled mechanism for DAG selection.
We name this score the "TPTN Ratio score", as it essentially checks the ratio of weighted true positive (TP) and true negative (TN) d-separation CI statements predicted by the model, while treating the p-values from the dataset as the soft ground-trum truth labels. Specifically, given the current sampled model parameter 0, obtain the weighted adjacency matrix via W (0), and then threshold to convert it to a binary graph Á, Å., = 1 (W..> > 0.5), Vx, y E [d]. Additionally, we would like to ensure that at every step of sampling, we are evaluating a DAG, as the graph À converted from the model parameters 0 found by the DI.P sampling algorithm may not fully satisfy the acyclicity constraint. This requirement can be achieved by pruning Á form an acyelic A. In DAGPA, we proceed by finding feedback are set (ARC) from Á and remove all the edges in construct A. We use the ARC method provided in the igraph package [7].
Then, setting the MeanEXx temperature (Equation (8)) very small value (e.g. a 1e 5), we use the differentiable d-separation scores (Definition 3.6) on A to obtain 5(0) and 5(1) Since a is small but not exactly 0, ()) and exp(5(()) close to but not exactly 0's or I's. Thus, threshold again to obtain the binary d-separation statements, 5 (a, v) = [expp (5 (10)
A [z, v) > 0.5] and ( (a,y y 2) = 1 (exp(S ()) (xc, y | =) > 0.5). We treat these as the binary d-separation statements predicted by the model.
27

---

Then, given the p-values pD from the data, we compute the true positives (TP), true negatives (TN), false positives (FP), and false negative (FN) scores as;
D)
Σ
S (I,yl (I,
Σ
(1,v v =) v
E.yE(d]
x.8.5E|1
TN(A,)) D)
Σ
(1
S (I,yl 2))(1 pD(2,
(1 (x, 2))(1 pD(2, 2))
x.yEld]
n.y.sE[d]
FP(A,D) Σ
( [1 y =) pD v))
Σ
S (2, 2) pp(2, =))
s.
FN(A,D)
E (1
S (z,y y 3))pp(z..)
Σ
(1
S (z,y 2))
z.v.c[dd
And then the TPTN Ratio score is computed as
TP A, D) + TN( (A))
TPTN-Ratio A, D) TP(A, D) TN(A, D) + A, D) (A,
(11)
which ranges in [0. 1] and the higher the score, the better the causal graph in matching the low-order CI statements found in the data.
C.4 Full algorithm of DAGPA
Combining all techniques together, DAGPA's full algorithm is given in Algorithm 3.
Algorithm DAGPA
Require: Data D, initial parameter 0o, number of steps T, number of best DAGs K, step size
1: for /= = 0: T do
2: We f (0.)
3: for (z, y) E [d]2 do
Can be done in parallel
compute (())
W. and (( as in Definition 3.6
4:
5: for V. 2) E 3 do
Can be done in parallel
a() 6: compute W. and ((1) W. as in Definition 3.6
7:
compute Lrv.. LTP.1. LIN., L.N.1, LDAG as in Definition 4.1
8:
U(O+) < Lrp- + LTP-1 + LrN + LIN-1 LDAG
compute - via PCGrad [43] (Algorithm 1) compute 0.+1 via DLP [50] (Algorithm 2)
9:
10:
11:
compute A++1 by converting 0++1 to discrete DAG (Appendix C.3)
12:
compute TPTN-Ratiol (A+1, D))
13:
return K DAGs from (Ar)T= with the Top-K K highest IPTN-Ration (atio((A,, D) score
C.5
Computation Complexity Analysis
The computational bottleneck of DAGPA is the p-values computation for all low-order CI statements and the reachability computation required by the differentiable d-separation formulae.
The complexity of p-values computation (which we GPU-accelerate) depends on the specific choice of statistical independence test. Take the Chi-squared test for example. Given data points, computing each single Chi-squared p-value takes O(n) as this is the time required for iterating over all data points and building the contingency. Thus, the overall time required is O(d')) as we obtain total of dl d number p-values for the unconditional and first-order conditional statements.
For the differentiable d-separation formulae: the reachability subroutine (Definition 3.3) involves maximum dl recursive steps, and each all-pairs Bellman-Food update step operates on ((d3) matrix entries (Equation (7)). Thus, the total computation time is O(d²). The resulting all-pairs reachability matrix will be cached and accessed by the cooneection/ddseprrattt score computation. Now, since in the lst-rrder d coomeecioon/d-sepaaatii scores (Equations (3) and (5)), we need reachability
28

---

the pode-deleted subgraph A-; for all conditioning node 2, this requires total of d + 1 all- reachability matrices, thus rendering the total runtime of this part ((d").
formulae (Equations (2) and (4)), computing (tD) ((, y) or ((O)) _A(x, y) for each (x, y) takes O(d) time, since it iterates over d possible common ancestors and accessing the reachability cache takes ((1). Thus, computing all-pairs Oth- order d-separation/d-conee condddennonoioonnnn scores take ((d)) time. Now, since later in the Ith-order formulae, need all-pairs Oth-order d peppaaatiooddcoonneettt scores for the node-deleted subgraphs A- for all conditioning node z, this adds another dimension. Thus, the total time is (dd)).
For the Oth-order deseparation/d-comeeceti
For the st-order dsseparation/dccoomcciiii formulae (Fquations (3) and (5)), similarly, computing the result for each query triple (x. y 2) takes ((d) time. Thus, getting the results for all triple of nodes takes total of (dd)) computations.
Thus in summary, computing all-pairs p-values takes (dd")) computations, computing reachability matrix (and cache for later use) takes O(d)) computations, and computing Oth-order and stt-rrder cewparaationdccoonettttt scores each take (dd*) computations. The total will be dominated by the reachability computation, which is ((d)).
Acceleration in Practice The total (ad) time complexity notwithstanding, we seize many promis- ing opportunities for acceleration in the code implementation of DAGPA.
First, we note that all computations of reachability and deeparatioo/ddcoonnctiio scores are matrix- vector operations, which can benefit from GPU vectorization. In practice, we leverage PyTorch GPU tensor library for all such computations, avoiding any explicit for-loops and significantly improving the speed of optimization.
Second, in practice can limit the reachability computation to only consider paths of a constant maximum length k; reducing the time complexity from (dd) to (d')),, rendering the entire pipeline Odd). In the experiments for larger graphs with d = 50 number of nodes (Appendix F.2), we limit the path length to k = 25. Observations during earlier method development suggest negligible performance degradation with significantly improved computational efficiency.
Third, the all-pairs p-values computation can be massively parallelized and then cached and reused. This is particularly useful in the case where one wishes to obtain multiple output causal graphs for the same input dataset to enhance solution diversity and exploration of the Markov equivalence class. In this case, the computation time for the all-pairs p- p-values step be amortized it can be reused in subsequent runs.
Finally, we expect future work to further examine how to incorporate sparsity assumptions into the graph further reduce the time complexity. For example, instead of assuming each node can connect to all d - 1 other nodes, one reasonable assumption is to restrict to a maximum of constant k degree. In that case, the reachability computation subroutine computation can drastically speed up.
D Experiment Details
D.1 Code and dataset release
We release our code and data in
https://githie thub. com//PurduuMIIIS/DDARR
D.2 Dataset creation and conversion
We generated synthetic binary datasets following the approach introduced in the k-PC codebase [18]. We simulated DAGs using both Rédōs-Rényi (ER) and Scale-Free (SF) graph models with d E {10, 50} nodes and expected edge-to-node ratios r E 4). For ER graphs, we set edge probability p r/d to achieve the desired are ratio. For SF graphs, we used the Barabási-Albert preferential attachment model with attachment parameter m = [r/2] to generate undirected graphs, then randomly oriented edges form DAGs.
For each generated DAG adjacency matrix A, we constructed a binary Bayesian network using pyAgrum [10]. Each variable X, takes values in {0, 1}. and we randomly generated conditional
29

---

probability tables (CPTs) for all nodes. For root nodes (variables with no parents), we sampled marginal probabilities P(X; = ~ Uniform (0.2, 0.8) avoid extreme probabilities. For non-root nodes with parent set (A))), we sampled conditional probabilities P(X; = PA(i) = c) ~ Uniform(0.2 0.8) independently for each parent configuration с E {0. | From these fully specified Bayesian networks, we drew n E {100, 1000, 10000, 100000} independent samples via ancestral sampling (forward sampling from the topological order). Since the data is binary, no normalization or scaling preprocessing was needed.
We also generated continuous data with simulated ER and SF DAG following the approach introduced in [4]. Specifically, for each DAG adjacency matrix A, we generated structural equation models of the form X, = ()EPA(A) Wyrf (X₃) + et. where (i( denotes the parent set of node i, edge weights w
sampled uniformly from [-2.0, -..5 U [0.5, 2.0]. and noise terms 6; ~ (O, 1) are independent Gaussian. For linear continuous data, we set (z) = 2. For nonlinear continuous data, we randomly
assigned each ʃj to one of four nonlinear functions: (x) E {a² _ .in(()) cos(x)) with equal probability. We generated datasets with the same sample sizes {100, 1000, 10000, }00000} and graph structures (ER and SF with d E {10, 50} nodes and are ratios r E {2, 4}) as the binary setting to enable direct comparison across data types.
For experiments on real-world dataset, we benchmarked our method and baselines over the Sachs dataset [28] and the LUCAS (LUng CAncer Simple set) dataset [14]. The Sachs dataset is a widely recognized benchmark in causal discovery research, consisting of protein signaling pathway data collected through flow cytometry experiments. It contains measurements of 11 phosphorylated proteins and phospholipids derived from thousands of individual primary immune system cells, gathered under various experimental conditions. What makes this dataset particularly valuable for benchmarking causal learning methods is that the pathways between these proteins are well- established in scientific literature, providing reliable ground truth causal graph against which algorithms can be evaluated. The dataset has been extensively used in numerous studies to assess the performance of causal discovery algorithms, including recent work that demonstrates how less restrictive modeling approaches can capture complex causal relationships in this data that traditional methods assuming additive noise often fail to identify. We obtained and did benchmark the subset of Sachs data containing approximately 800 samples with no perturbation. We used the standard discretized version preprocessed using the Hartemink discretization method, which converts the original continuous protein concentration measurements into 3-1evel categorical variables (low, average, high) while preserving the underlying dependence structure. The LUCAS dataset is an artificially generated benchmark specifically designed to evaluate causal discovery algorithms under different conditions. It features binary variables in causal Bayesian network structure and comes in several variants that present increasing levels of challenge. Among all experiments introduced in LUCAS, we used LUCASO (baseline unmanipulated data) of size 2000.
D.3 Evaluation metrics
The primary metric we adopt in this work is the Conditional Independence Matthews Correlation Coefficient (CI-MCC), which we describe in detail in Section 5 with visual illustration in Figure 2.
Additionally, we evaluate our method and baselines using the following standard graph structure-based metric:
CPDAG Arrowhead Fl: We compare predicted causal graphs with the ground-truth CPDAG derived from the ground-truth DAG. Since different methods produce varying graph types, we standardize all outputs by converting them to CPDAGs. For methods like NOTEARS [51] and DAGMA [4] that return DAGs, we convert their outputs to CPDAGs using the utility method from the causal- learn package [53]. For k-PC [18], which returns k- k-essential graphs containing cirele-marked edges (in addition to directed and undirected edges), we treat undirected edges in the CPDAG conversion. This treatment is appropriate because circle-markes edges in k-essential graphs denote edges whose directions cannot be determined from low-order CI statements alone. Once both predicted and ground-truth graphs are converted to CPDAGs, we compute the arrowhead F1 score exclusively over the directed edges in the CPDAGs.
edges as
CPDAG Skeleton Fl: Following the same standardization process as CPDAG Arrowhead F1, we convert all predicted causal graphs to CPDAGs. We then further convert all directed edges in both
30

---

predicted and ground-trutt CPDAGs to undirected edges, creating graph skeletons. The Fl score is computed on these undirected skeleton graphs.
CPDAG Structural Hamming Distance (SHD): After standardizing predicted causal graphs by converting them to CPDAGs, we compute the structural Hamming distance (SHD) between predicted and ground-truth CPDAGs. Each unit of SHD corresponds to a single edge difference between the two CPDAGs, including: missing edges, extra edges, incorrect edge orientations (directed vs. undirected),, or incorrect edge directions.
DAG F1: We also evaluate methods that directly output DAGs by comparing them with the ground-truth DAG. Since constraint-based methods like GES [6], PC [33], and k-PC [18] return CPDAGs or k-essential graphs rather than DAGs, this metric only applies to score-based methods: approach, NOTEARS (linear and nonlinear) [51, 52], and DAGMA [4].
We note that these graph structure- atructure-based metrics may pose an inherent disadvantage to our method compared to CI-MCC. Through empirical observation, have identified cases where DAGPA produces DAGs with similarly high CI-MCC scores but vastly different performance on structure- based metrics. For instance, one sampled DAG may achieve both high CI-MCC and high structure- based scores, while another DAG with comparable CI-MCC performance may score poorly on structure-basse based metrics. This suggests that while DAGPA consistently samples DAGs that align well with conditional independence patterns (as measured by CI-MCC), this alignment does not necessarily translate high performance on traditional graph structure metrics. We provide detailed analysis of this phenomenon in Appendix F.I.
D.4 DAGPA model details and hyperparameters
The most important and sensitive hyperparameter in DAGPA is the DLP sampling step size (Equation (9)). To this end, we first find values for all other hyperparameters through preliminary experimentations then fix them, and only vary the step size for the experiments on the synthetic binary dataset and the real-world datasets.
Some of the other important hyperparameters and their values:
DLP support logit set D:: This hyperparameter controls the support logits that the model parameter θ can take during sampling. We use D [-2.0, 0.0, 2.0].
LogMeenExx temperature This hyperparameter controls the approximation accuracy of the t-conorm operator (Equation (8)) and thus the accuracy of differentiable d-separation scores. Setting this value too large will lose approximation accuracy, while setting too low will induce unstable and unsmooth gradients. Thus, during training or parameter update cycles, we use a Orrain = 0.01, while during evaluation when computing the DAG selection score, use a Oreval le 5.
DAGMA's [4] log-det acyclicity constraint hyperparameter s: This hyperparameter controls the valid region of M-matrices in which the log-det acyclicity loss is well-defined. Setting this value too large will risk model parameters stepping out of this region and causing undefined gradients, whereas setting this value too large will cause the gradient to have very small norm. In our experiments, for small graphs (include n = 10 synthetic binary dataset and both Sachs [28] and Lueas [14]) we use s 3.0, while for large graphs (n 50 synthetic binary dataset) we use s 8.0.
Finally, for the DLP step size 3, for each dataset we choose a different range to run hyperpa- rameter search and choose the best sampled DAGs therein according to the DAG selection score (Appendix C.3). The specific value range of 3 is chosen according to the acceptance rate in the DLP's Metropolis-Hastings acceptance-rejection step. Specifically, we choose the lowest value to be the one that can roughly achieve 0.8 acceptance rate, and the highest 3 value to be the one that can roughly achieve 0.2 acceptance rate.
n 10 graphs in synthetic binary dataset: ,3 E (0.5, 0.7, 0.9, 1.0, 1.1,
n 50 graphs in synthetic binary dataset: 3 E (0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38). Sachs and Lucas: 3 e {0.76, 0.78, 0.80, 0.82, 0.84, 0.86, 0.88, 0.90)
31

---

D.5 Baseline model details and hyperparameters
We tested our method against 5 methods from the table in the following section. For PC and kPC, used chi-square independence test and chose significance level threshold 0.05, and tested with both first and second order conditional independence tests for kPC. We used the causal- learn implementation(53] for PC and GES algorithms. For GES we used the local BIC score. For NOTEARS and DAGMA linear mode, we used the 12 loss for non-linear mode, and we followed the default dimensions for the MLP layer as used in the authors' original code. Rest of parameters all followed default settings.
D.6 Compute resources used
For the baselines, we ran PC, kPC, GES and linear version of NOTEARS, DAGMA on a 64-core AMD Epyc 7662 "Rome" processor with 16 CPU cores and 32 GB memory requested. The non-linear version of NOTEARS and DAGMA were run on one A30 with same CPU and memory requirement. Every experiment is completed in 4 hours.
For DAGPA, run all experiments on an AMD GPU cluster, equipped with 32GB MI108 and 64GB MI210 and EPYC 7V13 cpu with 64 cores. Experiments on small graphs (n 10 synthetic binary, Sachs, and Lucas) terminate within 1.5 hours, whereas experiments on large graphs 50) terminates within 12 hours.
E
Related Work
The methodological evolution of causal discovery has progressed through two primary paradigms--constrainte constraint- based and score-based approaches- recent advances bridging their complementary strengths.
Constraint-Based Methods Pioneered by [34], constraint-based methods leverage conditional independence (CI) tests to reconstruct causal graphs. The PC algorithm [33] established core principles: (1) infer conditional independencies via statistical tests, (2) eliminate edges violating d-separation rules, and (3) orient edges using collider detection. While effective in principle, finite- sample reliability suffered from error propagation in high-order CI tests [9]. Extensions like FCI [46] addressed latent confounding through more sophisticated separation criteria, and Kernel CI Test (KCIT) [48] enabled nonparametric testing via Hilbert space embeddings. Modern variants like LOCI [42] demonstrated that low-order CI statements suffice for structure recovery under specific faithfulness conditions, while k-PC [18] provided theoretical guarantees for bounded-order testing.
Score-Based Methods Score- based methods reformulated causal discovery as a combinatorial optimization problem, maximizing score functions (e.g., BIC [29]) over DAG spaces. These methods can be broadly divided into those that search the discrete graph space and those that leverage continuous optimization. Among discrete methods, GES [6] advanced this paradigm through greedy cquivalence class search, Another major family consists of Bayesian methods, such as BayesDAG [1], which define prior distribution over DAG structures and likelihood function. By applying Bayes' rule, they compute posterior distribution over DAGs, typically using sampling methods like MCMC explore the structure space and average over models. While powerful for uncertainty quantification, these methods, like greedy search, are computationally intensive and constrained by the discrete search space. The field transformed with [51]'s NOTEARS framework, which redefined acyelicity through a differentiable constraint:
tr(e WoW -d = 0,
where W the weighted adjacency matrix. This enabled continuous gradient-based optimization. spawning derivatives like DAGMA [4] with improved learning stability and enhanced gradient behavior via a log-determinant characterization,
log (W) + d dloo
0,
and GOLEM [23] using likelihood-hased objectives. Nonlinear extensions [52] incorporated neural networks to model complex functional relationships while preserving acyclicity.
32

---

F Full Experiment Results
E1
Analysis of CI-MCC versus graph structure metries
In this work we primarily showease the conditional independence Matthews Correlation Coefficient (CI-MCC) metric, along with auxiliary graph- seaph-sstucture-bassee metrics like CPDAG Arrowhead F1, CPDAG Skeleton F1, CPDAG SHD, and DAG FI. We notice, however, that the graph-structure-bass metrics may pose unfair challenging to DAGPA. In particular, we found that the graph-structure metrics may not always align with DAGPA'S 's objective of matching the causal DAG's predicted d-separation statements with the low-order CI statements found in the dataset. There are many cases where the DAGs returned by DAGPA made few mistakes in aligning the CI statements, yet are still scored badly by the graph-structtur-biie based metrics. We provide concrete evidence in this section.
Here, we pick one experiment run of DAGPA on one of the syntheric binary dataset gener- ated by an ER graph with n 10 nodes with in total N 10000 samples. This particular run performed 1000 sampling steps, thus giv- ing us 1000 sampled DAGs. For each sampled DAG, we compute the DAG selection score (Ap- pendix C.3) and the CI-MCC metric score, as well as the scores of all graph-structure- based metrics. To recall, the DAG selection score is used by DAGPA to select the best DAGs out of all that were sampled during training and requires only the input data to compute. It is not test metric. Figure 8 shows the scatter plots of comparing CI-MCC to graph-structure- based metrics, Figure 7 shows the scatter plots of comparing the DAG selection score to graph- cructure-based method, and finally Figure 6 shows the relationship between the DAG selection score with CI-MCC.
0.9
0.8
800
Sampling Step
0.7
CI-MOC
600
0.6
0.5
100
0.4
0.3
200
0.2
0.2
0.3
0.4
0.5
0.6
0.7
DAG Selection Score
Figure 6: DAG selection score versus CI-MCC. There strong positive correlation.
From Figure 8, we can observe that, although the graph-structure- based metrics have generally positive correlation with CI-MCC, there are many samples concentrated in the bottom-right the region where these sampled DAGs achieves high CI-MCC but low graph-structure based metric values. Moreover, especially on regions of high CI-MCC values, the samples vary a lot on their graph-strructurebaasse metric values. This implies that there is misalignment between CI-MCC and the graph-structtre-based metric, where performing well on matching model's d-separation with data's low-order CI statements do not translate similarity on the graph structure.
Figure reveals a similar story. In this case, the x-axis is the score actually used by DAGPA during training. The similar pattern shows that, even if DAGPA discovers high-quality DAG with high DAG selection score, which in turn translates to making minimal false positive and false negative mistakes the low-order CI statements (Appendix C.3), it may still have drastically different structure than the ground-trutt DAG or CPDAG, yielding a very low graph-structure based metric value such as a low CPDAG Arrowhead Fl value. We leave the problem of proposing alternative DAG selection score that may yield much more positive correlation with graph-structure-basee metric to future research.
As samity check, Figure 6 shows the relationship between DAGPA's 's DAG selection score and the CI-MCC metric. Here, can finally observe strong correlation, demonstrating the validity of our proposed DAG selection score it is correctly doing what is designed to do, i.c., checking alignments of low-order CI statements between the model and input data.
Finally, throughout all three figures, we visualize the sampling step numbers in color scale, where the darker colors corresponds to early stages in the sampling, and lighter colors corresponds to later stages. We can observe that, in general, there exists a gradual transition from darker color to lighter color when going from bottom-left to top-right. This pattern suggests that DAGPA is indeed gradually approaching the regions in the weight space that have better and better performance, across all types of metrics.
33

---

CIOND SRAUSTITT
(a) CPDAG Arrowhead F1
(h) CPDAG Skeleton Fi
(c) CPDAG SHD
(d) DAG FI
Figure 7: Scatter plots of CI-MCC metric versus standard graph-structurr-baat based metries of DAGs sampled by exemplar DAGPA run a 10 synthetic binary data. The graph-structure- based metries are not strongly positively correlated with CI-MCC. Many DAGs have similarly good CI-MCC values, but vary significantly in their graph-structure-bant metric values. Thus, the graph- structure-based metries pose an unfair challenge to DAGPA that focuses on aligning low-order CI statements.
DAC Scadics Scar
57 AAD Stintee Sact
DAG Sww.ciee...ooo Sooe
DAG Delectio Score
(a) CPDAG Arrowhead F1
(b) CPDAG Skeleton Fl
(c) CPDAG SHD
(d) DAG FI
Figure 8: Scatter plots of DAGPA's 's DAG selection score (Appendix C.3) versus standard graph- structure-ban based metrics DAGs sampled by an exemplar DAGPA run on 10 synthetic binary data. Similar to the case of CI-MCC, many DAGs with similarly good DAG selection score may have vastly different values on the graph-streture--bsse score.
F.2 Results on synthetic binary data
We present the full suite of results on the synthetic binary dataset in Tables I to 4.
F.3 Results on real-world data
We show the results with additional metries on Sachs [28] dataset and on an additional real-world dataset Lucas [14] in Figure 9.
G Future Work Discussions
G.1 Integration with Score-Based Causal Discovery Methods
Our differentiable d-separation framework provides novel optimization objective based on con- ditional independence constraints, which can potentially be integrated with existing score-based methods like NOTEARS and DAGMA structured regularization technique. This section briefly discusses potential strategies for integrating our approach with existing score- based methods, particu- larly NOTEARS [51] and DAGMA [4]. We focus on the key technical challenge: ensuring parameter compatibility between score- score-based based and Cl-constraint objectives.
The challenge is ensuring compatibility between the model parameters 0 used in score-based methods and the probabilistic edge weights W E [0, 1]ddd required by our framework. For example, for NOTEARS linear models, applying sigmoid activation a (0) naturally produces valid probability measures. However, for nonlinear variants where the weighted adjacency matrix is derived from neural network parameters (e.g., (8U))TT)) for the first MLP layer), additional normalization is needed map unnormalized negatiee values to [0. 1] while preserving meaningful probabilistic semantics. Once this parameterization is established, the score-based prediction loss (e.g., MSE reconstruction loss) can be added as an additional task our multi-task optimization framework, combining functional relationship learning with explicit CI constraint enforcement. While this
34

---

Table 1: Experimental results on Synthetic binary ER, r 2
(a) Synthetic binary ER, r =
Medhod
Dataset Settings
d=100, no0
d=10, n=lk
d=10,
n==00n==000 n=
d=10,
m= 700k
d=500, n00
d=50, n=/k
d=500, n=/0k
d=300 m= 100
DAGMA (Linear) [4]:
0,2319 Looes
0,2259 Juin
0.2465 oonc
0,2380 Lone
0.0925 comss
0.0867 Lome
0.0850 1000s
0.0904 100150
DAGMA (Noolinear) [4]:
0,5133 anies
0.3544 nnees
0.5225 -0300
0.4707 .0.433
0.2668 a3131
011335 4000
0.1725 400652
0.199) a00
NOTEARS (Linear) [51]:
0.2626 souer
0.2276 sume
0,2689 -0083
0.2586
1.1069 sams
01881 -wmm
0.0915 10007s
0,983
NOTEARS (Noulinear) (52)
0,1949
0.1810
0.2123
0.1901
0.0758
0.0727
0.0718
0.0734
[6]
0,2983
0.8546-
0.7783
0.0178
0.1844
[38]: APC (=R)
0,2845
0.86666
0.7922
0.001
0.1073
0.3135
0,7585
0.8923
0.055
(k-2) [18]:
0.2943 J.omes
0,7100
40:154
0,5662
00004
0.1009
1000000
0.1679
000990
0.1139 411440
DAGPA (Ours)
0.4672 15
0.6848
-00011
0.8223
0000
0.797 -008
0.0604 +000
0.2156 -00ss
0.0601 1000
0.3844 w00
(b) Synthetic binary ER, r = 2, CPDAG Fl Arrowhead
Method
Dataset Settings
d=10, s=1000
d+ 11
d=10, ntl
d=10, 1000
di-50) n= 200
df-50, w=/k
d-36, 2-/
(fx53)) A= 1004
DAGMA (Lincar) [4]:
0,4273 40000
0:0000
120021
0.0000+1.000
0.0085
0,0214
0.0143 110
0.0146 4000
DAGMA (Nonlinerr [4]
0.2143
0,2344 11083
0.1921
0.3670 =
0.1791
0,3194 -1972
0.2749
NOTEARS (Lincur)
0,0296
0,0000
0.0000
0.0696
0,0213
0,0239
0.0251
NOTEARS (Ncollnncar) [52).
0.3424
0:0000
0.0133
0.0000
0.0519
0.0102
[6]:
10:1994
210113
240002
18353
0.6911 -1993
1331:
0,1025 -03018
0,3607
402211
0.5673
40:14
0.7620 ..154
0.1242 L3365
0.4024
0,592N -1.000
VPC (8-1) [18]
0.0524-6111
0.1991
0:2435
10:100
0.3359
+3385
0.2255
0.3395
0.7777
+0015
APC (k=2) [18]
0,0912 ...000
0.2298
110177
0,5012
10.825
0.5329
00993388228
0.2788
2.802
0.4137
-83009
0.556]
100ess
DAGPA (Ours)
0.1052 .0276
0,1852 80006
0,2041
100912
0.206444000
0.0005 18304
0.0047
22820
0,0015 10o0
0.054]
10003
(c) Synthetic binary ER, r = 2, CPDAG F1 Skeleton
Method
Dataset Settings
d= 10 1=000
d= mm
d= 70k
d= 10, 1000
5=50 m= 700
#=500 1m
d=500 mz 10k
d0,50 1006
DAGMA (Linear)
0.2921
0.3034
0.2082.4.123
0.2417
0.3431
DAGMA (Nonlinctr)
0,5418
0.5362
0.6114
0.6341
0.3667
0.5542
0.5701
NOTEARS (Linear)
0,3535 _ea0
0.2252
0.357418299
0,3057
NOTEARS (Noollttar) [52]:
0.1415
+0.133
0.1991 .012%
011124
0.1347 -1942
0.1319 .1881
0.1047
GES [6]:
0,3845 -0136
0:2013
10080
0.8018
20084
0.271 +0200
0,4003 8550
-
0.8588
8627
PC 1391:
0,6202 -03254
0,8121
0,9314
L0042
0.711
0
82272
0.7401 2
0
-2204
0,9626
10021
KPC (K-1) (18)
014270 -000r
0.7831
0.7793
0.8044
-1355
0
-3300
0.7488 -180
0.7560
0.688)
+0022
XPC (R-2) [18)
114302 0.0254
8.8022
00000
0.8968
00032
+337 0
0.4044
83275
IL74P6 2.11
018622
09110
00032
DAGPA (Ours)
0,4062 .0080
0,4273 100813
0,6642
0.44221.900
0.0062
0,0826 2.85
0.0079_1800
0,1043
10003)
(d) Synthetic binary ER, r = 2, CPDAG SHD
Method
Daltaed Settings
HL, nc Z00
40.00 a:: 10B
d= 100
dc,00 12
di:30, ez 10M2
DAGMA (Linear) (4):
17,0000 1.97
27.5000) 11803
17,1000 +3.35
15, 15.1000=15500 1.539
27.5000011159
20.2000-22010
57,2000 11
20,100012.400
DAGMA (Nostinear) [4]:
17,4000 23.70
10. 40000 47.770
15050000 15351
10. 1490
15. 4000 24:210
51.1000-11.01
NOTEARS (Lanee) 1511
17,7000 11500
17,5000
17,3000 1180
15.1000
35,3000 12000
96,6000 43214
05:9000 .3.3/3
(Naslincar) (53)
13,5000
10.1000
17,5000
9,0000
4,9000
95,0000
70,1000
49,1000
31,0000
KPC (5-1) (18):
17,5000
19,5000
94,90000
77,4000
134. 150.800022000
APC (em2) [18]
17,4000
12,5000
12,0000
10, 9000
94,5000
60.2000
$5.5000
50.3009
DAGPA (Oum))
21.4600 +3 m0
23.6200 3130
25, +1.5 8m
24 14000
100,6000 11
190. 180.5400-+3-00
106 2000
220.7500
(e) Synthetic binary ER, r = 2, DAG F1
Method
Dataset Sertings
70 100
d= 10:10 mm12
d= n= 70B
d= 10, 1000
4=30 m= 700
d=50 1m
d=,00 am 10k
d= nm 1006
DAGMA (Linear)
0.1706
0.1225
0.2137
0.1677 441313
0.2020
011002
0.1766
0.179]
DAGMA (Nonlinctr)
0.3243
0.2734
0.3150
0.2243 -1452-
0.3087
0.3452
0.9820
NOTEARS (Linear)
0.2076
0.2301 10:4000
0.1905
0.2304
0.2027 _c390
0.209] 1000ss
NOTEARS (Noolintar) 152]
0.0059 ao00d
B:179
40007
1301664
0.0827
0.0665
0.0725
0,0507
DAGPA (Ours)
0.1557
0.2333
(.2176
c0871
0.26113
63117
0.0372
0.0091
-18443
0,0599
40055
35

---

Table 2: Experimental results on Synthetic binary ER, r =
(a) Synthetic binary ER, r = 4, CI-MC
Method
Dutaset Settings
d=10,x=100 x-108
d-10, n=ik
d-10 R-10k
d-16, 1=1000
d-50, x=100
d-50,g-1k
d-50, n-/Ok
d-50
n=100k
DAGMA (Lincar) [4]:
0,2398 40235:
0.2085
0.2120
2040
0.2225 aosss
0.0428 10005
0.0398
0.039
.1007
0,0055
DAEMA (Noulincts) (4]:
0.5724 -0:09
04450
+0.50
0.5007
11:33
0.4624 -1
0.2658 +00815
0.8991
D4773
-13151
0.003-9001
0.0623
-00030
NOTEARS (Lince) [51]:
0,2758 -ocoss
0.2197 400m
0.2344 1.565
0.2294 -000ts
01053656000 4n00a1
0.0419 4c00
0/3434
-1896
0.0423
000007
NOTEARS (Nonlinerr) [52)
0.1882 -m:m
0.1804 souess
0.2091 13
0.1894 souers
0.0369 stuse
0.0350 stoms
0.0342
:1
0.0061
-0:005
GES 16
0,3341 28:212
0.5932 A01352
0:9135
0.9135.1:007
(.9251 Lure
0.3863 10:331
0.1514
-242277
-154
0,0000
23,000
PC 1331:
0.3252 18352
0.3071 103734
0,5397 1.335
0.9317 10084
0.2514 Joanss
0.0181 10004
0,0049
-1915
0,0100
190155
EPC (1-1)
0.4053
0 en10s
0.9469 11875
0.9872 -0N0
4L1670
-00012 c051
0.0353 -18147
150000-0.000
KPC (k=2) [18]:
0.3507
0.5565
0,9469
0.9378
0.2243
0.0036
0,0054
0,0000
DAGPA (Ours)
6.6214 10:002
0.5299
0,9493
0.9484
0,331
0,3292
0.2669
0.7209
(b) Synthetic binary ER, 4, CPDAG Fl Arrowhead
Method
Settings
d-10, x=1007
d-100
d-10,
m-10
d-10 1=1000
J-50 m-700
5f-50 w-12
4-50 4-50, a-10
d-50
DAGMA inaar) [4]
-
0:0000-0100
0.0333
00000-1300
0.004
0.0051
0
0.0040
+0045
DAGMA (Nonlinear] [4]:
0,1844 .034
0.1668 Ju.n
0,0857
10.175
0.1991 141344
0.1193
0.0613
0.1158 18864
0.0976
NOTEARS (Linear) [51]:
0.3490 -0077
8:0000..c0 60000
0,0533
40.554
0,000 +0000
0.01901
0.0061 2000
0.0130
0,0020
NOTEARS (Nonlinear) [53]:
0.00000 -one
0.0000 -0100
0.0847
+0.58
00314-1800
0.0109 -0315
0.0010 -10
0.0101 -105
0.0070
10022
GES [6]:
0,0685 -0870
0,2752
40.7%
0.3504
102241
0.4993 111712
0.0684
0.2896 11040
0.5302 -1855
0.1253
10207
PC [33].
0. 1.4049 -0000s
0.30% -0.551
114223
-0.155
0.459 -013
0.0819-8306
012572
0.3632 -18:04
0.4153
LPC (K=1) [18):
0.1014
1965
0,2402
0.2566
0.0377
0,0064
0,0063
0.0200
APC IMM
0,1034
0,0532
0.0684
DAGPA (Oues)
031
0,2733
0.2589x0:29
0.0424
010572
0,0892
0.0645
(c) Synthetic binary ER, 4, CPDAG F1 Skeleton
Method
Datasel Settings
d=10, x=:100
d=10.
d=10, m=/
d=10, 1=0k
d=550 n=/00
d=330 m=1k
d=50, n=/ok
d=300 n= 105k
DAGMA (Lincar) [4]:
0.2706 -6 --d
0.2260 oooe
0.2055 c000
0.2292 9000
0L0982 +03341
0.0732
14291
0.0052 11555
0.075] .0009
DAGMA (Nonlincary [4]:
0.5460 c00ma
0,5399 strec
0,5137
20:32
0.5488 2:1224
0.2296
0.1997 :1:00
0.2653
1189
0.2288
11048)
NOTEARS (Lincar) [51]
0.3598 a0.3i
0.2533 4000s4
0,2439 40.127
0.2722 43364
0.1388
0,0915 11824
0.1193
0J960
400077
NOTEARS (Nonlinear) [52)
D. 0.14535-1
B.0848 -crex
0,1777 -0.15
0. -1492-1189)
0.0188
010408
000352.0000
GES [6]:
0,3903
0,6616
0,8163
0.7991
08:820
0.4322
0,7155
0,1574
[39]:
0.4194
0,8299
0,8080
0.5469
0,7493
0.7816
APC [18):
0.7724
0.6504
0425
(k-2) [18]
0.419
0,7178
0.8127
100030
0.7558 L.3002
02118888200
0,5530 13340
0.7416 -1910
0,7458
100035
DAGPA (Ours)
0,4217 .00002
0,4450 000SS
0,5309 4008
0.5036
0.1057
0,0637 -202
(d) Synthetic binary ER, r = 4. CPDAG SHD
Uniaaet Settings
d-10 n-/
d-10) 1t-i
2-100 A-10k
d-10 /w/
d--00
6-30, R-1
d-50 n-10k
df 8-10
DAGMA (Linear) (4)
21,1000
20,4000
22,0000
21,2000
197,2000
197,7000
196,9000
DAGMA (Noalinean [4]:
21,1000
18,3000
21,0000
19,0000
191,7000
183,4000.450
NOTEARS (Lanzar) [51]:
19,9000
20.4000
21. 60.60000-0
21. 21.1000+:300
195.20 160.1000-00
197.50 0075000-10
196.1000
197:0000+.220
NOTEARS (Nualincari [32]
21,9000
21,3000
22. 12.1000
21 62.0000 422/02
195 5000
159.1000
198
198.9000 -11/2
GIOS (6
21,2000
1000 1122
19 0000000 -03751
4000 17362
190. 2411
154.50 -1755000-170
130 20
25 0000 1157577
PC 1331
20,7000 23.005
17,5000 a 06/0
17,3000 4122
17.0000 +1305
203 0000 15
177.3 (77.5000 1060
164.4000 a (3)
171 1000 143913
APC (3-1) (188:
21,1000 223054
17,7000 1808
22. 5000 4444
22. 4000 13080
203 00000000000 44200
1 17712
211,6000 212001
252 7000 28
KIC (em2)) [18)
20,5000
600
00000
4000
30:0
179.7000
19% 9:00 14112
(Our))
26,6900:1
27,
14
-00
(c) Synthetic binary ER, 4, DAG FI
Method
Dataset Settings
d=10 1=
d=10.
d=10,
nz/
d=10 1=
d=30 n=
d= m=lk
d=50, nok
d=330
m= 100k
DAGMA (Lincad) [4]:
0.1661 -0:
0.1428
0.1365.a100
0.071 4337
0,0569
0.0618-100
0.0452 40010
DAGMA (Nonlincar) [4]:
0.2813 sermo
0,3242
20130:
0,2824
0.3234 291329
0.1280
:15
0.1739
0.144] 2000
NOTEARS (Linear) (51]
0.2563 200874
0.1508
40327e
0.1517
400820
0.1492
0.0905
000688 11000
0,0761
0,0512
100224
NOTEARS (Nevlinear) (53]
0,0772-6000
B.038 -omer
0,1082
0.0730
0.0279
0.0079
0.0195
0.0196
DAGPA (Ours)
0.2017 mec
8,2056-4665
0.2874 -00m
0.3427
0.0483
0.0688 -
0.0433
0.0681 4005s
36

---

Table 3: Experimental results on Synthetic binary SF, r = 2
(a) Synthetic binary SE, r CI-MC
Method
Dattset Settings
d-10, 100
d-10, 11
d-10 1010
d-19) am /00L
d-50,
nm 200
d-:0) nm/k
d-50, n0A
d-50 Am TOKA
DAGMA (Linear) (4)
0.2279 Lovers
0.2247 10342
0,2256
0.2167 Jounzs
0,1113 L
0.14K3
0.1475 1o0md
0.1560 LOUIT
DAGMA (Nenlinear) [4]:
0.4979 suass
0.3630 sume
0,4339
10125
0.4195
201210
0.3170 10.0
0.2208
0,2464
210010
0.2867 20uEs
NOTEARS (Linear) [5]]:
0.2716 a0077
0,2422 a00
0.2548
0.2556 4005
0.3433 030233
0.1551 2000
0.1614 40000
0.1659 10019
NOTEARS (Noulinear) [52]
0.2018 -ome:
0.1919 -omm
0.1940-0006
0.1847 -00770
0.126-aame
0.125-6007
0.1263 4008l
011341
GES [6]:
0,3268
0.6257
0.7965
0.8394
0.0051
0.1955
0.2143
0,1430
[38]:
0,3559
0.8635
-0.0206
0.0422
APC (1R]
0.8032
0.8687
(k-2) [18]:
0,3435 400120
0,7805
200152
0.8199
.03322
0.868] 400022
0.0261 4
0.1487
101257
0,0895 100930
DAGPA (Ours)
0.5491 00.401
0.7731 401272
0.8684 0.112
0.865 400024
0,0728 4.022
0.1055 400001
(b) Synthetic binary SF, 2, CPDAG Fl Arrowhead
Method
Dataset Settings
d=100 s=:100
d=10, m=1k
d=10, m= /0k
d=10 d=100,==000 1=00
d=50 no
di:50 m=
d=500 n=
d=50, m= 1066
DAGMA (Lincar) [4]:
0,0000 m.000
0.0125 4000
1,4000 000000
0.0087
(10000
0.0044
0.00S4 4007
DAGMA (Nonlincar) [4]:
0.1955 00136
0.2072 40.300
0.2012
40.235
0.1725
0.3465
[L1891 18870
0.2219 18851
0,2210
NOTEARS incer) [51]
0.0495
0.0211 -01e
0,0000 -0000
0.0000
0.0273+
0.0148 -1874
0.0132
0.0149 +000
NOTEARS (Nealinear) (52]:
0,0451
0,0000
0,0000
0.0190
0.0284
0,0210
0,0106
0.0167
[6]:
0,3086
0,5625
0,5590
0.4695
0.7163
0.7741
[39]:
0,1330= 20:194
20.17
04159
0,4583
KPC (8=1) [18]:
0.1301 0,000
0,1593 1831336
0.2905 100330
0.2411 49107
0.0157 4
0.1542
0.2068 -13310
02871 A
XPC (8--2] [18)
0.1131 -0000
111612-6100
0.3378
0.2496-11244
0.0729
83341
0.1670-1854
0.2671
-88554
0:4013
DAGPA (Ours)
0,1363 02042
0.1919 100012
0,2266
10.126
0.2374
0.0037 83305
0.0186 12813
0,0078 -1000
0.0387
100075
(c) Synthetic binary SF, r = 2, CPDAG FI Skeleton
Method
Dataset Settings
10 k=0
d= m=
d=10, mz 70k
d= 10, 1000
d=50 m= 700
d=500 wm
d=500 mz 10k
d=:30 mm 1006
DAGMA (Lincar) [4]:
0.2159
0,1933 11085
0.2072 s0.170
0.2296
0.2210
0.2187
0,2148
0.1912 100035
DAGMA (Nonlincur) [4]:
0,5808
0,4768
0,5342
0.5854
0.3076
015338
0.5412
0,5433
NOTEARS
0.2536-6190
0,2358
0.2663
0.3141
0.2857-4
0.2748
0.2745
NOTEARS (Noolinear) [52]:
0.0776
0.1063
0,0560 1344
[o]:
0.3601
.00040
118113 400427
0.819 aam
0.0875
0.8637 .1817)
03839 40022
PC [33]:
0,4474
0.7446 -00m
0.8168 -10010
0.8004 -8
00:42.5:991
0.7079 -1:04
0,8377
0.800 -000
APC (k=1) (18):
0,6635 Leace
0,7454
0,7774
1000M
0.7737 LA0M
0.3689 2
0.6872 2.834
0.7782
0.750] 10005s
KPC (R-2) (18)
1.4474-0000
0.7404
0.8120
00000
0.7894 40375
0.3612
017025
20374
018231
03336
4000
DAGPA (Ours)
01,0003 60000
0.4466 4007
0:4839 CUIEE
0.48-6
+8361
0,0889 --853
0.0357
0.0851
40005
(d) Synthetic binary SF, r = 2, CPDAG SHD
Method
Dhstoet Settings
des 10M
dc 20, 6:
Jo 10, am 10k
d= (0, H002
di:S0 nc (10)
J=50 12
d5c d5,50 x: HOME
DMGMA (Linerr) [4]:
1M 21190
19.20 19.2000204475 2755
19 220000115500
18. 18:20000:10 11.101
97.2000 1530
56, 56,503 1220
55.9 55.900022000 22375
2
DAGMA (Noninearr (4):
18 2000 20m5
17:00:0 1393
17. (0000) 1588
15.99 11.5000-2100 -21055
147.2000-11000 11:000
56 1000 -450
$5.000-1125 1125
35 35.20000xxm0 mm
NOTEARS (Lineus) (50)
400044
18,7000 333233
98,9000 11003
1980
94.7000 2045
05,9000 (1633
(Nenlinar)) (52)
18
4.0000xsiM
18:3000
14,2000
60000000000
31,7000
71.
KPC (5-1) 118)
180000
16,0000
1000
12000000
911
105
LPC (c=2) [18)
18.1000 10:04
16.0000
15 15,3000
17.5 17,5000
118,0000
53. 1000
75. 750000
74, 74,3000
DAGRA (Ours)
2447900.2445
25.7000.220 24452
25 25.8000-20 -2392
119.7500.1554
124 .6955
130.70000.00
44335
192,3400 23:300
(c) Synthetic binary SF, r = 2, DAG F1
Method
Dataset Settings
d=10, 1=0100
d=10. n-Jh
d=10,
d=10, 100k
d-50, 200
di350 n-1k
d=30, z=
d=550
105k
DAGMA (Linear)
0.1383
0:09778
D.1 1.1899
0.1697
0.1529
0.3451
0.1372
DAGMA (Nonlinetr)
0.3603-0
0.3587
0.3766-
0.3446
0.3534
NOTEARS (Linear)
0,2089
0.1872 1000
0.1648
0.1801
0.1849 _c30t
NOTEARS (Ncolineur)122 152]:
0.0788
400453
0,0265
0.0402
8,0450
0.0281
-1.00
0.0458 AUISS
DAGPA (Ours)
0,1689 20280
0,2409 1008
0,2562
0.2682 10000
0.0061
18383
0.0005
0.0163
c800
0.0410
37

---

Table 4: Experimental results on Synthetic binary SF, r 4
(a) Synthetic binary SE, r = 4, CI-MC
Method
Dataset Settings
d=10, 100
d-10 nm/A
d=10 10k
d-20, 20, 2650
dl-50 100
2-50, 2m/A
de50, p0
J--00
100K
DAGMA (Linan) [4):
0.2029
0.1665
0.1764 1:05
0.1886 Liwc
0.0527
0,0516
10a90
0.0654
L1a5
0,0543
DAGMA (Nonlinear) (4]:
0.4897
0.3381
0.3735
0.3296
0.3055
0.0811
0,0708
0
NOTEARS (Linzar) 55];
0,2459
0.1940
0.2024
0.0613
0,0545
0,4494
0,0577
NOTEARS (Nonlincar) (52)
0.1762
0,155
0.1734
0.1763
0.0494
0.0456
0.0417
0.0462
0,2012
0.3465
019938
.0.1133
0.0046
[35]:
0.4374
0.9950
0.040
00,0552
(K-1)1188
0.4708 462203
3.907
<.1443 440024
0.0175 4005
-0.0044
(18)
0.4528 J827
0.9331
0.1893 J00030
0.0266
0.0014
100000
DAGPA (Ours)
6.7190 11:06
0 19699
0,9609 2:300
0.9927 0005
0.0907 +0000
0.10M
0.1851 0010
0.3492 (114
(b) Synthetic binary SF, = 4, CPDAG F1 Arrowhead
Method
Dataset Settings
d=10 s=:100
d=10, m=1k
d=10, m=
d=10 d=100,==0000 1=00
d=50 n=
di:50 m=
d=50. x=/0
d=50 m= 1066
DAGMA (Lincar) [4]:
0.0000..000
0,0000.0000
124000
000000
0.0000
0.0251 433111
0.0020 1:00
0.0000.1000
0.0010 4
DAGMA (Nonlinerr [4]
0.1100
0.114|
0,0755
0.0713+ames
0.1244
0.0727
0.0043 -1943-
0.0954
NOTEARS (Lincur)
0,0125
0,0000
0,0000
0.0000
0.0191
0,040
0.0011
0.0051
NOTEARS (Noollnncar) (52)
0.0105
0.003-
8.080
0.0051
0.0000
[6]:
0,0253
0.2427 1t073)
241329
:1897
0.5297
PC 1331:
10.34
0.3093 .0.0000
0.3147 491200
0.0330x82222
0.2456 48852
0,3968 -1350
XPC (8-1) [18]
0.0311 -0000
0.0578
-00601
0:1R13
0.1209
0.0360
0.0333-58144
0.0160 -1854
0,011
+001
APC (k=2) [18]
0,0159 _0350
0,0452
0,1128
10083
0.1637
140004
0.0100 2:3157
0.0127 L28150
0,0430_EE22
0.0590
DAGPA (Ours)
0,1558 0200
0,2338
0,2589
10.100
0.2131 141222
0.0172
0.0276
0,0031 -C
0.0476
100077
(c) Synthetic binary SF, r = 4, CPDAG Fl Skeleton
Method
Dataset Settings
10 100
10 m=
d= nz 70k
d= 10, 1000
d= m= 700
d=500 1m
d=500 mz 10k
d= am 1006
DAGMA (Linear)
0.1182
0,0846
0.0957
010722
0.0637
DAGMA (Nonlinctr)
0.3630
0.2010-011
0.2995
0.3001
0.2179
0.3140
(Linear)
.00144
Juae
1osets
1
0.1010
NOTEARS (Noolineur) [52]:
100884
0.0791
2tum4
29.97
18:244
0.0293 21539
0.0294-6000
10003
[o]:
0.1416
400000
1.7679
000014
0.4424
43300
0.1747
0/410 23005
0.7002
0-4622
PC [33]:
0.3141
0,5606
-00631
-00000
0.7840
0.5763-158
0.7387
0.7675
0004)
APC (k=1) (18):
0,3237 _e.000
0,6298
0,8036
100074
0,8551
0.2113
0.5345
11825
0.6921
0.6553
100Z14
KPC (R-2) (18)
0.3141 -0.156
0.5836
60060
0,7533
0.3036
2005 8220
015344
0.7402
0000
0.7385
00035
DAGPA (Ours)
0,5007 .03912
0,5839 100094
0,6645 110096
0.6512 193560
0.0135
0.0677 118775
0.0726_1150
1100
0.0956
L0003S
(d) Synthetic binary SF. r = 4, CPDAG SHD
Method
Datusel Settings
d=10 100
d=10, 12
2=10 10
2-10 om00
d-000 cm 100
dw:0 1=
d-35 mw 70k
20550 7=00
DAGMA (Linze) [c]:
30,0000
10:0
00.00000
195. (0000 1120
6000 12006
DAGMA (Noulinsar) [4]
37,1000 1100
36 300000 25476
36,1000
43.000
221-4000-00 11009
190. 500.5000 61.0
6000
4200 43300)
Linerr
NOTEARS
(52):
30.5000.12011
39.5000
39 50.500
39.500 13,5000
199
199 199,4000
199.4000
GRS
$1.30 11.30300
$5.00.01-700
15 40-0000-000
1000
197.0000
167.4000
1ST (0000
74 11000000000000
PC [33]
10 19:0000::1:0
$4.7000
$30000 1196
(2.0000)
2066 0000
151 151.4000
155 60000-1110000
172.7000 12:7000--11000
LPC (1-1)
39 .130 1131
3.00 33.9000
32.7000 47150
31, 1000 .4300
207 8000 21109
155.55 158,50000 1550
190 5000 13.185
270,2000 -20200
APC (18)
30.0000 1m00
$4.3000
mm
53 53.5000 64:0
32. 2000 1904
203 7000 2101
151 151.8000 mm2
167. 9000 10:00
201 50.0 m110
DAGPA (Qhrs)
47.1250
96,0000
1000
-400
233.0000 11196
344.4900-
272.6644 -71971
(e) Synthetic binary SF r 4, DAG
Method
Datasrt Settings
de 10 1=
d=10, m=12
d=10 m: 70k
d=10,
100k
d=50 700
d0 w=
d=550, ez
d= am 1066
DAGMA (Lincar) [4]:
0.0082 Le55s
0.4433
400271
120517
400000
0.0651 433443
0.0673
(10330
20125
013406 1006
OLU66E and
DAGMA (Nonlinear) [4]
0.2486 eres
1.1938
0.1950 -m
0.1977 -1961
0.1565-0000
0.1362 -1644
0.1359
01455
NOTEARS (Linear) (51)
0.1008
0,0832
0,0508
0.0879 14894
0.0845 18:219
0.0738
0.0628
0.0763
NOTEARS (Nculiaear) (52).
0,0369
0.0006
0,01276
0.014
8.012
0.0118-1.0
DAGPA (Oues)
0.2219
0.2752
0.3358
0.0235+63141
00034-54171
0.0532
38

---

integration direction is promising and could leverage complementary strengths of both paradigms, the normalization challenge for nonlinear models and comprehensive empirical evaluation remain important future work.
G.2 Alternative Conditional Independence Measurements and Tradeoffs
Our framework uses p-values from statistical conditional independence (CT) tests as soft measures of independence strength. This section discusses the motivation for CI-based measurements and compares different CI measurement approaches.
G.2.1 Motivation: Robustness and Interpretability
Robustness to Data Preprocessing. Many statistical CI tests offer inherent robustness due their scale-invariance properties. Standard CI tests remain reliable across different preprocessing pipelines: Fisher-z test statistics invariant under línear transformations since partial correlations PXY preserved when variables are standardized; chi-square tests operate on frequency counts which are unaffected by scaling; kernel-based tests use normalized embeddings providing ceale-invariance by design. We empirically verified this; DAGPA maintains CI-MCC scores with low variance of se 0.12 across unnormalized, standardized, min-max normalized, and robust-scaled versions the same dataset (n = 100, d = 10 ER-2 synthetic data), demonstrating stability that is valuable when data sources have heterogeneous scales or different preprocessing conventions.
Interpretability via Direct Hypothesis Testing, Cl-based measurements provide direct interpretabil- ity through testable hypotheses. Each structural decision corresponds specific CI tests with p-values, enabling transparent explanations. For example, when removing edge X + Y, we can justify: "Fdge removed because unconditional test shows X Y with p-value 0.001 and conditional test shows X Z with p-value 0)82,, indicating X influences Y likely through but does not have direct causal link Y Domain experts can independently verify these claims by examining raw data, running alternative tests, or checking whether distributional assumptions hold. This trans- parency is particularly valuable in high-stakes applications (medical diagnosis, policy making) where algorithmic decisions require human oversight and regulatory approval.
G.2.2 Overview of CI Measurement Methods
Different CI tests present tradeoff' across multiple dimensions. Table 5 summarizes key characteris- tics:
Table Comparison of measurement methods.
"Small-n" indicates performance with n < 200.
Method
Assumptions
Complexity
Small-n
Nonlinearity
p-value
Fisher-z [11]
Gaussian, linear
(nn)
Excellent
Poor
Native
Chi-square
Categorical
((n)
Good
Arbitrary
Native
KCIT [48]
Nonparametric
O(n²)
Poor
Excellent
A Approximate
RCIT [35]
Nonparametric
On
Moderate
Good
Approximate
CMI (k-NN)
Nonparametric
Onllog log
Poor
Good
Requires transform
Fisher-z Test: Tests whether partial correlation pxvız equals zero using transformation 2
|/n 7 3log 1-6 1 1). Strengths include excellent sample efficiency (reliable with 2 50 for 7 S 1), (n) computational cost, and native p-values. Limitations: designed for Gaussian data and primarily detects linear relationships. Our default choice for continuous data in low-sample regime.
Chi-square Test: Tests independence via contingency table using x² Σ., (O:-E.;)² Strengths
include distributional assumptions, detection of arbitrary associations, and good small-sample performance. Limitations: only applicable to categorical variables, requires sufficient cell counts (Eij 2 5), and exponential complexity with |Z]. Our default for binary/categorical data.
Kernel CI Test (KCTT): Tests independence via kernel embeddings using Hlilbert-Schmidt Inde- pendence Criterion. Strengths include detecting arbitrary nonlinear dependencies and no parametric
39

---

assumptions. Limitations: O(n²) complexity. requires >> 500 for reliable results, and kernel band- width selection affects performance. Recommended for large datasets (n 1000) with suspected nonlinear relationships.
Randomized CI Test (RCIT):
Approximates KCTT using random Fourier features, reducing plexity to (n n).. Provides middle ground between Fisher-z (fast, linear) and KCIT (slow, nonlinear) for moderate samples (500 < n 5000).
via k-NN or kernel density estimation. Provides ceemmmtiiontthhorreeee interpretation but lacks native p-values (requires transformation like p = exp(-3 3 Π)) and needs n > 1000 for reliable density estimation.
Conditional Mutual Information (CMI): Measures (X; = E log P(X.VIZ)
p(XT2)p((12)
Tradeoffs and Selection Guidance. For causal discovery with d variables, total CI testing requires (dF) order-O tests and O(d)) order-1 tests. With d 50, n 1000: Fisher-z completes in ~30 seconds, RCIT in ~5 minutes, KCIT in 445 minutes. Sample efficiency also varies: Fisher-z achieves high power with < 200 for linear relationships, while kernel methods require n 500 but eventually match or exceed Fisher-z for nonlinear cases. Our framework is measurement- agnostic- any function producing [0, 1] scores can be substituted via our API:
model (AAFF(ccittestt''iiheerr
Default
model DAGPA (ci test=t
For nonlinear data Custom implementation
model DAGPA MAGPA(ci_testtee (ec__ttest-tyyycussoteetttt
We recommend Fisher-z or chi-square for most applications due to sample efficiency and com- putational scalability, with kernel methods reserved for large datasets with confirmed nonlinear relationships. Practitioners should conduct sensitivity analysis across multiple tests when possible, as edges appearing consistently provide higher confidence regardless of specific test assumptions.
H Licenses
In this work, we evaluated our method on two publicly available causal discovery benchmark datasets:
Sachs Dataset: The Sachs dataset contains simultancous measurements of 11 phosphorylated proteins and phospholipids derived from thousands of individual primary immune system cells, subjected to both general and specific molecular interventions. This dataset was originally published by Sachs et al. (2005) and is widely used benchmark in causal discovery research. The dataset is publicly available through multiple repositories including bnlearn [30], other causal discovery toolboxes.
Availability: The dataset is publicly accessible for research purposes through various causal discovery software packages and repositories.
The LUCAS (LUng CAncer Simple set) dataset [14] is synthetic benchmark dataset consisting of 12 binary variables and 2000 instances, representing different causal rela- tionships in medical diagnosis problem for identifying patients with lung cancer. This dataset was created as part of the Causality Workbench project to provide standardized benchmarks for testing causal discovery algorithms.
LUCAS Dataset:
The dataset is publicly available for research and educational purposes through the Causality Workbench repository.
Availability:
Usage Declaration: Both datasets were used in accordance with their respective terms of for academic research purposes. No additional permissions were required for their use in this study.
I Societal Impact Statement
Our differentiable d-separation framework for causal discovery has potential positive impacts in healthcare (identifying causal factors disease), public policy (evaluating intervention effective- ness), scientific discovery (understanding complex systems), and algorithmic faireess (distinguishing causation from correlation in decision systems). However, several risks must be acknowledged: (1) incorrect causal discoveries could lead to harmful interventions if implemented without domain expert validation; (2) causal discovery in sensitive domains may raise privacy concerns, necessi- tating differential privacy techniques; (3) when applied to historically biased datasets, discovered
40

---

relationships might reflect and perpetuate these biases rather than ground truth; and (4) computational requirements might limit accessibility well-resoureed institutions. We recommend multiple mitiga- tions: returning diverse candidate structures rather than single models, requiring expert validation before implementation, implementing
techniques with sensitive data, examining discovered relationships for bias, and developing more efficient implementations to improve acces- sibility. Causal discovery tools require responsible application and domain expertise, especially in high-stakes domains where incorrect inferences could lead to harmful consequences.
41

---

NeurIPS Paper Checklist
1. Claims
Question:
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction claim that contribute novel framework for causal discovery that bridges constraint- based and mnnimmous-optimization score-based methods with promising empirical results. These claims accurately reflect our work's scope and contributions as evidenced in: (1) Section 3, where we develop the novel differentiable d-separation framework with theoretical guarantees and iyvoolattoo-bsee probabilistic interpretations; (2) Section 4, where we demonstrate one practical instantiation of this framework into the model we named DAGPA; and (3) Section 5, where empirically validate that our approach achieves competitive performance compared to established methods across multiple datasets and settings.
Guidelines:
The answer NA means that the abstract and introduction do include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
it fine to include aspirational goals motivation as long as it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We provide a dedicated Limitations paragraph in Section 7 that thoroughly discusses multiple aspects of our work's ' constraints and outlines specific future rescarch directions to address them. Throughout the paper, we are transparent about the scope our empirical validation and the provisional nature of our current implementation compared to the theoretical framework's s potential.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g., independence assumptions, noiseless settings, model well-speciffcation, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets with few runs In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might that reviewers discover limitations that aren't it acknowledged in the paper. The authors should use their best judgment and recognize
42

---

that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: In our paper, all assumptions are explicitly stated before each theorem and lemma, and the complete proofs are provided in Appendix A, with every item properly numbered and consistently cross-referenced throughout the paper. In addition, with intuitive explanations connecting the formal statements to the broader framework, and for the most important results, such as Theorem 3.2, we provide proof sketches in the main text that highlight the critical steps and insights.
precede each theoretical result
Guidelines:
The answer NA means that the paper does include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross-referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experi- mental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We provide complete reproducibility through detailed algorithm pseudocode for DAGPA (Appendix C), code and data release (Appendix D.1), dataset creation and preprocessing procedures (Appendix D.2), and comprehensive hyperparameter specifications for both our method (Appendix D.4) and baselines (Appendix D.5). Additionally, we document the computational resources used (Appendix D.6) and fully explain our evaluation metrics (Appendix D.3), ensuring all aspects of our experimental pipeline can be independently reproduced.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible important, regardless of whether the code and data are provided not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary either make it possible for others replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access to a hosted model (c.g., in the case of a large language model), releasing of model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
43

---

(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to reproduce the model (e.g,, with an open-source dataset instructions for how construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which authors welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users). but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We provide complete code and data access through an anonymous GitHub repository (linked in Appendix 1. that includes all implementation files, datasets (both synthetic generators and real-world data), and evaluation scripts necessary to reproduce our results. The repository contains detailed documentation with exact environment setup instructions (via Conda environment files), step-by-step commands to run all experiments, data preprocessing pipelines, and scripts to generate all figures and tables presented in the paper.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please the NeurIPS code and data submission guidelines
Caleess/CoSSubbiie.c issi onPolicy for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed the results. See the NeurIPS code and data submission guidelines (https / /nwppsippo/o//// cttpssccccpbb/c///ppp.c guides/ htte//ooeeuubiiii weabbbmmsssooopoiiioooooooieeeee onPolicy) for more details.
reproduce
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification:
The main text (Section 5) provides essential experimental details including dataset characteristics, evaluation metries, baselines, and the overall experimental design, while Ap- pendix D contains comprehensive information on hyperparameters, model selection criteria, data creation and preprocessing, and computational resources. We explicitly describe our parameter selection process and provide justification for key design choices to ensure results can be fully understood and contextualized.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core of the paper to level of detail that necessary to appreciate the results and make sense of them.
44

---

The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: Our primary results are presented empirical cumulative distribution functions (CDFs) over multiple (10) independent dataset instances per configuration, fully displaying perfor- mance variability across different data realizations and allowing for direct statistical comparison between methods. For tabular results, we report mean values with standard deviation error bars (1-sigma), clearly labeled as such, and calculated using standard statistical formulas across the independent trials. The sources of variability (different random DAG structures and data samples) explicitly described in Section 5, and all statistical comparisons are appropriately referenced in the text.
Guidelines:
The answer NA means that the paper does include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, least for the experiments that support the main claims the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1 I-sigma error bars, but one should state it. The authors should preferably report 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed reproduce the experi- ments?
Answer: [Yes]
Justification: Yes. We detail the computation resources used for each experiment settings in Appendix D.6.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e-g., preliminary or failed experiments that didn't 't make it into the paper).
Code of ethics
Question: Does the research conducted in the paper conform, every respect, with the NeurIPS Code of Ethics https: //werr/////.p...eeeee http://s......pceecaicccccciiaia
45

---

Answer: [Yes]
Justification:
Yes. This research fully conforms to NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there a special consideration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: We provide a comprehensive Broader Impact Statement in Appendix that discusses both potential positive impacts (in healthcare, public policy, scientific discovery, and algorithmic fairness) and potential negative impacts (misplaced trust in discovered causal relationships, privacy concerns, potential for reinforcing biases, and accessibility challenges). For each potential negative impact, we also suggest specific mitigation strategies that practitioners can implement when applying our method.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no socictal impact why the paper does not address socictal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g.. disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deploy- ment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to par- ticular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, is legitimate to point out that improvement in the quality of generative models could be used generate deepfakes for disinformation. On the other hand, is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put place for responsible release of data or models that have high risk for misuse (e.g. pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Our work presents a methodological framework for causal discovery that does not pose the high-risk misuse concerns typically associated with generative models or scraped datasets. We only synthetic datasets (generated with provided code) and well-established, properly cited benchmark datasets (e.g. Sachs) from trusted scientific sources. Our method does not learn representations that could be misused for generating harmful content, creating fake profiles, or similar high-risk applications that would necessitate special safeguards beyond standard research best practices.
Justification:
Guidelines:
46

---

The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse dual-use should be released with necessary safeguards to allow for controlled use the model, for example by requiring that users adhere usage guidelines or restrictions access the model implementngg safery filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We properly cite all original papers for datasets and code libraries used in work. In Appendix H, we provide detailed license information for the Sachs dataset (the only external real-world dataset used) and all third- parrd-patty libraries/packages employed in implementation, including version numbers and URLs to their repositories. For synthetic datasets, we specify that they are generated by our code which is released under MIT license.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi wwwode..eed...... com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.
For existing datasets that are re-packaged, both the original license and the license the derived asset (if it has changed) should be provided.
this information not available online, the authors are encouraged reach out the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We release our implementation of DAGPA via an anonymized GitHub repository (linked Appendix D.I) with comprehensive documentation including: installation instructions, usage examples. API documentation, parameter descriptions, trainingoppimizzttoonsssm/pppppp procedures, and known limitations. The repository includes a structured README file, detailed code comments, and clear MIT license statement.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their sub- missions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
47

---

Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
This work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may he required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applica- ble), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, non- standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: This work does not involve LLMs as any component our core methodology, algorithms, experiments. LLMs were used solely for writing assistance, editing, and LaTeX formatting, which per the NeurIPS LLM policy does not impact the scientific methodology or originality of our work and therefore does not require declaration beyond this checklist.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs any important, original, or non-standard components.
Please refer our L.LM policy lattps://neurif: http:://eerrrp..cc//offfeeeeeeee should or should not be described.
for what
48

---

Out
GB5
NOTEARS ILIWANO) NOTEARS (hivetioar) PC
DAGMA (lireer)
(k-11 MPC
Ouea
GBS
NOTEARS (Lizert
DAGNA (Noalmor??
(bm21 APC
BAGMA hheeerr
KNC [Rx
NOTEARS (Noctreem
DAGMA (Nertiinnn
2010/10
Enginnca
Iregitial
XCC (CLMMC)
oeder Statement (Ci
1t
(a) Sachs CI-MCC
(b) Lucas C1-MCC
DAGNA (Noattoar
0em21
BAGMA (Nonitteac)
APC (km
Tregirscal
CPDMU Artwssadd
Low-order CI Statement MOC (CI)EL)
(c) Sachs CPDAG Arrowhead
(d) Lucas CPDAG FI Arrowhead
GES
NOFFRARS (Liveer
Dere
GBS
NOTTARS (Urrvrr
DaCAS ILazeut
MPC (um))
NOTTARS /Nosllaouut
DACMA clincut
APC (ue
NOTEARS (Nvclise
(NonEwwwr)
21
BAGMA (Nenitzeart
B0C 31
Engginial CDF
Ereginool
(e) Sachs CPDAG FI Skeleton
( Lucas CPDAG F1 Skeleton
NOTEARS
Expitooal CDF
DAG
DAG
(g) Sachs DAG Fl
(h) Lucas DAG Fl
Figure 9: Full results Sachs [28] and Lucas [14] real-world dataset.
49


Reviews:
[{"review_id": "bk9Rb4qxv0", "reviewer": "Reviewer_TmCS", "summary": "This paper studies a new \"style\" of gradient-based causal discovery approach, where the gradient based optimization incorporates CI constraints as in constraints-based CD methods. This differs from the usual, existing score based approach, which encodes DAG in the constraint. The paper derives theoretical sounds showing soundness. The paper also demonstrates empirically that the method developed works particularly well in low sample settings, showing it is more robust in such settings.", "questions": "1. Can the authors comment more about the convergence of the (joint) optimization of the proposed objective function? This will make the paper more complete.\n\n2. Are there other replacements for the p-values, which are used as heuristic conditional independence measures?\n\n3. Is it correct to say that the approach developed in this paper can complement any existing score based approach? Specifically it can be added on top of any gradient-based method in that one can jointly optimize this paper's criterion and the criterion of the existing score based method. If yes, would this be better and if no, why not?", "limitations": "Please see weakness section.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength: The paper's motivation is interesting, aiming to produce a differentiable, constraint based causal discovery method. This could potentially lead to new lines of work just like how NOTEAR opened up a lot of directions and areas for improvement. I thought the paper did a nice job building up from d-separation and S_A to a main theoretical result (Lemma 4.2). The experiments do compare fairly against known baselines like NOTEARS, GES and PC (but the experiments are specifically evaluating for recovery of d-separation statements implied by the true underlying causal structure that generated the data, which is what the method optimizes).\n\nWeakness: With that said, I think more could be said about the optimization guarantees when optimizing the five multi-task loss functions (which was done in NOTEARS). I think the paper could also say more to motivate the necessity of having explicitly encoded constraint based approaches (what are some applications where this is important and NOTEARS is not sufficient).", "was_revised": true, "final_justification": "I think the paper is new and interesting, but there is some unknown about its convergence properties so it's a bit unclear if it can actually converge to the desirable optimum. As I said before, my overall rating is positive and I think the paper should be acepted.", "rebuttal": "We sincerely thank Reviewer TmCS for their positive and insightful feedback. We are encouraged that the reviewer found our motivation \"interesting\" and believes our work \"could potentially lead to new lines of work just like how NOTEARS opened up a lot of directions.\" We appreciate the thoughtful questions, which will help us significantly improve the clarity and completeness of our paper.\nBelow, we address the specific weaknesses and questions raised.\n\n> Q1. Can the authors comment more about the convergence of the (joint) optimization of the proposed objective function?\n\n**A1.** Thank you for this insightful question. We can provide a convergence analysis for our proposed objective function under certain conditions, which we can include in the revised paper. Below, we outline the key theoretical and practical arguments supporting the convergence of DAGPA.  \n\nTheoretically, DAGPA leverages the **Discrete Langevin Proposal (DLP)** sampler [1] to perform posterior sampling of DAGs from the true posterior distribution, whose energy function is defined as the sum of our five multi-task losses. Under the DLP framework, we can establish asymptotic convergence to the true posterior, provided that gradients of the objectives are added rather than processed via PCGrad. This convergence result stems from two key properties of DLP:  \n1. **Metropolis-Hastings (MH) rejection step**: Ensures that the sampler respects the true posterior distribution.  \n2. **Ergodicity of the Markov chain**: Guarantees that every possible DAG has a non-zero probability of being sampled, satisfying the conditions for asymptotic convergence as per [1].  \n\n\nAdditionally, **Lemma 4.2** in our paper establishes consistency, showing that as the sample size increases and the soft-logic approximation becomes exact, the true underlying DAG minimizes each loss component and thus has the highest posterior probability.  \n\nPractically, we adopt several strategies to enhance mixing (but of which we cannot evaluate their effectiveness):  \n1. **Acyclicity constraint**: We use the log-determinant characterization from DAGMA [2], which offers better-behaved gradients compared to NOTEARS, particularly when cycles are present. This choice ensures stable optimization even in complex graph structures.  \n2. **Gradient conflict resolution**: PCGrad [3] is employed to navigate conflicting gradients across tasks. Unlike naive loss combination methods, PCGrad projects gradients onto a consensus direction, provably improving convergence in common scenarios [4]. This is important for multi-task objectives, where gradients may pull the model in opposing directions.  \n3. **Enhanced exploration**: DLP enables parallel edge updates during sampling, significantly accelerating exploration of the parameter space. This is particularly valuable for causal discovery due to the non-uniqueness of solutions (e.g., Markov equivalence classes).  \n\nWhile quantifying an exact convergence rate under our multi-task losses remains a theoretical challenge, we address this by combining rigorous Bayesian guarantees with battle-hardened optimization practices. We will expand on these points in the revised paper and encourage future work to further analyze the complex loss landscape of CI-statement losses. We appreciate the reviewer’s feedback and look forward to incorporating these discussions.  \n\n\n> Q2. I think the paper could also say more to motivate the necessity of having explicitly encoded constraint based approaches\n\nA2. We thank the reviewer for their insightful feedback. We will revise the paper to better emphasize the unique advantages of our approach, as outlined below.\n\n1. **Differentiable Learning of k-Markov Equivalence Class**: Currently, there is no differentiable-learning-based approach similar to NOTEARS that can learn the k-Markov equivalence class in the low-sample regime, where we can only rely on low-degree CI tests since high-order ones are unreliable. We argue that this is a significant shortcoming, since real-world tasks often have a limited number of samples, and the k-Markov equivalence class theoretically quantifies the ceiling of what can be learned under this scenario.\n\n2. **Interpretability via Conditional Independence (CI) Tests**: Our method leverages CI tests as its core mechanism, which can offer greater interpretability in some scenarios when compared to other score-based methods that rely on likelihoods. This could be particularly valuable in real-world applications where interpretability is critical.\n\n3. **Robustness in Low-Sample Regimes**: Unlike traditional score-based approaches, which often struggle with overfitting or underfitting when fitting functional models to scarce data, our method (DAGPA) incorporates \"soft\" CI information derived from p-values. This design choice enhances robustness in low-sample settings, as evidenced by its superior performance at $n = 100$ compared to baselines (Figure 3a). This highlights a practical advantage of DAGPA in data-limited scenarios.\n\nWe appreciate the reviewer’s emphasis on the paper’s significance and will strengthen the introduction to explicitly articulate these points, thereby improving the overall motivation and clarity of the work.\n\n\n\n> Q3. Are there other replacements for the p-values, which are used as heuristic conditional independence measures?\n\nYes, indeed. Our framework happily accepts any CI measurement $p_{\\mathcal{D}}$ from input data as long as it satisfies the semantic that $p_{\\mathcal{D}} \\approx 1$ indicates (likely) independence, and $p_{\\mathcal{D}} \\approx 1$ indicates (likely) dependence. $p$-values fits this requirement; an alternative choice is the **conditional mutual information (CMI)** with some adequate conversion. \n\nSpecifically, mutual information and conditional mutual information are always non-negative quantities, with a value of 0 indicating independence and values greater than 0 indicating dependence. To convert to the aforementioned definition of the CI measurement $p_{\\mathcal{D}}$, one could adopt some form of indicator function. One candidate choice is the following function:\n$$\nf(x) = 1 - e^{-x / l} ,\n$$\nwhere $l$ is a hyperparameter controlling the sharpness such that when $l \\to 0$, $f(x) \\to 1[x > 0]$ approximates the hard indicator function that checks if input $x$ is greater than 0. Thus, given a raw CMI measurement $\\text{CMI}\\_{\\mathcal{D}}(x, y, \\mid z)$ for a triple of variables $x, y, z$, $f(\\text{CMI}_{\\mathcal{D}}(x, y, \\mid z))$ becomes a valid CI measurement ready to be used in our proposed framework.\n\nAlthough conceptually such an approach makes sense, during earlier method development we found that these converted CMI values are somewhat unreliable, as the values corresponding to the true independent statements vary a lot. Using these converted CMI values as input to the algorithm makes the method fragile and decreases the chances of sampling good DAG solutions. In addition, it introduces yet another hyperparameter $l$ that one needs to tune. In contrast, the p-values varied less and were more stable, and we were working with one less hyperparameter. Thus, we decided to adopt p-values for the paper. But CMI remains a promising avenue that we intend to explore further.\n\nWe thank the reviewer for bringing up this point, and we will add this discussion in the revised Appendix.\n\n> Q4. Is it correct to say that the approach developed in this paper can complement any existing score based approach?\n\nYes, we believe such a hybrid implementation is possible. The main challenge to address is defining how a weighted adjacency matrix $\\mathbf{W} \\in [0, 1]^{d \\times d}$ with edge weights in the range $[0, 1]$ can be derived from the model parameter $\\theta$ in a way that is compatible with both our differentiable d-separation framework and the score-based method. Particularly to our method, we require the edge weights to fall in $[0, 1]$ as it needs to be valid probability measures.\n\nFor instance, to convert the model parameter to our framework’s probabilistic edge weights, it should be possible to apply a sigmoid activation $\\sigma(\\cdot)$ for the NOTEARS linear model’s parameter $\\theta$. For the nonlinear version of NOTEARS though, its graph weighted adjacency matrix, obtained by $\\text{sqrt}({\\theta^{(1)}\\_j}^T \\theta^{(1)}\\_j)$ ($\\theta^{(1)}$ the first MLP layer weight) consists of unnormalized, nonnegative values. Thus, one needs to find a proper way to normalize such a weight into range $[0, 1]$ to satisfy the probabilistic semantics. The rest is to simply add the prediction loss, such as the MSE loss $\\mathcal{L}_{\\text{MSE}}(\\theta, X) = \\frac{1}{2n} \\| X - XW \\|^2_2$ for the linear NOTEARS model, as an additional task loss and apply the same multi-task optimization routine of DAGPA.\n\nDue to the limited time of the rebuttal period, we were not able to conduct comprehensive analysis of a concrete implementation of such a hybrid method. We will however emphasize this discussion in the revised main text and call for future work to investigate this promising research direction.\n\nReferences:\n\n[1] \"A langevin-like sampler for discrete distributions.\" PMLR 2022\n\n[2] \"Dagma: Learning dags via m-matrices and a log-determinant acyclicity characterization.\" NeurIPS 2022.\n\n[3] \"Gradient surgery for multi-task learning.\" NeurIPS 2020.\n\n[4] \"Revisiting scalarization in multi-task learning: A theoretical perspective.\" NeurIPS 2023.\n\n[5] \"Learning sparse nonparametric dags.\" PMLR 2020"}, {"review_id": "2UfJFTqBM3", "reviewer": "Reviewer_pEMd", "summary": "Differentiable Constraint‐Based Causal Discovery introduces a novel hybrid framework that marries the rigorous conditional‐independence testing of constraint‐based methods with the flexibility of gradient‐based optimization. The authors first show how low‐order d‐separation and d‐connection statements can be encoded as first‐order logical formula over graph reachability, then relax these discrete predicates into continuously differentiable “percolation” measures using probabilistic soft logic (LogLTN operators). This yields differentiable lower‐bounds on expected d‐separation/d‐connection that serve as losses alongside a log‐determinant acyclicity regularizer, forming the DAGPA algorithm. Empirical results on synthetic Erdős–Rényi and Scale‐Free graphs and the real‐world Sachs protein‐signaling dataset demonstrate DAGPA’s strong robustness in small‐sample regimes—surpassing both traditional constraint‐based (PC, k‐PC) and score‐based (NOTEARS, DAGMA) baselines—and competitive performance when ample data are available.", "questions": "I recommend adding the final optimization objective to the main text—showing concretely how the differentiable d-separation/d-connection losses and the acyclicity regularizer combine into a single formula will greatly improve readability.\n\nIt would be helpful if the authors could analyze the computational complexity of DAGPA and report empirical runtimes on representative benchmarks in the main paper. Are there opportunities to accelerate the method (e.g.\\ parallelizing percolation calculations, exploiting sparsity in the graph, or using approximate reachability)?\n\nThe experimental metric plotted in Figures 2 and 3 (“Low-order CI Statement MCC vs. Empirical CDF”) is not self-explanatory. Could you share more details on this part.", "limitations": "YES", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "**Strengths**\n\n1. The paper is exceptionally well written. Despite introducing many novel concepts and intricate notation, the authors clearly explain each idea and present their results in a coherent manner. To further aid comprehension, I suggest including a concrete example alongside each formal definition—particularly Definition 3.1—which currently feels under‐explained.\n\n2. The motivation and core idea are both compelling and original. The authors demonstrate how to embed the d-separation/d-connection criterion as a differentiable scoring mechanism within a gradient-based optimization framework. This novel integration clearly merits an accept.\n\n**Weaknesses**\n\nAs the authors acknowledge, the current approach only handles d-separation/d-connection for singleton conditioning sets. Extending this to larger conditioning sets would require exponentially many additional constraints, leading to severe computational challenges both in computing each individual score and in managing the objectives within the multi-task optimization.", "was_revised": true, "final_justification": "I find the paper presents a very interesting idea—personally, it caught my eye. The authors’ comprehensive rebuttal has addressed my concerns. This paper is solid 5, but I’m raising my score to a 6 as encouragement.", "rebuttal": "We are grateful to Reviewer pEMd for their positive and thorough review. We are delighted that the reviewer found our work \"exceptionally well written\" and the core idea \"compelling and original,\" and we sincerely appreciate their championing of the paper. The reviewer suggestions are all constructive and will certainly help us improve the paper's clarity.\nBelow, we address the reviewer's questions and suggestions.\n\n> Q1. I recommend adding the final optimization objective to the main text—showing concretely how the differentiable d-separation/d-connection losses and the acyclicity regularizer combined into a single formula will greatly improve readability.\n\nA1. We appreciate the constructive suggestion to enhance readability. In the revised manuscript, we will introduce a dedicated paragraph at the end of Section 4 to explicitly summarize all components of the multi-task loss function. While our method does not rely on a single unified formula due to its multi-objective nature, we will consolidate the individual loss terms into a centralized location for clarity. If space constraints prevent full inclusion in the main paper, we will provide explicit pointers to the detailed equations in the Appendix.\n\n> Q2. It would be helpful if the authors could analyze the computational complexity of DAGPA and report empirical runtimes... Are there opportunities to accelerate the method?\n\nA2. Thank you for this valuable suggestion. We will add a detailed analysis of complexity, runtimes, and acceleration opportunities in the Appendix. We will also emphasize the scalability limitation and point to promising avenues for future work in the revision.\n\n**Empirical Runtime Comparison**: The table below presents the wall-clock comparison of DAGPA against the baseline methods on the synthetic $d=10$ and $d=50$ experiment settings, which was the time these methods took to yield the results shown in the paper. We note that the wall-clock time of our method DAGPA shown below is an **overestimate** of the actual time it took DAGPA to reach the best resulting solution. This is because DAGPA uses Bayesian sampling (specifically Discrete Langevian Proposal (DLP) sampler) to continue exploring the DAG space. Rather than taking the last sampled DAG as the algorithmic output, we kept the algorithm running, kept track of all DAGs sampled at each step, evaluated them using the metric using only the input data (Appendix C3), and returned the best one. In almost all cases, the algorithm could have been terminated much earlier to return DAGs of equal quality to those shown in the paper’s experimental results.\n\nGPU: 24GB A10\n|  | PC | kPC (k=1) | kPC (k=2) | GES | NOTEARS (Linear) | NOTEARS (Nonlinear) | DAGMA (Linear) | DAGMA (Nonlinear) | **DAGPA (Ours)** |\n|---|---|---|---|---|---|---|---|---|---|\n| $d=10$ | 0.60 | 0.15 s | 0.32 s | 1.80 s | 1.36 s | 10.98 s | 1.26 s | 172.92 s | 1.5 h (1000 steps) |\n| $d=50$ | 2.98 s | 2.60 s | 1.97 s | 921.14 s | 19.66 s | 44.23 s | 4.35 s | 801.98 s | 10.5 h (500 steps) |\n\n\n\n**Computational Complexity:** The computational bottleneck of our approach (DAGPA) is the p-values computation for all low-order CI statements and the reachability computation required by the differentiable d-separation formulae.\n\nThe complexity of p-values computation (which we GPU-accelerate) depends on the specific choice of statistical independence test. Take the Chi-squared test for example. Given $n$ data points, computing each single Chi-squared p-value takes $O(n)$ as this is the time required for iterating over all data points and building the contingency. Thus, the overall time required is $O(d^3 n)$ as we obtain a total of $d^2 + d^3$ number of p-values for the unconditional and first-order conditional statements.\n\nFor the differentiable d-separation formulae: the reachability subroutine involves a maximum of $d$ recursive steps, and each all-pairs Bellman-Ford update step operates on $O(d^3)$ matrix entries (equation (7)). Thus, the total computation time is $O(d^4)$. The resulting all-pairs reachability matrix will be cached and accessed by the d-connection/d-separation score computation. Now, since in the 1st-order d-connection/d-separation scores (equation (3) and (5)), we need reachability on the node-deleted subgraph $\\mathbf{A}_{-z}$ for all conditioning node $z$, this requires a total of $d+1$ all-pairs reachability matrices, thus rendering the total runtime of this part $O(d^5)$.\n\nFor the 0th-order d-separation/d-connection formulae (equation (2) and (4)), computing $S^{(0)}\\_{\\mathbf{A}} (x, y)$ or $C^{(0)}\\_{\\mathbf{A}} (x, y)$ for each (x, y) takes $O(d)$ time, since it iterates over $d$ possible common ancestors and accessing the reachability cache takes $O(1)$. Thus, computing all-pairs 0th-order d-separation/d-connection scores take $O(d^3)$ time. Now, since later in the 1th-order formulae, we need all-pairs 0th-order d-separation/d-connection scores for the node-deleted subgraphs $A_{-z}$ for all conditioning node $z$, this adds another $d$ dimension. Thus, the total time is $O(d^4)$.\n\nFor the 1st-order d-separation/d-connection formulae (equation (3) and (5)), similarly, computing the result for each query triple $(x, y \\mid z)$ takes $O(d)$ time. Thus, getting the results for all triple of nodes takes a total of $O(d^4)$ computations.\n\nThus in summary, computing all-pairs p-values takes $O(d^3 n)$ computations, computing reachability matrix (and cache for later use) takes $O(d^5)$ computations, and computing 0th-order and 1st-order d-separation/d-connection scores each take $O(d^4)$ computations. The total will be dominated by the reachability computation, which is $O(d^5)$.\n\n**Acceleration Opportunities:** That being said, there are ample promising opportunities for acceleration.\n\n- First, we have to note that all the computations of reachability and d-separation/d-connection scores are matrix-vector operations, which can benefit from **GPUs**. This is what we did in our implementation, where we leverage PyTorch GPU tensor library for all such computations, avoiding any explicit for-loops, significantly improving the speed of our algorithm. \n\n- Second, in practice one can limit the reachability computation to only consider **paths of a constant length $k$**, reducing the time complexity from $O(d^5)$ to $O(d^4 k)$, **rendering the entire pipeline $O(d^4)$**. This is what we did in the experiments for the larger graphs with $d = 50$ number of nodes, where we limit the path length to $k = 10$. Our earlier observation during earlier method development suggest negligible performance degradation with significantly improved computational efficiency. \n\n- In addition, the all-pairs p-values computation can be massively parallelized and then cached and reused. This is particularly useful in the case where one wants to obtain multiple output causal graphs for the same input dataset to enhance solution robustness and exploration of the Markov equivalence class. In this case, the computation time for the all-pairs p-values step can be amortized as it can be reused in subsequent runs.\n\n- Finally, we expect future work to further examine how to incorporate sparsity assumptions into the graph to further reduce the time complexity. For example, instead of assuming each node can connect to all $d-1$ other nodes, one reasonable assumption is to restrict to a maximum of constant $k$ degree. In that case, the reachability computation subroutine computation can drastically speed up.\n\n\n\n\n> Q3. The experimental metric plotted in Figures 2 and 3 (“Low-order CI Statement MCC vs. Empirical CDF”) is not self-explanatory.\n\nA3. We appreciate the reviewer’s feedback regarding the clarity of the metric plotted in Figures 2 and 3. We recognize that this visualization may not be immediately intuitive to all readers. In the revised manuscript, we will provide a more detailed explanation of this metric.\nThe **Low-order CI Statement Matthews Correlation Coefficient (CI-MCC)** metric, illustrated in **Figure 1**, evaluates how well a model's predicted graph aligns with the ground-truth graph in terms of the low-order conditional independencies they imply. The process is as follows:\nFor a predicted graph $G_{\\text{pred}}$​ and the ground-truth graph $G_{\\text{true}}$​, we generate two lists containing all possible 0th- and 1st-order d-separation statements (e.g., $A \\perp\\perp B$, $A \\perp\\perp C \\mid D$, etc.).\nWe treat this as a binary classification task: the goal is for $G_{\\text{pred}}$ to correctly classify which d-separation statements hold true in $G_{\\text{true}}$​.\nWe then compute the standard **Matthews Correlation Coefficient (MCC)** on the resulting confusion matrix (TP, TN, FP, FN). MCC is a robust metric for binary classification that is well-suited for cases with imbalanced classes (as the number of true independencies vs. dependencies can vary greatly in practice).\nThe plots in Figures 2 and 3 show the **Empirical Cumulative Distribution Function (ECDF)** of these CI-MCC scores over multiple runs, which visualizes the distribution of performance for each method.\nWe will revise the first paragraph of Section 5 and the captions for Figures 2 and 3 to include this clear, step-by-step explanation of the CI-MCC metric and the ECDF plots.\n\n\nReferences:\n\n[1] \"Gradient surgery for multi-task learning.\" NeurIPS 2020."}, {"review_id": "5Asd2gpXoo", "reviewer": "Reviewer_JM3s", "summary": "This paper proposes a novel method to incorporate \"(conditional) independence statements\" into a differentiable causal discovery framework. The authors first express d-separation $S_A(x,y|z)$ using logical notation (Definition 3.1) and then introduce a continuous relaxation of this definition  $\\tilde{S}_W(x,y|z)$ through soft logic. They establish a connection between the soft version and the \"expected\" d-separation statements, enabling the integration of causal inference with differentiable optimization techniques.", "questions": "How many CI statements are needed? If you need test all pairs, PC would be much faster? I can’t imagine a practical scenario where testing all pairs is feasible, especially when the number of nodes is large. This seems to contradict the motivation behind methods like NOTEARS, which aim to reduce computational complexity.\n\n\nComments:\n1. The notation $\\independent_D$ is somewhat unfamiliar to me. Typically, independence is discussed in the context of a joint distribution $P$, so this notation may need further clarification.\n2. Is Theorem 3.2 a theorem or a definition？As far as I know， d-separation can be directly read off from the graph $A$ and is denoted by $\\independent_A$. From my perspective, Thm.3.2 and Def.3.1  are the same thing.\n3. Lemma 3.4 appear to be redundant in this draft and could be omitted to improve clarity and conciseness.\n4. It took me some time to understand the exact method for learning the graph. It would be helpful to highlight the key method more clearly at the beginning of Section 3.\n5. The mathematical notation is dense and may be difficult to follow for readers unfamiliar with the specific formalism used. Simplifying or breaking down the notation could improve readability.\n6. Despite these points, I appreciate the idea of bridging constraint-based and score-based methods, which is a promising direction in causal discovery.", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n- Novelty: This is the first time I’ve seen a natural way to combine constraint-based and score-based methods, which is a significant contribution to the field.\n\nWeaknesses\n- The method description could be presented in a more straightforward way.\n- Learning Time: The paper does not report or compare the learning time, which is crucial for evaluating the method's efficiency, especially in large-scale settings.\n- Computational Cost: The method requires performing CI tests on all pairs of variables, which can be computationally expensive and may not scale well to large networks.\n- Code is not available", "was_revised": true, "final_justification": "I've read the authors' response and will keep my score.", "rebuttal": "We thank Reviewer JM3s for their positive evaluation and for recognizing the novelty of our work as a \"significant contribution to the field.\" We are grateful for the detailed and constructive comments on how to improve the paper's clarity and presentation. We agree with the reviewer's suggestions and will incorporate them into our revision.\nBelow, we address the specific weaknesses and questions raised.\n\n**On Computational Cost, Scalability, and Learning Time (Weaknesses 2, 3 & Question 1)**\n\nA1. The reviewer raises an important point about the computational cost and scalability of our method, particularly the need to handle all low-order CI statements.\n\n1. **Empirical Runtime Comparison**: The table below presents the wall-clock comparison of DAGPA against the baseline methods on the synthetic $d=10$ and $d=50$ experiment settings, which was the time these methods took to yield the results shown in the paper. We note that the wall-clock time of our method DAGPA shown below is an **overestimate** of the actual time it took DAGPA to reach the best resulting solution. This is because DAGPA uses Bayesian sampling (specifically Discrete Langevian Proposal (DLP) sampler) to continue exploring the DAG space. Rather than taking the last sampled DAG as the algorithmic output, we kept the algorithm running, kept track of all DAGs sampled at each step, evaluated them using the metric using only the input data (Appendix C3), and returned the best one. In almost all cases, the algorithm could have been terminated earlier to return DAGs of equal quality to those shown in the paper’s experimental results.\n\nGPU: 24GB A10\n|  | PC | kPC (k=1) | kPC (k=2) | GES | NOTEARS (Linear) | NOTEARS (Nonlinear) | DAGMA (Linear) | DAGMA (Nonlinear) | **DAGPA (Ours)** |\n|---|---|---|---|---|---|---|---|---|---|\n| $d=10$ | 0.60 | 0.15 s | 0.32 s | 1.80 s | 1.36 s | 10.98 s | 1.26 s | 172.92 s | 1.5 h (1000 steps) |\n| $d=50$ | 2.98 s | 2.60 s | 1.97 s | 921.14 s | 19.66 s | 44.23 s | 4.35 s | 801.98 s | 10.5 h (500 steps) |\n\n2. Number of CI Statements & Comparison to PC: Our method considers all 0th-order ($O(d^2)$) and 1st-order ($O(d^3)$) CI statements, for a total of $O(d^3)$ constraints. While the PC algorithm may be faster on sparse graphs by adaptively pruning its search, our goal is different. DAGPA is designed for robustness in low-sample regimes, where PC's reliance on binary CI test results, obtained via hard-thresholding p-values, can be challenging. Our method's strength lies in using \"soft\" CI information from all low-order statements to guide a differentiable search, which, as our experiments show, leads to superior accuracy when data is limited. The trade-off is a higher computational cost for greater robustness.\n\n3. Motivation vs. NOTEARS: We thank the reviewer for the opportunity to clarify our method’s motivations versus NOTEARS. While we do not claim to address the computational complexity of causal discovery problems as directly as NOTEARS did, we believe our method indeed has the following practical advantage in real-world settings.\n\n    -  **Differentiable Learning of k-Markov Equivalence Class**: Currently, there is no differentiable-learning-based approach similar to NOTEARS that can learn the k-Markov equivalence class in the low-sample regime, when we can only rely on low-degree CI tests since high-order ones are extremely unreliable. We argue that this is a significant shortcoming, since real-world tasks often have a limited number of samples, and the k-Markov equivalence class theoretically quantifies the ceiling of what can be learned under this scenario.\n\n    - **Interpretability via Conditional Independence (CI) Tests**: Our method leverages CI tests as its core mechanism, which can offer greater interpretability in some scenarios when compared to other score-based methods that rely on likelihoods. This could be particularly valuable in real-world applications where interpretability is critical.\n \n    - **Robustness in Low-Sample Regimes**: Unlike traditional score-based approaches, which often struggle with overfitting or underfitting when fitting functional models to scarce data, our method (DAGPA) incorporates \"soft\" CI information derived from p-values. This design choice enhances robustness in low-sample settings, as evidenced by its superior performance at $n = 100$ compared to baselines (Figure 3a). This highlights a practical advantage of DAGPA in data-limited scenarios.\n\n\n4. Code Availability: The code is available. Please refer to Appendix D.1 in the supplementary zip file, in which we provide a link to an Anonymous repository containing all the code for DAGPA. Upon acceptance, we will replace this repository with a Github link.\n\n\n\n**On Method Description and Clarity (Weakness 1 & Comments)**\n\nA2. We appreciate the reviewer's detailed feedback on how to make our method easier to understand. We will implement the following changes in our revision:\n\n- **Comment 1: Notation for Independence:** Thank you for pointing this out. For simplicity of presentation, we used the graphical independence notation $\\perp\\perp$, which is used in works such as Pearl, (2009); Koller & Friedman (2009), to also express independence observed from the dataset $\\mathcal{D}$. In order to improve readability, we will add a note on the equivalence between $P(x \\mid y, z) = P(x \\mid z)$ and the notation in the paper. \n\n- **Comment 2: Theorem 3.2 vs Definition 3.1:** We thank the reviewer for highlighting this important aspect of the presentation. Theorem 3.2 serves to establish the correctness of our newly proposed first-order logical characterization of low-order d-separation (Definition 3.1), which is a foundational contribution of this work. The proof demonstrates that the output values derived from the logical formulae in Definition 3.1 precisely match those obtained via traditional procedural combinatorial algorithms for d-separation. This equivalence is important, as it enables the continuous, differentiable relaxation introduced in Section 3.2—a capability not readily achievable with traditional, non-differentiable combinatorial methods.\n\n\n- **Comment 3: Redundancy of Lemma 3.4:** We agree that it may interrupt the flow of the text, so we will move it to Appendix in the revised paper. Although the correctness of the reachability FOL formulae in Definition 3.3 may seem intuitively straightforward, Lemma 3.4 is still a necessary building block for a mathematically rigorous argument, since we believe it formally justifies the adoption of the reachability predicate $R_{\\mathbf{A}}$​ in our d-separation formulae.\n\n- **Comment 4 & 5: Highlighting the Method and Simplifying Notation:** We agree that the flow of Section 3 and a description and break-down of the key mathematical notation can be further improved. We will improve the \"roadmap\" paragraph at the beginning of Section 3 to make the method description more intuitive, and also review the entire section to add more intuitive explanations for our notation, such as a concrete running example, and break down complex formulae to improve readability.\n\nWe are confident that these revisions will improve the clarity and accessibility of our paper. We thank the reviewer again for their supportive assessment and valuable suggestions."}, {"review_id": "zcRA2wdF9z", "reviewer": "Reviewer_523A", "summary": "In this paper, the authors introduce DAGPA (DAG Percolation Apartness), a novel causal discovery method formulated as an optimization problem, which directly employs results from conditional independence tests. In detail, DAGPA computes the $p$-value of marginal independence tests between any two distinct variables $(x, y)$ and of conditional independence tests between any three distinct variables $(x, y, z)$. Then, the authors train the model by relating conditional independence statements with $d$-separation conditions. They inscribe this matching in a differentiable optimization problem by (i.) formulating $d$-separation as a first-order-logic formula and (ii.) applying a soft logic framework to optimize the DAG edges. The authors test empirically their work against differentiable score-based baselines (NOTEARS, DAGMA) and constraint-based baselines (PC).", "questions": "- In Appendix, the method underperforms in terms of F1 and SHD. Can the authors justify this behaviour? While I agree that MCC is an interesting measure, the paper should also discuss more in detail its relation with more widely adopted metrics in the main body.\n- It seems that to compute the loss, you need to precompute all possible conditional independence tests. Is this correct? How would this compare with the number of tests performed by PC? How does DAGPA compare to PC and score-based methods in terms of computational time?\n- Differentiable score-based methods, such as NOTEARS and DAGMA, are proven to work remarkably well on non-normalized data, while struggling to recover the DAG after rescaling data (\"Beware of the Simulated DAG! Causal Discovery Benchmarks May Be Easy To Game.\", Reisach et al., 2021). How is the pre-processing handled in your experimental pipeline?\n- Can this approach be integrated as a sort of \"regularization\" technique to existing score-based differentiable methods? For instance, in the linear case, I would have imagined to also minimize the reconstruction MSE loss $\\|\\|\\mathbf{W}^\\top\\mathbf{X}-\\mathbf{X}\\|\\|_2^2$ adopted by NOTEARS et similia. Is there some fundamental reason to avoid doing so?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths.\n\n- The method is designed to improve causal discovery performance in the small sample size regime, and experimental results hint at better results in terms of Matthews Correlation Coefficient (MCC) compared to baselines in this particular scenario.\n- The overall method is novel, and it is interesting the idea of mixing conditional independence-based methods with optimization approaches.\n- While the use of first-order-logic formulas lead to a quite convoluted presentation, the authors' effort in describing step-by-step most pivotal passages is appreciable.\n\nWeaknesses.\n\n- The method requires computing conditional independence among all possible pairs and triples of distinct variables to define the loss function of the model. This might significantly affect performance in terms of computational time.\n- Widely accepted metrics for causal discovery, such as F1 and SHD, are relegated to the Appendix where (i.) the method underperforms compared to the baselines (ii.) I would have expected the score-based baselines to have better results. However, this second point might not be a problem depending on the data-generating process adopted by the authors — please refer to the question.\n- Not really a weakness, but please notice two typos \"pioneerd\" (line 96) and \"Hammindg\" (line 295).", "was_revised": true, "final_justification": "The authors addressed my comments concerning the evaluation metric, the pre-processing pipeline, and the integration with linear causal structure learning algorithms. I am concerned by the computational time required by the method. However, the paper opens up an interesting novel direction and might be worth accepting.", "rebuttal": "We thank Reviewer 523A for their detailed feedback and constructive criticism. We appreciate that the reviewer found our method \"novel\" and the core idea \"interesting.\" The questions raised are important, and they will help us significantly improve the paper's discussion on scalability, evaluation, and its relationship with existing methods. We will also correct the typos noted.\n\nBelow, we address the specific weaknesses and questions.\n\n\n> W2 & Q1. The method underperforms in terms of F1 and SHD. Can the authors justify this behaviour? The paper should also discuss more in detail its relation with more widely adopted metrics in the main body.\n\n\nA1. We appreciate the reviewer's thoughtful question. To address this concern, we direct the reviewer to **Appendix E.1**, where we provide a detailed analysis of the relationship between CI-MCC and traditional metrics such as F1 and SHD. This section explicitly motivates our choice of CI-MCC as the primary evaluation metric.  \n\nThe core objective of our method is to learn a DAG whose **low-order conditional independence (CI) statements** align with those observed in the data. CI-MCC directly measures this alignment, whereas F1 and SHD assess structural similarity, which is not always equivalent to CI accuracy. As discussed in **Appendix E.1**, structurally distinct DAGs can imply identical sets of low-order CI statements. Our empirical analysis demonstrates numerous cases where DAGs sampled by our method achieve high CI-MCC scores (indicating few CI errors) but exhibit widely varying F1/SHD scores. This underscores that while our method successfully achieves its intended objective, it does not necessarily optimize for structural metrics like SHD or F1.  \n\nThis behavior is inherent to causal discovery, as different graphs within the same *k*-Markov equivalence class can yield identical CI statements despite structural differences. **Furthermore, recent work by Petersen (2025) highlights limitations of F1 and SHD, showing they often fail to distinguish model performance with sufficient statistical significance.** Given these challenges, we prioritized CI-MCC as the main evaluation metric while retaining F1/SHD results in the Appendix for completeness.\n\nWe acknowledge that our discussion of these metrics was confined to the Appendix in the initial submission. Thanks to the reviewer’s feedback, we will incorporate a concise paragraph in the main experimental section (Section 5) to explicitly highlight this distinction and direct readers to the detailed analysis in **Appendix E.1**.\n\nPetersen, A. H. (2025). *Are You Doing Better Than Random Guessing? A Call for Using Negative Controls When Evaluating Causal Discovery Algorithms*. UAI 2025.  \n\n\n> W1 & Q2: The method requires computing conditional independence among all possible pairs and triples of distinct variables to define the loss function of the model. This might significantly affect performance in terms of computational time. How does DAGPA compare to PC and score-based methods in terms of computational time?\n\nA2. The reviewer raises an important point about scalability, and it is correct that the DAGPA implementation requires pre-computing a bounded number of conditional independence tests $O(d^3)$. While this can still be computationally intensive, using parallel computation we observe in practice that our method has manageable runtime in our experiment setting. Below, we present the wallclock empirical comparison of our method against the baselines.\n\nWe note that the wall-clock time of our method DAGPA shown below is an **overestimate** of the actual time it took DAGPA to reach the best resulting solution. This is because DAGPA uses Bayesian sampling (specifically Discrete Langevian Proposal (DLP) sampler) to continue exploring the DAG space. Rather than taking the last sampled DAG as the algorithmic output, we kept the algorithm running, kept track of all DAGs sampled at each step, evaluated them using the metric using only the input data (Appendix C3), and returned the best one. In almost all cases, the algorithm could have been terminated earlier to return DAGs of equal quality to those shown in the paper’s experimental results.\n\n\nGPU: 24GB A10\n|  | PC | kPC (k=1) | kPC (k=2) | GES | NOTEARS (Linear) | NOTEARS (Nonlinear) | DAGMA (Linear) | DAGMA (Nonlinear) | **DAGPA (Ours)** |\n|---|---|---|---|---|---|---|---|---|---|\n| $d=10$ | 0.60 | 0.15 s | 0.32 s | 1.80 s | 1.36 s | 10.98 s | 1.26 s | 172.92 s | 1.5 h (1000 steps) |\n| $d=50$ | 2.98 s | 2.60 s | 1.97 s | 921.14 s | 19.66 s | 44.23 s | 4.35 s | 801.98 s | 10.5 h (500 steps) |\n\nIn addition, we would like to highlight several practical techniques that can greatly enhance the computation efficiency in real-world deployment:\n\n- First, our differentiable d-separation framework is amenable to efficient matrix-vector operators that can benefit from **vectorizations**. This is what we did for the DAGPA implementation, where we leverage PyTorch tensor library to implement all forward passes (i.e. computation of differentiable d-separation scores) as tensor operations without explicit for-loops. This allows us to utilize CUDA GPUs to significantly improve the computation speed, similar to how the nonlinear version of NOTEARs may leverage GPUs.\n\n\n- Second, the all-pairs p-values computation can be massively parallelized and then cached and reused. This is particularly useful in the case where one wants to obtain multiple output causal graphs for the same input dataset to enhance solution robustness and exploration of the Markov equivalence class. In this case, the computation time for the all-pairs p-values step can be amortized as it can be reused in subsequent runs.\n\n\n- For future work, another possibility is to use minibatching. That is, for very large graphs where computing all CI statements is infeasible, our framework can be extended so that at each step, we could sample a random subset of CI statements to estimate the gradients, trading a small amount of gradient accuracy for a significant reduction in per-step cost. \n\n\nIn the paper revision, we will add the empirical runtime analysis to the Appendix and in the main paper discussion section, emphasize which directions can yield scalability improvement worthy of future work to explore.\n\n\n> W2 & Q3. How is the pre-processing handled in your experimental pipeline?\n\n\nA3. Thank you for this interesting question regarding the work of Reisach et al. (2021). This allows us to highlight a key strength of our framework.\n\n1. **Experimental Pipeline:** For our synthetic experiments, we used binary data, where normalization/scaling is not needed. For the real-world Sachs dataset, we used the standard discretized version preprocessed using the Hartemink discretization method, which converts the original continuous protein concentration measurements into 3-level categorical variables (low, average, high) while preserving the underlying dependence structure.\n\n2. **Robustness of Our Method:** Our framework should be inherently more robust to pre-processing issues than other score-based methods. The issue highlighted by Reisach et al. is that methods like NOTEARS can \"game\" the benchmark by exploiting variance differences as a proxy for causality, a trick that fails when data is normalized. Our approach, however, relies on CI tests. Standard CI testers (like the Chi-squared test for binary data or Fisher's Z-test for continuous data) are robust to these kinds of scaling transformations. Therefore, DAGPA does not rely on such data artifacts and is expected to be more reliable across different pre-processing pipelines.\n\nWe will add a paragraph to the experimental setup section clarifying our pre-processing steps and explaining why our constraint-based approach is less susceptible to these known pitfalls of score-based methods.\n\n\n> Q4. Can this approach be integrated as a sort of \"regularization\" technique to existing score-based differentiable methods?\n\n\nA4. Yes, absolutely. Our formulation of the differentiable d-separation framework should be general enough to support such a potential hybrid method, combining the strength of both constraint-based scores and likelihood score (such as NOTEARS’ reconstruction MSE loss). Below we provide a more detailed explanation of how such a hybrid method might be implemented.\n\n\nThe main challenge to address is defining how a weighted adjacency matrix $\\mathbf{W} \\in [0, 1]^{d \\times d}$ with edge weights in the range $[0, 1]$ can be derived from the model parameter $\\theta$ in a way that is compatible with both our differentiable d-separation framework and the score-based method. Particularly to our method, we require the edge weights to fall in $[0, 1]$ as it needs to be valid probability measures.\n\nFor instance, to convert the model parameter to our framework’s probabilistic edge weights, it should be possible to apply a sigmoid activation $\\sigma(\\cdot)$ for the NOTEARS linear model’s parameter $\\theta$. For the nonlinear version of NOTEARS though, its graph weighted adjacency matrix, obtained by $\\text{sqrt}({\\theta^{(1)}\\_j}^T \\theta^{(1)}\\_j)$ ($\\theta^{(1)}$ the first MLP layer weight) consists of unnormalized, nonnegative values. Thus, one needs to find a proper way to normalize such a weight into range $[0, 1]$ to satisfy the probabilistic semantics. The rest is to simply add the prediction loss, such as the MSE loss $\\mathcal{L}_{\\text{MSE}}(\\theta, X) = \\frac{1}{2n} \\| X - XW \\|^2_2$ for the linear NOTEARS model, as an additional task loss and apply the same multi-task optimization routine of DAGPA.\n\n\n\n\nIt is our belief that our work opens up the opportunity of an interesting research direction, and merits further study on possible ways to integrate with various existing score-based methods."}]

Author rebuttal:

Meta-review: The paper introduces a new route for causal discovery that optimizes conditional-independence (CI) constraints with gradients by making d-separation differentiable. The authors (1) give formal first-order logic characterizations of 0th- and 1st-order d-separation and prove correctness, (2) relax these logical formulas into continuously differentiable quantities using LogLTN soft logic, yielding lower-bound estimators of expected d-separation over random DAGs, and (3) instantiate the idea in DAGPA, which learns a weighted adjacency by minimizing multi-task losses that align soft d-separation/connection with CI evidence while enforcing acyclicity. Empirically, the method is robust in low-sample regimes and outperforms popular baselines.

Strengths：
All reviewers agree that the paper is interesting and could potentially lead to new lines of work. 
Reviewer 523A further emphasizes that the better results obtained in the paper compared to baselines. Reviewer web4 appreciates that the paper presents the results in a coherent and step-by-step manner.  


Weaknesses：
1）Relies on CI p-values as heuristics. The training signal uses p-values from Fisher-z tests (even on nonlinear data like Sachs), which can inject bias/mismatch between the statistical test and the ground-truth mechanism.
2）Scope limited to low-order conditioning (|Z| ≤ 1). This is practical and well motivated, but it narrows identification power. Scaling the relaxation to higher-order CI would require exponentially many additional constraints, leading to severe computational challenges. 

The most important reasons for the decision：
1）A clear, principled bridge between constraint-based and continuous optimization methods via a novel, well-justified relaxation of d-separation;
2）A practical, consistent learning framework that operationalizes the idea with acyclicity and multi-task CI losses, validated on synthetic and real data.

Summarization of the discussion and changes during the rebuttal period: 
Some reviewers raised concerns about the computational complexity and empirical runtimes of the method, which the authors addressed with additional experiments and analysis in the rebuttal. 

Overall comments: 

After the rebuttal phase, the paper received one strong accept, one accept, and two borderline accepts, resulting in unanimous acceptance. The paper introduces a new way to combine constraint-based and score-based causal discovery methods. All reviewers agree that the proposed methods are interesting. Overall, the quality of the work is high and solid, and it is likely to be of interest to the community. Therefore, I suggest to accept this paper as a poster.
Decision: Accept (poster)