Differentiable Cyclic Causal Discovery Under Unmeasured Confounders
Authors: Muralikrishnna Guruswamy Sethuraman
Faramarz Fekri
Venue: NeurIPS 2025 spotlight ()
Keywords: causal discovery
confounders
feedback loops
generative modeling
Primary area: probabilistic_methods
TL;DR: We propose DCCD-CONF, a scalable, theoretically grounded framework for learning nonlinear cyclic causal graphs with unmeasured confounders using interventional data.

Abstract: Understanding causal relationships between variables is fundamental across scientific disciplines. Most causal discovery algorithms rely on two key assumptions: (i) all variables are observed, and (ii) the underlying causal graph is acyclic. While these assumptions simplify theoretical analysis, they are often violated in real-world systems, such as biological networks. Existing methods that account for confounders either assume linearity or struggle with scalability. To address these limitations, we propose DCCD-CONF, a novel framework for differentiable learning of nonlinear cyclic causal graphs in the presence of unmeasured confounders using interventional data. Our approach alternates between optimizing the graph structure and estimating the confounder distribution by maximizing the log-likelihood of the data. Through experiments on synthetic data and real-world gene perturbation datasets, we show that DCCD-CONF outperforms state-of-the-art methods in both causal graph recovery and confounder identification. Additionally, we provide consistency guarantees for our framework, reinforcing its theoretical soundness.

Differentiable Cyclic Causal Discovery Under Unmeasured Confounders
Muralikrishnna G. Sethuraman School of Electrical & Computer Engineering Georgia Institute of Technology mural i kgg@ggtechhhc edu
Faramarz Fekri
School of Electrical & Computer Engineering
Georgia Institute of Technology faramarz. fekr @ece gatech. edu
Abstract
Understanding causal relationships between variables is fundamental across sci- entific disciplines. Most causal discovery algorithms rely on two key assump- tions: (i) all variables observed, and (ii) the underlying causal graph is acyelic. While these assumptions simplify theoretical analysis, they are often violated in real-world systems, such as biological networks. Existing methods that account for confounders either assume linearity or struggle with scalability. To address these limitations, we propose DCCD-CONF, a novel framework for differentiable learning of nonlinear cyclic causal graphs in the presence of unmeasured con- founders using interventional data. Our approach alternates between optimizing the graph structure and estimating the confounder distribution by maximizing the log- likelihood of the data. Through experiments on synthetic data and real-world gene perturbation datasets, we show that DCCD-CONF outperforms state-of-the-an methods in both causal graph recovery and confounder identification. Additionally, provide consistency guarantees for our framework, reinforcing its theoretical soundness.
1 Introduction
Modeling cause-effeett relationships between variables is fundamental problem in science [1] 31. as it enables the prediction of system's 's behavior under previously unseen perturbations. These relationships typically represented using directed graphs (DGs), where nodes correspond to variables, and directed edges capture causal dependencies. Consequently, causal discovery reduces to learning the structure of these graphs.
Existing causal discovery algorithms can be broadly classified into three categories: (i) constraint- based methods, (ii) score-based methods, and (iii) hybrid methods. Constraint-based methods, such as the PC algorithm constraints observed in the data. However, since the number of conditional independence tests grows exponentially with the number of nodes, these methods often struggle with scalability. Score-based methods, such as the GES algorithm 77] (8), learn graph structures by maximizing penalized score function, such as the Bayesian Information Criterion (BIC), over the space of graphs. Given the vast search space, these methods often employ greedy strategies reduce computational complexity. significant breakthrough came with Zheng et al. 9 who introduced a continuous constraint formulation to restrict the search space acyclic graphs, inspiring several extensions (10) [ [2]
4 5 (), search for causal graphs that best satisfy the independence
13 14 [5] that frame causal discovery as a continuous optimization problem under various model assumptions. Hybrid methods 16
17 18 integrate aspects of both constraint-based and score-based approaches, leveraging independence constraints while optimizing score function.
Most causal structure learning methods assume (i) directed acyclic graph (DAG) with no directed cycles and (ii) complete observability, meaning no unmeasured confounders.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

While these ussumptions simplify the search space, they are often unrealistic, as real-world exhibit feedback loops and hidden confounders [19]. Enforcing these constraints can also in- crease computational complexity, particularly in ensuring acyclicity, which often requires solving challenging com- binatorial constrained optimization problems. These limitations hinder the practical applicability of existing methods in settings where such violations are unavoid- able.
in
X2
Xs
(a)
(b)
Figure 1: (a) Example of directed
mixed graph g, where the bidirectional edges represent hidden confounders, with Jij indicating their correspond- ing strengths; (b) Mutilated graph, do (G). resulting from the interven- tional experiment I& {X3}. where all incoming edges (including bidirectional edges) to X3 are removed.
Several approaches have been developed to address the challenge of feedback loops within causal graphs. Early work by Richardson [20 extended constraint-basee ap- proaches for Directed Acyelic Graphs (DAGs) to accom- modate directed cycles. Another key contribution came from Lacerda et al. [21] who generalized Independent Component Analysis (ICA)-based causal discovery to han- dle linear non-Gaussian cyclic graphs. More recently, a growing body of research has focused
ccore-based methods for learning cyclic causal graphs 24 25I. Additionally, some approaches leverage interventional data to improve structure recovery in cyclic systems. For instance, Hyttinen et al. [26 and Huetter and Rigollet (22] introduced frameworks that explicitly incorporate interventions to refine eyclic graph estimation. Sethuraman al. 27 further advanced this line of research by introducing differentiable framework for learning nonlinear cyclic graphs. Unlike differentiable DAG learners that enforce acyclicity through aug- solvers, their approach sidesteps these constraints by directly modeling the data likelihood, enabling more efficient and flexible learning of cyclic causal structures. However, their method assumes the absence of unmeasured confounders, which limits its applicability in real-world settings where hidden confounders are often present.
22
23
mented Lagrangian-baan
Causal discovery in the presence of latent confounders has seen limited development, with most existing approaches grounded constraint-basee methodologies. Extensions of the PC algorithm, such the Fast Causal Inference (FCI) algorithm (28). construct Partial Ancestral Graph (PAG) to represent the equivalence class of DAGs in the presence of unmeasured confounders, and can accommodate nonlinear dependencies depending
the chosen conditional independence tests. However, standard FCI does not incorporate interventional data, prompting extensions such as JCI-FCI 29 and related approaches [30] that combine observational and interventional settings. Jaber et al. [3 further advanced this line of work by allowing for unknown interventional targets. Additionally, Forré and Mooij [32 introduced a s-eparationn, a generalization of d-separation, enabling constraint-basee causal discovery in the presence of both cycles and latent confounders. Suzuki and Yang (33] introduce LiNGAM-MMI, generalization of the ICA-based LiNGAM (34] that quantifies and mitigates confounding via a KL divergence minimization. Integer-programen programming based formulations have also been proposed for causal discovery under latent confounding, such as [35 36
A few recent approaches, such as Bhattacharya et al. (37, have explored continuous optimization frameworks using differentiable constraints, though these methods are currently limited to linear settings. In parallel, several works exist on causal inference under latent confounding-nent Abadie et al. [38], Chernozhukov et al. (39 propose doubly robust estimators that integrate outcome modeling, weighting, and fross-fitting for reliable effect estimation. Overall, unified framework capable of handling nonlinearity, cycles, latent confounders, and interventions remains largely absent.
notably
Contributions. In this work, we tackle three key challenges in causal discovery: directed cycles, nonlinearity, and unmeasured confounders. Our main contributions are:
We introduce DCCD-CONE, a novel differentiable causal discovery framework for learning nonlinear eyclic relationships under Gaussian exogenous noise, with confounders modeled as correlations the noise term.
We show that exact maximization of the proposed score function results in identification of the interventional cquivalence class of the ground truth graph.
We conduct extensive evaluations, comparing DCCD-CONF with sate-of-the-an causal discovery methods on both synthetic and real-world datasets.

---

Organization. The paper is structured as follows: Section 2 introduces the problem setup. In Section 3 we present DCCD-CONE, our differentiable framework for nonlinear cyclic causal discovery with unmeasured confounders. We then evaluate its effectiveness on synthetic and real- world datasets in Section 4 Finally, Section concludes the paper.
2 Problem Setup
2.1 Structural Equations for Cyclic Causal Graphs
Let 9 = EV,E, B) represent a possibly cyclic directed mixed graph (DMG) that encodes the causal dependencies between the variables in the vertex set V [d]. where [d] {L..... d]. & denotes the set of directional edges of the form i + in G. and B denotes the set of bidirectional edges of the form i ++ j in G. Each node associated with a random variable X; with the directed edge + j E representing causal relation between X, and Xj, and the bidirectional edge i ++ j E B indicates the presence of hidden confounder between X; and X,. Following the framework proposed by Bollen 40 and Pearl [1], we use structural equations model (SEM) to algebraically describe the system:
X; = =F(XX =F((Xu Zz), i == ,,...dd
(i)(t),
where pag(1) = fje [d] j + i E E} represents the parent set of X, in G. and Xps denotes the components of (XI,..., Xa) indexed by the parent set pag(i). We exclude self-loops (edges of the form X; - X,) from G. as their presence can lead to identifiability challenges [42). The function Fú, referred to as causal mechanism, encodes the functional relationship between X; and its parents
X (a)) and the exogenous noise variable Z,.
The collection of exogenous noise variables Z = -Z...... Za) account for the stochastic nature as well as the confounding observed in the system. We make the assumption that the exogenous noise vector follows a Gaussian distributions: Z ~ NoO Σz). Notably, if (Ez)i * 0. then variables X; and X3 are confounded, i.e., + j E B. other words, confounding correlations in the exogenous noise variables. Intuitively, if X; and X, share hidden cause, their unexplained variation (the part not accounted for by their observed parents) will tend move together. By allowing the noise terms Z; and Z, to be correlated, this shared influence can be effectively captured. This formulation generalizes prior work by allowing cycles, extending both nonlinear cyclic models that assume independent noise terms [27. and acyclic models without confounders 43].
modeled through
By collecting all the causal mechanisms into the joint function F (Fi,..., Fa). we can then combine II over i = 1,.... d to obtain the equation
X = F(X,))
(2)
We will use ( to represent the causal system due to its simplicity for subsequent discussion. The observed data represents a snapshot of dynamical process where the recursive equations in ( define the system's 's state at its equilibrium. Thus, in our experiments we assume that the system has reached the equilibrium state. For given random draw of Z, the value of X is defined as the solution to (2). To that end, we assume that 2 admits a unique fixed point for any given Z. We refer to the map f. X + Z as the forward map, and f. Z + X the reverse map. In Section 3.1] we show that the chosen parametric family of functions indeed guarantees the existence of a unique fixed point. Under these restrictions, the probability density of X well defined and given by
Pg (X) pz (X)) det (Jr, (X)),
where Jr, (X) denotes the Jacobian matrix of the function f, at X.
2.2 Interventions
In our work, we consider surgical interventions [4], also known as hard interventions, where all the incoming edges to the intervened nodes are removed from G. Given a set of intervened upon nodes (also known as interventional targets), denoted as / C V. the structural equations in (T) are modified as follows
X; [Ca (X (1) Ze),
if X, I
if $ I,
3

---

where C2 is random variable sampled from a known distribution, i.e.. C2 ~ Pr (C.). We denote (G) be the mutilated graph under the intervention (see Figure D.. Note that X; no longer confounded if is intervened on.
We consider family of K interventional experiments I = {Tx (TR) where Ix represents the interventional targets for the k-th experiment. Let Us E {0, 1 denote diagonal matrix with (Us)a 1 if i & JR. and (UR)a: = 0 if i E Iu. Similar to the observational setting, 4 can be vectorized to obtain the following form
(5) where C = (CI....,, Ca) is a vector with C. ~ pr (C.) if E In. and C, = 0 otherwise. For the interventional targets 44 I. let (f denote the forward map. Similar to the observational setting, we make the following assumption on the set of interventions.
X = U,F(X, +)
Assumption 1 (Interventional stability). Let I [ee] be family of interventional targets. For each 1 E I, the structural equations in ( admits a unique fixed point given the exogenous noise vector Z.
Thus, the probability distribution of X for the interventional targets 4a is given by
(ao(Ia)(s)( = pr(C)pz ( (E (X) (x)) det (J (X)))
(6)
where U4 = ( T} denotes the index of purely observed nodes, and pz the marginal distribution of the combined vector Z, restricted to the components indexed by Uy.
[r. (x) (x))mm) is
Given a family of interventions I. our goal is to learn the structure of the DMG by maximizing the log-likelihood of the data, in addition to identifying the variables that are being confounded by the unmeasured confounders Z. The next section presents our approach addressing this problem.
3 DCCD-CONE: Differentiable Cyclic Causal Discovery with Confounders
In this section, we present our framework for differentiable learning of cyclic causal structures in the presence of unmeasured confounders. We start by modeling the causal mechanisms, then define the score function used for learning, followed by theorem that validates its correctness. Finally, we outline the algorithm for estimating the model parameters.
3.1 Modeling Causal Mechanism
We model the structural equations in 2 using implicir flows [44], which define an invertible mapping
between æ and z by solving the root of function G( 2) 0, where G 1g - Rd. Specifically,
we take G(æ, x) = x F F(x,= z). General implicit mappings, however, do not guarantee invertibility
or permit efficient computation of the log-determinant required for evaluating ( To balance
expressiveness with tractability, we adopt the structured form proposed by Lu et al. for the causal mechanism:
F F(x, 2) = -g (x) + g=(2) z,
where ga and g: are restricted to be contractive functions. A function g: Rd - Rd is contractive there exists a constant I. < such that [[g(>)) - g(y) S LIIa - yll for all æ, y E Rdd This contractiveness ensures that the associated implicit map is uniquely solvable and invertible (see Theorem in [44)). In other words, contractivity ensures that the process defined by the SEM converges to an equilibrium state.
Under this formulation, the forward map takes the form (x (id + +g.) --gg (id + (x), where id denotes the identity map. Given æ (or 2), the corresponding value of z (or x) can be computed via a moot-finding procedure, i.e., z RootFind(2: F(x, +)), specifically, we employ a quasi-Newtor method (i.e., Broyden's method 45 to find the root. To capture more complex nonlinear interactions between the observed variables X and latent confounders Z, multiple such implicit blocks can be stacked. This is true since f, is highly nonlinear and by suitably parameterizing g.. and g. any nonlinear interaction between æ: and can be modeled. For simplicity, we focus on single implicit flow block for subsequent discussion.

---

We parameterize the functions 8N and g: using neural networks. The adjacency matrix of the causal graph encoded as a binary matrix M E (0,1) representing the presence of directed edges and serving as mask on the inputs to ga- 'hhe diagonal entries of M are explicitly enforced to be zero to prevent self-loops. Similarly, the identity matrix is used to mask the inputs g:- Consequently, the causal mechanism is defined as:
[X 2)), = [- NN(M + X O₂)) NN 02) + Z:];
(8)
where NN(. 9) denotes a fully connected neural network parameterized by 0, O denotes the Hadamard product, and M.;- is the i-th column of M The contractivity of g. and ga can be enforced by rescaling their weights using spectral normalization [46]. Moreover, the contractive nature of the causal mechanism facilitates efficient computation of the score function used for learning causal graphs, as discussed in Section 3.2
While the contractivity assumption may seem restrictive, it ensures stability and well possedessss in the presence of directed cycles. If the causal graph is known to be acyclic, this assumption can be
relaxed (see Appendix C.I.
3.2 Score function
Given family of interventions I = kE(R) we would like to learn the parameters of the structural equation model, i.e., causal graph structure, causal mechanism, and confounder distribution. To that end, similar to prior work [27 47 [5] in this domain we employ regularized log- aigliikelhoodd of the observed nodes as the score function to be maximized. That is,
K
ST( = sup
log (ogga (X) A|9| X~ppik
E
(9)
0,5g
where p(8) is the data generating distribution for the k-th interventional experiment 7g. Σz is the parameter (covariance matrix) governing the confounder distribution pz. 0 = (0z, 0.) is the combined causal mechanism parameters, and [ denotes a sparsity enforcing regularizer on the edges of G, and (a)( (X) is given by (6).
We now present the main theoretical result of this paper. The following theorem establishes that, under appropriate assumptions, the graph & estimated by maximizing (9) belongs to the same general directed Markov equivalence class (introduced by (42) as the ground truth graph g* for each interventional setting I4 E I, denoted as 9 = G*, see Appendix A. Due to space constraints we provide the proof sketch below, see Appendix AA.3 complete proof of Theorem
Theorem 2. Let I {I be a family of interventional targets, let G* denote the ground truth directed mixed graph, let pik) denote the data generating distribution for It. and & := arg maxç S()). Then, under the Assumptions A.13 A.14,
A.5 and for suitably chosen X O, we have that 9 =1 G" That is, 9 is I-Markov equivalent to G*.
Theorem rests on three key assumptions. Assumption A. ensures that the data-generating distribu- tion lies within the model class, while Assumption A.14 guarantees that every statistical independence in the data corresponds to a a-separation in the ground-truth graph. Finally, Assumption | prevents the score function from diverging infinity.
Proof (Sketch). Building the characterization of general directed Markov equivalence class by Bongers et al. [42, extended to the interventional setting, we show that any graph outside this equivalence class has strictly lower score than the ground truth graph G". This follows from the fact that certain independencies present in the data are not captured by graphs outside the equivalence class. Combined with the expressiveness of the model class, this prevents such graphs from fitting the data properly.
If the intervention set consists of all single-node interventions, I = (1x) with I = (k), Hyttinen et al. (26] showed that the ground truth DMG can be uniquely recovered in the linear setting. Moreover, in the absence of cycles and confounders, this result extends to the nonlinear case, as demonstrated by Brouillard et al. [15). However, determining the necessary conditions on
5

---

interventional targets for perfect recovery in general DMGs with cycles and confounders remains an open problem. Nonetheless, practice, find that observational distribution combination with single-node interventions across all nodes lead to perfect recovery of the ground truth, even in the nonlinear case, as shown in Section
3.3 Updating model parameters
In practice, we use gradient based stochastic optimization to maximize (9) For this purpose, following Sethuraman al. [27] and Brouillard et [15), the entries of adjacency matrix Ma are modeled as Bemoulli random variable with parameters bij- grouped into the matrix o (B). We denote M ~ o (BB) to indicate that Mij Bern(b;y) for i,j E [d]. In this formulation, the sparsity regularizer is |M which is computationally intractable and thus we use the f1-nomm, ||M]] proxy. Consequently, the score function in (
is replaced by the following relaxation:
S7
sup E
ΣΣ log Pau(I, (e)
- A -
(10)
0, 0,E₂ M
where we replace the expectation with respect to data distribution in (9) with sum over the finite samples, (1.A) represents the i-th data sample in the k-the interventional setting. We note that, since Pz (0, Σz), the covariance of the exogenous confounder vector, Σz, implicitly embedded within Pao(Tr)(()) (ix()) in the score function.
The optimization of the score function carried two steps. First, we optimize S(B) with respect to the neural network parameters 0 and the graph structure parameters B. Next, we optimize S(B) with respect to the parameters of the exogenous noise distribution, Σz. However, maximizing (B) presents two main challenges: (i) computing log (X) is computationally expensive due to the presence of det((J g ))}, which requires (dd) gradient calls, and (ii) updating Σz via stochastic gradients could lead to stability issues as Σz may loose its positive definiteness.
We now describe how these challenges are addressed, along with the specific procedures for updating the individual model parameters.
3.3.1 Computing log determinant of the Jacobian
is a significant challenge in maximizing the score function S(B). To address this, we utilize the unbiased estimator of the log-determinant of the Jacobian introduced by Behrmann
As discussed earlier, computing log [J J (X) |
46), which is based on the power series expansion of
log(1 + z). Since r!!x) (x) = (id U&g:) (id (x)
log det (J (J)) = log det (I + Jung. (x)) log det (I + Ju.g. (2))
1-1)m+1
Tr Tr.Tu (X) TT(JU.. (2))
(11)
m
where E Rdxdd denotes the identity matrix, U.g. represents the Jacobian matrix raised to the m-th power, and Tr denotes the trace of matrix. The series in (T]) guaranteed to converge if the causal functions gx and ga are contractive [48].
In practice, the power series is truncated finite number terms, which may introduce bias into the estimator. To mitigate this issue, we follow the stochastic approach of Chen et al. [49]. Specifically, we sample random cut-off point n ~ Ps((n) for truncating the power series and weight the i-term in the finite series by the inverse probability of the series not ending at i. This yields the following unbiased estimator
(-1)m+1
Tr (X) Tr (Z)}
E log det (J (J (J)))()))] a>pp(())
Σ
(12)
m
pa(e 2 m)
m=1
The gradient calls can be reduced even further using the Hutchinson trace estimator Appendix B for more details.
50 see

---

3.3.2 Updating neural network and graph parameters.
In the first step of the parameter update, keeping Σz fixed, the parameters of the neural network O and the graph structure B are updated using the backpropagation algorithm with stochastic gradients. The gradient of the score function (B) with respect to B is computed using the Straight- Straight-Through Gumbel estimator. This involves using Bernoulli samples in the forward pass while computing score, and using samples from humbel-Softmas distribution in the backward pass to compute the gradient, which can be differentiated using the reparameterization trick [51]
3.3.3 Updating the confounder-nois distribution parameters
In second parameter update step, we fix the value of 0 and B and focus on the confounder-nois noise distribution parameter Σz. First, consider the case where no interventions are applied, i.e, I& Note that the dependence of ((B on Σz arises solely from pz, which is embedded within (C) (X). Therefore, we can thus ignore the remaining terms in (BB and focus exclusively on Pz. Let {2x(0) ( =i=1 denote the observational data. From the forward map, we have 2(1) f, (x:)).). Given that pz (0, Σz), the relevant parts of S(B) with respect Σg, denoted as [(Ix), are expressed as:
N
1 2
(x(1) Σg (0)
((1))
sup Σz
log
Σz|
(13)
==1
Simplifying 13 yields a more convenient form:
E(Ik)
sup -T((SEZ)) log Σz/,
(14)
Σz
where S
1 ΣE' 1=1 | (z))T is the sample covariance of Z.
Maximizing 14 directly using backpropagation and stochastic gradients results in stability issues as Σg may lose its positive definiteness. However, Friedman et al. [52] demonstrated that the sparsity-regularized version of (14) is a concave optimization problem in Σz that can be efficiently solved by optimizing the columns of Σz individually. This is achieved by formulating the column recovery as a lasso regression problem. We adopt this strategy while updating the Σz during the maximization of (B).
Let W Σz be the estimate of the covariance matrix. We reorder W such that the column and row being updated can be placed at the end, resulting in the following partition
W 11
W12
S 1
812
W
(15)
W12
W22
12
8:2
Then, shown by Friedman al. [2] regression problem, denoted as lasso W 11, P):
w12
W B3 where B is the solution the following lasso
min B
w 1/2
y + +|1B111.
(16)
where y = WT 812, and p is the regularization constant that promotes sparsity Σz
In an interventional setting 4,, the dependence Š(B) on Σg arrises from the marginal distribution of Z restricted to components indexed by Us. i.e., purely observed nodes. Since Z follows Gaussian distribution, Zue also follows Gaussian distribution with Zue (O, Σn, ). From the properties of Gaussian distribution
533 we have Σ/. EEzlu. u Consequently, for the interventional
setting It. [4] becomes
((Tx) sup
-Tr(S (E2)it log (Ez)u Ez
(17)
where Sr. = E 24, S the sample covariance of Z corresponding to the purely observed nodes. In this case, we set W = (E2),,, and the rest of the update procedure remains the same. The overall parameter update procedure summarized in Algorithm in Appendix
7

---

Graph Recowary (Nerlinoar SEM
Confounder Recovery (Noninear SEM
Graph Recovery (Nonlinear SEM
Confounder Rocovery (Nominear SEM
Nome und Sieg
0.77
0.75
0.25
0.25
Num confounders
Num
Nam
Nu.m
-
dood
-
Todays
+
10
-M dagin
- zdig
-
dood
-
modags
10
- dagma
-H LOV1
(a) Varying number of confounders
(b) Varying number of cycles
Graph (Nerlinoar SEM
Confoundor (Norinner 8EM
Recovary (Norinneur
Confoundar (Nominear SEM
1.00
Norre und SMa
0.75
0.76
0.50
0.25
0.5
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.90
0.75
1.00
0.0
2.5
5.0
7.5
10.0
0.0
2.5
5.0
10.0
Nart
Narms
inderverdion
nodaga
dagma
dood
rodage
dagma
(c) Varying number of degree of nonlinearity
(d) Varying number of interventions
Figure 2: Performance of causal graph and confounder recovery under varying problem dimensions. In all the cases the number of observed variables is fixed at = 10. (Top row, left column) number of latent confounders ranges from 2 to 8, (top row, right column) number of cycles ranges from 0 to 8. (Bottom row, left column) the degree of nonlinearity is varied between 0 and 1, (bottom row, right column) the number of training interventions is varied between O and 10.
4
Experiments
The code for DCCD-CONF is available the repository:
decd, conf
We evaluated DCCD-CONF on both synthetic and world datasets, comparing its performance against several state-of-the-ar baselines: NODAGS-FIoW (27), LLC [26), DAGMA [54], and the linear ADMG recovery method proposed by Bhattacharya al. (37 (which we refer as ADMG). NODAGS-FIOW learns nonlinear cyclic causal graphs but does not model unmeasured confounders.
LLC accounts for confounders but is limited linear cyclic SEMs. DAGMA handles nonlinearity under causal sufficiency while being limited to acyclic graphs. ADMG handles confounding but limited to acyclic graphs and linear SEMs. Note that both DAGMA and ADMG do not natively support interventional data and hence we use these models in combination with the Joint Causal Inference (JCI) framework (29) and treat interventions as multiple contexts. We also include a comparison between DCCD-CONF and two constraint based models LiNGAM-MMI JCI-FCI 29 in the Appendix (see Appendix[C)
33 and
4.1 Synthetic data
In all synthetic experiments, the cyclic graphs were generated using Erdös-Rényi (ER) random graph model with the outgoing edge density set to 2. We evaluated DCCD-CONF and the baselines on both linear well as nonlinear SEMs described in Section 2 Our training data set consists of observational data and single-node interventional over all the nodes in the graph, i.e. I (0. {1).... {d}} (unless stated otherwise), with N₂ = 500 samples per intervention. Furthermore. in all the experiments presented here, the SEM was constrained to be contractive. However, we also compare the performance of DCCD-CONF to the baselines on non-contractive SEMs in the appendix. For causal graph recovery (directed edges), we use the normalized structural Hamming distance (SHD) as the error metric. SHD counts the number of operations (addition, deletion, and reversal) needed to match the estimated causal graph to the ground truth, and normalization is done with respect to the number of nodes the graph (lower the better). For confounder identification (bidirectional edges), we compare the non-diagonal entries of the estimated confounder-non covariance matrix to those of the ground truth. We use score as the error metric (higher the better). More details regarding the experimental setup is provided in Appendix B
8

---

Impact of confounder count. We evaluate the performance of DCCD-CONE and the baselines using the previously defined enror metrics, varying the confounder ratio (number of confounders divided by the number of nodes) from 0.2 to 0.8. In this case, the number of nodes in the graph is set to 10. The results, summarized in Figure Za show that DCCD-CONF consistently achieves lower SHD across all confounder ratios in both linear and nonlinear SEMs. Notably, in nonlinear SEMs, DCCD-CONF outperforms all baselines in causal graph recovery. Additionally, it demonstrates competitive results in confounder identification, highlighting its robustness in both tasks.
Impact of number of cycles. With d = 10 nodes and confounder ratio of 0.3, we vary the number of cycles in the graph from 0 Figure compares the performance of DCCD-CONF with the baselines under this setting. As shown, increasing the number of cycles does not lead to any noticeable degradation in performance for either directed or bidirected edge recovery.
Impact of degree of nonlinearity. In this experiment, we vary the degree of nonlinearity in the SEM by adjusting between and I, where
x - ()W x) + B tanh W x + 2).
The SEM is fully linear when 3 = 0 and fully nonlinear when 8 1. Figure 2c summarizes the results. As shown, DCCD-CONE attains the highest performance as 3 approaches one, for both directed and bidirected edge recovery. When 3 is small (i.e., the system is more linear), LLC slightly outperforms DCCD-CONE, with both models performing comparably around / = 0.25.
Impact of number of interventions. In this section, we evaluate graph recovery performance as the number of training interventions K varies from 0 to d, with d = 10 fixed. The case K = 0 corresponds the observational dataset. Results for the nonlinear SEM setting are presented in Figure 2d As illustrated, with fewer interventions all DCCD-CONF and the baselines tend to exhibit similar performance (less then interventions). As the numbre of interventions increase, the performance gap widens with DCCD-CONF dominating all of the baselines. It is also worth noting that L.C. cannot operate in the purely observational setting (K 0).
Scaling with nodes. We compare the per- formance of DCCD-CONF and the baselines the number of nodes (d) varies from 10 to
Recovery 3EM
Confoundar (Noclinaar SEM)
Nomra aod SHD
0.75
80, with results summarized in Figure 3 The number of confounders is set to 0.3d. As the number of nodes increases, SHD rises across all methods, reflecting the increased difficulty of causal graph recovery in larger graphs. How- ever, DCCD-CONF consistently outperforms the baselines in many cases, achieving lower SHD and higher Fl score, suggesting superior scalability with increasing graph size.
0.50
0.00
+ DOOD-CONE
NODAGS-Foo
-H
LLC
Figure 3: Performance comparison between DCCD-CONF and the baselines and d is varied between 10 and 80.
4.2 Real World data
We evaluate DCCD-CONF
Table Results on The table presents the average Negative Log-Likelihood (NLL) on the test set, averaged over multiple trials (standard deviation is reported within parentheses).
[53]
gene perturbation dataset.
on learning the causal graph structure of a gene regula- tory network from real-world gene expression data with ge- netic interventions. Specif- ically, we use the Perturb- CITE-seq dataset [55], which contains gene expression data from 218,331 melanoma cells across three conditions: (i)
Method
Control
Co-Culture
IFN-y
DCCD-CONF
1.375 (0.103)
1.245 (0.039)
1.235 (0.338)
NODAGS
1.465 (0.015)
1.406 (0.012)
1.504 (0.009)
LLC
1.385 (0.039)
1.325 (0.029)
1,430 (0.048)
DCDI
1.523 (0.036)
1.367 (0.018)
1.517 (0.041)
control, (ii) co-culture, and (iii) IFN-y. Due to computational constraints, we restrict our anal- ysis subset of 61 genes from the 20,000 genes in the genome, following the experimental setup

---

of Sethuraman et al. [27] (see Appendix B for details). Each cell condition is treated as a separate dataset consisting single-node interventions on the selected 61 genes.
Since the dataset does not provide a ground truth causal graph, SHD cannot be used for direct performance comparison. Instead, assess DCCD-CONF and the baselines based predictive performance over unseen interventions. To evaluate performance, we split each dataset 90-10 using the smaller portion as the test set, and measure performance using negative log-likelihood (NLL) on the test data after model training (lower the better). The results are presented in Table [ From Table we can see that DCCD-CONF outperforms all the baselines across all the three cell conditions, showcasing the efficacy of the model and prevalence of confounders in real-world systems. Additionally, we also report the performance of DCCD.CONF and the baselines with respect MAE on the test data error metric in Table in Appendix C with two additional baselines: DCDFG 47 and Bicycle [56]
Additionally, we also provide results in Appendix C for the following settings: (i) performance comparison on non-contractive SEMs when the underlying graph is restricted to DAGs, (ii) performance comparison function of training data size, (iii) performance comparison as a function noise variance, (iv) performance comparison as a function outgoing edge density, and (v) performance comparison between DCCD-CONF and additional baselines: JCI-FCI and LiNGAM-MML
Additional experiments.
Discussion
In this work, we introduced DCCD-CONE, novel differentiable causal discovery framework that handles directed cycles and unmeasured confounders, assuming Gaussian exogenous noise. It models causal mechanisms via neural networks and learns the causal graph structure by maximizing penalized data likelihood. We provide consistency guarantees in the sarge-sample regime and demonstrate, through extensive synthetic and real-world experiments, that DCCD-CONE outperforms state-of-the-art methods, maintaining robustness with increasing confounders and graph size. On the Perturb-CITE-seq dataset, our model achieves superior predictive accuracy.
While the focus of this work is limited to Gaussian exogenous noise, we plan to investigate other noise distributions for future research. Other future directions include supporting missing data, and relaxing interventional assumptions by incorporating soft interventions and unknown interventional targets.
Acknowledgment
This material is based upon work supported by the National Science Foundation under Grant No. CCF-2007807 and 2502298.
References
[1] Karen Sachs, Omar Perez, Dana Pe'er, Douglas A. Lauffenburger, and Garry P. Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308(5721): 523-529, 2005.
[2] Eran Segal, Dana Pe'er, Aviv Regev, Daphne Koller, Nir Friedman, and Tommi Jaakkola. Learning module networks. Journal of Machine Learning Research, 6(4), 2005.
[3] Bin Zhang, Chris Gaiteri, Liviu-Gabria Bodea, Zhi Wang, Joshua McElwee, Alexei A. Podtelezhnikov, Chunsheng Zhang, Tao Xie, Linh Tran, and Radu Dobrin. Integrated sys- tems approach identifies genetic nodes and networks in late-onset Alzheimer's disease. Cell, 153(3):007-220, 2013.
[4] Peter Spirtes, Clark Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
[5] Soñaa Triantafillou and Ioannis Tsamardinos. Constraint-based causal discovery from multiple interventions over overlapping variable sets. The Journal of Machine Learning Research, 16(1): 2147-2205, 2015.
10

---

[6] Christina Seinze-Demll Jonas Peters, and Nicolai Meinshausen. Invariant causal prediction for nonlinear models. Journal of Causal Inference, 6(2), 2018.
[7] Christopher Meek. Graphical Models: Selecting causal and statistical models. PhD thesis,
Carnegie Mellon University, 1997.
[8] Alain Hauser and Peter Bühlmann. Characterization and greedy learning of interventional markov equivalence classes of directed acyclic graphs. The Journal of Machine Learning Research, (3(11:2409-243 224444944444 2012.
[9] Xun Zheng. Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. DAGs with NO TEARS: Continuous optimization for structure learning. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Dssa-Biacchi, and R. Garnett, editors, Advances in Neural Information Process- ing Systems, volume 31, 2018. URL http:: RtLLitttss//proceeddnge www.nnnuiiss.. wttpppsscccppppp/ppeeeepppeeee /018/
Paper pdf
[10] Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. International Conference on Machine Learning, pages 7154 7163. PMLR, 2019.
[11] Ignavier Ng. AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and DAG constraints for learning linear dags. Advances in Neural Information Processing Systems, 2020.
[12] Ignavier Ng. Shengyu Zhu, Zhuangyan Fang, Haoyang Li, Zhitang Chen, and Jun Wang. Masked gradient- graddent-based causal structure learning. In Proceedings of the 2022 SIAM International Conference on Data Mining (SDM), pages
SIAM, 2022.
[13] Xun Zheng. Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning sparse nonparametric DAGs. In Silvia Chiappa and Roberto Calandra, editors, Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108, 26-28 Aug 2020.
pages
[14] Hao-Chih Lee, Matteo Danieletto, Riccardo Miotto, Sarah T Cherng, and Joel T Dudley. Scaling structural learning with NO-BEARS infer causal transcriptome networks. In Pacific Symposium on Biocomputing 2020, pages 391-402 World Scientific, 2019.
[15] Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin. Differentiable causal discovery from interventional data. Advances in Neural Information Processing Systems, 53:21865-2187,, 2020.
[16] Ioannis Tsamardinos, Laura E Brown, and Constantin F Aliferis. The max-min hill-climbing bayesian network structure learning algorithm. Machine learning, 65(1):31-78, 2006.
[17] Liam Solus, Yuhao Wang, Lenka Matejovicova, and Caroline Uhler. Consistency guarantees for permutation-based causal inference algorithms. arXiv preprint arXiv:1702.033530 2017.
[18] Yuhao Wang, Liam Solus, Karren Yang. and Caroline Uhler. Permutation based causal inference algorithms with interventions. Advances in Neural Information Processing Systems, 30, 2017.
[19] Jacob W. Freimer, Oren Shaked, Sahin Naqvi. Nasa Simnott-Armstonng, Arwa Kathiria, Christian M. Garrido, Amy F. Chen, Jessica T. Cortez, William J. Greenleaf, Jonathan K. Pritchard, and Alexander Marson. Systematic discovery and perturbation of regulatory genes in human T cells reveals the architecture of immune networks. Nature Genetics, pages 1- July 2022. ISSN 1546-1718. doi: 10. 1038//41588-022-000 106-y. URL https://www http://ww......eeee//a//aaaa nature. com/ http:/arriclsssseeee wttp://iccceee//111-000----11110
[20] Thomas Richardson. A discovery algorithm for directed eyclic graphs. In Proceedings of the Twelfth international conference on Uncertainty in artificial intelligence, pages 454 461, 1996.
[21] Gustavo Lacerda, Peter Spirtes, Joseph Ramsey, and Patrik O. Hoyer. Discovering cyclic causal models by independent components analysis. In Proceedings of the Twenty-Furth Conference on Uncertainty in Artificial Intelligence, UAI008, page 366-374, Arlington, Virginia, USA, 2008. AUAI Press. ISBN 0974903949.
[22] Jan- -hnrihisttann Huetter and Philippe Rigollet. Estimation rates for sparse linear eyclic causal models. In Jonas Peters and David Sontag, editors, Proceedings of the 36th Conference on Uncertuinty in Artificial Intelligence (UAI), volume 124 of Proceedings of Machine Learning Research, pages 1169-1178. PMIR, 03-06 Aug 2020. URL. http:htt/:///ptt htttps//pproceednngs. .mr press/v1244 http://///eeettaaaeeeommmm

---

[23] Carlos Améndola, Philipp Dettling, Mathias Drton, Federica Onori, and Jun Wu. Structure learning for eyelic linear causal models. Conference on Uncertainty in Artificial Intelligence. pages 999-1008. PMLR, 2020.
[24] Joris M Mooij and Tom Heskes. Cyclic causal discovery from continuous equilibrium data. In Uncertainty Artificial Intelligence, 2013.
[25] Mathias Drton, Christopher Fox, and Y. Samuel Wang. Computation of maximum likelihood estimates in cyclic structural equation models. The Annals of Statistics, 47(2):663 690, 2019. doi: N0.1214/17-ASS6600 URL https: //w.i. org/oo../10 1214/11111111 AOS1602
[26] Antti Hyttinen, Frederick Eberhardt, and Patrik O Hoyer. Learning linear cyclic causal models with latent variables. The Journal of Machine Learning Research, 13(1):3387-3499 2012.
[27] Muralikrishnna G Sethuraman, Romain Lopez, Rahul Mohan, Faramarz Fekri, Tommaso Biancalani, and Jan-Christian Hütter. Nodags-flow: Nonlinear cyclic causal structure learning. In International Conference on Artificial Intelligence and Statistics, pages 6371-6387, PMLR, 2023.
[28] Peter Spirtes. An anytime algorithm for causal inference. Thomas Richardson and Tommi Jaakkola, editors, Proceedings of the Eighth International Workshop Artificial Intelligence and Statistics, volume R3 of Proceedings of Machine Learning Research, pages 278-285. PMLR, 04-07 Jan 2001. URL https: httttps://proceedinggss mlr wttp:peeeerr//rpreeeeeeeeee html Reissued by PMLR on 31 March 2021.
[29] Joris M Mooij, Sara Magliacane, and Tom Claassen. Joint causal inference from multiple contexts. Journal of machine learning research, 21 2020.
[30] Murat Kocaoglu, Amin Jaber, Karthikeyan Shanmugam, and Elias Bareinboim. Characterization and learning of causal graphs with latent variables from soft interventions. Advances in Neural Information Processing Systems, 32, 2019.
[31] Amin Jaber, Murat Kocaoglu, Karthikeyan Shanmugam, and Elias Bareinboim. Causal discovery from soft interventions with unknown targets: Characterization and learning, Advances in neural information processing systems, 33:9551-9561, 2020.
[32] Patrick Forré and Joris M Mooij. Constraint-based causal discovery for non-linear structural causal models with cycles and latent confounders. arXiv preprint /rXiv:/807. 1807.03024, 2018.
[33] Joe Suzuki and Tian-Le Yang. Generalization of lingam that allows confounding. In 2024 IEEE International Symposium on Information Theory (ISIT), pages 3540-3545, 2024. doi:
10.1
[34] Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. A linear non- non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006.
[35] Petr Ryšavý, Pavel Rytīf Xiaoyu He, Georgios Korpas, and Jakub Marečck. Exmag: Learning of maximally ancestral graphs. arXiv preprint arXiv:2503. 082455 2025.
[36] Rui Chen, Sanjeeb Dash, and Tian Gao. Integer programming for causal structure learning in the presence of latent variables. In International Conference on Machine Learning, pages 1550--5600. PMLR, 2021.
[37] Rohit Bhattacharya, Tushar Nagarajan, Daniel Malinsky, and Ilya Shpitser. Differentiable causal discovery under unmeasured confounding. In International Conference on Artificial Intelligence and Statistics, pages 2314-2322 PMLR, 2021.
[38] Alberto Abadie, Anish Agarwal, Raaz Dwivedi, and Abhin Shah. Doubly robust inference in causal latent factor models. arXiv preprint arXiv:2402. 1122.11652 1652, 2024.
[39] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. Double/debiased machine learning for treatment and causal parameters. Technical report, 2024.
[40] Kenneth A Bollen. Structural equations with latent variables, volume 210. John Wiley & Sons, 1989.
[41] Judea Pearl. Causality. Cambridge University Press, 2 edition, 2009. doi: 10. 10.1017/ #B0978051180316
12

---

[42] Stephan Bongers, Patrick Forré, Jonas Peters, and Joris M Mooij. Foundations of structural causal models with cycles and latent variables. The Annals of Statistics, 49(5):2885-2915, 2021.
[43] Jonas Peters and Peter Bühlmann. Identifiability of gaussian structural equation models with equal error variances. Biometrika, 101(1):219-225 2014.
[44]
Cheng Lu, Jianfei Chen, Chongxuan Li, Qiuhao Wang, and Jun Zhu. Implicit normalizing flows. In International Conference on Learning Representations, 2021. /ttpp/ttps//ooeerrvveee
[45] Charles G Broyden. class of methods for solving nonlinear simultaneous equations. Mathe- matics of computation, 19992):577-593, 1965.
[46] Jens Behrmann, Will Grathwohl, Ricky TQ Chen, David Duvenaud, and Jörn-Henrik Jacobsen. Invertible residual networks. In International Conference on Machine Learning, pages 573-582.
PMLR, 2019.
[47] Romain Lopez, Jan-Christian Hütter, Jonathan K Pritchard, and Aviv Regev. Large-scale differentiable causal discovery of factor graphs. In Advances in Neural Information Processing Systems, 2022.
[48] Brian C. Hall. Lie Groups, Lie Algebras, and Representations, pages 333-366. Springer New York, New York, NY, 2013. ISBN 978-14461--711655. doi: UtLpstss:
org/ 10. 1000000978--411111
[49] Ricky TQ Chen, Jens Behrmann, David K Duvenaud, and Jörn-Henrik Jacobsen. Residual flows for invertible generative modeling. Advances in Neural Information Processing Systems,
32, 2019.
[50] Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in 1076, 1989.
and Computation, 18(3): 1059
[51] Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with gumbel-softmax In International Conference on Learning Representations, 2017. URL /tperppp
[52] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432 2008
[53] Michael Jordan. An introduction to probabilistic graphical models, 2003.
[54] Kevin Bello, Bryon Aragam. and Pradeep Ravikumar. Dagma: Learning dags via m m-matrices and log-determinant acyclicity characterization. Advances in Neural Information Processing Systems, 35:8226-82333 8239, 2022.
[55] Chris J Frangieh, Johannes C Melms, Pratiksha Thakore, Kathryn R Geiger-Schuller, Patricia Ho, Adrienne M Luoma, Brian Cleary, Livnat Jerby-Arnon, Shruti Malu, Michael S Cuoco, al. Multimodal pooled Perturb-CITE-sqq screens in patient models define mechanisms of cancer immune evasion. Nature genetics, 53(3):332-341, 2021.
[56] Martin Rohbeck, Brian Clarke, Katharina Mikulik, Alexandra Pettet, Oliver Stegle, and Kai Ueltzhöffer. Bicycle: Intervention- causal discovery with cycles. In Causal Learning and Reasoning, pages 209-242. PMLR, 2024.
[57] Patrick Forré and Joris M Mooij. Markov properties for graphical models with cycles and latent variables. arXiv preprint arXiv:1710.087755 2017.
[58] Peter L Spirtes. Directed cyclic graphical representations of feedback models. arXiv preprint 1302.4332 2013.
[59] Thomas S. Richardson. A factorization criterion for acyclic directed mixed graphs. In Proceed- ings of the Twenty-Fitte Conference on Uncertainty in Artificial Intelligence, UAI '09, page 022-470, Arlington, Virginia, USA, 2009. AUAI Press. ISBN 9780974903958.
[60] Steffen L. Lauritzen and Finn Verner Jensen. Local computation with valuations from a commu- tative semigroup. Annals of Mathematics and Artificial Intelligence, 21:51-69, 1997.
[61] Peter Spirtes and Thomas S. Richardson. A polynomial time algorithm for determining dag equivalence in the presence of latent variables and selection bias. In David Madigan and Padhraic Smyth, editors, Proceedings of the Sixth International Workshop on Artificial Intelligence and
13

---

Statistics, volume RI of Proceedings of Machine Learning Research pages 489-500. PMLR, 04- 07 Jan 1997. URL htttps: by PMLR on 30 March 2021.
mlr w.l..ppees/////p2/rte//// htmll html Reissued
[62] R. Ayesha Ali, Thomas S. Richardson, and Peter Spirtes. Markov equivalence for ancestral graphs. The Annals of Statistics, 37(5B), October 2009. ISSN 0090-5364. doi: 10.1214) 08-aos626. URL htLLhttt:://x: doi org/10. 10. 10144/08AASS620
[63] Jiji Zhang. Causal reasoning with ancestral graphs. Journal of Machine Learning Research,
(7), 2008.
[64] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In Yoshua Bengio and Yann LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015. URL http:://are arxiv. otgg/ba///1...o.. 1212 6980
14

---

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 main claim and contribution of the paper is clearly stated in the abstract and at the end of introduction.
Guidelines:
The answer NA means that the abstract and introduction do not 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 is fine to include aspirational goals as motivation long 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]
We discuss the limitations and potential directions for addressing them in the discussions section.
Justification:
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 sell-speciiiatiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in 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 few datasets or 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 low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails 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 faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
15

---

Answer:
[Yes]
Justification: The full set of assumptions are provided in Appendix A and assumption is stated in Section 2. All the proofs are provided in Appendix A.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proof' in the paper should be numbered and cross referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in 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
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: All the details regarding the experimental setup and the model configurations are provided in Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to 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 a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are 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 have some path reproducing or verifying the results.
5. Open access to data and code
16

---

reproducible, they Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The code is provided the supplementary materials and will be made public upon publication of the paper.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.c cc/
public/guides/
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "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 to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/
Cooee//dddeeebmmmcia ssi onPolicy) for more details.
The authors should provide instructions 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 subset of experiments 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 to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All the details regarding the experimental setup and model configurations provided in Appendix B
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer:
[Yes]
Justification: The standard deviation over repeated trials reported in each plot and table.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of 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, or overall with given experimental conditions).
17

---

The method for calculating the error bars should be explained (closed form formula, call a 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 the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative rates).
show in tables or
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: All the compute resource for the synthetic and real-world data experiments are included in Appendix
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, 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 make into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification: This paper conforms with every aspect of NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
We do not expect the paper to have any negative social impact.
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 societal impact why the paper does not address societal impact.
18

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment 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 particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that 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 intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional 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 in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This work poses no such risk.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety 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 is 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
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?
Question:
Answer: [Yes]
Justification: See Appendix B for details on the baselines.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, 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 use in package should be provided. For popular datasets, paperswithcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
19

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer:
[NA]
Justification:
No new assets are released in this work.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions 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 used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: Human subjects were not involved in this work.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be 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 by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: Human subjects were not involved in this work.
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 be 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 to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
20

---

16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or 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, originality of the research, declaration is not required.
Answer: [NA]
Justification: LLMs are not part of this work.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to L.LM policy htte:p///eerriipc.cc/o/oootrrers for what should or should not be described.
21

---

Appendices
The appendices are structured as follows: Appendix A presents the theoretical foundations of differentiable cyclic causal discovery in the presence of unmeasured confounders, including the proof of Theorem 2 and a characterization of the equivalence class of DMGs that maximize the score function. Appendix BeaddixxBiroooe provides implementation details DCCD-CONF and the baselines. Finally, Appendix C provides additional experimental results comparing DCCD-CONF with the baselines.
A Theory
In this section, we establish the theoretical foundations of differentiable cyclic causal discovery in the presence of unmeasured confounders. We begin by reviewing key definitions and results from prior work that are essential for proving Theorem starting with fundamental graph terminology.
A.1 Preliminaries
Consider a directed mixed graph 9 = B)) A path π between nodes i and j is a sequence (io, f1,,,,,..,,, En. ), where (io...., _,. C V and {e1. ( C & U B, with in i and in path is directed if each edge EK follows the form i + in for all & E [n]. A cycle through node i consists of directed path from to some node j and an additional edge + i. For any node V, the ancestor set is defined as ang := {j E V a directed path from j to exists in G). while the descendant set is given by deg (1) == [j E V directed path from to j exists in 9). The spouse set of node i is defined as SPg (i) == [j E j + i E B). If i is both a spouse and an ancestor of j, this creates almost directed cycle. A mixed graph is called ancestral if it contains neither directed or an almost directed cycle. The strongly connected component of i, denoted seg(1), is the intersection of its ancestors and descendants: seg (i) ang (i) deg The district of node V is defined as (isg(i) = {j + ... ++ i E G or i j}. We can apply these definitions to subsets u S V by taking union of over the items of the subset, for instance, ang (24) = ( A vertex set A C V is said to be barren if i E A has no descendants in G that are in A. however, may have descendants in not A, that is, barreng (A) {i i E A: de n {i]). A subset C V is ancestrally closed if A contains all of its ancestors. We define A(G) = {A ang (A4) = A} as the set of ancestrally closed sets in G.
Definition A.1 (Collider). For a directed mixed graph G = (V,E, node is E V in path (io, E1, d1, E2..... 1,,11, En is called collider k + 0, (non-endpoint) and the two edges Eje Ek+1 have their heads pointed at i, i.e., the subpath (is-1, i Ek41, (441))4+1) is of the form i - ig <- ik+1, ig-1 <> i +- is+1,, ik-1 > ig + i 1, 4--1 ++ is +> i441. The node ig is called a non-collider if i is not a collider.
Note that the end points of walk are always non-colliders. We now define the notion of d-separation extended to DMGs.
Definition A.2 (d-separation). Let G EV,E,B) be directed mixed graph and let C C V be subset of nodes. A path π (10,E1, E1, 11, E2. in-1, Ens in) is said to be d-blocked given C if
1. π contains a collider in & ang
2. π contains a non-collider& in E C.
The path π is said to be d-open given C if is not d-blocked. Two subsets of nodes A, B CVV is said to be d-separated given C if all paths between a and b, where a A and b E B, is d-blocked given C, and is denoted by
A ALBIC.
If the underlying graph acyclic, d-separation implies conditional independence. That is, for subsets of nodes A, B, V.
d A ALBIC - XA Pg
Xc:
g
where pg denotes conditional independence, and Pg denotes the observational distribution. This is known the directed global Markov property of G [57]. However, in general, cyclic graphs do
22

---

X1
X2
X1
X2
X
X
g
acy(g)
Figure 4: (Left) Illustration a directed mixed graph that disoheys directed global Markov property.
(Right) The graph on the right represents the graph 9 after the acyclification process.
not obey the directed global Markov property as shown by the counterexample below taken from 42 58
Example A.3. Consider the SEM given by:
f1(X Zı, fi Zz.
Z3 = - + Za,
and pz is the standard normal distribution. One check that X1 is not independent of X2 given {X3, XA). However, the X1 and X2 are d-separated given X4} the graph corresponding to the SEM (see Figure 7.
Forré and Mooij [57] introduced o-separation as a generalization d-separation to extend the directed global Markov property cyclic graphs. This concept motivated by applying d-separation the acyclified version of the DMG. Before delving into o-separation, we first define the acyclification procedure of directed mixed graph, following 42
Definition A.4 (Acyclification of directed mixed graph). Let G = EV,E, B) denote directed mixed graph, the acyclification of G maps to the acyclified graph acy(G) (V, E,, B), where - i E if and only if E pag(sog(( (i)) s0g (i), and <> j E B if and only if there exists i' scg(() and j' E scg(i) such that i' = j' or i' + j e B.
It is important to note that the existence of a acylified graph for SEM relies on the solvability of the SEM over all the strongly connected components of the DMG corresponding to the SEM. This say that have solution for X (1) given X paassc(())))es
and Zzcp(i). This is indeed the case as we assume that the forward map f, is invertible for the all the SEMs under consideration, see Bongers al. 42 for more details. Figure 4 iillustrates the acyclification process for the graph corresponding to Example A.3
Definition A.5 (o-separation). Ler G = EV,E, be a directed mixed graph and let C C V be a subset of nodes. A path π (fo, E1, i1, E2...., Ens in) is said to be a-blocked given C if
1. the first node of π,, in E C' or its last node in E C, or
2. n contains collider in & ang(C)
3. x contains a non-colle non-collider is C that points towards neighbor that is not in the same strongly connected component as is in 9, i.e, such that i + in in π and it-1 & (ig() ig + i in and ig+1 & scg(is).
The path π is said to be o-open given C' if it is not a o-blocked. Two subsets of nodes A, B C is said to be o-separated given C if all paths between a and b, where and b B, is a-blocked given C, and is denoted by
A AÎBIC
Note that o-separation reduces to d-separation for acyclic graphs, that is, when scg (i) {i) for all i E V. The following result in (57) relates o-separation and d-separation.
Proposition A.6 (57)). Let g = (V,E,B) be a directed mixed graph, then for A, B, B,C
= A d B C C
B
g
23

---

X1
X2
X1
X2
X1
X2
X1
X2
C1
C2
X3
X3
X
X
G
do((Xa))
do((X)))(
G
Figure 5: Illustration of the augmented graph G corresponding to the set of interventional targets (0, {X3). {x.]). do((Xa)) and do((X3)) corresponds the graph obtained after hard interventions on X3 and X4 respectively. The augmented graph here is the union of the graphs G,
do({Xa)), do (Xa)) along with the context variables.
Using --separation we can now define the general directed global Markov property.
Definition A.7 (General directed global Markov property (57). Let G = (V, B) be directed mixed graph and pg denote the probability density of the observations X. The probability density Pg satisfies the general directed global Markov property for A, B, C.C
A G
XA
XB
Xc,
Pg
that is, XA and XB B are conditionally independent given Xc.
A.2 Joint Causal Modelling and Markov properties
In order to incorporate multiple interventional settings into single causal modeling framework, we follow the joint causal model introduced by [29), where we augment the system with set of context variables C ( (CI,.... Ck) each corresponding to a non-empty interventional setting, In this case, C = 0 for all & = 1. ..... K corresponds the observational setting. We construct augmented graph, denoted by G consisting of both the system variables X and the context variables C such that the chg (C&) I.. and no context variable has any parent or spouse. Figure 5 illustrates the augmented graph for the graph from Example A.3 and the intervention sets {0, (X²), [Xa)). The new system containing both the observed variables and the context variables is called the meta system. Finally, given a family of interventional targets = {4 and the corresponding context variable C₁, the structural equations of the meta systems governing the observations X and C has the following form:
Pi(X (a)))))) m) Z,) = (C&b [X (1) Zz), otherwise.
if E [K] s.t. Ok + 0, and X; E IN.
We call the distribution over the context variables PCC)) the context distributions and as noted by Mooij et al. [29]. the behavior of the system is usually invariant the context distribution. We
assume access to the context distribution as the interventional settings are known apriori. Note that, the observational distribution corresponds to PGz (X C1
Ck 0). Similarly, the
interventional distribution for the interventional setting I corresponds to (x C = Ere, C=+
9), i.e.,
Pgz (X) C₁ = Elg. C- = 0) = PG)
Furthermore,
Pg= (C ;x) Pgz pp= X C2).
(18)
Recall that for the interventional setting 4x. the probability density function governing the observations X is given by (6), which we repeat here for convenience
PG) (X) pz(C)pz
(x)] (EL")(()](L det (J ( (X))),
Definition A.S. Let G = (V, B) be a directed mixed graph, and I = { with Lo = be a family of interventional targets. Let Mg( denote the set of positive densities Pgz : 1d2d + R such that pgz is given by (18 for all F R2dd - Rdd, with F;(X; Z) = F. pag (i), Zi), such that the resulting forward map f, is unique and invertūble, and Σz > and (Ez)a 3))s kg # 0 if and only if i ++ j B.
24

---

Proposition A.9. For directed mixed graph G = (V,E, B) and family of interventional targets I {Ix] such that Io 9, let p E Mr(9), then p satisfies the general directed global Markov property relative to G².
Proof. For an DMG G and choice of F I2dd + R such that f, is unique and invertible and Σ > O, the structural equations are uniquely solvable with respect to each strongly connected component of G. Morover, the addition of context variables in the augmented graph does not introduce any new cycles. Therefore the meta system forms simple SCM. Thus, from Theorem A.21 in [42]. the distribution Pgz is unique and it satisfies the general directed global Markov property.
We now define the notion interventional Markov equivalence class for DMGs based on the set of distribution induced by them.
Definition A.10 (I-Markov Equivalence Class), Two directed mixed graphs G1 and G2 are I-Markov equivalent if and only if M1( Mz))),, denoted as 91 =I G2. The set of all directed mixed graphs that are I-Markov equivalent to 91 is the I-Markov equivalence class of 91. denoted as I-MEC(G1).
From Proposition A.6 for DMG 9 and family of interventional targets I ( any o-separation statement in GZ translates to d-separation statement in the acyclified graph (cy( Consequently, the acyclified graph acy(g²) is equivalent to the augmented graph G². Furthermore, by the results of [59], Pgz admits a factorization, as formalized in the theorem below.
Theorem A.11 (59).. J. A probability distribution p obeys the directed Markov property for an acyclic directed mixed graph if and only for every A A(G),
p(XA)= II p(XH (
(19)
HE|Ale
where [A]g denotes a partition of A into sets (HI...., Hx)-
Each term in the factorization above is of the form p( X XT), H,T and H = 0. Following Richardson (59). Lauritzen and Jensen (60] we refer to H as the head of the term p(XH XT), and I' the tail. An ordered pair of sets (H,,)) form the head and tail of the factor
associated with G if and only if all of the following conditions hold:
1. H barreng (aan((())))
2. If every nodes h E H is connected via a path in the graph obtained by removing all the directed edges the graph 9 when restricted to the nodes ang(H) (H), and
3. T (disang (H) H) U pag (d))
Proposition A.12. Let G = (V,E, B) be directed mixed graph and I = { be a family of interventional targets. The set of interventional distributions Pgz Mz(g)) if and only if pgz admits factorization of the form given by (19
Proof. Since any Pgz E (g) also Markov to acy acy(g*), the proposition above is a direct implication of applying of Theorem 19 on (g))
A.3 Proof of Theorem 2
We now present the main result of this paper. Recall the score function introduced in Section 3.2
K
Sr(9) = sup
E x
log logPao(L.M0) A(gg,
where (80) is the data-generating distribution for la E I. and & Σz) represents the set of all model parameters. In the context of the meta system, since we assume access to the context distribution, the score function above is equivalent the following score:
Sr) = sup E &
log Pg (X, C (Φ))= \وgkk
(X (X,C)~pt
25

---

where Pgz (X.C C Φ) is given by (18) for a specific choice of Ø. and pz denotes the joint ground-
truth distribution for the observed and the context variables. We define Pr(G) as the of all
distributions pgz (X,C' C' () that can be expressed by the model specified by equations (5) and (8) That is,
Pz(G) = 3 s.t.
- P-
(20)
From the above definition it is clear that Pz(9) S Mz(9). Theorem relies on the following set of assumptions. The first one ensures that the model is capable of representing the ground truth distribution.
Assumption A.13 (Sufficient Capacity). The joint ground truth distribution py is such that pzz E Pz()"), where 9 is the ground truth graph.
In other words, there exists & such that pi = Pgz (- Φ). The second assumption generalizes the notion of faithfulness assumption to the interventional setting.
Assumption A.14 (Z-h-faathfullness) Let (X, f. any subset of nodes A, B, VUCC, and 4 E I
A
B C = VA VB Vc. gz
PGz
The above assumption implies that any conditional independency observed in the data must imply o-separation in the corresponding interventional ground truth graph.
Assumption A.15 (Finite differential entropy). For {Ix]
Epp log (X, C)I
The above assumption ensures that the hypothetical scenario where S(g*) and S(g) are both infinity is avoided. This is formalized in the lemma below taken from [15].
Lemma A.16 (Finiteness of the score function 15)). Under assumptions 00.
A.13 and A.15 [Sx(9)
<
From the results of [15], the minimization KL diverengence plus the difference the regularization terms.
express the difference in score function between 9 and G as
Lemma A.17 (Rewritting score function [5]). Under assumptions (3.3 and A.75 15 we have
S(g*) S(G) inf DEL (- ()) +
We will now prove the following technical lemma (adapted from proving Theorem 2
15 which we will be used in
Lemma A.18. Let (V,E,B) be directed mixed graph, for of interventional targets H and pa & MM(()), then
inf D(p" (p) > pEM,(9))
Proof. Let V = (X,C4 from theorem A.ll any E Mr(G))
admits factorization of the form
p(V)
II p(V VT). HE(V]g
Let us define new distribution follows:
p(V)=
II P(V Vr).
HelVle
where
p" (VH, VT)
Φ(H VT)
1² (V+)
26

---

From proposition A.12 we see that P E MM) and hence + P. We will show that
p = arg pEM,((9)
min (K (p)))
For an arbitrary p E MM())), consider the following:
Ep- log
(VV
Exp
E
log
D H
Vr)
(21)
P(V)
p VH Vx)
HeVVa
p P'(HH VT)
Σ
Ep
log
(22)
p(v
V+)
HEVI
In the equation above, leverage the linearity of expectation, which holds under Assumption A.15 ensuring that we don't sum infinities of opposite signs. We now show that each term in the right hand side of the above equality is expectation of KI. divergence which is always in [0, 00).
Egg log p² (Vn
VT) PVH VT)
/ 1² p'(n -
p Q"(nn VT)
VT log
dVndV
(23)
P(VH VT)
| Dv I V dVr.
VT)
(24)
Thus, Ep- log P(V) p(V) E [0, oo).
We now show that p = arg M(dd((((() (ki(p"le2
p" (V) P(V)
(KI(P" (lp) Ryp log
(25)
((V) P(V)
p (V)
p0 V)
Ep log
+ Ep- log
(26)
p(V)
p (V)
p V)
DxI(p" (1)) + Exp log
(27)
p( (V)
2 Dri(p" (1)
(28)
Since the expectations in 26 are both in [0, oo), splitting the expectation is valid. The very last inequality holds since p" p. Thus,
inf (0M
D(p 1p) 2 (KL (10) > 0.
This proves the lemma.
We are now ready prove Theorem
2
Recall,
Theorem 2 Let = {Ix R=o
{ be a family of interventional targets, let G denote the ground truth directed mixed graph, ((t)) denote the data generating distribution for Itp and G := arg maxg S(G). Then, under the Assumptions/ A.13 A.14 = g" That is, I-Markov equivalent to g*
15 and for a suitably chosen X >> O, we have that
Proof. It is sufficient to show that for G & I-MEC(G*). the score function of 9 is strictly lower than the score function of g", i.e., S(G*) > S(G). Since I-MEC(G") and pz MMz(g*) (by Assumption A.13), there must exist subsets of nodes B, C CCVCCC such that either:
A A C and 9
A B C g.
or
A ALBIC and A AÍC.C.
no such subsets exist, then G and 9* impose the same -separation constraints and thus induce the same set of distributions. This would imply that G E I-MEC(G"). contradicting assumption.
27

---

Since p1 E Mr(())), it must be true that VA 4p00) VB Ve (Assumption A.14 Therefore pI* doesn't satisfy the general directed Markov property with respect to G and hence p & Mr) For convenience, let
(g) == inf DKL (PLlpilppgz( Φ)).
Note that
(g)) = inf Dx (pi Pp= |-1 Φ)) 2 peM)((0
inf DKL (p (k) [lp) > 0,
6
where we use Lemma A.18 18888
for the final inequality. Thus, from Lemma A.17
S(G") - S(G) 7(9) X( (
(29)
Following 15
we now show that by choosing À sufficiently small, the above equation is stictly positive. Note that if [g] 2 | then S(g*) - S(G) O. Let G == {9 | < |g*1}. Choosing
such that 0 < A < minger+ 7(
we see that:
()(
min GEa
(30)
9- |
(9))
4G E G+
(31)
9- G]
+ X(69) -91) < (() vg E G+
(32)
- 0 (9) + A(19) (g = S(g)) Ag e G+
(33)
Thus, every graph outside of the general directed Markov equivalence class of (g*)² has a strictly lower score.
A.4 Characterization of Equivalence Class
Let 9 BV) be directed mixed graph, and consider family of interventional targets I I= I=== with 0. From Proposition A.6 for any graph Gi I-MEC(G), the corresponding augmented graph of is equivalent to the acyclification acy(g²) of G². Several prior works have studied the characterization of equivalence classes of acyclic directed mixed graphs (ADMGs),, including [61] 30, 62). We now provide graphical notion of the I-Markov equivalence class of DMG G. A graph G is said to be maximal if there exists no inducing path (relative the empty set) between any two non-adjacent nodes. An inducing path relative to subset 7. is a path on which every non-endpoint node É L is a collider on the path and every collider is an ancestor of an endpoint of the path. A Maximal Ancestral Graph (MAG) is one that is both ancestral and maximal. Given an ADMG aCy( it is possible to construct a MAG over the variable set V (X, C that preserves both the independence structure and ancestral relationships encoded in acy (g see [63 for details. We denote M (Gca (G=) )) to mean MAG that is constructed from acy Therefore all the independencies encoded in acy (g²) is also present in M (acy Aa(cy((2²) Before present the condition for MAG equivalence, we introduce the notion of discriminating path. A path in) in (C) is called discrininating path for in-1 if (1) x includes at lest three edges; (2) 1 is non-endpoint node on म, and is adjacent to in on π.; and (3) in and in not adjacent, and every node in between in and i is collider æ and is parent of i,. The following theorem from Spirtes and Richardson [61] characterizes the equivalence of MAGs.
(fo. E..... 1,-1,
Theorem A.19 (Spirtes and Richardson 61 only if:
Two MAGs 91 and G2 are Markov equivalent and
1. 9i and G2 have the same skeleton;
2. G1 and G2 have the same unshielded colliders; and
3. if n forms a discriminating path for i in 91 and Gz, then i a collider on π if and only it is collider on in G2-
Therefore, from Proposition and Theorem A.19 two DMGs 9i and Ga are equivalent and only M AG (GF)) and M AG[acy(G3) satisfying the conditions of Theorem i.e., M AGaacy and M (a)y(()))) )): (i) have the same skeleton, (ii) same unshielded colliders, and (iii) same discriminating paths with consistent colliders.
28

---

B Implementation Details
this section we provide the implementation details of DCCD-CONF and the baseline models along with the details of the experimental setup.
B.1 Hutchinson trace estimator for computing log determinant of the Jacobian
Computing the log determinant of the Jacobian matrix present in (6) poses significant challenge. However, following Behrmann et [46], in Section 3.3.1 showed that the log-determinant of the Jacobian can be estimated using the following estimator
a
(-1) Tr (X)) - Tr (Z)}
log det (J (J [x))()))]
E (epp((N)
m
pr(( 2 m)
The estimator above still has major drawback: computing the Tr(Ju.g) still requires (d gradient calls compute exactly. Fortunately, Hutchinson trace estimator [50) can used to stochastically approximate the trace of the Jacobian matrix. This then results in the following estimator that can be computed efficiently via reverse-mode automatic differentiation
1)
log det (J (X))
E
V
U.g. (X) V TT.J..g. V
(34)
prlle 2 m)
B.2 Parameter update via score maximization
As described in Section the model parameters are updated in two stages. In the first stage, the param- of the neural networks and the Gumbel-Softman distribution, used sample adjacency matrices, are updated via backpropagation using stochastic gradient descent. Since (5) forms an implicit block of an implicit normalizing flow, following 44], we directly estimate the gradients of S(B) with re- spect to x and Φ = (0, B). The gradient computation involves two terms: Ø(-) log det(I+Jge (æ, Φ))
05 0z and dz .(.) where (.) is placeholder for and Φ. From estimators for the gradients:
46
the following unbiased
i log del(I JR. (x, Φ))
(
1-1)* P(N k)
OJ g. (a, Φ)) v.
E V
VTJ²)
(35)
0(-)
8(.)
On the other hand,
BS(B) oz o
can be computed according the implicit function theorem as follows:
D$(B) Oz
S(B)
G(x,z, 0(.)
8z Jc. (z)
(36)
Oz 0())
where G=(2) (z. z, and recall that G(a, 2, Φ) = g. Φ) + x + +z
Φ) + z. See
44 for more details. The procedure SGUPDATE shown in Algorithm [ performs the gradient
computation in 35 and (36).
In the second stage, the entries of the covariance matrix of the endogenous noise distribution are updated column-wise by solving a sequence of Lasso optimization problems. The complete parameter update procedure is summarized in Algorithm
B.3 DCCD-CONF and the baselines code details
DCCD-CONE. We implemented our framework using the libraries Pytorch and Scikit-learn in Python and the code used running the experiments can be found in the following Github repository: https: : thub. com www.murwwiggs////oooooooooooeeee hccpi//sdddddcodooltl
29

---

Algorithm PARAMETER UPDATE
Require: Family of interventional targets ) 1==1" interventional dataset {x(i.k) regularization coefficients À and p.
J... 1==1,==1'
Melk
Ensure: Learned neural network parameters 0, graph structure parameters B, confounder-noise distribution parameters Σz.
1: Initialize the parameters: ((0) ~ po((), B(0) pB((),, and Σg =
2: Iteration counter: t
3: while NOT CONVERGED do for k to K do
4:
5:
1 < +
6:
W < --E% Ju
Compute score function (()()) (()) W, Tx) ((0-1),0(0+)) SGUPDATE( Bio, 010) for j- to do
7:
8:
9:
10:
Push j-th row and column in W to the end
11:
3 - lasso (W 11. $12, p)
12:
W12 <- W
13:
end for (
14:
lu uu W
15:
end for
16:
end while
17: 0, B, Σg + ((1) E(1)
return 0, O,
Starting with an initialization of the model parameters (0(") ( E()) we iteratively alternate between maximizing the score function with respect to (0(1) B()) and {tt) as described in Algo- rithm Standard stochastic gradient updates are used for (o()) ( while coordinate gradient descent, implemented via the Scikit-leare library, applied E/ For modeling the causal function ga. we follow the setup of Sethuraman et al. [27). employing neural networks (NNs) with dependency masks parameterized by Gumbel-softme distribution. The log-determinant of the Jacobian is computed using a power series expansion combined with the Hutchinson trace estimator. To mitigate bias from truncating the power series expansion, the number terms sampled from Poisson distribution, as detailed in Section 3.3 and Appendix B.I The final objective is optimized using the Adam optimizer (64]
The learning rate in all our experiments was set to 10-2, The neural network models used in our multi-layer perceptron layer. No nonlinearities were added the neural
experiments contained networks for the linear SEM experiments. We used tanh activation for the nonlinear SEM experi- the Berturb-CTTE--sq data set. The graph sparsity regularization
ments and for the experiments constant X was set to 10-2 for all the experiments. The sparsity inducing regularization constant for the inverse covariance matrix of the confounder distribution, P. was set 10 in all the experiments. The models were trained and evaluated on NVIDIA RTX6000 GPUs.
For NODAGS-Flow, we used the code provided by authors 27 available at |https //github. woam//ennntecch/ee wwwp:hhhh//ooagaoooooooooo The default values were set for the hyperparameters. We implemented the LLC algorithm based on the details provided in (26). The implementation can be found within the codes/ baselines folder in the supplementary materials. For FCI, we used the implementation that is available in the causallearn python library /teetps:/ggitie For DCDI, we used the codebase provided by the authors [[S]. available at httphttps:s http:/gtttttttt/.e/ com/ s http://slcaaaalleeeeeaaaa /cdii The default hyperparameters were used while training and evaluating model. For DAGMA and ADMG, we used the codebase provided by the authors [54] and [37], available at wntttts:///. //github.. com htte:/evinnsbllllon. dagma and https //gitlab. wtt.w/owtttt.coc/eeeeeeeo respectively. We implemented LiNGAM-MMI based on the details
Baselines.
com/py- http://ywwwyycuusa--leeeeee causal-
provided by the authors [33].
30

---

B.4
Experimental setup
In this section, we describe how the data sets were generated for the various experiments conducted.
B.4.1 Synthetic Experiments
We begin by sampling directed graph using the Erdös-Rényi (ER) random graph model with an edge density of 2 unless specified otherwise, which determines the directed edges in the DMG G. Next, we generate a random matrix and project it onto the space of positive definite matrices to obtain the confounder covariance matrix Σg. setting the maximum exogenous noise standard deviation to 0.5 unless specified otherwise. The nonzero off- off-diagonal entries of Σg correspond to the bidirectional edges in G. For the linear SEM, edge weights are sampled uniformly from Unif((-0.. (-9.9, -..) U (0.2, 0.9)). In all experiments except those on non-contractive SEMs, the edge weight matrix rescaled ensure Lipschitz constant of less than one. For nonlinear SEMs, we apply a tanh nonlinearity to the linear system defined by the edge weights, i.e.,
ac Lanh
2);
where W the weighted adjacency matrix. In all the experiments, the training data consisted of 500 samples per interventional setting (unless specified otherwise).
Impact of Confounder Count is fixed 10. Training data consists of combination of observational data and single-node interventions over all nodes, i.e., I = U {{i} i E [d]]. The confounder ratio (number of confounders divided by the number of nodes) is varied from 0.2 to 0.8.
In this experiment, the number of observed nodes in the graph
Impact of Cycles
We fix the number of nodes the graph = 10, The number of cycles in the graph is varied between 0 and 8, with the confounder ratio set to 0.3. The training data consists of observational data as well as single node interventions over all the nodes in the graph.
Impact of Nonlinearity In this setting, the degree of nonlinearity (controlled by () is varied between and 1. That is, the data is generated from the following SEM:
- B
=) + tanh(W
+ z).
Here, 3 implies the SEM purely linear and 3 1.0 implies the data purely nonlinear. The confounder ratio is set to 0.3. The number of cycles is randomly set.
Scaling with Number of Nodes We fix the confounder ratio at 0.4. The total number of nodes in the graph is varied from 10 to 80. As in the previous setup, training data consists of combination of observational data and single-node interventions across all nodes, with 500 samples per interventional setting.
Scaling with Interventions We fix d 10 and set the confounder ratio to 0.4. The number of interventions during training varies from to Zero interventions corresponds to observational data. When fewer than interventions provided, the intervened nodes selected arbitrarily. Each interventional setting consists of 500 samples.
In this case, explicitly enforce non-contractive causal mechanism F by rescaling edge weights to ensure that the Lipschitz constant of the edge weight matrix exceeds one. We set 10 and provide observational data and single-node interventions across all nodes, with 500 samples per intervention. The confounder ratio varies between 0.2 and 0.8.
Non-Contractive SEM
Scaling with Training Samples To examine the sample requirements of DCCD-CONE, set the confounder ratio to 0.4. Training data consists of observational data and single-node interventions over all nodes, while the number of samples per intervention is varied from 500 to 2500.
Scaling with outgoing edge density In this case, outgoing edge density of the ER random
graphs is varied from I to 4. The confounder ratio is set to 0.4 and the number of nodes d = 10. The training data consists of observational data and single node experiments over all the nodes in the graph.
31

---

Scaling with noise standard deviation In this setting, we vary the maximum noise standard deviation between 0.2 and 0.8. The confounder ratio is set to 0.4 and the number of nodes d 10. The training data consists of observational data and single node experiments over all the nodes in the graph.
Evaluation Metrics Across all experiments, we use Structural Hamming Distance (SHD) to evaluate the accuracy of the estimated directed edges relative to the ground truth. SHD measures the number of modifications (edge additions, reversals, and deletions) required to match the estimated graph to the ground truth. For DCCD-CONT and NODAGS-FUS we fix threshold value of 0.8 for the estimated adjacency matrix. The recovery of bidirectional edges assessed using the score, which is defined as:
precision recall
Fl score
2 precision + recall
where
TP
TP
precision
TP+ + FP recall
TP + FN
and TP, FP, FN denote true positives, false positives, and false negatives, respectively. use threshold of 0.01 for the estimated covariance matrix to identify the bidirectional edges.
Additionally, we also measure the performance of DCCD-CONF and the baselines using Area Under Drecision-Realle Curve (AUPRC) as the eπror metric. AUPRC computes the area under the precision-recat curve evaluated at various threshold values (the higher the better).
B.4.2 Gene Perturbation Data set
The dataset was obtained from the Single Cell Portal the Broad Institute (accession code SCP1064). Following the experimental setup of Sethuraman al. (27 filtered out cells with fewer than 500 expressed genes and removed genes expressed in fewer than 500 cells. Due to computational constraints, selected subset of 61 perturbed genes (Table 2 from the full genome. The three experimental conditions cenditions-cocultuuee and models were trained and evaluated on each condition independently.
IFN-y, and control-wert partitioned into separate datasets,
Table 2: The list of chosen genes from
dataset
155
ACSL3
ACTA2
B2M
CCNDI
CD274
CD58
CD59
CDK4
CDK6
CDKNIA
CKSIB
CST3
CTPSI
DNMTI
BIF3K
EVAIA
FKBP4
FOS
GSEC
GSN
HASPIN
HLA-A
HLA-B
HLA-C
HLA-E
IENGRI
IFNGR2
ILF2
IRF3
JAKI
JAK2
LAMP2
LGALS3
MRPL47
MYC
P2RX4
PABPCI
PAICS
PETIOO
PTМA
PUP60
RNASEHZA
RRSI
SATI
SEC1IC
SINHCAF
SMAD4
SOX4
SP100
SSR2
STATI
STOM
TGFBI
TIMP2
TM4SFI
TMEDIO
TMEMI73
TOPIMT
TPRKB
TXNDC17
VDAC2
C Additional Results
C.1 Additional synthetic experiments
Experiments on Non-contractive DAGs. We evaluated the performance of DCCD-CONE and baseline methods on non-contractive SEMs, where the ground truth DMG acyelic. While DCCD- CONF assumes contractive causal mechanism, we adapt it for non-contractive settings using the preconditioning trick proposed by Sethuraman diagonal preconditioning matrix A, transforming the causal mechanism as follows:
(27. This approach introduces learnable
6x A gx A.
where gx as defined K8) remains contractive (see Sethuraman al. 27 details). We vary the confounder ratio, and the results are summarized in Figure 6al As shown in Figure, DCCD-CONF effectively learns the ADMG even in non-contractive SEM settings, demonstrating competitive performance against the baselines.
32

---

Graphi Recovery Nun aa
Confouunerr Racovaty (Non Non-coottrecive
Graph Recovory (Noninncar
Confounder Rocovery (Noniinar SEM
1.00
1.5
0.75
0.75
6.-
0.00
0.30
0.25
0.25
0.20
0.00
920
1500
2000
1000
2500
- DCC0-GONI
-
MEDAGGS-TW
-
LLC
-W
DC31
+ DCCD-CONT
- NODAGS- Fm
-
LLC
- DGBI
(a) Noncontractive SEM
(b) Performance V. sample size
Graph Recovery (Nonliinarr SEMI
Confounder Recovery (Noninnear SEM
Graph Recovary (Noniineur SEM
Confoundor (Nortnnaa SEM
20
Nomazze 9t
1.5
60
00
0.0
-+
DOCD-CONF
NDDAGS-FOW
-
DODI
DOCD-CONF
NODAGS-FIo
LLC
-
DODI
(c) Performance v. Noise deviation
(d) Performance outgoing edge density
Figure 6: Performance comparison between DCCD-CONF and baseline methods on causal graph recovery and confounder identification, evaluated across varying model parameters. Each subplot details specific experimental setting.
Performance comparison vs. sample size We also assess the sample requirements DCCD- CONF. Figure 6b summarizes the results obtained by varying the number of samples per intervention. As shown in the figure, when the confounder ratio 0.4, DCCD-CONE achieves low SHD even with 500 samples per intervention and attains near-perfect accuracy from 2000 samples onward. However, performance declines slightly as the confounder ratio increases.
Performance comparison max endogenous noise st. deviation In this setting, compare DCCD-CONE with the baseline by varying the maximum standard deviation of the endogenous noise terms between 0.2 and 0.8, the results are summarized in Figure 6cc DCCD-CONE outperforms the baselines for all noise standard deviations. However, the performance of the models does deteriorate slightly as the noise standard deviation increases.
Performance comparison vs. outgoing edge density In this case, the expected number of outgoing edges from each node is varied between I and 4. This affects the sparsity of the resulting graph. The results are summarized in Figure 6d) As seen from Figure 6d. DCCD-CONF still outperforms the baselines, even though the performance of all the models worsens as the edge density increases.
C.2 Additional performance metries
In addition to SHD for directed edge recovery and Fl score for bi-directional edge recovery. We also AUPRC to compare DCCD-CONF and the baseline for all of the experimental settings stated in Appendix B.4.1. The results are summarized in Figure Overall, DCCD-CONF performs better than LLC on nonlinear SEMs across all the settings, while achieving perfect AUPRC scores in several cases.
C.3 Comparison with FCI-JCI
Here, compare the performance of DCCD-CONF with FCI-JCI 22),, which is an extension of FCI algorithm that is capable of handling multiple contexts (in this case interventional settings). FCI-JCI outputs Partial Ancestral Graph (PAG), which graph structure that represents the equivalence class of MAGs. We define modified SHD score in order to check if the DMG estimated by DCCD-CONF belongs to the same equivalence class of the ground truth DMG. To that end, we convert the ground-truth DMG and the estimated DMG to their augmented DMGs and then construct the MAG of the acyclified version of the augmented DMGs. The modified SHD score then computes the discrepancies in the conditions of Theorem A.19 i.c., we count: (i) the
33

---

Graph Recovery
Confounder Recovery
Graph Recovery
Confounder Reaovrry
1.00
1.00
1.00
1.00
AUPRO
0.76
0.75
0.76
0.70
0.50
0.2
0.6
0.8 0.2
0.6
0.8
10
40
00
10
20
40
Confounder
Confounder
DGUD GONF
DGED GONF
LLG
(a) Performance V. number of confounders
(b) Performance V. number of nodes
Graph Recovery
Confounder Reaovery
Graph Recovery
Conlounder Reaevery
1.00
1.00
AUFRO
0.75
0.75
0.75
0.75
0.50
0.50
0.50
0.50
0.26
0.85
0.26
0.85
0.00
Interventions
Interventions
samplex
samplex
-W DCCD.CONE CONF
DCCD-CONF
- LLC
(c) Performance v. number of interventions
(d) Performance sample size
Graph Recovory
Contounder Recovery
Graph Recovery
Confoundor Recovery
1.00
1.00
AUFRC
0.50
0,60
0.50
0.25
0.55
025
0.25
0.00
0.00
0.00
0.00
0.2
0.4
0.6
0.8 0.2
0.4
0.6
0.8
Max nolce dov
Max noice dev
Edge density
Edge density
- DOCID-CONF
-M--
LLC
+ DOCD-CONF
LLG
(e) Performance V. max noise st. dev
(1) Performance v. edge density
Figure 7: Comparison of DCCD-CONE and bascline methods on causal graph recovery and con- founder identification, measured using AUPRC across varying model parameters.
intervvntiooss
2 interventions
10.0
Modntd SHD
Modifeed
5.0
2-
2.5
0-26
0
0.2
0.4
0.6
0.8
0.2
Confounder nallo
Confounder ratio
- DCCD-CONF
- FCNJCI
Figure 8: Performance comparison between DCCD-CONF and FCI-JCI with respect modified SHD as the еттоr metric. The number of observed nodes was set to d = 5. For the left plot, all single node interventions along with the observational data were provided as training data. For the right plot, observational data and interventions over two of the nodes (randomly chosen) were used as training data.
number of extra edges (N₂) using the skeletons of the estimated MAG and ground truth MAG, (2) number of mismatched unshielded colliders (N2), and (3) discrepancies in the discriminating paths (Ng). Similarly, for FCI-JCI, we count the disagreements between the ground -ruth MAG and the estimated PAG, i.e., mismatch in skeleton (N1), mismatch in unshielded colliders (N2), invariant edge orientation discrepancies (Na). Finally, the modified SHD = N₁ + N₂ N₂- We compare DCCD-CONF and JCI-FCI over two different settings: (i) the training data consists of observational data and single node interventions over the nodes in the graph, and (ii) the training data consists of observational data and interventions over 2 nodes (randomly chosen) in the graph. Due to the
34

---

complexity of computing the modified SHD (as it involves iterating over the discriminating paths and inducing paths) we fix the number of observed nodes be d 5. In the both the cases, the confounder ratio is varied between 0.2 and 0.8, and nonlinear SEM is used to generate the data. The
results are summarized in Figure g
As seen from Figure 8 DCCD-CONF outperforms FCI-JCI in the both the settings. However, the performance does decrease the number of training interventions reduces. We attribute this to the increase sample requirements the number of training interventions goes down.
Graph Recovery (Noninear SEM)
Normalzedd SHD
C.4 Comparison with LINGAM-MMI
0.2
0.4
0.6
Confounder ratio
We compare the performance of DCCD-CONE with LiNGAM- MMI [33], which extends the ICA-based LiNGAM 34 to handle hidden confounding. Since LiNGAM- requires iterating over all possible node permutations, its computational cost grows rapidly with the number of nodes; hence, we restrict our comparison d 5. The training data include both observational samples and single-node interventions for all nodes, with 500 samples per interventional setting. We vary the confounder ratio (i.c., the ratio of confounders observed nodes) between 0.2 and 0.8, As shown in Figure DCCD-CONF consistently outperforms LiNGAM-MMA across all tested confounder ratios.
-X- dood - Irgaa-mmi
Figure 9: Performance comparison between DCCD- CONF and LiNGAM-MMI on d = 5 node graphs.
C.5 Hyperparameter Sensitivity
We evaluate the sensitivity of DCCD CONF its hyperparameters: (i) the directed edge sparsity reg- ularization coefficient Ae, and (ii) the bidirected edge sparsity regularization coefficient p. Figure [0] summarizes the results for Åes م E 0.1, 0.01, 0.001 graphs with d = 10 nodes and a confounder ratio of 0.3. The model was trained using both observational data and all single-node interventions. As shown in the figure, DCCD-CONE remains fairly robust to hyperparameter variations, achieving normalized SHD values below 0.3 and F1 scores above 0.65 for most combinations of X and
Illustration of DCCD-CONF performance for various choices of hyperparameters A (directed edge sparsity regularization coeff.) and (bidirected edge sparsity regularization coeff.).
Figure 10:
C.6 Training Time Comparison
Figure I1 compares the training times of DCCD-CONE and baseline methods. Unlike the other algorithms, LLC does not require stochastic gradient-based training, it relies solely on solving a series of linear regressions. Consequently. LLC is considerably faster, shown in Figure [T Excluding LLC, NODAGS-FIUL Flow is the most efficient in terms of runtime; however, cannot account for confounders within its framework. DCCD-CONF achieves training times comparable to ADMG and DCDI while simultaneously handling confounders, cycles, and nonlinear dependencies. The training times are computed on d = 10 node graphs with training dataset consisting of observational and all single-node interventions. DCCD-CONF was trained for 200 epochs. All models were run on RTX6000 GPUs.
C.7 Additional Results on Perturb-CITE-seq Dataset
In addition to test-set NLL., we evaluate the performance of DCCD-CONF and the baselines on the Perturb-CITE CITE-seg dataset (55) using Interventional Mean Absolute Error (I-MAE) the evaluation
35

---

Training tim
1990
600
400
rsdege
dggaa
Figure 11: Training time comparison between DCCD-CONF and the baselines
metric. I-MAE computed as the mean of x + g=(x))l over all observations æ in the held-out test set. Beyond the baselines discussed in Section also include Bicycle (56) and DCDFG [47] for comparison. Bicycle supports nonlinear and cyclic structures, whereas DCDFG restricts the search space acyclic graphs. The results, summarized in Table show that DCCD-CONF remains competitive with state-of-the-ar methods under the I-MAE metric.
Table 3: Results Perturb-CITE-seq 55 gene perturbation dataset. The table presents the average Mean Absolute Error (MAE) on the test set, averaged over multiple trials (standard deviation is reported within paranthesis).
Method
Control
Co-Culture
IFN-7
DCCD-CONT
0.781 (0.037)
0.765 (0.046)
0.843 (0.035)
NODAGS
0.847 (0.018)
0.762 (0.018)
0.861 (0.023)
Bicycle
0.782 (0.042)
0.735 (0.036)
0.883 (0.028)
DCDFG
0.845 (0.066)
0.774 (0.038)
0.891 (0.041)
Additionally, report the adjacency matrix of the recovered causal graph for the cell condition "Co-Culture" in Figure [2] DCCD-CONF identified 38 feedback cycles. This number validates prior work showing that gene regulatory networks rich in feedback loops [55]
Figure 12: Adjacency matrix of learnt by DCCD-CONF for "Co-Culture" cell condition of Perturb- CITE-seq dataset.
36

---

C.8 Protein Signaling Dataset
We further evaluate DCCD-CONE on biological dataset for protein signaling network discovery [ which is widely used as benchmark for causal discovery algorithms.
The dataset contains continuous measurements of multiple phosphory- lated proteins and phospholipid components, in human immune system cells, with the corresponding network capturing the ordering interac- tions among pathway components. Based on n. 7466 samples across m 11 cell types, Sachs al. identified 20 edges in the underly- ing graph. Using the consensus network from [ as ground truth, we evaluate performance using the Structural Hamming Distance (SHD) as the error metric. The results, summarized in Table 4 show that DCCD-CONF performs comparably to the baselines. The recovered
Table 4: Performance com-
parison on Sachs et al. [
protein signaling dataset.
Method
SHD
DCCD-CONF
18
DAG-GNN TO
19
DAGMA
21
NOTEARS
22
directed graph is visualized in Figure 13
True Graph (Directed)
DCCD-CONF: Estimated Graph (Diested)
Vic
644
PP2
PP0
46-
46
Mor
nr2
nro
Ei
Figure 13: (Left) Consensus graph from Sachs al. (1] (right) graph learned by DCCD-CONE.
37


Reviews:
[{"review_id": "Eob4q6zIov", "reviewer": "Reviewer_LryH", "summary": "The authors correctly point out that much of the work on causal learning is rather limited, in the sense that it focusses on models that are directed acyclic graphs (DAGs), often without considering confounders and non-linearities. The authors' focus on cycles without considering time-series data on the input, however, obscures the fact that dynamic Bayesian networks (where cycles in the intra-slice graph are replaced by relationships of the inter-slice kind, and thus \"cycles\" are allowed only with some delay) or some other \"dynamic\" causal models may be more appropriate in many cheminformatics and bio-medical applications. Then, they proceed to formulate a variant of causal learning, where deep learning learns a model with cycles and with confounders. In terms of optimization, this can be easier to optimize over directed graphs or mixed graphs than dealing with the cycle-avoidace in DAGs. Unfortunately, there is little analysis and only a very limited empirical study.", "questions": "Q1: what was the research question you tried to answer with the empirical study? What was the reasoning for picking NODAGS-Flow, LLC, DCDI as the methods you compare against? \n\nQ2: could you estimate the effects of non-linearity, effects of considering cycles, and effects of working with confounders on the statistical performance? \n\nQ3: what would be the statistical performance on any of the commonly considered benchmarks? For example DREAM4, https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0013397, is used across many causal learning methods, including those disregarding cycles etc.", "limitations": "The claims are not particularly well supported by the empirical evidence.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths:\n\nS1: The method has a number of desirable features (allowing for confounders, non-linearity, and cycles), but still seems to scale better than many methods for learning DAGs (other than FCI and their variants). \n\nWeaknesses: \n\nW1: The empirical comparison seems pooly designed. The authors present a comparison only against NODAGS-Flow (a method for causal learning with cycles from UAI 2023), LLC (a method for causal learning with cycles and latent variables, JMLR 2012), DCDI (a differentiable approach focussing on working with interventional data, NeurIPS 2020), with focus on Perturb-CITE-seq, a little known dataset. There is no ablation study in the sense of comparing what are the effects of the non-linearity, what are the effects of considering cycles, and what are the effects of confounding, nor any comparison against the best-cited methods in the field that disregard confounding and cycles, but do consider non-linearity (e.g., DAGMA, NeurIPS 2022, or variants of FCI, UAI 1995), or methods that diregard non-linearity, but consider confounding. \n\nW2: I appreciate that Perturb-CITE-seq presents an important challenge in the processing of genomics data, but it has not been used only in a couple of papers on causal learning, yet, which makes it hard to evaluate the performance of the method from the results on Perturb-CITE-seq. Having said that, including the results of DCD-FG (NeurIPS 2022), Bicycle (ICLR 2024), and DAG-BLR (Nature Machine Intelligence 2023), which have tested on Perturb-CITE-seq, would be fabulous, considering that they are already available. \n\nW3: confounding has been studied in many papers in causal learning, including the doubly-robust inference (e.g., https://arxiv.org/abs/1608.00060, https://arxiv.org/abs/2402.11652), differentiable methods (e.g., https://arxiv.org/abs/2010.06978), methods based on LinGAM (https://arxiv.org/abs/2401.16661), or integer-programming methods for learning mixed graphs (e.g., https://arxiv.org/abs/2102.03129, https://arxiv.org/abs/2503.08245). Most of these methods are not cited, nor they are compared against. Considering that many of these methods are exact, it would be easy to estimate the effect of working with confounders.", "was_revised": true, "final_justification": "I have been well impressed by the additional experiments conducted by the authors and raised my score to 5. Thank you!", "rebuttal": "We thank the reviewer for their thoughtful and constructive feedback. We greatly appreciate the recognition of our method’s ability to jointly handle confounders, non-linearity, and cycles while remaining computationally scalable compared to many DAG-based approaches, as well as the acknowledgment of the relevance of our problem setting in challenging real-world applications. Below, we address specific comments and questions raised. \n\n## Regarding Q1 and the first half of comment W1\n\nWe thank the reviewer for this question. Our empirical study was designed to evaluate DCCD-CONF in settings where **confounding, cycles, and interventions** jointly occur, scenarios where most existing methods are inapplicable or have limited support. Accordingly, we selected baselines that isolate specific factors:\n- **LLC** [1]: handles cycles and confounding but is linear, highlighting the value of nonlinearity.\n- **NODAGS-Flow** [2]: supports nonlinearity and cycles but no confounding, isolating the effect of latent confounders.\n- **DCDI** [3]: handles interventions but assumes acyclicity and no confounding, representing the standard differentiable DAG-based approach.\nWe also compare with **JCI-FCI** [4] in the appendix as a constraint-based method for confounding. Further, we systematically analyzed the impact of confounders (Fig. 2a,b), nodes (Fig. 2c,d), interventions (Fig. 3), and robustness factors (sample size, noise, edge density, and non-contractive SEMs in the appendix) to identify where DCCD-CONF offers the greatest benefit.\n\n## Regarding Q2 and the second half of the comment W1\n\nWe thank the reviewer for this question. To directly estimate the individual effects of nonlinearity, cycles, and confounders on statistical performance (Q2), we conducted controlled experiments that systematically vary each factor while holding the others fixed. In addition to the baselines used in the main paper (LLC, NODAGS-Flow, DCDI), we include DAGMA [5], Bhattacharya et al. [6], and LiNGAM-MMI [7], each integrated with the JCI framework to incorporate interventional data.\n\n**Effects of nonlinearity.**\nWe vary the degree of nonlinearity in the data-generating process using\n$$X = (1 - \\beta)\\cdot (W^\\top X + Z) + \\beta\\cdot\\tanh(W^\\top X + Z),$$ where $\\beta$ controls nonlinearity.\n\nDirected edge recovery (Error metric: SHD)\n\n| **Beta**               | **0** | **0.25** | **0.5** | **0.75** | **1** |\n| ---------------------- | ----- | -------- | ------- | -------- | ----- |\n| DCCD-CONF              | 4.8   | 4.4      | 3.3     | 3.1      | 2.6   |\n| JCI-Bhattacharya et al | 13.4  | 20.6     | 20.4    | 21       | 21.1  |\n| JCI-Dagma              | 12.3  | 15.6     | 12.6    | 16       | 14.7  |\n| LLC                    | 3.2   | 5.6      | 6.6     | 8.2      | 7.7   |\n\nBidirectional edge recovery (Error metric: F1 score)\n\n| Beta                   | 0     | 0.25  | 0.5   | 0.75  | 1     |\n| ---------------------- | ----- | ----- | ----- | ----- | ----- |\n| DCCD-CONF              | 0.752 | 0.804 | 0.83  | 0.87  | 0.88  |\n| JCI-Bhattacharya et al | 0.249 | 0.29  | 0.213 | 0.202 | 0.174 |\n| LLC                    | 0.164 | 0.164 | 0.163 | 0.164 | 0.164 |\n\n_Interpretation:_ As nonlinearity increases, DCCD-CONF’s error decreases, while linear baselines (e.g., Bhattacharya et al.) degrade sharply, demonstrating the benefit of DCCD-CONF’s nonlinear modeling capacity.\n\n---\n\n**Effects of cycles.**\nWe vary the number of cycles while fixing $\\beta = 1.0$.\n\nDirected edge recovery (Error metric: SHD)\n\n| **Num. Cycles**        | **0** | **2** | **4** | **6** | **8** |\n| ---------------------- | ----- | ----- | ----- | ----- | ----- |\n| DCCD-CONF              | 0.6   | 2.5   | 2.7   | 4     | 4.1   |\n| JCI-Bhattacharya et al | 5.7   | 9.6   | 10.3  | 14.6  | 20    |\n| JCI-Dagma              | 12.1  | 12.3  | 12.6  | 21.5  | 27.2  |\n| LLC                    | 2.8   | 4.3   | 5.8   | 13.7  | 17.6  |\n\nBidirectional edge recovery (Error metric: F1 score)\n\n| Num. Cycles            | 0     | 2     | 4     | 6     | 8     |\n| ---------------------- | ----- | ----- | ----- | ----- | ----- |\n| DCCD-CONF              | 0.93  | 0.84  | 0.82  | 0.81  | 0.79  |\n| JCI-Bhattacharya et al | 0.25  | 0.18  | 0.2   | 0.23  | 0.18  |\n| LLC                    | 0.184 | 0.161 | 0.173 | 0.167 | 0.151 |\n\n_Interpretation:_ DCCD-CONF maintains low error even with many cycles, while acyclic baselines (DAGMA, Bhattacharya et al.) degrade significantly, confirming that explicit cycle modeling is crucial.\n\n---\n\n**Effects of confounders.**\n\nWe vary the number of confounders by adjusting the noise covariance structure.\n\nDirected edge recovery (Error metric: SHD)\n\n| **Num. Confounders** | **2** | **4** | **6** | **8** |\n| -------------------- | ----- | ----- | ----- | ----- |\n| DCCD-CONF            | 1.2   | 1.7   | 3.2   | 4.7   |\n| JCI-Bhattacharya     | 17.1  | 18    | 15.3  | 21.6  |\n| JCI-Dagma            | 11    | 13.3  | 15.6  | 16    |\n| LLC                  | 10.2  | 11.8  | 12    | 9     |\n\nBidirectional edge recovery (Error metric: F1 score)\n\n| Num. Confounders | 2    | 4    | 6     | 8    |\n| ---------------- | ---- | ---- | ----- | ---- |\n| DCCD-CONF        | 0.93 | 0.9  | 0.89  | 0.89 |\n| JCI-Bhattacharya | 0.21 | 0.19 | 0.17  | 0.13 |\n| LLC              | 0.1  | 0.21 | 0.289 | 0.36 |\n\n_Interpretation:_ DCCD-CONF is substantially more robust to increasing confounders compared to non-confounder-aware methods (DAGMA, Bhattacharya et al.), demonstrating the importance of jointly modeling latent confounding.\n\n---\n\n**Comparison with LiNGAM-MMI.**\n\nDue to its computational cost, we compare LiNGAM-MMI (as it iterates over all permutations of the observations) on 5-node graphs:\n\nDirected edge recovery (Error metric: SHD)\n\n| **Num. Confounders** | **1** | **2** | **3** | **4** |\n| -------------------- | ----- | ----- | ----- | ----- |\n| DCCD-CONF            | 0     | 0.2   | 0.2   | 0.2   |\n| LiNGAM-MMI           | 4     | 7     | 10    | 8.5   |\n\n_Interpretation:_ Even in small graphs, DCCD-CONF outperforms LiNGAM-MMI, highlighting that differentiable modeling scales more effectively with confounding.\n\n---\n\nThese experiments confirm that the performance gains of DCCD-CONF arise from its ability to jointly handle nonlinearity, cycles, and confounders. We will integrate these results into the revision along with a more rigorous ablation.\n## Regarding Q3\nWe thank the reviewer for this question. We compare the performance of DCCD-CONF and several state-of-the-art causal discovery methods on the protein signaling dataset of Sachs et al. [8], a widely used benchmark in causal discovery. We use the consensus network from [8] as ground truth. We selected this dataset instead of DREAM4 because DREAM4 provides only a single steady-state measurement per intervention, which is not conducive to differentiable likelihood-based methods like DCCD-CONF that require multiple samples per interventional setting.\n\n| **Method**   | **SHD** |\n| ------------ | ------- |\n| DCCD-CONF    | 19      |\n| DAG-GNN [9]  | 19      |\n| DAGMA        | 21      |\n| NOTEARS [10] | 22      |\n\nAs shown above, DCCD-CONF matches the best-performing model on this standard benchmark, demonstrating its competitiveness on real-world datasets commonly used in causal discovery.\n\n## Regarding comment W2\n\nWe thank the reviewer for this comment. We agree that Perturb-CITE-seq is not yet a widely adopted benchmark in the causal discovery community, which makes it difficult to evaluate methods solely based on this dataset. To address this concern, we have (i) added a comparison on the well-established protein signaling dataset (Sachs et al. [8]) in the previous section, and (ii) extended our Perturb-CITE-seq results by including DCD-FG, Bicycle, and NODAGS-Flow for a fairer evaluation.\n\nWe report the I-MAE scores, which measure the reconstruction error in unseen interventional conditions for genes not directly targeted:\n\n| **Condition** | **Cocult** | **IFN** | **IL-4** |\n| ------------- | ---------- | ------- | -------- |\n| DCCD-CONF     | 0.761      | 0.837   | 0.779    |\n| Bicycle       | 0.735      | 0.890   | 0.782    |\n| DCD-FG        | 0.774      | 0.883   | 0.845    |\n| NODAGS-Flow   | 0.762      | 0.861   | 0.847    |\n\nAs seen from the table, DCCD-CONF is competitive with state-of-the-art methods on this challenging dataset. We will include these results in the final version of the manuscript. We also attempted to include DAG-BLR in our comparison but were unable to obtain a workable implementation within the rebuttal timeframe; we will add this comparison in the camera-ready version.\n\n## Regarding comment W3\n\nWe thank the reviewer for pointing out these relevant works. We will include a discussion of these methods in the related work section of the revised manuscript to better situate our contribution within the broader literature.\n\nAs part of our rebuttal, we have added preliminary comparisons with two of the suggested methods: the differentiable approach of Bhattacharya et al. [6] and LiNGAM-MMI [7]. We did not include doubly-robust inference methods, as they primarily focus on causal effect estimation given (partial) causal graphs, which differs from our goal of full causal structure discovery.\n\nFor the remaining methods, including the integer-programming approaches, we were unable to obtain workable code within the rebuttal timeframe. However, we will make every effort to incorporate these comparisons into the final version of the manuscript.\n\n**References**\n\n[1] Hyttinen et al., JMLR 2012\n\n[2] Sethuraman et al., AISTATS 2023\n\n[3] Brouillard et al., NeurIPS 2020\n\n[4] Mooij et al., JMLR 2020\n\n[5] Bello et al., NeurIPS 2022\n\n[6] Bhattacharya et al., AISTATS 2021\n\n[7] Suzuki & Yang, ISIT 2024\n\n[8] Sachs et al., Science 2005\n\n[9] Yu et al., ICML 2019\n\n[10] Zheng et al., NeurIPS 2018"}, {"review_id": "5bQLACt0Um", "reviewer": "Reviewer_9ULc", "summary": "The paper provides an approach to learn cyclic causal graphs in the presence of (unmeasured) confounders. It defines cyclic causal graph using structural equations models with Gaussian noise. \nIt proposes a special parametrization of the structural equations that enables the differentiable learning of the structural equation as well as the Gaussian covariance structure (the confounders).\nThrough a combination of numerical approximation and mathematical tricks, the authors are able to train their model and show that it works competitively on simulated data.\nThe authors finally validate their method on real data from gene expression.", "questions": "- Can you clarify in section 2 what formulations and equations are standard/from prior work? Just to avoid confusion between prior work and your contributions. \n- Can you reconcile your theorem 2 (which has no condition on the family {I_k} and the sentence: “If the intervention set consists of all single-node interventions, Hyttinen et al. [26] showed that the ground truth DMG can be uniquely recovered in the linear setting.” Is it because “ground truth DMG” here is more restrictive than I-markov class? Isn’t by definition an I-markov class composed oof graphs that are indistinguishable from data from interventions I?\n- Lines 230:  “First, we optimizeˆS(B) with respect to the neural network parameters θ and the graph structure parameters B.” — the notations are a bit inconsistent: either you include both parameters S(B, theta) or you include them in the argmax, but not one in each —  or you need to rephrase the sentence.\n- Eq (13): consider writing explicitly L(Ik = \\emptyset)\n- Line 281: typo arrises", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strength:\n- I found the paper to be clear and well written\n- The ideas are simple/efficient and build elegantly on previous works (notably regarding the numerical tricks and approximations)\n- The method seems to work well on simulations, and the authors made an effort to evaluate on real data.\n\nWeakness\n- The results seem to rely on the parametrization of eq (7). Can you explain/justify/motivate the practical properties of such a parametrization. The contractivity does not bother me (it seems important to have an equilibrium), but the pseudo additivity of x and z seems a bit at odds with the claim of non-linearity. It seems that such a parametrization removes all nonlinear relations between x and z. Is this limiting? You mentioned “multiple such implicit blocks can be stacked.” I am not sure if this answers the question, but this statement was vague and should be clarified or removed.\n- The training of DCCD-CONF seems to rely on many approximations and tricks with potentially unstable training. How robust is the training to the hyperparameters (e.g. lambda)? By experience, those differentiable causal discovery methods are brittle with respect to it.\n- Lack of computational power comparison (e.g. training time).\n- Thank you for providing a real-data evaluation. However, it does not contain graph ground-truth and I am not convinced by the NLL. The NLL measures the quality of the generative model, and I would not be surprised that a great generative model with a wrong causal graph could beat a worse generative model with correct causal graph. Since you have “only” 61 genes, can you analyze the output and show how many cycles are learned, do they make sense biologically? Are we finding expected directionality between major genes (e.g. transcription factors?).", "was_revised": true, "final_justification": "The authors addressed my questions. The paper is well written and propose interesting new methodology. I have my own disagreements with the field of differentiable causal discovery field, but among this field, I think the paper is good.", "rebuttal": "We thank the reviewer for their thorough and constructive feedback, as well as for recognizing the clarity, simplicity, and practical value of our approach. We appreciate the positive assessment of our contributions and the helpful suggestions regarding the parametrization in Eq. (7), training robustness, computational considerations, and the real-data analysis. We address each of these points in detail below.\n## Regarding weakness pt 1\n\nWe thank the reviewer for raising this point. Our parameterization in Eq. (7) follows the structured implicit flow formulation of Lu et al. [1], chosen to balance expressiveness and tractability. Specifically, the pseudo-additive form\n$$F(X, Z) = -g_x(X) + g_z(Z) + Z$$\nensures a contractive mapping, which guarantees a unique equilibrium and invertibility, while enabling efficient computation of the log-determinant via the power series expansion (Eq. 11).\n\nIn practice, this structure offers several advantages:\n- **Stable and efficient inference:** The contractive map leads to fast convergence of root-finding (typically <10 iterations with Broyden’s method).\n- **Tractable likelihood estimation:** As detailed in section 3.3.1 in the manuscript, the additive decomposition of $g_x$ and $g_z$ simplifies Jacobian factorization, allowing unbiased log-determinant estimation without costly matrix operations. The main benefit comes from the fact that the power series expansion of the log-determinant of the Jacobian converges given our parameterization. \n- **Expressive modeling:** While a single block already captures nonlinear dependencies within $X$ and between $X$ and $Z$, stacking multiple implicit blocks (as suggested in [1]) increases expressiveness while preserving stability. Single layer of the implicit layer above already models nonlinear interaction between $X$ and $Z$ through the function $g_z$. Furthermore, note that the mapping between X and Z is defined implicitly as\n$$X = (id + g_x)^{-1}\\circ (id + g_z)(Z),$$where $g_x$ and $g_z$ are universal function approximators. By stacking multiple implicit normalizing flow layers, our framework can capture highly nonlinear dependencies between X and Z.\n\nWe will clarify these points in the revised manuscript and include runtime statistics for root-finding and log-determinant estimation to further support this choice.\n\n## Regarding weakness pt 2\n\nWe thank the reviewer for this insightful comment. We agree that differentiable causal discovery methods can be sensitive to hyperparameters, and thus we evaluated the robustness of DCCD-CONF with respect to its two main sparsity parameters: (i) $\\lambda$, which controls directed edge sparsity, and (ii) $\\rho$, which controls bidirected edge sparsity.\n\n  Using the same experimental setup described in Section 4, we report the average SHD for directed edge recovery (10 trials):\n\n**$\\rho$ fixed to $10^{-2}$**\n\n| **λ** | **1E-5** | **1E-4** | **1E-3** | **1E-2** | **1E-1** | **1** |\n| ----- | -------- | -------- | -------- | -------- | -------- | ----- |\n| SHD   | 4.2      | 3.6      | 0.4      | 0.1      | 16.8     | 19.7  |\n\n**$\\lambda$ fixed to $10^{-2}$**\n\n| **ρ** | **1E-4** | **1E-3** | **1E-2** | **1E-1** | **1** | **10** |\n| ----- | -------- | -------- | -------- | -------- | ----- | ------ |\n| SHD   | 8.9      | 7.5      | 0.1      | 0.8      | 2.2   | 2.4    |\n\nThese results indicate that DCCD-CONF is **stable across a broad range of hyperparameter values**, with an optimal region around $(\\lambda, \\rho) \\approx (10^{-2}, 10^{-2})$. Outside this range (e.g., $\\lambda > 10^{-1}$), performance degrades gracefully rather than collapsing, suggesting that the method is not excessively brittle.\n\nWe will add these results (and corresponding plots in the appendix) to clarify that DCCD-CONF can be tuned reliably in practice.\n\n## Regarding weakness pt 3\n\nWe thank the reviewer for this comment. We will include a detailed computation time analysis of the proposed method. Here, we provide preliminary comparison of the training time for graphs of size $d=10$ with 1100 training samples for DCCD-CONF and the baselines. \n\n| Model                  | Training time (sec) |\n| ---------------------- | ------------------- |\n| DCCD-CONF              | 50.31               |\n| LLC                    | 0.095               |\n| NODAGS-Flow            | 11.45               |\n| JCI-Bhattacharya et al | 73.859              |\n| JCI-Dagma              | 120.25              |\n| Lingam-MMI             | DNF                 |\n\n## Regarding weakness pt 4\n\nWe thank the reviewer for this valuable suggestion. We agree that NLL alone cannot validate causal correctness and thus the learned graph structure. For the Perturb-CITE-seq. dataset, DCCD-CONF identified **38 feedback cycles**. While we do not claim to provide definitive biological validation, this number is consistent with prior work showing that gene regulatory networks are rich in feedback loops (e.g., [2]).\n\nWe further note that gene CKS1B has the most outward going edge, while gene B2M has the most incoming edges. We will include a detailed list of top-ranked edges in the appendix to facilitate expert inspection.\n\nAlthough a full biological validation is outside the scope of this work, we believe these results provide supporting evidence that DCCD-CONF recovers biologically plausible structure, beyond what is reflected in the NLL. We will explicitly note this in the revised manuscript and highlight this as an important direction for future collaboration with domain experts.\n\n## Regarding question 1\n\nWe thank the reviewer for their suggestion. To clarify, Equations (1)–(4) and the accompanying definitions of structural equations, exogenous noise, and surgical interventions follow standard formulations from prior work on causal graphs [3, 4, 5]. Our contributions in this section are to (i) extend these formulations to possibly cyclic directed mixed graphs, (ii) model hidden confounders through correlated exogenous noise, and (iii) enable nonlinear interactions between exogenous and endogenous variables, in contrast to the additive noise models that typically assume linearity. Importantly, our formulation supports these extensions while maintaining efficient likelihood computation. We will revise the manuscript to make these distinctions explicit.\n\n## Regarding question 2\n\nWe thank the reviewer for this question. It is correct that our theorem does not guarantee perfect recovery of the ground-truth graph. Rather, it establishes that, for any given family of interventions, the recovered graph lies in the same interventional Markov equivalence class as the ground truth.\n\nIn the linear setting, however, Hyttinen et al. [6] showed that the edge weights can be identified by solving a collection of linear systems—one per interventional setting. They further proved that if the intervention targets satisfy the _pair condition_ (i.e., for every pair $(X_i, X_u) \\in V \\times V$, there exists an intervention where $X_i$ is intervened upon while $X_u$ is passively observed), then the interventional Markov equivalence class collapses to the ground-truth graph. A simple example of such a family of interventions is the set of all single-node interventions.\n\nWe will clarify this distinction in the revision and note that investigating analogous conditions for the nonlinear setting is an important direction for future work.\n\n**Regarding minor comments and suggestions**\n\n- **Line 230 (notation):** We thank the reviewer for pointing this out. We will revise the sentence for consistency. Specifically, we change the notation to $S(B, \\theta)$ explicitly to avoid confusion.\n- **Eq. (13):** We will explicitly write $\\mathcal{L}(I_k = \\emptyset)$ as suggested, to improve clarity.\n- **Line 281 (typo):** We will correct the typo “arrises” to “arises.”\n\n**References**\n\n[1] Lu et al., ICLR 2021\n\n[2] Frangieh et al., Nature genetics 2021\n\n[3] Bollen, Structural equations with latent variables (1989)\n\n[4] Pearl, Causality (2009)\n\n[5] Peters and Buhlmann, Annals of Statistics (2014)\n\n[6] Hyttinen et al., JMLR (2012)"}, {"review_id": "BGVtCjswQK", "reviewer": "Reviewer_mFwQ", "summary": "The paper presents the DCCD-CONF framework for causal discovery under presence of latent confounders, learning non-linear cyclic causal relationships among variables using interventional data with Gaussian noise assumption. The work demonstrates its applicability through experiments on real-world and synthetic datasets.", "questions": "* Suggestions\n1. Perhaps you can state assumptions A.16 ~ 18 in the main paper, at least in high-level on what those represent? The main theoretical result (Theorem 2) is a crucial part of the results, but it is hard to tell what those assumptions imply, from the main paper.", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "* Strengths\n1. The paper addresses the limitations of prior work on causal discovery, i.e., nonlinear cyclicity of causal structures, presence of unmeasured confounders, and interventions, simultaneously.\n2. The work provides a guarantee that a graph induced from the maximum score of the proposed score function results in a Markov equivalence class of the ground truth graph, under certain assumptions (Theorem 2).\n3. Experiments performed over real-world data, and synthetic data testing across various qualitative/quantitative metrics add concrete empirical support to the results.\n\n* Weaknesses\n1. If possible, it could be more useful to have a showcase of the resulting graph, semantically analyze some of the present edges generated by DCCD-CONF, and compare with baseline methods. For example, an edge A <-> B outputted by DCCD-CONF makes sense under certain domain knowledge, but some other methods output either A -> B, A <- B, or none, which does not necessarily match the domain knowledge.", "was_revised": true, "final_justification": "I appreciate that the authors provided detailed explanation regarding other concerns (e.g., hyperparameter settings and comparing performance with other methods), in addition to certain assumptions that I was looking for its clarification.", "rebuttal": "We thank the reviewer for their positive assessment of our contributions, particularly regarding our unified treatment of nonlinear cyclic structures, latent confounders, and interventional data, as well as the theoretical guarantees and extensive experiments. Below, we address the specific points raised by the reviewer. \n\n## Regarding weakness\n\nWe thank the reviewer for this helpful suggestion. We agree that showcasing the learned graph and analyzing edges in relation to domain knowledge would strengthen the paper. In the revised manuscript, we will add an appendix with a semantic comparison of representative edges recovered by DCCD-CONF against baseline methods for the gene perturbation dataset, as well as a similar analysis on synthetic data to better illustrate the differences between DCCD-CONF and the baselines.\n\n## Regarding the assumption\n\nWe thank the reviewer for their comment. We agree that providing a high-level explanation of the technical assumptions will improve the readability of the main theoretical result. We will add such a summary in the main text.\n- **Assumption A.13 (Sufficient capacity):** Ensures the data-generating distribution lies within the model class.\n- **Assumption A.14 (Interventional faithfulness):** Guarantees that any independence in the data corresponds to a $\\sigma$-separation in the ground truth graph.\n- **Assumption A.15 (Finite differential entropy):** Prevents the score function from diverging to infinity.\n\nFrom these assumptions, we derive the following lemmas:\n- **Lemma A.16 (Finite score):** The score remains bounded (follows from A.15).\n- **Lemma A.17:** The score difference can be expressed as an infimum of KL divergence plus a sparsity penalty.\n- **Lemma A.18:** Any distribution outside the interventional Markov equivalence class has strictly positive KL divergence from the ground truth.\n\nThese results form the basis of Theorem 2."}, {"review_id": "pizqbVAkQ9", "reviewer": "Reviewer_Ctrj", "summary": "This paper presents DCCD-CONF, a novel differentiable framework for learning nonlinear cyclic causal graphs under unmeasured confounding, leveraging both observational and interventional data. The approach models causal mechanisms using neural networks with implicit flows, allowing cycles and Gaussian noise-based confounder modeling. It introduces a score function based on data likelihood regularized for sparsity, with a two-step optimization alternating between causal graph structure and confounder distribution estimation.", "questions": "Please answer my concerns in W1 and W3.", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "- S1: the problem of cyclic causal discovery with latent confounders is underexplored. DCCD-CONF explicitly allows cycles and latent confounders.\n- S2: the authors provide consistency guarantees and show that the method recovers the I-Markov equivalence class under reasonable assumptions.\n- S3: The method outperforms strong baselines such as NODAGS-Flow, LLC, and DCDI across synthetic and real-world (gene perturbation) datasets, particularly in nonlinear and confounded settings.\n\n- W1: the characteristics of a \"cyclic\" graph are vague and not well clarified in the formualtion in the sense that it’s unclear what \"the equilibrium state\" refers to in the context of equations (1) and (2). an illustrative example would be helpful.\n- W2: line 74: the authors claim \"A few recent approaches, such as Bhattacharya et al. [33], have explored continuous optimization frameworks using differentiable constraints, though these methods are currently limited to linear settings,\" which is inaccurate. For instance, N-ADMG [1] provides a differentiable method to learn an ADMG from nonlinear data.\n- W3: line 113: the authors claim \"This formulation generalizes prior work by allowing cycles, extending both nonlinear cyclic models without unmeasured confounders that assume independent noise terms [ 27 ], and acyclic models without confounders [37].\" but I cannot identify non-trivial differences from the one formulated in Bhattacharya et al [33] (which is not even mentioned in the text). if I understand correctly, the correlated noise term (i.e., $(\\Sigma_Z)_{ij}\\neq 0$) is the same as what is defined in [33], except for replacing the additive noise in SCM with a general function $F_i(pa_i, z_i)$. Then, in the concrete formulation in Eqn (7), it appears that the noise term is still additive.\n- W4: given the differentiable nature of the proposed method, a more reasonable baseline would be JCI+Bhattacharya et al [33] or JCI+N-ADMG, instead of JCI+FCI. \n- W5: the runtime is not reported.\n\n---\n[1] Causal reasoning in the presence of latent confounders via neural ADMG learning", "was_revised": true, "final_justification": "The authors have addressed most of my concerns and confusion. Therefore, I raised my rating.", "rebuttal": "We thank the reviewer for the thoughtful and constructive feedback. We are pleased that the reviewer recognized the novelty of addressing **cyclic causal discovery with latent confounders (S1)**, the theoretical consistency guarantees (S2), and the strong empirical performance of DCCD-CONF across synthetic and real-world datasets (S3). We address each of the raised concerns (W1–W5) below and clarify how we will revise the manuscript accordingly.\n\n**Regarding concern W1**\n\nIn the acyclic case, we can establish a causal ordering $\\pi$ that allows us to generate samples iteratively: $X_{\\pi_1} = F_{\\pi_1}(Z_{\\pi_1})$, $X_{\\pi_2} = F_{\\pi_2}(X_{pa_\\mathcal{G}(\\pi_2)}, Z_{\\pi_2})$, …, $X_{\\pi_d} = F_{\\pi_d}(X_{pa_\\mathcal{G}(\\pi_d)}, Z_{\\pi_d})$, where $\\mathrm{pa}_\\mathcal{G}(\\pi_i) \\subseteq \\pi_{<i}$.   \n\nIn contrast, for cyclic graphs such an ordering does not exist. Here, the causal model represents a process that evolves over time and is assumed to reach an **equilibrium state** before observation. Concretely, starting from an initial value $X(0)$ and sampling one vector of disturbances $Z$, we iteratively update:\n$$X(t) = F\\big(X(t-1), Z\\big),$$\nuntil the system converges to a fixed point $X^\\ast$. Equation (2) in the paper corresponds exactly to this fixed-point relation. A sufficient condition for convergence is that $F$ is contractive in $X$ (as in prior work on cyclic SEMs, e.g., [1], [2]).\n\n**Illustrative example:** Consider two variables with mutual regulation:\n$$X_1 = \\tanh(0.7 X_2 + Z_1), \\quad X_2 = \\tanh(0.5 X_1 + Z_2).$$\nStarting from $X(0) = (0, 0)$ and iterating the update converges to a unique fixed point $(X_1^\\star, X_2^\\star)$, which represents the observed equilibrium. For instance, if $(Z_1, Z_2) = (-0.67, 1.11)$, then the fixed point (also the observation) is the given by $(X_1^\\ast, X_2^\\ast) = (-0.12, 0.78)$. \n\nThus, “equilibrium” refers to the steady-state solution of this iterative process, at which point the time dependence is removed and the model matches equation (2).\n\n**Regarding concern W2**\n\nWe thank the reviewer for pointing out this oversight. We agree that our original statement was inaccurate, as N-ADMG [3] indeed provides a differentiable approach for learning ADMGs from nonlinear data. We will revise the manuscript to reflect this and include N-ADMG both in the related work discussion and as a baseline in our experiments. We also note that our method differs from N-ADMG in that we focus on cyclic graphs (directed mixed graphs) rather than acyclic mixed graphs, which we will clarify in the revision.\n\n**Regarding concern W3**\n\nWe thank the reviewer for this helpful comment. The statement in line 113 was intended to clarify that our framework extends the problem formulations of [2] (cyclic graphs without confounding) and [4] (acyclic graphs without confounding), while also handling interventional data. We agree that Bhattacharya et al. [5] model confounding in a similar way (via correlated noise, i.e., $(\\Sigma_Z)_{ij}\\neq 0$), but their approach is limited to **linear acyclic models** and does not explicitly incorporate interventional data. Our work generalizes this by supporting **nonlinear mechanisms**, **feedback cycles**, and interventions, making it a non-trivial extension.\n\nMoreover, while Eq. (7) is expressed in an additive form, the mapping between X and Z is defined implicitly as\n$$X = (id + g_x)^{-1}\\circ (id + g_z)(Z),$$\nwhere $g_x$ and $g_z$ are universal function approximators. By stacking multiple implicit normalizing flow layers, our framework can capture highly nonlinear dependencies between X and Z, far beyond the linear additive structure of Bhattacharya et al. [5].\n\n**Regarding concern W4**\n\nWe thank the reviewer for this valuable suggestion. We agree that JCI+Bhattacharya et al. [5] and JCI+N-ADMG are appropriate baselines for comparison. We will update our empirical study to include both methods in the revised manuscript. Due to the lack of publicly available, executable code for N-ADMG within the rebuttal timeframe, we instead provide a preliminary comparison with JCI+Bhattacharya et al. and will incorporate N-ADMG in the camera-ready version.\n\nPerformance comparison between DCCD-CONF and JCI+Battacharya et al as a function of number of confounders. The synthetic graphs consists of ER random graphs with $d=10$ nodes. The SEM is nonlinear with the number of confounders specified in the top row. The training data consists of observational data as well as single node interventions across all the nodes in the graph. \n\nDirected edge recovery (Error metric: SHD)\n\n| Num. Confounders       | 2    | 4    | 6    | 8        |\n| ---------------------- | ---- | ---- | ---- | -------- |\n| DCCD-CONF              | 1.25 | 1.76 | 3.2  | 4.78     |\n| JCI-Bhattacharya et al | 17.1 | 18   | 15.3 | 21.6 |\n\nBidirectional edge recovery (Error metric: F1 score)\n\n| Num. Confounders       | 2    | 4    | 6    | 8    |\n| ---------------------- | ---- | ---- | ---- | ---- |\n| DCCD-CONF              | 0.93 | 0.9  | 0.89 | 0.89 |\n| JCI-Bhattacharya et al | 0.21 | 0.19 | 0.17 | 0.13 |\n\nPerformance comparison between DCCD-CONF and JCI+Battacharya et al as a function of number of cycles in the graph. The synthetic graphs consists of ER random graphs with $d=10$ nodes. The SEM is nonlinear with the number of confounders specified in the top row. The training data consists of observational data as well as single node interventions across all the nodes in the graph. \n\nDirected edge recovery (Error metric: SHD)\n\n| Num. Cycles            | 0   | 2   | 4    | 6    | 8   |\n| ---------------------- | --- | --- | ---- | ---- | --- |\n| DCCD-CONF              | 0.6 | 2.5 | 2.7  | 4    | 4.1 |\n| JCI-Bhattacharya et al | 5.7 | 9.6 | 10.3 | 14.6 | 20  |\n\nBidirectional edge recovery (Error metric: F1 score)\n\n| Num. Cycles            | 0    | 2    | 4    | 6    | 8    |\n| ---------------------- | ---- | ---- | ---- | ---- | ---- |\n| DCCD-CONF              | 0.93 | 0.84 | 0.82 | 0.81 | 0.79 |\n| JCI-Bhattacharya et al | 0.25 | 0.18 | 0.2  | 0.23 | 0.18 |\n\n**Regarding concern W5**\n\nWe thank the reviewer for this comment. We will include a detailed computation time analysis of the proposed method. Here, we provide preliminary comparison of the training time for graphs of size $d=10$ with 1100 training samples for DCCD-CONF and the baselines. \n\n| Model                  | Training time (sec) |\n| ---------------------- | ------------------- |\n| DCCD-CONF              | 50.31               |\n| LLC                    | 0.095               |\n| NODAGS-Flow            | 11.45               |\n| JCI-Bhattacharya et al | 73.859              |\n| JCI-Dagma              | 120.25              |\n| Lingam-MMI             | DNF                 |\n\n\n**References**: \n\n[1] Bongers, Stephan, et al. \"Foundations of structural causal models with cycles and latent variables.\" _The Annals of Statistics_ 49.5 (2021): 2885-2915.\n\n[2] Sethuraman, Muralikrishnna G., et al. \"Nodags-flow: Nonlinear cyclic causal structure learning.\" _International Conference on Artificial Intelligence and Statistics_. PMLR, 2023.\n\n[3] Ashman, Matthew, et al. \"Causal reasoning in the presence of latent confounders via neural ADMG learning.\" _arXiv preprint arXiv:2303.12703_ (2023).\n\n[4] Brouillard, Philippe, et al. \"Differentiable causal discovery from interventional data.\" _Advances in Neural Information Processing Systems_ 33 (2020): 21865-21877.\n\n[5] Bhattacharya, Rohit, et al. \"Differentiable causal discovery under unmeasured confounding.\" _International Conference on Artificial Intelligence and Statistics_. PMLR, 2021."}]

Author rebuttal:

Meta-review: **Summary.** The paper introduces a differentiable causal structure learning algorithm for learning non-linear (NL), cyclic (Cyc) causal models with unmeasured confounders (UC); the algorithm also incorporates interventional data (Int). They validate their algorithm on both simulated and real data, showing strong performance against baselines such as NODAGS-Flow (NL+Cyc), LLC (Cyc+UC), and DCDI (NL+Int). The appendix also includes JCI-FCI (NL+UC+Int)

**Strengths.** The paper is well-written and builds elegantly on previous work. The reviewers appreciate that the method simultaneously handles many issues (nonlinearity, cyclicity, causal insufficiency, and interventions) and found the experimental results to be convincing.

**Weakenesses.** 
1. *Ablation* One of the main weakness pointed out by the reviewers (e.g., LryH and Ctrj) was that the experiments were somewhat poorly designed, not performing a proper ablation over the four relevant aspects of the approach (NL, Cyc, UC, Int) by considering methods that leave out only one of the aspects (e.g., the experiments did not contain any method of the form NL+Cyc+UC). In the rebuttal, the authors added three new baselines, using the Joint Causal Inference (JCI) approach to adapt existing methods to the interventional setting. In particular, the authors ran JCI-DAGMA (NL+Int), JCI-Bhattacharya et al. (NL+UC+Int), and JCI-LINGAM-MMI (UC+Int). Overall, the reviewers were satisfied with these additional experiment, though I still believe that the experiment design could be more carefully executed.
2. *Unfamiliar real dataset* Reviewer LryH points out that the Perturb-CITE-seq dataset is not commonly used in causal structure learning, and thus it is difficult to evaluate the model's performance. In the rebuttal, the authors address this by including the results from other methods which have used this dataset, including Bicycle and DCD-FG, and also add results for the more commonly used Sachs dataset. This additional context will substantially help readers interpret the results.
3. *Other* The reviewers raised a variety of other minor weaknesses, including (a) missing citations to other methods for causal structure learning under unobserved confounding (Reviewer LryH), (b) the stableness of training wrt hyperparameters (Reviewer 9ULc), (c) missing comparisons on computational power (Reviewers 9ULc and Ctrj), (d) the need to give a high-level explanation of technical assumptions (Reviewer mFwQ), and (e) the need to better discuss the contractivity condition and the pseudo-additivity parameterization for cyclic SCMs (Reviewers 9ULc and Ctrj). The authors addressed these concerns in the rebuttal: for (a), they will include a discussion of the more methods in their revision, for (b), they demonstrated the performance of their method across a range of hyperparameters, for (c), they added the comparison of computational power, and for (d) and (e), they provided high-level explanations which will be added to the text.

**Conclusion.** Overall, the paper received solid ratings and the concerns of the reviewers were thoroughly addressed. It looks like the revision will involve several changes in terms of both presentation and results, but these changes are not so major as to disqualify the paper. The paper addresses an important gap in the field (a scalable method combining all four aspects mentioned), it is technically solid, and it is well-written. Thus, I suggest accepting the paper.
Decision: Accept (spotlight)