Continuous Bayesian Model Selection for Multivariate Causal Discovery
Authors: Anish Dhir
Ruby Sedgwick
Avinash Kori
Ben Glocker
Mark van der Wilk
Venue: ICML 2025 poster (2025)
Keywords: bayesian model selection
gaussian process
causal discovery
Primary area: general_machine_learning->causality
TL;DR: Bayesian model selection for multivariate causal discovery allows for using more flexible models which leads to better performance.

Abstract: Current causal discovery approaches require restrictive model assumptions in the absence of interventional data to ensure structure identifiability. These assumptions often do not hold in real-world applications leading to a loss of guarantees and poor performance in practice. Recent work has shown that, in the bivariate case, Bayesian model selection can greatly improve performance by exchanging restrictive modelling for more flexible assumptions, at the cost of a small probability of making an error. Our work shows that this approach is useful in the important multivariate case as well. We propose a scalable algorithm leveraging a continuous relaxation of the discrete model selection problem. Specifically, we employ the Causal Gaussian Process Conditional Density Estimator (CGP-CDE) as a Bayesian non-parametric model, using its hyperparameters to construct an adjacency matrix. This matrix is then optimised using the marginal likelihood and an acyclicity regulariser, giving the maximum a posteriori causal graph. We demonstrate the competitiveness of our approach, showing it is advantageous to perform multivariate causal discovery without infeasible assumptions using Bayesian model selection.

Continuous Bayesian Model Selection for Multivariate Causal Discovery
Anish Dhir
Ruby Sedgwick * 12 Avinash Kori!
Ben Glocker'
Mark van der Wilk 3
Abstract
Biology (Sachs et al., 2005) to Medicine (Feuerriegel et al., 2024) to Economics (Hicks et al., 1980). There are two main approaches for learning the ground truth causal structure from data. The first approach assumes a restricted model class for the data generating process (Peters et al., 2017, Ch. 4), but its guarantees often fail when these assumptions do not hold. The second approach requires interventional data for all variables, which can be costly, ethically chal- lenging, even impossible (Li et al., 2019). Consequently, learning unique causal structure often relies on impractical assumptions, limiting its real-world applicability.
Current causal discovery approaches require re- strictive model assumptions in the absence of in- terventional data to ensure structure identifiability. These assumptions often do not hold in real-world applications leading to a loss of guarantees and poor performance in practice. Recent work has shown that, the bivariate case, Bayesian model selection can greatly improve performance by ex- changing restrictive modelling for more flexible assumptions, at the cost of a small probability of making an error. Our work shows that this approach is useful in the important multivariate case as well. We propose a scalable algorithm leveraging a continuous relaxation of the discrete model selection problem. Specifically, we employ the Causal Gaussian Process Conditional Den- sity Estimator (CGP-CDE) a Bayesian non- parametric model, using its hyperparameters to construct an adjacency matrix. This matrix then optimised using the marginal likelihood and an acyclicity regulariser, giving the maximum a posteriori causal graph. We demonstrate the com- petitiveness of our approach, showing it is advan- tageous to perform multivariate causal discovery without infeasible assumptions using Bayesian model selection.
Recent work has shown that, the absence of interventional data and the bivariate case, Bayesian model selection can more accurately infer the unique causal structure of system (Dhir et al., 2024). Unlike previous methods, including other Bayesian methods which impose constraints on the func- tion class, this framework allows the use of more flexible model classes. This comes at a cost of loosened identifia- bility guarantees but with the ability to posit more realistic assumptions. However, Dhir et al. (2024) only study the bi- variate case, while most applications require understanding causal relationships between multiple variables. It is thus an open question whether the advantages of this method hold in the multivariate case. Furthermore, the approach in Dhir et al. (2024) requires computing and comparing the posteri- ors of all possible causal graphs. As the number of graphs grows super-exponentially with the number of variables, this becomes computationally intractable for scalability.
A potential solution lies in continuous optimisation causal discovery approaches (Zheng et al., 2018: Nazaret et al., 2024). These approaches interpret causal discovery as a single continuous optimisation problem that is amenable to gradient based methods. The emphasis of this line of work has been either on scaling to a larger number of vari- ables loosening model assumptions. Our contribution tackles the latter. We argue that the potential advantages of Bayesian methods enable them to outperform their non- Bayesian counterparts.
1. Introduction
In many systems, such as protein signalling networks, vari- ables causally relate to each other, since changing variable only modifies certain variables (Sachs et al., 2005). Uncover- ing these unique underlying causal structures from data can allow us to gain new insights in a wide range of fields, from
'Equal contribution 'Department of Computing, Imperial Col- lege London, London, UK *Xyme, Oxford, UK Deparrment of Computer Science, University of Oxford, Oxford, UK. Correspon- dence to: Anish Dhir <anish.dhir 3 imp.rial.ac.uk>..
We investigate whether and how well Bayesian model se-
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).
'These methods are usually concerned with inferring a distri- bution instead of single causal structure.

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
lection can recover the ground truth causal structure in mul- tivariable systems. We use model based on the Gaussian Process conditional density estimator (GP-CDE) (Titsias & Lawrence, 2010; Lalchand et al., 2022) that is more flexible than traditional identifiable models that rely on restrictive model assumptions. Interpreting the model hyperparameters as a graph, we optimise over them using a Bayesian model selection-based loss and regulariser to enforce acyelic- ity, ultimately producing the maximum a posteriori (MAP) causal graph. Our experiments address two critical ques- tions: (1) How much performance is lost compared to enu- merating every causal graph? (2) How does our scalable Bayesian model selection compare to existing multivariate causal discovery methods in terms of performance? Our results show that Bayesian model selection outperforms methods enforeing strict identifiability even at larger scales. While it introduces small error probability, the flexibility afforded enables more realistic and practical multivariable causal discovery, advancing toward real-world applications.
2.2. Learning Causal Structure
Formally, the task is to recover a DAG & given dataset of N samples with D variables, X E RNXD, Without any in- terventional data, only the Markov equivalence class (MEC) of a DAG can be recovered with the faithfulness assump- tion (Pearl, 2009) and with infinite data (Appendix A). To recover a unique DAG, and differentiate within an MEC, we require additional assumptions.
2.3. Learning Causal Structure with Functional Restrictions
One way to recover a DAG is to impose a priori restrictions on the allowable conditional distributions CijPA These restrictions are chosen so that the joint distributions implied by causal model cannot be expressed by any other causal model specifically, for P E Fg it holds that P ş FH for any H + 9 (Guyon et al., 2019, Ch. 2). Thus, if data is sampled from one of the causal models, the causal structure can be identified. As an example, additive noise models (ANM) (Hoyer et al., 2008) restrict all CiJPAg() to the form I (())a)()) + e for some arbitrary non-linear function f. and arbitrarily distributed e. If we consider different graph H but the same ANM restrictions on C())A it is not possible to approximate the joint induced by the original set of conditional distributions Cg- = Ille CiIPAo()) using Ch = HEy CjjPA This allows the graph G to be identi- fied. For lot of these restrictions, the maximum likelihood score can identify the causal structure (Zhang et al., 2015; Peters & Bühlmann, 2014; Immer al., 2023).
2. Background
Here, we outline the assumptions relevant the task of causal discovery. Additional background is in Appendix A.
2.1. Causal Model
We assume that data X E RDxDD can be explained by a Causal Model, which we refer to Mg. We make the common assumption that there are no hidden confounders.
Methods relying on restricted model classes also restrict the datasets they can model. Specifically, assume that data was generated by a distribution II, then it may be that II & Fg for all G. Here, none of the causal models can approximate the true data- generateeernnggg distribution and guarantees about causal identifiability no longer hold. That is, the models are misspecified. A solution would be loosen the restrictions on Cy)PAp()) providing greater flexibility in approximating distributions. However, this could also allow causal models with different DAGs to express the true data-generating distribution, thereby losing the ability to recover the true DAG (Zhang et al., 2015, Lemma 1).
Definition 2.1. (Consistent with Pearl (2009) but using no- tation from Dhir et al. (2024)) A Causal Model Mg == (G, Cg) is defined as a tuple containing a DAG 9 with vertex set V and edge set E, along with a set of conditional distribu- tions for each variable (i) where (1) denotes the parent index set of node index in G. The set of all possible conditionals is the Cartesian product Cg := Ievv C))PAg()) The set of all possible conditional distributions induces set of joint distributions denoted Fg.
An element of Co induces joint distribution over D ran- dom variables X = (X1...., XD} when its constituent members are multiplied together. That is, for an element (P. : E V)= E Cg the joint is II:2 (. Hence Cg induces set of joint distributions Fg (II2_ Px(Xx(XPAA(t)) (P; i E V) E Ce} (Dhir et al., 2024). We also assume likelihood modularity (Geiger & Heckerman, 2002), that is for variable X, if (i) = PAg (1) for 9 # G', then CipAA(()) = C, (A) This means that for graphs where variables have the same par- ents, we assume the same set of conditional distributions.
2.4. Learning Causal Structure with Bayesian Model Selection
Restricted function classes are required for identifiability within an MEC as it implies that, with infinite data, only one of the causal models can approximate the data distribution and thus achieve higher score such as the maximum likeli- hood. Without restrictions, multiple causal structures can yield the same maximum likelihood score, making it impos- sible to distinguish between them. However, using such an
2

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
unrestricted causal model with Bayesian model selection- based score can restore the ability differentiate between causal structures. This is because, while multiple structures may share the same maximum likelihood score, they typi- cally yield distinct scores under Bayesian model selection. This insight was demonstrated for the two-variable case in Dhir et al. (2024), building on prior efforts (Friedman & Nachman, 2000; Stegle et al., 2010). While this approach no longer guarantees strict identifiaiility and admits small probability of error, the error probability can be empirically estimated and remains low with appropriate models (Dhir et al., 2024). The key advantage of this framework lies in its flexibility. By enabling the use of broader model classes, mitigates the risk of misspecification, which is a com- mon limitation where the assumptions and guarantees of restricted models fail (Section 2.3). Our work builds on these insights.
Multiple Bayesian causal models may achieve the same likelihood score given infinite data, but the prior and the encoded ICM assumption play a key role in yielding distinct marginal likelihood scores. The ICM assumption, which encodes the causal factorisation in the prior, typically does not hold in any other factorisation (Dhir et al., 2024). To see this, consider the model with graph X - Y, with Cx parametrised by 0 and Cr|x parametrised by é. Here, placing prior on Cx-x - is equivalent to placing a prior on the parameters 0 and Φ. ICM implies that the priors over 0, Φ factorise, which means the joint factorises as
x (0))(()P(X1)))P() o).
(3)
We can find the reverse factorisation of the variables for the same model using Bayes' rule. Here, in general, it is not possible to parametrise Cy and Cxly with independent parameters (Dhir et al., 2024)
In Bayesian model selection, the evidence for causal model Mg (as defined in Section 2.1) is quantified by the posterior
Equation (3) = (0) (X))) (4)
=
(5)
P(Mg|X) & P(X
(1)
where a, are some reparametrisations. Hence, if T(a, 3) * (T)))(()) the model with the graph X + Y only satisfies the ICM assumption, as specified by factorised priors on the set of conditionals, in the causal factorisation of the joint (P P. E Cx. E Cy1x).
where P(Mg) is the prior over the causal model. As we assume no preference over causal direction between vari- ables, our prior belief is equal over graphs in the same MEC. The posterior within MEC is then dependent only on the term P(X]Mg), known as the marginal likelihood. Given that our task is to infer the single ground truth causal struc- ture, we select the causal model Mÿ based on the highest posterior probability (Kass & Raftery, 1995)
Given a prior, the marginal likelihood is
D P(XIMg) = /
(6)
My = argiaxx P(Mg)X).
Mg
The marginal likelihood score has an automatic Occam's 's razor effect that balances model fit along with complexity (MacKay, 2003, Ch. 28). This means that although mul- tiple causal models may have parameter settings that fit a given distribution, some models may provide a simpler ex- planation than others. Specifically, those whose prior, and encoded ICM assumption, align with properties of the data generating process. This allows for distinguishing between causal models even with loosening the restrictions discussed in Section 2.3.
To calculate the marginal likelihood, a prior (denoted т)) over the distributions in the causal model must be specified. This defines Bayesian Causal Model.
Definition 2.2. (Dhir et al., 2024) A Bayesian Causal Model, denoted as (Mg. T), is defined as a causal model Mgg (Definition 2. with a prior distribution ı over Cg.
Each Bayesian causal model should encode the common and foundational principle that change in one causal mech- anism should not affect any others (Aldrich, 1989; Pearl, 2009). This can be formalised by the independent causal mechanism (ICM) assumption (Janzing & Schölkopf, 2010). The ICM assumption implies that, in the causal factori- sation, knowledge about any distribution P. E C)) should not inform any other P3- E C)IPA9()) for A i. In the Bayesian approach this implies defining factorised pri- ors on each set of variable distributions where ,. E ((C))Ae()) is the prior on the distributions for variable X,, factorised according to the causal structure g, and P is the set of all distributions over an object. This is the only constraint we place on the prior.
Dhir et al. (2024) only consider Bayesian model selection for the two variable case. It is not clear how this approach will perform when the number of variables is increased, especially in contrast with competing methods. We study the multiple variable case where we address the following challenges: First, we use the theory in Dhir et al. (2024) to clearly show that Bayesian model selection can be applied to causal discovery with multiple variables. Secondly, the number of variables increases, the number of possible DAGs
π = Il Ti
Exceptions exist, notably normalised linear Gaussian models which are known be unidentifiable.
3

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
grows super-exponentially, making direct comparison of marginal likelihoods across all DAGs infeasible. To address
The advantage of the Bayesian approach is apparent when considering cases where the underlying causal model is unidentifiable with maximum likelihood. This is the case for graphs within an MEC, where functional restrictions have not been made.
this, we propose a scalable approach for applying Bayesian model selection to multiple-variable causal discovery.
3. Distinguishability with Multiple Variables
Remark 3.3. Suppose that the ICM assumption only holds in the causal factorisation of the Bayesian causal model. Bayesian model selection can then distinguish within an MEC (Theorem B.7).
We use theory from Dhir al. (2024) to show that in the multiple variable case, Bayesian causal models find the correct DAG using only priors with the ICM condition, rather than hard model restrictions. We proceed by first showing that Bayesian causal models can be identified up to an MEC. Then, we show that Bayesian causal models can differentiate within an MEC, without functional restrictions, but with some probability of error. We only make remarks on key results in this section, but the full set of assumptions, theorems, and proofs is in Appendix B. Note that these definitions hold in the population setting.
Thus, without restrictions and if the ICM condition only holds in the causal factorisation, Bayesian model selection can be used to distinguish Bayesian causal models, in that the probability of error is lower than a uniformly random chosen graph. The model class we choose only satisfies the ICM condition in the causal factorisation (Dhir et al., 2024, Appendix D.3). The probability of error for our chosen model class was shown to be very low for the two variable case in Dhir et al. (2024, Section 4.3). We estimate the probability of error for the multivariate case and show that it is low in Section 6.1.
First, we define the probability of making an error, given datasets from the chosen Bayesian causal model.
Definition 3.1. (Dhir et al., 2024) Given a score S, such as the marginal likelihood, we define the probability of error for model Mg as
As is common with identifiability results, the probability of error estimate assumes that the data is generated from one of the candidate causal models. However, this assumption may not hold in practice. Dhir et al. (2024, Section 4.4) demon- strate that the model's 's estimated probability of error and the true probability of error (from the actual data-generating process) are bounded by the total variation distance between the model and the true distribution. Consequently, mild violations of model assumptions, such as the prior, can still yield accurate error probability estimates.
P(E\Mg) P(XIMg)dX, JR
where R {X: S(XIM)) > Mg) for H. * G). and S(X ((X)M)) is the score achieved by model MM on dataset X.
That is, given datasets X from a Bayesian causal model with DAG G, the integral over datasets where the wrong causal model chosen by the score. For Bayesian causal models, this quantifies the overlap in the posteriors. Note that, simi- lar to previous work, this definition assumes that the data is sampled from one of the models under consideration. The probability of error can thus be estimated empirically for Bayesian causal models by sampling a graph, then sampling datasets from a chosen Bayesian causal model, and infer- ring the causal graph using the chosen score (Dhir et al., 2024). An identifiable model implies probability of error of zero in the infinite data limit (Dhir et al., 2024; Guyon et al., 2019). In contrast, an unidentifiable model will have a probability of error equal to that of uniformly random chosen graph. Similarly, we say a model is distinguishable if the probability of error is less than that of uniformly chosen random graph.
4. CGP-CDE: Causal Gaussian Process Conditional Density Estimators
Although the posterior allow for more accurate causal discovery with fewer restrictions, computing and comparing the posterior of every causal structure becomes intractable when the number of variables is increased. We propose a method that scales up Bayesian model selection based causal discovery by using the insights of Zheng et al. (2018) (Appendix A.2). We first describe the model that we use. then we show how the model can continuously parameterise the space of graphs. This is followed by loss function that maximises the posterior probability and ensures that the final graph is a DAG.
Remark 3.2. With the faithfulness assumption, Bayesian model selection can identify a Bayesian causal model up to the MEC of its graph (Theorem B.3).
4.1. The GP-CDE Model
The Bayesian model selection approach to causal structure learning described in Section 2.4 requires defining distribu- tions and priors for a causal model. To take advantage of the Bayesian framework, we wish to use a model for each Ci)PAg(0) that is more flexible than previous attempts. For
Although this is well known for certain parametric models (Heckerman, 1995; Chickering, 2002), we show that it is true for the general non-marametric case.
4

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
this, we use version of the Gaussian process conditional density estimator (GP-CDE) model (Dutordoir et al., 2018). Here, (1) is a class of nonparametric functions and the prior over the conditionals is a Gaussian process prior.
variable Xy. This implies an adjacency matrix
{ 0
Osj i
(10)
otherwise.
The input to the estimator for each variable is determined by the causal graph of the causal model. For a given causal DAG g. the joint p(X,F f. W/A, Mg) of our model is
A value of 0.g 0 implies that changing X, does not change X, and results in an absence of the edge X, + X; in A.
D N II
The Gaussian process prior (f, X ([f][PAAall (1), Fi 0i, w.) is then parameterised as follows
plane|f (PPp(())))))) As, wns)p(wne), (8)
where W E RND are the latent variables, & E R are the likelihood variances, A; are the kernel hyperparame- ters, and Mg denotes the causal model with DAG g. Each variable x; is modelled as the output of a function f, with the parents of the variable in G. (t), and a latent vari- able w; being the inputs. The inclusion of w; ~ I) allows for each variable to be described with non-Gaussian and heteroscedastic noise, greatly increasing the expres- sivity of the model (Dutordoir et al., 2018)². The term p(f, XpAa()):A; A,, (ins) is the Gaussian process prior (Ras- mussen, 2003), and equal
(O, .o.., ((X w.). (X w::'));
(11)
where PA (i) denotes the parents of X; in G (the graph implied by the adjacency A), and X a denotes all variables except X₁. Thus, for a variable X,, all the other variables are inputs, and the values of the hyperparameters 0, control the dependence of the rest of the variables X.,, on X,, and hence define the parents of X,, The hyperparameters the above model thus continuously parameterise the space of graphs. The exact kernel we use is in Appendix C.I.
4.3. Priors on Graphs
(0, KA, ((XPAo()) ws), (X), w.))),,
As discussed in Eggeling al. (2019), the lack of a prior on graphs leads to higher weight on denser graphs due the larger number of dense graphs versus sparse graphs. Hence, to even out this effect, we place a prior on graphs. By intro- ducing priors on the hyperparameters 0 in our model, we can effectively encode priors on the graph in Equation (10). We thus place Gamma prior with parameters chosen to pre- fer small values of 0, P(O)- = Gamma(n, 3). Note that our prior is symmetric between graphs within the same MEC. That is, we do not impose preference over the causal di- rection between variables, only over the number of edges in the graph. Other priors on graphs could be considered (Eggeling et al., 2019), but we leave that for future work.
where k(,,) -) is a chosen kernel matrix parameterised by hyperparameters A.
The above approach requires defining a separate model for each causal DAG. Next, we show how to interpret the hyper- parameters of this model as an adjacency matrix. This will allow for continuously parameterising the space of graphs.
4.2. Continuous Relaxation: The CGP-CDE Model
The key insight to continuously parameterise graphs (not DAGs) in the above model is that the hyperparame- ters control dependence between variables (Williams & Rasmussen, 1995). In our Gaussian process prior for a variable X,, we use kernels with different hyperpa- rameters Ai for each dimension j- In general, ker- nels of this form factorise as KA, ((X, w:), (X, w,)') = KA,e( (W,, w!) I1y:: KA (x),x)). Certain hyperparame- ters of the kernels can then be used to control the variability of the function with respect to specific inputs. We denote these hyperparameters as O with 0 C A and all other hyper-
4.4. Score
As stated in Section 2.4, our decision rule is to pick the causal graph with the highest posterior probability. For this, we need to calculate the log marginal likelihood for every variable in our model (Equation (6)). To find the log marginal likelihood, we must integrate over the priors of f and W, log p(x, (r) o,, O., oz) is given as
parameters as o A \ O. Specifically
log / [ pix./fe
0f. ax,
XPAgA Xi) Wi, o, Oc)p(( (w; wdw)dwdddd dw,
0 0
O.
Calculating this marginal likelihood is analytically in- tractable due to the non-linear dependence of f, on the latent term w. We thus use variational inference, with variational posteriors q. to optimise a tractable lower bound to the marginal likelihood (details and derivation
Hence, if the kernel hyperparameter value for Oij is near zero, the function for variable X; is constant with respect to
Ouur model can be thought of as Bayesian conditional varia- tional autoencoder.
5

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
as well handling of other hyperparameters are in Ap- pendix C.2). We denote the full lower bound for all variables as LELBO(q, 0, o, )). The optimisation problem is
Brouillard al., 2020; Ke al., 2020). Since we assume access only to observational data, we focus on the former.
Continuous optimisation methods cast causal discovery as single unconstrained optimisation problem (Zheng et al., 20. NOTEARS (Zheng et al., 2018) first introduced a regulariser for learning acyclic causal graphs in linear causal models. uollow-up work extended these methods to more flexible models such as neural networks, or graph neural networks (Lachapelle et al., 2019; Yu et al., 2019; Zheng et al., 2020). At the same time, works such as NOBEARS (Lee et al., 2019) and SDCD (Nazarct et al., 2024) have worked on improving the scalability of the acyclic regularis- ers. However, the above methods rely on restricted models to recover a DAG, which can limit performance when as- sumptions are not met (Section 2.3). When these models are not restricted, their guarantees no longer hold.
max LELBO 0 or, +) log p(O) s.t. G DAGS. 0,0.4.4
Maximising o, Φ, and gives accurate estimate of the marginal likelihood while maximising & finds the graph that maximises the posterior probability.
To enforce the condition that G^ E DAGS we use the spectral acyclicity constraint method (Lee et al., 2019)
((A))) = I)a(Aa))
(12)
where (A())) denotes the largest eigenvalue of Ae. which has been shown to be more stable and scalable than other acyclic constraints (Nazaret et al., 2024). Following Nazaret et al. (2024), a penalty method is used to allow for continu- ous optimisation, making the final loss
Other methods first search for node orderings and then prune to ensure a DAG structure. CAM (Bühlmann et al., 2014) uses sparse regression to learn node neighbourhoods and searches for orderings that maximise a score. SCORE (Rol- land et al., 2022) learns the ordering based on the score function of ANM models with Gaussian noise and then prunes using sparse regression. NOGAM (Montagna et al., 2023) extends SCORE to arbitrary noise distributions. The above methods specifically rely on the ANM assumption.
LALBO(y, O, a (),)))= +ogr wh((A))
(13)
where the weighting "/t is increased by P at each epoch. A solution A is found when h(A') < T for some T > O. See Appendix A.2 for more details.
4.5. Computational Cost
The Bayesian framework has previously been used for causal discovery but has typically focused on computing a posterior over causal structures to quantify uncertainty (Friedman & Koller, 2003; Heckerman et al., 2006). Geiger & Heckerman (2002) use Bayesian linear Gaussian mod- els but construct their model such that graphs in the same Markov equivalence class receive the same score. Similarly, Cundy et al. (2021) also only consider posteriors over lin- ear Gaussian models. DiBS (Lorch et al., 2021) assumes Bayesian ANMs and parametrises the space of graphs using latent variables. Annadani al. (2023) also use sampling methods to sample from the posterior over ANM causal models, but use permutation based parametrisation of the graph (Charpentier et al., 2021; Yu et al., 2021).
The cost of computing the loss for N samples in a Gaussian Process is usually ( We use the uncollapsed induc- ing point formulation that introduces M < NN inducing points allowing for mini-batching (Appendix C.2)(Hessma et al., 2013). This changes the cost to (M*)). As we have D variables, the total cost is ((DM)). The cost of computing the acyclicity regulariser is O (D²2. Hence, the cost is dominated by M. The choice of M depends on the problem, with more inducing points leading to better ap- proximations (Bauer et al., 2016). However, it is possible to use LELBO to find the optimal number of inducing points (Burt et al., 2020). Gaussian processes, and Bayesian meth- ods in general, can be more computationally expensive than their non- Bayesian counterparts. However, Bayesian model iidection-based causal discovery provides key advantages, including the ability to use more flexible model classes.
Previous Bayesian methods have also relied on the restric- tive assumptions common in likelihood based methods. In contrast, our work demonstrates that the advantage of Bayesian model selection is the allowance of fewer restric- tions compared to likelihood-basse methods (Stegle et al., 2010; Dhir et al., 2024; Cooper & Herskovits, 1992). This is strongly related to the idea of using independence of causal mechanisms (and the Kolmogorov complexity of causal models) to distinguish causal models (Janzing & Schölkopf, 2010) (see Dhir et al. (2024, Appendix B)). While previous studies demonstrated this with two variables, we show its feasibility and advantages in the important multi variable regime.
5. Related Work
Certain approaches only use information about indepen- dences between variables to discover the causal structure, recovering only an MEC, and not the whole DAG (Spirtes et al.. 2000; 1995; Chickering, 2002; Huang et al., 2018). Identifying a unique DAG requires additional assumptions, typically involving restrictions on model classes or access to interventional data (Lippe et al., 2021; Ke et al., 2022;
6

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
SIID | for 3 Variable
SID for 3 Variable
F1 for 3 Variable
1.00
0.75
4
0.50
2
2
0.25
0.00
NOTEARS
SDCD
CAM NOGAM SCORE
DiBS CGP CDE DCP
CDE
Random
NOTEARS
SDCD
CAM NOGAM
SCORE DIBS
CDE
CDE
Random
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBs CGP-
CDE DGP CDE Random
CGP DCP
Figure 1: Boxplot of metrics for the three variable dataset. The plots show the distribution of the metrics across the different datasets for each experiment. The centre line is the median, the black dot is the mean, and the white dots represent outliers. Lower is better for SHD and SID, and higher is better for F1.
6. Experiments
the ground truth adjacency matrix. SID: interventional distance counts the number of interventional distributions that are incorrect if the predicted graph is used to form the parent adjustment set instead of the ground truth graph (Peters & Bühlmann, 2015). F1: The F1 score is the harmonic mean of the precision and recall, where an edge is considered the positive class.
The structural
We compare our approach to various baselines that perform multivariate causal discovery to recover a DAG. We test our method on synthetic data generated from our model (Section 6. I). and data not generated from our model (Sec- tion 6.2). Then we test model on common semi- synthetic benchmark (Section 6.3). Additional experi- ments and results are in Appendix I. See Appendix C.6 for full details of the CGP-CDE implementation. Code can be found at: https: / //gtthbg.thh/b///e com/An AAiiss14ii/////// Cont inntinnoosssSSSSS ructurreerrttt gitt.
6.1. Synthetic 3 Variables
Using more flexible estimators may introduce some proba- bility error (Section 3). In the two variable case, Dhir et al. (2024) showed that the probability of error for our model class is low; we show that this is true for the multivariate case as well. To estimate the error only due to overlap in posteriors of different causal models, we perform discrete model comparison (labelled DGP-CDE, see Appendix D. for more details). This involves enumerating and computing the posterior of every causal model. We study the 3 variable case, where it is possible to enumerate every possible causal structure, 25 in total. Here, as required by the probability of error, data is generated from our GP-CDE model (details in Appendix F). We generate five datasets of 1000 samples for each of the six graphs that are unique up to a permutation of the variables.
was recently found that causal discovery methods can be sensitive to variance information that is an artefact of synthetic causal datasets (Reisach et al., 2021; Ormaniec et al., 2024). Hence we standardise every variable, including during data generation (see Appendix E).
Details for all baselines are in Appendix G. We compare against non- Bon-Bayeiian counterparts that use acyclic- ity regularisers, such as NOTEARS (Zheng al., 2018), and SDCD (Nazaret et al., 2024). SDCD uses the same spectral acyclicity regulariser as the CGP-CDE, the only difference being that the model uses a neural network and a likelihood score to learn the data distribution (Lachapelle et al., 2019). Thus, comparison against this method provides direct evi- dence for the usefulness of the Bayesian approach. We also compare against methods that first estimate topological order of the graph and then prune to get the final graph SCORE (Rolland et al., 2022), NOGAM (Montagna et al., 2023), and CAM (Bühlmann et al., 2014). These all assume ANM. Finally, we compare against the MAP estimate of the Bayesian model DiBS. Other Bayesian methods are sample based and hence do not allow for easy MAP calculation.
Baselines:
The results for the 3 variable case can be seen in Figure 1. Here, we can see that the discrete comparison DGP-CDE achieves very good performance. This shows that in the GP-CDE model, the probability of error in the multivariable case is low and that our inference does not contribute greatly to errors. Thus, if the GP-CDE model is good descriptor of the dataset at hand, we can expect the MAP value to identify the true causal direction with high probability. We also compare our continuous relaxation model against the discrete model that enumerates each graph type. In Figure 1, the performance of CGP-CDE shows that we can expect higher error due to optimisation of the continuous relaxation
Metrics: We use the following metrics to compare the methods. SHD: The structural Hamming distance is the Hamming distance between the predicted adjacency and
7

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
SHD J for 50 Variable ER4
SID for 50 Variable ER4
F1 t for 50 Variable ER4
800
2250
600
2000
0.6
400
1750
0.4
1500
200
1250
0.2
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM SCORE
DiBS
CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
CGP
(8)
SHD for Syntren
SID | for Syntren
for Syntren
300
100
0.4
80
200
60
0.2
100
40
0.0
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CGP-CDE CDE
(b)
Figure 2: Box plots of metrics for different graph types: (a) Erdos-Renyi (ER) graphs with four expected edges per variable, (b) Syntren data. Black dot is the average metric. Lower is better for SHD and SID, higher is better for F1.
in the CGP-CDE. However, the CGP-CDE still performs competitively compared to other baselines.
datasets in Appendix I.1.
The full results are in Appendix I, and we show the 50 vari- able ER4 results in Figure 2a. Compared to the DGP-CDE results of Section 6.1, we can expect higher number of errors due to the continuous relaxation and the shift in the data generation process. Nevertheless, the CGP-CDE out- performs all other baselines, especially in terms of SHD and F1 scores, and the increase in performance compared to the baselines greater for these larger graphs. While both SDCD and CGP-CDE employ flexible function ap- proximators and the same acyclic regulariser, the superior performance of CGP-CDE suggests that Bayesian methods are more effective at distinguishing causal structures. In con- trast, the poor performance of DiBS (MAP estimation) can be attributed to poor inference in Bayesian neural networks.
6.2. Synthetic Higher Variable
Next, we analyse the performance on a higher number of variables with varying graph density. Here, performing discrete model comparison is too costly as the number of DAGS is prohibitively large, laxation, labelled CGP-CDE. We generate five datasets for each graph size 50 using Erdos-Rényi (ER) and scale-free (SF) (Barabási & Albert, 1999) graph generation schemes, with expected total edge counts of 1D and 4D, respectively. In ER graphs, each edge is generated independently, while SF graphs promote few nodes with high number of edges. To test our claim that Bayesian model selection allows for more flexible models, ensure that the data is not gener- ated from our model. The data is generated from randomly initialised neural networks with noise included as an input. Performing well in such cases showcases the advantage of the added flexibility in the Bayesian approach. Data genera- tion details are in Appendix F. To show that our method also performs well with lower number of variables and samples, show results for graph size of 20 in Appendix 1.3 and Appendix 1.2. We also perform ablation study with ANM
continuous
6.3. Syntren
Syntren is gene-regulatory network simulator that gener- ates gene expression data from real biological graphs. There are datasets of 20 nodes with 500 samples each. The results on Syntren can be seen in Figure 2b. Here, the CGP. CDE outperforms other methods, especially in terms of the SID and Fl scores.
8

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
7. Conclusion
interventional data. Advances in Neural Information Pro- cessing Systems, 33, 2020.
In this work, we have shown why Bayesian model selec- tion works, and that highly effective approach for causal discovery multivariable settings. Then, we pro- pose scalable model that outperforms existing Bayesian and non-Ban Ba-eBBay methods across multiple settings. The consistent performance of our approach highlights the ad- vantages of framing causal discovery as Bayesian model selection problem. We believe this perspective addresses a critical limitation in causal cereovery-naeelyy the inability use flexible models, paving the way for more practical and widely applicable causal inference in real-world scenarios.
Bühlmann, P., Peters, J., and Ernest, J. Cam: Causal additive models, high- digeedimensional order search and penalized regression. The Annals of Statistics, 42(6), 2014.
Burt, D. R., Rasmussen, C. E., and Van Der Wilk, M. Con- vergence of sparse variational inference in gaussian pro- cesses regression. Journal of Machine Learning Research, 21(131), 2020.
Charpentier, B., Kibler, S., and Günnemann, S. Differ- entiable dag sampling. In International Conference on Learning Representations, 2021.
Impact Statement
Chickering, D. M. Optimal structure identification with greedy search. Journal of machine learning research, (Nov), 2002.
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none of which we feel must be specifically highlighted here.
Cooper, G. and Herskovits, E. A bayesian method for the induction of probabilistic networks from data. Machine learning, 9, 1992.
Acknowledgements
R.S. was supported by the Wellcome Trust [222836/2/21/7)] B.G. acknowledges the support of the UKRI AI programme, and the Engineering and Physical Sciences Research Coun- cil, for CHAI- EPSRC Causality in Healthcare AI Hub (grant EP/Y028856/1). A.K. is supported by UKRI (grant number EP/S023356/1) for Doctoral Training in Safe and Trusted AI. A.D. acknowl- edges the support of G-Research's 's March 2025 grant.
Cundy, C., Grover, A., and Ermon, S. Bed nets: Scalable variational approaches for bayesian causal discovery. Ad- vances in Neural Information Processing Systems, 34, 2021.
part of the UKRI Centre
Cvetković, D. M., Doob, M., and Sachs, H. Spectra of graphs: Theory and application. Academic Press, 1980.
Dhir, A., Power, S., and Van Der Wilk, M. Bivariate causal discovery using Bayesian model selection. In Proceed- ings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learn- ing Research. PMLR, 2024.
References
Aldrich, J. Autonomy. Oxford Economic Papers, 41(1), 1989.
Annadani, Y., Pawlowski, N., Jennings, J., Bauer, S., Zhang, C., and Gong, W. BayesDAG: Gradient-based posterior inference for causal discovery. In Thirth-seventh Confer- ence on Neural Information Processing Systems, 2023.
Dutordoir, V., Salimbeni, H., Hensman, J., and Deisen- roth, M. Gaussian process conditional density estimation. Advances in neural information processing systems, 31, 2018.
Barabási, A.-L. and Albert, R. Emergence of scaling in random networks. science, 286(5439), 1999.
Eggeling, R., Viinikka, J., Vuoksenmaa, A., and Koivisto, M. On structure priors for learning bayesian networks. In The 22nd International Conference on Artificial Intelligence and Statistics. PMLR, 2019.
Bauer, M., Van der Wilk, M., and Rasmussen, C. E. Under- standing probabilistic sparse gaussian process approxima- tions. Advances in neural information processing systems, 29, 2016.
Erdos, P., Rényi, A., et al. On the evolution of random graphs. Publ. math. inst. hung. acad. sci, 5(1):17-60. 1960.
Bello, K., Aragam, B., and Ravikumar, P. Dagma: Learn- ing dags via m-matrices and a log-determinant acyclicity characterization. Advances in Neural Information Pro- cessing Systems, 35:8226-822 82396,82222222 2022.
Feuerriegel, S., Frauen, D., Melnychuk, V., Schweisthal, J., Hess, K., Curth, A., Bauer, S., Kilbertus, N., Kohane, S.S. and van der Schaar, M. Causal machine learning for predicting treatment outcomes. Nature Medicine, 30(4): 958-968, 2024.
Brouillard, P., Lachapelle, S., Lacoste, A., Lacoste-Julien S., and Drouin, A. Differentiable causal discovery from

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
Friedman, J. H. Another approach to polychotomous classifi- cation. Technical Report, Statistics Department, Stanford University, 1996.
Janzing, D. and Schölkopf, B. Causal inference using the algorithmic markov condition. IEEE Transactions on Information Theory, 56(10), 2010.
Friedman, N. and Koller, D. Being bayesian about network structure. a bayesian approach to structure discovery in bayesian networks. Machine learning, 50, 2003.
Kass, R. and Raftery, A. E. Bayes factors. Journal of the american statistical association, 90(430), 1995.
Ke, N. R., Bilaniuk, O., Goyal, A., Bauer, S., Larochelle, H., Schölkopf, B., Mozer, M. C., Pal, C., and Bengio, Y. Learning neural causal models from unknown interven- tions. stat, 1050:23, 2020.
Friedman, N. and Nachman, I. Gaussian process networks. In Proceedings of the Sixteenth conference on Uncertainty in artificial intelligence, 2000.
Ke, N. R., Chiappa, S., Wang. X., Bornschein, J., Goyal, A., Rey, M., Weber, T., Botvinick, M., Mozer, M. C., and Rezende, D. J. Learning to induce causal structure. In International Conference Learning Representations, 2022.
Geiger, D. and Heckerman, D. Parameter priors for directed acyclic graphical models and the characterization of sev- eral probability distributions. The Annals of Statistics, 30 (5), 2002.
Guyon, I., Goudet, O., and Kalainathan, D. Evaluation methods of cause-effect pairs. In Cause Effect Pairs in Machine Learning, chapter Springer, 2019.
Kingma, D. P. and Ba, Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6988, 2014.
Lachapelle, S., Brouillard, P., Deleu, T., and Lacoste-Julien, S. Gradient- based neural dag learning. In International Conference on Learning Representations, 2019.
He, K., Zhang, X., Ren, S., and Sun, J. Delving deep into rectifiers: Surpassing human-level performance imagenet classification. In Proceedings of the IEEE inter- national conference on computer vision, 2015.
Lalchand, V., Ravuri, A., and Lawrence, N. D. Generalised gplvm with stochastic variational inference. In Proceed- ings of The 25th International Conference on Artificial Intelligence and Statistics, 2022.
Heckerman, D.
A bayesian approach to learning causal networks. In Proceedings of the Eleventh conference on Uncertainty in artificial intelligence, pp. 285-295, 1995.
Lee, H.-C., Danieletto, M., Miotto, R., Cherng, S. T., and Dudley, T. Scaling structural learning with no-bears to infer causal transcriptome networks. In Pucific Sympo- sium Biocomputing 2020, pp. tific, 2019.
Heckerman, D., Meek, C., and Cooper, G. A bayesian approach to causal discovery. Innovations in Machine Learning: Theory and Applications, 2006.
World Scien-
Hensman, J., Fusi, N., and Lawrence, N. D. Gaussian processes for big data. In Proceedings of the Twenty- Ninth Conference on Uncertainty in Artificial Intelligence, 2013.
Li, H., Wu, K., Ruan, C., Pan, J., Wang, Y., and Long, H. Cost-reduction strategies in massive genomics experi- ments. Marine Life Science & Technology. 1, 2019.
Hicks, al. Causality in economics. Australian National University Press, 1980.
Lippe, P., Cohen, T., and Gavves, Efficient neural causal discovery without acyclicity constraints. In International Conference on Learning Representations, 2021.
Hoyer, P., Janzing, D., Mooij, J. M., Peters, J., and Schölkopf, B. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21, 2008
Lopez, R., Hütter, J.-C., Pritchard, J., and Regev, A. Large- scale differentiable causal discovery of factor graphs. Ad- vances Neural Information Processing Systems, 35: 19290-19303, 2022.
Huang, B., Zhang, K., Lin, Y., Schölkopf,, B., and Glymour, C. Generalized score functions for causal discovery. In Proceedings of the 24th ACM SIGKDD international con- ference knowledge discovery & data mining, 2018.
Lorch, I.., Rothfuss, J., Schölkopf, B., and Krause, Dibs: Differentiable bayesian structure learning. Advances in Neural Information Processing Systems, 34, 2021.
MacKay, D. J. Comparison of approximate methods for handling hyperparameters. Neural computation, 11(5), 1999.
Immer, A., Schultheiss, C., Vogt, J. E., Schölkopf, B., Bühlmann, P., and Marx, A. On the identifiability and estimation of causal location-scale noise models. In- ternational Conference on Machine Learning. PMLR, 2023.
MacKay, D. Information theory, inference and learning algorithms. Cambridge university press, 2003.
10

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
Magnus, J. R. On differentiating eigenvalues and eigenvec- tors. Econometric theory, pp. 179-191, 1985.
Spirtes, P., Meek, C., and Richardson, T. Causal inference in the presence of latent variables and selection bias. In Proceedings of the Eleventh conference Uncertainty in artificial intelligence, 1995.
Montagna, F., Noceti, N., Rosasco, L., Zhang, K., and Lo- catello, F. Causal discovery with score matching addi- tive models with arbitrary noise. In Conference on Causal Learning and Reasoning. PMLR, 2023.
Spirtes, P., Glymour, C. N., and Scheines, R. Causation, prediction, and search. 2000.
Nazaret, A., Hong, J., Azizi, E., and Blei, D. Stable dif- ferentiable causal discovery. In Forty-first International Conference on Machine Learning, 2024.
Stegle, O., Janzing, D., Zhang, K., Mooij, J. M., and Schölkopf, B. Probabilistic latent variable models for distinguishing between cause and effect. Advances in neural information processing systems, 23, 2010.
Ormaniec, W.. Sussex, S., Lorch, L., Schölkopf, B., and Krause, A. Standardizing structural causal models. arXiv preprint arXiv:2406. 11601, 2024.
Titsias, M. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics. PMLR, 2009.
Pearl, J. Causality. Cambridge university press, 2009.
Titsias, M. and Lawrence, N. D. Bayesian gaussian process latent variable model. In Proceedings of the thirteenth international conference on artificial intelligence and statistics. JMLR Workshop and Conference Proceedings, 2010.
Peters, J. and Bühlmann, P. Identifiability of gaussian structural equation models with equal error variances. Biometrika, 101(1), 2014.
Peters, J. and Bühlmann, P. Structural intervention distance for evaluating causal graphs. Neural computation, 27(3), 2015.
Van den Bulcke, T., Van Leemput, K., Naudts, B., van Re- mortel, P., Ma, H., Verschoren, A., De Moor, B., and Marchal, K. Syntren: generator of synthetic gene pression data for design and analysis of structure learning algorithms. BMC bioinformatics, 7, 2006.
Peters, J., Janzing, D., and Schölkopf, B. Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
Rasmussen, C. and Ghahramani, Z. Occam's razor. Ad- vances in neural information processing systems, 13, 2000.
Walker, S.. Damien, P., and Lenk, P. On priors with kullback-leibler property. Journal of the American Sta- tistical Association, 99(466), 2004.
Rasmussen, C. E. Gaussian processes in machine learning. In Summer school on machine learning. Springer, 2003.
Walker, S. G. Modern bayesian asymptoties. Statistical Science, 2004.
Reisach, A., Seiler, C., and Weichwald, S. Beware of the simulated dag! causal discovery benchmarks may be easy to game. Advances in Neural Information Processing Systems, 34, 2021.
Wei, D., Gao, T., and Yu, Y. Dags with fears: A closer look at continuous optimization for learning bayesian networks. Advances in Neural Information Processing Systems, 33:3895 3906, 2020.
Rolland, P., Cevher, V., Kleindessner. M., Russell, C., Janz- ing, D., Schölkopf. B., and Locatello. F. Score matching enables causal discovery of nonlinear additive noise mod- els. In International Conference on Machine Leurning. PMLR. 2022.
Williams, C. and Rasmussen, C. Gaussian processes for regression. Advances in neural information processing systems, 8, 1995.
Williams, C. K. and Rasmussen, C. E. Gaussian processes
for machine learning, volume 2. MIT press Cambridge, MA, 2006.
Sachs, K., Perez, O., Pe'er er, D., Lauffenburger, D. A., and Nolan, G. P. Causal protein- siggiinnggn networks derived from multiparameter single-cell data. Science, 308(1))1), 2005.
Yu, Y., Chen, J., Gao, T., and Yu, M. Dag-gnn: Dag structure learning with graph neural networks. In International Conference on Machine Learning. PMLR, 2019.
Salimbeni, H., Eleftheriadis, S., and Hensman, J. Natural gradients in practice: Non-conjugate variational inference in gaussian process models. In International Conference on Artificial Intelligence and Statisties. PMLR, 2018.
Yu, Y., Gao, T., Yin, N., and Ji, Q. Dags with no curl: An efficient dag structure learning approach. In International Conference on Machine Learning. Pmlr, 2021.

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
Zhang, K., Wang, Z., Zhang, J., and Schölkopf, B. On estimation of functional causal models: general sults and application to the post-nonlinear causal model. ACM Transactions on Intelligent Systems and Technology (TIST), 7(2), 2015.
Zheng, X., Aragam, B., Ravikumar, P. K., and Xing, F. P. Dags with no tears: Continuous optimization for structure learning. Advances in Neural Information Processing Systems, 31, 2018.
Zheng. X., Dan, C., Aragam, B., Ravikumar, P., and Xing, E. Learning sparse nonparametric dags. In Interna- tional Conference on Artificial Intelligence and Statistics. PMLR, 2020.
12

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
A. Additional Background
A.1. Causal Model Basics
The causal model satisfies the Independent Causal Mechanism (ICM) assumption by construction (Janzing & Schölkopf, 2010). Given a graph 9, this means that changing term in the causal factorisation of the joint I1 (i(X) does not change any of the other terms. This is achieved by assuming that the choice of conditional distributions for each variable X; from (1) is independent of the choice for other variables Xj from C))PAc(s) with j # i. This property is assumed to hold in the causal factorisation but does not necessarily hold in any other factorisation of the joint (Peters et al., 2017), This can be seen in the two variable system with graph X + Y. If the terms in ((YI)) are chosen independently, the choice of the terms in (())((X Y) are tied to each other. Changing P(X) or P(Y[X) in general will require changes in both P(Y) to maintain equality.
As variables are dependent only on their parents a graph by construction, a causal model satisfies the Markov assumption. This states that d-separations in & between variables implies conditional independence between the same variables (Pearl, 2009). We also assume faitffunnes:: any conditional independences between variables represented by d-separations in G (Pearl, 2009). Together, faithfulness and Markov assumptions imply a one-to-one relationship between independences in the distributions and d-separations in the graph (Peters et al., 2017). These assumptions are enough to recover a causal structure from data up to a Markov equivalence class (MEC), that is, the class of causal structures that have the same d-separations.
A.2. Continuous Optimisation for Learning Causal Structure
Given a score (such as the likelihood or marginal likelihood) S, causal structure learning can be defined as an optimisation problem over graphs (Zheng et al., 2018)
g argmax S(G) such that G E DAGs.
(14)
That is, finding the graph that optimises the score under the constraint that the graph must be a valid DAG. This converts causal structure learning into a single optimisation problem, tackling the scalability issue. However, ensuring that the graph acyclic is non-trivial with gradient based methods. Zheng et al. (2018) first solve this by encoding graph as a weighted adjacency matrix A E 12 D and using a measure of acyclicity h (AA) as constraint. The function h A - Eg0 is a measure of the number of weighted directed walks that allow for returning to starting node. Hence, h(A) = 0 implies that there are no directed cycles in the adjacency A and that the graph implied by A is valid DAG.
Numerous variations on this approach have been proposed to improve the scalability and stability of optimisation (Lee et al., 2019; Bello et al., 2022; Wei et al., 2020; Nazaret et al., 2024). We use the spectral acyclicity constraint method in Equation (12), which is equivalent to the largest eigenvalue magnitude of A, due to its stability and scalability (Lee al., 2019; Nazaret et al., 2024). This constraint uses the fact that A is acyclic if and only if all its eigenvalues are zero (Cvetković al., 1980).
The gradients of the spectral acyclicity constraint can be easily calculated as
Vh(A)=
(15)
U/ ud
where ua and Va the right and left eigenvectors associated with >a(A) (Magnus, 1985). We follow the example of Nazaret et al. (2024) and use the power iteration method to estimate these eigenvectors in O(d)) rather than finding them exactly, which costs (dd))
Equation (14) can now be written as an optimisation problem with the differentiable constraint h(A) = 0. Nazaret et al. (2024) use penalty method, writing the objective as:
argmax S(GA) -A),
(16)
A
where the weighting % of the penalty ((A) is increased by p each epoch. A solution A is found when h(A*) < T for some convergence threshold 0.
13

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
B. Proof of Distinguishability with Multiple Variables
Bayesian model selection uses the marginal likelihood (equivalent to the posterior over models) to distinguish between Bayesian causal models with different causal graphs. Thus, following Dhir et al. (2024), we will provide conditions such that two Bayesian causal models will have different distributions over marginal likelihood values. This implies that there will always exist datasets such that the marginal likelihood values of two Bayesian causal models differ. We will do this by defining two Bayesian causal models as Bayesian distribution-equiialent if they have the same distribution over marginal likelihood values, and then provide conditions such that two Bayesian causal models cannot be Bayesian distribution-equivallee Note that this can be viewed as a version of observation equivalence (Pearl, 2009), but for Bayesian causal models. If models are not observationally equivalent, they can be learnt from observational data (for example, by using methods such as PC (Spirtes et al., 2000)). A similar reasoning holds for Bayesian distribution equivalence.
We split the proof into two parts, first, we will show that Bayesian model selection can effectively identify up to a Markov equivalence class. This is known for certain classes of models, for example, curved exponential models (Chickering, 2002), but it has not been shown for the general non- naraaaaaeetric case. Second, we will provide conditions such that two Bayesian causal models within the same Markov equivalence class are not Bayesian distribution-equivllenn Here, there may be some probability of error depending on the overlap of the posteriors that can be computed empirically. This second result is implied by the result in Dhir et al. (2024, Proposition 4.7).
We define Bayesian distribution eeuiivaleee models. These are the Bayesian counterparts of models that are unidentifiable. That is, given any dataset, these models will achieve the same marginal likelihood score.
Definition B.1. (Dhir et al., 2024, Definition 4.4) Given two Bayesian causal models (Mg, mg). (Mx, Tx). say they are Bayesian distributiinn-quiialenn if P(- Mg) = P(- Mx), i.e. for all N E N, and for all (xx, ) ) holds that p(xN, y MM) = r(x y Mx)-
(x*, E (X x y) it
We proceed by showing that two Bayesian causal models with graphs in different MECs cannot be Bayesian distribution equivalent. Further, if the estimate of the Bayes factor (Kass & Raftery, 1995) (ratios of marginal likelihood) is consistent, it can identify the correct graph.
Lemma B.2. Assume faithfulness and two Bayesian causal models (Mg: ng), and (MH; πH) Ifg and H in separate Markov equivalence classes (MEC), then Fg O Fy. 0, where Fg and Fu are the sets of all joints expressible by Mg and MR respectively.
X; Lg X3 Xg, where 1p denotes distributional independence, and 1g signifies d-separation in graph G. Hence, as G and H are in separate MECs, there exists some X; dg Xj|X& such that X; K Nx. This implies that for every P E Fg. we have that X; Lp Xx, and for every Q E FH we have that X; La X,lXx. Hence P & FH.
Proof. The Markov and Faithfulness assumptions together means that X, Lp X,
Theorem B.3. Assume faithfulness and two Bayesian causal models (M₂, πG), and (M).TH). IfG and H are in separate Markov equivalence classes (MEC), then they are not Bayesian distribution equivalent for any choices of priors Tg, TH- Assume N samples are generated from one of the models. If the estimate P(X P(X* consistent as N + 00, the probability of error between these two models is O.
(MX))M
Proof. The first claim follows directly from Lemma B.2 and Dhir et al. (2024, Proposition 4.5). Given that Fg nF = 0. |
only one of the models can be the true model. Consistency implies P(**MMx) - oo if the data generating distribution
P [Ms) P E Fg. and P(XN [Mx) + 0 if the data generating distribution P E Fu (Walker, 2004). Hence, the marginal likelihood is higher for the correct model for each distribution in the union of Fg Fx leading to probability of error of O.
Consistency in the above requires certain conditions on the prior (Walker et al., 2004, Theorem 1). Suppose the data is generated from XN ~ Po. we require that for the correct model, suppose it is Mg. the prior puts neg-nggle mass around Po: TT((P: [P] < e}) > 0 for all E > 0 (Walker et al., 2004). Second, the posterior for both models cannot paradoxically concentrate in a region of negligible prior mass: lim inf N [L Pv 2 c for all c > 0 where A = [P: : PL > e). PNA is the posterior (over N samples) computed from the prior but restricted to the set of densities A, which are e far from the true density Po as measured by the KI--divergence (Walker, 2004). This condition
14

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
stops the wrong model (which has negligible prior mass around the true density) from concentrating around the true density for large enough sample size.
The condition for two causal models that F& n Fxx = 0 also holds when the causal models are identifiable (Guyon et al., 2019, Ch. 2). Hence, Bayesian model selection can recover the causal structure if the causal models are identifiable.
Corollary B.4, Assume two Bayesian causal models (Mg. mg), and (Mu, Tn) where the undertying causal models are identifiable: Fg FH = 9. If the estimate P(X~ [Mg], consistent as N + 5o, the probability of error between these two models O.
((XNMM)
Next, we consider the case of two Bayesian causal models where the two graphs are in the same MEC. Without restrictions, it is well known that in this case there may exist distributions in two Bayesian causal models that can express given dataset (Pearl, 2009). Thus, the result from Lemma B.2 does not hold. The key idea that allows for distinction here is that the ICM assumption (as represented by factorised priors) may only hold in the causal factorisation of Bayesian causal model. Further, the distributions of marginal likelihood values are sensitive to the ICM assumption. As the causal factorisation of the two Bayesian causal models differs, they will have a different distribution over marginal likelihood values.
We first formalise the ICM assumption as the separability of the priors with respect a certain factorisation of the joint. Definition B.5. (Dhir et al., 2024, Definition 4.2) Given a Bayesian causal model (Mg, mg), a prior is separable with respect to Cg if it factorises liv T, such that T, E ((C)))A)()) where P is the set of all distributions over an object.
That is, if π is separable with respect to Cg. the joint can be written as ILev (X) (T) (dPP), with P. E Ci) Separability of the prior (and hence the ICM condition) may hold in multiple factorisations. This is the case with normalised linear Gaussian models (Dhir et al., 2024, Appendix D.2). To formalise this, we reiterate the notion of separability- compatibility.
Definition B.6. (Dhir et al., 2024, Definition 4.6) Given two Bayesian causal models (Mg, mg). and (Mx, Tx). Denote Cg - Cu bijection such that for any P E Cg such that Q = 8(P) E Cu, there holds an equality of joint measures The two Bayesian causal models are separable-compaible if: i) the
Y:Cg
Iliey P:(X:)XPAo(t) = Iliev pushforward mig o is separable with respect to Cy, ii) THe o y is separable with respect to Cg.
The operation ô is simply the transformation by Bayes' rule that transforms the factorisation according to G the factorisation according to H while ensuring the joint is the same. Separability separability-compatibiii compatibility checks whether the prior separates according to multiple factorisations. If it does, ICM holds in multiple factorisations.
Theorem B.7. Given two Bayesian causal models (Ma, πG). and (MH,TH). I and H are in the same Markov equivalence classes (MEC), and the Bayesian causal models are not separable compatible, they camnot be Bayesian
distribution-eqiivllenn Further, assuming data is generated from one of the models, then the probability of error is less than 50%.
Proof. The first statement follows directly from Dhir et al. (2024, Proposition 4.7). The probability of error in the two model case can be written as (Dhir al., 2024)
1
P(E) = 2 -1- (VPP((T)))
(17)
If the two models are not Bayesian distributionn-equivalenn, then (VP(( > 0. This implics that the probability of error is less than 1
2'
Theorem B.7 shows that if the ICM assumption only holds in the causal factorisation of the model (prior is only separable with respect to the causal factorisation), then there will always exist datasets that give different marginal likelihood values for the two Bayesian causal models. The above holds for any combination of Bayesian causal models with graphs in the same MEC. Hence, the total probability of error for the comparison within an MEC will be less than that of a graph randomly chosen from within the MEC (Friedman, 1996).
C. CGP-CDE Details
For our Bayesian causal model, we use the Gaussian process conditional density estimator (Dutordoir et al., 2018). We use a sum of commonly used kernels Rasmussen (2003, Ch. 5), as outlined in Appendix C.I. We also apply variational inference
15

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
optimise the lower bound, as is standard for latent variable Gaussian process models (Dutordoir al., 2018; Titsias & Lawrence, 2010) and use inducing points with stochastic variational inference to improve data scalability (Hensman et al., 2013). We outline the variational inference in Appendix C.2 and detail our optimisation schedule in Appendix C.3.
C.1. CGP-CDE Kernel
The ultimate aim of the approach is to find the causal graph that maximises the posterior probability. The choice of kernels will be problem dependent and it is possible to carry out Bayesian model selection over the kernels themselves (Rasmussen, 2003, Ch. 5). To allow for good fits of a range of datasets, we use sum of kernlls that can express functions with range of lengthscales and roughness well.
The kernel for function fę is
k; Kinn, +
Rage, + Km12, km32,i + kig.is
(18)
where each the kernels defined below (Williams & Rasmussen, 2006, Ch. 4). For all the kernels, 8₁ is the hyperparameter for the latent dimension.
i is the linear kernel:
Σ
iy + ww,dd, jtt
(19)
where fjn,ki indicates the jth dimension of hyperparameter in Esge, is the squared exponential kernel
(xg - x2y}
+ exp (_22 (wi w()² -
Rugei
022 oggo
Σ cxp j
-82 -0pp,i
(20)
2
where Aage is the precision parameter.
km12, and k are the Matérn 12 and Matérn32 kernels with // y== and v = respectively, with the general form
Kmu, o Σ 21-" (Dmu.d) 201z, (2mi.
- adi)) Kr (Omu.s V0mniyV2ule V2 rryl)
I(v)
141
21-v + mv, T(U) (Pmeits V 22/w: - 220(uu -ori: B, (0 - wil). (21)
where [(v) is the gamma functions and B, is a modified Bessel function. Similarly to the squared exponential kernel, Am
the precision parameter.
king is the rational quadratic kernel, which is equivalent the sum of many squared exponential kernels with different precision hyperparameters (Williams & Rasmussen, 2006, Ch. 4),
02 (x, x1)2
-a. + +1+ 02 (wi w))
-0;
kiq.i
Σ jæL
1
(22)
20d
2a;
where a a hyperparameter that learned.
The kernel variance terms o 2n and o mean individual kernels can be "switched off" by setting this term to zero they do not contribute improving the evidence lower bound.
In all these kernels, the hyperparameter denoted by 0 controls the variability of the function, discussed in subsection 4.2. This means as they tend to zero, the input and output of the function become decorrelated. To construct the adjacency matrix in Equation 10, we sum these hyperparameters, so that
0.j Oiinii Osgris + 0m1224 0m32,i5 + Omg
(23)
16

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
C.2. CGP-CDE Lower Bound
We use variant on the GP-CDE (Dutordoir et al., 2018) define our causal model and define prior over distributions. Intractability in these models is commonly solved using variational inference (Lalchand et al., 2022; Dutordoir et al., 2018).
Each variable x; is modelled as a function of all other variables X and latent variable wa plus some Gaussian noise, parameterised by noise variance o²
X; f((-i, w:) t e~ NO of).
(24)
The presence of the latent variable allows learning of non-Gaussian and heteroscedastic noise (Dutordoir et al., 2018). A Gaussian process prior placed on f.
p f XPA (1), (t), A,, w;) = (O (().,,, (X w,))),
(25)
where XPAça (6) are the parents of x; A, and K is the kernel with set of hyperparameters denoted A,. We denote kernel hyperparameters that control dependence and are included in the adjacency matrix collectively 0, while the rest o, with A {0, ).
Different prior beliefs of the form of f can be expressed in the choice of kernel. To find the causal DAG with the highest posterior, we need to maximise the marginal likelihood of each variable. This is written as
p(x; XPAaa(0)))))) (1): A. A.) = log / / ( (PA As, w.lp(w.)dw,dff.
However, the latent variables w, appear non-linearly in p(f, XPp (i) Ag, w.) making the calculation of (eapp(xXXPAAAAA PAAA (0), A,, Φ,) analytically intractable.
log p(x:
The addition of inducing points u, with corresponding inducing inputs Zi, helps both with the intractability of 1(x, XPA (f), A,, Φ,) and scaling to large datasets (Titsias, 2009; Titsias & Lawrence, 2010). We can write the joint between and u, as:
f. (()== ((.
0 (1)-
Kerr Kg,ta
N
(26)
Ki,uf
Ks,uu
where Ke is the covariance matrix between training data, Kim is the covariance matrix between inducing points and K and Kut are the covariance matrices between the two. Suppressing the hyperparameters and Z; for brevity, the marginal likelihood can then written as:
log (XAga(())))))) log /
PAgA (1), Wi, (wwp)/)))))),))dddddd..
(27)
where p(u:) (O, Ki,u)). This is still not tractable, we use variational inference to define an evidence lower bound (ELBO) to the marginal likelihood (Titsias, 2009). This is done by introducing a variational distribution qff;, w,, u.) and rewriting the marginal likelihood as (suppressing conditioning terms for neatness):
log p(x;) = logg / pix.lt, (f) u, XPA (0) w:)p(u,)p(w))
(fi,, Wi, U:) (ff,, w,, u.)
du,dw,df,
The variational distribution takes the form (Titsias & Lawrence, 2010; Dutordoir et al., 2018):
qff, Wi, u:) p(f/ww, u: (w))).
(28)
where q(u,) = (u,,,,,,,,... Sw.)) ) is the Gaussian variational distribution of inducing points u, and g(w:) = -N(,, E,) is
the variational distribution of w..
Rearranging and using Jensen's incquality, we can then get the lower bound for variable d as (Dutordoir et al., 2018)
log,(x;)> 2 ((logp(x;(f)) (f.))) g(w.)) KLLg(u)lp((u)] KL[g(w:))lpt [p(w:)).
(29)
17

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
where ((.)) - / p(f, u, w.)a(u)du. We denote the lower bound for variable X; as LELBO, (.). The full lower bound then is [q, A,, 0:) :- - LELBO,i(g; A., 0i), remembering A = =00;,00 o.).
By incorporating the optimal Gaussian variational distribution q(u), which can be calculated in closed form, it is possible to collapse the bound by integrating out the inducing variables (Titsias, 2009; Titsias& & Lawrence, 2010). However, the computational complexity of this approach is (N (NV D)D), meaning for large values of N training becomes infeasible (Lalchand et al., 2022). To allow method to scale to large datasets, and keep (u)) uncollapsed so that stochastic variational inference can be used. This means the data can be mini-batched in training, reducing the computational complexity to O(M To further improve computational efficiency, we use an encoder to learn a function 9q.i (Xhmch) + (1i,, Σi) for the variational distribution g(w,) = N(u,, E,), where Xbarch E x is batch of data points, 14 E 12'h and M, E IR for batch size b (Dutordoir et al., 2018). This is more efficient than the alternative of learning N parameters. We use multi-layer pereeptron for the encoder and list the hyperparameters
follow the approach of Hensman et al. (2013)
used in Table 1.
Table 1: Hyperparameters for the variational encoder 9q. : (Xnatch) + (Mi, E.).
Parameter
Value
hidden layer size
128
number of layers
activation function
ReLU
Other hyperparameters: We follow common practice and maximise the lower bound with respect to the hyperparameters a and Φ as well (Rasmussen, 2003) (along with 0 and q). This is justified based on the fact that the posterior for these hyperparameters tends to be highly peaked in practice (MacKay, 1999). Thus, a Laplace approximation can be approximated by the maximum value.
C.3. Optimisation schedule
We split the training schedule up into three steps. The warm-up phase allows the model to learn conditional independences present in the data. In the acyclic constraint phase, where a penalty method is used to enforce the DAG constraint, the model then learns the direction of any edges. The learnt DAG is then enforced, and the ELBO is optimised in the cool-down phase to find the ELBO of the final model, which can then be used for comparing random restarts.
Warm-up phase: The warm-up phase allows the model to learn any conditional independences from the data. It also ensures the hyperparameters and variational parameters of the Gaussian process are reasonable values at the beginning of the acyclic phase. Although the number of iterations needed will depend on the problem hand, we found To = 25, 000 iterations for this phase sufficient to ensure convergence.
Acyelic constraint phase: After the warm-up phase has concluded, a lot of the edges corresponding to conditional independences are switched off. However, for variables that are correlated, the causal direction has not been decided. This part of the training forces the constraint that the final graph has be acyelic. The loss is
a, LELBO 0 (a. x) = Lg. 0 の) + ((g)(0) 7ah(A),
(30)
optimised, using the acyelic regulariser proposed by Nazaret et al. (2024). The weighting 7. of the acyclic penalty term ((AA) increased linearly each epoch by p = 50, starting from 0. As we use minibatching, one epoch is full pass through the dataset. For some of the 50 variable datasets, is because the unconstrained loss scales with the number of variables, so for the 50 variable dataset, higher p ensures a higher ratio between the acyclic constraint and the rest of the loss, ensuring faster convergence.
increased this to p 250 to reduce the time required to train. This
We use 50 power iterations to approximate the eigenvectors required to calculate the gradients of the acyclic constraint, as described in Appendix A.2. The acyelic constraint phase is terminated when h(Ae.) < T where 7 = 0.005, or the number of iterations exceeds maximum value of Ta 50, 000.
18

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
Acyclic thresholding: After the end of the acyclic phase, the values of the hyperparameters are sometimes not exactly 0 due to numerical precision, find the adjacency matrix with extra thresholding step. Following Lachapelle et al. (2019), we threshold weights to 0, starting with the lowest weight, until the adjaceney is exactly acyclic (h(Ao,) 0).
Cool-down phase: We found that random restarts helped tackle local optima in optimisation Appendix C.4. To allow for the comparison of ELBO values across restarts (see Appendix C.4) we continue training the found acyclic DAG for Ty iterations. This ensured that the found ELBO approximates the ELBO of the final resultant DAG. As some hyperparameters may have been pushed to extreme values by the acyelic threshold if competing with an anti-causal edge, we reinitialise the hyperparameters of the edges that are still active to their initial value. This reduces the number of iterations needed for the cool-down phase. We found that f 25, 000 cool-down iterations sufficed here.
Final matrix thresholding: Finally, we perform an extra thresholding step to remove edges that do not encode correlations. As our data is normalised, we check the kernel evaluation of points -3 and which should bound 95% of the data. If the function is highly correlated between these points, it is effectively constant for 95% of the data. We thus threshold the linear variances to 10-4 and the graph parameters (the rest of the kernels) to 0.05, which corresponds to lengthscale of 20. These values ensure that points between -3,3 are highly correlated.
C.4. Random Restarts
The main principle behind our method is to find the graph that maximises the marginal likelihood (or a lower bound to it). Maximising the ELBO with respect to the hyperparameters A in Gaussian process models is known to suffer from local optima issues (Rasmussen, 2003, Ch. 5). The final hyperparameters found, and hence the adjacency matrix, can be dependent on the initialisation of the hyperparameters. To mitigate this, we use a widely adopted technique of performing random restarts (Dhir al., 2024). We optimise the loss starting from multiple initialisations (this can be done in parallel) N, times. Then we pick the graph that achieves the highest CELBO as the candidate for the most likely graph. Adding more random restarts for the CGP-CDE should improve performance, as shown for the DGP-CDE Appendix D.I.
Algorithm 1 Optimisation procedure for the Causal GP-CDE.
Input: Data X, number of random restarts N,, acyclic penalty weighting update p. threshold т, acyclic penalty weighting = 0 initial A ()
Result: Most likely adjacency matrix A*
Initialise empty list graphscore
for i to N, do
for j = I Update A by maximising L (Equation (30))
do
end
1 0-0
while h(Aa,) and < Ta do Update A by maximising L (Equation (30))
1 1- 1+1
epoch ended then
" x- 7t +
end
end
end
Set "t. 0 Threshold Ae, starting from lowest weight until h(Ao.)
for k = to T' do I Update q, A by maximising L (Equation (30))
end
Append (Ag, L) to graphscore.
A +- Ag from graphscore with the maximum L.
19

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
C.5. Hyperparameter Priors
As discussed in Section 4.3, we can place priors on the graphs by placing priors on the hyperparameters 0 of the CGP- CDE. For our implementation,
place a prior 0 that favours sparser graphs. Specifically, we use the Gamma prior P(0) Gamma(n, B), where 7) is the shape parameter and 3 is the rate parameter. For all experiments, we set 17 1 and B 10. Our hyperparameters for the priors were heuristically chosen and give log probabilities that are a fraction of the rest of the loss.
C.6. Implementation Details
In this section, we outline the details of our implementation.
Hy pyperaaraeeeer initialisations We initialise the hyperparameters 0 ~ Uniform(0.00 1), except Ojin = 0.25. All kernel variances o are initialised to so they all have the same initial weighting and likelihood variance initialised as Ø² 6² where K ~ Uniform(50, 100). The a term in the rational quadratic kernel is initialised as a ~ Uniform(0.1, 10).
Variational parameters We use 400 inducing points for all datasets as we found this reasonable trade off between computation time and accuracy. We initialise our inducing point locations at subset of the data inputs. We use a multilayer perception for the latent variable encoder, as described in Appendix C.2, with hyperparameters listed in Table 1. The encoder weights are initialised from a truncated normal distribution centred on 0 with a standard deviation of V hidden layer size (He al., 2015).
Loss calculation We use a minibatch size of 256 to calculate the loss. We take 50 Monte Carlo samples to integrate out the latent variables w, N(s, Ex).
Acyelic constraint penalty The acyelic constraint penalty has three hyperparameters: the coefficient to linearly increase the weighting of the penalty term by, p. the threshold value for the acyclicity, т, and the number of power iterations used to approximate the eigenvectors for the gradients of the acyclic constraint. We tuned these variables using five datasets generated from 10 variable ER1.5 graphs. We chose 50 power iterations, as it led to better performance than smaller numbers of iterations and after 50 we saw minimal improvement. For the 10 variable ER1.5 datasets, we found p = 50 and r = 0.005 leads to the best performance across SHD, SID and F1, while still converging in a reasonable amount of time. However, as the number of variables scaled, this value of p took prohibitively long, except for the Syntren dataset. This is because the ELBO scales with the number of variables, LELBO A,, di) := E CELBO. (q: A., Øz). Therefore, for the 3, 20 and 50 variable datasets we scale P with the number of variables, such that p = 5D. This ensures the relative magnitude of the ELBO and the acyelic penalty remains roughly the same, meaning the time for the acyelic constraint phase to converge is not dependent on the number of variables. We did not find scaling affected the performance, but it did speed up the acyclic phase.
Optimisers Each optimisation step consists of two steps. First, optimise the parameters of the variational distribution of inducing points q(u;) = (N [ Su.;) using natural gradients, which have been shown to significantly improve training time for variational Gaussian processes (Salimbeni et al., 2018). We linearly increase the natural gradient step size from 0.0001 to 0.1 for the first five iterations, and then use a step size of 0.1. Second, we optimise the Gaussian process hyperparameters using Adam with learning rate of 0.05 (Kingma Ba, 2014).
Random restarts We did as many random restarts as we had computational resources for. This means we did restarts for the three-variable dataset, restarts for each of the 50 variable datasets and the Syntren dataset, and restart for the 20 variable datasets.
D. DGP-CDE Details
In Section 6.1 we introduce the discrete Gaussian process conditional density estimator (DGP-CDE). which is used to enumerate every possible causal structure for the 3 variable case. For this, a separate DGP-CDE is trained for each of the 25 possible causal graphs, with the Gaussian process prior for each variable defined in Equation (9). The DGP-CDD differs from the CGP-CDE in that only takes Xp) as inputs, whereas the CGP-CDE takes X as inputs and then learns the
20

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
adjacency matrix.
As the adjacency matrix is fixed for the DGP-CDE there is no acyclic penalty, the loss function is just the lower bound Bo. A, Φ) := Σ CFLBO.s((5 A,, o.). When all 25 separate DGP-CDE models have been trained, the most likely graph is selected by choosing the model with the highest marginal likelihood.
D.I. Implementation Details
Loss Calculation As the DGP-CDE is limited in the number of variables it can model by the need to enumerate over all graphs there is no need to make the lower bound scalable using the methods discussed in Appendix C.2 (Hensman et al., 2013). Instead, we use the collapsed version of the lower bound, which can be calculated in closed form (Titsias, 2009).
Kernels We use a sum of a linear kernel defined in Equation (19) and a squared exponential kernel defined in Equation (20). This kernel means the kernel expectations needed to calculate the lower bound can be computed in closed form.
The kernel variances are initialised as o, = 1, the likelihood variance is randomly sampled as of ~ (00formm100+11000 and the precision parameter 0, Uniform(1,100)
Hyperparameters
Variational parameters As the enumeration over graphs is computationally expensive, we use 200 inducing points, which we found be sufficient. As we are not using stochastic variational inference (Hensman et al., 2013) for the DGP-CDE, we can no longer mini-batch, so don't 't use multilayer perceptron for the latent variables. Instead, we initialise the latent variable mean μ4 0.1x, and standard deviation is randomly sampled Σ, Uniform(O, 0.1).
Optimisation Schedule The use of the collapsed bound for the DGP-CDE requires different optimisation schedule than for the CGP-CDE. This consists of a two part optimisation scheme. Due to the loss function being highly non-convex and suffering from local optima, the first step of the optimisation scheme uses Adam (Kingma & Ba, 2014) with a learning rate of 0.05 to get into a good region of the decision space. Once the loss function reaches the value a noise model would have, the optimisation scheme switches to the FeaueennFletcher-GGllaaaatt for the final part of the optimisation.
(BFGS) algorithm to perform gradient descent
Random Restarts We did 10 random restarts for each DGP-CDE, selecting the one with the best lower bound for our final result. We analysed the effect of the random restarts by permuting the order of the random restarts and plotting, in Figure 3, the mean and standard deviation of the SHD, SID and Fl across permutations the number of restarts increases. From the plot, it is clear that the number of random restarts increases, the mean value of the metrics, average, improves. This is because more random restarts allow more thorough exploration of the loss function. It also shows that higher value of LELBO does lead to better causal structure discovery.
E. Data Standardisation
For additive noise models, it has been shown that when data is sampled from simulated DAG and not normalised, the marginal variance tends to increase along the causal order (Reisach et al., 2021). When this is the case, it is possible to, least partially, determine the causal order by ranking the marginal variance of variables. Reisach et al. (Reisach et al., 2021) demonstrate that commonly used synthetic datasets often have this property, meaning a simple baseline method based on variance sorting and regression can perform as well as some recent continuous structure learning algorithms. This variance can be manipulated by rescaling (such as changing measurement units), which would then give a different causal order. As such, results that take advantage of the variance may be sensitive to measurement scales. Such effects may also persist in real-world datasets due the measurement scale used.
To get rid of these effects, standardise all datasets before applying any of the causal discovery methods discussed in this paper. For the synthetic datasets, we also normalise each variable during the data generating process. This has been shown to create graphs that are not Var-sortable or R² -orttbbble, meaning the causal order cannot be determined using variances and correlation artefacts (Ormaniec et al., 2024). This means our results for the baselines don't 't match the results reported in the original papers for similar datasets (Montagna et al., 2023; Lachapelle et al., 2019; Rolland et al., 2022) (note that some of these works also define SHD of an anti-causal edge as where we define it as 2).
21

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
SHD 1
SID
F1 1
0.6
0.55
0.90
0.50
0.88
0.5
0.45
0.86
0.4
0.40
0.84
0.35
0.82
0.3
0.80
2
4
6
0
6
8
0
2
8
Number of Random Restarts
Number of Random Restarts
Number of Random Restarts
Figure 3: Effect of adding random restarts for the DGP-CDE. We record the metries (SHD, SID and F1) for the restart with the lowest loss out of growing set of restarts, starting with one restart and adding one restart time. We randomly permute the order in which the random restarts are added 50 times and plot the mean and standard deviation for each of the metrics. This shows that, on average, increasing the number of random restarts improves performance.
F. Data Generation Details
3 Variable Synthetic Data:
For each of the six distinct causal structures for three variables (up to permutations of the variables), we generated five datasets, each with 1000 samples. The data was generated by sampling a GP-CDE with a sum of linear kernel and one other kernel randomly selected out of Matérn12, Matérn32, Matérn52, squared exponential and rational quadratic. The kernel hyperparameters are sampled from o, Uniform(1, 100) and 0, ~ Gamma( (1.mma(1.5,, 1) while the other parameters are sampled as Φ; Uniform(0.01, 1) and w; (0 1). The data is then sampled from this GP-CDE.
20 & 50 Variable Synthetic Data: For both 20 and 50 nodes, create synthetic graphs using the ER (Erdos et al., 1960) and SF (Barabási & Albert, 1999) sampling schemes, with expected edges of one and four to simulate sparse and dense graphs. For each graph type, we generate five random graphs. Data is then generated for each node X; using
X; = f. (XPAç()) 9))
(31)
where (i) the parents of X, and 6 ~ (O, 1) is sampled from 6 ~ WN 1). NNN is randomly initialised neural network with two layers, 128 units and ReLU activation functions. For each graph, we sample 1000 data points. This data generating scheme ensures that the final data is generated from a complicated distribution which does not directly correspond to any of the models.
Syntren Dataset: Syntren is synthetic data generator that approximates real transcriptional regulatory networks (Van den Bulcke et al., 2006). Networks are selected from previously described transcriptional regulatory networks. Relationships between the genes are based on the Michaelis-Menten and Hill enzyme kinetic equations. The kinetic equations contain some biological noise and some lognormal noise is added on top. We use the data generated by Lachapelle et al. (2019). This contains 10 datasets of 500 samples, each with 20 nodes. The data was generated using the syntren data generator (Van den Bulcke et al., 2006), using the coli network with the default parameters except for the probability for complex 2-regulator interactions, which was set to (Lachapelle al., 2019).
G. Details about Baselines
SCORE NoGAM: For SCORE (Rolland et al., 2022) and NoGAM (Montagna et al., 2023) we use the implementation in the dodiscover package https: //ittuu com/ www/pymww//ccccccccc dodi scover. SCORE and NoGAM have very similar hyperparameters, and for both methods, we use the values of ridge regression suggested by Montagna et al., which they
22

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
tuned to minimise the generalisation error on the estimated residuals (Montagna et al., 2023). We also use the values suggested by the authors for both methods for all the hyperparameters for the Stein gradient, Hessian estimators and CAM pruning step except for the cutoff value for the CAM pruning step which we set to 0.01 is standard for CAM and follows Montagna et al. (Montagna et al., 2023). NoGAM has the additional hyperparameter of number of cross validation models, for which we use the preset value of five.
CAM: For CAM, we preliminary neighbourhood selection and DAG pruning method with cutoff value of 0.001 for the edge pruning as suggested by the authors (Bühlmann et al., 2014). We used the implementation in the https: / //github. terithu.com//a com/ wtw../ooo/okooonne wwwowaaan/rraaatt DAG repository.
For NOTEARS, we used the implementation provided in https::///gtt wtw.uuu..ooaa.iiai wttp:////.///aiiiaa.aa.a which uses binary search to find the adjacency matrix threshold that gives the largest possible DAG. This procedure has been found to give better performance than the default threshold of 0.3 in the original NOTEARS paper (Nazaret et al., 2024; Lopez et al., 2022).
NOTEARS:
SDCD: Nazaret et al. did not tune their hyperparameters for every dataset, instead using the same hyperparameter values, which were found to perform well empirically, for all their experiments (Nazaret et al., 2024, Appendix C.2). We use these same hyperparameter values and their implementation provided in ht https: http//ggthhb.///e// com/ www./occciaiiiat.c miizz/iadasddd...
We use the non-linear model provided by the authors in https: graph priors, we use the ER prior for the ER graphs and SF prior for the SF graphs, Syntren and the three variable baseline. This is because SF more closely matches the structure of the Syntren and the three variable datasets. For the other hyperparameters we use the values fine-tuned by the authors, including an observational noise of 0.1 and 3000 optimisation steps. The authors fine-tuned different values for the bandwidths of the kernel and slope of the linear schedule depending on the number of variables. We use these values, which are listed in (Lorch et al., 2021, Appendix E.3). DiBS typically returns a posterior over DAGs using Stein variational gradient descent. When a single particle is used, this becomes the MAP estimate of the graph. Therefore, as wish to compare particle for DiBS, and perform five random restarts, selecting the one with the highest marginal likelihood value.
DiBS:
/iiss. For
method, which returns MAP estimate, we use one
For the three variable dataset, we include uniform random sampling of DAGs. To compute this, for each ground truth graph, we enumerate all 25 possible graphs and evaluate their performance metrics. This exhaustive enumeration ensures that our results approximate those obtained in the limiting case where the number of sampled DAGs tends to infinity.
Random:
H. Code Availability & Computational Resources
Code and data to replicate the experiments in this paper can be found at https: //github. com/ Anish114 44/ ContinuousMMSEE wwwtutettttuuuuuuuuuueeee The experiments in this paper were run on 4100 and RTX 4090 GPUs.
I. Additional Experiments
I.1. 20 Variable ER4 Comparison with Additive Noise Model
Our model is more expressive than previous methods, many of which rely on the additive noise assumption. Here, we answer the following question: If the dataset is generated from an ANM model, how does our more expressive model perform?
To test this, we introduce the CGP, which is the same as the CGP-CDE model (Section 4), except that it is restricted to an additive noise model. This is done by simply removing the latent variable w in Equation (8). The likelihood then is a Gaussian, hence an ANM model. Removal of the latent variable also means that inference more accurate in this model than the CGP-CDS however the cost of a loss of flexibility. All the settings of the CGP model are exactly the same as in Appendix C.
Figure 4a shows the performance on five datasets generated from 20 variable ER4 graphs with Gaussian process functional relations. Thus, the datasets were generated from an additive noise model. It is worth noting that the CGP outperforms other methods that make the ANM assumption, such as CAM, SCORE, NOGAM, DiBS, and SDCD. Despite being more expressive than the CGP model, the CGP-CDE model performs on par with the CGP model. This shows that with the
23

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
SHD J for 20 Variable ER4 ANM
SID for 20 Variable ER4 ANM
F1 ? for 20 Variable ER4 ANM
120
350
0.6
100
300
80
0.4
250
60
0.2
40
NOTEARS
SDCD
CAM
NOGAM
SCORE
CGP CDE
DiBS
CGP
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CDE
CGP
NOTEARS
SDCD
CAM
NOGAM
SCORE
CGP CDE
DiBS
CGP
CGP CGP-CDV
(a)
SHD for 20 Variable ER4
SID for 20 Variable ER4
for 20 Variable ER4
100
350
0.8
80
300
0.6
60
250
0.4
40
200
20
0.2
NOTEARS
SDCD
CAM
NOGAM
SCORE
CGP CGP-CDE
DiBS
CGP
NOTEARS
SDCD
CAM
NOGAM
SCORE
CGP CDE
DiBS
CGP
NOTEARS
SDCD
CAM
NOGAM
SCORE
CGP CDE-CDE CDE
DiBS
CGP
(b)
Figure 4: 20 variable ER4 graphs, with data generated from (a) additive noise model (Gaussian process) and (b) non-additive noise model (neural network) described by Equation (31). Black dot is the average metric. Lower is better for SHD and SID, and higher is better for Fl. Despite the CGP-CD being more expressive than the CGP, it is still able to perform comparably for the ANM dataset, demonstrating the added flexibility of the CGP-CDE does not hinder the performance on datasets generated from restricted models.
added flexibility our model, we do not lose performance on datasets generated from restricted models. This is due to the Occam's 's razor effect of Bayesian model selection, which allows for expressing simpler models when the data requires (Rasmussen & Ghahramani, 2000).
For comparison, Figure 4b shows results for five 20 variable ER4 non-additive noise datasets generated using neural network described in Equation (31). Here, the CGP-CDE outperforms the CGP due to its more flexible functional assumptions. Also, the CGP outperforms the other ANM methods, showing the benefit of the Bayesian model selection approach even when restrictive functional assumptions are made.
I.2. 20 Variable with Different Numbers of Samples Results
To investigate the effect of the number of samples on the performance of the models, we performed experiments on five 20 variable ER4 graphs, with 100, 500 and 1000 random samples. These results are shown in Figure 5. The CGP-CDE generally outperforms the other methods for metrics and numbers of samples, except for with 100 samples, where CAM performs slightly better and NoGAM and SCORE perform comparatively, however, the CGP-CDE outperforms these methods on SID and F1.
24

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
SHD |. for 20 Variable ER4
SID for 20 Variable ER4
F1 for 20 Variable ER4
350
0.8
125
100
300
SHDD
75
250
0.4
50
200
25
0.2
100
500
1000
100
500
1000
100
500
1000
No. Sumples
No. Samples
No. Samples
Model
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP-CDE
Figure 5: 20 variable ER4 graphs with data generated for 5 graphs described by Equation (31). The axis is the number of samples per graph, the circular markers are the mean and the vertical lines show the standard error across graphs. Lower is better for SHD and SID, higher is better for F1.
1.3. 20 Variable Results
To show the CGP-CDE performs well with a smaller number of variables on a range of graph types and densities, we performed experiments on 20 variable graphs with the following graph types and densities: ER1 (Figure 6a), ER4 (Figure 4b), SF1 (Figure 6b), and SF4 (Figure 6c). These results show the CGP-CDE consistently performs well, although it outperformed some metrics for the ER1 and SF1 graphs, partieularly by NOGAM, which seems to perform better sparser graphs, perhaps because it is an ordering based method. However, the CGP-CDE is more consistent in performance than NOGAM across all our experiments with varying number of variables and graph types.
1.4. Full 50 Variable Results
To demonstrate the performance of the CGP-CD variable graphs with the following graph types and densities: ER1, ER4, SF1, and SF4. The results of these experiments can be seen in Figure 7. These results show the CGP-CDE consistently performs well. For the ER graphs, the CGP-CDE outperforms all the baselines on SHD, SID and F1. For the SF graphs, where a few variables cause the rest of the variables, some of the baselines outperform the CGP-CDE on SID. We hypothesise this is because the structure of the SF graphs benefits methods that use topological search, such as NOGAM and SCORE. However, the CGP-CD greatly outperforms these methods on SHD and F1 for both the SF graph densities.
range of graph types and densities, performed experiments on 50
25

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
SHD for Variable FRI
STD 20 Variahle
F1 for Variable ERI
0.8
100
30
60
0.4
40
10
0.2
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDP-CDE
CDE
NOTEARS
SDCD
CAM
NOGAM SCORE
DIBS
CGP CDE
(a)
SHD for 20 Variable SF1
STD for 20 Variahle SFI
Fl for 20 Variable SFI
40
200
0.6
30
150
100
10
50
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DIBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
(b)
SHD for 20 Variable SF4
SID for 20 Variable SF4
F1 for 20 Variable S14
0.8
350
80
300
0.6
60
250
0.4
40
200
0.2
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DIBS
CGP CDE
(0)
Figure 6: (a) Erdos-Renyi (ER) graphs with expected edge per variable. (b) Scale-Free (SF) graphs with expected edge per variable. (d) SF graphs with four expected edges per variable. For results for Erdos-Renyi (ER) graphs with four expected edges per variable, see Figure 4b. Black dot is the average metric. Lower is better for SHD and SID, higher is better for FI.
26

---

Continuous Bayesian Model Selection for Multivariate Causal Discovery
SHD J for 50 Variable BRI
SID for 50 Variable ERI
Fl 1 for 50 Variable ERI
0.6
400
200
300
0.4
100
200
0.2
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP - CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CDE
NOTEARS
SDCD
CAM
NOGAM SCORE
DiBS
CGP CDE
CGP CGP-
@
SHD J for 50 Variable ER4
SID for 50 Variable ER4
Fl 1 for 50 Variable ER4
800
2250
600
2000
400
1750
0.4
1500
200
1250
0.2
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DIBS
CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DIBS
CGP CDE
CGP
(b)
SHD for Variable SF1
SID J
50 Variable SF1
Fl for Variable SFI
600
1200
0.6
1000
400
800
0.4
200
600
0.2
400
NOTEARS
SDCD
CAM
NOGAM
SCORE
DIBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM SCORE
DIBS
CDE
NOTEARS
SDCD
CAM
NOGAM SCORE
DiBS
CGP CDE
CGP
(c)
SHD for 50 Variable SF4
SID | for 50 Variable SF4
Fl for 50 Variable SF4
600
2000
0.6
400
1500
0.4
200
1000
0.2
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DiBS
CGP CDE
NOTEARS
SDCD
CAM
NOGAM
SCORE
DIBS
CGP CGP-CDE
CDE
(d)
Figure 7: (a) Brdos-Remi (ER) graphs with one expected edge per variable. (b) ER graphs with four expected edges per
variable. (c) Scale-Free (SF) graphs with one expected edge per variable. (d) SF graphs with four expected edges per variable. Black dot is the average metric. Lower is better for SHD and SID, higher is better for F1.
27


Reviews:
[{"review_id": "sVJbK84tCA", "reviewer": "Reviewer_apg4", "summary": "The paper proposes a new method called CGP-CDE for (Bayesian) causal model discovery that allows for less restrictive model assumptions and can be applied to higher dimensional systems as well. It is based on a GP approach to obtain a nonparametric conditional density estimator for each node given its parents in the causal DAG, ultimately capturing the likelihood of a graph. The search for the best fitting model is turned into a continuous optimization problem by adding a familiar acyclicity constraint with penalty on the weight matrix representing the graph. The method is evaluated on synthetic/realistic data and found to compare equal or favourably to other alternatives.\n\n## update after rebuttal\nI thank the authors for their reply, and, having read the other reviews & rebuttals as well, I am happy to leave my score at '4: accept'.\n\nOne final remark: I understand the temptation to go for maximal informativeness, but the authors will know that typically the score differences within a MEC are much smaller than between MECs, and therefore the extra information from orienting the full graph tends to be much less reliable, in turn making the entire output less trustworthy = less useful in practice. Given that any directed graph has a unique mapping to a MEC one does not need to 'remove' the preference within a MEC. Instead, we can easily show both outputs, but with a much higher reliability for anything implied by the MEC representation, giving practitioners an intuitive way to distinguish between more and less reliable conclusions in the output.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "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\"**.\n\n> The paper initially suggests that it will solve the problem of restrictive / unrealistic model assumptions encountered when tackling real world data\n\nPrevious 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). \nWe 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). \n\n> Yes it ensures identifiability, but you are essentially finding back what you put in and assume/hope that it matches reality.\n> ...using an in practice hard to justify/assess prior\n\nOur 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).\n\nThe specific priors over functional mechanisms are also important.\nOur 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).\n\nWe 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.\n\n[1] Hjort et al., eds. Bayesian nonparametrics. Vol. 28. Cambridge University Press, 2010.\n\n> 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.\n\nThis 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.\n\n> 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...\n\nThe papers you mention are important but only consider simple linear models.\nHowever, 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].\n\n[1] Dhir et al., \"Bivariate Causal Discovery using Bayesian Model Selection.\" ICML, 2024.\n\n[2] Friedman et al., \"Gaussian process networks.\" Proceedings of the Sixteenth conference on Uncertainty in artificial intelligence. 2000.\n\n[3] Stegle et al. \"Probabilistic latent variable models for distinguishing between cause and effect.\" Advances in neural information processing systems 23 (2010).\n\n> I would have preferred a modest but more robust aim for e.g. the MEC.\n\nThis 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.", "claims_and_evidence": "The paper initially suggests that it will solve the problem of restrictive / unrealistic model assumptions encountered when tackling real world data, but this is of course nonsense. It still starts from the causal sufficiency and acyclicity assumptions (which is not how the world works), and relies on a rather arbitrary prior to go from Markov equivalence class (MEC) to unique model. Yes it ensures identifiability, but you are essentially finding back what you put in and assume/hope that it matches reality.\nHowever, I do really like the GP approach to modelling the likelihood, which indeed does provide a significant improvement over existing parametric approaches to Bayesian/likelihood based inference.", "methods_and_evaluation_criteria": "Experimental evaluation is limited but ok, and demonstrates the potential of the core contribution.", "theoretical_claims": "No explicit theoretical claims in the main paper, but overall approach is sound.", "experimental_designs_or_analyses": "As mentioned, experimental evaluation is limited but ok. 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. But I will not hold that against the authors :)", "relation_to_broader_scientific_literature": "Relevant (recent) work is discussed.", "essential_references_not_discussed": "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 subsequent work by Heckerman showed how to combat this (rather undesirable) behaviour by introducing a score that ensures MEC-equivalence, but the current paper now seems to suggest that they came up with the idea of selecting models *within* a MEC by using a non-equivalent score …", "other_strengths_and_weaknesses": "As stated: paper suggests more than it does, and some aspects (like aiming for unique identifiability from using an in practice hard to justify/assess prior) actually weakens the output, and I would have preferred§ a modest but more robust aim for e.g. the MEC.\nAlso starting from the causal sufficiency assumption (and to some degree also the acyclicity assumption) is `unforgivable’ when the aim is to make a causal discovery method more suitable to real-world applications. Yes it makes everything easier, but that is not a good enough reason to do / keep doing it.\n\nOn the plus side: paper is clear and well written, and the GP approach at the heart of the method is a significant and promising contribution, and for that reason I will recommend accept.", "other_comments_or_suggestions": "No", "questions_for_authors": "No", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ffpBXrY6Ob", "reviewer": "Reviewer_Ry2n", "summary": "Recent work shows that in the bivariate case Bayesian model selection can be used for structure identification under more flexible assumptions at the cost of a small probability of error. This paper extends the previous result to the multivariate case. The authors empirically validate the method by comparing to multiple baselines under extensive experimental settings.\n\n## update after rebuttal\n\nThank you for the author's response. After reading all the review comments I decide to keep my rating unchanged.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "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.\n\n> do we need to know the correct prior in advance?\n\nIt 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]).\n\n[1] Hjort et al., eds. Bayesian nonparametrics. Vol. 28. Cambridge University Press, 2010.\n\n> what will happen under misspecification of prior?\n\nWe 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]. \n\nWe 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**. \n\n[1] Dhir et al., \"Bivariate Causal Discovery using Bayesian Model Selection.\" ICML, 2024.\n\n> 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.\n\nDhir 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?\nWe 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.\nWe then show that this approach outperforms competing methods in the multivariate case.", "claims_and_evidence": "Yes.", "methods_and_evaluation_criteria": "It makes sense to me that for two distribution equivalent graphs we can rely on independent priors to distinguish them. Yet, it is unclear to me that do we need to know the correct prior in advance? If yes, what will happen under misspecification of prior?", "theoretical_claims": "No.", "experimental_designs_or_analyses": "Yes. The experimental designs look sensible to me. The setting is extensive including both synthetic setting with multiple variables and real-life experiments.", "relation_to_broader_scientific_literature": "Key contributions of the paper compared to prior work are properly discussed.", "essential_references_not_discussed": "Related works are properly discussed.", "other_strengths_and_weaknesses": "Strengths \n\n1. The proposed method and theory of using bayesian model selection for multivariate structure learning look novel to me.\n\n2. The experimental setting is extensive including both synthetic and real-life data.\n\n3. It is interesting to see methods that can distinguish between graphs that are distribution-equivalent.\n\nWeaknesses\n\n1.  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.", "other_comments_or_suggestions": "N.A.", "questions_for_authors": "Please refer to my question in Methods And Evaluation Criteria part.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "kGYmeQsJqZ", "reviewer": "Reviewer_2axY", "summary": "This paper presents a multivariate causal discovery approach based on Bayesian model selection. It builds on the work of Dhir et al. (2024), who proposed to use Bayesian model selection to identify causal direction in the bivariate case. The Bayesian model selection framework allows for a trade-off between a model's goodness of fit and complexity. The original work showed that Bayesian Model Selection can discriminate causal directions even when Maximum Likelihood (ML) methods fail due to distribution equivalence, thanks to the marginal likelihood, which, unlike the likelihood, is not symmetric (they consider the whole function space instead of the best-fitting function).\nThey use the Causal Gaussian Process Conditional Density Estimator (CGP-CDE) to impose no restrictions on the model (e.g. linear model, additive noise...), which is a significant contribution of this work. This makes it possible to relax the assumptions, one of the most critical challenges in causal discovery. On the other hand, their method does not come with a strict identifiability guarantee due to potential overlaps in the posterior of different causal models. They propose a way to estimate the error probability using a sampling strategy.\nThe causal identifiability proof is based on the Independent Causal Mechanism (ICM) assumption, which the authors use to claim that well-chosen Bayesian priors for the distribution of a cause and the distribution of an effect given the cause should be independent.\nIn practice, to avoid comparing all possible DAGs and their posteriors, the model is identified by continuously optimizing the hyperparameters of the Gaussian process priors using a variational autoencoder (VAE) with a Bayesian model selection-based loss and acyclicity penalty. They also define priors on the resulting directed acyclic graphs (DAGs) to favor sparse models. As with most existing methods, their approach only works if all confounders are observed (causal sufficiency assumption).", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "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\"**. \n\n> whether the found model respects the Markov condition \n\n**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). \n\nWe 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.\n\n[1] Dawid et al., \"Beware of the DAG!.\" Causality: objectives and assessment. PMLR, 2010.\n\n[2] Janzing et al., \"Causal inference using the algorithmic Markov condition.\" IEEE Transactions on Information Theory (2010).\n\n> It is unclear whether the continuous optimization method finds well-chosen priors\n\nOur continuous optimization does not find the prior, but the *posterior* given a prior (L255 RHS). \n\n> literature on the information-theoretic view of causality (Janzing, Schölkopf, 2010) \n\nThe 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.\n\n[1] Dhir et al., \"Bivariate Causal Discovery using Bayesian Model Selection.\" ICML, 2024.\n\n> There is no comparison with other score-based approaches\n\nWe 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. \n\n> The causal model is not strictly identifiable.\n\nThis 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.\n\n> How do you interpret the poor performance of CGP-CDE in the 3-variable experiment?\n\nWe 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.\nThe 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).\n\n> Could another initialization reduce the number of iterations needed?\n\nWe 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.\n\n> In the 3-variable experiment... How do you explain this?\n\nThe 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.\n\nNote 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.\n\nThe reviewer refers to questions 4 and 5 but the questions only go up to 3.", "claims_and_evidence": "Not all claims are theoretically substantiated. The authors claim to find a causal model. Their continuous optimization method yields a direct acyclic graph, but it is not made clear whether the found model respects the Markov condition (each node is conditionally independent of its non-descendants, given its parents), which is required for a DAG to be causal. It is unclear whether the continuous optimization method finds well-chosen priors, which is required for identifiability.", "methods_and_evaluation_criteria": "The synthetic data experiment is performed with different generation settings (neural network as function, noise additive or not...). Syntren is a simulated data set. There is no experiment on real data. SHD and SID are commonly used metrics for causal discovery.", "theoretical_claims": "I read the proofs in the appendix and found no errors.", "experimental_designs_or_analyses": "The experimental design includes several settings. Results on synthetic data generated with a neural network show that the proposed method largely outperforms competitors, which was expected since their assumptions (e.g. additive noise) are violated. However, the author(s) also provide results for data generated with additive noise and their performance is still competitive with other methods, showing that the relaxation of the assumptions does not come at the expense of overall performance.", "relation_to_broader_scientific_literature": "This paper contributes to the field of causal discovery and, more specifically, to the score-based and continuous optimisation-based lines of research. It extends the work of Dhir et al. (2024), who proposed using Bayesian model selection for causal edge orientation. This paper's novelty is its adaptation and proof for the multivariate case, the use of CGP-CDE for flexible modeling, and the continuous optimization of model parameters to find a DAG.", "essential_references_not_discussed": "The authors do not discuss the literature on the information-theoretic view of causality (Janzing, Schölkopf, 2010) in the related works although it is conceptually very close -- instantiating Occam's razor to distinguish between models. There is no comparison with other score-based approaches that also aim to balance complexity and goodness of fit, such as GES (Chickering, Maxwell, 2002), which uses the Bayesian Information Criterion (BIC), and GLOBE (Mian, 2021), which uses a two-part Minimum Description Length (MDL) score.\n\nReferences:\n- Chickering, David Maxwell. \"Optimal structure identification with greedy search.\" Journal of machine learning research 3.Nov (2002): 507-554.\n- Mian, Osman A., Alexander Marx, and Jilles Vreeken. \"Discovering fully oriented causal networks.\" Proceedings of the AAAI Conference on Artificial Intelligence. Vol. 35. No. 10. 2021.\n- Janzing, Dominik, and Bernhard Schölkopf. \"Causal inference using the algorithmic Markov condition.\" IEEE Transactions on Information Theory 56.10 (2010): 5168-5194.", "other_strengths_and_weaknesses": "Strengths:\n- The paper is well written and well motivated.\n- The use of a flexible model and VAE allows for fewer assumptions than usual and therefore makes the model applicable to a wider range of applications.\n- The theory is well covered and convincing.\n- The extensive experiments show that the proposed methods significantly outperform the evaluated competitors.\n\nWeaknesses:\n- The causal model is not strictly identifiable.\n- Continuous optimisation may not be efficient (see questions 1, 2 and 3).\n- It is unclear whether the continuous optimisation approach maintains the theoretical guarantees (see questions 4 and 5).", "other_comments_or_suggestions": "NOTEARS is a structure learning approach, not a causal discovery algorithm. That is, it finds a DAG, but there are no guarantees this a causal one. This should be mentioned. I have the same concern about the proposed method.\n\n\n### update after rebuttal ###\n\nI thank the authors for their answers, but I'm afraid they did not take away my concerns regarding the identifiability of the approach and therewith whether this is a method that is guaranteed to return causal structures.", "questions_for_authors": "1) How do you interpret the poor performance of CGP-CDE in the 3-variable experiment?\n\n2) In the continuous optimization, the warm-up phase seems computationally expensive (about 25,000 iterations for 1,000 samples?). Is this a typical number of iterations for an initialization phase? Could another initialization reduce the number of iterations needed? What about the number of iterations for the cooling phase (also 25,000 iterations for 1,000 samples)?\n\n3) In the 3-variable experiment, where the data are generated using a Bayesian causal model, DGP-CDE acts as a sanity check, as it considers all models and thus includes the true one used to generate the data.  However, the poor performance of the continuous optimization is worrying given the small search space (25 graphs) and the fact that the data was generated using a Bayesian causal model (unlike the opponent whose assumptions were violated). How do you explain this?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bQ7CD34PdV", "reviewer": "Reviewer_7Mfs", "summary": "This paper studies structure learning for observational data using Bayeisan model selection. It falls into the category of score based learning and uses model evidence as score to select DAG. It shows the existing work on bivariate case [Dhir et al., 2024] can be extended to multivariate case, and applies a flexible model called Causal Gaussian Process Conditional Density Estimator (CGP-CDE) with continuous reparametrization and variational inference to learn the DAG. Experiments are conducted to show the effectiveness.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "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. \n\n>However, the proof of this theorem only says it follows directly from [Dhir et al., 2024] without any details.\n\nOur 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.\n\n> 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.\n\nWe 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.\n\n> I wonder for two models that are not Bayesian equivalent, why can they have positive probrobability of error\n\nThe 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**.\nHence, 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).\n\n[1] Dhir et al., \"Bivariate Causal Discovery using Bayesian Model Selection.\" ICML, 2024.\n\n[2] Hoyer et al. \"Nonlinear causal discovery with additive noise models.\" Advances in neural information processing systems 21 (2008).\n\n> 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?\n\nThe 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.\n\n>Or why do we have identifiability and distinguishabllity and what are their differences?\n\nIdentifiability 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.", "claims_and_evidence": "The claims on improvement are supported via experiments.", "methods_and_evaluation_criteria": "The evaluation criteria via SHD, SID, F1 are a commonly used critera.", "theoretical_claims": "One of the proofs on the main results --  Theorem B.6 -- concerns me. One of the main contribution of this paper is to extend the existing identifiability result to multivariate case. However, the proof of this theorem only says it follows directly from [Dhir et al., 2024] without any details.", "experimental_designs_or_analyses": "The experimental design makes sense.", "relation_to_broader_scientific_literature": "Causal discovery is widely used in scientific research for exploratory analysis. Developing identifiability and scalable method is important for broader scientific literature.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "- The paper derives from the ICM assumption for identifiability and allows for much flexible class of identifiable Bayesian network\n- The proposed method leverages the powerful tools for modelling and learning nonparametric DAGs in a scalable manner.\n- The experiments show competitive performance with many benchmarks.", "other_comments_or_suggestions": "- 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.", "questions_for_authors": "- The definition of Bayesian equivalence (identifiability and distinguishabllity) involves data $X$, also relates to the definition of probrobability of error. I wonder for two models that are not Bayesian equivalent, why can they have positive probrobability of error. 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? Or why do we have identifiability and distinguishabllity and what are their differences?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper primarily extends Dhir et al. (2024) from the bivariate to the multivariate setting, using the old idea (as one reviewer points out, this dates back at least to Cooper & Herskovitz, 1992) of letting priors select causal between models that would otherwise be equivalent. The key contribution is to propose a new differentiable NOTEARS-style model under a nonlinear Gaussian process model for the conditional densities (CGP-CDE) that is more flexible. The authors show that their method outperforms recent methods based on differentiable causal discovery.

An important part of this line of work is to make the methods and evaluation open-source. During the discussion, the authors clarified that they "plan to make the Github repo containing the code public upon publication", which I expect the authors will do, with prominent links to the repo in the camera ready.
Decision: Accept (poster)