We are grateful for your engagement and the constructive feedback on our submission.

Weaknesses

Identifiability. The author did not discuss what kind of DAGs the proposed model can identify. Since the VI procedure requires knowledge of the likelihood, it seems that the noise distribution has to be known to guarantee identifiability, which is a limitation that needs to be declared.

Thank you for highlighting this point. Indeed, identifiability depends on both the data-generating SEM (including the noise distribution) and the user-specified likelihood. For instance:

• If the likelihood is correctly specified (e.g., Gaussian likelihood applied to Gaussian noise), the true DAG (or its equivalence class) is theoretically identifiable. In our experiments (e.g., Figure 3), we observe this behavior as the posterior increasingly concentrates on the true structure as $n$ grows; and
• If the likelihood is misspecified (e.g., Gaussian likelihood applied to non-Gaussian noise), full identifiability may no longer hold. In this scenario, ProDAG will reflect the resulting uncertainty in its posterior.

Thus, ProDAG captures and reports the inherent uncertainty determined by the combination of the true SEM and the chosen likelihood. We will explicitly clarify this subtlety in the revised manuscript.

Some experiments are not sufficient. For example, the experiment for nonlinear data is limited to $d=10$, and those for non-Gaussian noise distributions are limited to $d=20$. Experiments on higher dimensions, e.g., $d=100$ would be helpful.

We added the following experiments with more nodes $p$ ($d$ in your notation) to address your concerns:

1. Nonlinear data with $p=20$ (double the previous dimension) in Table R3;
2. Gumbel noise data with $p=100$ in Table R4; and
3. Exponential noise data with $p=100$ in Table R5.

Due to NeurIPS image restrictions this year, we report a snapshot of the experiments at $n=100$.

Table R3 Larger nonlinear experiment with $p=20$.

Table R4 Larger Gumbel noise experiment with $p=100$.

Table R5 Larger exponential noise experiment with $p=100$.

The results align with our earlier findings, confirming that ProDAG scales to larger, more complex settings. Note that DiBS/DiBS+ ran out of memory in the big nonlinear experiment so are not included. Likewise, BOSS is not reported in Table R3 because it does not accommodate nonlinear models, and Tables R4 and R5 because it encountered errors.

Time comparison with other Bayesian methods is not provided.

We added a time comparison in Table R6 for linear and nonlinear DAGs (corresponding to the experiments in Figures 3 and 5 of the paper). These timings show that ProDAG is faster than BayesDAG and slower than DiBS. However, ProDAG's superior performance in accuracy and uncertainty quantification makes it a strong choice, even with the higher computational cost compared to DiBS.

Table R6 Run times (seconds) against Bayesian baselines.

Questions

Do you assume knowledge of the true noise distribution for experiments on non-Gaussian data?

No, the model is always fit with a Gaussian likelihood $x=W^\top x+\epsilon,\epsilon\sim N(0,\Sigma)$ independent of the true noise distribution in the simulated data.

Is ProDAG robust if the log-likelihood function is misspecified? For example, using a Gaussian likelihood for data generated by non-Gaussian noise distributions, or an equal variance likelihood for a heterogeneous non-Gaussian noise distribution?

Yes, as in our non-Gaussian experiments, when we deliberately apply the Gaussian likelihood to data generated with Gumbel or exponential noise, ProDAG's performance is roughly similar to the Gaussian noise case, indicating it is resilient to moderate likelihood misspecification.

How sensitive are the results to the choice of $\mathbb{P}$?

Any continuous prior on $\tilde{W}$ is pushed through the same projection $\tilde{W}\mapsto W$, so the induced prior is largely governed by the $\ell_1$ ball and acyclicity constraints, not by the exact density of $\mathbb{P}$. We therefore use an isotropic Gaussian, which is quite general. Having said that, some priors such as Laplace, which shrink towards zero, are not well-suited as they can combine with the $\ell_1$-ball constraint to induce double shrinkage.

I thank the author for the detailed reply, which addressed most of my concerns. Thus, I will raise my score to 4. However, I still have a follow-up question on the noise variance. Do you assume equal noise variance in the experiments? If that is the case, I would recommend that the author include experiments on column-standardized data [1] or varying variance, as recent works [1,2,3] have identified limitations of differentiable causal discovery methods on data with unequal variance. Clarifying this would further strengthen the paper.

[1] Reisach, A., Seiler, C., & Weichwald, S. (2021). Beware of the simulated dag! causal discovery benchmarks may be easy to game. Advances in Neural Information Processing Systems, 34, 27772-27784.

[2] Ng, I., Huang, B., & Zhang, K. (2024). Structure learning with continuous optimization: A sober look and beyond. In Causal Learning and Reasoning (pp. 71-105). PMLR.

[3] Berrevoets, J., Raymaekers, J., Van der Schaar, M., Verdonck, T., & Yao, R. (2025). Differentiable Causal Structure Learning with Identifiability by NOTIME. In International Conference on Artificial Intelligence and Statistics (pp. 3115-3123). PMLR.

Thank you for taking the time to engage with our work and for your thoughtful feedback.

Weaknesses

The method inherits standard limitations including cubic complexity, which limits its scalability (experiments only go to p=100). The use of variational inference does not guarantee good approximation to the true posterior. The non-convexity of the acyclicity constraint means optimization can only reach stationary points rather than global optima.

It is correct that non-convexity means we can only attain a local optimum, though this property is true of all differentiable DAG learning methods (e.g., DAGMA). Regarding scalability, we now include larger experiments (see our response below).

While generally well-written, some important experimental details are deferred to appendices rather than being in the main text, and the combination of multiple existing techniques (though novel in combination) makes it hard to work out which components drive the performance improvements.

With the extra page allowed in the camera ready version, we will move additional experimental details (e.g., data generation and hyperparameter tuning) from the appendix to the main text. Regarding the source of performance gains, please see our response below concerning component attribution and ablations.

The performance improvements over baselines are consistent but small, and the scalability limitations to networks with around 100 variables may limit practical applicability to larger real-world problems where the method's benefits would be most valuable.

We agree that practical applicability depends on scalability. As noted in our response below, we now include results on the Andes network with $p=223$ in Table R2 (see response to Reviewer VF9w), which demonstrate that ProDAG scales beyond $p=100$ while maintaining its performance advantages.

While the overall framework is novel, the individual technical components such as the acyclicity constraints, L1-ball projections, and variational inference) are all existing techniques, so the contribution is primarily in their creative combination rather than fundamental algorithmic innovations. Nonetheless there are some significant innovations in the work.

The individual tools (acylclicity characterizations, $\ell_1$ balls) are indeed known, but their composition to define a projected prior is, as you say, new. No prior work samples exact sparse DAGs by projecting continuous draws and then instantiates that construction in a variational framework. This step is the key innovation that drives ProDAG's gains.

Questions

Can you demonstrate the method on larger networks (p=200-500), or propose any algorithmic improvements that could significantly increase the scaleability of the approach?

We now include experiments on the $p=223$ Andes network from bnlearn (Table R2). On this larger, real‑world‑structured benchmark, ProDAG continues to delver excellent performance and represents one of the largest Bayesian DAG evaluations reported to date.

Can you determine which components are responsible for the performance improvements? Some ablations could help with this, such as comparing against a simpler baselines with the same acyclicity constraints without projected distributions.

We believe that DAGMA, a key baseline in our experiments, reflects the ablation you propose: it enforces the same acyclicity constraints but omits projected distributions. DAGMA can also be interpreted as a MAP version of ProDAG. The consistent improvements over DAGMA underscore the value of modeling full posteriors directly over DAGs, made possible by our projected variational framework. In practice, these gains trace to two components: (i) projected DAG distributions that place mass on sparse, exactly acyclic structures, and (ii) variational inference that leverages this distributional support to propagate uncertainty during learning.

Can you provide recommendations about the settings in which your approach should be chosen over existing methods, given the modest improvements and higher computational cost?

Choose ProDAG when one of these situations arise:

1. You need well-calibrated posterior uncertainty (e.g., decision making, active experimentation);
2. The number of observations is limited, so precise sparsity and exact zeros aid recovery; or
3. Exact acyclicity and hard sparsity matter for interpretability (e.g., policy, biology).

Can you provide more analyses as to when the variational approximation works well? This could include diagnostics or an experiment comparing against exact posteriors on a small synthetic problem.

While exact posterior comparisons via MCMC are nontrivial, we report Brier score, SHD, F1 score, and AUROC as direct assessments of the variational approximation’s structure learning capacity. To further stress-test this approximation, we added new higher-dimensional experiments on nonlinear and non-Gaussian data (Tables R3–R5), and the Andes network. In all cases, ProDAG maintains strong performance, demonstrating robustness of the approximation across a range of challenging settings.

Thank you for your thoughtful review and engagement with the paper.

Weaknesses

The following works also considers a continuous constraint optimization problem The authors should compare their methods to the following works qualitatively and use some of them as benchmarks.

Thank you for these references. We agree they relate to our work but target different settings and all are non-Bayesian. Brouillard et al. (2020), Faria et al. (2022) and Lopez et al. (2022) rely on interventional data, whereas ProDAG assumes purely observational data. Bhattacharya et al. (2021) also uses observational data but learns acyclic directed mixed graphs (graphs with bidirected edges), not DAGs. He et al. (2021) adds a residual‑independence regularizer to frequentist NOTEARS‑type models, but it is unclear how to embed this within ProDAG’s Bayesian projection prior. We will cite the papers and discuss the extension of ProDAG to interventional data as future work.

The authors should provide more details on the “path-following algorithm”.

We agree that the path-following algorithm (presented in Appendix B) would benefit from additional explanatory text to complement the algorithmic details. We will revise the appendix to include a more descriptive explanation of the algorithm’s steps.

The authors should provide some intuition on “implicit function theorem.” and “Karush-Kuhn- Tucker (KKT) conditions of the projection” in the main paper.

We agree. While the full details were placed in Appendix C due to space limitations, we will include descriptions in the main text in the final version.

The role of $l_1$ ball is not clear. The author should discuss how large or small $l_1$ ball affects the algorithm performance. Also it is not clear how $l_1$ is chosen.

The $\ell_1$ ball constraint is included to sample sparse graphs; a larger $\ell_1$ ball mean more edges while a smaller $\ell_1$ ball means fewer edges. In our numerical experience, the algorithm's speed is not particularly sensitive to the size of the ball, though smaller balls tend to converge slightly faster than larger balls. The size $\lambda$ of the ball is chosen using a validation set (see Appendix E).

It is not clear how the authors obtain $q_\theta(\tilde{W})$ with variational inference in their algorithm. I would suggest the authors discussing this in detail.

The distribution $q_\theta(\tilde{W})$ is the variational family chosen by the user (see the last paragraph of Section 5.1) such as a multivariate Gaussian, and its parameters $\theta$ are optimized via gradient descent to maximize the ELBO. The projected distribution over DAGs, $q_\theta(W)$, is implicitly defined by applying the projection to samples from $q_\theta(\tilde{W})$.

The authors should show more experimental results with real-world data. They might check bnrepository for this purpose.

Thank you for the reference to this great repository. We have now added Table R2, which reports results on the Andes network ($p=223$) from the repository, drawing $n=100$ samples from its DAG with our usual data generation pipeline. This large, real‑world‑structured benchmark shows that ProDAG performs strongly on very high-dimensional DAGs, achieving the lowest Brier score and SHD. While DAGMA scores well on F1 and AUROC, ProDAG provides superior calibration and structure recovery.

Table R2 Experiments on real-world Andes DAG with $p=223$ nodes and $s=338$ edges. We report 8 runs here due to time constraints and will extend to 10 in the final version, but results are already clearly differentiated.

To our knowledge, this experiment represents one of the largest Bayesian DAG evaluations reported to date. Aside from ProDAG, DAGMA was the only method to successfully run on this large dataset. The other Bayesian baselines encountered difficulties: DiBS and DiBS+ produced nonsensical results (Brier scores and SHDs in the tens of thousands), indicative of convergence issues, while BayesDAG ran out of memory (24 GB available). BOSS encountered crashes.

Figure 2 might be replaced with a table and the space can be used for describing other concepts in detail.

If space is a concern in the final version, we will move Figure 2 to the appendix to make room for additional explanations in the main text.

Questions

How do the author set the value for $\lambda$.

The parameter $\lambda$ is chosen by maximizing the expected likelihood on a validation set. See Appendix E for details.

The authors said that the projection (2) is unique. Why can there not be two different $\tilde{W}$ with equal distance from $W$?

Uniqueness means that for almost every draw of $\tilde{W}$ the projection (2) has one minimizer, so $\operatorname{pro}_\lambda(\tilde{W})$ is single‑valued for each fixed $\tilde{W}$. This result does not imply injectivity: different $\tilde{W}$ can still project to the same $W$, but no single $\tilde{W}$ has more than one optimal projection. We will clarify this point in the text.

The authors mentioned that step 1 projects the matrix $\tilde{W}$ onto the set of acyclic matrices through equation 4. How does the algorithm have access to $W$ to calculate the subtraction here?

The $W$ in (4) appearing inside $f_\mu(W;\tilde{W})=\frac{\mu}{2}\|\tilde{W}-W\|_F^2+h(W)$ is the current optimization variable. The optimizer repeatedly updates this variable while it computes the difference $\tilde{W}-W$ until $f$ is minimized. After convergence we reuse the same symbol to denote the minimizer, so the duplicate notation is purely stylistic and we never assume prior knowledge of the solution.

We appreciate your time and attention in reviewing our work.

Weaknesses

Regarding weaknesses, while the empirical evaluation shows strong performance, I would like to have seen some more insight into why ProDAG performs better than existing approaches (see questions below).

ProDAG places both the prior and variational posterior directly on DAG space, enforcing exact acyclicity through projection. Competing methods (i) operate in relaxed spaces with soft penalties that cannot guarantee acyclicity, or (ii) introduce discrete variables (e.g., permutations) that complicate optimization. See the discussion in Section 2. By pushing irrelevant edges to exact zeros while keeping the ELBO differentiable, ProDAG yields sharper uncertainty estimates and lower SHD scores.

Additionally, the empirical evaluation compares to variational approaches such as BayesDAG and DiBS that employ augmented spaces. However, for completeness there should also be comparison to other types of Bayesian methods, in particular (i) permutation-based variational methods such as BCDNets; and (ii) MCMC-based approaches such as Gadget [1] which may exhibit qualitatively different behavior (e.g. consistency of MCMC). The related work section should also mention some of these MCMC-based approaches.

For permutation-related variational methods such as BCDNets, we note that we have, in fact, compared to BayesDAG which is representative of the state-of-the-art in this class of methods and has been shown in its original paper to exceed BCDNets's performance. See Section 3.2 of the BayesDAG paper for the equivalence to permutation approaches and Section 6 for the comparisons against BCDNets. For MCMC-based methods, we have now added comparisons against Gadget with $p=20$ and $n=100$. Table R1 shows it substantially underperforms ProDAG. This underperformance is not due to Gadget’s configuration, as we evaluated many settings and report the best-performing one. We will expand the related work section to include a discussion of Gadget and MCMC-based methods.

Table R1 Comparisons against Gadget.

Finally, a seeming limitation is that ProDAG requires the sparsity hyperparameter $\lambda$ to be set in advance, and due to the hard projection constraint this may inevitably lead to errors even in the large-data limit (i.e. if one sets $\lambda<d(d-1)/2$ then one could never learn the complete graph even if it is the true graph).

We treat $\lambda$ purely as a regularisation knob and apply a grid search, so the tuning procedure (see our response below) selects the $\lambda$ best supported by the data. The grid search includes $\lambda=\infty$ so the true graph, however dense, is never ruled out. Only an intentionally small $\lambda$ would exclude it.

Questions

Intuitively, why is it possible to derive the gradients of the optimal solution (Apx C) in closed form but not the optimal solution itself? Given that the iterative procedure only produces an approximation solution to the projection, is there a significant discrepancy between the exact autograd gradients for the approximate solution and the gradients of the optimal solution?

The projection has no closed-form solution because it solves a high-dimensional non-convex problem, but once a stationary point $W^\star$ is attained its KKT conditions are algebraic relations that can be implicitly differentiated, giving the compact gradient formulas in Proposition 2 (in short, local sensitives are easier to obtain than global minimizers). The gradient formulas are evaluated at the $W^\star$ returned by the projection, so they exactly coincide with the gradient autograd would produce at the same point, i.e., there is no discrepancy.

Why did the authors choose a multivariate Gaussian variational posterior, given that the true posteriors can be multimodal? (unless the projection step induces multimodality). It would be interesting to see some examples in low dimensions to visualize the effect of the projection step.

We choose a Gaussian in the unconstrained space $\tilde{W}$ for differentiable reparameterization, closed-form KL, and efficient optimization. We stress, however, that the projection $\tilde{W}\mapsto W$ does indeed induce a non-Gaussian, multimodal posterior over DAGs. We will clarify this point explicitly in the revised text. We are not allowed to include visualizations in the rebuttal, but Figure 1 in the paper roughly illustrates the effect of the projection.

Did the authors try learning a variational posterior with a Laplace-style prior instead of the hard sparsity projection? If so and it performed worse, are there some qualitative insights into why the hard projection step should be preferred?

Yes, during prototyping we evaluated the DAG projection under a Laplace prior and saw noticeably higher SHD scores because Laplace only shrinks edges toward zero and does not set them exactly to zero, so the resulting graphs were too dense. In contrast, the $\ell_1$ ball prior assigns positive probabilities on exact zeros and thus yields crisper sparsity patterns and better edge recovery. See Xu & Duan (2023, arXiv:2006.01340) for a detailed comparison of $\ell_1$ ball versus Laplace.

How did the authors tune the sparsity hyperparameter $\lambda$ on the validation data? How sensitive is the performance of the method to the choice of $\lambda$ (i.e. would it be reasonable to simply choose a single $\lambda$ value and use it across all experiments)?

We sweep $\lambda$ over ten values from \lambda_\min=0 to \lambda_\max (\lambda_\max is the average $\ell_1$-norm of $W$ from the learned posterior when $\lambda=\infty$) and pick the value that maximizes the expected likelihood on a held-out 10% validation set (see Appendix E). Because $\lambda$ directly controls graph sparsity, performance drops if it is set too low (underconnected) or too high (overconnected), so exploring a few points on a grid is preferable to a single value. We apply the same grid‑search to every method with a sparsity parameter.