Near-Optimal Multi-Perturbation Experimental Design for Causal Structure Learning
Year: 2021

Abstract:
Causal structure learning is a key problem in many domains. Causal structures can be learnt by performing experiments on the system of interest. We address the largely unexplored problem of designing a batch of experiments that each simultaneously intervene on multiple variables. While potentially more informative than the commonly considered single-variable interventions, selecting such interventions is algorithmically much more challenging, due to the doubly-exponential combinatorial search space over sets of composite interventions. In this paper, we develop efficient algorithms for optimizing different objective functions quantifying the informativeness of a budget-constrained batch of experiments. By establishing novel submodularity properties of these objectives, we provide approximation guarantees for our algorithms. Our algorithms empirically perform superior to both random interventions and algorithms that only select single-variable interventions.

heading: 1 Introduction
text: The problem of finding the causal relationships between a set of variables is ubiquitous throughout the sciences. For example, scientists are interested in reconstructing gene regulatory networks (GRNs) of biological cells [11]. Directed Acyclic Graphs (DAGs) are a natural way to represent causal structures, with a directed edge from variable X to Y representing X being a direct cause of Y [33].
Learning the causal structure of a set of variables is fundamentally difficult. With only observational data, in general we can only identify the true DAG up to a set of DAGs called its Markov Equivalence Class (MEC) [35]. Empirically, for sparse DAGs the size of the MEC grows exponentially in the number of nodes [17]. Identifiability can be improved by intervening on variables, meaning one perturbs a subset of the variables and then observes more samples from the system [8, 16, 40]. There exist various inference algorithms for learning causal structures from a combination of observational and interventional data [16, 36, 40, 28, 34]. Here we focus on the identification of DAGs that have no unobserved confounding variables.
Performing experiments is often expensive, however. Thus, we are interested in learning as much about the causal structure as possible given some constraints on the interventions. In this work, we focus on the batched setting, where several interventions are performed in parallel. This is a natural setting in scientific domains like reconstructing GRNs. Existing works propose meaningful objective
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
functions for this batched causal structure learning problem and then give algorithms that have provable guarantees [3, 13]. However, these works focus on the setting where only a single random variable is perturbed per intervention. It is an open question as to whether there exist efficient algorithms for the multiple-perturbation setting, where more than one variable is perturbed in each intervention.
For the task of reconstructing GRNs, it is now possible for experimenters to perturb multiple genes in a single cell [2, 7]. Figures 1 b) and c) illustrate a specific example where a two-node intervention completely identifies a DAG in half as many interventions as single-node interventions. In general, it is possible for a set of q-node interventions to orient up to q-times more edges in a DAG than single-node interventions (see the supplementary material for a more general example). While multiperturbation interventions can be more informative, designing them is algorithmically challenging because it leads to an exponentially larger search space: any algorithm must now select a set of sets.
Our main contribution is to provide efficient algorithms for different objective functions with accompanying performance bounds. We demonstrate empirically on both synthetic and GRN graphs that our algorithms result in greater identifiability than existing approaches that do not make use of multiple perturbations [13, 40], as well as a random strategy.
We begin by introducing the notation and the objective functions considered in this work in Section 2, before reviewing related work in Section 3. In Section 4 we present our algorithms along with proofs of their performance guarantees. Finally, in Section 5 we demonstrate the superior empirical performance of our method over existing baselines on both synthetic networks and on data generated from models of real GRNs.

heading: 2 Background and Problem Statement
text: Causal DAGs Consider a causal DAG G = ([p], E) where [p] := {1, ..., p} is a set of nodes and E is a set of directed edges. Let (i, j) ∈ E iff there is an edge from node i to node j. Each node i is associated with a random variable Xi. In the GRN example, Xi would be the measurement of the gene expression level for gene i. An edge from i → j would represent gene i having a causal effect on the expression of gene j. The functional dependence of a random variable on its parents can be described by a structural equation model (SEM).
The probability distribution over X = (X1, . . . , Xp) is related to G by the Markov property, meaning each variable Xi is conditionally independent of its non-descendants given its parents [33]. From conditional independence tests one can determine the MEC of G, a set of DAGs with the same conditional independancies between variables. All members of the MEC share the same undirected skeleton and colliders [35]. A collider is a pair of edges (i, k), (j, k) ∈ E such that (i, j), (j, i) /∈ E. The essential graph of G, Ess(G), is a partially directed graph, with directed edges where all members of the MEC share the same edge direction, and with undirected edges otherwise [5]. Ess(G) uniquely represents the MEC of G. These MECs can be large, so we seek to perform interventions on the nodes to reduce the MEC to a smaller set of possible DAGs.
Interventions We use the term intervention to refer to a set I ⊂ [p] of perturbation targets (variables). We assume all interventions are hard interventions, meaning intervening on a set I removes the incoming edges to the random variables XI := (Xi)i∈I and sets their joint distribution
to some interventional distribution PI [9]. In the GRN reconstruction example this corresponds to, for example, running an experiment where we knockout all genes in set I . Some of our results extend easily to the alternative model of soft interventions [25], as we discuss in the supplementary material.
We use I = 2[p] to refer to the set of all possible interventions. Our goal will be to select a batch ξ of interventions, where ξ is a multiset of some I ∈ I. For practical reasons, we typically have constraints on the number of interventions, i.e., |ξ| ≤ m and on the number of variables involved in each intervention |I| ≤ q, ∀I ∈ ξ. Namely, there are at most m interventions per batch and each intervention contains at most q nodes. A constraint on the number of perturbations per intervention is natural in reconstructing GRNs, since perturbing too many genes in one cell will leave it unlikely to survive. We refer to the set of ξs satisfying these constraints as Cm,q . The observational distribution (no intervention) is given by ξ = ∅. For any set of interventions ξ and DAG G, there is a set of ξ-Markov equivalent DAGs. These are the set of DAGs that have the same set of conditional independencies under all I ∈ ξ and under the observational distribution. This set of DAGs is no larger than the MEC of G and can similarly be characterized by an essential graph Essξ(G) [16].
We will always assume that there exist no unobserved common causes of any pair of nodes in G. We also assume that the distribution of the random variables satisfies faithfulness with respect to G [33].
Choosing optimal interventions We seek to maximize an objective function F that quantifies our certainty in the true DAG. In general our goal is to determine
argmax ξ∈Cm,q F (ξ).
A natural choice for F is given by Agrawal et al. [3]. They assume that there exist parameters θ that determine the functional relationships between random variables. For example, this could be the coefficients in a linear model. Given existing data D, we try to choose ξ that maximizes Objective 1 (Mutual Information (MI)).
FMI(ξ) = EG|DEy|G,θ̂,ξ [ ŨM.I(y, ξ;D) ] , (1)
where y is the set of samples from the interventions, θ̂ is the current estimate of the parameters given D and G, and ŨM.I(y, ξ;D) is the mutual information between the posterior over G and the samples y. Each intervention produces one sample in y. The use of mutual information means the objective aims to, in expectation over all observed samples and true DAGs, minimize the entropy of the posterior distribution over DAGs. There already exist a number of algorithms for determining the posterior over DAGs from observational or experimental data [40, 36, 16].
Infinite sample objectives Finding algorithms that optimize the MI objective is difficult because we have to account for noisy observations and limited samples. To remove this complexity, we study the limiting case of infinitely many samples per unique intervention. The constraints given by Cm,q still stipulate that there can be only m unique interventions, but each intervention can be performed with an infinite number of samples. We also assume that an essential graph is already known (i.e., we have infinite observational samples and infinite samples for any experiments performed so far). For objectives with infinite samples per intervention, we treat ξ as a set of interventions, not a multiset, since there is no change in objective value for choosing an intervention twice. Consider ξ′ to be the set of interventions contained in our dataset before our current batch. In this setting, maximizing Objective 1 reduces to maximizing Objective 2 (Mutual info. inf. samples (MI-∞)).
F∞(ξ) = − 1 |G| ∑ G∈G log2|Essξ∪ξ ′ (G)|, (2)
where Essξ∪ξ ′ (G) refers to the updated essential graph after performing interventions in ξ. The objective aims to, on average across possible true DAGs, minimize the log of the essential graph size after performing the interventions. The derivation of this objective is given in the supplementary material.
Ghassami et al. [13] study a different objective in the infinite-sample setting. The objective seeks to, on average across possible G given the current essential graph, orient as many edges as possible. Let
R(ξ,G, ξ′) be the set of edges oriented by ξ if the true DAG is G, and the essential graph is given by the ξ′-MEC.
Objective 3 (Edge-orientation (EO)).
FEO(ξ) = 1 |G| ∑ G∈G |R(ξ,G, ξ′)|. (3)
The function R is computed as follows. Firstly, ∀I ∈ ξ, orient undirected edge i− j in Essξ ′ (G) if i ∈ I but j /∈ I or vice-versa. Secondly, execute the Meek Rules [26], which allow inferring additional edge orientations (discussed in the supplementary material) on the resulting partially directed graph. Finally, output the set of all edges oriented. Agrawal et al. [3] show that this objective is not consistent; however, this is because they fix G to be the MEC instead of using the most up-to-date essential graph for each batch. In the supplementary material, we show that the version of the objective we work with is indeed consistent. We will drop the dependence of R on ξ′ for readability.
Below, we provide algorithms with near-optimality guarantees for Objectives 2 and 3, while motivating a practical algorithm for Objective 1.

heading: 3 Related Work
text: Causality has been widely studied in machine learning [30, 31]. Here we focus on prior research that is most relevant to our work.
Agrawal et al. [3] and Ghassami et al. [13] give near-optimal greedy algorithms for Objectives 1 and 3 respectively. Ahmaditeshnizi, Salehkaleybar, and Kiyavash [4] present a dynamic programming algorithm for an adversarial version of Objective 3, optimizing for the worst case ground truth DAG in the MEC. However, all these algorithms only apply to single-perturbation interventions. Both of these works use the submodularity of the two objectives. In this paper we address the exponentially large search space that arises when designing multi-perturbation interventions, a strictly harder problem.
Much existing work in experimental design for causal DAGs is focused on identifying the graph uniquely, while minimizing some cost associated with doing experiments [9, 18, 32, 21, 23]. When the MEC is large and the number of experiments is small, identifying the entire graph will be infeasible. Instead, one must select interventions that optimize a measure of the information gained about the causal graph.
Lindgren et al. [23] show NP-hardness for selecting an intervention set of at most m interventions, with minimum number of perturbations, that completely identifies the true DAG. This, however, does not directly imply a hardness result for our problem.
Gamella and Heinze-Deml [10] propose an approach to experimental design for causal structure learning based on invariant causal prediction. While our approach has guarantees for objectives relating to either the whole graph or functions of the oriented edges, their work is specific to the problem of learning the direct causes of one variable.
Acharya et al. [1] consider testing between two candidate causal models. However, the setting differs from ours: they assume the underlying DAG is known but allow for unobserved confounding variables.
Designing multi-perturbation interventions has been previously studied in linear cyclic networks, with a focus on parameter identification [14]. Here we focus on causal graph identification in DAGs.

heading: 4 Greedy Algorithms for Experiment Design
text: All of our algorithms follow the same general strategy. Like in previous works on single-perturbation experimental design [3, 13], we greedily add interventions to our intervention set. We add I maximizing F (ξ ∪ {I}) where ξ is the currently proposed set of interventions. This greedy selection is justified because our objectives are submodular, a property we define formally later. For singleperturbation experimental design, this is algorithmically simple since there are only p possible interventions. However, for multi-perturbation interventions even selecting greedily is intractable at scale since we have ( p q ) possible interventions. Therefore we provide ways to find an intervention
that is approximately greedy, i.e, an intervention with marginal improvement in objective that is close to that of the greedy intervention.
A further challenge with the greedy approach is that it involves evaluating the objective, which for our objectives is a potentially exponential sum over members of an essential graph. Each of the two algorithms we give has a different strategy for overcoming this.
In Section 4.1, we present DOUBLE GREEDY CONTINUOUS (DGC) for optimizing Objective 3, the edge-orientation objective. For greedily selecting interventions to maximize an exponential sum, we employ the stochastic continuous optimization technique of Hassani et al. [15].
In Section 4.2, we present SEPARATING SYSTEM GREEDY (SSG) for optimizing Objective 2, MI-∞. To greedily select interventions, we use the construction of separating systems (SS) [37, 32, 23], to create a smaller set of interventions to search over. Collectively, the interventions in the SS fully orient the graph. To handle tractably evaluating the objective, we use the idea of Ghassami et al. [13] and Agrawal et al. [3] to optimize an approximation of the objective constructed using a limited sample of DAGs.
To give near-optimality guarantees for these algorithms, we will use two properties of the objectives: monotonicity and submodularity. Definition 1. A set function F : 2V → R is monotonically increasing if for all sets I1 ⊆ I2 ⊆ V we have F (I1) ≤ F (I2). Definition 2. A set function F : 2V → R is submodular if for all sets I1 ⊆ I2 ⊆ V and all v ∈ V \I2 we have F (I1 ∪ {v})− F (I1) ≥ F (I2 ∪ {v})− F (I2).
Submodularity is a natural diminishing returns property, and many strategies have been studied for optimizing submodular objectives [22]. In both the above definitions, V is called the groundset, the set that we can choose elements from. In the single-perturbation problem, the groundset is just [p], whereas in our case it is all subsets of up to q nodes.
We show that DGC achieves an objective value within a constant factor of the optimal intervention set on Objective 3. SSG does not achieve a constant-factor guarantee, but for both infinite sample objectives we obtain a lower bound on its performance.
All of our algorithms run in polynomial time; however, they assume access to a uniform sampler across all DAGs in the essential graph. This exists for sampling from the MEC [39] but not for essential graphs given existing interventions. In practice, we find that an efficient non-uniform sampler [13] can be used to achieve strong empirical performance.

heading: 4.1 Optimizing the Edge-orientation Objective
text: In the following, we develop an algorithm for maximizing Objective 3. In fact, the algorithm we provide has a near-optimality guarantee for a more general form of FEO, namely
FEO(ξ) = ∑ G∈G a(G) ∑ e∈G w(e)1(e ∈ R(ξ,G)),
where ∀e, w(e) ≥ 0 and ∀G, a(G) ≥ 0. The weights a(G) can be thought of as corresponding to having a non-uniform prior over the DAGs in the essential graph, whilst the weights w(e) can be thought of as assigning priority to the orienting of certain edges. The inner sum above is a weighted coverage function [22] over the set of edges.
We will first show that FEO is monotone submodular over groundset I. This generalizes a result by Ghassami et al. [13] who showed the same result for groundset [p] (single perturbation interventions). Lemma 1. FEO is monotone submodular over the groundset I.
Proof. All proofs are presented in the supplementary material unless otherwise stated.
As mentioned, we cannot use a greedy search directly since the groundset I is too large. Instead, we develop an algorithm for selecting an intervention with near-maximal utility compared to the greedy choice. In particular, our strategy is to prove a submodularity result over the function F with modified domain. Consider the set function F ξEO(I) = FEO(ξ ∪ {I}) for fixed ξ.
Lemma 2. F ξEO is non-monotone submodular over the groundset [p].
The Non-monotone Stochastic Continuous Greedy (NMSCG) algorithm of Mokhtari, Hassani, and Karbasi [27] can therefore be used as a subroutine to select, in expectation, an approximately greedy intervention to add to an existing intervention set. The algorithm uses a stochastic gradient-based method to optimize a continuous relaxation of our objective, and then rounds the solution to obtain a set of interventions. The continuous relaxation of F ξEO is the multilinear extension
fξEO(x) = ∑ I∈I F ξEO(I) ∏ i∈I xi ∏ i/∈I (1− xi)
with constraints ∑ i xi ≤ q, 0 ≤ xi ≤ 1 for all nodes i. The multilinear extension can be thought of as computing the expectation of F ξEO(I), when input x is a vector of independent probabilities such that xi is the probability of including node i in the intervention. The sum over I in fξEO and the sum over DAGs in F ξEO make computing the gradient of this objective intractable. Therefore, we compute an unbiased stochastic approximation of the gradient ∇fξEO(x) by uniformly sampling a DAG G from G and intervention I from the distribution specified by x. Define
f̂ξEO(I,G) = |R(ξ ∪ {I}, G)|. (4) Mokhtari, Hassani, and Karbasi [27] show that an unbiased estimate of the gradient of f(x) can be computed by sampling G and I to approximate
∂
∂xi fξEO(x) = E G,I|x
[ f̂ξEO(I,G; Ii ← 1)− f̂ ξ EO(I,G; Ii ← 0) ] , (5)
where Ii ← 0 means that if i ∈ I , remove it. The use of a stochastic gradient means that FEO can be efficiently optimized despite it being a possibly exponential sum over G. After several gradient updates we obtain a vector of probabilities x that approximately maximizes fEO. To obtain an intervention I one uses a ROUND function, for example pipage rounding which, on submodular functions, has the guarantee that E[F ξEO(ROUND(x))] = f ξ EO(x) [6].
Theorem 1 (Mokhtari, Hassani, and Karbasi [27]). Let I∗ be the maximizer of F ξEO. NMSCG with pipage rounding, after O ( p5/2/ 3 ) evaluations of R, achieves a solution I such that
E [ F ξEO(I) ] ≥ 1 e F ξEO(I ∗)− .
The original result measures runtime in terms of the number of times we approximate the gradient in Equation 5 with a single sample (in our case a single G, I tuple). From Equations 4 and 5 we can see that the number of gradient approximations is a constant factor of the number of evaluations of R. Hence, we measure runtime in terms of number of evaluations of R. The bottleneck for evaluating R is applying the Meek Rules, which can be computed in time polynomial in p [26]. Note that the NMSCG subroutine can be modified to stabilize gradient updates by the approximation of a Hessian, in which case the same guarantee can be achieved in O ( p3/2/ 2 ) [15]. We use this version of NMSCG for the experiments.
Our main result now follows from the fact that selecting interventions approximately greedily will lead to an approximation guarantee due to lemma 1.
Theorem 2. Let ξ∗ ∈ CN,b be the maximizer of Objective 3. DGC will, after O ( m4p5/2/ 3 ) evaluations of R, achieve a solution ξ such that
E[FEO(ξ)] ≥ ( 1− 1
e1/e
) FEO(ξ ∗)− .
We highlight this is a constant-factor guarantee with respect to the optimal batch of interventions. Our bound requires compute that is low order in p with no dependence on q. Without even accounting for computing the possibly exponential sum in the objective, merely enumerating all possible interventions for fixed q is O(pq). Ghassami et al. [13] give a similar constant-factor guarantee for batches of single-perturbation interventions (q = 1). Their bound is within 1− 1e of the optimal single-perturbation batch. In the
Algorithm 1 DOUBLE GREEDY CONTINUOUS(DGC)
Input: essential graph G, constraints Cm,q, objective FEO Init ξ ← ∅ while |ξ| ≤ m do I ← ROUND(NMSCG(F ξEO, q)) ξ ← ξ ∪ {I} end while Output: a set ξ of interventions
Algorithm 2 SEPARATING SYSTEM GREEDY(SSG)
Input: essential graph G, constraints Cm,q, objective F̃∞ Init I ← ∅ S ← SEPARATE(q,G) while |ξ| ≤ m do I ′ ← argmaxI F̃∞(I ∪ {I}) ξ ← ξ ∪ {I ′} end while Output: a set ξ of interventions
supplementary material, we show that the optimal multi-perturbation intervention can orient up to q times more edges than the optimal single perturbation intervention. In Section 5 we experimentally verify the value of multi-perturbation interventions by comparing DGC to the algorithm presented in Ghassami et al. [13]. If we allow for soft interventions, the 1− 1e guarantee can also be obtained for multi-perturbation interventions. See the supplementary material for details.

heading: 4.2 Optimizing the Mutual Information Objective
text: We now consider an algorithm for maximizing Objective 2. First, we note that computing the sum over G and the size of ξ ∪ ξ′-essential graphs is computationally intractable. The computation of |Essξ∪ξ ′ (G)| makes F∞ a nested sum over DAGs. Like Agrawal et al. [3], we optimize a computationally tractable approximation to F∞. First, we uniformly sample a multiset of DAGs G̃ from G to construct our approximate objective
F̃∞(ξ) = − 1
|G̃| ∑ G∈G̃ log2|Ẽss ξ∪ξ′ (G)|,
where Ẽss ξ∪ξ′ is the submultiset of G̃ consisting of elements in Essξ∪ξ ′ . A submodularity result
similar to lemma 1 can also be proven for F̃∞.
Lemma 3. F̃∞ is monotone submodular.
We first show that an approach similar to that used by DGC does not so easily give a near-optimal guarantee. Similarly to above, define F̃ ξ∞(I) = F̃∞(ξ ∪ {I}) for fixed ξ. Proposition 1. There exists G, G̃ such that F̃ ξ∞ is not submodular.
Hence we cannot use existing algorithms for submodular optimization to construct near-greedy interventions. We instead take a different approach. Suppose we can reduce I to some set of interventions S much smaller than I, such that F̃∞(S) has the maximum possible objective value. Due to lemma 3, we can obtain a guarantee by greedily selecting m interventions from S . A method for constructing the set S comes from separating system constructions. Definition 3. A q-sparse G-separating system of size N is a set of interventions S = {S1, S2..., SN} such that |Si| ≤ q and for every undirected edge (i, j) ∈ G there is an element S ∈ S such that exactly one of i, j is in S [32].
A separating system of G completely identifies the true DAG, and hence obtains the maximum possible objective value 0 without necessarily satisfying the constraints Cm,q. As an example, in Figure 1(b,c) we see 1 and 2–sparse separating systems respectively.
We will make use of an algorithm SEPARATE(q,G) which efficiently constructs a q-sparse separating system of G. Wegener [37] and Shanmugam et al. [32] give construction methods that are agnostic to the structure of G (it will identify any DAG with p nodes). Lindgren et al. [23] give a construction method that depends on the structure of G. Since the separating system obtains the maximum objective value, we can greedily select m interventions from this set and obtain a lower bound on the objective value due to submodularity (lemma 3).
Theorem 3. For q ≤ bp/2c,m ≤ |S|, SSG outputs ξ ∈ Cm,q with objective value
F̃∞(ξ) ≥ (1− m
dp/qedlog pe )F̃∞(∅)
in O(m|G̃|pq log p) evaluations of R, when using SEPARATE as in Shanmugam et al. [32].
Increasing q does not necessarily increase F̃∞(ξ). However, the bound we give becomes more favourable as q increases because the upper bound on |S| decreases. In practice, we run SSG ∀q′ ≤ q and pick the intervention set with the highest objective value on F̃∞.
Note that SSG can also be used with a similar guarantee for an analogous approximation of FEO.

heading: 5 Experiments
text: To evaluate our algorithms we consider three settings. Firstly, randomly generated DAGs, using infinite samples per intervention. Secondly, randomly generated DAGs with linear SEMs, using finite samples per intervention. Finally, subnetworks of GRN models, using infinite samples per intervention. Full details on all of our experiments can be found in the supplementary material. For code to reproduce the experiments, see https://github.com/ssethz/multi-perturbation-ed.
Infinite samples We evaluate our algorithms using Objective 3. We consider selecting a batch of experiments where only the MEC is currently known. We vary the type of random graph and the constraint set Cm,q . The following methods are compared:
• RAND: a baseline that for m interventions, independently and uniformly at random selects q nodes from those adjacent to at least one undirected edge;
• GREEDY: greedily selects a single-perturbation intervention as in Ghassami et al. [13];
• DGC: our stochastic optimization approach;
• SSG-A: our greedy approach using the graph agnostic separating systems of Shanmugam et al. [32];
• SSG-B: as above, using the graph-dependent separating system constructor of Lindgren et al. [23].
Since there are infinite samples per intervention, the exact SEM used to generate data is not relevant. We plot the mean proportion of edges identified and error bars of 1 standard deviation over 100 repeats. Noise between repeats is due to randomness in the graph structure and in the algorithms themselves.
In Figure 2 a) we display the results for Erdös-Renyí random graphs with edge density 0.1 (ER(0.1)) and 40 nodes. To prevent trivial graphs and large runtimes, we only consider graphs with MEC sizes in the range [20, 200]. The observations given here were also found for denser Erdös-Renyí graphs and tree graphs, in addition to graphs with less nodes.
For all constraint values, all the algorithms improve greatly over RAND and GREEDY. SSG-B outperforms SSG-A, likely because the graph-sensitive separating system construction tends to return a groundset of more effective interventions.
SSG-B achieves similar objective value to DGC. However, DGC behaves most robustly when the graph structure is chosen adversarially. For example, consider Figure 2 b). Here we plot the proportion of identified edges on a p = 5 fully connected graph. On this graph, the separating system construction of Lindgren et al. [23] will always return the set of all single-node interventions. Therefore, its performance does not improve with q, whilst DGC’s does. An adversarial example for SSG-A is constructed in Figure 2 c): an MEC that consists of 3 disconnected star graphs with 7, 7 and 6 nodes. In this case, DGC and SSG-B can orient most of the graph in a single intervention, whereas SSG-A will likely not contain such an intervention in the separating system it constructs.
Finite samples We use linear SEMs, with weights generated uniformly in [−1,−0.25] ∪ [0.25, 1]. Measurement noise is given by the standard normal distribution. The underlying DAGs are generated in the same way as the infinite sample experiments. Before experiment selection, we obtain 800 observational samples of the system. 3 samples are obtained for each intervention selected by our algorithms. Each perturbation fixes the value of a node to 5. We approximate Objective 1 using the methods of Agrawal et al. [3]. In particular, an initial distribution over DAGs is estimated by bootstrapping the observational data and using the techniques of Yang, Katcoff, and Uhler [40] to infer DAGs. Each DAG in the distribution is weighted proportionally to the likelihood of the observational data given the DAG and the maximum likelihood estimate of the linear model weights. The posterior over DAGs after interventions is computed by re-weighting the existing set of DAGs based on the likelihood of the combined observational and interventional data. To ensure the distribution has support near the true DAG, we include all members of the true DAG’s MEC in the initial distribution.
With finite samples, our methods do not have guarantees but can be adapted into practical algorithms:
• GREEDY: greedily optimize Objective 1 with single perturbation interventions (Agrawal et al. [3]);
• DGC-∞: optimizes Objective 3, with the summation over DAGs being a weighted sum over the DAGs in the initial distribution;
• SSG-B: optimizes FMI, greedily selecting from the separating system of Lindgren et al. [23];
• SSG-B-∞: approximates the objective using Objective 2 and optimizes with SSG.
For each algorithm, we record FMI of the selected interventions, over 200 repeats. In Figure 2 d) the objective values obtained by each algorithm are shown for varying batch size and q = 5. DGC-∞ performs worse than SSG-B and SSG-B-∞, perhaps because it is optimizing FEO which is not totally aligned with FMI. Between SSG-B and SSG-B-∞, there was a small benefit to directly optimizing FMI as opposed to its infinite sample approximation F∞. Accounting for finite samples may lead to greater improvements when there is heteroscedastic noise or a wider range of weights.
Performance on FMI corresponds closely with performance on a downstream task as shown in Figure 2 e). For each algorithm, we compute the posterior over DAGs given the selected interventions. Then, we independently predict the presence of each edge in the true DAG. Figure 2 e) plots the average F1 score for each algorithm. In finite samples, our approaches outperform both RAND and GREEDY.
DREAM 3 networks We evaluate our algorithms under infinite samples using subgraphs of GRN models. From the “DREAM 3 In Silico Network” challenge [24], we use the 5 subgraphs with p = 50 nodes. Here, we present the results for “Yeast1", which is the graph requiring the most interventions for our methods to orient. Results for the other graphs are in the supplementary material.
We compare the same algorithms as considered in the other infinite sample experiments. In Figure 2 f) we record the proportion of unknown edges that were oriented by each algorithm. Our methods and RAND all have intervention sizes of 5. For each method we perform 5 repeats. On Yeast1, our methods all perform similarly and outperform RAND and GREEDY. This is found for the other DREAM 3 graphs too, except on one subgraph where GREEDY performs similarly to our methods.

heading: 6 Conclusions
text: We presented near-optimal algorithms for causal structure learning through multi-perturbation interventions. Our results make novel use of submodularity properties and separating systems to search over a doubly exponential domain. Empirically, we demonstrated that these algorithms yield significant improvements over random interventions and state-of-the-art single-perturbation algorithms. These methods are particularly relevant in genomics applications, where causal graphs are large but multiple genes can be intervened upon simultaneously.
Acknowledgments and Disclosure of Funding
This research was supported in part by the Swiss National Science Foundation, under NCCR Automation, grant agreement 51NF40 180545. Caroline Uhler was partially supported by NSF (DMS1651995), ONR (N00014-17-1-2147 and N00014-18-1-2765), IBM, and a Simons Investigator Award.
Thank you to Raj Agrawal for a helpful discussion regarding experiments.
Experiments were performed on the Leonhard cluster managed by the HPC team at ETH Zürich

Review summaries:
The paper considers a challenging problem of designing multi-perturbation experiments for causal discovery. They focus on two objective functions and develop one algorithm for each objective. The approximation accuracy is studied by showing submodularity of the proposed algorithms. Empirical studies are used to demonstrate the superiority of the proposed algorithms over random design and single-perturbation design.

This paper tackles the problem of causal structure learning from interventional data. Specifically, they seek a series of interventional experiments which can rapidly identify the correct causal graph, or reduce some measure of uncertainty about the true causal structure. In contrast to previous work, this paper considers the more general setting in which interventions can be made on more than one variable simultaneously (within one experiment). To select an optimal sequence of subsets on which to intervene, the authors begin by discussing 3 objectives functions that could be used to select “good” sets of interventions: two based on MI, one based on orienting edges. All objectives are grounded in existing literature. The authors propose two different algorithms to actually solve the experimental design problem. First, they propose Double Greedy Continuous. As the name implies, this method is greedy both with respect to adding interventions to the experiment sequence and with respect to adding variables to intervene on in a single interventional experiment. The authors provide strong theoretical guarantees for this method with respect to their edge-orienting objective. For the MI objective, a suitable sub-modularity property does not hold for double greedy. Instead, the authors propose a new algorithm, Separating Systems Greedy. This is greedy with respect to adding interventional experiments, but to choose sets of variables to intervene on, this algorithm uses an exhaustive search over a separating system. The authors prove guarantees relating this algorithm to their infinite data MI objective.

The paper studies experimental design for causal structure learning. Unlike existing works that focus on experiments with single interventions, it focuses on the experiments that can possibly simulnatenously intervene on multiple variables. The search space becomes much larger because the power set of all variables has exponential cardinality. The paper establish submodularity properties for a two optimization objectives, which allows them to employ/adapt existing variants of greedy algorithm to perform active causal structure learning. The paper also establish theoretical guarantees that shows the algorithm can achieve a constant factor of the optimal objective value with the greedy algorithms.