Continuous Bayesian Model Selection for Multivariate Causal Discovery

Thank you for your positive and encouraging feedback on our work. We appreciate your acknowledgement that the proposed method "allows for learning nonparametric DAGs in a scalable manner" and that our "experiments show competitive performance with the benchmarks". We address your comments in the following.

However, the proof of this theorem only says it follows directly from [Dhir et al., 2024] without any details.

Our construction and proofs earlier in the sections lead us to write the first part of our theorem exactly as Prop. 4.7 in Dhir et al., 2024. Hence, we can directly refer to the proof. We will make this clearer by rewriting the theorem from Dhir et al., (2024) and explicitly stating where we invoke the result.

Many mathemetical details are hidden in the appendix. Especially for Section 3 on the theoretical extension, there is no one theoretical statement in this section while everything is burried in the appendix.

We would love to put more details in the main paper. Our theory requires us to define several concepts not relevant to the rest of the paper, but that allow us to state accurate and precise theorems. As this takes up a lot of space, we thought it best to do this in the Appendix. We do provide a general intuition for our theoretical results in the main paper (L174 LHS). We hope you agree with this approach.

I wonder for two models that are not Bayesian equivalent, why can they have positive probrobability of error

The posteriors may overlap even in the infinite data setting. We can show this with a very simple example. Normalised linear Gaussian models are not identifiable (in the population setting) [1, Appendix D.2]. If a chosen model can approximate a normalised linear Gaussian model, and we sample from the model, there is a non-zero probability that we sample a normalised linear Gaussian dataset. Thus the chosen model must have non-zero probability of error. Hence, while non-linear additive noise models are identifiable in general, additive noise models (also containing linear functions) will have a positive probability of error (but not completely unidentifiable).

[1] Dhir et al., "Bivariate Causal Discovery using Bayesian Model Selection." ICML, 2024.

[2] Hoyer et al. "Nonlinear causal discovery with additive noise models." Advances in neural information processing systems 21 (2008).

Should it be consistent to use the model evidence to find the true DAG? I guess the question is does the sample size in the definition of Bayesian equivalence goes to infinity?

The definition of Bayesian distribution-equivalence holds for any sample size. However, our theorems (B.3 and B.6) hold in the population setting. We do state that in B.3, but we will state that in theorem B.6. Thank you for pointing this out.

Or why do we have identifiability and distinguishabllity and what are their differences?

Identifiability is where the probability of error is zero, whereas distinguishability is where the probability of error is less than random uniform. We will make it clear that these are in the population setting where we define these concepts in L200 LHS.

Thank you for your insightful review. We appreciate your recognition of our contribution to the field of causal discovery. We are glad you think the "theory was well covered and convincing", the paper is "well written and well motivated" and our "extensive experiments significantly outperform the evaluated competitors".

whether the found model respects the Markov condition

Our model itself satisfies the causal Markov condition at the end of training by construction. This happens in the acyclicity thresholding phase of the training (L935). The causal Markov assumption follows trivially from the assumption of no hidden confounders (L82 LHS) and that the noise terms are independent (eq 8).

We are aware of the larger discussion of whether the causal Markov condition implies conditional independences in all generality [1]. However, the same work states that the added assumption of "Modularity" (Property 7.3 in [1]) provides a "fully quantitative solution to the problem of inferring causality from observational data". The ICM assumption we make (L145 LHS) implies a form of modularity [2]. We note that works like NOTEARS do not make this assumption, and the only similarity to our work is that we also use a continuous optimisation scheme.

[1] Dawid et al., "Beware of the DAG!." Causality: objectives and assessment. PMLR, 2010.

[2] Janzing et al., "Causal inference using the algorithmic Markov condition." IEEE Transactions on Information Theory (2010).

It is unclear whether the continuous optimization method finds well-chosen priors

Our continuous optimization does not find the prior, but the posterior given a prior (L255 RHS).

literature on the information-theoretic view of causality (Janzing, Schölkopf, 2010)

The basic principle that allows for distinguishing causal structure is very related to the work you have mentioned [1, Appendix B.2]. We will include a discussion on this.

[1] Dhir et al., "Bivariate Causal Discovery using Bayesian Model Selection." ICML, 2024.

There is no comparison with other score-based approaches

We have run GES and GLOBE for all experiments, results can be viewed here: https://anonymous.4open.science/r/Additional-Results/. GES performs poorly compared to CGP-CDE and the other baselines, except for competitve performance on the 50 var ER1 dataset. This is because GES performs well on sparse graphs, but struggles with denser graphs as it gets stuck in local optima. GLOBE didn't perform well on most of the datasets, except getting the best SHD (but poor SID and F1) for Syntren.

The causal model is not strictly identifiable.

This is true, but we show that for our model the probability of error is small. Our motivation is that restrictions made to gain strict identifiability can be unrealistic. When the data is generated from a different model, identifiability does not hold anyway. We argue that tolerating a small probability of error for gaining more flexibility can be beneficial in these cases. We show exactly this in our experiments.

How do you interpret the poor performance of CGP-CDE in the 3-variable experiment?

We discuss this in Line 368 RHS, although we note the wrong figure was referenced, which we will fix. The good performance of the discrete case shows that our objective empirically identifies the correct causal graph in the multivariate case. The continuous approximation uses (noisy) gradients to find the most likely causal structure, which introduces errors. The comparison to the continuous case thus quantifies the error introduced by the continuous approximation. We believe this shows that the continuous approximation scheme can be improved on. Nevertheless, we emphasise that we clearly outperform baselines in datasets with more variables (Appendix I).

Could another initialization reduce the number of iterations needed?

We were overly conservative with the number of iterations for the warm up and cool down phases. You are correct, better initialisation will reduce the number of iterations required. This can be monitored simply by looking at the training curve.

In the 3-variable experiment... How do you explain this?

The point of this experiment was to show the price we pay for scalability. The difference in performance is due to the continuous optimisation, which introduces errors (see above). These results show that while our Bayesian principle is correct, large improvements can be made in the continuous approximation scheme. As we experiment on all possible 3 variable graphs, this experiment also contains a more dense graphs (relative to the number of variables), which may be harder to find.

Note we also provide other experiments with more variables where the data isn't generated from our model and the ANM experiments where the baseline models' assumptions aren't violated. In both cases, we outperform the baselines.

The reviewer refers to questions 4 and 5 but the questions only go up to 3.

Thank you for your positive and encouraging feedback. We are glad you found the paper "clear and well written", and appreciate your comment that our method is a "significant and promising contribution".

The paper initially suggests that it will solve the problem of restrictive / unrealistic model assumptions encountered when tackling real world data

Previous methods regularly make restrictive functional assumptions that may not hold in practice. In this paper, we address relaxing these restrictive functional modelling assumptions (L12 RHS, L172 LHS). We wholeheartedly agree causal sufficiency and acyclicity can be unrealistic assumptions and are worth relaxing. These are very common in causal discovery algorithms, and the baselines we compare against also make these assumptions. We hope we didn't overclaim on this point and will make clear that this paper is only a step towards more realistic causal discovery (L70 LHS).

Yes it ensures identifiability, but you are essentially finding back what you put in and assume/hope that it matches reality. ...using an in practice hard to justify/assess prior

Our work shows that a relatively simple assumption - independent causal mechanisms (see L148 LHS) - can allow for distinguishing causal structure within a Markov equivalence class (L129 RHS, Theorem B.6). This assumption can be encoded in the prior by ensuring the priors factorise appropriately (L153 LHS).

The specific priors over functional mechanisms are also important. Our approach for this was to ensure, as much as possible, that the prior does not put zero probability on any dataset (see for example [1, Section 1.2]). The model/prior that we chose is not arbitrary, but close to known identifiable models (non-linear additive noise models), as shown by the good performance on additive noise datasets (Appendix I.1). The difference to previous methods is that we do not make hard restrictions - our model allows us to approximate more than just additive noise datasets (L232 LHS).

We note that assumptions made in previous methods that a-priori restrict functional form are also unverifiable. Any unverifiable assumption requires empirical verification. This is exactly what we do in sections 6.2 and 6.3. Here, the data is generated from mechanisms that are different to all the models. We vary other factors as well (graph types and density). Our method outperforms previous methods (full results in Appendix I). We believe this shows the usefulness of the Bayesian approach.

[1] Hjort et al., eds. Bayesian nonparametrics. Vol. 28. Cambridge University Press, 2010.

One of the most striking findings is that the discrete version (DGP-CDE) is vastly superior (FIg.1), but that the authors still choose to focus on the (currently in vogue) continuous approximation, even though it suffers from the exact same problems w.r.t. finding the global optimum in anything with higher dimensions.

This is indeed very interesting. The cost of the discrete version, which requires enumerating over all possible graphs, is too high for more than a few variables. We included this result because the difference in performance between the CGP-CDE and the DGP-CDE clearly shows that, although our principle is correct, there is room for improvement in the continuous relaxation. However, CGP-CDE allows us to scale to a larger numbers of variables.

Cooper & Herskovitz (1992): this old but seminal work already contains a Bayesian score-based method that relies on the prior to select an optimal posterior causal DAG...

The papers you mention are important but only consider simple linear models. However, it has been known that with more complicated models, Bayesian models tend to have an opinion within an MEC [1,2,3]. The main reason (as shown in Appendix B) is because to create equivalent models, the ICM assumption has to hold in multiple factorisations (L129 RHS, Theorem B.6). While it is relatively simple to construct models that do this with linear models [1, Appendix D.1, D.2], it is not clear whether this is possible with more complex models [1, Appendix D.3].

[1] Dhir et al., "Bivariate Causal Discovery using Bayesian Model Selection." ICML, 2024.

[2] Friedman et al., "Gaussian process networks." Proceedings of the Sixteenth conference on Uncertainty in artificial intelligence. 2000.

[3] Stegle et al. "Probabilistic latent variable models for distinguishing between cause and effect." Advances in neural information processing systems 23 (2010).

I would have preferred a modest but more robust aim for e.g. the MEC.

This would obscure the fact that the model (due to the ICM assumption in the prior) has a preference for certain causal structures within an MEC over others (Line 138 RHS and line 720). This preference within an MEC can only be removed by breaking the ICM assumption (Theorem B.6). We thus make use of this preference.

Thank you for your detailed feedback. We appreciate your acknowledgement of the novelty of our approach in applying Bayesian model selection for multivariate structure learning and your recognition of the thorough discussion of our contributions relative to prior work. We are also pleased you found the experimental design sensible and extensive. Additionally, we are glad you highlighted the importance of our method’s ability to distinguish between distribution-equivalent graphs, which is a key strength of our approach. We address your specific comments in the following.

do we need to know the correct prior in advance?

It is not necessary to know the correct prior in advance. A key reason that allows for distinguishing causal graphs is the ICM assumption that is encoded in the prior (L153 LHS, "separable compatible" in Theorem B.6). The prior over functional mechanisms is also important. Our approach for this was to try and choose a model/prior that ensures as much as possible we do not put zero probability on any dataset (see for example [1, Section 1.2]).

[1] Hjort et al., eds. Bayesian nonparametrics. Vol. 28. Cambridge University Press, 2010.

what will happen under misspecification of prior?

We discuss this in L183 RHS. The prior defines what datasets are likely under our model (data distribution). This data distribution can be used to define a probability of error, which is a measure of how distinguishable causal graphs under our model are (eq 7). Given datasets from a separate distribution (data generating process), the accuracy of the estimated probability of error depends on how far the model's data distribution is from the data generating process. Small differences in the model and the data generating process do not result in complete invalidity of the estimated probability of error. As with any method, large variations will mean the estimated probability of error differs from the true probability of error (of the data generating process). For more discussion on this see [1, Section 4.4].

We note that the a-priori functional restrictions imposed by previous methods are also unverifiable assumptions. As with any unverifiable assumption, it is necessary to empirically validate how well it works in practice. We test exactly this in our experiments (Section 6.2, 6.3), where the data generating processes don't match our model prior and our model outperforms the baseline methods.

[1] Dhir et al., "Bivariate Causal Discovery using Bayesian Model Selection." ICML, 2024.

I am not very familiar with bayesian model selection so I am not sure about the significance of extending from bivariate (Dhiretal.,2024) to the multivariate scenario. I will defer to other reviewers regarding this point.

Dhir et al.,2024 consider the bivariate case, but we answer the question: does the theory hold in the multivariate case, and how does the performance scale with number of variables? We theoretically show Bayesian model selection works for multivariate datasets and propose a method to effectively scale to large numbers of variables, overcoming the costs which would be super-exponential with number of variables. We then show that this approach outperforms competing methods in the multivariate case.