Differentiable Cyclic Causal Discovery Under Unmeasured Confounders

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.

Regarding concern W1

In 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}$.

In 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: $X(t) = F\big(X(t-1), Z\big),$ until 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]).

Illustrative example: Consider two variables with mutual regulation: $X_1 = \tanh(0.7 X_2 + Z_1), \quad X_2 = \tanh(0.5 X_1 + Z_2).$ Starting 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)$.

Thus, “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).

Regarding concern W2

We 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.

Regarding concern W3

We 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.

Moreover, while Eq. (7) is expressed in an additive form, the mapping between X and Z is defined implicitly as $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, far beyond the linear additive structure of Bhattacharya et al. [5].

Regarding concern W4

We 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.

Performance 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.

Directed edge recovery (Error metric: SHD)

Bidirectional edge recovery (Error metric: F1 score)

Performance 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.

Directed edge recovery (Error metric: SHD)

Bidirectional edge recovery (Error metric: F1 score)

Regarding concern W5

We 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.

References:

[1] Bongers, Stephan, et al. "Foundations of structural causal models with cycles and latent variables." The Annals of Statistics 49.5 (2021): 2885-2915.

[2] Sethuraman, Muralikrishnna G., et al. "Nodags-flow: Nonlinear cyclic causal structure learning." International Conference on Artificial Intelligence and Statistics. PMLR, 2023.

[3] Ashman, Matthew, et al. "Causal reasoning in the presence of latent confounders via neural ADMG learning." arXiv preprint arXiv:2303.12703 (2023).

[4] Brouillard, Philippe, et al. "Differentiable causal discovery from interventional data." Advances in Neural Information Processing Systems 33 (2020): 21865-21877.

[5] Bhattacharya, Rohit, et al. "Differentiable causal discovery under unmeasured confounding." International Conference on Artificial Intelligence and Statistics. PMLR, 2021.

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.

Regarding weakness

We 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.

Regarding the assumption

We 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.

• Assumption A.13 (Sufficient capacity): Ensures the data-generating distribution lies within the model class.
• Assumption A.14 (Interventional faithfulness): Guarantees that any independence in the data corresponds to a $\sigma$-separation in the ground truth graph.
• Assumption A.15 (Finite differential entropy): Prevents the score function from diverging to infinity.

From these assumptions, we derive the following lemmas:

• Lemma A.16 (Finite score): The score remains bounded (follows from A.15).
• Lemma A.17: The score difference can be expressed as an infimum of KL divergence plus a sparsity penalty.
• Lemma A.18: Any distribution outside the interventional Markov equivalence class has strictly positive KL divergence from the ground truth.

These results form the basis of Theorem 2.

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.

Regarding weakness pt 1

We 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 $F(X, Z) = -g_x(X) + g_z(Z) + Z$ ensures 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).

In practice, this structure offers several advantages:

• Stable and efficient inference: The contractive map leads to fast convergence of root-finding (typically <10 iterations with Broyden’s method).
• 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.
• 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 $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.

We will clarify these points in the revised manuscript and include runtime statistics for root-finding and log-determinant estimation to further support this choice.

Regarding weakness pt 2

We 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.

Using the same experimental setup described in Section 4, we report the average SHD for directed edge recovery (10 trials):

$\rho$ fixed to $10^{-2}$

$\lambda$ fixed to $10^{-2}$

These 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.

We will add these results (and corresponding plots in the appendix) to clarify that DCCD-CONF can be tuned reliably in practice.

Regarding weakness pt 3

We 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.

Regarding weakness pt 4

We 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]).

We 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.

Although 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.

Regarding question 1

We 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.

Regarding question 2

We 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.

In 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.

We will clarify this distinction in the revision and note that investigating analogous conditions for the nonlinear setting is an important direction for future work.

Regarding minor comments and suggestions

• 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.
• Eq. (13): We will explicitly write $\mathcal{L}(I_k = \emptyset)$ as suggested, to improve clarity.
• Line 281 (typo): We will correct the typo “arrises” to “arises.”

References

[1] Lu et al., ICLR 2021

[2] Frangieh et al., Nature genetics 2021

[3] Bollen, Structural equations with latent variables (1989)

[4] Pearl, Causality (2009)

[5] Peters and Buhlmann, Annals of Statistics (2014)

[6] Hyttinen et al., JMLR (2012)

Dear Reviewer 9ULc,

Thank you for the feedback on our rebuttal, and we clarify below our hyperparameter selection procedure and plans for additional real-world biological validation.

For weakness pt 2 (hyperparameter tuning)

We did not cherry-pick the values of $\lambda$ and $\rho$. These values were selected via a grid search using a validation set composed of 2-node interventions over randomly chosen nodes, with the negative log-likelihood (NLL) as the selection criterion. The grid search was conducted on graphs with $d=10$ nodes. While these values may vary with graph size, one can approximate suitable hyperparameters by optimizing for the NLL on held-out interventional data, which does not require ground-truth graphs.

For weakness pt 4 (real-world biological validation)

We agree that real-world validation of causal discovery methods remains challenging. As future work, we plan to provide additional biological validation using Perturb-CITE-seq data. In response to reviewer LryH, we have also evaluated our method on the Sachs et al. protein-signaling dataset [1], which includes a consensus network treated as ground truth, and we report the SHD comparison with the baselines below:

The table shows that DCCD-CONF matches the best-performing baseline (DAG-GNN). We will include this evaluation in the revised manuscript along with a qualitative comparison of the learned graph against the consensus network.

References:

[1] Sachs et al., Science 2005

[2] Yu et al., ICML 2019

[3] Bello et al., NeurIPS 2022

[4] Zheng et al., NeurIPS 2018

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.

Regarding Q1 and the first half of comment W1

We 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:

• LLC [1]: handles cycles and confounding but is linear, highlighting the value of nonlinearity.
• NODAGS-Flow [2]: supports nonlinearity and cycles but no confounding, isolating the effect of latent confounders.
• DCDI [3]: handles interventions but assumes acyclicity and no confounding, representing the standard differentiable DAG-based approach. We 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.

Regarding Q2 and the second half of the comment W1

We 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.

Effects of nonlinearity. We vary the degree of nonlinearity in the data-generating process using $X = (1 - \beta)\cdot (W^\top X + Z) + \beta\cdot\tanh(W^\top X + Z),$ where $\beta$ controls nonlinearity.

Directed edge recovery (Error metric: SHD)

Bidirectional edge recovery (Error metric: F1 score)

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.

---

Effects of cycles. We vary the number of cycles while fixing $\beta = 1.0$.

Directed edge recovery (Error metric: SHD)

Bidirectional edge recovery (Error metric: F1 score)

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.

---

Effects of confounders.

We vary the number of confounders by adjusting the noise covariance structure.

Directed edge recovery (Error metric: SHD)

Bidirectional edge recovery (Error metric: F1 score)

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.

---

Comparison with LiNGAM-MMI.

Due to its computational cost, we compare LiNGAM-MMI (as it iterates over all permutations of the observations) on 5-node graphs:

Directed edge recovery (Error metric: SHD)

Interpretation: Even in small graphs, DCCD-CONF outperforms LiNGAM-MMI, highlighting that differentiable modeling scales more effectively with confounding.

---

These 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.

Regarding Q3

We 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.

As 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.

Regarding comment W2

We 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.

We report the I-MAE scores, which measure the reconstruction error in unseen interventional conditions for genes not directly targeted:

As 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.

Regarding comment W3

We 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.

As 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.

For 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.

References

[1] Hyttinen et al., JMLR 2012

[2] Sethuraman et al., AISTATS 2023

[3] Brouillard et al., NeurIPS 2020

[4] Mooij et al., JMLR 2020

[5] Bello et al., NeurIPS 2022

[6] Bhattacharya et al., AISTATS 2021

[7] Suzuki & Yang, ISIT 2024

[8] Sachs et al., Science 2005

[9] Yu et al., ICML 2019

[10] Zheng et al., NeurIPS 2018