ABCD-Strategy: Budgeted Experimental Design for Targeted Causal Structure Discovery
Raj Agrawal Chandler Squires Karren Yang Karthik Shanmugam Caroline Uhler







Targeted Causal Structure Discovery






Agrawal, Squires, Yang, Shanmugam, Uhler

[

ABCD-Strategy: Budgeted Experimental Design 
for Targeted Causal Structure Discovery

 Raj Agrawal  Chandler Squires  Karren Yang  Karthikeyan Shanmugam    Caroline Uhler 



 MIT MIT MIT MIT-IBM Watson AI Lab 
 IBM Research NY MIT



]


Determining the causal structure of a set of variables is critical for both scientific inquiry and decision-making. However, this is often challenging in practice due to limited interventional data. Given that randomized experiments are usually expensive to perform, we propose a general framework and theory based on optimal Bayesian experimental design to select experiments for targeted causal discovery. That is, we assume the experimenter is interested in learning some function of the unknown graph (e.g., all descendants of a target node) subject to design constraints such as limits on the number of samples and rounds of experimentation. While it is in general computationally intractable to select an optimal experimental design strategy, we provide a tractable implementation with provable guarantees on both approximation and optimization quality based on submodularity. We evaluate the efficacy of our proposed method on both synthetic and real datasets, thereby demonstrating that our method realizes considerable performance gains over baseline strategies such as random sampling. 




§ INTRODUCTION
 
Determining the causal structure of a set of variables is a fundamental task in causal inference, with  widespread applications not only in artificial intelligence but also in scientific domains such as  biology and economics <cit.>. 
One of the most common ways of representing causal structure is through a directed acyclic graph (DAG), where a directed edge between two variables in the DAG represents a direct causal effect and a directed path indicates an indirect causal effect <cit.>. 





Causal structure learning is intrinsically hard, since a DAG is generally only identifiable up to its Markov equivalence class (MEC) <cit.>. Identifiability can be improved by performing interventions <cit.>, and several algorithms have been proposed for structure learning from a combination of observational and interventional data <cit.>. Since experiments tend to be costly in practice, a natural question is how principled experimental design (i.e., selection of intervention targets) can be leveraged to maximize the performance of these algorithms under budget constraints.

Seminal works by <cit.> and <cit.> showed that experimental design can improve structure recovery in causal DAG models. 
However, these methods assume a basic framework in which experiments are performed one sample at a time. 
In practice, experimenters often perform a batch of interventions and collect samples over multiple rounds of experiments; and they must also factor in budget and feasibility constraints, such as on the number of unique interventions that can be performed in a single experiment, the number of experimental rounds, and the total number of samples to be collected. 

In genomics, for instance, genome editing technologies have enabled the collection of batches of large-scale interventional gene expression data <cit.>. 

An imminent problem is understanding how to optimally select a batch of interventions and allocate samples across these interventions, over multiple experimental rounds in a computationally tractable manner. 

Since the initial works by <cit.> and <cit.>, there have been a number of new experimental design methods under budget constraints <cit.>. These methods suffer from two drawbacks: (1) poor computational scaling <cit.> or (2) strong assumptions including the availability of infinite observational and/or interventional data from each experiment <cit.>.   

Since it is difficult to learn the correct MEC in a limited sample setting, it is desirable to use interventional samples not only to improve identifiability but also to help distinguish between observational MECs. 



Generalizing the frameworks in <cit.>, we assume the experimenter is interested in learning some function f(G) of the unknown graph G. Returning to gene regulation, one might set f(G) to indicate whether some gene X is downstream of some gene Y, i.e. if X is a descendant of Y in G. Using targeted experimental design, all statistical power is placed in learning the target function rather than being agnostic to recovering all features in the graph. In addition, we also explicitly take into account that only finitely many samples are allowed in each round, and work under various budget constraints such as a limit on the number of rounds of experimentation.

We start by reviewing causal DAGs in <ref> and then propose an entropy-based score function that generalizes the one by <cit.> and <cit.> in Section <ref>. Since optimizing this score function is in general computationally intractable, we propose our ABCD-Strategy consisting of approximations via weighted importance sampling and greedy optimization in Section <ref>. We also provide guarantees for this algorithm based on submodularity. 
Further, in contrast to earlier score functions, we show that our proposed score function is provably consistent. Finally, in <ref> we demonstrate the empirical gains of the proposed method over random sampling on both synthetic and real datasets.  



§ PRELIMINARIES




Causal DAGs:   Let G = ([p], A) be a directed acyclic graph (DAG) with vertices [p] := {1,…,p} and directed edges A, where (i,j)∈ A represents the arrow i→ j. A linear causal model is specified by a DAG G and a corresponding set of edge weights θ∈ℝ^|A|. Each node i in G is associated with a random variable X_i. Under the Markov Assumption, each variable X_i is conditionally independent of its nondescendants given its parents, which implies that the joint distribution factors as
∏_i=1^p ℙ(X_i  |  𝖯𝖺_G(X_i)), where 𝖯𝖺_G(X_i) denotes the parents of node X_i <cit.>. This factorization implies a set of conditional independence (CI) relations; 
 the Markov equivalence class (MEC) of a DAG G consists of all DAGs that share the same CI relations <cit.>. The 
 essential graph Ess(G) is a partially oriented graph that uniquely represents the MEC of a DAG by placing directed arrows on edges consistent across the equivalence class and leaves the other edges undirected <cit.>. 


Learning with Interventions:   Let  intervention I ⊆ [p] be a set of intervention targets. Intervening on I  removes the incoming edges to the random variables X_I (X_i)_i ∈ I in G and sets the joint distribution of X_I to a new interventional distribution ^I. The resulting  mutilated graph is denoted by G^I. 
A typical choice of ^I is the product distribution ∏_i ∈ I f_i (X_i), where each f_i(X_i) is the probability density function for the intervention at X_i. We denote by ^* {I_1, ⋯, I_K } the set of all K ∈ℕ allowed interventions and by ⊆^* the subset of selected interventions. An intervention I = ∅ indicates observational data. We assume that ^* is a conservative family of interventions, i.e., for any i ∈ [p], there exists some I_j ∈^* such that i ∉ I_j <cit.>. Given a conservative family of targets , two DAGs G_1 and G_2 are -Markov equivalent if they are observationally Markov equivalent and for all I ∈, G_1^I and G_2^I have the same skeleta <cit.>. The set of -Markov equivalent DAGs can be represented by the -essential graph Ess^(G), a partially directed graph with at least as many directed arrows as Ess(G) <cit.>.


Bayesian Inference over DAGs:   In various applications, the goal is to recover a function f(G) of the underlying causal DAG G given a mix of n independent observational and interventional samples D = {(X_mi, I^(m)): I^(m)∈^*, m∈ [n], i∈ [p]}. For example, we might ask whether an undirected edge (i,j) is in A, or we might wish to discover which nodes are the parents of a node i. We can encode our prior structural knowledge about the underlying DAG through a prior ℙ(G).
The likelihood (D | G) is obtained by marginalizing out θ:


    (D | G)    = ∫_θ( D, θ| G)   dθ
       = ∫_θ(D |θ, G) (θ| G)   dθ


and can be computed in closed-form for certain distributions <cit.>. Applying Bayes' Theorem yields the posterior distribution (G | D) ∝(D | G) (G), which describes the state of knowledge about G after observing the data D. Given the posterior, we can then compute _(G | D) f(G), the posterior mean of some target function f(G). Note that when f is an indicator function, this quantity is a posterior probability.





§ OPTIMAL BAYESIAN EXPERIMENTAL DESIGN


Our goal is to learn some feature f(G) of the unknown graph through experimental design under budget constraints such as limited number of experimental rounds. In principle, this question can be answered using optimal Bayesian experimental design, namely 
by selecting the experiment that maximizes the expected value of some utility function U, where the expectation is with respect to hypothetical data generated according to our current beliefs <cit.>. Here, the expected utility function U is a function defined on multisets of ^*:     

 
The expected utility U^f(ξ; D) of a multiset of interventions ξ∈ℤ^^* for learning a function f(G) given currently collected data D is given by


    U^f(ξ; D)    = _y ∼(y | D, ξ) U^f(y, ξ; D) 
       = _G, θ| D_y | G, θ, ξ U^f(y, ξ; D),    y ∈^|ξ|,


where U^f(y, ξ; D) ∈ is a function measuring the utility of observing additional samples y from a proposed design ξ and |ξ| ∑_I ∈^* |# times I in ξ|. The optimal Bayesian design ξ^* under a set of design constraints C is given by


    ξ^* ∈ξ∈ℤ^^*∩  C U^f(ξ; D).


We denote samples collected from such an optimal strategy ξ^* by D_ξ^* 


In <ref>, y is distributed according to our current beliefs (y | D, ξ) = _G, θ| D[ (y | G, θ, ξ) ], a mixture distribution over (G, θ), and the utility function U^f(y, ξ; D) is averaged over this distribution. There are many potential choices for U^f(y, ξ; D), a popular one being mutual information. <cit.>, <cit.> and <cit.> propose optimizing mutual information for the problem of recovering the full graph. More precisely, they consider the problem where f(G) = G in the active learning setting, where the experimenter can adaptively collect one sample at a time. 
We here extend their framework to general functions f(G) and the 
batched setting, where multiple samples are collected at once and the total number of batches is fixed by the experimenter. Hence, U^f must be defined on multisets instead of elements of ^* since multiple samples (i.e., interventions of the same type) may be collected in each batch. Note that the difficulty in solving <ref> stems from the constraint set C, which renders this optimization problem combinatorial.  




Recently, <cit.> proposed a Bayesian experimental design method to work in the batched setting. The authors proposed a utility function based on the expected number of additional edges that could be oriented by performing a particular intervention given the observational MEC. This function is similar to the one proposed by <cit.> and <cit.>, in which interventions are chosen that fully identify the causal network given the MEC. Unfortunately, the algorithm in <cit.> has factorial dependence on the size of the batch; in addition, we prove in <ref> that their proposed utility function is, in general, not consistent; see <ref> for a definition of  consistency.  

We therefore follow the approach taken by <cit.> and <cit.> and consider the utility function U^f(y, ξ; D) to be given by mutual information. Maximizing the mutual information is equivalent to picking the set of interventions that leads to the greatest expected decrease in entropy of f(G). The mutual information utility function is given by 


    U_M.I.^f(y, ξ; D)  H(f | D) - H(f | D, y=y, ξ),


where the entropy H(f | D) equals


    ∑_e: f(G) = e-(f(G) = e | D) log(f(G) = e | D), 
     and (f(G) = e | D) = _(G | D)(f(G) = e), 
    (G | D) ∝∫_θ(D | G, θ) (θ| G) (G).


To better understand the behavior of U_M.I.^f, we prove the following proposition, which  highlights the behavior of U_M.I.^f in the limit of infinite samples per intervention; this is the setting studied by <cit.> and <cit.>. 

 
Suppose that the Markov equivalence class 𝒢 of G^* is known and the goal is to identify the underlying true DAG G^*. 
Furthermore, assume a uniform prior over 𝒢, infinite samples per intervention I ∈, and at most K unique interventions per batch as in <cit.>. Then, U_M.I. selects the interventions 


    _M.I.∈|| ≤ K 1/|𝒢|∑_G ∈𝒢log_2 |Ess^(G)|








where |Ess^(G)|  |{G^'∈𝒢: G^'∈Ess^(G) }|. 
















This result (proof in Appendix) shows that in the limiting case, mutual information selects interventions that lead to the finest expected log -MEC sizes. This limiting behavior of mutual information parallels what graph-based score functions do, such as the ones considered by <cit.>, <cit.> and <cit.>, that invoke the Meek Rules <cit.> to select interventions that orient the most number of edges in the -essential graphs (in expectation).   




A score function based on mutual information is particularly appealing since it not only has desirable properties in the infinite sample setting, but also does not require the MEC to be known, naturally handling the case of finite sample sizes. In particular, a score function based solely on Meek rules will not pick the same intervention twice by definition, since repeating the same intervention does not improve identifiability. As a result, adapting graph-based score functions in the finite sample regime requires first constructing an intervention set and then allocating samples instead of jointly picking and allocating samples. Mutual information, on the other hand, can pick the same intervention twice; for example if a particular intervention is very informative, selecting it twice and allocating more samples to it might lead to a greater expected decrease in entropy than a new intervention. 





 §.§ Budget Constraints



So far we have not specified the constraint set C in <ref>. To this end, we assume that the experimenter has a total of N samples to allocate across B batches. While one could try to optimize the partition of N samples across batches, in this work we study the simpler case where each batch b, 1 ≤ b ≤ B, receives a pre-specified amount of samples N_b with ∑_bN_b=N. For simplifying notation we assume throughout that N_b = N/B. 
We leave the study of adaptive batch sizes N_b for future work. The constraint set then equals,


    C_N, b   {ξ∈ℤ^^*: |ξ| = N_b },


where the subscripts on C emphasize the dependence on N and b.
Then, the optimal design in batch b is obtained by solving the following combinatorial optimization problem:


    ξ^*_b ∈ξ∈ℤ^^*∩ C_N, b U(ξ; D_b-1),


where D_b-1 [D_ξ^*_1, ⋯, D_ξ^*_b-1] is all the data collected at the start of batch b and U(ξ; D_b-1) could, for example, be the mutual information defined in <ref>. Notice that while a particular form of U(ξ; D_b-1) is provided in <ref>, U(ξ; D_b-1) need not necessarily be a Bayesian utility function to fit within the framework of <ref>.  

We now define a natural notion of consistency for any experimental design method that can be cast as an optimization routine in the form of <ref>. Since the consistency of a utility function should not depend on a specific constraint set such as C_N, b, <ref> assumes the constraint set is arbitrary and set by the practitioner.  

  
Suppose f(G) is identifiable in Ess^^*(G^*), where G^* is the true unknown DAG. Let C_N, b, 1 ≤ b ≤ B denote the constraints in batch b. 
A utility function U(ξ) is budgeted batch consistent for learning a target feature f(G) if


    ( f(G) | D_B) (f(G) = f(G^*)),


as N, B →∞, where μ^* is the law determined by the true unknown causal DAG (G^*, θ^*) 








 
U_M.I.^f is budgeted batch consistent for single-node interventions, i.e., when ^* = {{1}, ⋯, { p}}.



While Theorem <ref> may not be surprising (proof in the Appendix), we found that various utility functions that seem natural and have been proposed in earlier work are not consistent in the budgeted setting. 
In particular, in <ref> we show that the utility function proposed by <cit.> is not consistent for single-node interventions. The main issue is that there are DAGs and constraint sets for which the same interventions keep getting selected, instead of selecting new interventions to fully identify f(G).    





§ TRACTABLE ALGORITHM
 


While <ref> provides a general framework for targeted experimental design, there are several computational challenges that we have not yet addressed. The first challenge is computing U_M.I.^f(ξ; D). This objective function requires summing over an exponential number of DAGs and marginalizing out the edge weights θ. In this section, we discuss how to approximate U_M.I.^f(ξ; D) by sampling graphs (either through MCMC or the DAG-bootstrap <cit.>) and using the maximum likelihood estimator of θ for each graph. Taken together, these approximations not only allow the mutual information score to be computed tractably but also lead to desirable optimization properties. In particular, we prove in <ref> that our approximate utility function is submodular. This property enables optimizing the approximate objective in a  sequential greedy fashion with provable guarantees on optimization quality.  




 §.§ Expectation over (G,θ)



A serious problem from a computational perspective is the expectation over (G, θ) in <ref>. Since the number of DAGs grows superexponentially with p, enumerating all possible DAGs is intractable. Instead, in each batch b, we propose to sample T graphs according to the posterior (G | D_b-1). This can be done using a variety of different Markov chain Monte-Carlo (MCMC) samplers; see for example <cit.>. An alternative that is often faster but still achieves good performance, is approximating the posterior via a high-probability candidate set of T DAGs Ĝ_T <cit.>. While there are many ways to build up this set, a popular approach is through the nonparametric DAG bootstrap <cit.>. The main idea is to subsample the data (with replacement) T times and fit a DAG learning algorithm to each of the generated datasets to construct Ĝ_T. Each G ∈Ĝ_T can then be weighted according to the ratio of unnormalized posterior probabilities,


    w_G, D(G) (D | G)/∑_G ∈Ĝ_T(G) (D | G)


to form an approximate posterior (G)  w_G, D(G ∈Ĝ_T). The DAG learning algorithm used for this purpose must be able to handle  a mix of observational and interventional data. Two recent methods that have been developed for this purpose are given in  <cit.> and <cit.>. We summarize constructing an approximate posterior via the DAG bootstrap in <ref>.




Given Ĝ_T, which can be constructed from <ref> or sampled from a Markov chain, we next discuss how to compute the expectation over θ. Recall that U_M.I.^f(ξ; D) is given by


    _G | D[ _y | G, ξ U^f_M.I.(y, ξ; D) ] 
       = _G | D[ _θ| G, D_y | G, θ, ξ U^f_M.I.(y, ξ; D) ].


Instead of carrying out the expensive expectation over θ| G, D in <ref>, we use the MLE of θ for each sampled G. This approximation is justified by the Bernstein-von Mises Theorem, which implies that  




    (θ| G, D)    → N(, 1/n I(θ_G)^-1),
       θ (D | G, θ).


Here, n is the number of datapoints in D, and I(θ_G) is the Fisher information matrix of the parameter θ_G, which is the asymptotic limit of the maximum likelihood estimator  <cit.>. Therefore, the posterior distribution θ| D, G concentrates around  at the standard O(1/√(n)) statistical rate. Hence, for moderate n (e.g., when a moderate amount of observational data is provided at the start of the experimental design), <ref> implies   


    _G | D_θ| G, D[ _y | G, θ, ξ U^f_M.I.(y, ξ; D) ] 
       ≈_G | D[ _y | G, , ξ U^f_M.I.(y, ξ; D) ].


In <cit.>, the authors provide a closed-form expression for  when y | G, θ is multivariate Gaussian. In this case,  is a simple function of the sample covariance matrix. 



 §.§ Approximating Mutual Information

While in the previous subsection we showed how to approximate the expectations in <ref>, computing U_M.I.^f(y, ξ) even for a fixed y is intractable since we must sum over all possible DAGs. Recall from <ref> that the mutual information utility function is    


    U_M.I.^f(y, ξ; D) = H(f | D) - H(f | D, y=y, ξ).


Note that H(f | D) is a constant and does not matter in the optimization over ξ. More care is required for  computing the second term in <ref>, since the posterior of G changes as a result of observing y, the realizations of the interventions specified by ξ. We therefore cannot immediately use the samples in Ĝ_T to approximate this term. To overcome this problem, we propose to use weighted importance sampling and approximate H(f | D, y, ξ) by a weighted average of DAGs in Ĝ_T. We define the importance sample weights for DAG G_i, 1 ≤ i ≤ T, by





    w_i (D, y | G_i, ξ)/(D | G_i).


In general, w_i is not equal to (y | G, ξ) since D and y are  dependent without conditioning on θ. 








While (D, y | G, ξ) can be computed in closed-form if the prior on θ| G belongs to one of the families described in <cit.>, the dependence on previous samples in the importance weights makes greedily building up the intervention set ξ expensive. In particular, since <ref> does not factorize, the importance weights must be recomputed with every new additional intervention, which again requires an integration over all parameters. 
Motivated by the approximation in <ref>, where the parameters of each sampled G ∈Ĝ_T are not marginalized out, we instead propose using the importance sample weights   


    ŵ_i    (D, y | G_i, θ̂^G_i_MLE, ξ)/(D | G_i, θ̂^G_i_MLE )
       = (y | G_i, ξ, θ̂^G_i_MLE);






ŵ_i has the natural interpretation of re-weighting each DAG  by the likelihood of the newly observed data y. 

Recall from <ref>, that U_M.I.^f(y, ξ; D) is based on weighting each DAG according to its posterior probability (G | D) ∝∫_θ(D | G, θ) (θ| G) (G). Using the importance sample weights ŵ_i translates into approximating the mutual information against a different posterior distribution in <ref>, namely 


    (G | D) ∝(D | G, ) (G),


which is a specific instance of an empirical Bayes approximation. In what follows, we denote the mutual information score based on the posterior in <ref> by Ũ_M.I.^f(y, ξ; D). 
















 §.§ Greedy Optimization


The cardinality constraint |ξ| = N_b makes our optimization problem a difficult integer program. In the following, we show how to overcome this final computational hurdle using a generalized notion of submodularity for multisets <cit.>. In particular, we prove that greedily selecting interventions provides a (1 - 1/e) guarantee on optimization quality. 



 Suppose f(G) = G i.e. the goal is to recover the full graph as in <cit.>. Then the difference between the global optimum


    v^*_b = max_ξ∈ℤ^^*∩ C_N, b _G | D_b-1_y | G, , ξ Ũ_M.I.^f(y, ξ; D)


and ṽ_b = (Ũ_M.I.^f, N_b, ^*), the output of <ref> in batch b,
satisfies ṽ_b ≥ (1 - 1/e)v^*_b, where C_N, b is defined as in <ref>.



We conjecture that <ref> holds for arbitrary functions f, but we currently only have a proof (see Appendix) for the case when f(G) = G. 


We conclude this section by summarizing the developed Active Budgeted Causal Design Strategy (ABCD-Strategy) in <ref> and then summarizing all the proposed approximations.  



In terms of approximations, <ref> implies

    _G, θ| D[_y | G, θ, ξŨ^f_M.I.(y, ξ; D) ] 
    
       ≈_G, | D[ _y | G, , ξ Ũ_M.I.(y, ξ; D) ] 
       ≈∑_t=1^T ∑_m=1^M Ũ^f_M.I.(y_tm, ξ; D), 
           s.t.  y_tmi.i.d∼ y | G_t, , ξ
       ≈∑_t=1^T ∑_m=1^M Û^f_M.I.(y_tm, ξ; D),  where



    Û^f_M.I.(y_tm, ξ; D)  H_1(f | D) - H_2(f | D), 
       _1(G | D)  w_G, D(G ∈Ĝ_T), 
       _2(G | D, y, ξ) w_G, D(y | G, ξ, )/∑_t=1^T w_G_t, D(y | G_t, ξ, θ̂^G_t_MLE),

where M is the number of synthetic datasets generated, H_1 and H_2 are the entropies induced by _1 and _2 respectively, and w_G, D is defined in <ref>. Note that Û^f_M.I. is based on the importance sample weights given in <ref>.
 
The total runtime of <ref> with input utility function Û^f_M.I. is O(pTκ^3 + |^*| MT^2N_bκ p), where κ is the maximum indegree of a graph in Ĝ_T. 


See <ref> for the proof of <ref>.


§ EXPERIMENTS
 

We begin by considering a simple case to demonstrate the behavior of our ABCD-strategy under easily interpretable conditions. Consider the chain graph on 2m -1 nodes,

    1 → 2 →…→ m →…→ p=2m-1.

The corresponding essential graph is completely undirected, and the MEC has 2m - 1 members, one with each node as the source. Assume that sufficient observational data is available to identify the MEC, and we are interested in fully identifying the DAG. Then, our ABCD-strategy selects interventions in order to minimize the expected entropy of the posterior over this MEC. Given a limit of one intervention per batch but infinite samples per batch, <ref> implies the expected entropy after intervening at node i or 2m - i, 1 ≤ i ≤ m, is 


    1/2m-1( ∑_j < ilog (i-2) + ∑_j > ilog (m-(i+2)) ),




which is minimized by choosing the midpoint i = m. Analogously, we see that the updated {∅, {m}}-MEC is of the same form, so in the second batch, the optimal intervention will be halfway through the remaining nodes. This process of bisection is illustrated in Figure <ref> and matches the behavior of our algorithm even in the finite-sample regime as described next.

Figure <ref> illustrates the performance of our ABCD-strategy on fifty 11-node chain graphs with random edge weights sampled from [-1, -.25] ∪ [.25, 1]. For comparison, we consider a random intervention strategy that uniformly distributes the samples in each batch to k interventions picked uniformly at random, where k is the maximum number of unique interventions allowed per batch. Whereas the median-performing random strategy barely reduces the entropy, the ABCD-strategy reduces the entropy significantly in all runs. When all k interventions are picked for the same batch, so that ABCD receives no feedback, the median-performing run of active learning still reduces the entropy as much as the best-performing runs of the random strategy. 



Having demonstrated the behavior of ABCD for a simple case, we now analyze the performance of our method on more general DAGs. The skeleton of each graph is sampled from an Erdös-Rényi model with density ρ=0.25. The edges of these graphs are directed by sampling a permutation of the nodes uniformly at random and orienting the edges accordingly. To avoid long runtimes when enumerating the MEC, we disposed of graphs with more than 100 members in their MEC.[From a sample of 10,000 graphs, only 54 had MEC size greater than 100. Based on the results by <cit.>, we expect the MECs to be typically small.] When the MEC is known, we may define a variant of the random strategy, Chordal-Random, which only intervenes on nodes that are in chordal components of the essential graph, i.e., nodes adjacent to at least one undirected edge. Since the Meek rules can only propagate by intervening within chordal components, Chordal-Random is a more fair baseline strategy for comparison than simple random sampling.   





Figure <ref> demonstrates the improvement in selecting interventions using the ABCD-strategy as compared to Chordal-Random when the number of unique interventions per batch is bounded by one. The entropy reduction for an interventional data set D_ξ is defined as H(G) - H(G | D_ξ)/H(G), and it is used as a metric so that MECs of different sizes are comparable. Since the number of total possible unique interventions is kB, an increase in the number of batches also increases the variability of the interventions, reflected in the increase of entropy reduction with batch size. Already with only 192 samples and 3 total batches, our ABCD-strategy is able to learn most graphs with complete certainty. The comparable performance of the Budgeted Experiment Design (BED) strategy <cit.> suggests that for the given experimental setup, the interventions that orient the most edges correspond well to those that most reduce entropy as we discussed in <ref>. Figure <ref> shows that the performance of the ABCD-strategy remains strong even when the MEC of the graph is not known. Specifically, up to 3 additional MECs were generated by randomly flipping non-covered edges that did not create cycles, and again only graphs for which the union of these MECs had cardinality less than 100 were kept. Note that we are not able to compare with BED since BED requires that the MEC is known.




DREAM4 Synthetic Dataset. Finally, we applied our experimental design strategy to gene expression data from the DREAM4 10-node in-silico network reconstruction challenge <cit.>. These data are generated from stochastic differential equations and simulate microarray data of gene regulatory networks. We constructed an observational dataset from the wild-type, multifactorial perturbation, and time 0 time-series samples (16 samples in total), and similarly, interventional datasets from the knockdown and knockout samples (2 samples each).

Previous work on experimental design 
applied to biological datasets <cit.> has focused on learning the entire network. In practice, practitioners may be specifically interested in performing experiments to elucidate a functional of the network, such as the pathway or local network surrounding a gene of interest. To emulate this setting, we applied our ABCD-strategy towards learning the downstream genes of select genes from the true network (Figure <ref>, top). Despite high variations in learning due to the small size of the dataset, we observed an improvement over the random strategy for several central genes (<ref>, bottom; <ref>). These results illustrate the promise of applying targeted experimental design for applications to genomics.



§ CONCLUDING REMARKS


We proposed Active Budgeted Causal Design Strategy (ABCD-Strategy), an experimental method based on optimal Bayesian experimental design with provable guarantees on approximation quality. Empirically, we demonstrated that ABCD yields considerable boosts over random sampling for both targeted and full causal structure discovery. 
Such experimental design strategies are particularly relevant for applications to genomics, where the number of possible experiments is huge due to the possibility of intervening on combinations of genes.




§ ACKNOWLEDGEMENTS


R. Agrawal was partially supported by IBM. K.D. Yang was supported by an NSF graduate fellowship and ONR (N00014-18-1-2765).  C. Uhler was partially supported by NSF (DMS-1651995), ONR (N00014-17-1-2147 and N00014-18-1-2765), IBM, and a Sloan Fellowship. 

abbrvnat 








§ PROOFS




 §.§ Proof of <ref>

Given infinite samples per intervention I ∈, G^* is recovered up to its -Markov equivalence class. Hence, the resulting entropy after placing an infinite number of samples at each intervention is equal to log_2 |Ess^(G)| when the true DAG is G. Since the true DAG is unknown, this entropy must be averaged over our prior distribution on 𝒢, which is uniform. Hence, the entropy after observing an infinite number of samples per intervention in  equals 1/|𝒢|∑_G ∈𝒢log_2 |Ess^(G)|. Minimizing this entropy over all possible interventions sets of size at most K completes the proof. 



 §.§ Proof of <ref>


Let


    ^∞{I ∈^*: ∑_b=1^∞ |Ĩ∈ξ_b: Ĩ = I | = ∞ μ^* a.s. },


where ξ_b denotes the interventions selected at batch b by U_M.I^f. Since ^* is finite, ^∞ is non-empty. When |^∞| > 1, ^∞ is a conservative family of targets since ^* is a family of single-node interventions. Hence, we identify the ^∞-MEC of G^* in the limit of an infinite number of batches and samples <cit.>. Assume |^∞| > 1. If f(G) is identifiable in Ess^^∞(G^*), then 


    ( f(G) | D_B) (f(G) = f(G^*)).


Hence, it suffices to show that the interventions U_M.I^f selects infinitely often identifies f(G) in the limiting interventional essential graph Ess^^∞(G^*). Suppose towards a contradiction that f(G) were not fully identifiable in Ess^^∞(G^*). 
By definition of almost sure convergence, there exists some b^* < ∞ such that any Ĩ∈^* ∖^∞ is never selected again after batch b^* with probability one since ^* is finite. Maximizing U^f_M.I. is equivalent to minimizing the conditional entropy,


    H_ξ^b(f | Y_ξ) _y ∼(y | D_b, ξ) H(f | D_b, Y=y).


If b > b^*, then  


    ξ∈ℤ^^*∩  C_b H_ξ^b(f | Y_ξ) = ξ∈ℤ^^∞∩  C_b H_ξ^b(f | Y_ξ)


since any batch b after b^* never selects an intervention in Ĩ∈^* ∖^∞. Since f is not identifiable in Ess^^∞(G^*), that implies


    lim_b →∞ H_ξ_∞^b(f | Y_ξ_∞) → L > 0.


Since ^* consists of all single-node interventions, ^* can identify f(G) <cit.>. Hence, there must be some Ĩ∈^* ∖^∞ and ϵ > 0 such that 


    lim_b →∞ H_ξ_∞∪Ĩ_∞^b(f) < L - ϵ,

where Ĩ_∞ denotes selecting Ĩ infinitely many times. But <ref> implies that there must exist some batch b > b^* such that the conditional entropy of the design ξ̃ = {Ĩ} is uniformly smaller than the conditional entropy of any ξ∈ℤ^^∞. But this is a contradiction because then Ĩ would be selected again after some batch b > b^* and <ref> would no longer hold. 


















For |^∞| = 1, we no longer have a conservative family of targets. However, a nearly identical argument works by noting that, in the limit, we learn the observational equivalence class of the ^∞ mutilated graph of G^*.



 §.§ Consistency Counterexample
 



Suppose we know the Markov equivalence class of G^* and the goal is to fully recover G^*. Suppose C_b = {ξ: ξ_0 = K }, where ·_0 counts the number of unique interventions in ξ. Since there is no constraint on the number of samples, only on the number of unique interventions, we may allocate an infinite number of samples per intervention within each batch. This constraint is equivalent to the one examined in <cit.>. The scores in both <cit.> and <cit.> select interventions by maximizing the expected number of oriented edges in the interventional Markov equivalence classes. In particular, the utility function in <cit.> is equivalent to maximizing,


    U(; D) = ∑_G ∈𝒢 A(Ess^(G)) (G),



where A(Ess^(G)) equals the additional number of edges oriented relative to the observational Markov equivalence class. Suppose G^* equals  the graph in <ref> and that K=1 unique interventions are allowed within each batch. Assume that ^* = {{1}, ⋯, {4}} and that we start with a uniform prior over 𝒢. Then, since all arrows are undirected in the observational Markov equivalence class, symmetry implies U({j}; ∅) = U({j}; ∅) for all i, j ∈ 1, ⋯, 4. Without any loss of generality suppose intervention one is selected in batch one. We show that every subsequent batch will select intervention {1}. If only {1} were selected, U(; D) would not be consistent since the {∅, {1 }}-MEC(G^*) contains two graphs, as shown at the bottom of <ref>. After batch one, the posterior is supported on these two graphs since an infinite number of samples are allocated to the intervention at node one. 

The utility function in <ref> scores interventions relative to the observational equivalence class, which causes the consistency issue. In particular, the posterior in batch two is only supported on the two DAGs given in the bottom box of <ref>.  The score of {1} equals 5 in batch two while the scores of interventions {2}, {3}, {4} equal 4, 3, 4, respectively. Hence, in batch two, intervention {1} will be selected again, but the posterior will remain the same since the {∅, {1 }} interventional Markov equivalence class of G^* is already known.

An easy way to fix <ref> (for this given counterexample) would be to only select interventions not selected in previous batches. This modification would fix the issue with the counterexample, namely prevent intervention one from being selecting infinitely often. However, when one can only allocate a finite number of samples per batch, this modification would not lead to a consistent estimator. In particular, if a certain intervention is done in some batch, and that intervention must be conducted in order to identify f, then only placing finitely many samples to that intervention in that batch and never placing any more samples in subsequent batches will not lead to a consistent method. 

























































































































































 §.§ Proof of <ref>









 <cit.>
Let E be a finite set. A function f: ℤ^E → is diminishing returns submodular (DR-submodular) if for x ≤ y 


    f(x + χ_e) - f(x) ≥ f(y + χ_e) - f(y),  x, y ∈ℤ^E


where e ∈ E and χ_e is the ith unit vector.  


 
Ũ_M.I.^f(ξ; D) is DR-submodular.


f(G) = G so we omit f in Ũ_M.I.^f to simplify notation. Since the sum of submodular functions is submodular, it suffices to show 


    _y | G, , ξ Ũ_M.I.(y, ξ; D)     = H(G) - H(G | Y_ξ) 
        = I((G, ), Y_ξ)


is DR-submodular, where I is the mutual information. Consider an A ⊆ B ∈ℤ^^*. Take any C ∈^*. Since entropy decreases with more conditioning,  


    H(Y_C | Y_A) - H(Y_C | (G, )) ≥
     H(Y_C | Y_B) - H(Y_C | (G, )).


By conditional independence, 

    H(Y_C | (G, ))     = H(Y_C | (G, ), Y_A) 
       =  H(Y_C | (G, ), Y_B).


Hence, <ref> may be rewritten as,


    I((G, ), Y_C | Y_A) = 
     H(Y_C | Y_A) - H(Y_C | (G, ), Y_A) ≥
     H(Y_C | Y_B) - H(Y_C | (G, ), Y_B) = 
    
            I((G, ), Y_C | Y_B).


<ref> implies 


    I((G, ), Y_A + Y_C) - I((G, ), Y_A) 
       ≥ I((G, ), Y_B + Y_C) - I((G, ), Y_B)


as desired.


The proof of <ref> then follows directly from <ref> and <cit.>.



 §.§ Proof of <ref>
 
For each graph G ∈𝒢_T, compute the associated edge weights . Computing each  takes O(pκ^3) time using the formula given in <cit.>. Since there are T DAGs, the total time to compute the MLE estimates of the edge weights of each DAG is O(Tpκ^3). 
Sampling from a multivariate Gaussian with bounded indegree with known adjacency matrix takes O(pκ) time. Û^f_M.I. requires a total of |^*|MN_bT^2 samples. Hence, the total computation time of sampling all the y_mt in <ref> is O(|^*|MN_b κ pT^2). Evaluating Û^f_M.I. takes O(M T^2) time using these samples, which is of lower computational complexity than computing Û^f_M.I.. Hence, the total runtime is O(pκ^3 + |^*|MN_b κ pT^2). 



 §.§ Constraint on the Number of Unique Interventions
 

If we are only allowed to allocate at most K unique interventions per batch, we modify <ref> by allocating N_b/K samples per intervention in <ref>. 
Once an intervention is selected, that intervention is removed from I^* and another one is greedily selected from the remaining set. With this strategy, <ref> will terminate after K iterations. Hence, there will be at most K unique interventions as desired. 



 §.§ DREAM4 Supplementary Figures

We applied our targeted experimental design strategy towards learning the downstream pathways of select genes from a 10-node network from the DREAM4 challenge. We observed a modest improvement over the random strategy for some central genes in the network (<ref>, top). However, the results are subject to high variations (<ref>, bottom), which we surmise to be due to the small size of the observational dataset. Nevertheless, these preliminary results illustrate the promise of applying targeted experimental design to real, large-scale biological datasets.



