Causal Discovery via Bayesian Optimization
Authors: Bao Duong
Sunil Gupta
Thin Nguyen
Venue: ICLR 2025 Poster (2025)
Keywords: causal discovery
causal structure learning
bayesian optimization
Primary area: causal reasoning
TL;DR: We propose a novel score-based causal discovery method based on Bayesian optimization, allowing us to reach higher-scoring DAGs with fewer trials and less time.

Abstract: Existing score-based methods for directed acyclic graph (DAG) learning from observational data struggle to recover the causal graph accurately and sample-efficiently. To overcome this, in this study, we propose DrBO (DAG recovery via Bayesian Optimization)—a novel DAG learning framework leveraging Bayesian optimization (BO) to find high-scoring DAGs. We show that, by sophisticatedly choosing the promising DAGs to explore, we can find higher-scoring ones much more efficiently. To address the scalability issues of conventional BO in DAG learning, we replace Gaussian Processes commonly employed in BO with dropout neural networks, trained in a continual manner, which allows for (i) flexibly modeling the DAG scores without overfitting, (ii) incorporation of uncertainty into the estimated scores, and (iii) scaling with the number of evaluations. As a result, DrBO is computationally efficient and can find the accurate DAG in fewer trials and less time than existing state-of-the-art methods. This is demonstrated through an extensive set of empirical evaluations on many challenging settings with both synthetic and real data. Our implementation is available at https://github.com/baosws/DrBO.

Published as a conference paper at ICLR 2025
CAUSAL DISCOVERY VIA BAYESIAN OPTIMIZATION
Bao Duong, Sunil Gupta, and Thin Nguyen Applied Artificial Intelligence Institute (A²I²), Deakin University, Geelong, Australia [b. duong, sunil. gupta, thin. nguyen @aaryeeneddekkii @deakiin edu.au
ABSTRACT
Existing score-based methods for directed acyclic graph (DAG) learning from observational data struggle to recover the causal graph accurately and sample- efficiently. To overcome this, in this study, we propose DrBO (DAG recovery via Bayesian Optimization)- a novel DAG learning framework leveraging Bayesian optimization (BO) to find high-scoring DAGs. We show that, by sophisticat- edly choosing the promising DAGs to explore, we can find higher-scoring ones much more efficiently. To address the scalability issues of conventional BO in DAG learning, we replace Gaussian Processes commonly employed in BO with dropout neural networks, trained in continual manner, which allows for (i) flexi- bly modeling the DAG scores without overfitting. (ii) incorporation of uncertainty into the estimated scores, and (iii) scaling with the number of evaluations. As a result, DrBO is computationally efficient and can find the accurate DAG fewer trials and less time than existing state-of-the-ar methods. This is demon- strated through an extensive set of empirical evaluations on many challenging settings with both synthetic and real data. Our implementation is available at https:/ / hwg.tb...h...gggggg com baosws Saaosws/DrBO. /DrBO.
INTRODUCTION
Learning directed acyclic graphs (DAGs) encoding the underlying causal relationships, also known as causal discovery, provides invaluable insights about interventional outcomes and counterfactuals, and thus significant research effort has been dedicated on this frontier. Our study focuses on the score-based framework-a major class of causal discovery methods which casts the DAG learn- ing problem optimization problem, maximizing over the space of DAGs predefined score function measuring how a DAG G fits the observed data D;
G arg max S (D,9). GEDAGs
(1)
There are several challenges associated with this formulation. First, this optimization problem is NP-hard in general (Chickering, 1996), due to the combinatorial search domain that scales super- exponentially with the graph size (Robinson, 1977) and the acyclicity condition that is nontrivial to maintain. The second challenge is computational cost, as the score evaluation can be expensive for complex models (Zhu et al., 2020; Wang et al., 2021), and hence, methods requiring too many trials will ineur heavy computational expense.
Bayesian Optimization (BO) has emerged as an effective approach for expensive black-box opti- mization thanks to its sample efficiency. Its applicability covers many domains due the pervasive- ness of optimization tasks in virtually every field. The central idea is that past evaluation data can reveal potential candidates better solutions using fewer evaluations. However, while BO has been employed for active causal discovery (Toth et al., 2022; Zhang al., 2024) to suggest cost-effective intervention strategies to discover causal graphs from active experiments, its power is yet to be harnessed for the problem of observational causal discovery, where no active intervention available.
evaluate next, and effectively exploiting them allows
arrive at
We have identified two potential difficulties preventing BO to optimize Eq. (1). The first is the limited scalability of BO in general, which usually restricted to only a few hundred dimensions and thousands of evaluations (Wang et al., 2023), while DAG learning often involves far more trials (Zhu et al., 2020; Wang et al., 2021; Duong et al., 2025). The second challenge is how to efficiently optimize the "acquisition function" measuring the potential of DAG solutions in BO. This itself is scsre DAG learning problem and is expected to be cheaper than the original problem as is repeated many times in the BO pipeline, thus requiring to be very efficient to be practical.

---

Published as a conference paper at ICLR 2025
Present study. In this study, we tackle these obstacles and bring forth the benefits of BO into causal discovery with the introduction of DrBO (DAG recovery via Bayesian Optimization) a novel causal discovery algorithm that leverages BO to find the highest- DAG. More partic- ularly, we solve the aforementioned challenges by several key design choices. (i) Inspired by Yu et al. (2021); Massidda et al. (2024); Duong et al. (2025), we devise a low-rank DAG representa- tion that relaxes the constrained optimization in Eq. (1) to an unconstrained optimization problem with a dimensionality growing linearly with the number of nodes (Sec. 4.1), enabling the use of BO with an amenable dimensionality. (ii) We replace the Gaussian processes (GPs) in conventional BO approaches that scale cubically w.r.t the number of evaluations with dropout neural networks (Srivastava et al., 2014; Guo et al., 2021), which offer expressive uncertainty-awar modeling ca- pabilities, as well as faster acquisition function calculation and optimization (Sec. 4.3). (iii) Our surrogate model learns the DAG score indirectly via node-wise local scores, allowing us to predict the DAG scores more accurately with enhanced training information, compared with learning a di- rect map from a DAG to its score and not exploiting the local scores (Sec. 4.4). (iv) Instead of being retrained every step with all data, our neural networks are trained continually, enabling our method to scale better with the number of trials (Sec. 4.5). Compared existing optimization strategies in causal discovery, like gradient-based -baed methods which do not attempt to exploit past exploration data to prioritize visiting the most promising DAGs, our method makes more informed decisions about which DAGs to investigate next, leading to fewer unnecessary trials to reach high-scoring DAGs.
Contributions. The main contributions of our study are summarized as follows:
1. To facilitate sample-efficiet causal discovery. we propose DrBO-a causal discovery method employing BO to optimize for the DAG score. DrBO is specifically designed for causal discovery, aiming to be not only accurate, but also computationally manageable. To our knowledge, this is the first score-based causal discovery method based on BO for purely observational data.
2.
We demonstrate the effectiveness of our method on a comprehensive set of experiments, showing that DrBO can consistently surpass sute-of-the-aa baselines on various condi- tions, with fewer evaluations and less time, signifying the sample efficiency of BO in our design. In addition, extensive ablation studies verify the significance of our design choices.
2 RELATED WORK
Causal discovery methods can be broadly categorized into two major classes, namely constraint- based and score-based methods. The former category (Spirtes et al., 2000; Colombo et al., 2012) involves performing a series of hypothesis tests to recover the undirected skeleton of the causal graph, before orienting the edges using graphical rules. On the other hand, score-based causal dis- covery recasts the problem as combinatorial optimization task. Classical methods (Chickering, 2002; Ramsey et al., 2017; Ramsey, 2015) greedily traverse the DAG space by adding and remov- ing edges one-at-a-time to maintain acyclicity. Recent advances include relaxing the combinatorial optimization to continuous optimization problem (Zheng et al., 2018; Yu et al., 2019; Zheng et al., 2020: Yu et al., 2021; Zhang et al., 2022; Bello et al., 2022; Annadani et al., 2023). In addition, meth- ods based on Reinforcement learning (RL) (Zhu et al., 2020; Wang et al., 2021; Yang et al., 2023a; Duong al., 2025) have recently emerged competitive search strategies. We also acknowledge interventional causal discovery studies (Hauser & Bühlmann, 2012; Brouillard et al., 2020; Lippe al., 2022) where interventional data is exploited to help identify the causal DAG, however, here we focus on the more challenging setting where no intervention is available.
Furthermore, Bayesian causal discovery studies are an intriguing direction (Deleu et al., 2022: Tran et al., 2023; Annadani et al., 2023), where the aim is to infer the posterior distribution over causal graphs given observed data, in order to quantify uncertainty DAG estimates. Meanwhile, our use of BO involves quantifying uncertainty in DAG scores, so despite sharing the term "Bayesian", " these studies are distant from us. Moreover, while Bayesian optimization has been employed for causal discovery in the active setting (Toth et al., 2022; Zhang et al., 2024), these methods utilize BO to suggest optimal interventions to quickly recover the causal DAG, necessitating the ability to perform active experiments. Meanwhile, our study utilizes BO to optimize for the score function that can be calculated purely from observational data. In addition, causal Bayesian optimization (Aglietti et al., 2020: 2021), which deals with optimizing variable of interest that is part of causal system with known DAG via a series of interventions suggested by BO, is also a different problem from ours, which focuses on finding the unknown DAG instead.
2

---

Published as a conference paper at ICLR 2025
3 BACKGROUND
Notations. In this paper, unless specifically indicated, normal lowercase letters indicate scalars (e.g., I, y) or functions (e.g.. f. 9). while bold lowercase letters represent vectors (e.g., x. y), and bold uppercase letters denote matrices (e.g.. X). Mcanwhile, subscripts and bracketed superscripts index dimensions and samples, respectively, e.g. xx? denotes the i-th dimension of the j-th sample.
3.1 BAYESIAN OPTIMIZATION
We provide here only the details necessary to understand our contributions. For a more compre- hensive review of BO, see Wang al. (2023), for example. Consider maximization of a function f that is expensive to evaluate: x* = arg maXxex' f (x). BO is a class of sequential optimization methods, which iteratively (i) proposes potential candidate(s) evaluate based on an "acquisition function", then (ii) evaluates said candidate(s), and (iii) updates its statistical model with the newly acquired observations. More specifically, BO defines a probabilistic "surrogate" model over the distribution f given observed data y, i.e., (J y),, to define the acquisition function.
Surrogate model. To model P y) Gaussian process (GP) the standard in BO since offers closed-form solution (Rasmussen, 2003). Assuming we have observed a dataset of n evalu- ations X = [x(").... x.f... x(")]" and y = .y(").... y(m) where x ~ X = IR and == f (x), GP assumes that y follows multivariate Gaussian distribution governed by mean function μ X R and positive-defintte covariance function k: X x X' + IR: y ~ N (A (x). Kx,x) where Kx.x == [ss (x(*), x()) i.j=1.....n is the n x n covariance matrix. Using Bayes' theorem, the posterior of the function value at a new location x is given analytically as: P (y x, x, - N (n (x), (x)) where μ (x) := xx,XxxXX and o² (x) == K (x,x) kx Kx xxxxxx x GPs scale poorly with n due to the need to invert Kx,x, so alternative statistical models like random forest (Hutter et al., 2011), Bayesian linear regression (Snoek al., 2015), and Bayesian neural network (Springenberg et al., 2016) have been employed as more scalable surrogate models.
Acquisition functions (AFs) BO judge how promising an arbitrary candidate x is based on the posterior inferred by the surrogate model, to make a more informed candidate suggestion. As an example, upper confidence bound (UCB) is a common choice, designed to minimize regret in the multi-armed bandit literature (Srinivas et al., 2010) and given by AF (x) := ft (x) + 3a (x), where > hyperparameter controlling the explotaation-eploratiio trade-off. Another promising acquisition is Thompson sampling (TS, Thompson, 1933), which is stochastic function that uses a random draw from the posterior as the potential indicator: AF (x) ~ (1 x, x)).
3.2 STRUCTURAL CAUSAL MODEL
JA- denote an i..dd dataset of n samples from P (x). The structural causal model (SCM, Pearl, 2000; 2009) among said variables can be described by (i) a DAG G = (V,E) where each node E V {1, (1,...,4 d} corre- sponds to a random variable xy. and each edge (j + i) E E C Vx Vimplies that æj is a direct cause of 24, (ii) set of functions {f dictating the causal mechanisms, and (iii) noise distribution
Let E IRd be the random vector capturing the system of interest, and D {x((a) )
(e). Together, these components define generative process 2; := (xpa Vi = 1, I....,d. d,
where pa. [j e V (j+ e E} the set of direct causes (a.k.a. structural parents) of node in G, and E ~ P (e) is the noise vector. Then, the observational causal discovery problem is con- cerned about recovering the DAG G from the observational dataset D. In addition, following (Zhu et al., 2020; Wang et al., 2021; Yang et al., 2023a;b; Duong al., 2025), we also assume: (i) causal sufficiency: there is no unobserved confounders among the variables; (ii) causal minimality: there no function fa that is constant to any of its argument (Peters et al., 2014); and (iii) identifiable causal models: this means G is the unique causal graph that can induce P (x), and thus it is possible to be recovered. For instance, while general (Spirtes et al., 2000), examples for identifiable causal models include linear-Gann Gaussian models with the equal-variance assumption (Peters al., 2014) and nonlinear additive noise models (ANMs) in general (Hoyer et al., 2008). Our experiments will adopt these identifiable models.
models are known to be unidentifiable
3.3 SCORE- BASED CAUSAL DISCOVERY
A critical component of score-based causal discovery is the proper specification of a scoring func- tion. With the proper scoring function, the optimization of the score is equivalent reaching to
3

---

Published as a conference paper at ICLR 2025
ground truth DAG. Consistent scoring functions (Chickering, 2002) are known to satisfy this quirement. In the main text, we demonstrate our method with the Bayesian Information Criterion (BIC, Schwarz, 1978), which is consistent sidered in numerous existing studies (Chickering, 2002; Zhu et al., 2020; Wang et al., 2021; Yang al., 2023a;b; Duong et al., 2025).
as shown by Haughton (1988) and is widely con-
More formally, let 0 = (ii P(e)} be the parameters of an SCM, then the BIC score is given
by SpIc (D,9) = 21p (D10.)) ||] where := arg max (D) 0, g) is the maximum- likelihood estimator of the causal model parameters. n is the sample size of D, and [g] denotes the number of edges G. The generality of BIC allows for its adaptation numerous causal models. In the main paper, we showcase our method with the BIC defined for the popular additive noise model (ANM, Hoyer et al., 2008): x; == Si (Xpa.) + & where e; N (0. of) In the general form, where the noise variances can be non-equal, the BIC for ANM is specified as follows:
SVICNV (D,9) = n Σ In MSE, (paç) In
where NV stands for "non-equal variance" and MSE, (pa n Σラ= f
is
pa.
the mean squared error after regressing 2; on xpa.. In addition, if we further assume that the noise variables have equal variances (Bühlmann et al., 2014) then the BIC yields:
El MSE; (pa.) d
SVICEV (D,9) == -nd In
[g] In n.
(3)
Eqs. (2) and (3) are common in prior studies (Zhu et al., 2020; Wang et al., 2021: Yang et al., 2023a;b; Duong et al., 2025), and we also provide their derivation in Appendix B.
3.4 PARAMETRIZED DAG GENERATION
For effective acquisition function optimization, we find it crucial for our method to be able to quickly generate candidate DAGs within specific regions determined by low-dimensional search space. This would help narrow down the regions of interest and improve the quality of suggested candi- dates. A potential approach towards this end is autoregressive DAG generation techniques (Wang et al., 2021; Deleu al., 2022; Yang et al., 2023a; Deleu et al., 2024), which break down the DAG generation into multiple sequential steps. However, each step of the generation process com- putationally involved due to their autoregressive nature (Duong et al., 2025), accumulating into significant DAG generation cost. Recently, constraint-free DAG representations with quadratic complexity are proposed in (Yu et al., 2021: Massidda al., 2024; Duong al., 2025), which in- troduce different maps from an unconstrained real-valued representation to the space of DAGs. For example, the Vec2DAG operator in (Duong et al., 2025) takes as input a continuous "node potential" vector Rd and strictly upper- trianulaa "edge potential" matrix E E dd to deterministically
create DAG: Vec2DAG (p, E) := H (grad @ H (E (T+ where H (x) 1, if > 0,
0, otherwise
is the entry-wise Heaviside step function, @ is the Hadamard product, and grad (p): = Pj ps is
the gradient flow operator (Lim, 2020).
0
2
4
Example. Consider system of 3 nodes with potentials p [-1,3,2] 2] and E
0
0
7
0
0
0
-
(
0
4
3
0
-
Then, Vec2DAG (p, E) = H
-4
0
-
H
2
0 7
0
0
3
-
This is the adjacency of the DAG 1 determined by the first and second terms of Vec2DAG, respectively.
+ 2 < 3. where the edge directions and connectivities are
Using this approach, it is simple to sample candidate DAGs whose representations are in the neigh- borhood around some values p and E of interest.
4

---

Published as a conference paper at ICLR 2025
Algorithm The DrBO method for causal discovery.
Require: Dataset D {x(i) [ j=1 of d nodes and observations, score function S (D,
DAG rank k, batch size B, no. preliminary candidates C, and total no. of evaluations T.
Ensure: A DAG that maximizes S gD,
1: Initialize empty experience H. := and node-wise dropout neural nets:
[DropoutNN,)
2: while [H| <
do
3: Generate random DAGs: (g (2())) where [-1,
Secs. 4.1 & 4.2.
(1)
DropoutNN, (pagto)
C
4: Sample local scores:
Sec. 4.3.
1=1
C
AFF = Combine (t)
5: Combine local scores:
Sec. 4.4.
6: Select top B candidates with highest AF values: J1.... jB == argtop g
AF
Sec. 4.2.
7: Evaluate these candidates and update experience: : HU HU s (D, g
j...j Sec. 4.5.
8: Update the neural nets H.
9: end while
10: Get highest-scoring DAG so far: 9 == arg maxgeR (D,g).
Prune I if needed.
Sec. 4.6.
4 DrBO: DAG RECOVERY VIA BAYESIAN OPTIMIZATION
An overview of our framework is illustrated in Algorithm 1. In the following, describe the proposed DrBO algorithm step-by-step.
4.1 SEARCH SPACE
To effectively utilize BO, the search space should be unconstrained. Thus, following et al. (2021); Massidda et al. (2024); Duong et al. (2025), transform the constrained combinatorial optimization problem in Eq. (1) to an unconstrained optimization task. However, the dimensionality of O (d²) of their search spaces ality, facilitating easier acquisition function optimization. For this purpose, similarly to Fang et al. (2023), we assume that the edge potential matrix in Vec2DAG (Duong et al., 2025) is low-rank and thus consider an adaptation that offers a search space that grows linearly with the number of nodes. Specifically, each node is now associated with dindenensional embedding vector r; E R* with k < d, and two nodes i and are connected if and only if (ri,r)) > O. The total dimensionality of this search space is thus only d (1 + k). More formally, given a node potential E Rd and an embedding matrix R. RRdk we define the following map
still be reduced to mitigate the effect of the curse of dimension-
T (p, R) ;= (grad (p)) @ H (R.R.).
(4)
The following Lemma ensures the acyclicity of the DAG corresponding to r (p, R). Lemma 1. For all d, k N e Ra and R E Raxk, let Rd x dd + {0, 13 axd be defined as Eq. (4). Then, 7 (p. R.) represents a binary adjacency matrix of a DAG.
The proof can be found in Appendix A.1. In addition, like Vec2DAG, our variation also exhibits the following scale-invariance property.
Lemma 2. For all d. k E N P Rd and R E Raxk let Rdd x Baxk - +0, 100,1aaxd be defined as in Eq. (4). Then, for all a > 7 Rp.
(op, oR).
The proof is provided in Appendix A.2. This insight allows us to restrict the search domain to fixed range (e.g., [-1, 1)) for numerical stability. For brevity, the remaining parts of this manuscript will use vector z of d' (1 + k) dimensions as the concatenation of and the flattened R, and we adopt the notation (z) 7 (p, R). short, we can now translate the original optimization problem in Eq. (1) to the following unconstrained optimization problem:
arg max S (D, (z)). (E
(5)
5

---

Published as a conference paper at ICLR 2025
Obviously, if & 2 d then R R is full-rank, Duong et al., 2025), and thus G := T (zz") is a maximizer of Eq. (1) for any maximizer z" of Eq. (5). While this may not hold for k < d. our empirical evaluations reveal that this representation suffices even for very complex graphs. Further, we find that lower ranks typically result in better sample-effciieny than higher ranks (Figure 3(a)) and provide an explanation in Appendix F.1.6.
represent any DAG possible (Theorem 1,
4.2 ACQUISITION FUNCTION OPTIMIZATION
To perform each step of the BO pipeline (Sec. 3.1), use the acquisition function obtained far to select a batch of most promising candidates to evaluate, also known as Batch BO (Joy et al,, 2020). This done by first generating C B preliminary candidates {x(1) | from hypercube centered at the current best solution z (see Appendix C.1 for more details). Subsequently, we evaluate the acquisition function of the DAGs induced by these candidates, and choose the top B candidates based on the acquisition function values. The quality of the candidates batch strongly depends C, i.e., if we evaluate the acquisition function of lot of candidates, then the top B candidates are likely to have higher values. However, this also increases computational cost, so acquisition function must scale well with the number of candidates C to mitigate this overhead, which leads us to the next point.
4.3 SURROGATE MODELING WITH DROPOUT NETWORKS
To overcome the scalability issues of standard BO due to the use of GPs as discussed earlier, we in- stead pursue neural networks, which well-known for their scalability and flexibility (Snock et al., 2015). Our networks must be able to model uncertainty to help the optimizer prioritize evaluating uncertain but promising candidates. Towards this end, we employ dropout activations (Srivastava et al.,, 2014), whose original purpose was to reduce overfitting in training neural networks, and were later found to be also useful an approximate Bayesian inference method (Gal & Ghahramani, 2016), and thus have been successfully applied to BO (Guo et al., 2021). We provide detailed discussion on this choice compared with other models in Appendix C.2.
Specifically, we devise a single-layer neural network with dropout activation as follows. Let p E (0, 1) be the dropout rate, d denotes the dimensionality of the input, is the number of hidden units, W1 E Rd and W2 E pax1 are weight matrices, b E I2t and b2 E R are biases. Our dropout networks are then defined as:
DropoutNN (x) W₂ BatchhNrmm RoLU (1-p (1- m) (W]x+ bb)))++))))))>))) + b2,
1
where we also follow common practice to employ Batch Normalization (loffe & Szegedy, 2015) for improved training efficiency, and m Bernoulli (p)* is drawn for every invocation of Dropout (x) in both train and test modes. By training this model on observed data X and y with the square loss, performing stochastic forward pass y Dropout (x) be interpreted as drawing from an approximate posterior ~ q (y x, X, (Gal & Ghahramani, 2016). Us- ing Thompson sampling, we do not need to characterize the whole posterior and this mere sample suffices for acquisition function optimization (Russo & Van Roy, 2014; Eriksson et al., 2019).
4.4 INDIRECT SURROGATE MODELING
A naive surrogate modeling approach for DAG learning is modeling direct map from a DAG to its score. However, DAG scores can typically be decomposed into independent node-wise components, which can be further exploited for better modeling. For example, the BIC scores in Eqs. (2) and (3)
involve the local components {MSE, {paf)}= observable after each DAG score invocation. To exploit these information to the fullest, we propose to learn local surrogate models predicting the node-wise scores, then combine them using the rule in Eq. (2) or (3), depending on the situation.
the evaluation data {paçf .SSE,
g()) (pa,")))
where gl) is the j-th evaluated DAG, to train dropout network DropoutNN, predicting In MSE, from paf. which is represented by binary vector of d dimensions. The models are independent among
Particularly, for each node i,
'Our AFs do not take as input z, but T (z) instead, many z's can produce the same DAG.

---

Published as a conference paper at ICLR 2025
all nodes instead of being shared to avoid spurious correlations. To summarize, for dataset of nodes, we jointly train local surrogate models [Dropout.NN,) i=- To sample the DAG score of a graph G. we first sample each local score I DropoutNN, (pa then combine them with AF (G) == CombineBIc-NY (41,...,ld) == n Ei= 4, | In if the non-equal variance BIC score is being considered, AF (g) CombineRIc-BE EV (11 (l,,.... ........ la) -nd Σin elt -I9] In for the case of equal variance BIC score, which resemble Eqs. (2) and (3), respectively.
4.5 CONTINUAL MODEL TRAINING
Upon acquiring batch of B of new evaluations, the neural networks require retraining to update their weights. avoid retraining with all data so far, which scales at least quadratically with the number of evaluations, we adopt continual training approach (Wang et al., 2024) that only updates the models with the new data and a small portion of past data to mitigate forgetting. Specifically, perform Rggrads gradient steps within each BO iteration, each of which is calculated on the B new datapoints and a replay buffer of Preplay past observations using Reservoir Sampling (Vitter, 1985).
4.6 FINALIZING THE RESULT
Pruning the resultant DAG is common practice to suppress the redundant edges (Bühlmann et al., 2014; Zheng et al., 2018; Wang et al., 2021; Bello et al., 2022; Duong et al., 2025), and is also employed in our framework. This can be done by thresholding the weight matrix for linear data (Zheng et al., 1)), or employing significance testing for nonlincar data using generalized additive model regression (CAM pruning, Peters al., 2014), conditional independence testing under the faithfulness assumption (Duong et al., 2025). More details are provided in Appendix C.3.
EXPERIMENTS
In this section, we verify our claim in the introduction: DrBO both more accurate and sample- efficient than existing approaches in score-based observational DAG learning. We show this by comparing our DrBO method with a number of the most recent advances in causal discovery that are based on sequential optimization, including gradient-based methods DAGMA (Bello et al., 2022), COSMO (Massidda et al., 2024). GOLEM (Ng et al., 2020), NOTEARS (Zheng et al., 2018) with TMPI constraint (Zhang al., 2022), well as RL-based approaches CORL (Wang et al., 2021) and ALIAS (Duong et al., 2025). We note that CORL, ALIAS, and DrBO directly optimize the BIC score, COSMO, DAGMA, and NOTEARS optimize the penalized least-square loss, while GOLEM optimizes a penalized log-likelihood. For gradient- methods, we consider a gradient update equivalent to one DAG evaluation. Additional information, including implementation details and metrics, are provided in Appendix D.
5.1 RESULTS SYNTHETIC DATA
We consider the standard Erdōs-Rényi (ER) graph model (Erdōs & Rényi, 1960) to generate data, where graphs with d nodes and de edges on average are referred to as dERe graphs (e.g, 10ER4).
5.1.1
LINEAR-GAUSSE DATA
After simulating DAG G, we sample edge weights with wys ~ 14 (1-2,-0.3) U [0.5,2) like prior studies (Zheng et al., 2018; Zhu et al., 2020; Wang et al., 2021). Then, we generate a dataset of 1,000 i.i.d. samples according to linear-Gaussin SCM := Σjepa, WydC, + Ex, where E N (0, 1). For fairness, we prune the DAGs returned by DrBO, along with ALIAS and CORL as prescribed in their papers, by thresholding the weight matrix obtained via linear regression at 0.3. This is not done for DAGMA and COSMO because their implementations already incorporated the same pruning scheme.
Dense graphs. We stress-test our method with complex structures in Figure 1(a). As depicted in the first column, our method is the only approach that can achieve absolute overall performance in all five metrics, surpassing the second- and third-best methods ALIAS and CORL by large margins of more than 20% in each metric, while gradient-based methods DAGMA and COSMO struggle with much worse performance. From the second and third columns, our DrBO approach can reach higher BIC scores, and thus lower SHDs, very sharply with the number of evaluations, highlighting
7

---

Published as a conference paper at ICLR 2025
1 - FDR
BIC score (higher is bette))
SHD (lower is betzer)
SHD (lower is better)
176.8
181.6
152.2
F
177.8
200
153.7
-2
145.4
130.8
105.2
116.4
Precision
Recall
1-6
2.0
20000 40000
25000 50000
0.0
2.5
5.0
DAG evaluations
DAG evaluations
Minutes
Dr80 (Ours)
ALJAS
CORL
COSMO
DAGMA
GOLEM
NOTRARS +
Dense graphs with Linear-Gaussian data (DAGs with 30 nodes and ==240 edges). For fairness,
summary metrics (first column) are calculated at 50,000 evaluations for all methods.
1 FDR
BIC score (higher is bette)
SHD (lower is better)
SHD (lower is bettes)
400
400
117.0
112.6 6
Fi
-1.00
104.9
2
-3200
11A
13.66
8.6
8..6
8.4
-4
-6.6
Precision
Recall
-14
1.4
20000 40000
25000 50000
20
DAG evaluations
DAG evaluations
Minutes
D:8O (Ours)
A11AS
CORL
COSMO
DAGMA
GOLEM
NOTEARS + TMMP
(b) Large graphs with Linear- Gaussian data (DAGs with 100 nodes and 200 edges). For
fairness, summary metrics (first column) are calculated at 50,000 evaluations for all methods.
FDR
BIC score (higher is better)
SHD (lower is better)
40
iis Anaate
TPR
F
25.6
29.2
20.0
17.0
8:4
-1.6
7.2
7.2
4.0
Precision Recaal
0
10000
20000
10000 20000
DAG evaluations
DAG evaluations
Minutes
DrBO (Ours)
AIJAS
CORL
COSMO
DAGMA
NOTEARS TMMP
(c) Nonlinear data with Gaussian processes (DAGs with 10 nodes and 4:00 edges). For fairness, summary metrics (first column) calculated at 20,000 evaluations for all methods.
Figure 1: DAG learning results on Synthetic data. First column: overall performance in terms of True Positive Rate (TPR, higher is better), Precision, Recall, and F1 score (higher better), as well as False Discovery Rate (FDR, lower is better). Second column: we track the best Bayesian Information Criterion (BIC. higher better) so far at every optimization step. Third and Fourth columns: we monitor the Structural Hamming Distance (SHD, lower is better) of the DAG whose best BIC so far at every optimization step. Shaded areas in the line plots indicate 95% confidence intervals over random datasets. NOTEARS+TMP! usually stops early before the time limit.
its sample-efficieney. Lastly, the fourth column shows that our method's SHD improvement over time is continuous and faster in terms of runtime than all other baselines.
Large-scale graphs. Next, we demonstrate the scalability of our method on high-dimensions graphs of 100 nodes, having up to 10,000 edges. Our results in Figure 1(b) shows that DrBO is still the leading method for high-dim data, where it obtains absolute overall performance along with the lowest SHD among all methods in limited runtime.
5.1.2 NONLINEAR DATA
Following prior works (Zhu et al., 2020; Wang et al., 2021; Yang et al., 2023a; Duong et al., 2025), demonstrate Lachapelle et al. (2020). Specifically, we employ the datasets generated according to an ANM 24 = Si (Xpa,) + E; where fi is drawn from Gaussian process prior with unit bandwidth RBF
method nonlinear datasets generated using Gaussian processes from
8

---

Published as a conference paper at ICLR 2025
kernel, and E .V (f) with non-equal noise variances of sampled uniformly in [0.4, 0.8]. ALIAS, CORL, and DrBO optimize SBICNVV, while COSMO, DAGMA, and NOTEARS+TMP! optimize the least-square objective. The empirical results reported in Figure I(c) confirms that our DrBO method also excels on complex nonlinear data in both accuracy and computational cost. This is evidenced by the leading overall performance (first column), along with vanishing SHD of only 0.4 at the end of the learning curve, surpassing other methods by a visible gap both in SHD and convergence speed.
5.2 RESULTS ON REAL DATA AND STRUCTURES
Table I: Causal Discovery Performance on Real-world Structures (Scutari, 2010). The per- formance is measured Structural Hamming Distance (SHD, lower is better). The numbers are mcan + std over 5 independent datasets with 1,000 observational samples. For faimesss all methods are limited 20,000 evaluations.
Dataset
Alarm
Asia
Cancer
Clold
Earthquake (5 modes, 4 edges)
Method
(37 uodes, 46 edges)
18 modes,, N edges)
(5 modes, 4 edges)
(20 uodes, 25 edges)
ALIAS [Doong at. 2025)
25.8 1. 7.8
0.2 0.3
0.0 0.0
288
0.0 0.0
CDRL (Wang et al., 2021)
19.8
0.0
1.6 2.3
0.0 0.0
COSMO (Massidla et al., 2024)
25.8 =: 4.0
4.3 1.6
2:2 0.S
11.64:2..
2:2 2 11S
DAGMA, (Bello et al., 2012)
25.0 +: 5.3
2.5 1.3
1.0 10:1:12
8.0 5.2
1,0 20:2:1
GOLEM (Ng et al. 2020)
4.8 1: 6.2
0.0 0.00
0.0 0.0.00
0.0 0.0
0.0 0:0
NOTEARS-TMPI (Zheng et al.,, 2018; Zhang et 2022)
7.0 7.9
0.4.1.0.9
0.0.1.00
0.0 0.0
0.0 0.0
DrBO (Ctrs)
1.0
0.0++0000
0.0+00.
0.0+0000
0.0++00
Benchmark data. We verify the performance of our method on real data using the popular benchmark flow cy- tometry dataset (Sachs et al., 2005), concerning a protein signaling network based expression levels of proteins and phospholipids. We employ the observational portion of the dataset containing 853 observations and known causal network with 11 nodes and 17 edges. We apply
ALIAS
DAGMA
CORL
NOTEARS-TMPI
COSMO
-
DrBO (Ours)
SHD (lower is better)
16
the similar settings in Sec. 5.1.2 for all methods. The evaluations shown in Figure verify the effectiveness of
our method on real data, where effortlessly achieves
10
lowest SHD of 9 with fewer evaluations compared with the competitors.
5000
10000
15000
DAG evaluations
Real-world structures. To further illustrate the capabili- ties of our approach on real-world scenarios, we conduct experiments on real structures provided by the Bnl-earn repository (Scutari, 2010). Each dataset contains 1,000 observational samples and a ground truth causal network belonging to real-world applications with varying size. Additional details regarding these datasets are given in Appendix D.3. The results are presented in Table 1, highlighting that our method is the only approach that can consistently achieve zero SHD on four out of five real-world structures. On the Alarm dataset, which appears to be most challeng- ing, our DrBO method still leads with much lower SHD compared with all other baselines.
Figure 2: Causal Discovery Per- formance on the Benchmark Sachs Dataset (Sachs et al., 2005). NOTEARS+TMPI stops early before the max no. of evaluations is reached.
5.3
SUPPLEMENTARY RESULTS
Extended Causal Discovery Settings. We investigate the performance of DrBO in extended sce- narios in Appendix F as follows. Varying sample sizes: we show in Figure 4 that method can achieve low SHDs even with limited data. Different graph models: in Table 3, DrBO achieves low SHDs and surpasses competitors for both ER and SF graphs, even on the dense graphs. Dif- ferent noise distributions: Table 4 shows that DrBO also outperforms the baselines under five different noise types. BGe score: Figure 5 confirms that our method can work with different score and match the score of ground truth graphs with low structural errors. Discrete Data: in Figure 6, we show that DrBO also obtains the highest scores and lowest SHDs compared with the baselines for non-continuous data. Standardized Data: we show in Figure 7 that our method is also robust against data standardization for both linear and nonlinear data. Large-scale Nonlinear Data: Fig- ure 8 demonstrates DrBO's competitive performance and efficiency for nonlinear data on 50- and 100-node graphs.
Additional baselines. In Appendix G, also examine other baselines that are not based on quential optimization, showing that DrBO also significantly outperforms these methods on linear, nonlinear, well as real data.

---

Published as a conference paper at ICLR 2025
250
Dropout Nets
Training time
Full rank (465 dims)
SHD (lower better)
Rank 32 (990 dims)
Exact GPs
Sampling time
200
Rank 12 (390 dims)
Approx. GPs
150
Rank (270 dims)
time
Cumulative runtime
Rank (150 dims)
Seconds per step
100
Rank dims)
300
Seconds
200
50
100
0-
2000
2000 DAG evaluations
DAG evaluations
DAG evaluations
(a) Effeet of Low-Rank Representation. We compare this for different ranks & with the Full-rank representation (Duong et al., 2025).
(b) Effect of Dropout Networks. We compare these with Exact GPs (Rasmussen, 2003) and Approximate GPs (AGP, Hensman et al.. 2015). Details in Appendix 1.211.
See our discussion in Appendix F.1.6.
Direct (2048 units)
Direct [128 units)
Direct (1024 unies)
Direct (64 units)
Full retraining
Continual training
Direct (512 unins)
Indirect (64 units)
BIC prediction RMSE
Training time
Cumulative runtime
Seconds per step
10
100
200
SHD
150
Minutes
10a
100
10a
10000 20000
0
10000
20000
20000
40000
0
20000
40000
DAG evaluations
DAG evaluations
DAG evaluations
DAG evaluations
(c) Effect of Indirect DAG Modeling, We compare this with learning Direct map from DAG to score, with varying no. hidden units.
(d) Effect of Continual Training. We compare this approach with Full retraining using all experienced data at every BO iteration. Details are in Appendix F.2.2.
Figure 3: Ablating our design choices. All configurations are evaluated on 5 linear- Gaussian datasets of 1,000 samples 30ER8 graphs. Shaded areas indicate 95% confidence intervals.
Runtime, In Appendix H, we also detail the numerical runtime among all methods.
5.4 ABLATIONS
In Figure 3, through ablation studies, we justify the design choices outlined in the Introduction, showing that every component contributes considerately to the accuracy and/or scalability of our method. The remaining hyperparameters of algorithm are also studied in Appendix F.2.
CONCLUDING REMARKS
This study presents DrBO, novel BO method search for high-scoring DAGs. We have shown that, by meticulously choosing promising DAGs to evaluate, we can find the optimal one more efficiently and cost-effectively. Our comprehensive experiments demonstrate that DrBO performs well even in many intricate settings like dense graphs, high-dimensional,, and nonlinear data.
Regarding limitations, our method, the surrogate model architecture is manually chosen and mains fixed through the course of optimization, and thus is prone underfitting at the end. To mitigate this, incremental neural architecture search techniques (Liu et al., 2018: Geifman & El-Yaniv, 2019) can be employed to facilitate autonomous architecture selection and scaling. In addition, tinual learning component simple and can benefit from more advanced techniques (Wang et al., 2024) to further improve the performance of surrogate models.
For future developments, it would be an interesting direction to combine our approach with ac- tive causal discovery methods recover the causal structure more efficiently. Morcover, method can be extended to solve causal discovery problems with hidden confounders, where the outputs longer DAGss
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Virginia Aglictti, Xiaoyu Lu, Andrei Paleyes, and Javier González. Causal Bayesian optmizzttion In Proceedings of the International Conference on Artificial Intelligence and Statistics, pp. 3155- 3164. PMLR, 2020.
Virginia Aglietti, Neil Dhir, Javier González, and Theodoros Damoulas. Dynamic causal bayesian optimization. In Advances in Neural Information Processing Systems, pp. 10549-10560, 2021.
Yashas Annadani, Nick Pawlowski, Joel Jennings, Stefan Bauer, Cheng Zhang, and Wenbo Gong. BayesDAG: Gradient-based posterior sampling for causal discovery. In Advances in Neural In- formation Processing Systems, 2023.
Albert-Lászin Barabási and Rćka Albert. Emergence of scaling in random networks. Science, 286: 509-512, 1999.
Ingo Beinlich, Henri Jacques Suermondt, R Martin Chavez, and Gregory Cooper. The ALARM monitoring system: A case study with two probabilistic inference techniques for belief networks. In Proceedings of the European Conference on Artificial Intelligence Medicine, pp- 247-256,
1989.
Kevin Bello, Bryon Aragam, and Pradeep Ravikumar. DAGMA: learning DAGs via M-matrices and a reg-determinant acyclicity characterization. In Advances in Neural Information Processing Systems, pp. 8226-82333 2022.
Philippe Brouillard, Sébastien I Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin. Differentiable causal discovery from interventional data. Advances Neural Information Processing Systems, 33:21 865-218877 2020.
Peter Bühlmann, Jonas Peters, and Jan Ernest. CAM: Causal additive models, high-dimensional order search and penalized regression. The Annals of Statistics, 42:2526 2556, 2014.
David Maxwell Chickering. Learning Bayesian networks is NP-complete. Learning from data: Artificial intelligence and statistics V, pp. 121-130, 1996.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3:507-554, 2002.
Diego Colombo, Marloes H Maathuis, Markus Kalisch, and Thomas S Richardson. Learning high- dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics, pp. 294-321, 2012.
Samuel Daulton, David Eriksson, Maximilian Balandar, and Eytan Bakshy. Multi-objective bayesian optimization over high-dimensional search spaces. In Uncertainty in Artificial Intelligence, pp- 507-517. PMLR, 2022.
Tristan Deleu, António Góis, Chris Emezue, Mansi Rankawat, Simon Lacoste-Julien, Stefan Bauer, and Yoshua Bengio. Bayesian structure learning with generative flow networks. In Proceedings of the Uncertainty in Artificial Intelligence, pp. 518-528, 2022.
Tristan Deleu, Mizu Nishikawa- Toomey, Jithendaraa Subramanian, Nikolay Malkin, Laurent Char- lin, and Yoshua Bengio. Joint Bayesian inference of graphical structure and parameters with a single generative flow network. In Advances in Neural Information Processing Systems, 2024.
Bao Duong, Hung Le, Biwei Huang, and Thin Nguyen. Reinforcement learning for causal discovery without acyclicity constraints. Transactions Machine Learning Research, 2025.
Paul Erdös and Alfréd Rényi. On the evolution of random graphs. Publications of the Mathematical Institute of the Hungarian Academy of Sciences, 1960.
David Eriksson. Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization. Advances in Neural Information Processing Systems, 32, 2019.

---

Published as a conference paper at ICLR 2025
Zhuangyan Fang, Shengyu Zhu, Jiji Zhang, Yue Liu, Zhitang Chen, and Yangbo He. On low-rank directed acyclic graphs and causal structure learning. IEEE Transactions on Neural Networks and Learning Systems, 35(4):4922-49)) 4937, 2023.
Yarin Gal and Zoubin Ghahramani. Dropout as a bayesian approximation: Representing model un- certainty in deep learning. In Proceedings of the International Conference on Machine Learning, pp. 1050-1059, 2016.
Yonatan Geifman and Ran EI- Yaniv.. Deep active learning with a neural architecture search. In Advances in Neural Information Processing Systems, 2019.
Dan Geiger and David Heckerman. Learning Gaussian networks. In Uncertainty Proceedings, pp- 235-243. 1994.
Dan Guo, Yaochu Jin, Jinliang Ding, and Tianyou Chai. Heterogeneous ensemble-based infill crite- rion for evolutionary multiobjective optimization of expensive problems. IEEE Transactions Cybernetics, 49(3):1012-102, 2018.
Dan Guo, Xilu Wang, Kailai Gao. Yaochu Jin, Jinliang Ding. and Tianyou Chai. Evolutionary optimization of dimensional multiobjective and many-objective expensive problems assisted by dropout neural network. IEEE transactions systems, man, and cybernetics: systems, 52: 2084-2097, 2021.
Dominique Marie-Annicl Haughton. On the choice of a model to fit data from an exponential family. The Annals of Statistics, pp. 342-355, 1988.
Alain Hauser and Peter Bühlmann. Characterization and greedy learning of interventional Markov equivalence classes of directed acyelic graphs. Journal of Machine Learning Research, 13(1): 2409-2464, 2012.
David Heckerman, Dan Geiger, and David M Chickering. Learning Bayesian networks: The com- bination knowledge and statistical data. Machine Learning, pp. 197-243, 1995.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian pro- cess classification. In Proceedings of the International Conference on Artificial Intelligence and Statistics, pp. 351-360. PMLR, 2015.
Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems, 2008
Frank Hutter, Holger H Hoos, and Kevin Leyton-Brown, Sequential model-based optimization for general algorithm configuration. In Learning and Intelligent Optimization: 5th International Conference, LION 5, Rome, Italy, January 17-21, 2. Selected Papers 5, pp. 507-523. Springer, 2011.
Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In Proceedings of the International Conference Machine Learning. 2015.
Tinu Theckel Joy, Santu Rana, Sunil Gupta, and Svetha Venkatesh. Batch Bayesian optimization using multi-scale search. Knowledge-Bass Based Systems, 187:104818, 2020.
Daphne Koller and Mehran Sahami. Toward optimal feature selection. Technical report, 1996.
Kevin B Korb and Ann E Nicholson. Bayesian artificial intelligence. CRC press, 2010.
Jack Kuipers, Giusi Moffa, and David Heckerman. Addendum the scoring of Gaussian directed acyclic graphical models. The Annals of Statistics, 42(4):1689 1691, 2014.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-basse neural DAG learning. In Proceedings of the International Conference on Learning Representa- fions, 2020.
12

---

Published as a conference paper at ICLR 2025
Steffen L Lauritzen and David Spiegelhalter. Local computations with probabilities graphical structures and their application to expert systems. Journal of the Royal Statistical Sociery: Series B (Methodological), 50(2):157-194, 1988.
Lek-Heng Lim. Hodge Laplacians on graphs. Siam Review, 62:685-715, 2020.
Phillip Lippe, Taco Cohen, and Efstratios Gavves. Efficient neural causal discovery without acyclic- ity constraints. In Proceedings of the International Conference on Learning Representations, 2022.
Chenxi Liu, Barret Zoph, Maxim Neumann, Jonathon Shlens, Wei Hua, Li-Jia Li, Li Fei-Fei, Alan Yuille, Jonathan Huang, and Kevin Murphy. Progressive neural architecture search. In Proceed- ings of the European Conference on Computer Vision, pp. 19-34, 2018.
Riccardo Massidda, Francesco Landolfi, Martina Cinquini, and Davide Bacciu. Constraint freee structure learning with smooth acyelic orientations. In Proceedings of the International Confer- ence on Learning Representations, 2024.
Michael D McKay, Richard J Beckman, and William Conover. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 42(1))55-61, 2000.
Francesco Montagna, Nicoletta Noceti, Lorenzo Rosasco, Kun Zhang, and Francesco Locatello. Scalable causal discovery with score matching. In Conference on Causal Learning and Reason- ing. pp. 752-771. PMLR, 2023.
Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and DAG constraints for learning linear DAGs. In Advances in Neural Information Processing Systems, pp. 17943-17954, 2020.
Art B Owen. Scrambling Sobol' and Niederreiter-Xing points. Journal of complexity, 14(4):466- 489, 1998.
Judea Pearl. Models, reasoning and inference. Cambridge University Press, 2000.
Judea Pearl. Causality. Cambridge University Press, 2009.
Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 2014.
Jonas Peters, Dominik Janzing, and Berhard Schölkopf. Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
Joseph Ramsey, Madelyn Glymour, Ruben Sanchez- Romero, and Clark Glymour. A million vari- ables and more: the Fast Greedy Equivalence Search algorithm for learning nigh-dimensional graphical causal models, with application to functional magnetic resonance images. Interna- tional Journal of Data Science and Analytics, 3:121-129, 2017.
Joseph D Ramsey. Scaling up greedy causal search for continuous variables. arXiv preprint 11Xiv:1507.07749. 2015.
Carl Edward Rasmussen. Gaussian processes in machine learning. In Summer school on machine learning, pp. 63-71. Springer, 2003.
Rommel G Regis and Christine A Shoemaker. Combining radial basis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization. Engineering Optimization, 45(5):529-555, 2013.
Alexander Reisach, Christof Seiler, and Sehastian Weichwald. Beware the simulated dag! causal discovery benchmarks may he easy to game. In Advances in Neural Information Processing Systems, volume 34, pp. 27772-27784, 2021.
Robert W Robinson. Counting unlabeled acyelic digraphs. In Combinatorial Mathematics V, pp- 28 43. Springer, 1977.
13

---

Published as a conference paper at ICLR 2025
Paul Rolland, Volkan Cevher, Matthäus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score matching enables causal discovery of nonlinear ad- ditive noise models. In Proceedings the International Conference Machine Learning, pp-
18741-18753, 2022.
Daniel Russo and Benjamin Van Roy. Learning to optimize via posterior sampling. Mathematics of Operations Research, 39(4):121-12243 2014.
Karen Sachs, Omar Perez, Dana Pe'er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 308:523-529, 2005.
Pedro Sanchez, Xiao Liu, Alison O' Neil, and Sotirios A. Tsaftaris. Diffusion models for causal discovery via topological ordering. In Proceedings of the International Conference on Learning Representations, 2023.
Gideon Schwarz. Estimating the dimension of model. The Annals of Statistics, pp. 161-46116666 1978.
Marco Seutari. Learning Bayesian networks with the balearn R package. Journal of Statistical Software, 35:1-22, 2010.
Edward Snelson and Zoubin Ghahramani. Sparse gaussian processes using pseudo-inputs. Advances in neural information processing systems, 18, 2005.
Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary, Mr Prabhat, and Ryan Adams. Scalable bayesian optimization using deep neural networks. In Proceedings of the International Conference Machine Learning, pp. 2171- 2180. PMLR, 2015.
David J Spiegelhalter. Learning probabilistic expert systems. Bayesian statistics, 4:447-46 1992.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT Press, 2000.
Jost Tobias Springenberg, Aaron Klein, Stefan Falkner, and Frank Hutter. Bayesian optimization with robust Bayesian neural networks. Advances in Neural Information Processing Systems, 29, 2016.
Niramjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger. Gaussian process opti- mization in the bandit setting: No regret and experimental design. In Proceedings of the Interna- tional Conference on Machine Learning, 2010.
Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout: a simple way to prevent neural networks from overfitting. Journal of Machine Learning Research,
2014.
William R Thompson. On the likelihood that one unknown probability exceeds another view of the evidence of samples. Biometrika, 25(3-4):285-294, 1933.
Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Pro- ceedings of the Artificial Intelligence and Statistics, pp. 567-574. PMLR, 2009.
Christian Toth, Lars Lorch, Christian Knoll, Andreas Krause, Franz Pernkopf Robert Peharz, and Julius Von Kügelgen. Active bayesian causal inference. In Advances Neural Information Processing Systems, pp. 16261-16275, 2022.
Quang-Duy Tran, Phuoe Nguyen, Bao Duong, and Thin Nguyen. Differentiable bayesian structure learning with acyclicity assurance. In Proceedings of the IEEE International Conference Data Mining, pp. 598-607, 2023.
Paul E Utgoff. Incremental induction of decision trees. Machine learning, 4:161-186, 1989.
Paul E Utgoff, Neil C Berkman, and Jeffery A Clouse. Decision tree induction based on efficient tree restructuring. Machine Learning, 29:5-44, 1997.
14

---

Published as a conference paper at ICLR 2025
Jeffrey S Vitter. Random sampling with reservoir. ACM Transactions Mathematical Software (TOMS), 11(1):37-57, 1985.
Liyun Wang, Xingxing Zhang, Hang Su, and Jun Zhu. A comprehensive survey of continual learning: theory, method and application. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2024.
Xiaoqiang Wang, Yali Du, Shengyu Zhu, Liangjun Ke, Zhitang Chen, Jianye Hao, and Jun Wang. Ordering-based causal discovery with reinforcement learning. In Proceedings of the International Joint Conference on Artificial Intelligence, pp. 3566-3573, 2021.
Xilu Wang, Yaochu Jin, Sebastian Schmitt, and Markus Olhofer. Recent advances in Bayesian optimization. ACM Computing Surveys, 55:1-36, 2023.
Zi Wang. Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka. Batched large-scale Bayesian optimization in high-dimensional spaces. In Proceedings of the International Conference on Artificial Intelligence and Statistics, pp. 745-754. PMLR, 2018.
Eric P Xing, Michael I Jordan, Richard M Karp, et al. Feature selection for high-dimensional genomic microarray data. In Proceedings of the International Conference on Machine Learning. volume 1, pp. 601-608. Citeseer, 2001.
Dezhi Yang, Guoxian Yu, Jun Wang, Zhengtian Wu, and Maozu Guo. Reinforcement causal struc- ture learning on order graph. In Proceedings of the AAAI Conference on Artificial Intelligence, pp. 10737-10744, 2023a.
Dezhi Yang. Guoxian Yu, Jun Wang, Zhongmin Yan, and Maozu Guo. Causal discovery by graph attention reinforcement learning. In Proceedings of the SIAM International Conference on Data Mining, pp. 28-36, 2023b.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In Proceedings of the International Conference Machine Learning, pp. 7154-7163, 2019.
Yue Yu, Tian Gao, Naiyu Yin, and Qiang Ji. DAGS with No Curl: An efficient DAG structure learning approach. In Proceedings of the International Conference on Machine Learning, pp.
12156-12111 12166, 2021.
Keli Zhang, Shengyu Zhu, Marcus Kalander, Ignavier Ng, Junjian Ye, Zhitang Chen, and Lujia Pan. gCastle: A python toolbox for causal discovery. arXiv preprint arXiv:2111.15155 2021.
Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kernel-based conditional independence test and application in causal discovery. In Proceedings of the Conference Un- certainty in Artificial Intelligence, pp.
2011.
Zeyu Zhang, Chaozhuo Li, Xu Chen, and Xing Xie. Bayesian active causal discovery with multi- fidelity experiments. In Advances in Neural Information Processing Systems, 2024.
Zhen Zhang, Ignavier Ng. Dong Gong, Yuhang Liu, Ehsan Abbasnejad, Mingming Gong, Kun Zhang, and Javen Qinfeng Shi. Truncated matrix power iteration for differentiable DAG learning. In Advances in Neural Information Processing Systems, pp. 18390-18402, 2022.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. DAGs with NO TEARS: Contin- uous optimization for structure learning. In Advances in Neural Information Processing Systems, pp. 9472-9483, 2018.
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. Learning sparse nonparametric DAGs. In Proceedings of the International Conference on Artificial Intelligence and Statistics, 2020.
Shengyu Zhu, Ignavier Ng, and Zhitang Chen. Causal discovery with reinforcement learning. In Proceedings of the International Conference on Learning Representations, 2020.
15

---

Published as a conference paper at ICLR 2025
A PROOFS
A.1 PROOF OF LEMMA 1
Proof. The acyclicity of Eq. (4) is ensured by the first term: H (grad (p)). This binary matrix representing the adjacency matrix of a directed graph. In this graph, the presence of an edge i - equivalent to P: < Pj- As such, cycle i - ... - i1, if exists, would lead to p; < Pis which contradictory, meaning this graph must be DAG. Multiplying this adjacency matrix with the second term H (R- RR). which also a binary matrix, has the effect of removing already existing edges in this graph, and thus cannot introduce any cycle. This concludes our proof.
A.2 PROOF OF LEMMA 2
Proof. For any > 0, we have H (p; -Pj) = H (a (p. - Pj)) = H (ap: app) and H RR. R = H R = H (aR) (aR) Thus, H (grad (p)) o (R (R-R =
H (grad (ap)) @ H (aR) (aR)
concluding our proof.
B DERIVING BIC SCORES
P(e)). . the general BIC is given by Recall that for causal model with parameters 0 = {(ft)
SEIC (D,9) = 2mp (D (0) | Inn,
where 0 : arg maxa (D) 9) is the maximum-likelihood estimator of the causal model param- eters, n is the sample size of D, and | denotes the number of edges in G. The first term in Eq. (6) is ligeliillhood objective similar to GOLEM (Ng et al., 2020), GraN-DAG (Lachapelle et al., 2020), etc., while the second term penalizes extra edges.
The flexibility of Eq. (6) allows for the adoption of BIC in various causal models by simply specify- ing the likelihood model. In the following, we derive the BIC scores for additive noise models with non-equal and equal variances, as well as logistic models.
B.
BIC FOR ANM WITH NON-QUAL VARIANCES
Let us consider an ANM defined by
I; == fo (xpaz) + Eis =i= I,..,.., d,
where the noise is assumed to be Gaussian with fixed variance: ~ N (0,0?) Vi = 1,... d. The likelihood of dataset D {x())* } under this model given by
(2 24
(1 pag
fo
C= -- (D19, 1 (A) (or)
In of constant.
o?
(8)
Taking its derivative and setting it to zero, we can solve for the parameters as follows:
03 = 1
1 fe (xx1)) pap ))
(9)
MSE,
where {i]] ==1 are least-square estimators obtained by minimizing (2.00) f1 (x(2)
paf ))
over
restricted hypothesis class /~ F. Following past studies (Zhu et al., 2020; Wang et al., 2021; Yang
16

---

Published as a conference paper at ICLR 2025
et al., 2023a;b; Duong et al., 2025), linear regression for linear data and Gaussian process regression for nonlinear data. However, any other regression method can be used as desired.
Subsequently, the first term in Eq. (8) cancels out and the maximum log- likeliikd is reduced to
d
11
In MSE; + constant.
(10)
Finally, the BIC score for ANMs with non-equal variances is obtained by:
SBIC-NV (D,9) -1 In MSE, - |g
(11)
B.2 BIC FOR ANM WITH EQUAL VARIANCES
By setting 0; Vi = -.... d, repeat the previous steps and obtain the following solution for the maximum likelihood estimators:
d
8²
2
x(()
pa.
(12)
d
MSE.
The maximum likelihood now is given by
nd
MSE,
È 2
constant.
(13)
d
Finally, we obtain the BIC score for ANMs with equal variances follows:
Σ MSE;
SBIC (D,9) -nd In
[ n.
(14)
d
B.3 BIC FOR BINARY DATA WITH LOGISTIC REGRESSION
From the formulation in Eq. (6), we can also adapt it to non-continuous data. For example, let us consider the logistic causal model governed by
E3 ~ Bernoulli (fe (xxpat
(15)
where fa models the conditional probability of given xpa.-
The log-likelihood of data under this model is determined by
d
- L (DIG.())))) = ΣΣ (1) In (x ) + (1-1(1) In (1 (1- fo (x()
pag )))
j=
Maximizing this objective and plugging into Eq. (6) gives us with the BIC for logistic model:
Spic-Logistic (D,9) 2ć |g] Inn.
(16)
The first term is equivalent to the logistic loss used in, e.g., DAGMA (Bello et al., 2022), while the second term punishes redundant edges as always.
17

---

Published as a conference paper at ICLR 2025
C ADDITIONAL DISCUSSIONS AND ALGORITHM DETAILS
C.I PRELIMINARY CANDIDATE GENERATION
In Sec. 4.2, we generate set of C preliminary candidates in hyperrectangle centered at the best solution so far z This hyperrectangle can be seen as a single "trust region" in BO (Eriksson al., 2019; Daulton et al., 2022). More specifically, 1)1,1)((_ and the hypercube of length centered at z' Following Eriksson et al. (2019), we adaptively update L according to the learning progress. This is done by maintaining success (and failure) counter that keeps track the number of consecutive BO iterations that improves (or fails to improve, resp.) the DAG score. After Tsnce consecutive successes, we enlarge L by two times in order to shift the focus to other regions, and after Niail consecutive failures, we shrink it by two times to zoom more into the current region. In all experiments, Aaai = 5, and is initialized with the value of and is clipped to be within [0.01, 2] after each update.
trust region is the intersection of the search space
the fixed values of Nluxe = 3,
To produce the preliminary candidates on which we generate Thompson samples, for the very first suggestions, following Eriksson al. (2019), we employ the Latin hypercube design (LHD, McKay et al., 2000), which is a space-filiigg method used to generate near-random samples from a multi- dimensional space, which ensures that each dimension of the hypercube is evenly covered, while random sampling can lead to an unevenly covered space. For subsequent suggestions, we follow the established procedure in (Eriksson al., 2019) first generate a scrambled Sobol sequence (Owen, 1998) within the current trust region, then we use the perturbation value in the Sobol sequence with ā[[+))} for each given candidate and dimension, and the value of the center z" otherwise. As noted in (Regis & Shoemaker, 2013; Eriksson al., 2019), perturbing only a few dimensions can lead to significant performance improvement for high-dim scenarios.
probability min {
C.2 SCALABLE SURROGATE MODELING
The literature of BO is vast and here we only discuss a few promising alternative approaches to scale up surrogate models in BO, and justify of our dropout neural network choice.
Bayesian neural networks (BNN) are a also natural replacement thanks to the flexibility of neural networks combined with the inherent ability to model uncertainty of the Bayesian ideology (Sprin- genberg et al., 2016). However, to stay as close as possible to a truly Bayesian treatment, i.e., char- acterizing the exact posterior distributions, they require stochastic gradient Markov Chain Monte Carlo (MCMC) to sample from the posterior, which necessitates several sampling steps to reach the desirable posterior (Springenberg et al., 2016). Meanwhile, in method, sacrifice the accuracy of the posterior inference, so that we can sample from the posterior faster with a single forward pass through the dropout neural networks. That being said, our DAG learning accuracy still strong, which justifies this sacrifice.
Sparse GPs (or Approximate GPs) have also been studied to enhance scalability of GPs (Snelson & Ghahramani, 2005; Hensman et al., 2015; Titsias, 2009) by reducing the number of data points to set of "inducing points" that are representatives of the true data, which can be chosen to be much smaller than the original data set, However, it is observed that these methods usually do not scale well with dimensionality (Wang et al., 2018).
Random forest (RF) has also been considered as an alternative surrogate model in BO (Hutter et al., 2011). An RF is composed of multiple decision trees, each of which is constructed from portion of the training set and a few random dimensions. Therefore, uncertainty modeling in RF is achieved via the variation of the individual trees' predictions. While RF known to be efficient, especially in tabular data, we find it more straightforward to train neural networks continually (Wang et al., 2024) compared with tree-based models (Utgoff, 1989; Utgoff et al., 1997). This is crucial for our BO framework to scale with the number of iterations aforementioned.
Ensembling methods (Wang et al., 2018: Guo et al., 2018) sidestep the scalability challenge of BO with the use of ensemble of models trained different partitions of the samples and dimensions. A similar idea is also proposed in (Eriksson et al., 2019), where local GPs are trained on multiple local trust regions of the search space. However, these methods require training multiple models,
18

---

Published as a conference paper at ICLR 2025
and managing them is less straightforward than maintaining a single neural network like our method, which also scales very well with a strong causal discovery performance.
C.3 PRUNING TECHNIQUES
Pruning is typically employed in causal discovery to reduce false positive estimates, and there are several approaches depending on the causal model as follows.
Linear data. For linear data, given the resultant DAG & that needs to be pruned, linear regression is used to find the linear coefficients {úsg (i--))e associating with the edges in G. Then, only the weights satisfying certain absolute strength is kept, while the remaining are removed. In other words, the edge (i + j) E 9 is kept iff [ > a. Typically, the weights in the true generative process for linear models are sampled from U ([-2, -0.5] U [0.5, 2)), and the common value for is 0.3 (Zheng et al., 2018; Zhu et al., 2020; Wang et al., 2021; Bello et al., 2022; Massidda et al., 2024; Duong al., 2025).
Nonlinear additive model. For nonlinear data under additive models, a common approach is based on feature selection using generalized additive model (GAM) regression (Bühlmann et al., 2014), also known as CAM pruning. Particularly, each node i is regressed on its parents Xpaî using GAM, then the significance test of covariates is conducted, and parent is kept if its p-value is lower than the significance level of 0.001. This is also done in (Zhu et al., 2020; Wang et al., 2021; Massidda et al., 2024; Rolland et al., 2022; Sanchez et al., 2023; Duong et al., 2025).
Non-additive models. For general non-additive models, conditional independence (CI) testing can also be employed to prune edges. To be more specific, the Faithfulness assumption (Peters et al., 2017) implics that any conditional independence observed in data reflects the corresponding d- separation in the causal graph, so if x; 1. Tj Xpa,(f) for some j pa,, then j is extra edge and needs to be removed. This follows the same idea of constraint-based causal discovery (Spirtes et al., 2000; Colombo et al., 2012) and feature selection via Markov Blankets (Koller & Sahami, 1996; Xing et al., 2001). As observed in (Duong et al., 2025), Cl-based pruning leads to better performance on the Sachs dataset compared with CAM pruning, and hence we also employ for all methods on the Sachs dataset. The specific CI test is the popular kernel-based method KCIT (Zhang et al., 2011) with a significance level of 0.001.
D EXPERIMENT DETAILS
D.1 DAG LEARNING METRICS
We evaluate the performance of each method using the following standard measures, each of which calculates the disparity between an estimated DAG and the ground truth DAG:
Structural Hamming Distance (SHD, lower is better): this is the most common metric in causal discovery, which counts the minimum number of edge additions, removals, and reversals to turn the estimated DAG into the true graph.
BIC score (higher is better): we also monitor the BIC score of the estimated DAG for all methods in the main text (despite the fact that some baselines do not optimize this score). For linear data, we use the BIC with equal variances and linear regression, while for nonlinear data, BIC with non-equal variance with Gaussian process regression is used, and BIC with logistic regression employed for binary data.
True Positive Rate (TPR, higher is better): this measures the ratio of correctly recovered edges over the true edges in the ground truth DAG.
False Discovery Rate (FDR, lower is better): this measures the proportion of incorrectly estimated edges over all estimated edges.
Precision, Recall, and Fj: this measures the binary classification performance by treating the binary adjacency matrix as set of individual binary classification tasks.
19

---

Published as a conference paper at ICLR 2025
D.2 SYNTHETIC CAUSAL GRAPHS
Our synthetic causal graphs involve one of the two well-known graph models:
Erdös-Rényi (ER, Erdös & Rényi, 1960): the edges in this type of graph are added inde- pendently with fixed probability. To generate a DAG of d nodes with an expected in-degree e, we first generate an undirected graph where edges are added with probability a(d then orient the edges using a random permutation over the list of nodes.
Scale-Free (SF, Barabási & Albert, 1999): these are graphs where the degree distribution follows power law, where few nodes have many connections while others have only few connections. To generate SF graphs with x de edges, we start with an empty graph then repeatedly grow it by attaching new nodes, each with & edges, that are preferentially attached to existing nodes.
Remark It is noteworthy that the majority of methods perform well only for graphs with up to = 4 (Zheng et al., 2018; Zhu et al., 2020; Ng et al., 2020; Yu et al., 2021; Wang et al., 2021; Bello et al., 2022; Yang et al., 2023b). Indeed, was noted (Bello et al., 2022) that ER4 and SF4 are the hardest settings. However, in this study, we have shown that our DrBO method is still very accurate on much denser ER8 and SF8 graphs
D.3 BNLEARN DATASETS
The BnLearn repository (Scutari, 2010) contains of Bayesian networks of varying sizes and complexities from different real-world domains. The chosen networks in our study include Alarm (Beinlich al., 1989). Asia (Lauritzen & Spiegelhalter, 1988), Cancer (Korb & Nicholson, 2010), Child (Spiegelhalter, 1992), and Earthquake (Korb & Nicholson, 2010). However, these datasets only describe the conditional probability distributions of discrete variables, while the baselines con- sidered in our method are mostly implemented for continuous data. Since the purpose of this ex- periment is to show that our method can correctly recover real structures, we use the networks from the BnLearn to generate continuous data. Specifically, we employ Sinaa--a SCMs as in Sec. 5.1.1 to generate synthetic data adhering to real causal networks, and each dataset contains only 1,000 observational samples.
D.4 IMPLEMENTATION AND PLATFORM
Implementations. We employ the following implementations for the considered baselines as fol- lows:
DAGMA (Bello et al., 2022): we use the official implementation provided by the authors at hthttps://gtthbb com/ wwe.kkseeellooooeeee /dagga.
COSMO (Massidda et al., 2024): we use the original implementation attached supple- mentary material at https:////// openreview wwtenreviewnnee//eeeeeeee net/ forum? id=Kw08LSUCW..
CORL (Wang et al., 2021): we employ the official implementation provided in the gCastle library (Zhang et al., 2021) https://gitie httppttththhooo..oooeooooeeeeete cwm/huawet -nnnaaeeiinhhh/// tzustworthyAI.
ALIAS (Duong et al., 2025): we reimplement their method by following the exact instruc- tions and librarics described.
GOLEM (Ng et al., 2020): we adopt the implementation provided by gCastle (Zhang et al., 2021) at https: : /tt////tttb///a com hwwwwi nww./trsstwortteettotoo
NOTEARS+TMPI (Zheng al., 2018; Zhang et al., 2022): we replace the DAG constraint in NOTEARS' implementation at https: /ttt:///githbb com/ wwnzhengggggguueeeeeeee hoainnngnnttrrrrrr with the TMPI constraint at htttts://ggiii com/ 22h.ng187/n
Platform. The majority of our experiments are conducted on a machine with Intel@ TM i9- 13900KF processor and NVIDIA RTX 4070 Ti GPU. The only exception is the case of CORL on
Data is publicly downloadable https: //ww.//bbo leann www..nccoooooooooommmmmmoooo.. wwwp/onnwwooooeoeeoooeooeoeeeeee tory
20

---

Published as a conference paper at ICLR 2025
large graphs (Figure 1(b)) which requires more than 16Gb of CUDA memory, and therefore these experiments are conducted on an NVIDIA A100 GPU with 40G of CUDA memory instead.
Additionally, in this study, we evaluate the performance of various methods with respect to both the number steps and runtime, addressing two independent questions: "How accurate can a method become given a fixed number of steps?" and "How accurate can it be within given runtime?". To ensure fair comparisons, we account for potential biases in measuring performance solely by the number of DAG evaluations. This is particularly important for methods like gradient-based approaches (e.g., DAGMA), which may require many steps but still exhibit low overall runtime.
To address this, we use runtime more equitable efficiency metric. Specifically, we a high number of steps for all methods (e.g., we use 'T' 800 iterations instead the default of only T= 5 for DAGMA on linear data) and disable early stopping if applicable, capture their progression over an extended period of time. We then truncate the tracking data, which contains performance metries and timestamp every step, either at a fixed number of steps or a specified runtime, as illustrated in Figure 1. This ensures that the results in the last column of Figure are not constrained by the number of steps. For instance, at the 5-minute mark in the last column of Figure la, DAGMA completes approximately 5 million steps compared to only 50,000 steps in the third column.
E HYPERPARAMETERS
We provide the specific set of hyperparameters for our DrBO method in Table 2. More details can be found in our published source code.
Table 2: Hyperparameters for DrBO. Unless specifically indicated, the default hyperparameters here are used for all experiments. More details can be found in our published source code. For nonlinear data with GPs, use GP regression with regularization noises have near-unit variances) and radial basis function (RBF) kemel, where the length scale is optimized over (I, 10²). FFor the Sachs dataset, we also use GP regression with the same kernel as above. In addition, as the noise variances are unknown, we set a really small value a 10-5 just to ensure positive definiteness of the covariance matrix, and following Wang et al. (2021); Duong et al. (2025), we also employ the median bandwidth heuristic for the kernel, by dividing the predictors by their median pairwise euclidean distance before GP regression is applied.
(because the additive
Experiment
Hyperparameter
Linear data
Nonlinear duta with GPs
Sachs data
Normalize data Scoring function
No Saic with lincar regression Linear pruning
No
Yes
Spic with GP regression" No pruning
Sprco with GP regressiont CIT pruning
Pruning method
Batch size B
64
DAG rank k
8
No. training steps rigrod
10
No. preliminary candidates
100,000
Optimizer
Adam
Learning rute
0.1
Replay buffer size Raplay
1,024
No. hidden units
64
Dropout rate
0.1
For ALIAS and CORL, apart from the recommended default hyperparameters, for fairness, we also use the same batch size of 64 as ours (except for CORL which requires the batch size to be at least the number of nodes, so for 100-node graphs we have to set the batch size to 100), as well as the same scoring function and pruning method as ours in each experiment. Regarding COSMO and DAGMA, we use the linear and nonlinear versions specific to each experiment and the hyperparameters are set recommended. We also use the same BIC score and pruner specific to each experiment to track the progress of all methods in Figure 1.
Regarding the number of evaluations, for all methods, we run more than needed then cut off at common thresholds, as mentioned in Appendix D.4. Remark Following the same line as in (Bello et al., 2022; Zheng et al., 2018) and many studies in this field, we intentionally avoid hyperparameter tuning. This is to prevent injecting spurious
21

---

Published as a conference paper at ICLR 2025
information of the dataset characteristics into the causal discovery results, which is against the idea of causality. Therefore, we use a fixed set hyperparameters for each method in each setting. The hyperparameters of baseline methods chosen recommended in the original manuscripts, while ours are chosen based on common practice and prior experience. Indeed, our ablation studies reveal that there are better hyperparameter choices that can further improve the performance of our method.
F ADDITIONAL CAUSAL DISCOVERY SETTINGS
F.1 DIFFERENT SAMPLE SIZES
In Figure 4, test our method with varying sample sizes linear-Gaussian data. The results indicate that our method is already very accurate with SHD s 0 at merely 500 samples.
Structural Hamming Distance
Cumulative runtime
Sample size n
100
80
500
Minutes
60
1000
40
5000
10000
20
2500
5000
7500
10000
0
2500
5000
7500
10000
DAG evaluations
DAG evaluations
Figure 4: Causal Discovery Performance with Varying Sample Sizes. We apply our DrBO
method linear-Gaussian data with 20ER4 graphs. Shaded areas represent 95% confidence interval over runs.
F.1.1 DIFFERENT GRAPH TYPES
We evaluate method different graph types, namely ER and SF, with varying densities, as shown in Table 3. Our method perfectly identifies the correct DAG in all cases for both graph models, even in the dense graphs ER8 and SF8.
22

---

Published as a conference paper at ICLR 2025
Table Causal Discovery Performance with Different Graph Types. The considered graph mod- els include Erdös-Rényi (Erdōs & Rényi, 1960) and Scale-Free (SF, Barabási & Albert, 1999), with expected in-degrees of and corresponding to the case of sparse and dense structures, respectively. We compare our DrBO method with ALIAS (Duong et al., 2025) and DAGMA (Bello et al., 2022) on
datasets with 20 variables and 10,000 samples. The numbers are mean std over
5 random datasets. For fairness, all methods are limited to 10,000 score evaluations.
Graph Expected In-degree
Method
SHD
FDR
TPR t
ALIAS
21.6 = 7.2
0.32 + 0.12
0.68 0.08
2
DAGMA
25.8 = 7.3
0.28 = 0.13
0.55 + 0.07
ER
DrBO (ours)
0.0 0.0
0.00 = 0.00
1.00 0.00
ALIAS
8.5
0.15 + 0.03
0.53 0.04
8
DAGMA
122.4 4.7
0.19 + 0.03
0.28 + 0.02
DrBO (ours)
0.0
0.00 = 0.00
1.00 = 0.00
ALIAS
13.6 2.7
0.18 + 0.09
0.71 + 0.02
DAGMA
23.4 9.8
0.25 + 0.19
0.53 0.13
SF
DrBO (ours)
0.0
0.00 0.00
1.00 0.00
ALIAS
68.8
0.34 0.08
0.52 0.06
8
DAGMA
94.4
0.49
0.2 0.06
DrBO (ours)
0.0 0.0
+ 0.00
1.00 0.00
F.1.2 DIFFERENT NOISE DISTRIBUTIONS
To examine the performance of our method compared with others under noise misspecification, in Table 4, we evaluate causal discovery performance various noise distributions. Specifically, we con- sider linear SCM F; = Σ jEpa, Mijiyt-EE where E; is drawn from one of the following distributions
Exponential noise: E, Exp (1), Vi I,... d, where I the scale parameter.
Gaussian noise: E; ~ N (0,1), = 1,.... d, where 0 is the mean and is the variance.
Gumbel noise: E, ~ Gumbel (0, (0,1),W) Vi 1,..., d, where is the location parameter and 1 is the scale parameter.
Laplace noise: 6; ~ Laplace (0,1),4 Vi = 1,...,d d, where O is the location parameter and the scale parameter.
Uniform noise: E; maximum value
(-1,11), 1,...,d, d, where the minimum value is -1 and
Remark 3. We note that our BIC score kept unchanged for different noises to show that our method can work well beyond the assumed Gaussian noise.
23

---

Published as a conference paper at ICLR 2025
Table Causal Discovery Performance with Different Noise Distributions. We consider 5 noise distributions: Exponential, Gaussian, Gumbel, Laplace, and Uniform. Our DrBO method is com- pared with ALIAS (Duong al., 2025) and DAGMA (Bello et al., 2022) on linear- -aaas-Ga datasets with 20ER4 graphs and 10,000 samples. The numbers are mean + std over random datasets. For
fairness, all methods are limited 20,000 score evaluations.
Noise Type
Method
SHD L
FDR
TPR
ALIAS
36.2 = 13.3
0.24 0.07
0.75 + 0.09 0.37 + 0.05
Exponential
DAGMA
56.6 6.8
0.22 0.10
DrBO (ours)
0.4 0.6
0.00 0.00
0.99 -: 0.01
ALIAS
34.2 11.4
0.23 0.08
0.77 0.05
Gaussian
DAGMA
57.2 7.6
0.25 0.09
0.39 + 0.09
DrBO (ours)
0.0 0.0
0.00 0.00
1.00 + 0.00
ALIAS
36.6 - 12.9
0.26 +
0.77 + 0.06
Gumbel
DAGMA
54.8 7.4
1.25+0000
0.44 0.08
BrBO (ours)
0.2 = 0.5
0.00 - 0.01
1.00 + 0.01
ALIAS
36.4
14.2
0.26 + 0.07
0.76 + 0.09
Laplace
DAGMA
56.0 8.8
0.25 + 0.08
0.40 + 0.09
BrBO (ours)
0.0 = 0.0
0.00 - 0.00
1.00 + 0.00
ALIAS
34.4
13.0
0.21 + 0.08
0.77 + 0.05
Uniform
DAGMA
59.0 7.3
0.26 + 0.08
0.39 + 0.09
DrBO (ours)
0.0
0.00 0.00
1.00 0.00
F.1.3 BGE SCORE FOR MARKOV EQUIVALENCE CLASS DISCOVERY
Here, show that method is not restricted to the BIC score and be extended to other scores as well. We consider the Bayesian Gaussian equivalent (BGe, Geiger & Heckerman, 1994; Heckerman et al., 1995) score for learning the Markov Equivalence Class (MEC) of DAGs in linear- Gaussian settings. The BGe score assigns equal scores to DAGs belonging to the same MEC, and be decomposed the of local scores follows:
d
SAGe (D,9) Σ LocalBGe; (pa
(17)
We refer readers to, for example, Kuipers et al. (2014), for the specific formula of the BGe score. To adapt our method with this score, we simply train each dropout network DropoutNN, to predict LocalBGe; from paf and combine them using Eq. (17). It can be seen that this adaptation does not involve changing any other component our method.
We report the results in Figure 5, where we compare DrBO with two popular baselines that are well-known to recover the MEC, namely PC (Spirtes et al., 2000) and GES (Chickering, 2002)." The results illustrate that our method can find the DAGs with the highest BGe scores, while the scores from PC and GES's estimations are well below those of the ground truths. As result, our recovered structures are more accurate than the baselines.
We use the gCastlc library (Zhang al., 2021) for their implementations, where hyperparameters are left default.
24

---

Published as a conference paper at ICLR 2025
DrBO (Ours)
PC
GES
-7.055
Extimated DAG's BGe
Estimated DAG's BGe
Estimated DAG's BGe
-7.4
-7.4
7.0
-7.6
1.00
0.10
-0.04
S1ID 0.61111
B1DD C 1.212.00
SLID C 19:26
-7.103
-7.100
-7.000
7.085
-7.0
Ground truth DAG's BGc
Ground truth DAG's
Ground truth DAG's BGe
Figure 5: BGe for Markov Equivalence Class Discovery. We compare the BGe score of ground truth DAGs and the estimations from DrBO with two popular baselines PC (Spirtes et al., 2000) and GES (Chickering, 2002). Each point corresponds to one of 50 random datasets with linear Gaussian data on ER graphs of nodes and 5 edges average. The Pearson correlation coefficient r between the scores the estimated and ground truth DAGs are included. addition, we also report the SHD-C metric, which measures the structural distance between MECs.
F.1.4 DISCRETE DATA
this section, we show that our method is not limited to continuous data either. To demonstrate, we consider binary data with logistic causal models:
~ Bernoulli (ft (xxa))
We adapt method to this situation by simply changing the BIC score to take into account logistic models. More particularly, we use the BIC score for logistic data as in Eq. (16), where the maximum log-likelihood can be decomposed into
= LocalMLL;
where LocalMLL, = E (z.6) (x() + +11-
pa. ))) with f being x0.)) In (1 (1- (x())
pag
maxinum-likelihood estimator found via logistic regression.
Then, we employ our local surrogate models DropoutNN, to model LocalMLL,, and the acquisi- tion function values are obtained by summing the local score samples. Again, this adaptation only involves changing the scoring function and does not modify any other component of our BO frame- work.
We report the causal discovery performance under binary data in Figure 6, where also compare DrBO with the state-of-th of-te-of-the-rr art DAGMA, which is now set to use the logistic loss supported. be seen that our method also performs well for binary data, where it consistently obtains the highest scoring graphs, resulting in a low structural error. Meanwhile, DAGMA usually finds only sub-optimal solutions which lead to high structural errors.
25

---

Published as a conference paper at ICLR 2025
DrBO (Ours)
DAGMA
X10
X30
-5.4
-5.50
Estimated DAG's BIC
-5.6
Estimated DAG's BIC
-5.75
-5.8
6.0
-6.00
6.2
-6.25
66.4
-6.50
-6.6
1.00
-6.75
0.92
-6.8
SHD-C 0.311.1
SHD-C 2.4
6.55
-6.0
-5.5
6.5
6.00
-5.5
X10-1
x10-
Ground truth DAG's BIC
Ground truth DAG's BIC
Figure 6: Causal Discovery performance on Binary Data. We compare our DrBO method using the BIC score for the logistic model with DAGMA (Bello et al., 2022) using the logistic loss. Each point corresponds to one of 50 random datasets with logistic data on ER graphs of 5 nodes and edges on average. The Pearson correlation coefficient between the scores of the estimated and ground truth DAGs are included. In addition, we also report the SHD-C metric, which measures the structural distance between MECs.
F.1.5 STANDARDIZED DATA
As discussed in Reisach al. (2021), marginal variances may contain crucial information about the causal ordering among the causal variables, and thus revealing the causal DAG by simply sorting the variables by increasing variances. For this reason, in this section, we investigate the causal discovery performance of the proposed method in comparison with the baselines under uninformative marginal variances, by standardizing the observed data to have zero mean and unit variance per dimension, before feeding it to causal discovery methods.
Linear-Gaussian data. Since the noise variances are non-equal after standardization, we employ the non-equal variance versions of the methods that support it, including ALIAS, CORL, GOLEM, and DrBO. In addition, as data is standardized, the usual threshold 0.3 for pruning is no longer appropriate because significant edge weights may be rescaled much smaller values after standardization, so this experiment, we increase the sample size to 100,000 to reduce weight estimation variance, and lower the pruning threshold 0.01. The results presented in Figure 7(a) confirm that our method is still robust for standardized data. Overall, while all methods obtain a non- zero SHD due the difficulty of standardized data, which render the equal- variance SCM unidentifiable, our method still outperforms other baselines significantly, where we achieve an SHD*33, while the second-best SHD nearly 20, highlighting DrBO's improved effectiveness over existing approaches in this intricate scenario. Additionally, even though the same pruning threshold is used for all methods, our method barely predicts any extra edge, while other baselines still suffer from high false discovery rate. Moreover, our method does not predict too many reverse edges, as opposed to most methods, showing that while data standardization negatively impacts causal discovery performance to some extent, the effect our method is minimal.
Nonlinear data with GPs. Figure 7(b) presents the results for nonlinear data with GPs with stan- dardization, showing that our method can achieve very low SHD and surpasses other methods considerably. This result is similar to Figure 1(c), where the same datasets employed are not stan- dardized, indicating that the performance of our method is not affected by data standardization in this case, which could be potentially thanks to the fact that nonlinear ANMs remain identifiable after standardization.
26

---

Published as a conference paper at ICLR 2025
SHD 4
Missing
Extra +
Reverse +
Ireo Alias
Buasss
Noneass Couls
centio
msoo (Dwn)
AlMA
coal.
couso
prto town) Aldds
crat.
crsuup
Alss
cas
Dassss
NOTEARS delt
crstio
NUTRASS-
(a) Linear- Gaussian Data (DAGs with 10 nodes and 20 edges). For fairness, all metrics are calculated at 20,000 evaluations for all methods.
SHD
Missing
Extra 1
Reverse 4
30
Dee (Oue
cons
NOTROss raay
Deto rown)
AuAs
camo
AYTARSS raay
Addss
copt
paraa
NOTEARS man
aseo (owe)
recass
ersso
NoT may
Urit
(b) Nonlinear Data with GPs (DAGs with 10 nodes and s: edges). For faimess,, all metrics are
calculated at 20,000 evaluations for all methods.
Figure 7: Causal Discovery Performance on Standardized Data. Performance metrics are Struc- tural Hamming Distance (SHD), number of Missing, Extra, and Reverse edges. Lower values are more preferable. Error bars indicate 95% confidence intervals over 5 simulations.
F.1.6 WHY LOW-RANK DAG REPRESENTATION TENDS TO PERFORM BETTER
In Figure 3(a), we have shown that the rank our DAG representation plays an important role in the sample-effciency of our method, where lower ranks clearly enable reaching the accurate DAGs carlier than their high-rank counterparts. This is because it is much more challenging to search in very high-dimensional space compared to lower-dimensional one. Specifically, for higher ranks, the search space much larger and sparser than the low-rank ones, and due to the curse of dimensionality, sampling the same number of random DAG candidates (step in Algorithm 1) in the higher-rank search spaces tends to lead to fewer unique candidates compared with a lower-rank one, reducing the chance meet the higher-scoring DAGS earlier.
To empirically verify this, we calculate the number of unique DAGs among 1,000 random 30-node DAGS generated with different ranks in Table 5. It can be seen that, typically, the lower the rank, the more unique DAGs we can pre-examine for exploration. For k = 2. almost every DAG among 1,000 generated DAGs is unique, and thus our method can reach SHD very quickly Figure 3(a), whereas the full-rank representation is dimees-di and can only generate fewer than half the unique DAGs. addition, for generate fewer than 10% of unique DAGs, explaining why this representation is the least sample- efficient one Figure 3(a).
32 d, where the dimensionality highest, we can only
Table 5: Effect of DAG Rank on Exploration Diversity. We generate 1,000 DAGs with 30 nodes using g := T (z) E [-1, (1+k)) with different The numbers are mean + std over 10 simulations.
Rank Eq. (4)
Number of dimensions
Number unique 30-node DAGs over 1,000 random DAGs
2
90
926.7 = 7.0
150
779.2 12.7
8
270
403.5 12.3
12
390
332.1 = 10.8
32
990
Full rank (Vec2DAG, Duong et al., 2025)
465
421.9 = 13.8
27

---

Published as a conference paper at ICLR 2025
F.1.7 LARGE-SCALL NONLINEAR DATA
data with nonlinearity, following Zhang et al. (2022), we evaluate our method 50ER2 and 100ER1 graphs with nonlinear SCM == B cos (x) + E, where the weights B are sampled uniformly in |-2 -0.5] U [0.5, 2] and & ~ N (0, This model is identifiable according to Bühlmann et al. (2014). To calculate the BIC score, we employ linear regression with cosine features. The results reported in Figure 8 demonstrate that our method is also competitive on large-scale nonlinear data, where it can outperform the baseline in both causal discovery perfor- mance and runtime by visible margin.
For
1 - FDR
BIC score (higher is better)
SHD (lower better)
SHD (lower is better)
100
TPRR
Fi
50
0.0
-6.4
4.2
Precision
Recall
-3
1.6
2.2
0
10000
20000
10000 20000
0
10
DAG evaluations
DAG evaluations
Minutes
DrBO (Ours)
ALIAS
(minutes) to
Method
SHD=20
SHD>110
SHD=5
SHDD=
ALIAS
0.6
5.1
8.3
9.1
DrBO (ours)
1.5
21=14
452
6.2
(a) 50 nodes and 100 edges. For fairness, summary performance metrics (circular plot)
calculated at 20,000 evaluations for both methods.
1 -
FDR
BIC score (higher is better)
SHD (lower better)
SHD (lower is better)
76
75
TPR
Fi
50
Precision
Recall
-4-
10000
20000
10000 20000
20
DAG evaluations
DAG evaluations
Minutes
DrBO (Ours)
ALLAS
(minutes)
Method
SHD=20
SHD-10
SHD=5
SHD=2
ALIAS
21.0
33.7
DrBO tours)
2.9 = 0.6
5.0 5.0=13 13
8.0 = 3.3
12.2
(b) 100 nodes and ~J 100 edges. For fairness, summary performance metrics (circular plot) are calculated at 20,000 evaluations for all methods.
Figure 8: DAG learning results on large-scale nonlinear data.
F.2 ABLATION EXPERIMENTS
F.2.1 EFFECT OF DROPOUT NETWORKS
In Figure 3(b), we compare our dropout networks as the surrogate model with exact GPs and approx- imate GPs. Approximate GPs learn a set of pseudo data points called inducing points and conduct inference via these points instead of the real data (Hensman et al., 2015). Here we number 100 inducing points for Approximate GPs, to see they can scale well with few inducing points. Due to the limited scalability and intensive memory requirement for GPs, can only use C' = 10,000 preliminary candidates on which we sample from the posteriors, while our default hyperparameter C 100,000 using dropout networks.
small
28

---

Published as a conference paper at ICLR 2025
F.2.2 EFFECT OF CONTINUAL TRAINING
In Figure 3(d), we compare the continual training approach with fully retraining using all data. For continual learning, meaning for each BO iteration, we perform 10 gradients update, each update calculated from 1,024 64 1,088 datapoints. To ensure faireess, also rlgrads 10 epochs and mini-batch size of 1,088 for the full retraining approach.
default hyperparameters repllyy 1,024, B 64, and Ngnds 10,
F.2.3 EFFECT OF EVALUATION BATCH SIZE (B)
In Figure 9, we show the influence of the evaluation batch size B onto the performance and scala- bility of our method. Overall, is clear that smaller batch sizes lead better SHD but much worse runtime since the surrogate model is updated more frequently, and vice versa. However, B 64 seems to achieve balance, where it enables SHD As 0 and lower runtime than smaller batch sizes.
Structural Hamming Distance
Minutes
Evaluation batch size (B)
16
64
256
512
Figure 9: Effect of Evaluation Batch Size B. We evaluate our method on Emear-Gaussin data with 20ER4 graphs and 1,000 observations. Error bars indicate 95% confidence intervals over runs. The number of evaluations is limited to 20,000.
F.2.4 EFFECT OF NUMBER OF PRELIMINARY CANDIDATES (C)
We show the variation our DrBO's performance w.r.t different numbers of preliminary candidates C' in Figure 10, showing that the best performance and runtime can be achieved at 10,000 candidates, and even with 10x more candidates, the runtime of our method increases only slightly.
Structural Hamming Distance
Minutes
No. preliminary candidates (C)
5000
10000
20000
100000
10
8.0
8.5
9.0
9.5
Figure 10: Effect of Number of Preliminary Candidates C. We evaluate our method on linear- Gaussian data with 20ER4 graphs and 1,000 observations. Error bars indicate 95% confidence intervals over runs. The number of evaluations is limited to 20,000.
29

---

Published as a conference paper at ICLR 2025
F.2.5 EFFECT OF NUMBER OF TRAINING STEPS PER BO ITERATION (NORADS)
We study the effect of the number of gradient steps in each BO iteration (Agradss) in Figure 11. In general, small values may lead underfitting and large values may be prone overfitting, so values in the middle are better for this hyperparameter.
Structural Hamming Distance
Minutes
No. training steps (rigras))
[e- 5
10-
[+-
20
HH
50
2
10
15
20
Figure 11: Effect of Number of Training Steps per BO Iteration Rigrads We evaluate our method on Amner-GGassiaal data with 20ER4 graphs and 1,000 observations. Error bars indicate 95% confi- dence intervals over 5 runs. The number of evaluations is limited to 20,000.
F.2.6 EFFECT OF REPLAY BUFFER SIZE (PREPLAY)
In Figure 12, we show that higher values for the replay buffer size significantly reduces. SHD but does not considerately influence the runtime.
Structural Hamming Distance
Minutes
No. retraining samples (Mreplay)
256
512
HH
2048
4096
2
8.2
8.4
8.6
Figure 12: Effect of Replay Buffer Size nreplay We evaluate method Ghuas-Ga data with 20ER4 graphs and 1,000 observations. Error bars indicate 95% confidence intervals over runs. The number of evaluations is limited to 20,000.
F.2.7 EFFECT OF LEARNING RATE
Figure 13 depicts that the learning rate has weak effect the performance and scalability of our method, where any value below can achieve the same level of SHD and runtime. The SHD only becomes large for high learning rate of 1.
30

---

Published as a conference paper at ICLR 2025
Structural Hamming Distance
Minutes
0.0001
Learning rate
0.001
0.01
1.0
8.46
6.47
8.48
8.49
8.50
Figure 13: Effect Learning Rate. evaluate our method on linear-Gaussian data with 20ER4 graphs and 1,000 observations. Error bars indicate 95% confidence intervals over 5 runs. The number evaluations is limited to 20,000.
F.2.8 EFFECT OF NUMBER HIDDEN UNITS
We present in Figure 14 that the number of hidden units in our dropout networks also has a visible effect method's performance, but not much the runtime. Specifically, moderate value of 32 achieves vanishing SHD with the equivalent runtime as others. Meanwhile, to many hidden units may challenge the training process performance may drop.
Structural Hamming Distance
Minutes
H
No. hidden units (h)
H
64
128 |-|1
512
[o]
8.2
8.4
Figure 14: Effeet of Number of Hidden Units h. We evaluate our method on data with 20ER4 graphs and 1,000 observations. Error bars indicate 95% confidence intervals over runs. The number of evaluations is limited to 20,000.
F.2.9 EFFECT OF NUMBER OF DROPOUT RATE
Figure suggests that the performance of our method improves with higher dropout rates, while the runtime does not vary significantly.
31

---

Published as a conference paper at ICLR 2025
Structural Hamming Distance
Minutes
0.01
Dropout rate (p)
0.05
0.1
0.3
8.42
8.44
6.46
8.48
8.50
Figure 15: Effect of Dropout Rate p. We evaluate our method on data with 20ER4 graphs and 1,000 observations. Error bars indicate 95% confidence intervals over 5 runs. The number evaluations is limited to 20,000.
G ADDITIONAL BASELINES
Apart from the score-based competitors compared so far, here we also consider additional baselines for comparison, such as conventional methods PC (Spirtes al., 2000) and GES (Chickering, 2002). In addition, ordering-bant based methods have also been gaining popularity (Reisach et al., 2021; Rolland al., 2022; Sanchez al., 2023; Montagna et al., 2023). These methods sidestep the acyelicity issue by first learning topological ordering of the causal DAG, then pruning the fully-connected DAG induced by the ordering to obtain the causal structure.
In Table 6, we provide comprehensive comparison between method and conventional methods PC and GES, as well as popular ordering baser method, including sortnregress (Reisach et al., 2021), SCORE (Rolland et al., 2022), and DiffAN (Sanchez et al., 2023), on linear, nonlinear, as well as real data. The empirical evaluations indicate that our method is also able to surpass these methods in all metries and settings.
32

---

Published as a conference paper at ICLR 2025
Table Comparison with Ordering-based Methods. We compare our method with conventional methods PC (Spirtes et al., 2000) and GES (Chickering, 2002), as well as popular ordering- based algorithms sortnregress (Reisach al., 2021), SCORE (Rolland al., 2022), and DiffAN (Sanchez et al., 2023). The orderings produced by those methods are first transformed into respective fully connected DAGs. Then, fairness, all DAGS, including ours, are pruned using the weight matrix thresholded at 0.3 for linear data, CAM pruning (Bühlmann et al., 2014) for nonlinear data, and KCIT (Zhang et al., 201 for the Sachs dataset. The figures are mean = over datasets, except for the Sachs dataset. Note that DrBO's result in Figure ((c) with SHD 0.4 is obtained without pruning, while using CAM pruning leads to a few missing edges as shown below.
Dataset
Method
SHD
FDR
TPR
F1 1
PC
10.8
5.1
0.29 - 0.1
0.49 +49
0.56 = 0.2
GES
11.6 6.4
0.40 =.2
0.61
0.58 = 0.2
Linear Data (1OER2 graphs. 1, 1,000 samples)
sortnregness
2.6
0.10
0.94
0.92 0.1
SCORE
0.0
0.0
1.00 0.0
1.00 0.0
DifFAN
16.6
0.54
0.4
0.48 + 0.1
DrBO (ours)
0.0 0.0
0.0 0.0
1.00 0.0
1.00 + 0.0
PC
227.0 + 8.6
0.44 0.0
0.10 +.0
0.17 + 0.0
GES
Failed to halt
Linear Data (30ER8 graphs,
sortnregress
101.4 21.8
0.0 0.76
0.80 + 0.1
0.78 = 0.0
SCORE
247.0 + 11.8
0.05 0.0
0.08 + 0.1
1, 1,000 samples)
DiffAN
236.0 + 3.7
0.58 0.1
0.16
0.23 + 0.1
DrBO (ours)
1.6 + 1.5
0.00 + 0.0
0.99 0.0
1.0 0.0
PC
32.2
0.50
0.21 0.1
0.29 - 0,1
GES
29.6
0.2
0.30
0.38 - 0.1
Nonlinear Data
(10ER4 graphs,
sortnregress
3.3
0.1
0.37 0.1
0.45 = 0.1
SCORE
9.4 4.6
0.11
0.83 0.1
0.86 - 0.1
1,000 samples)
DiffAN
18.6
2.663388 3.8
0.34
0.61 + 0.1
0.63 0.1
DrBO (ours)
4.2 1.3
0.00 0.0
0.90 = 0.0
0.95 0.0
PC
11
0.25
0.35
0.39
Sachs et al. (2005), (11 nodes, 17 edges, 853
GES
11
0.25
0.35
0.39
sortnregress
13
0.44
0.29
0.38
SCORE
12
0.33
0.35
0.46
samples)
DiffAN
16
0.75
0.12
0.16
DrBO (ours)
9
0.11
0.47
0.62
H
TIME COMPARISONS
Table Runtime comparison on 30ER8 Gaeeaa- data (experiment in Figure 1(a)).
Max number of evaluations
Method
10,000
20.000
50,000
SHD
Runtime (mins)
SHD
Runtime (mins)
SHD
Runtime (mins)
ALIAS
208.8 + 13.1
0.5
165.0 18.0
105.2
2.4 = 0.0
CORL
1.1
162.0 13.4
145.4 29.3
5.2 0.1
COSMO DAGMA
213.6 + 17.4
0.6 0.0
203.4 + 21.0
1.1 +
195.8 10.8
2.8 + 0.0
222.2 4: 11.2
0.0 0.0
218.8 4: 10.6
0.0 0.0
181.6 22.4
0.0 -:: 0.0
DrBO (ours)
2.0 + 1.2
4.6 0.0
2.0 + 1.4
9.2 0.0
1.6 1.5
22.9 + 0.1
33

---

Published as a conference paper at ICLR 2025
Table 8: Runtime comparison 100ER2 Finear-Gaussie data (experiment in Figure 1(b)..
Max number of evaluations
Method
10,000
20,000
50,000
SHD
Runtime (mins)
SHD
Runtime (mins)
SHD
Runtime (mins)
ALIAS
230.4 +: 32.7
0.3
136.8 +: 30.2
5.5 4.
32.0
16.7
13.5 0.6
CORL
148.0 = 42.9
11.1
120.2 + 18.1
22.0 0.2
81.0 = 19.9
54.2
COSMO
111.0 = 10.9
0.8
111.8 12.2
1.5
112.6 = 13.4
3.7
DAGMA
124.8 17.4
0.0
93.6 +
6,6 4.0
0.3 0.3:00
DrBO (ours)
29.2 = 16.7
12.7
3.4
25.4 +
1.4 1.1
62.7 2..7:...5
Table 9: Runtime comparison on 10ER4 nonlinear data with Gaussian processes (experiment in
Figure I(c)).
Max number of evaluations
Method
1,000
2,000
20,000
SHD
Runtime (mins)
SHD
Runtime (mins)
SHD
Runtime (mins)
ALIAS
15.8 - 4.4
4.0 + 0.2
12.6 + 4,2
5.8 = 0.2
4.0 2.3
8.0 0.3
CORL
1.4
20.9 00.933 3.0
9.2 + 2.3
26.5 + 3.1
8.4
29.9 2.4
COSMO
2.6
0.3 0.3
33.0 3.3
0.4
25.6 4.1
1.8 0.1
DAGMA
40.4
0.0
38.4 2.7
0.1
34.2 3.5
0.8 0.1
DrHO (ours)
3.2+1.6
3.7 0.8
3.9 0.9
0.4 + 0.5
6.1 1.0
34


Reviews:
[{"review_id": "xgUlQatpON", "reviewer": "Reviewer_gMGJ", "summary": "The authors propose DrBO (DAG recovery via Bayesian Optimization)—a novel DAG learning framework leveraging Bayesian optimization (BO) to find high-scoring DAGs.  To address the scalability issues of conventional BO in DAG learning,  the authors replace Gaussian Processes commonly employed in BO with dropout neural networks, trained in a continual manner. DrBO is computationally efficient and can find the accurate DAG in fewer trials and less time than existing state-of-the-art methods. This is demonstrated through an extensive set of empirical evaluations on many challenging settings with both synthetic and real data.", "questions": "Could the authors give more details on how to ensure the binary adjacent matrix is a DAG in the optimization steps?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Learning DAG from data using BO is novel and interesting.  The authors overcome the scalability issue of conventional BO by leveraging dropout in neural networks. Experimental results show that the proposed method is effective and can achieve improved results. The paper was written with technical details.", "weaknesses": "N/A", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "RWL0M8yzfO", "reviewer": "Reviewer_UyS6", "summary": "This paper proposes an efficient causal discovery algorithm with Bayesian optimization. In particular, the authors consider a variant of Vec2DAG (Duong et al., 2024) for the DAG constraint, and use dropout neural networks and continual training scheme to optimize the adjacency matrix. Experiments show the efficiency and effectiveness of their method.", "questions": "- In Eq.(4), is the matrix $R$ strictly upper-triangular?\n- In Figure 1(b), did the authors still use 1000 samples for the large-graph experiments? $n=1000$ for the graph with 100 nodes?\n- In Figure 3(a), why in general smaller $k$ could obtain higher performances, compared with the full-rank cases?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- This paper is written clearly, with clear and detailed descriptions of their method and experiments.\n\n- They performed extensive experiments for validation.", "weaknesses": "- While there exist some causal discovery algorithms with Bayesian optimization, it seems not proper to state “To our knowledge, this is the first score-based causal discovery method based on BO ”. I think it should be corrected.\n\n- Throughout the paper, from the experiments, it is demonstrated that the proposed method can give better performances in both accuracy, sample-efficiency, and scalability, compared with other SOTA baselines. Generally, such a great method needs more assumptions or conditions to be satisfied. But intuitively, I cannot find these assumptions or conditions. Did this method have some other implied assumptions or conditions (like more hyperparameters)? \n\n- In experiments, it would be good to compare some Bayesian causal discovery methods (Deleu et al., 2022: Tranet al., 2023; Annadani et al., 2023), since they are all causal discovery methods. Or explain the reasons why not comparing with them.\n\n- The code is not available for reproduction.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "ZMHkc1xO7f", "reviewer": "Reviewer_VYj5", "summary": "This paper develops a Bayesian optimization method for score-based causal discovery. Several design choices are adopted, which include (1) developing a low-rank DAG representation, (2) replacing Gaussian process in conventional Bayesian optimization with dropout neural networks, (3) learning the DAG score indirectly via node-wise local scores, and (4) training in a continual way. Empirical studies are provided.", "questions": "- Is there a reason why the paper considers only identifiable models? That is, why general linear Gaussian model cannot be learned by the BIC-NV score?\n- Although BIC-NV is given, it seems that the experiments focus on equal variances. I would suggest adding experiments for different variances as well.\n- Baselines: For linear case, GOLEM may also be included. Also, adding the results for more conventional search methods, such as GES/FGES, may also be helpful.\n- Why does DAGMA-MLP performs so poorly for nonlinear data? The TPR is close to 0. If the reason is due to instability in optimization, the paper may consider adding NOTEARS-MLP that may be more stable.\n- For Section 5.2, specifically Sachs data, did the paper use linear or nonlinear version of the method?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "- The paper is well written and easy to follow.\n- Developing effective search procedure for score-based causal discovery is an interesting and important topic. The proposed method adopts various design choices and is practical.\n- The search method is reasonable.\n- The empirical studies demonstrate that the proposed method considerably outperforms existing methods.", "weaknesses": "- Some of the baselines considered are not adequate.\n- Some of the results may seem too good to be true. For example, achieving a SHD of 1.6 with only 1000 samples across 30 nodes and 240 edges seems highly challenging due to finite sample error. This concern is especially relevant when dealing with nonlinear data. (I look forward to the authors' clarification/explanation on this, and please correct me if I misunderstood anything.)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "VQIzwzikwv", "reviewer": "Reviewer_ZehJ", "summary": "This paper introduces DrBO, a Bayesian Optimization-based framework for efficient and accurate DAG learning from observational data. By leveraging dropout neural networks instead of Gaussian Processes, DrBO addresses scalability issues while integrating uncertainty in score estimation. Empirical results demonstrate DrBO's improved efficiency and accuracy over existing methods across synthetic and real datasets.", "questions": "1. To my knowledge, the CBO series of papers assumes that the DAG structure is known and primarily focuses on optimizing policies with this prior knowledge. Therefore, these papers may not be directly relevant to the active causal discovery literature. Please revise this in the introduction to reflect the distinction.\n\n2. The synthetic dataset, as first used in NOTEARS, is known to be relatively easy to learn, making pursuit of very low SHD scores less meaningful in recent research, especially the very simple linear gaussian case. How does your method perform on this dataset after standardization, as described in the paper *Beware of the Simulated DAG*?\n\n3. The proposed method appears to be limited to ANMs in causal discovery, which restricts the scope of the paper. It may be more accurate to frame the task as DAG structure learning or Bayesian structure learning rather than causal discovery.\n\n4. In Section 4.1, the authors mention that their method incorporates a low-rank adaptation of Vec2DAG. Does this imply an assumption about the data’s structure, as discussed in *On Low Rank Directed Acyclic Graphs and Causal Structure Learning*? Additionally, what would occur if \\( k < d \\)?\n\n5. Replacing the Gaussian Process in Bayesian Optimization with Dropout is not uncommon, so it may not warrant being highlighted as a novel contribution in this paper.\n\n6. While many prior works employ CAM as a pruning method, I believe this approach may lack justification here. Why would score-based search methods, including this paper, attempt to prune under nonlinear conditions? It's unusual for newly proposed methods to rely on post-processing from an older method.\n\n7. Please compare this baseline method, *Truncated Matrix Power Iteration for Differentiable DAG Learning*, to your approach.\n\n8. In Figure 1, several methods being compared fail to converge. For the tabulated results, have you ensured that all comparison methods have converged? Additionally, consider moving the running time details from the appendix to the main content, as the high time complexity is a notable limitation of the proposed method.\n\n9. Please provide results for nonlinear functions with datasets of 50 and 100 nodes, detailing both performance metrics and running time.", "limitations": "", "rating": 6, "confidence": 5, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- Using BO in causal discovery is interesting and novel. \n- The paper is very well written.", "weaknesses": "see questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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:

Decision: Accept (Poster)