Sample Efficient Bayesian Learning of Causal Graphs from Interventions

Authors: zihan  zhou
muhammad qasim elahi
murat  kocaoglu
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 21:36:56

Abstract: Causal discovery is a fundamental problem with applications spanning various areas in science and engineering. It is well understood that solely using observational data, one can only orient the causal graph up to its Markov equivalence class, necessitating interventional data to learn the complete causal graph. Most works in the literature design causal discovery policies with perfect interventions, i.e., they have access to infinite interventional samples. This study considers a Bayesian approach for learning causal graphs with limited interventional samples, mirroring real-world scenarios where such samples are usually costly to obtain. By leveraging the recent result of Wienöbst et al. [2023] on uniform DAG sampling in polynomial time, we can efficiently enumerate all the cut configurations and their corresponding interventional distributions of a target set, and further track their posteriors. Given any number of interventional samples, our proposed algorithm randomly intervenes on a set of target vertices that cut all the edges in the graph and returns a causal graph according to the posterior of each target set. When the number of interventional samples is large enough, we show theoretically that our proposed algorithm will return the true causal graph with high probability. We compare our algorithm against various baseline methods on simulated datasets, demonstrating its superior accuracy measured by the structural Hamming distance between the learned DAG and the ground truth. Additionally, we present a case study showing how this algorithm could be modified to answer more general causal questions without learning the whole graph. As an example, we illustrate that our method can be used to estimate the causal effect of a variable that cannot be intervened.

Main paper: 
 introduction: 
 Causal discovery refers to learning the underlying causal graph of a data generating process using a combination of observational and interventional data. This is a fundamental problem with applications across various areas including economics, genomics, meteorology, could computing, and etc. [Pearl, 2009 , Hoover, 1990 , King et al., 2004 , Ikram et al., 2022 , Runge et al., 2019] . In recent decades, the technological progress has facilitated the accumulation of vast quantities of observational data, i.e., data collected without perturbing the underling causal mechanisms. However, with only observational data, one can only recover the causal graph up to its Markov equivalence class (MEC) [Verma, 1991 , Andersson et al., 1997] . In general, interventional data, i.e., data collected after a perturbation in Additionally, the performance of our algorithm is stable across different settings (order and density of the random DAGs).• We discuss how we can modify our algorithm to answer more general causal questions without learning the whole graph using DAG sampling from the MEC. We show how we can estimate the posterior of a set to be the backdoor adjustment set of a causal query given the essential graph and some interventional data as an example. We further show through simulated experiments that we can modify our algorithm to estimate causal effects of variables that cannot be intervened.Outline of the paper: In Section 2, we summarize the related works. In Section 3, we state the preliminary notations and background knowledge. In Section 4, we describe how we setup the initialization steps of the main algorithm. In Section 5, we describe the algorithm with details and analysis on the convergence to the ground truth. In Section 6, we discuss how our algorithm can be modified to answer more general causal questions with an example. In Section 7, we present the experiment results on simulated datasets and compare with baselines. In Section 8, we conclude with discussion of limitations and future extensions. 
related works: 
 We can roughly divide the causal discovery methods into 3 categories: non-adaptive, adaptive, and Bayesian. Non-adaptive methods design intervention targets for each experiment before retrieving interventional samples. Since no information is shared between experiments, they can be conducted in parallel. Eberhardt [2007] discussed non-adaptive methods as fixed searching strategies and they show that in the worst-case scenaris O(log n) size n 2 interventions are necessary to orient the whole causal graph, regardless of the adaptivity of the algorithm. Hyttinen et al. [2013] show the connection between orienting causal graph via intervention and the concept of separating system from combinatorics. Ghassami et al. [2018] studied the problem of orienting the maximum expected number of edges with a budgeted number of size 1 interventions. Kocaoglu et al. [2017] , Lindgren et al. [2018] studied the problem of non-adaptive causal graph learning while minimizing intervention costs when each vertex is assigned with a different intervention cost. Many of the early works in causal discovery are non-adaptive, and therefore they construct the theoretical basis of causal discovery algorithm. Separating system is also used in our proposed method.Adaptive methods make the decision of next intervention based on the results of previous experiments. Less interventions are usually required than non-adaptive counterparts. Eberhardt [2007] compares adaptive strategies with non-adaptive ones and propose a heuristic algorithm. He and Geng [2008] proposed MinMaxMEC and MaxEntropy algorithms that select the next intervention target that minimizes the maximum I-MEC size or maximizes the intervention entropy respectively. Hauser and Bühlmann [2014] proposed OptSingle which calculates the target that minimizes the maximum number of unoriented edges in the interventional graph. Squires et al. [2020] show that any algorithm would take at least ⌈ ω 2 ⌉ size 1 interventions to recover the whole graph and proposed a two-stage active learning algorithm. In the first stage, it finds the directed clique tree representation of the unoriented graph and orient the residuals separately in the second stage. Choo et al. [2022] show that a size 1 invervention set can orient the whole graph if and only if it cuts all the covered edges in the ground truth DAG, and propose a learning algorithm that intervenes on a 1/2-separator in each experiment to iteratively decrease the size of the currently unoriented components.The aforementioned works all assume the existence of an infinite intervention oracle. The Bayesian approach is usually used when this assumption is untenable. Heckerman et al. [1997] investigated learning Structural Causal Models (SCMs) under the assumption that each local likelihood is a combination of multinomial distributions and all variables are discrete. Similarly, Buntine [1994] discussed simple linear local likelihoods involving both discrete and continuous variables. Tong and Koller [2001] proposed an active learning algorithm for Bayesian network learning. Previous studies typically employed Markov Chain Monte Carlo (MCMC) for sampling Directed Acyclic Graphs (DAGs) to estimate posterior distributions over DAGs and function parameters, which grow exponentially with the size of the model. Acharya et al. [2018] studied the problem of distinguishing two causal models on the same graph with limited interventional samples. Nishikawa-Toomey et al.[2022] uses variational inference and GFlowNets [Bengio et al., 2023] to learn and parameterize the joint posterior distribution over DAGs and mechanisms while assuming linear model and equal noise variances. Recent researches [Lorch et al., 2021 , Charpentier et al., 2022 , Cundy et al., 2021 , Toth et al., 2022 , Hägele et al., 2023] started leveraging gradient information for more efficient inference, while they suffer from poor inference quality. These studies all make specific parametric assumptions regarding functions or noises in the SCM, which lead to incorrect outcomes when these assumptions are violated. Kuipers et al. [2022] combines constraint-based methods and MCMC methods to efficiently search DAGs, but no guarantee of convergence is provided. To circumvent such issues, Greenewald et al. [2019] directly update the posterior of a specific subgraph containing the root when the graph is a tree without any constraint on the SCM. However in practice, causal graphs are usually more complicated than trees, which motivates our work. 
preliminaries: 
 A causal graph D = (V, E) is a directed acyclic graph (DAG) where the vertex set V represents a group of random variables. In the context of D, a directed edge (X, Y ) ∈ E from variable X to variable Y (denoted as X → Y ) signifies that X serves as an immediate parent of Y . We denote the parent set of variable Y by pa(Y ). The cut at vertex set X, denoted as E[X, V \ X], constitutes the set of edges between X and some vertex in the set V \ X. According to the Markov assumption, the joint distribution can be decomposed as P (v) = n i=1 P (v i |pa(X i )). A causal graph entails specific conditional independence (CI) relationships among variables via d-separation statements. The d-separation serves as a criterion to determine whether a set of variables X is independent of another set Y, given a third set Z. This approach links dependence to connectedness, meaning the existence of a connecting path, while it associates independence with the absence of connection or separation. A set of DAGs is deemed Markov equivalent when they entail identical conditional independence (CI) statements. All Markov equivalent DAGs must have the same adjacencies and unshielded colliders (or v-structures). An unshielded collider is a graph structure of the form A → B ← C where A and C are non-adjacent vertices. The set of all Markov equivalent DAGs is called a Markov Equivalence Class (MEC). We denote the set of all DAGs that is Markov equivalent to D as [D] . Definition 1 (Faithfulness [Zhang and Spirtes, 2012] ). If the population distribution exhibits a conditional independence relation only when there exists a corresponding d-separation statement in the causal graph, then we say that the population distribution is faithful to the causal graph.The positivity assumption is fundamental for making causal inferences. It asserts that, in theory, every individual possesses a non-zero chance of being both exposed and unexposed [Hernan and Robins, 2020] . In many cases, we have access to abundant observational data, which allows for precise estimation of the true observational distribution. The no positivity violation assumption is merely needed for theoretical guarantees to hold while not required by the algorithm. We assume access to the observational distribution, and the observational distribution is faithful to the true causal graph, with no positivity violations.A partially directed acyclic graph (PDAG) is a partially directed graph free from directed cycles. All Markov equivalent DAGs can be represented by a completed partially directed acyclic graph (CPDAG), denoted by C. A DAG D can be represented by a CPDAG when they share the same set of adjacencies and unshielded colliders, and every oriented edge in the CPDAG is also present in D [Meek, 2013] . CPDAGs are chain graphs with chordal chain components [Andersson et al., 1997] . In graph theory, a chordal graph is one in which cycles of four or more vertices always contain an additional edge, called a chord. We denote the set of all chain components of a PDAG P as CC(P).An intervention on a subset of variables W ⊆ V, denoted by the do-operator do(W = w), involves setting each W j ∈ W to w j . For every intervention, we have an induced post-interventional graph denoted by D S with incoming edges to vertices in S removed. We denote the interventional and observational distributions as P D w (v) and P D (v) respectively for a given DAG D. Using the truncated factorization formula over the post-interventional graph (D S ), we have the following:EQUATIONWhere v must be consistent with the intervention do(W = w). (V) to represent the interventional distribution in E(D) with the cut configuration C k (S i ). Also we use the notation C * (S i ) for the cut configuration in the true DAG D * . Given an intervention set S i , if we assume perfect intervention, one can recognize all the edges adjacent to vertices in S i and further apply Meek Rules. The resulting PDAG is called an interventional essential (I-essential) graph of D denoted as E S i (D). When S i is empty, it is the observational essential graph. I-essential graphs are also known to be maximally oriented partially directed acyclic graphs (MPDAGs). We denote an MPDAG as M and the I-essential graph with C k (S i ) as M C k (Si) . MPDAGs are also chain graphs with chordal components [Hauser and Bühlmann, 2012] . We denote the KL divergence between two distributions p, q as D KL (p||q). 
algorithm initializations: 
 Assume the true DAG is D * , our algorithm aims to return a most 'probable' causal graph D given N interventional samples. With the access to infinite observational data, we can retrieve the essential graph G = E(D * ) by PC algorithm [Spirtes et al., 2001] and joint distribution P (V). Here we describe how we compute the separating system of G and the causal effect of a given intervention set.As mentioned in Section 2, separating system plays the key role for non-adaptive intervention design of causal discovery. Roughly speaking, a separating system on a set of elements is a collection of subsets such that for every pair of elements from the set, there exists at least one subset which contains exactly one element from the pair. Consider an undirected complete graph G of n vertices indexed as 1, ..., n, a separating system S on [n] would cut every edge of G. In the worst case, a separating system is needed to learn the causal graph with G as its essential graph [Shanmugam et al., 2015] . If we further bound the size of each S ∈ S such that |S| ≤ k, k < n 2 , the resulting set is called an (n, k)-separating system. We provide the formal definition below: Definition 2 ((n, k)-Separating System [Katona, 1966 , Wegener, 1979] ). An (n, k)-separating system on [n] is a set of subsets S = {S 1 , S 2 , ..., S m } such that |S i | ≤ k and for every pair i, j there is a subset S ∈ S such that either i ∈ S, j / ∈ S or i / ∈ S, j ∈ S. We discuss the detailed steps of how we construct the (n, k)-separating system in Appendix D.To use a Bayesian approach, we need to construct a set of disjoint events and track their posteriors.Here we consider the events as interventional causal effects when we intervene on a set of vertices. By the assumption of faithfulness, Lemma 1 shows that the post-interventional distribution is determined by the edge configurations that are adjacent to vertices in the intervention set. Perkovic [2020] proposed a formula is to identify any causal effect in an MPDAG. One can enumerate all valid edge cuts and use the identification formula to calculate the post-interventional distributions. Here we propose a simple Algorithm 3 to enumerate through all possible configurations of a given set S and calculate the post-interventional distribution P C(S) s (V) via DAG sampling. For each candidate configuration C k (S), we check if it contains invalid structures, i.e. unshielded colliders or cycles. If it is valid, we apply Meek Rules to get the MPDAG M C k (S) which is an I-essential graph E S (D * ) for some DAG D C k (S) ∈ [D * ] that is consistent with C k (S). The chain components of M are chordal graphs and can be oriented independently [Hauser and Bühlmann, 2012] . To efficiently calculate the interventional distribution, we use a DAG sampler [Wienöbst et al., 2023] to sample a DAG for each chain component and replace the edges in M C k (S) with the arcs in the sampled DAG to get a fully oriented graph D. The sampling process takes linear time. Given the DAG D, the interventional distribution could be calculated by using the Equation 1 over the DAG. 
algorithm design and analysis: 
 In the previous section, we discuss two key steps for the algorithm initialization, which include the construction of (n, k)-separating systems and enumerating causal effects for all possible cutting edge configurations for all the intervention targets in the separating system. For a given intervention target S i , we utilize the following result that under the faithfulness assumption, shows that we have a unique interventional distribution for every cutting edge orientation C j (S i ).Lemma 1. [Elahi et al., 2024] Assume that the faithfulness assumption holds and D * is the true DAG. For any DAG D 1 ̸ = D * , if P D1 s = P D * s for some S ⊆ V, they must share the same cutting edge orientation C(S).Lemma 1 allows us to uniquely orient cutting edges for every intervention target in the separating system, which in turn orients every edge and gives us the true DAG. We use the notation D si ab := D KL ( PCa(Si) si (V)||P C b (Si) si (V) ) > 0 ∀a ̸ = b , ∀a, b ∈ [n Si ] and D si := min ∀a̸ =b , ∀a,b∈[n S i ] D siab . The idea is to use the interventional data do(S i = s i ) to determine the true cutting edge configuration C * (S i ) and repeat this for all intervention targets in the separating system. Suppose we have access to m si i.i.d. samples from the intervention do(S i = s i ), which we denote as Data do(si) = {v 1 , v 2 , ..., v ms i }. We define the posterior probabilities of C j (S i ) for all possible cutting configurations at S i , i.e., for all j ∈ [n Si ] where n Si ≤ 2 kdm as follows:Definition 3. (Posterior Probabilities of Cutting Edge Configurations) Consider an intervention target S i and a collection of all possible cutting edge configurations C j (S i ) for all possible cutting configurations at S i , i.e., for all j ∈ [n Si ] and the interventional dataset of i.i.d. samples Data do(si) = {v 1 , v 2 , ..., v ms i }. We define probabilities of Cutting Edge Configurations as:P (C j (S i ) | Data do(si) ) = P Cj (Si) si (v 1 , v 2 , ..., v ms i ) p j n S i a=1 P Ca(Si) si (v 1 , v 2 , ..., v ms i ) p a ∀j ∈ [n Si ] , ∀S i ∈ S (2)where p a for all a ∈ [n Si ] are the set of priors such thatn S i a=1 p a = 1.Given the fact that the interventional samples are i.i.d we can rewrite the posterior probabilities asEQUATIONAssumption 1. We assume access to the observational distribution that is faithful to the true causal graph. Furthermore, we have logP Ca (S i ) s i (v) P C b (S i ) s i (v)≤ β for every intervention set S i ∈ S and for alla, b ∈ [n Si ].The second half of the assumption can be seen as a version of the positivity assumption commonly used in causal discovery and many other applications. We show that for the problem of learning the true cutting edge configuration, the posterior is consistent, i.e., as the number of samples m si → ∞, the posterior of the true cutting edge configuration given the data converges to 1 with high probability.Lemma 2. (Posterior Consistency) Consider an intervention target S i ∈ S and the corresponding true cutting-edge configuration C * (S i ). As the number of samplesm si → ∞ in Data do(si) = {v 1 , v 2 , . . . , v ms i }, the posterior of the true cutting-edge configuration P (C * (S i ) | Data do(si) )converges to 1 with high probability. More precisely, we have the following high probability lower bound on the posterior probability of the true cutting-edge configuration.EQUATIONWhere α 1 and α 2 are constants depending on the prior and the problem instance. Thus, for any small choice of the probability δ, with a sufficiently large number of samples m si , the posterior of the true cutting-edge configuration P (C * (S i ) | Data do(si) ), converges to 1 with a probability at least 1 -δ.Note that the randomness in the high-probability result in Equation 4arises from the finite sample size. Lemma 2 guarantees that with a sufficiently large number of samples from interventions on S i , the posterior of the true cutting edge configuration P (C * (S i ) | Data do(si) ) converges to 1 with high probability. However, in scenarios where we want to know ahead of time how many interventional samples we need from every target S i ∈ S to ensure that the posterior of the corresponding true cutting orientations P (C * (S i ) | Data do(si) ) ≥ 1 -γ with high probability, we need to determine the required number of samples. The parameter γ ≥ 0 represents the error tolerance. We extend Lemma 2 to provide the minimum number of samples m si required for an intervention on the target set S i to ensure that the posterior of the corresponding true cutting orientations is greater than some threshold. Theorem 3. Given that the Assumption1 hold, consider an intervention target S i ∈ S such that |S i | ≤ k and the corresponding true cutting edge configuration C * (S i ). If the number of samples m si in Data do(si) = {v 1 , v 2 , ..., v ms i } satisfies the following:m si = 2β 2 (D si ) 2 ln 2 (k+1)dm δ + 2 D si ln 2 kdm (1 -γ)(1 -p * ) p * γ (5)where p * is the prior assigned to the true cutting-edge configuration C * (S i ), then we haveP (C * (S i ) | Data do(si) ) ≥ 1 -γ with a probability at least 1 -δ.Since we have multiple intervention targets in the (n, k)-separating system of the formS = {S 1 , S 2 , ..S p } such that |S i | ≤ k for all i ∈ [p].We can have a total of p bad events one for every target set S i ∈ S where posterior of corresponding true cutting edge configuration C * (S i ) is not greater than the desired threshold 1 -γ. Thus we extend the Theorem 3 as follows:Corollary 4. Given that Assumption1 hold, consider a separating system of the formS = {S 1 , S 2 , ..S p } such that |S i | ≤ k for all i ∈ [p].If the number of samples m si in Data do(si) = {v 1 , v 2 , ..., v ms i } for every target set S i ∈ S satisfies the following:m si = 2β 2 (D si ) 2 ln 2 (k+1)dm δ ′ + 2 D si ln 2 kdm (1 -γ)(1 -p * ) p * γ (6)where p * is the prior assigned to the true cutting edge configuration C * (S i ) and δ ′ = δ p , then we have P (C * (S i ) | Data do(si) ) ≥ 1 -γ with probability at least 1 -δ for every S i ∈ S.To ensure that for every target set S i ∈ S, we have the posterior probability of the corresponding true cutting-edge configuration P (C * (S i ) | Data do(si) ) ≥ 1 -γ, we need m si samples (as in Equation 6) from every target. Therefore, for the causal discovery problem with p targets in a separating system, we require a total of p • m si samples.Based on the analysis, we propose Algorithm 1 as the main algorithm. We start by calculating a (n, k)-separating system S using the labeling procedure described in Appendix D for the essential graph G = E(D * ), and them identify all valid configurations C(S), ∀S ∈ S using Algorithm 3. For each valid configuration, one can simply assume that their priors are uniform. In our algorithm, we assume that each possible DAG in [D * ] are equally likely to be the true DAG. Therefore we can calculate the prior more accurately using the MEC counting algorithm [Wienöbst et al., 2023] . Specifically, for each valid configuration C k (S), we find the MPDAG M C k (S) that matches with it by applying Meek Rules. M C k (S) is a chain graph with chordal components that can be oriented independently. Thus, the interventional MEC size of M C k (S) could be calculated by:|[M C k (S) ]| = H∈CC(M C k (S) ) |[H]| (7)The prior of each configuration could then be estimated by:EQUATIONIf we are allowed to collect N total interventional samples, we can randomly choose a target set from the (n, k)-separating system for each sample. Algorithm 3 iterates all valid configurations to calculate the likelihood P (v do(S) |C k (S)) by sampling DAG from the interventional MEC. After calculating all the likelihoods, we can update the posteriors and priors. For the last step, we need to combine all the configurations of each (n, k)-separating set to return a DAG. We consider a greedy approach. At each step, we choose the configuration with the highest posterior across all unvisited target sets that is compatible with chosen configurations. Given the anytime nature of our algorithm, we can propose a candidate DAG after each interventional sample. With a large enough sample number, according to Theorem 3, our algorithm will return the true DAG with a high probability.P C k (S) S (C k (S)|v do(S) ) ← P (C k (S))×P (v do(S) |C k (S)) C k (S)∈C(S) P (C k (S))×P (v do(S) |C k (S)) ; P (C k (S)) ← P C k (S) S (C k (S)|v do(S) ) D ← G, Svisit = ∅; while S ̸ = Svisit do C k (S) ← arg max C j (S)∈C(S),S∈S P C k (S) S (C k (S)|v do(S) ); if C k (S)6 Case Study: Estimating Causal Effects of Non-Intervenable VerticesHere we show how to use DAG sampler to estimate the causal effect when some vertex cannot be intervened in the causal graph. The detailed problem setup is as follow. D = (V, E) is an undirected underlying causal graph, X, Y ∈ V are non-adjacent vertices. We want to estimate p(y|do(x)) using interventional data while X is not intervenable in the graph. Thus we cannot directly use the method in Perkovic [2020] to estimate the causal effect. With the DAG sampler, we can first intervene on X's neighborhood N e(X) and then sample DAG to estimate the likelihood of each configuration being the true one in the causal graph. When we have large sample size, the posterior of the true configuration will go to 1 with a high probability according to Lemma 2. Then, we iterate through each configuration and their posterior to calculate the average divergence D between the estimated and true causal effect. The divergence D here can be KL divergence or Total Variation Distance (TVD). If the ground truth causal effect is p * , the interventional data is Data do(N e(X)) , then D is formally defined as:EQUATIONIn this experiment, we randomly create 50 causal graphs with n = 5, 6, 7, ρ = 0.3, 0.6 using a similar approach described in Section 7. Then we randomly choose a pair of non-adjacent vertices X, Y . We then intervene on X and collect 100,000 samples. We plot D KL and D T V D between estimated causal effect and the ground truth given different number of interventional samples. The mean and standard deviation are plotted in Figure 1 . We can see that as the number of interventional samples increases, the estimated causal effect gets close to the ground truth and the variance is also decreasing. The difference decreases sharply w.r.t. the number of samples, showing the efficiency of our proposed approach. We generate random connected moral DAGs with order n and density ρ using a modified Erdős-Rényi sample approach, similar to the process in Squires et al. [2020] . We start by generating a random ordering τ over vertices. Then, for the i th vertex, we sample its in-degree as X i = max(1, Bin(n -1, ρ)). The vertices precede it in the ordering are uniformly assgined as its parents. In the last step, we apply the elimination algorithm in Koller and Friedman [2009] to chordalize the graph. The elimination algorithm uses an ordering that is the reverse of τ . Based on the generated graph, we randomly sample the conditional probability table (CPT) that is consistent with the graph and strictly positive. The essential graph of the generated DAG is then fed into the causal discovery algorithms together with the Bayes network that is consistent with the DAG.To measure the performance of the algorithms, we plot the Structural Hamming Distance (SHD) between the ground truth and learned DAGs with respect to the number of collected interventional samples. SHD counts the number of edge insertions, deletions, or flips to transform from one DAG into another. For each setting of n and ρ, we sample 50 random DAGs and calculate the mean and standard deviation of SHD by each causal discovery algorithm. The mean is plotted as a curve and the shaded area around the curve represents the standard deviation. We use Chi-Square independence test from the Causal Discovery Toolbox [Kalainathan et al., 2020] to perform statistical tests with limited samples.The results in Figure 2 and 3 show that our algorithm outperforms other causal discovery algorithms.Figure 2 shows the performance of causal discovery algorithms on complete graphs with order 5, 6, and 7 respectively. Our algorithm uses significantly less samples to reach a low SHD and the number of samples required does not increase much with the order of the graph compared to the other algorithms. Figure 3 shows the results of graphs with the same order but varying densities. Similarly, our algorithm uses comparable number of samples across different densities. In contrast, the baseline algorithms require much more samples in denser graphs to reach the same SHD. To wrap up, in both cases, our algorithm reaches a low SHD using much fewer samples than the baselines and is more stable across different settings. The results on larger graphs (n = 20) are provided in Appendix H.1. The algorithm code is provided at https://github.com/CausalML-Lab/ Bayesian_SampleEfficient_Discovery.  
conclusion: 
 In this work, we discuss the problem of learning causal graphs with minimum number of interventional samples. We propose an algorithm that solves this problem in a Bayesian approach. Specifically, we keep track of the posteriors of interventional distributions of each target set in the (n, k)-separating system and then combine the configurations that have high posteriors to propose a DAG. We show in theory that given enough samples, our algorithm can return the true causal graph with a high probability. Also, according to the Bayesian nature of our algorithm, we can stop the algorithm anytime and check the output graph. Experiments on simulated datasets show that compared with the baselines, our algorithm use significantly fewer interventional samples to achieve the same SHD.Additionally in the case study, we demonstrate that we can modify our algorithm to answer special causal queries with DAG sampling. For future extensions, it is of great interest to decrease the space of posteriors being tracked, since on large dense graphs, it will be intractable to track all interventional distributions of a target set. Besides, we can try to remove/weaken the assumptions like causal sufficiency and faithfulness while maintaining theoretical guarantees. Our work can potentially have some societal consequences including ethical considerations related to performing interventions and the risk of biased or partial understandings which may cause misled decision-making in real-world cases. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: We do not use any dataset or model that has this concern.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We cite the assets and mention the license explicitly in the paper and Readme of the code.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [Yes] Justification: We provide the descriptions in the Readme file. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: A Proof of Lemma 2:Consider an intervention target S i and a collection of all possible cutting edge configurations C j (S i ) for all possible cutting configurations at S i , i.e., for all j ∈ [n Si ] and the interventional dataset of i.i.d. samples Data do(si) = {v 1 , v 2 , ..., v ms i }. We revisit some notations from the main paper. We denote the interventional distribution with cut configuration C j (S i ) by P Cj (Si) si. we have n Si different cut configurations which we denote by C 1 (S i ), C 2 (S i ) .... C n S i (S). Under faithfulness assumption we have PCa(Si) si ̸ = P C b (Si) si for all a ̸ = b , a, b ∈ [n Si ]. Also note that for with (n, k) separating system we have n Si ≤ 2 kdm . Without loss of generality assume C 1 (S i ) is the true cutting edge configuration. We use the notation D si ab := D KL ( PCa(Si) si (V) | | P C b (Si) si (V) ) > 0 ∀a ̸ = b , ∀a, b ∈ [n Si ] and D si := min ∀a̸ =b , ∀a,b∈[n S i ] D si ab .Let us consider P C1(Si) si and Ponly for now.EQUATIONSince interventional dataset is composed of i.i.d. samples we have:EQUATIONEQUATIONEQUATIONSince C 1 (S i ) is the true cutting edge configuration so the interventional samples are actually sampled from P C1(Si) si. So we have the following:E[L j ] = E vj ∼P C 1 (S i ) s i -log P C2(Si) si (v j ) P C1(Si) si (v j ) (15) E[L j ] = E vj ∼P C 1 (S i ) s i log P C1(Si) si (v j ) P C2(Si) si (v j ) (16) E[L j ] = D KL ( P C1(Si) si (V) | | P C2(Si) si (V) ) = D si 12 > 0 (17) Under the Assumption | L j | ≤ β ∀j ∈ [m si ] for all L j = -log P C 2 (S i ) s i (vj ) P C 1 (S i ) s i (vj ). using hoeffding inequality with m si samples we have:EQUATIONNow, suppose we assign a non-zero prior p j > 0 to the all cut configuration C j (S i ) ∀j ∈ [n S I ] such that n j=1 p i = 1. The posterior distribution can be written as:EQUATIONFrom the result in Equation 24 with probability at least 1 -δ we have:EQUATIONWith probability at least 1 -δ we have:P (C 1 (S i ) | Data do(si) ) ≥ p 1 p 1 + n S i j=2 e -ms i D s i 1j -2β 2 ms i ln 2×2 kdm δ p j (27)Since n Si ≤ 2 kdm and D si := min ∀a̸ =b , ∀a,b∈[n S i ] D si ab . With probability at least 1 -δ we have:P (C 1 (S i ) | Data do(si) ) ≥ p 1 p 1 + (1 -p 1 )2 kdm e -ms i D s i -2β 2 ms i ln 2×2 kdm δ (28) P (C 1 (S i ) | Data do(si) ) ≥ 1- 1 1 + p1 (1-p1)2 kdm e ms i D s i -2β 2 ms i ln 2×2 kdm δ (w.p. at least 1-δ) (29)which can equivalently written as follows:EQUATIONWhere α 1 and α 2 are constants depending on the priors used and the problem instance. Note the in our proof the true cutting edge configuration C * (S i ) = C 1 (S i ). This completes the proof of Lemma 2.The Proof of the Theorem 3 is the continuation of the proof of Lemma 2. Suppose we want P (C 1 (S i ) | Data do(si) ) ≥ 1 -γ we have:EQUATIONm si D si -2β 2 m si ln 2 × 2 kdm δ ≥ log 2 kdm (1 -γ)(1 -p 1 ) p 1 γ (33)Solving the above equation for number of interventional samples m si we have:EQUATIONIn order to simplify the expression we select the number of interventional samples m si as follows:m si = 4β 2 ln 2 (k+1)dm δ +4D s i ln 2 kdm (1-γ)(1-p 1 ) p 1 γ + 4β 2 ln 2 (k+1)dm δ +4D s i ln 2 kdm (1-γ)(1-p 1 ) p 1 γ 2 4(D s i ) 2 (35) m si = 4β 2 ln 2 (k+1)dm δ + 4D si ln 2 kdm (1-γ)(1-p1) p1γ + 2(D si ) 2 (36) m si = 2β 2 (D si ) 2 ln 2 (k+1)dm δ + 2 D si ln 2 kdm (1 -γ)(1 -p 1 ) p 1 γ (37)for the above choice of number of samples m si we have the following result:EQUATIONNow, suppose p * is the prior assigned to the true cutting-edge configuration C * (S i ). Consequently, for the following choice of number of samples m si :m si = 2β 2 (D si ) 2 ln 2 (k+1)dm δ + 2 D si ln 2 kdm (1 -γ)(1 -p * ) p * γ (39)We have P (C * (S i ) | Data do(si) ) ≥ 1 -γ with probability at least 1 -δ. Note the in our proof the true cutting edge configuration C * (S i ) = C 1 (S i ). This completes the proof of Theorem 3.C Proof of the Corollary 4: 39for number of samples m si then for a particular S i ∈ S we haveP (C * (S i ) | Data do(si) ) ≥ 1 -γ with probability at least 1 -δ p .Taking union bound we have the following result:If the number of samples m si in Data do(si) = {v 1 , v 2 , ..., v ms i } for every target set S i ∈ S satisfies the following:m si = 2β 2 (D si ) 2 ln 2 (k+1)dm δ ′ + 2 D si ln 2 kdm (1 -γ)(1 -p * ) p * γ (40)where p * is the prior assigned to the true cutting edge configuration C * (S i ) and δ ′ = δ p . Then we have with P (C * (S i ) | Data do(si) ) ≥ 1 -γ with probability at least 1 -δ for every target set S i ∈ S. This completes the proof of Corollary 4.D.1 (n, k)-Separating System Lemma 5 ( [Shanmugam et al., 2015] ). There exists a labeling procedure that gives distinct labels of length l for all elements in [n] using letter from the integer alphabet {0, 1, ..., a}, where l = ⌈log a n⌉. Furthermore, in every position, any integer letter is used at most ⌈ n a ⌉ times. The labeling method in Lemma 5 from Shanmugam et al. [2015] is described as: Labelling Procedure: Let a > 1 be a positive integer. Let x be the integer such that a x < n ≤ a x+1 .x + 1 = ⌈log a n⌉. Every element j ∈ [n] is given a label L(j) which is a string of integers of length x + 1 drawn from the alphabet {0, 1, ..., a} of size a + 1. Let n = p d a d + r d and n = p d-1 a d-1 + r d-1 for any integers p d , p d-1 , r d , r d-1 , where r d < a d and r d-1 < a d-1 . Now, we describe the sequence of the d-th digit across the string labels of all elements from 1 to n:1. Repeat the integer 0 a total of a d-1 times, and then repeat the subsequent integer, 1, also a d-1 times from {0, 1, ..., a -1} till p d a d . 2. Following this, repeat the integer o a number of times equal to ⌈ r d a ⌉, and the repeat the integer 1 ⌈ r d a ⌉ times, continuing this pattern until we reach the nth position. It is evident that the nth integer in the sequence will not exceed n -1. 3. Each integer that appears beyond the position a d-1 p d-1 is incremented by 1.Once we have a set of n string labels, we can easily construct a (n, k)-separating system using the following Lemma: Lemma 6 ( [Shanmugam et al., 2015] ). Consider an alphabet A = [0 :⌈ n k ⌉] of size ⌈ n k ⌉ + 1 where k < n 2 .Label every element of an n element set using a distinct string of letters from A of length l = ⌈log ⌈ n k ⌉ n⌉ using the labeling procedure in Lemma 5 with a = ⌈ n k ⌉. For every 1 ≤ a ≤ l and 1 ≤ b ≤ ⌈ n k ⌉, we choose the subset I a,b of vertices whose string's a-th letter is b. The set of all such subsets S = {s a,b } is a k-separating system on n elements and |S| ≤ (⌈ n k ⌉)⌈log ⌈ n k ⌉ n⌉.While an (n, k)-separating system separates every edge in the complete graph of order n, usually the causal graphs are sparse with much fewer edges. Thusly, we also provide a simple algorithm to construct a separating system that cuts all the edges in a graph.Definition 4 (G-Separating System, Definition 3 of Kocaoglu et al. [2017] ). Given an undirected graph G = (V, E), a set of subsets I ⊆ 2 V is a G-separating system if for every edge {u, v} ∈ E, there exists I ∈ I such that either (u ∈ I i and v / ∈ I i )or (u / ∈ I i and v ∈ I i ). Kocaoglu et al. [2017] show that an intervention set I learns every DAG entailed by G if and only if it is a G-separating system. For a graph with n vertices, a G-separating system S = S 1 , S 2 , ..., S m , such that|S i | ≤ k; ∀i ∈ [m], is called a (n, k)-separating system. There are many ways to get a G-separating system. Algorithm 2 gives a simple algorithm that uses vertex coloring. We first find the perfect elimination order (PEO) L of G. PEO can be found by maximum cardinality search [Tarjan and Yannakakis, 1984] . Then based on the reverse PEO, we apply greedy coloring to assign each vertex with a number i ∈ [ω] where ω is the clique number of G. ω colors are guaranteed to color G since G is a perfect graph. In each step, we assign a vertex with the minimum color number that is not used by its visited neighbors, thus no adjacent vertices will have the same color. For each vertex, the union of it and its visited neighbors induce a clique, meaning that its color number is bounded by ω. Consequently, Algorithm 2 returns a G-separating system with ω intervention sets as it cuts every edge in G. We use f (•) to represent a function that maps a set of vertices to a set of colors that are used by the vertices. N (v) denotes the neighborhood of a vertex v ∈ V.Algorithm 2: G-separating system of a given unoriented graph GData: Input UCCG G = (V, E), |V| = n Result: Output the G-separating system S Initialization: Find the PEO L and ω of G by Maximum Cardinality Search; L ← Reverse L, Si = ∅, ∀i ∈ [ω], f (V) = {0}; for i = 1 to n do ci ← min{[ω] -f (N (Li))}; f (Li) ← {ci}; Sc i ← Sc i ∪ {Li}; S = ∪ i∈[ω] {Si}; return S;E Algorithm for Enumerating Causal Effects For Bayesian approach of learning causal graph, we mainly care about the space complexity because for run time, all the steps (e.g. calculating likelihood, updating posteriors) are efficient. For dense graphs, a fully Bayesian method and our algorithm are both computationally expensive. For example, for a complete graph with order n, Fully Bayesian method needs to store O(n!) distributions because there are such many DAGs in the MEC. Our method needs to track O(2 n-1 ) posteriors, which is less than fully Bayesian but still exponentially many and would be intractable when n is large. However, our method would be a lot better when the given graph is sparse, i.e. the maximum degree d in the UCCG G is around log n. If we consider atomic interventions, the valid configurations around a vertex v ∈ V, |C(v)| < 2 d = n. While the MEC size of a UCCG with order n is at least n, i.e., |[G]| ≥ n. The equality holds if and only if G is a tree. Thus, our method is at least as good as a fully Bayesian method. However, even for sparse graphs, its MEC size could be exponentially many in n. Let's consider the following example. Given a sparse UCCG G with order n. G consists of [ n [log n] ] small cliques. Each small clique contains [log n] vertices. The cliques are connected together to form a 'clique line'. The MEC size of G would be:EQUATIONAccording to Stirling's Approximation, we have:EQUATIONWhich is approximately O(n log[log n] ). If we consider an interventional target of size k for our algorithm, there are at most2 kd = n k configurations to track. If n is large enough, k < log[log n], ∀k > 0.Thus, the MEC size for G is greater than any polynomial complexity, while our algorithm tracks at most polynomially many configurations.In Section 4.2, we mentioned that the Causal Identification Formula for MPDAG from Perkovic [2020] could also be used when enumerating the causal effects. Here we show that our use of DAG sampler provides a more flexible and general approach of causal discovery with special queries that do not rely on the cut configurations of the MPDAG. In many real-world settings, our focus is on specific queries rather than understanding the entire graph. For instance, when a cloud computing system encounters an error, we aim to efficiently identify the root cause. Similarly, when a patient exhibits symptoms of cancer, we want to determine if certain attributes of the patient are direct causes of the cancer. Here, we consider the scenario mentioned in Malinsky [2024] , where the goal is to estimate the posterior probability of a set of variables being the adjustment set of a causal query given the essential graph. The adjustment set is crucial in estimating the causal effect. Given the essential graph G, denote the variable of a set S being the adjustment set as A S . A S = 1 if S is the adjustment set and 0 otherwise. We want to estimate P (A S |v), where v could be any data collected, both observational and interventional, Denote the true DAG as D * , the posterior could be expressed as: P (A S |vAccording to Hoeffding Inequality, the difference between S and E S could be bounded byP ( S M -E S ≥ ϵ) ≤ exp(-2M ϵ 2 ) (49)By estimating the posteriors for both P (A S = 1|v) and P (A S = 0|v), we can normalize them to get the exact posterior. For large and dense graphs which are hard to iterate through each DAG, we could estimate efficiently by sampling DAGs according to intended accuracy. With enough DAG samples, it shows that the estimation converges to the true posterior with a high probability.We further test our proposed algorithm against the baseline methods on larger random chordal graphs.The results are shown in Figure 4 . We can see that our algorithm converges very fast to a low SHD.The performance of our algorithm does not differ much from that in small graphs. The baselines, on the other hand, require a lot more samples to reach the same SHD as in small graphs.For extremely large graphs, it would be intractable to store the whole observational distribution, instead we could all the conditional distributions of the factorizations. Given large enough observational samples, the conditional distributions will be accurate and then we can perform the Algorithm 1 for Bayesian learning. In this experiment, we compare our proposed algorithm with a Bayesian causal discovery method, AVICI, proposed in Lorch et al. [2022] . For AVICI, we fine-tuned the pretrained models, scm-v0 and neurips-rff, with 50 random complete graphs, each with 1000 observational samples. The results are plotted in Figure 6 . We used SHD to measure our algorithm and expected SHD for AVICI. Here, since AVICI can not take in too many samples, we just show the performance with 1000 interventional samples. At each step for both methods, we choose the interventional target randomly. For AVICI, we feed in 1000 observational samples together with the interventional samples. The results show that our algorithm performs better, while AVICI failed to learn the causal structure with high expected SHD. AVICI's bad performance in this setting may be caused by its parametric model on the SCM or noise. Also, there is no theoretical guarantee on AVICI's performance and convergence to the ground truth. Furthermore, our algorithm benefits has the advantage of being able to update the posterior given any number of samples, while AVICI would fail when the number of samples is large. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 


Decision: This work presents Bayesian approaches to the challenging problem of active learning of causal graphs via interventions. Reviewers were overall positive and I recommend acceptance. That said, reviewers did raise issues with various aspects of the writing - authors please do incorporate the recommendations and address the comments as promised in the rebuttal.

Reviews: review: paper_summary: paper_summary: In this paper, the authors present a Bayesian method to learn the causal graph from observational data and limited interventional data. The authors assume that there are plenty of observational data, which can be used to learn a ground-truth CPDAG. Then, by using the efficient DAG enumeration method to sample DAGs and calculate the posterior, one can find the DAG that is most consistent with the interventional data.
strength_weakness: strength_weakness: Using Bayesian method to learn causal relations is not novel. There are many existing studies in the literature. Further, from the viewpoint of considering limited samples, it is more sensible to take the uncertain in learning the CPDAG into account.

The writing could be improved. There are a lot of contents in introduction that are better to appear in Related works. Besides, it is not quite clear that what role Section 4.1 play in Section 4: Algorithm Initializations.

It is better to distinguish \citet and \citep.
questions: questions: It seems that Algorithm 3 could be improved further. If the authors orient the edges of $\mathbf{S}$ one by one and update the graph with Meek rules (due to the completeness of Meek rules in incorporating background knowledge), the complexity will be reduced and there is no need to detect unshielded colliders or cycles. Am I right?
limitations: limitations: No.
paper_summary: paper_summary: The paper considers the problem of learning causal graphs using limited interventional samples through the a Bayesian perspective. An algorithm is proposed which returns the most probable causal graph given a limited set of samples. The approach is empirically evaluated and code is given.
strength_weakness: strength_weakness: The experiments feels a little "small scale". This is probably due to the need to enumerate all possible configurations via Algorithm 3 for this approach to work (please correct me if I am mistaken).
questions: questions: - Lines 39-41 discusses the trade-off between non-adaptive and adaptive methods. This was explored in [1], a paper you have already cited. You may want to consider adding a reference to [1] in this section of the introduction in your revision.
- Line 104: Typo of $O(\log n)$?
- Line 104: Is it size 1 or size $n/2$ interventions? For example, see Page 58 of [2] which you cited.
- Line 243 and 514: Do you mean $D_{KL}$? Also, you should define KL divergence properly somewhere in the preliminaries.
- Equation (2): Do you mean $p_j$ in the numerator?
- Line 605: Typo of "Enumerating"?

[1] Davin Choo and Kirankumar Shiragur. "Adaptivity complexity for causal graph discovery." Uncertainty in Artificial Intelligence. PMLR, 2023.

[2] Frederick Eberhardt. "Causation and intervention." Unpublished doctoral dissertation, Carnegie Mellon University 93 (2007).
limitations: limitations: Nil
paper_summary: paper_summary: This paper proposes a Bayesian approach for learning causal graphs with limited interventional samples. The proposed algorithm first constructs a separating system to design intervention targets and then enumerates the causal effects for all possible cutting edge configurations for each target, and tracks their posteriors. The authors provide theoretical guarantees on the convergence to the true causal graph with sufficient interventional samples, and the experiments on simulated chordal graphs demonstrate that the proposed method requires significantly fewer interventional samples than baselines to achieve low SHD.
strength_weakness: strength_weakness: Major:
1. The proposed method should not be classified strictly as "Bayesian learning of causal graphs," because it does not calculate the posterior distribution over causal graphs directly. Instead, it computes posterior probabilities of specific cutting edge configurations within the graph. Based on these probabilities, the proposed method updates the output DAG accordingly. It would be better if the authors clarified that in the paper to avoid any confusion.

2. It is unclear to me how the proposed method selects intervention targets at each step. I wonder if the authors could further clarify whether the proposed method also designs the intervention targets or the interventions are performed randomly.

3. The authors did not discuss the computational complexity in the paper. The sizes of nodes in the experiments are also relatively small. I wonder if the authors could discuss the complexity of the proposed method and whether it could scale up to graphs with a larger number of nodes.

4. In the experiments, the authors only consider simulated chordal graphs, which naturally fit the proposed method. I think the authors should conduct more experiments on different types of graphs (e.g., scale-free graphs) and perhaps some semi-synthetic graphs (e.g., graphs simulated using realistic simulators [1]).

5. The authors only consider 3 baselines in the experiments, and I wonder if the authors could compare the proposed method with Bayesian causal discovery methods that can handle interventional data (e.g., [2], [3]).


Minor: 
1. Missing related work: [3].

2. Please use \citet and \citep correspondingly rather than only using \citet. For example, in line 72: Meek Rules Meek [1995] --> Meek Rules [Meek, 1995] 


[1] Dibaeinia, P., & Sinha, S. (2020). SERGIO: a single-cell expression simulator guided by gene regulatory networks. Cell systems, 11(3), 252-271.

[2] Lorch, L., Sussex, S., Rothfuss, J., Krause, A., & Schölkopf, B. (2022). Amortized inference for causal structure learning. Advances in Neural Information Processing Systems, 35, 13104-13118.

[3] Hägele, A., Rothfuss, J., Lorch, L., Somnath, V. R., Schölkopf, B., & Krause, A. (2023, April). Bacadi: Bayesian causal discovery with unknown interventions. In International Conference on Artificial Intelligence and Statistics (pp. 1411-1436). PMLR.
questions: questions: 1. Please see the questions in the Weaknesses part.

2. How is the value of $k$ of the $(n, k)$-separating system determined in the experiments? 

3. How is the prior of each configuration defined?
limitations: limitations: The authors adequately addressed the limitations of their work.
paper_summary: paper_summary: This paper introduces a way of doing causal discovery from interventional samples using bayesian inference. They prove some theoretical results of convergence of the task at hand. Finally, they test their approach against common benchmarks on synthetic data.
strength_weakness: strength_weakness: In my opinion, the related work section is not written very well. It reads like a list with one sentence of what they are doing and not putting your research in the context of previous research. The only exception is line 124.
questions: questions: I’m still unsure about the overall significance of the paper. Yes, there is section 6 where the authors give a case study for estimating the probability of a set being a valid adjustment set. Furthermore, the authors briefly explain the case of cloud computing and cancer research. However, I don’t understand how that relates to querying interventional data from the system. Could the authors please expand on this? I have not reduced the score from this because I currently believe the strength of the paper are the theoretical results and the approach itself.
limitations: limitations: See questions

Comments: title: Final Thoughts
Thank you
Official Comment by Reviewer i7ee
Official Comment by Reviewer i7ee
Answer to response
comment: We thank all reviewers for providing thoughtful and constructive feedback. In response, we have carefully considered each of the points raised and made significant revisions to the manuscript to address the concerns. Specifically, we added extra experiments required by Reviewer i7ee to test the performance of proposed algorithm on scale-free graphs and compare our proposed algorithm with AVICI proposed in **[1]**. In the revision, we will do the following changes suggested by the reviewers:

- We will also fix the typos and citation formats raised by the reviewers. 
- We will take the suggestions on writing by Reviewer uDQi and HSDU.
- We will clarify the objective of this work as suggested by Reviewer i7ee.
- We will improve the performance of AVICI, add experiments to validate section 6, and include experiments using real-world simulators as suggested by Reviewer i7ee.

Here, we briefly conclude our work. In this work, we propose a novel sample-efficient causal discovery algorithm using a Bayesian approach. It is well-known that a fully non-parametric Bayesian approach would be intractable as the MEC size can be large. Previous Bayesian approaches in the literature have some parametric assumptions such linear SCM or additive Gaussian noise. Such methods will fail when the assumption do not hold. To handle this, our proposed algorithm uses $(n, k)$-separating system as interventional targets and further partition the DAGs according to the edge configuration of the target set, which makes the Bayesian learning tractable. Our algorithm also has the anytime property, i.e., predict an optimal graph with any amount of interventional samples. Additionally, we show in theory that our proposed algorithm is guaranteed to converge to the true graph both asymptotically and non-asymptotically.  In section 6, we briefly show that our method could be extended to estimate more general causal queries efficiently using DAG sampler proposed in **[2]**.

We believe these changes will strengthen the manuscript and clarify the contributions of our work. We appreciate the opportunity to improve our paper and are confident that the revised version better reflects the significance and rigor of our research.

Thank you again for your time and effort in reviewing our submission. We look forward to your final decision.

 **[1]** *Lorch, Lars, Scott Sussex, Jonas Rothfuss, Andreas Krause, and Bernhard Schölkopf. "Amortized inference for causal structure learning." Advances in Neural Information Processing Systems 35 (2022): 13104-13118.*

**[2]** *Wienöbst, Marcel, Max Bannach, and Maciej Liśkiewicz. "Polynomial-time algorithms for counting and sampling Markov equivalent DAGs with applications." Journal of Machine Learning Research 24, no. 213 (2023): 1-45.*
We thank the reviewer for the increasing the score. We want to briefly clarity the concerns below.

Although Bayesian causal discovery is not new, a fully non-parametric Bayesian approach is not tractable. Our work introduces a novel idea by using the $(n, k)$-separating system, followed by an enumeration of causal effects to partition the set of all possible DAGs into tractable sets. At the prediction stage, we combine the configurations with high posterior to return a DAG. Furthermore, we show in theory that our method is guaranteed to converge to the true DAG as sample size $m\rightarrow \infty$ in Lemma 1. Additionally, we show that when the sample size is large, the predicted graph would be the true graph with a high probability as shown in Theorem 3. This is the first work to show such guarantees in the context of sample-efficient causal discovery problems under mild faithfulness assumptions.
Thank the authors for the response. After reading the rebuttal of the authors and the other reviewers' comments, I think the theoretical result in Thm. 3 is somewhat interesting, thus I increase my score to 4. However, causal discovery using a Bayesian method is not firstly studied. And I agree that using the result of Wienöbst et al. is a nice step, but it is not contributive enough. Hence I insist on my negative score.
We thank the reviewer for increasing our score but would like to briefly address the concerns highlighted by reviewer.

**AVICI experiment details**

In our experiment, we did not fine-tune the pre-trained models using observational samples. We directly feed them with the interventional samples. We will try to fine-tune the pre-trained model first with the observational samples.

**Question of real-world application**

In real-world applications, if we have access to large number of observational samples, we can apply our algorithm on a set of candidate CPDAGs, We will also add a few simple experiments to validate the case study in section 6.

**Experiments with real-world simulators**

We will add experiments using real-world simulators.
Thank you for the explanation and response. I now understand the main objective of the proposed method.

1. Regarding the experimental details, I would like to clarify whether "For AVICI, we feed 1000 observational samples to the pre-trained models" means that you fine-tuned the pre-trained AVICI using these 1000 observational samples. If the pre-trained AVICI was directly applied to the simulated graphs without fine-tuning, the poor performance could be explained by a significant difference between the distributions of the simulated graphs and the pre-trained graphs. However, if the pre-trained AVICI was fine-tuned with the observational samples, it is surprising that AVICI's performance is still so poor (a SHD of 8 or 9, which suggests the predictions were almost entirely incorrect). 

2. I am still not very clear about when and how the proposed method can be applied in real-world scenarios. From my understanding, the proposed method is best suited for inferring causal graphs when there are limited intervention samples with known intervention targets. However, as the proposed method requires the observational distributions to be known, I am not sure whether this assumption can be easily fulfilled in real-world applications (e.g., gene regulatory networks). I feel like the case study in Section 6 might be a good application of the proposed method, but the authors did not perform experiments for validation. 

3. I still feel that the experiments should include real-world (or realistic simulated) experiments to demonstrate the effectiveness of the proposed method.

I appreciate the authors' effort and time during the rebuttal, therefore I raised my score by 1. However, I still feel that this paper needs more preparation, especially in describing the main objective, applying the theoretically sound method to real-world applications, and relaxing the assumptions to allow the causal discovery for a larger number of variables.
We thank the reviewers for the constructive reply. Below we address each of the points mentioned.

**Question of objective**

In our setup, we begin by constructing a set of targets to intervene on. We use an $(n, k)$-separating system here because it is guaranteed to learn every graph in the MEC. During the discovery process, we randomly choose a target from the set, so it is different from the experimental design works where interventional targets are adaptively selected based on previous interventions. The basic purpose of this work is to come up with an anytime algorithm that predicts an optimal graph given limited interventional samples. We will clarify this in the revised manuscript.

**Question metrics in the additional experiment**

For the results in Figure 2 of extra experiments, we just use SHD of the prediction as done in Avici's tutorial. More specifically, a threshold of $0.5$ is used to filter the predicted adjacency matrix. We agree that expected SHD would be a better metric to use here for Avici. We can modify this part in revision. 

**Questions about experiment details**

**Intervention targets:** For both methods, we choose target randomly for each sample. For Avici, we feed 1000 observational samples to the pre-trained models before feeding interventional samples.

**Graph generation:** We generate the graph following the steps described in section 7, second paragraph. More specifically, here, we first generate a random ordering $\tau$ of $[5]$. Then, we orient $a \rightarrow b$ if $a$ is previous to $b$ in $\tau$ for $a \neq b, a, b \in [5]$. Since the complete graph is chordal, we do not need to perform the last chordalizing step. A DAG $\mathcal{D}$ randomly generated in this way is guaranteed to have a chordal skeleton $G$ and $\mathcal{D} \in [G]$. We do notice that the performance of Avici is highly related to the pre-trained model it is using. We can try to improve Avici's performance by modifying the training part.

**Initial SHD:** For a complete graph with $5$ vertices, there will be $10$ edges, thus SHD is bounded by $10$. Our method has an initial SHD of $5$, which is basically a random guess of all the edges. However, Avici's metric uses a threshold of $0.5$, which probably masked out most entries in the adjacency matrix, resulting in a high SHD.

**Question about complexity:**

We agree with the reviewer that one can obtain observational distributions from real-world simulators, which is similar to synthetic generation. We briefly discuss the computational complexity of our proposed algorithm. Our approach has two major steps. In the first step, we initialize the targets and compute the interventional distributions. In the second step, we merely sample from the intervened Bayes net, put the sample into the interventional distributions, and update the prior and posteriors. The most computationally expensive part is saving the observational and interventional distributions. Since we are using binary variables in our simulation, the joint distribution of $50$ vertices would be a table of $2^{50}$ float numbers, which is intractable. If there is a compact way to represent the joint distributions, our algorithm will still work. We will add this clarification in the revised manuscript.

We hope we have addressed the questions and would be please to discuss further if reviewer has further concerns.
Thank you for the detailed response and the additional experiments. However, I still have some questions regarding the proposed method and the experimental setup.

1. From my understanding, the proposed method consists of two steps: intervention design and causal discovery. It is unclear whether the method is intended to generate higher-quality interventions or to infer the causal graph more accurately with limited interventional samples compared to other causal discovery methods. If the focus is on the former, the method may be more aligned with causal experimental design rather than causal discovery. Could you please elaborate on this point?

2. Regarding the additional experiments, since AVICI is a Bayesian causal discovery method, did you use the expected SHD as the measure? The figures only show SHD, so it would be helpful if you could provide more details about the experimental setup.

3. When comparing with AVICI, did you use the same interventions for both methods? Specifically, how were the interventions for AVICI generated—were they random or derived from the proposed method? Additionally, how did you generate random complete graphs? If pre-trained AVICI models were used, the prediction accuracy could drop significantly if the distribution of the generated graphs differs from the pre-trained set. From Figure 2 in the newly added document, it seems that the SHD of the proposed method is much smaller than that of AVICI at the initial stage. Could you clarify this?

4. Finally, while complexity is a crucial factor in real-world scenarios, the proposed method seems computationally expensive. Additionally, I believe we can obtain observational distributions from real-world simulators, which is similar to synthetic generation (e.g., scale-free). Could you address the computational aspects of your method in this context?

I am willing to adjust my score if these concerns are addressed.
Thank you for your responses. I am still not convinced by the size of the experiments but I do see value in the proposed approach. For instance, I felt it was interesting that they used the method of Wienöbst for efficiency. As such, I will keep my positive score as it is.
I thank the authors for taking the time to answer to my questions and for taking the suggestion to improve the related work section. Give my current understanding of the paper and the authors' responses I will keep my score as is.