Causal Discovery via Bayesian Optimization

In Figure 1, several methods being compared fail to converge. For the tabulated results, have you ensured that all comparison methods have converged?

We would like to clarify that, to demonstrate the sample-efficiency of our approach, it is our deliberate choice to show that the baselines fail to converge before our method. This is because the number of optimization steps, or more generally, the computational budget, is an important hyperparameter in score-based methods that strongly influences the causal discovery performance, but is usually overlooked. Therefore, our empirical study aims to control for this factor, and thus our Figure 1 is intended to show that our method can converge to low SHDs using fewer steps than other baselines, when they have not converged and their SHDs are still high.

To further show that other methods can indeed converge when given more optimization steps, the Table below contains the converged performance of all methods for the case of linear data with 10 nodes (mean±std over 10 simulations), in which 5/6 baselines can converge to a near-zero SHD. COSMO converges to a non-zero SHD, which is in accordance with the reported results in their paper.

Additionally, consider moving the running time details from the appendix to the main content, as the high time complexity is a notable limitation of the proposed method.

Regarding runtime, in the main paper we already have the runtime analysis wrt. performance (last column of Figure 1), showing that our method can arrive at near-zero SHD with less time than other methods, and the runtime details in Appendix H are merely numerical supplementary.

In addition, we believe that when assessing a causal discovery method, running time should be accompanied with a performance metric like SHD, which is the primary measure of causal discovery accuracy, since a method can be very fast and still achieve poor performance. Our results in Figure 1, especially the last two columns, have demonstrated that within a given budget (either number of DAG evaluations or runtime), our method usually achieves better performance than other methods, and given more budget, our method can arrive at very accurate DAGs, while many methods still struggle. For example, in 5 minutes, while our method has achieved a very low SHD=2 for the highly challenging graphs 30ER8, the second-best method is still at SHD=16, and the remaining baselines are still very far behind with SHD>100. This trend can also be observed for nonlinear data, when our method can reach SHD≈0 in 5 minutes, when other methods are still at SHD>6.

Please provide results for nonlinear functions with datasets of 50 and 100 nodes, detailing both performance metrics and running time.

We are conducting these experiments and will update once they are completed.

In the meantime, we look forward to further discussions with you to resolve any outstanding issues.

In Section 4.1, the authors mention that their method incorporates a low-rank adaptation of Vec2DAG. Does this imply an assumption about the data’s structure, as discussed in On Low Rank Directed Acyclic Graphs and Causal Structure Learning?

Yes, our study also assumes the low-rank structure, but implicitly. Specifically, like the sparsity assumption that is usually implicitly imposed on the causal structure in existing methods, we did not explicitly outlined it as an assumption in our presentation. Moreover, our experiments did not enforce the low-rank assumption on the synthetic data, showing that our method is not restricted to low-rank structures. The high empirical performance of our method on general graphs of different types (ER, SF, dense, large, and real structures) suggests that our low-rank representation is robust to many kinds of graphs.

What would occur if ( k < d )?

When $k<d$, the search space is much lower dimensional, allowing us to potentially reach higher-scoring DAGs faster, as shown in Figure 3a. At the same time, the set of DAGs representable by the low-rank representation is reduced from the set of all DAGs, potentially excluding the ground truth DAG. However, our experiments using only rank $k=8$ have shown that our method can still attain a very low SHD even for much larger $d\in\\{30,100\\}$ (see our Figure 1a and 1b). In addition, in Figure 3a, we have empirically analyzed the effect of different ranks and found that even a low rank $k=2$ can still result in a near-zero SHD for complex 30ER8 graphs that have dense connections.

Furthermore, due to the curse of dimensionality, lower ranks, which associate with lower-dim representations, allow for generating more unique candidates compared with higher ranks, so we can more likely meet the higher-scoring DAGs earlier during exploration. We empirically test this by calculating the number of unique DAGs among 1000 random DAGs generated using different ranks in the Table below (the numbers are mean±std over 10 simulations).

This shows that the lower the rank, the more unique DAGs we can consider for exploration. For $k=2$, almost every DAG among 1000 generated DAGs is unique, whereas the lowest number of unique DAGs of less than 10% of all DAGs is obtained using $k=32\approx d$.

Replacing the Gaussian Process in Bayesian Optimization with Dropout is not uncommon, so it may not warrant being highlighted as a novel contribution in this paper.

We would like to clarify that, as outlined in the "Contributions" part of the Introduction, we did not highlight replacing GPs with dropout networks as a novelty in our paper. Rather, we mentioned it as one of the "key design choices", highlighting it as a neccessary leverage to enable BO in causal discovery, which is our novel contribution in the context of causal discovery.

While many prior works employ CAM as a pruning method, I believe this approach may lack justification here. Why would score-based search methods, including this paper, attempt to prune under nonlinear conditions? It's unusual for newly proposed methods to rely on post-processing from an older method.

We would like to clarify that we employed only the pruning step from CAM, not the whole "CAM as pruning method". In CAM pruning, each variable is regressed against its parents in the raw estimated DAG, which may contain many extra edges due to overfitting, by using generalized additive model regression, then insignificant parents are removed. This is itself a reasonable method for pruning extra edges in additive models because it can approximate the data generation process, and thus is widely employed in not just score-based methods (GraN-DAG, RL-BIC, ALIAS, etc.), but also ordering-based methods (e.g., CORL, SCORE, DiffAN, etc.), where pruning is required to remove extra edges in the fully-connected DAGs induced by the returned causal orderings.

Please compare this baseline method, Truncated Matrix Power Iteration for Differentiable DAG Learning, to your approach.

Following your suggestion, we have added this baseline to the main experiments in our revision. We used NOTEARS's implementation combined with TMPI DAG constraint and call the method 'NOTEARS+TMPI' in our experiments. This method's performance is relatively strong to some other baselines, especially in sparse graphs (Figure 1b and Table 1 of the revised manuscript) and nonlinear data (Figure 1c and 2 of the revised manuscript), however, it is still surpassed by our method in all cases.

We are greatful for Reviewer ZehJ's invaluable comments. We try our best to address your concerns and we look forward to your responses.

To my knowledge, the CBO series of papers assumes that the DAG structure is known and primarily focuses on optimizing policies with this prior knowledge. Therefore, these papers may not be directly relevant to the active causal discovery literature. Please revise this in the introduction to reflect the distinction.

Thank you for your kind suggestion. We would like to clarify that we mentioned them to highlight that our application of BO is entirely different from existing causal discovery studies, since it is easily mistaken that our method is another CBO method, and we have revised our manuscript to be more precise.

The synthetic dataset, as first used in NOTEARS, is known to be relatively easy to learn, making pursuit of very low SHD scores less meaningful in recent research, especially the very simple linear gaussian case.

This is not neccessarily the case. We would like to clarify that our dataset involves much denser graphs, and achieving a low SHD in this case is very challenging for existing methods. Specifically, for dense graphs 30ER8, most baselines, especially the sortnregress approach introduced in Beware of the Simulated DAG, obtains a very high SHD of 100+, while our method's SHD is only 1.6±1.5 (see our Figure 1a and Table 5).

Indeed, this can be confirmed by calculating the Varsortability from Beware of Simulated DAG. If Varsortability=1 then the causal structure can be recovered simply by sorting the nodes by increasing variances, i.e., using sortnregress. Below, we show the Varsortability (mean±std over 100 simulations) for 30-node graphs with varying densities. The results suggest that while sortnregress may be useful for very sparse graphs, it becomes invalid very quickly for denser graph.

In addition, since Varsortability=0 for 30ER8 graphs, it may be tempting to think that sorting the nodes by decreasing variance can help. To test this, we again apply sortnregress with both sorting directions in the Table below (mean±std over 100 simulations).

This result indicates that our synthetic data is not easy and our method's ability to reach low SHDs in such cases is significant.

How does your method perform on this dataset after standardization, as described in the paper Beware of the Simulated DAG?

Following your suggestion, we are conducting experiments with standardized data and will update you once they are completed.

The proposed method appears to be limited to ANMs in causal discovery, which restricts the scope of the paper. It may be more accurate to frame the task as DAG structure learning or Bayesian structure learning rather than causal discovery.

We would like to clarify that our method is not limited to ANMs. We have demonstrated in Appendix F.1.4 that our method can also handle logistic nonlinearity. In general, our method can be applied to other causal models, as long as a suitable scoring function is provided.

Dear Reviewer ZehJ,

Thank you for bringing this to our attention!

We assure you that our comparison remains fair.

We would like to clarify that, we did not limit DAGMA to just 50k steps in the last column of Figure 1, but in fact, we limited it to millions of evaluations, and we removed its early stopping condition to allow it to run for an arbitrary amount of time, so its cut-off runtime in the last column in Figure 1 can be equal to all other methods. Without this, we cannot answer the question "what is the performance of each method if it is run for x steps / y minutes?". We have also mentioned this in Appendix E: "Regarding the number of evaluations, for all methods, we run more than needed then cut off at common thresholds."

Specifically, for DAGMA, we run it for around 24 million steps by setting T = 800 (the default is T = 5), so that we can either use a threshold at 50k evaluations or at 5 mins with the full learning progression of the method. In deed, DAGMA can reach 50k steps in just a few seconds as usual, but costs around 5 million steps at 5 mins for 30ER8 data.

We strongly believe that our comparison is reasonable for controling for both the number of evaluations and runtime, and we will most certainly clarify this better in the revision.

In light of this clarification, we hope you can reconsider your assessment, and we are open to further discussions regarding any outstanding concern. Thank you again for your time!

Dear Reviewer ZehJ,

In our last response, we have clarified that DAGMA runs for minutes in our study instead of just a few seconds is simply because we allowed it to, so that we can answer the question of How good a method can become within a given runtime?.

In addition to that, we have uploaded a revision clarifying your concerns.

Particularly, in Appendix D.4, we have included:

Additionally, in this study, we evaluate the performance of various methods with respect to both the number of steps and runtime, addressing two independent questions: “How accurate can a method become given a fixed number of steps?” and “How accurate can it be within a given runtime?”. To ensure fair comparisons, our second question accounts for potential biases in measuring performance solely by the number of DAG evaluations. This is particularly important for methods like gradient-based approaches (e.g., DAGMA), which may require many steps but still exhibit low overall runtime.

To address this, we use runtime as a more equitable efficiency metric. Specifically, we set a high number of steps for all methods (e.g., we use T=800 iterations instead the default of only T=5 for DAGMA on linear data) and disable early stopping if applicable, to capture their progression over an extended period of time. We then truncate the tracking data, which contains performance metrics and timestamp at every step, either at a fixed number of steps or a specified runtime, as illustrated in Figure 1. This ensures that the results in the last column of Figure 1 are not constrained by the number of steps. For instance, at the 5-minute mark in the last column of Figure 1a, DAGMA completes approximately 5 million steps compared to only 50,000 steps in the third column.

We have taken every possible measure to ensure fairness in our evaluations and strongly discourage any form of misconduct. In light of our clarification, we would greatly appreciate it if you could reconsider your evaluation.

Thank you sincerely for your time and effort in assessing our paper!

We thank Reviewer VYj5 for the valuable insights. We address your concerns in the rebuttal below.

Some of the baselines considered are not adequate

Following your suggestion, we have enhanced the baselines considerably. Our revised manuscript has incorporated several additional baselines, including GOLEM and NOTEARS with TMPI constraint (see Figure 1, 2, and Table 1). These gain better performance than other baselines in several cases, e.g., they both achieve slightly better performance than DAGMA for large graphs (Figure 1b), and much lower SHDs compared with other baselines for real-world structures (Table 1). In addition, NOTEARS+TMPI also obtains the third-best performance with SHD≈7.2 for nonlinear data with GPs (Figure 1c), and we have also improved DAGMA's performance significantly in this setting, from SHD≈28 to SHD≈17 (Figure 1c).

Some of the results may seem too good to be true. For example, achieving a SHD of 1.6 with only 1000 samples across 30 nodes and 240 edges seems highly challenging due to finite sample error. This concern is especially relevant when dealing with nonlinear data. (I look forward to the authors' clarification/explanation on this, and please correct me if I misunderstood anything.)

Thank you for recognizing the hardness of our data setting. We firmly assure you that this accuracy is possible with our method, thanks to BO's ability to effectively optimize the DAG score. Additionally, this level of accuracy has also been achieved in ALIAS, even though with far many more DAG evaluations than ours. To further prove this, we have attached our source code in the Supplementary, so that you can reproduce our results.

Is there a reason why the paper considers only identifiable models? That is, why general linear Gaussian model cannot be learned by the BIC-NV score?

We would like to clarify that we considered both identifiable and unidentifiable models in our experiments. While we examined identifiable models in the main paper, we have also studied two unidentifiable causal models in the Appendix, namely the general linear models (without equal noise variance) in Appendix F.1.3 and logistic models in Appendix F.1.4, where our method still works well and can find the highest scores with lowest structural errors most of the time.

Although BIC-NV is given, it seems that the experiments focus on equal variances. I would suggest adding experiments for different variances as well.

We would like to clarify that our experiments examined both equal and non-equal variances settings. Specifically, our non-linear experiments (Figure 1c and 2) are for non-equal variances, in which BIC-NV is used for scoring and optimization. We have mentioned in Section 5.1.2 that the noise variances for GP data are sampled uniformly in $[0.4, 0.8]$, and for the real Sachs dataset, the noises are unlikely to have equal variances. Our revised Section 5.1.2 has further clarified this.

Baselines: For linear case, GOLEM may also be included. Also, adding the results for more conventional search methods, such as GES/FGES, may also be helpful.

Thank you for your kind recommendation. We have incorporated GOLEM results for linear data in Figure 1 and Table 1 of the revised manuscript. Regarding conventional search methods like GES/FGES, since they are usually shown to perform poorly in many recent studies, and to avoid cluttering the presentation with too many baselines, we did not include them. However, following your suggestion, we will add them in the next revision.

Why does DAGMA-MLP performs so poorly for nonlinear data? The TPR is close to 0. If the reason is due to instability in optimization, the paper may consider adding NOTEARS-MLP that may be more stable.

Thank you for your kind suggestion. We have run NOTEARS-MLP for nonlinear data and obtained an SHD of 11.8±2.59. We found that NOTEARS-MLP set the $\ell_1$ and $\ell_2$ weights to zero, while DAGMA-MLP used $\ell_1=0.02$ and $\ell_2=0.005$, leading it to predict very sparse graphs and thus resulting in low TPR. Therefore, we have updated our Figure 1c with $\ell_1=\ell_2=0$ for DAGMA, in which it achieves a significantly higher performance with SHD≈17. Nevertheless, we have also included NOTEARS-MLP in conjunction with TMPI constraint to achieve an SHD≈7 (see our Figure 1c of the revised manuscript).

For Section 5.2, specifically Sachs data, did the paper use linear or nonlinear version of the method?

As mentioned in our Appendix E, we used the nonlinear version of our method with GP regression and BIC-NV scoring for the Sachs data, as with all considered baselines.

We hope this rebuttal sufficiently addresses your concerns, and we look forward to further discussions to resolve any outstanding issues.

We thank Reviewer UyS6 for your positive evaluation. Your concerns are addressed in the rebuttal below.

While there exist some causal discovery algorithms with Bayesian optimization, it seems not proper to state “To our knowledge, this is the first score-based causal discovery method based on BO ”. I think it should be corrected.

Thank you for your kind comment. We have revised the manuscript to highlight that our study is the first score-based causal discovery method based on BO for purely observational data.

Throughout the paper, from the experiments, it is demonstrated that the proposed method can give better performances in both accuracy, sample-efficiency, and scalability, compared with other SOTA baselines. Generally, such a great method needs more assumptions or conditions to be satisfied. But intuitively, I cannot find these assumptions or conditions. Did this method have some other implied assumptions or conditions (like more hyperparameters)?

We would like to clarify that our method does not require any additional assumptions compared with existing ones. Specifically, our assumptions outlined in Section 3.2 include causal sufficiency, causal minimality, and identifiable models. These assumptions are standard and are similar to many studies, e.g., NOTEARS, DAGMA, RL-BIC, CORL, etc. Our scoring functions and evaluation data are also the same as RL-BIC, CORL, ALIAS, etc.

Our method is better than the baselines thanks to the effectiveness of BO, in which we predict the scores of DAG candidates to prioritize the DAGs that are most likely to have higher scores, before actually exploring/evaluating them. This helps us avoid examining non-informative DAGs, e.g., DAGs that we are certain to have low scores, and at the same time reveals higher-scoring DAGs earlier. Meanwhile, most existing methods do not effectively take into account past exploration data to make an informed candidate selection, thus resulting in more unneccessary trials than our method.

In addition, our code is now available in the Supplementary material, which can be used to confirm our strong performance.

In experiments, it would be good to compare some Bayesian causal discovery methods (Deleu et al., 2022: Tranet al., 2023; Annadani et al., 2023), since they are all causal discovery methods. Or explain the reasons why not comparing with them.

We would like to explain that we did not compare with them for several reasons. Firstly, our study focuses on score-based causal discovery with an emphasis on sequential optimization, which revolves around solving the optimization problem in Eq. (1), so Bayesian causal discovery methods do not directly fall into the same setting as ours. Secondly, it is not very common in the literature for a point-estimate causal discovery study to compare with Bayesian causal discovery studies, so we simply followed common practice. However, following your suggestion, we will certainly add them in the next revision.

The code is not available for reproduction.

We have made our code available for reproduction in the Supplementary material.

In Eq.(4), is the $R$ matrix strictly upper-triangular?

No, $R$ does not need to be strictly upper-triangular.

In Figure 1(b), did the authors still use 1000 samples for the large-graph experiments? $n=1000$ for the graph with 100 nodes?

Yes, we used 1000 samples for all main experiments. We hope this clarification helps you find the performance of our method significant. You can confirm our results with the code provided.

In Figure 3(a), why in general smaller $k$ could obtain higher performances, compared with the full-rank cases?

In short, this is because it is much more challenging to search in a very high-dimensional space compared to a lower-dim one. Specifically, for full-rank cases, the search space is much larger and sparser than the low-rank ones. Due to the curse of dimensionality, sampling the same number of random DAG candidates in the full-rank search space tends to lead to fewer unique candidates compared with a low-dim one, reducing the chance to meet the optimal solution earlier.

To empirically verify this, we calculate the number of unique DAGs among 1000 random 30-node DAGs generated with different ranks in the Table below (the numbers are mean±std over 10 simulations).

It can be seen that, typically, the lower the rank, the more unique DAGs we can pre-examine for exploration. For k=2, almost every DAG among 1000 generated DAGs is unique, whereas the full-rank representation is higher-dim and can only generate fewer than half the unique DAGs.

We hope this rebuttal addresses your concerns, and are open for further discussions to resolve your remaining concerns.

We wholeheartedly thank Reviewer gMGJ for the positive assessment. Below we try our best to address your question.

Could the authors give more details on how to ensure the binary adjacent matrix is a DAG in the optimization steps?

In short, the binary adjacency matrix is always ensured to be a DAG thanks to our DAG representation explained in Section 4.1.

In more details, let us first recall our DAG representation:

$\tau(\bf{p},\bf{R})=H(\mathrm{grad}(\bf{p}))\odot H(\bf{R}\cdot \bf{R}^\top)$

The acyclicity of this binary matrix is enforced through the first term $H(\mathrm{grad}(\bf{p}))$. This term is an adjacency matrix of a directed graph where there is an edge $i\rightarrow j$ if and only if $p_i < p_j$. By contradiction, assuming there is a directed cycle $i_1\rightarrow i_2\rightarrow\ldots\rightarrow i_1$ in this graph, then it must be that $p_1 < p_1$, which is always impossible. Therefore, there cannot be any directed cycle in this graph, rendering it a DAG. The second term in the equation above, $H(\bf{R}\cdot \bf{R}^\top)$, plays the role of a selection matrix that chooses some edges in the DAG induced by the first term to include in the final result, so there is no new directed cycle and thus the obtained binary matrix is a DAG.

Dear all reviewers,

In our recent revision, we have included a new section (Appendix F.1.5) to delve deeper into how our method performs on standardized data, which may makes causal discovery less trivial, as discussed in the study Beware of the Simulated DAG!. We have found that our method continues to perform really well on both linear and nonlinear standardized datasets, further solidifying its reliability:

• For linear-Gaussian data, standardization renders the causal model unidentifiable, thus making it impossible to consistently recover the correct causal structure. Yet, the empirical results reveal that our method still significantly outperforms other baselines with very low errors (SHD≈3 for 10ER2 graphs, while the second-best SHD is ≈20). Specifically, our method barely produces any missing or extra edge and only mis-orients a few edges due to the unidentifiability of the causal model.
• For nonlinear ANM data, since the causal model remains identifiable after standardization, the causal discovery result is similar to our experiment on the same datasets without standardization. More particularly, our method can still accurately recover the causal structures with an SHD≈0.

For the detailed analyses, please refer to Appendix F.1.5 of our latest revision.