Budgeted Experiment Design for Causal Structure Learning
AmirEmad Ghassami∗† , Saber Salehkaleybar† , Negar Kiyavash∗† , Elias Bareinboim‡
∗

arXiv:1709.03625v1 [cs.LG] 11 Sep 2017

†

Department of ECE, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA.
Coordinated Science Laboratory, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA.
‡
Department of Computer Science, Purdue University, West Lafayette, Indiana, USA.
†
{ghassam2,sabersk,kiyavash}@illinois.edu, ‡ eb@purdue.edu

Abstract
We study the problem of causal structure learning when the
experimenter is limited to perform at most k non-adaptive experiments of size 1. We formulate the problem of finding the
best intervention target set as an optimization problem, which
aims to maximize the average number of edges whose directions are resolved. We prove that the objective function is submodular and a greedy algorithm is a (1 − 1e )-approximation
algorithm for the problem. We further present an accelerated
variant of the greedy algorithm, which can lead to orders of
magnitude performance speedup. We validate our proposed
approach on synthetic and real graphs. The results show that
compared to the purely observational setting, our algorithm
orients majority of the edges through only a small number of
interventions.

Introduction
The problem of learning the causal relations among a set
of variables is of interest in many scientific domains and
AI. Causal structures are commonly represented by directed
acyclic graphs (DAGs), where the vertices are random variables and a directed edge from variable X to Y indicates that variable X is a direct cause of Y (Pearl 2009;
Spirtes, Glymour, and Scheines 2000). In order to uncover
the causal relations among a set of variables, if the experimenter is restricted to work with only observational data
from the variables, by running conditional independence
based algorithms, he is usually left with some (or in some
cases many) causal relations unresolved. On the other hand,
sufficient experiments involving interventions can identify
the underlying causal graph completely.
An interventional inference algorithm consists of a set
of experiments. In each experiment, the experimenter intervenes on a subset of the variables and collects data from
the intervened system. In this setting, two natural questions
arise:
1. What is the smallest required number of experiments in
order to fully learn the causal graph?
2. For a fixed number of experiments (budget), what portion
of the causal graph can be learned?
The first problem has been addressed in the literature under different assumptions (see Related Work section), while

the second question has received less attention. We address the second question herein. Specifically, we consider
a setup with k experiments, each containing exactly one intervention. The reason we consider single-intervention experiments is that in many applications, such as some experiments in biology, performing simultaneous interventions in
multiple variables may not be feasible.
The majority of the work in the field of causality separates
the design of algorithms from dealing with the statistical
concerns. The reason is that the statistical tests required are
often well-studied in statistics. Many of the seminal works
in the field of causal inference (Verma and Pearl 1991;
Pearl and Bareinboim 2011; Shpitser and Pearl 2006;
Shpitser,
VanderWeele,
and
Robins
2012;
Hyttinen, Eberhardt, and Hoyer 2013; Eberhardt, Glymour, and Scheines 2006; Schölkopf et al. 2012;
Hoyer et al. 2009), and almost entire books (Spirtes,
Glymour, and Scheines 2000) and (Pearl 2009) decouple
the statistical from the algorithmic considerations. We
focus on the algorithm design side of the problem and
formulate it as an optimization problem. We propose an
approximation algorithm which allows an experimenter to
take high advantage of a limited budget of experiments.
Contributions. In our interventional inference algorithm,
first an observational test, such as PC algorithm (Spirtes,
Glymour, and Scheines 2000), is performed on the set of
variables. This test learns the skeleton as well as the orientation of some of the edges of the causal graph. Based on the
result of the initial stage, the complete set of k experiments
is designed in a non-adaptive manner. The advantage of a
non-adaptive approach is that having the complete set of experiments enables the experimenter to perform the interventional experiments in parallel. The formal description of the
problem of interest is provided in the Problem Description
Section. We present the problem of finding the best intervention target set as an optimization problem which aims to
maximize the average number of edges whose directions are
resolved. We prove that this objective function is monotonically increasing and submodular. This implies that a general greedy algorithm is a (1 − 1e )-approximation algorithm.
Since computing the objective function is intractable for a
given set, we present an unbiased estimation of this function. Additionally, we propose a fast estimator with a close

performance to the unbiased estimator (Proposed Approach
Section). In Improved Algorithm Section, we implement an
accelerated variant of the general greedy algorithm through
lazy evaluations, originally proposed by Minoux (Minoux
1978). This algoirthm can lead to orders of magnitude performance speedup (Leskovec et al. 2007). Using synthetic
and real data, in Experiments Section, we show that the proposed approach recovers a significant portion of the edges
by only a few number of interventions in many graphs of
interest.

Related Work
The generic algorithms for purely observational setting are
IC and IC∗ (Pearl 2009), PC and FCI (Spirtes, Glymour, and
Scheines 2000) and Greedy Equivalence Search (GES) algorithm (Meek 1997; Chickering 2002). Such purely observational approaches reconstruct the causal graph up to
the Markov equivalence class, and hence, the direction of
some of the edges may remain unresolved. Albeit under
some stronger assumptions, full structure learning using
merely observational data is feasible (Shimizu et al. 2006;
Hoyer et al. 2009; Peters and Bühlmann 2012). There is a
growing body of research on learning causal structures using
interventional data (Eberhardt 2007; Hauser and Bühlmann
2012; Cooper and Yoo 1999; He and Geng 2008). Specifically, Hauser and Bühlmann (Hauser and Bühlmann 2012)
extended the notion of Markov equivalence to the interventional case and proposed an interventional generalization of GES algorithm. Note that performing an intervention could be viewed as optimally designing a change in the
environment. In (Peters, Bühlmann, and Meinshausen 2016;
Ghassami et al. 2017) the case where data comes from different unknown changes in the environment is studied.
Regarding the first question discussed earlier, (Eberhardt,
Glymour, and Scheines 2005) considered the complete
graph as the underlying causal structure to obtain the worst
case bounds on the number of required experiments under the assumption of causal sufficiency. In that work, both
cases of experiments containing bounded and unbounded
number of interventions were studied. In (Eberhardt 2007;
Hyttinen, Eberhardt, and Hoyer 2013), it was shown that the
problem of designing a proper set of experiments for causal
inference is connected to the problem of finding a separating
system in a graph. Hence, through this connection, the authors derived fundamental bounds on the number of required
experiments for fully learning the graph.
Adaptive algorithms for experiment design were proposed
in (Hauser and Bühlmann 2014; Shanmugam et al. 2015). In
(Hauser and Bühlmann 2014), two algorithms that minimize
the number of experiments in the worst case were proposed.
These algorithms are adaptive, and in the one with polynomial complexity, the size of experiments can be as large
as half the order of the graph, which may not be practical
in many real-life applications. In (Shanmugam et al. 2015),
the authors present information-theoretic lower bounds on
the number of required experiments for both deterministic
and randomized adaptive approaches. They also proposed
an adaptive algorithm that allowed to learn chordal graphs
completely. In (Kocaoglu, Dimakis, and Vishwanath 2017),

the authors considered costs for intervening on each variable
and derived an experiment design algorithm with minimum
total cost that reconstructs the whole structure. They also
considered the case in which the number of experiments is
limited. However, each experiment in their setup is allowed
to include intervening on arbitrary number of variables.
The budgeted experiment design problem is similar to the
influence maximization problem. The goal in the latter is to
find k vertices (seeds) in a given network such that under a
specified influence model, the expected number of vertices
influenced by the seeds is maximized. In (Kempe, Kleinberg, and Tardos 2003) the authors showed that selecting the
most influential vertices is NP-hard and provided the first
provable approximation guarantees for efficient algorithms
based on submodularity of the objective function. There
are many improvements of (Kempe, Kleinberg, and Tardos
2003) such as (Leskovec et al. 2007) and (Chen, Wang, and
Yang 2009). Besides the interpretative differences, an important distinction between the two problems is that in maximum influence problem, the goal is to spread the influence
to the vertices of the graph while in budgeted experiment
design problem, the goal is to pick the initial k vertices in
a way that leads to discovering the orientation of as many
edges as possible. As such the optimal solution to these two
problems for a given graph can be quite different (see the
appendix for an example).

Problem Description
Definition 1. A mixed graph G = (V, E, A) consists of set
of vertices V , set of undirected edges E, and set of directed
edges A. A directed graph G = (V, A) consists of vertex set
V and set of directed edges A. Directed graph G is a DAG
if it is a finite graph with no directed cycles.
A DAG G is called causal if its vertices represent random
variables V = {X1 , ..., Xn } and a directed edge (Xi , Xj )
indicates that variable Xi is a direct cause of variable Xj .
In our model, we assume no restrictions on the underlying
causal DAG. We consider a structural equation model (Pearl
2009) on the variables, which is a collection of n equations Xi = fi (P AXi , Ni ), i ∈ {1, ..., n}, where P AXi
denotes the parents of Xi in G, and Ni s are jointly independent noise variables. Furthermore, we assume faithfulness assumption on the probability distribution PV and
causal sufficiency (Spirtes, Glymour, and Scheines 2000;
Pearl 2009).
Definition 2. Two causal DAGs G1 and G2 over V are
Markov equivalent if every distribution that is compatible
with one of the graphs is also compatible with the other. The
set of all graphs over V is partitioned into a set of mutually
exclusive and exhaustive Markov equivalence classes, which
are the set of equivalence classes induced by the Markov
equivalence relation (Koller and Friedman 2009). The graph
union1 of all DAGs in the Markov equivalence class of a
1

Graph union of a set of mixed graphs over a skeleton, is a
mixed graph with the same skeleton which contains directed edge
(X, Y ), if the direction of the edge between X and Y is either from
the former to later, or not specified in all members of the set. The
rest of the edges in the union remain undirected.

DAG G is called the essential graph of G and is denoted by
Ess(G).
Using purely observational data (referred to as the null
experiment by (Eberhardt, Glymour, and Scheines 2005)),
one can learn the causal structure up to Markov equivalence.
This could be done by performing a “complete” conditional
independence based algorithm which we define as follows.
Definition 3. A complete conditional independence-based
(CCI) algorithm includes first finding the skeleton and vstructures2 of the graph by performing conditional independence tests on the given observational data, and then, applying the Meek rules.
A CCI algorithm leads to learning the Markov equivalence class of the ground truth DAG. After performing CCI
on the ground truth DAG G∗ , we denote the set of directed
edges revealed by the procedure as A(Ess(G∗ )) and the set
of undirected edges by E(Ess(G∗ )).
Interventions can enable us to differentiate among the different causal structures within a Markov equivalence class.
We use the same notion of ideal intervention as in (Eberhardt, Glymour, and Scheines 2005). That is, for an intervention I on variable X ∈ V , the influence of all the variables on the target variable X is removed, and X is randomized by forcing values from an independent distribution
on it. This intervention changes the joint distribution of all
variables in the system for which X is a direct or indirect
cause. We denote this intervention by I = X. An interventional inference algorithm consists of a set of k experiments3
E = {E1 , E2 , ..., Ek }, where each experiment Ei contains mi
(i) (i)
(i)
interventions, i.e., Ei = {I1 , I2 , ..., Imi } for 1 ≤ i ≤ k.
An inference algorithm may be adaptive, in which the experiments are performed sequentially and the information
obtained from the previous experiments is used to design
the next one, or passive, in which all the experiments are
designed beforehand.
The approach we take in our algorithm is to first perform
a CCI algorithm to learn the skeleton and the direction of
the edges in A(Ess(G∗ )), and then design the experiments
in a passive manner4 under a budget constraint. This approach gives the experimenter the ability to perform the experiments in parallel without the need to wait for the result
of one experiment to choose the next one. For example, in
the study of gene regulatory networks (GRNs), for the cases
that the GRN of all cells are the same, experiments can be
performed on different cells simultaneously.
We consider a setup in which a budget of k interventions
is given and the goal is to design a set of experiments to learn
the directions of as many edges in E(Ess(G∗ )) as possible.
2
A v-structure is a structure containing two converging directed
edges whose tails are not connected by an edge.
3
Note that in most of the other work in this area, each intervention is what we refer to as experiment here and hence, each
intervention can contain as many as n variable be randomized.
4
This approach is referred to as the passive setup by (Shanmugam et al. 2015), while (Eberhardt, Glymour, and Scheines
2005) uses the term passive for a setting in which the interventions
are selected without performing the null experiment.

We focus on single-intervention experiments setup in which
for all experiments, we have mi = 1. Therefore, our exper(1)
(2)
(k)
iment set is of the form E = {{I1 }, {I1 }, ..., {I1 }}.
To simplify the notation, we donate the set of intervention
(i)
targets as I = {I1 , ..., Ik }, where Ii = I1 . Moreover,
as shown in (Eberhardt 2007), observing the result of the
null experiment, one can find the orientation of the edge
{Xi , Xj } ∈ E(Ess(G∗ )), if there exists El ∈ E such that
(Xi ∈ El , Xj ∈
/ El ) or (Xj ∈ El , Xi ∈
/ El ). On the
other hand, if for all experiments El ∈ E we have both
Xi ∈ El and Xj ∈ El , the orientation of {Xi , Xj } cannot
be learned. An experiment in which both Xi and Xj are intervened on is called a zero information experiment for Xi
and Xj . Having mi = 1, for all i ∈ {1, ..., k}, prevents such
zero information experiments.
So formally, we focus on the following problem: If the experimenter is allowed to perform k experiments, each of size
1, what portion of the graph could be reconstructed on average? We formalize the problem statement as follows. For
any given set of directed edges A from the underlying DAG
G∗ , define R(A, G∗ ) as the set of edges which do not belong to A(Ess(G∗ )) and their direction can be learned by
applying Meek rules starting from the set of directed edges
A ∪ A(Ess(G∗ )); i.e., the set of edges whose direction is
learned as a result of new information in A. As mentioned
earlier, intervening only on variables I ⊆ V , we learn the
orientation of all the edges intersecting with members of
(I)
I. Denoting the set of these edges by AG∗ , intervening on
variables in I results in learning the direction of edges in
(I)
(I)
R(AG∗ , G∗ ). That is, R(AG∗ , G∗ ) is the set of edges whose
direction is learned as a result of intervening on I and ap(I)
(I)
plying Meek rules on AG∗ . Let D(I, G∗ ) = |R(AG∗ , G∗ )|,
and consider G as the set of all DAGs in the Markov equivalence class of G∗ . Since we do not know the ground truth
we define the average number of edges whose direction can
be learned as a result of intervention on target set I as
D(I) =

1 X
D(I, Gi ).
|G|
Gi ∈G

Therefore, the problem of interest could be stated as finding
intervention target set I ⊆ V of cardinality k which maximizes D(I), that is
max D(I) s.t. |I|= k.

I:I⊆V

(1)

This problem is hard to solve for two reasons: First, finding the optimum I requires a combinatorial search. Second,
even for a given set I, computing D(I) requires considering
all the DAGs in the Markov equivalence class of G∗ , which
in the worst case, say, if G∗ is a complete graph, requires
super-exponential computation time.

Proposed Approach
Definition 4. A set function f : 2V → R is monotonically
increasing if for all sets I1 ⊆ I2 ⊆ V , we have f (I1 ) ≤
f (I2 ).

Definition 5. A set function f : 2V → R is submodular if
for all subsets I1 ⊆ I2 ⊆ V and all v ∈ V \ I2 5 ,

Lemma 1. The set function D is monotonically increasing,
i.e., for sets I1 ⊆ I2 , we have
D(I2 ) ≤ D(I1 ).
Proof. First we show that for a given directed graph Gi ∈ G
the function D(I, Gi ) is a monotonically increasing function of I. In the proposed method, intervening on elements
(I)
of I, we first discover the orientation of the edges in AGi ,
and then applying the Meek rules, we possibly learn the orientation of some extra edges. Having I1 ⊆ I2 implies that
(I )
(I )
AGi1 ⊆ AGi2 . Therefore using I2 , we have more information about the direction of edges. Hence, in the step of
applying Meek rules, by soundness and order-independence
of Meek algorithm, we recover the direction of more extra
(I )
(I )
edges, i.e., R(AGi1 , Gi ) ⊆ R(AGi2 , Gi ), which in turn implies that D(I
P1 , Gi ) ≤ D(I2 , Gi ). Finally, from the relation
1
D(I) = |G|
Gi ∈G D(I, Gi ), the desired result is immediate.
The following lemma plays a fundamental role in the
proof of submodularity of the set function D.
Lemma 2. For sets I1 , I2 ⊆ V ,
(I )

({v})

, Gi )|

, Gi )|. Since I1 ⊆ I2 , as seen
(I )

in the proof of Lemma 1, R(AGi1 , Gi ) ⊆ R(AGi2 , Gi ).

In their seminal work, Nemhauser et. al. showed that if
f is a submodular, monotonically increasing set function
with f (∅) = 0, then the set Î with |Î|= k found by the
greedy algorithm satisfies f (Î) ≥ (1 − 1e ) maxI:|I|=k f (I)
(Nemhauser, Wolsey, and Fisher 1978). That is, the greedy
algorithm is a (1 − 1e )-approximation algorithm. We will use
this result in our proposed approach. We will show that the
set function D is monotonically increasing and submodular,
and hence, the greedy algorithm is a (1 − 1e )-approximation
algorithm to the maximization problem (1).

(I ∪I2 )

(I )

−|R(AGi2 , Gi ) ∩ R(AGi

(I )

f (I1 ∪ {v}) − f (I1 ) ≥ f (I2 ∪ {v}) − f (I2 ).

R(AG∗1

({v})

Similarly, D(I2 ∪ {v}, Gi ) − D(I2 , Gi ) = |R(AGi

(I )

, G∗ ) = R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ ).

(I )

({v})

Therefore, −|R(AGi1 , Gi ) ∩ R(AGi

(I )

, Gi )|≥ −|R(AGi1

({v})

, Gi )∩R(AGi , Gi )|, which implies that D(I1 ∪{v}, Gi )−
D(I1 , Gi ) ≥ D(I2 ∪ {v}, Gi ) − D(I2 , Gi ). This together
with the fact that the function D(I, Gi ) is a monotonically
increasing function of I (observed in the proof of Lemma
1) shows that D(I, Gi ) is aP
submodular function of I. Fi1
nally, we have D(I) = |G|
Gi ∈G D(I, Gi ). Since a nonnegative linear combination of submodular functions is also
submodular, the proof is concluded.
As mentioned in Problem Description Section, another issue regarding solving the optimization problem (1) is the
computational intractability of calculating D(I) for a given
intervention target set I. We propose running Monte-Carlo
simulations of the intervention model for sufficient number of times to obtain an accurate estimation of D(I). The
pseudo-code of a baseline estimator is presented in Subroutine 1. In this subroutine, we first orient the edges in
E(Ess(G∗ )) in one of the possible directions independently according to Bernoulli( 12 ) distribution to obtain a directed graph Ĝ0i . We keep those DAGs which are in the
same Markov equivalence class as G∗ in a multiset G 0 .
To check whether Ĝ0i is in the same Markov equivalence
class as G∗ , it suffices to check if they have the same set
of v-structures (Verma and Pearl 1991). Note that G 0 is a
multiset in which repetition is allowed, and operator ] in
the pseudo-code indicates the multiset addition. Finally, we
calculate the estimated value D̂(I) on G 0 instead of G as
P
D̂(I) = |G10 | G0 ∈G 0 D(I, G0i ).
i

Theorem 2. For any I ⊆ V , D̂(I) obtained from Subroutine 1 is an unbiased estimate of D(I), i.e., E[D̂(I)] =
D(I).
See the appendix for the proof.

See the appendix for the proof.
Theorem 1. The set function D is a supermodular function.
Proof. Due to Lemma 1, it suffices to show that for I1 ⊆
I2 ⊆ V , and v ∈ V , we have D(I1 ∪ {v}) − D(I1 ) ≥
D(I2 ∪ {v}) − D(I2 ). First we show that for a given directed graph Gi ∈ G the function D(I, Gi ) is a submodular
(I ∪{v})
function of I. From Lemma 2, we have R(AGi1
, Gi ) =
(I )

({v})

R(AGi1 , Gi ) ∪ R(AGi

, Gi ). Therefore,

D(I1 ∪ {v}, Gi ) − D(I1 , Gi )
(I ∪{v})

= |R(AGi1

(I )

(I )

, Gi )|−|R(AGi1 , Gi )|
({v})

= |R(AGi1 , Gi ) ∪ R(AGi
({v})

= |R(AGi
5

(I )

, Gi )|−|R(AGi1 , Gi )|

(I )

({v})

, Gi )|−|R(AGi1 , Gi ) ∩ R(AGi

, Gi )|

We can relax it to v ∈ V if f is monotonically increasing.

Our experiments showed that in some cases, out of T
generated graphs in Subroutine 1, only a few of them are
in the same Markov equivalence class as G∗ . Therefore,
we propose an efficient estimator for D(I). The pseudocode of the proposed estimator is presented in Subroutine 2.
In this subroutine for the given mixed graph E(Ess(G∗ )),
we generate N DAGs from the Markov equivalence class
of G∗ as follows: We consider all subsets of size 3 from
V in a random uniformly random order (achieved by uniformly shuffling the labels of elements of V ). For each subset {Xi , Xj , Xk }, we orient the undirected edges among
{Xi , Xj , Xk } independently according to Bernoulli( 21 ) distribution. If the resulting orientation on the induced subgraph on {Xi , Xj , Xk } became a directed cycle or a new
v-structure which was not in Ess(G∗ ), we redo the orienting. We keep checking all the subsets of size 3 until the induced subgraph on all of them are directed and none of them

Subroutine 1 Baseline D(I) Estimator Subroutine
Input: Ess(G∗ ), intervention target set I, and T .
initiate: G 0 = ∅, i = 1
for j = 1 to T do
Orient each undirected edge {u, v} ∈ E(Ess(G∗ )) as
(u, v) or (v, u) independently according to Bern( 21 ) to
get directed graph Ĝ0j .
if Ĝ0j is a DAG and is in the Markov equivalence class
of G∗ then
G0i = Ĝ0j
G 0 = G 0 ] G0i
i=i+1
end if
end for
P
Output: D̂(I) = |G10 | G0 ∈G 0 |D(I, G0i )|
i

is a new v-structure, which did not exist in Ess(G∗ ), or a directed cycle.
Lemma 3. Each generated directed graph G0i in D(I) Estimator Subroutine belongs to the Markov equivalence class
of G∗ .
Proof. We require the following claim for the proof:
Claim 1. A chordal 6 graph has a directed cycle only if it
has a directed cycle of size 3.
To prove Claim 1, if the directed cycle is of size 3 itself,
the claim is trivial. Suppose the cycle Cn is of size n > 3.
Relabel the vertices of Cn to have Cn = (v1 , ..., vn , v1 ).
Since the graph is chordal, Cn has a chord and hence we
have a triangle on vertices {vi , vi+1 , vi+2 } for some i. If
the direction of {vi , vi+2 } is (vi+2 , vi ), we have the directed
cycle (vi , vi+1 , vi+2 , vi ) which is of size 3. Otherwise, we
have the directed cycle Cn−1 = (v1 , ..., vi , vi+2 , .., vn , v1 )
on n − 1 vertices. Relabeling the vertices from 1 to n − 1
and repeating the above reasoning concludes the claim.
All the components in the undirected subgraph of
Ess(G∗ ) are chordal (Hauser and Bühlmann 2012). Therefore, by Claim 1, to insure that a generated directed graph
is a DAG, it suffices to make sure that it does not have any
directed cycles of length 3, which is one of the checks that
we do in the proposed procedure. For checking if the generated DAG is in the same Markov equivalence class as G∗ ,
it suffices to check if they have the same set of v-structures
(Verma and Pearl 1991), which is the other check that we do
in the proposed procedure.
We add the resulted generated DAG to a multiset G 0 . Finally, we calculate the estimated value D̂(I) on G 0 instead
P
of G as D̂(I) = |G10 | G0 ∈G 0 D(I, G0i ).
i
Next we will consider the required cardinality of the set
G 0 to obtain a desired accuracy in estimation D(I). We use
Chernoff bound for this purpose.

Subroutine 2 D(I) Estimator Subroutine
Input: Ess(G∗ ), intervention target set I, and N .
initiate: G 0 = ∅
for i = 1 to N , generate G0i as follows: do
Uniformly shuffle the order of the elements of V .
while the induced subgraph on any subset of size 3 of
the variables is not directed, or a directed cycle, or a
v-structure which was not in Ess(G∗ ) do
for all {Xi , Xj , Xk } ⊆ V do
Orient the undirected edges among {Xi , Xj , Xk }
independently according to Bern( 21 ) until it becomes a directed structure which is not a directed
cycle or a v-structure which was not in Ess(G∗ ).
end for
end while
G 0 = G 0 ] G0i
end for
P
Output: D̂(I) = |G10 | G0 ∈G 0 D(I, G0i )
i

Proposition 1 (Chernoff Bound). Let X1 , ..., XN be independent random variables such that for all i, 0 ≤ Xi ≤ 1.
PN
Let µ = E[ i=1 Xi ]. Then
P (|

N
X

Xi − µ|≥ µ) ≤ 2 exp(−

i=1

Theorem 3. For an estimator with E[D(I, G0i )] =
D(I), given set I and , δ > 0, if N = |G 0 |>
|E(Ess(G∗ ))|(2+)
ln( 2δ ), then
2
D(I)(1 − ) < D̂(I) < D(I)(1 + ),
with probability larger than 1 − δ.
D(I,G0 )

Proof. Define Xi = |E(Ess(Gi∗ ))| , for i ∈ {1, ..., N }. By
1
the assumption of the theorem, E[Xi ] = |E(Ess(G
∗ ))| D(I)
(Note that we saw in the proof of Theorem 2 that this assumption is satisfied). Using Chernoff bound we have
N
X
P (|
Xi −
i=1

N
N
D(I)|≥ 
D(I))
|E(Ess(G∗ ))|
|E(Ess(G∗ ))|

N 2
D(I))
|E(Ess(G∗ ))|(2 + )
N 2
≤ 2 exp(−
).
|E(Ess(G∗ ))|(2 + )
PN
Therefore, P (| N1 i=1 D(I, G0i ) − D(S)|≥ D(I)) ≤
N 2
2 exp(− |E(Ess(G
∗ ))|(2+) ). Hence, P (|D̂(I) − D(I)|<
≤ 2 exp(−

2

N
D(I)) > 1 − 2 exp(− |E(Ess(G
∗ ))|(2+) ).
∗

6

A chord of a cycle is an edge not in the cycle whose endpoints
are in the cycle. A hole in a graph is a cycle of length at least 4
having no chord. A graph is chordal if it has no hole.

2
µ).
2+

Setting N > |E(Ess(G2 ))|(2+) ln( 2δ ), upper bounds the
right hand side with 1 − δ and concludes the desired result.

Algorithm 1 General Greedy Algorithm
Input: Joint distribution over V , and budget k.
Obtain Ess(G∗ ) by performing a CCI algorithm.
initiate: I0 = ∅
for i = 1 to k do
vi = arg maxv∈V \Ii−1 D̂(Ii−1 ∪ {v}) − D̂(Ii−1 )
Ii = Ii−1 ∪ {vi }
end for
Output: Î = Ik

Our general greedy algorithm is presented in Algorithm
1. Defining the marginal gain of variable v when the previous chosen set is I as ∆v (I) = D(I ∪ {v}) − D(I), the
algorithm iteratively adds a variable which has the largest
marginal gain to the intervention target set until it runs out
of budget. Either D(I) estimators can be used in calculating
an estimation of the marginal gains.
Theorem 4. For any 0 , δ 0 > 0, there exists , δ > 0, such
that if for any set I, D(I)(1 − ) < D̂(I) < D(I)(1 + )
with probability larger than 1 − δ, then Algorithm 1 is a
(1− 1e −0 )-approximation algorithm with probability larger
than 1 − δ 0 .
See the appendix for the proof.

Improved Algorithm
Submodularity of function D can be exploited to implement an accelerated variant of the general greedy algorithm
through lazy evaluations, originally proposed by Minoux
(Minoux 1978). In each round of the general greedy algorithm, we check the marginal gain ∆v (I) for all remaining
vertices in V \I. Note that as a consequence of submodularity of function D, the set function ∆v is monotonically
decreasing. The main idea of the improved greedy algorithm is to take advantage of this property to avoid checking all the variables in each round of the algorithm. The
idea is as follows: Suppose for vertices v1 and v2 , in the
i-th round of the algorithm we have obtained marginal gains
∆v1 (Ii ) > ∆v2 (Ii ). Now, if in the (i + 1)-th round we
have calculate ∆v1 (Ii+1 ) and have ∆v1 (Ii+1 ) > ∆v2 (Ii ),
then by monotonic decreasing property, we can conclude
that ∆v1 (Ii+1 ) > ∆v2 (Ii+1 ), and hence, there is no need
to calculate ∆v2 (Ii+1 ).
The detailed description of the improved greedy algorithm is given in Algorithm 2. The idea can be formalized
as follows (Leskovec et al. 2007): We define a profit parameter pv for each variable v and initialize the value of all of
them with ∞. Also, for all variables we define an update
flag updv , which will be set to false at the beginning of every round of the algorithm and will be switched to true if
we update pv with the value of the marginal gain of vertex
v. In each round of the algorithm, we pick vertex v ∈ V \I
with the largest profit, update its profit with the value of the
marginal gain of v, and set updv to true. This process is repeated until the vertex with the largest profit is already updated, i.e., its update flag is true. Then we add this vertex to
I and end the round. For example, if in a round, the vertex
v has the highest profit and after updating the profit of this

Algorithm 2 Improved Greedy Algorithm
Input: Joint distribution over V , and budget k.
Obtain Ess(G∗ ) by performing a CCI algorithm.
initiate: I0 = ∅, and pv = ∞, ∀v ∈ V .
for i = 1 to k do
updv =false, ∀v ∈ V \Ii−1
while true do
v ∗ = arg maxv∈V \Ii−1 pv
if updv∗ then
Ii = Ii−1 ∪ {v ∗ }
break;
else
pv∗ = D̂(Ii−1 ∪ {v ∗ }) − D̂(Ii−1 )
updv∗ = true
end if
end while
end for
Output: Î = Ik

vertex, pv is still larger than all the other profits, we do not
need to evaluate the marginal gain of any other vertex and
we add v to I.
The correctness of the improved algorithm follows directly from submodularity and Theorem 4 holds for Algorithm 2 as well. This algorithm can lead to orders of magnitude performance speedup (Leskovec et al. 2007).

Experiment
Synthetic Graphs
In this subsection, we evaluate the performance of the proposed improved greedy algorithm on synthetically generated chordal graphs. We use randomly chosen partial elimination ordering (PEO)7 of the vertices to generate our
underlying chordal graphs (Hauser and Bühlmann 2014;
Shanmugam et al. 2015). For each graph, we pick a random
ordering of the vertices. Starting from the vertex v with the
highest order, we connect all the vertices with lower order
to v with probability inversely proportional to the order of
v. Then, we connect all the parents of v with directed edges,
where each directed edge is oriented from the parent with
the lower order to the parent with the higher order. In order to make sure that the generated graph will be connected,
if vertex v is not connected to any of the vertices with the
lower order, we pick one of them uniformly at random and
set it as the parent of v.
To evaluate the performance of the proposed algorithm,
for any underlying graph, we consider the ratio of the number of edges whose directions are discovered merely as a result of interventions to the number of edges whose directions
were not resolved from the observational data. Note that due
to our specific graph generating approach, the orientation of
none of the edges is learned from the observational data. We
7
A perfect elimination ordering {v1 , v2 , ..., vn } on the vertices
of an undirected chordal graph is such that for all i, the induced
neighborhood of vi on the subgraph formed by {v1 , v2 , ..., vi−1 }
is a clique.

1

1

Discovered edge ratio

Discovered edge ratio

0.95
0.9
0.85
0.8
0.75

0.8
0.6
0.4
0.2
0
Ecoli1 Ecoli2 Yeast1 Yeast2 Yeast3

0.7
1

2

3

4

5

Figure 3: Discovered edge ratio in five GRNs from DREAM
3 challenge.

Budget (k)

Figure 1: Discovered edge ratio versus budget for n = 20.

Discovered edge ratio

1

0.95

0.9

0.85
10

15

20

25

30

Graph order (n)

Figure 2: Discovered edge ratio versus the order of graph for
k = 3.
generated 100 instances of chordal DAGs of order 20. In all
experiments, the improved greedy algorithm is utilized and
Subroutine 2 is used for estimating D(I). Figure 1 depicts
the discovered edge ratio with respect to the budget k. As
seen in this figure, three interventions suffices to discover
the direction of more than 90% of the edges whose direction
was unknown prior to performing interventions. Further to
investigate the effect of the order of the graph on the performance of the proposed algorithm, we evaluated the discovered edge ratio for budget k = 3 on graphs with order
n ∈ {10, 15, 20, 25, 30} in Figure 2.
Furthermore, to compare the performance of the proposed
algorithm with the optimal solution, we generated 100 instances of chordal DAGs of order 10 and performed a brute
force search to find the optimal solution of (1) for budget
k = 2. The discovered edge ratio was 0.9 and 0.916 for our
proposed algorithm and the optimal solution, respectively.

Real Graphs
We evaluated the performance of the proposed algorithm in
gene regulatory networks (GRN). GRN is a collection of biological regulators that interact with each other. In GRN, the
transcription factors are the main players to activate genes.
The interactions between transcription factors and regulated
genes in a species genome can be presented by a directed

graph. In this graph, links are drawn whenever a transcription factor regulates a gene’s expression. Moreover, some of
vertices have both functions, i.e., are both transcription factor and regulated gene.
We considered GRNs in “DREAM 3 In Silico Network”
challenge, conducted in 2008 (Marbach et al. 2009). The
networks in this challenge were extracted from known biological interaction networks. Since we know the true causal
structures in these GRNs, we can obtain Ess(G∗ ) and give
it as an input to the proposed algorithm. Figure 3 depicts the
discovered edge ratio in five networks extracted from GRNs
of E-coli and Yeast bacteria with budget k = 5. The order of
each network is 100. As it can be seen, the discovered edge
ratio is at least 0.65 in all GRNs.

Conclusion
We studied the problem of experiment design for causal
inference when only a limited number of experiments are
available. In our model, each experiment consists of intervening on a single vertex, which makes results in taking
the most advantage from a finite budget of interventions and
makes the model suitable for applications in which intervening on several variables simultaneously is not feasible. Also,
in our model, experiments are designed merely based on the
result of an initial purely observational test, which enables
the experimenter to perform the interventional tests in parallel. We addressed the following question: “How much of the
causal structure can be learned when only a limited number
of experiments are available?” We formulated the problem
of finding the best intervention target set as an optimization problem which aims to maximize the average number
of edges whose directions are discovered. We showed that
the objective function is monotonically increasing and submodular. Consequently, the greedy algorithm is a (1 − 1e )approximation algorithm for this problem. Moreover, we
proposed estimation methods in order to compute the objective function for a given set of intervention targets. We
further presented an accelerated variant of the greedy algorithm which can lead to orders of magnitude performance
speedup. We examined our proposed improved greedy algorithm on synthetic as well as real graphs. The results show
that a significant portion of causal graphs can be learned by
only a few number of interventions.

Appendices
Example of Comparison with the Influence
Maximization Problem
Suppose k = 1. Figure 4 depicts a graph for which the optimal solution to the influence maximization problem is different from the optimal solution to the budgeted experiment
design problems. Clearly, influencing vertex v1 leads to influencing all the vertices in the graph, and hence, this vertex
is the solution to the influence maximization problem. But,
intervening on v1 leads to discovering the orientation of only
3 edges, while intervening on, say v2 , leads to discovering
the orientation of 5 edges.

Proof of Lemma 2
(I )

(I )

The direction R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ )
⊆
(I ∪I )
R(AG∗1 2 , G∗ ) is proved in the proof of Lemma 1.
Also, as observed in the proof of Lemma 1, we have
(I ∪I )
(I )
(I )
R(AG∗1 2 , G∗ ) ⊆ R(R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ ), G∗ ).
(I ∪I )
Therefore, in order to prove that R(AG∗1 2 , G∗ ) ⊆
(I1 )
(I
)
R(AG∗ , G∗ ) ∪ R(AG∗2 , G∗ ), it suffices to show
(I )
(I )
that
R(R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ ), G∗ )
⊆
(I1 )
(I2 )
∗
∗
R(AG∗ , G ) ∪ R(AG∗ , G ), for which it suffices to
(I )
(I )
show that if e 6∈ R(AG∗1 , G∗ ) and e 6∈ R(AG∗2 , G∗ ), then
(I )
(I )
e 6∈ R(R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ ), G∗ ).
(I )

Proof by contradiction. Let e 6∈ R(AG∗1 , G∗ ) and
(I )
e 6∈ R(AG∗2 , G∗ ), but its orientation is learned in the
(I )
first iteration of applying Meek rules to R(AG∗1 , G∗ ) ∪
(I2 )
R(AG∗ , G∗ ) ∪ A(Ess(G∗ )). Then, we have learned the orientation of e due to one of Meek rules (Verma and Pearl
1992):

• Rule 3. e = {a, b} is oriented as (a, b) if ∃c, d s.t. e1 =
(I )
(I )
(c, b) ∈ R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ ) ∪ A(Ess(G∗ )),
(I )
(I )
e2 = (d, b) ∈ R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ ) ∪
∗
∗
A(Ess(G )), {a, c} ∈ skeleton of G , {a, d} ∈ skeleton
of G∗ , and {c, d} 6∈ skeleton of G∗ .
• Rule 4. e = {a, b} is oriented as (a, b) and e = {b, c} is
(I )
oriented as (c, b) if ∃d s.t. e1 = (d, c) ∈ R(AG∗1 , G∗ ) ∪
(I2 )
R(AG∗ , G∗ ) ∪ A(Ess(G∗ )), {a, c} ∈ skeleton of G∗ ,
{a, d} ∈ skeleton of G∗ , and {b, d} 6∈ skeleton of G∗ .
In what follows, we show that the orientation of e cannot
be learned due to any of the Meek rules unless it belongs to

𝑣2

𝑣4

𝑣3

Figure 4: Example of comparison with the influence maximization problem.
(I )

(I )

R(AG∗1 , G∗ ) or R(AG∗2 , G∗ ).
Rule 1.
(I )
Without loss of generality, assume e1 ∈ R(AG∗1 , G∗ ) ∪
∗
A(Ess(G )). Therefore, we should have the condition of
rule 1 satisfied when only intervening on I1 as well, which
(I )
implies that e ∈ R(AG∗1 , G∗ ), which is a contradiction.
Rule 2.
(I )
If both e1 and e2 belong to R(AG∗1 , G∗ ) ∪ A(Ess(G∗ ))
(I )
(or R(AG∗2 , G∗ ) ∪ A(Ess(G∗ ))), then we should have the
condition of rule 2 satisfied when only intervening on I1 (or
(I )
I2 ) as well, which implies that e ∈ R(AG∗1 , G∗ ) (or e ∈
(I2 )
R(AG∗ , G∗ )), which is a contradiction. Therefore, it suffices to show that the case that e1 belongs to exactly one of
(I )
(I )
R(AG∗1 , G∗ ) ∪ A(Ess(G∗ )) or R(AG∗2 , G∗ ) ∪ A(Ess(G∗ ))
and e2 belongs only to the other one, does not happen. To
this end, it suffices to show that there does not exist interven(I)
tion target set I such that e1 ∈ R(AG∗ , G∗ ) ∪ A(Ess(G∗ )),
(I)
and e, e2 6∈ R(AG∗ , G∗ ) ∪ A(Ess(G∗ )), i.e., there does not
exist intervention target set I that has structure S0 , depicted
in Figure 5, as a subgraph of Ess(G∗ ) after applying the ori(I)
entations learned from R(AG∗ , G∗ ).
𝑐
𝑒1

• Rule 1. e = {a, b} is oriented as (a, b) if ∃c s.t. e1 =
(I )
(I )
(c, a) ∈ R(AG∗1 , G∗ )∪R(AG∗2 , G∗ )∪A(Ess(G∗ )), and
∗
{c, b} 6∈ skeleton of G .
• Rule 2. e = {a, b} is oriented as (a, b) if ∃c s.t. e1 =
(I )
(I )
(a, c) ∈ R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ ) ∪ A(Ess(G∗ )),
(I )
(I )
and e2 = (c, b) ∈ R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ ) ∪
∗
A(Ess(G )).

𝑣1

𝑎

𝑆0

𝑏

Figure 5: Structure S0
(I)

If e1 ∈ AG∗ , then a ∈ I or c ∈ I, which implies
(I)
(I)
e ∈ AG∗ or e2 ∈ AG∗ , respectively, and hence, e ∈
(I)
(I)
R(AG∗ , G∗ ) or e2 ∈ R(AG∗ , G∗ ), respectively. Therefore,
(I)
in either case, e ∈ R(AG∗ , G∗ ), and S0 will not be a sub(I)
graph. Therefore, e1 6∈ AG∗ , and hence, e1 was learned by
applying one of the Meek rules. We consider each or the
rules in the following:
• If we have learned the orientation of e1 from rule 1, then
we should have had one of the structures in Figure 6
as a subgraph of Ess(G∗ ) after applying the orientations
(I)
learned from R(AG∗ , G∗ ). In case of structure S1 , using
rule 1 on subgraph induced on vertices {v1 , a, b}, we will
also learn (a, b). In case of structure S2 , using rule 4, we
will also learn (b, c). Therefore, we cannot learn only the
direction of e1 and hence, S0 will not be a subgraph.

𝑐

𝑐

𝑒1

𝑒1
𝑏

𝑎

𝑏

𝑎

𝑣1

𝑣1

𝑆1

𝑆2

Figure 6: Rule 1
• If we have learned the orientation of e1 from rule 3, then
we have had one of the structures in Figure 7 as a subgraph of Ess(G∗ ) after applying the orientations learned
(I)
from R(AG∗ , G∗ ). In case of structures S3 and S4 , using
rule 1 on subgraph induced on vertices {v2 , c, b}, we will
also learn (c, b). In case of structure S5 , using rule 3 on
subgraph induced on vertices {b, v2 , c, v1 }, we will also
learn (b, c). Therefore, we cannot learn only the direction
of e1 and hence, S0 will not be a subgraph.
𝑐

𝑐

𝑣2

𝑎

𝑏

𝑣2

𝑣1

𝑎

𝑒1

𝑆4

𝑣1

𝑎

𝑆5

𝑣1

Figure 7: Rule 3
• If we have learned the orientation of e1 from rule 4, then
we have had one of the structures in Figure 8 as a subgraph of Ess(G∗ ) after applying the orientations learned
(I)
from R(AG∗ , G∗ ). In case of structures S6 , using rule 1
on subgraph induced on vertices {v1 , c, b}, we will also
learn (c, b). In case of structure S7 , using rule 1 on subgraph induced on vertices {v2 , v1 , b}, we will also learn
(v1 , b), and then using rule 4 on subgraph induced on vertices {b, a, v2 , v1 }, we will also learn (a, b). In case of
structures S8 , using rule 4 on subgraph induced on vertices {b, v2 , v1 , c}, we will also learn (b, c). Therefore, we
cannot learn only the direction of e1 and hence, S0 will
not be a subgraph.
𝑣1

𝑣1

𝑣1

𝑣2

𝑐 𝑣2

𝑐

𝑒1

𝑎

𝑆6

𝑏

𝑣2

𝑐
𝑒1

𝑒1

𝑎

𝑆7

𝑏

𝑎

𝑆8

𝑏

Figure 8: Rule 4
• If we have learned the orientation of e1 from rule 2, then
we should have had one of the structures in Figure 9
as a subgraph of Ess(G∗ ) after applying the orientations
(I)
learned from R(AG∗ , G∗ ). In case of structure S9 , using
rule 1 on subgraph induced on vertices {v1 , c, b}, we will

𝑣1

𝑐
𝑒1
𝑎

𝑏

𝑒1

𝑆3

𝑣1

𝑐

𝑏 𝑣2
𝑒1

also learn (c, b) and hence, S0 will not be a subgraph. In
case of structure S10 , if v1 ∈ I, then the direction of the
edge {v1 , b} will be also known. If the direction of this
edge is (v1 , b), then using rule 2 on subgraph induced on
vertices {a, v1 , b}, we will also learn (a, b); otherwise, using rule 2 on subgraph induced on vertices {b, v1 , c}, we
will also learn (c, b). Therefore, v1 6∈ I. Also, as mentioned earlier, a 6∈ I. Therefore, we have learned the orientation of (a, v1 ) from applying Meek rules.
In the triangle induced on vertices {v1 , b, a}, we have
learned only the orientation of one edge, which is (a, v1 ).
But as seen in structures S1 to S9 , all of them lead to
learning the orientation of at least 2 edges of a triangle. In
the following, we will show that a structure of form S10 ,
does not lead to learning the orientation of only (a, v1 )
and making S10 a subgraph either.
𝑐
𝑒1
𝑎

𝑏

𝑆9

𝑏

𝑆10

Figure 9: Rule 2
Suppose we had learned (a, v1 ) via a structure of form
S10 , as depicted in Figure 10(a). Using rule 4 on subgraph
induced on vertices {v2 , v1 , c, b}, we will also learn (b, c).
Therefore, we should have the edge {v2 , c} too. Also, using rule 2 on triangle induced on vertices {v2 , v1 , c}, the
orientation of this edges should be (v2 , c). Therefore, in
order to have S10 as a subgraph, we need to have the
structure depicted in Figure 10(b) as a subgraph. As seen
in Figure 10(b), we again have a structure similar to S10 :
a complete skeleton K5 , which contains (vj , c), (a, vj ),
{vj , b}, for j ∈ {1, 2} and (v2 , v1 ), with a triangle on
vertices {v2 , b, a}, in which we have learned only the orientation of (a, v2 ).
𝑣1

𝑣1

𝑣2

𝑐

𝑣2

𝑐

𝑒1

𝑎

(𝑎)

𝑏

𝑒1

𝑎

(𝑏)

𝑏

Figure 10: Step of the induction.
We claim that this procedure always repeats, i.e., at step
i, we end up with skeleton Ki , which contains (vj , c),
(a, vj ), {vj , b}, for j ∈ {1, ..., i} and (vk , vj ), for 1 ≤
j < k ≤ i, with a triangle induced on vertices {vi , b, a},
in which we have learned only the orientation of (a, vi ).
We prove this claim by induction. We have already proved
the base of the induction above. For the step of the induction, suppose the hypothesis is true for i − 1. Add vertex
vi to form a structure of form S10 for (a, vi−1 ). vi should
be adjacent to vj , for j ∈ {1, ..., i − 2}; otherwise, using
rule 4 on subgraph induced on vertices {vi , vi−1 , vj , b},
we will also learn (b, vj ). Moreover, using rule 2 on triangle induced on vertices {vi , vi−1 , vj }, the direction of

{vi , vj } should be (vi , vj ). Also, using rule 4 on subgraph induced on vertices {vi , vi−1 , c, b}, we will also
learn (b, c). Therefore, we should have the edge {vi , c}
too.
We showed that S0 is a subgraph only if S10 is a
subgraph, and S10 is a subgraph only if the structure in
Figure 10(b) is a subgraph, and this chain of required
subgraphs continue. Therefore, since the order of the
graph is finite, there exist a step where since we cannot
add a new vertex, it is not possible to have one of the
required subgraphs, and hence we conclude that S0 is not
a subgraph.

∗

Since we can orient E(Ess(G∗ )) in 2|E(Ess(G ))| ways and
noting that we used the notation G0 for those DAGs which
are in the Markov equivalence of G∗ , we have
1
X
|E(Ess(G∗ ))|
D(I, G0i ) 2 |G|
E[D(I, G0 )] =
G0i ∈G

=

Rule 4.
(I )
Without loss of generality, assume e1 ∈ R(AG∗1 , G∗ ) ∪
∗
A(Ess(G )). Therefore, we should have the condition of
rule 4 satisfied when only intervening on I1 as well, which
(I )
implies that e ∈ R(AG∗1 , G∗ ), which is a contradiction.
The argument above proves that there is no edge e such
(I )
(I )
that e 6∈ R(AG∗1 , G∗ ) and e 6∈ R(AG∗2 , G∗ ), but e ∈
(I )
(I )
R(R(AG∗1 , G∗ ) ∪ R(AG∗2 , G∗ ), G∗ ).

Therefore, with N := |G 0 |, we have
N
N
1 X
1 X
D(I, G0i )] = EN [E[
D(I, G0i )|N ]]
N i=1
N i=1

= EN [

N
N
1 X
1 X
E[D(I, G0i )]] = EN [
E[D(I, G0i )]]
N i=1
N i=1

= EN [

N
1 X
1
D(I)] = EN [ N D(I)] = EN [D(I)] = D(I).
N i=1
N

Proof of Theorem 4
∗

Let I
have

=

{v1∗ , ..., vk∗ }

(a)

+

k
X
∗
[D(Ii ∪ {v1∗ , ..., vj∗ }) − D(Ii ∪ {v1∗ , ..., vj−1
})]
j=1

(b)

Proof of Theorem 2

E[D(I, Ĝ0 )|Ĝ0 ∈ G]
X
=
D(I, Ĝ0 = Ĝ0i )P (outputing Ĝ0 = Ĝ0i |Ĝ0 ∈ G)
Ĝ0i ∈Gtotal

D(I, Ĝ0 = Ĝ0i )P (outputing Ĝ0 = Ĝ0i |Ĝ0 ∈ G)

Ĝ0i ∈G

+

X

D(I, Ĝ0 = Ĝ0i )P (outputing Ĝ0 = Ĝ0i |Ĝ0 ∈ G)

Ĝ0i 6∈G

=

X

0

D(I, Ĝ =

Ĝ0i )P (outputing

0

Ĝ =

Ĝ0i |Ĝ0

∈ G)

Ĝ0i ∈G

+

X

=

D(I, Ĝ0 = Ĝ0i )

Ĝ0i ∈G

=

X
Ĝ0i ∈G

[D(Ii ∪ {vj∗ }) − D(Ii )],

(2)
where (a) follows from Lemma 1, and (b) follows from Theorem 1. Define D̂i,v,1 and D̂i,v,2 as the first and second calls
of subroutine in i-th step for variable vv, respectively. By
the assumption of the theorem we have
D(Ii ∪ {vj∗ }) − D(Ii ∪ {vj∗ }) < D̂i,vj∗ ,1 (Ii ∪ {vj∗ }),
with probability larger than 1 − δ. Therefore
D(Ii ∪ {vj∗ }) < D̂i,vj∗ ,1 (Ii ∪ {vj∗ }) + D(I ∗ ),
with probability larger than 1 − δ. Similarly
D̂i,vj∗ ,2 (Ii ) < D(Ii ) + D(Ii )

w.p. > 1 − δ,

⇒ −D(Ii ) < −D̂i,vj∗ ,2 (Ii ) + D(I ∗ )

w.p. > 1 − δ,

Therefore,
D(I, Ĝ0 = Ĝ0i ) × 0

D(Ii ∪ {vj∗ }) − D(Ii ) < D̂i,vj∗ ,1 (Ii ∪ {vj∗ })

Ĝ0i 6∈G

X

k
X
j=1

Let Gtotal be all directed graphs with the same skeleton as
G∗ and Ĝ0 ∈ Gtotal .

X

∈ arg maxI:I⊆V,|I|=k D(I). We

D(I ∗ ) ≤ D(I ∗ ∪ Ii ) = D(Ii )

≤ D(Ii ) +

=

1 X
D(I, G0i ) = D(I).
|G| 0
Gi ∈G

E[D̂(I)] = E[
Rule 3.
Since edges e1 and e2 form a v-structure, they should
appear in A(Ess(G∗ )) as well. Therefore, we should have
the condition of rule 3 satisfied when only intervening on
(I )
I1 as well, which implies that e ∈ R(AG∗1 , G∗ ), which is a
contradiction.

2|E(Ess(G∗ ))|

D(I, Ĝ0 = Ĝ0i )

0

P (outputing Ĝ =
P (Ĝ0

Ĝ0i , Ĝ0

∈ G)

P (outputing Ĝ0 = Ĝ0i )
P (Ĝ0 ∈ G)

∈ G)

− D̂i,vj∗ ,2 (Ii ) + 2D(I ∗ )

w.p. > 1 − 2δ.

(3)

Also, by the definition of the greedy algorithm,
D̂i,vj∗ ,1 (Ii ∪ {vj∗ }) − D̂i,vj∗ ,2 (Ii )

.

≤ D̂i,vi+1 ,1 (Ii ∪ {vi+1 }) − D̂i,vi+1 ,2 (Ii )
= D̂i,vi+1 ,1 (Ii+1 ) − D̂i,vi+1 ,2 (Ii ),

(4)

and similar to (3), we have
D̂i,vi+1 ,1 (Ii+1 ) − D̂i,vi+1 ,2 (Ii ) < D(Ii+1 )
− D(Ii ) + 2D(I ∗ )

w.p. > 1 − 2δ.

(5)

Therefore, from equations (3), (4), and (5) we have
D(Ii ∪{vj∗ })−D(Ii ) < D(Ii+1 )−D(Ii )+4D(I ∗ ), (6)
with probability larger than 1 − 4δ. Plugging (6) back in (2),
we get
D(I ∗ ) < D(Ii ) +

k
X
[D(Ii+1 ) − D(Ii ) + 4D(I ∗ )]
j=1

= D(Ii ) + k[D(Ii+1 ) − D(Ii )] + 4kD(I ∗ ),
with probability larger than 1 − 4kδ. Therefore,
D(I ∗ ) − D(Ii )
< k[D(I ∗ ) − D(Ii )] − k[D(I ∗ ) − D(Ii+1 )] + 4kD(I ∗ ),
with probability larger than 1−4kδ. Defining ai := D(I ∗ )−
D(Ii ), and noting that a0 = D(I ∗ ), by induction we have
ak = D(I ∗ ) − D(Ik )
< (1 −

k−1
X
1 k
1
) D(I ∗ ) + 4D(I ∗ )
(1 − )j
k
k
j=0

1
< [ + 4k]D(I ∗ )
e
It concludes that
1
D(Ik ) > (1 − − 4k)D(I ∗ )
e
0

w.p. > 1 − 4k 2 δ.

w.p. > 1 − 4k 2 δ.
0


δ
Therefore, for  = 4k
and δ = 4k
2 , Algorithms 1 is a (1 −
1
0
−
)-approximation
algorithm
with
probability larger than
e
1 − δ0 .

References
[Chen, Wang, and Yang 2009] Chen, W.; Wang, Y.; and
Yang, S. 2009. Efficient influence maximization in social
networks. In Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining, 199–208. ACM.
[Chickering 2002] Chickering, D. M. 2002. Optimal structure identification with greedy search. Journal of machine
learning research 3(Nov):507–554.
[Cooper and Yoo 1999] Cooper, G. F., and Yoo, C. 1999.
Causal discovery from a mixture of experimental and observational data. In Proceedings of the Fifteenth conference
on Uncertainty in artificial intelligence, 116–125. Morgan
Kaufmann Publishers Inc.
[Eberhardt, Glymour, and Scheines 2005] Eberhardt,
F.;
Glymour, C.; and Scheines, R. 2005. On the number of
experiments sufficient and in the worst case necessary to
identify all causal relations among n variables. 178–184.
[Eberhardt, Glymour, and Scheines 2006] Eberhardt,
F.;
Glymour, C.; and Scheines, R. 2006. N-1 experiments
suffice to determine the causal relations among n variables.
In Innovations in machine learning. Springer. 97–112.

[Eberhardt 2007] Eberhardt, F. 2007. Causation and intervention. Unpublished doctoral dissertation, Carnegie Mellon University.
[Ghassami et al. 2017] Ghassami, A.; Salehkaleybar, S.;
Kiyavash, N.; and Zhang, K. 2017. Learning causal
structures using regression invariance.
arXiv preprint
arXiv:1705.09644.
[Hauser and Bühlmann 2012] Hauser, A., and Bühlmann, P.
2012. Characterization and greedy learning of interventional
markov equivalence classes of directed acyclic graphs. Journal of Machine Learning Research 13(Aug):2409–2464.
[Hauser and Bühlmann 2014] Hauser, A., and Bühlmann, P.
2014. Two optimal strategies for active learning of causal
models from interventional data. International Journal of
Approximate Reasoning 55(4):926–939.
[He and Geng 2008] He, Y.-B., and Geng, Z. 2008. Active
learning of causal networks with intervention experiments
and optimal designs. Journal of Machine Learning Research
9(Nov):2523–2547.
[Hoyer et al. 2009] Hoyer, P. O.; Janzing, D.; Mooij, J. M.;
Peters, J.; and Schölkopf, B. 2009. Nonlinear causal discovery with additive noise models. In Advances in neural
information processing systems, 689–696.
[Hyttinen, Eberhardt, and Hoyer 2013] Hyttinen, A.; Eberhardt, F.; and Hoyer, P. O. 2013. Experiment selection for
causal discovery. Journal of Machine Learning Research
14(1):3041–3071.
[Kempe, Kleinberg, and Tardos 2003] Kempe, D.; Kleinberg, J.; and Tardos, É. 2003. Maximizing the spread
of influence through a social network. In Proceedings
of the ninth ACM SIGKDD international conference on
Knowledge discovery and data mining, 137–146. ACM.
[Kocaoglu, Dimakis, and Vishwanath 2017] Kocaoglu, M.;
Dimakis, A. G.; and Vishwanath, S. 2017. Cost-optimal
learning of causal graphs. arXiv preprint arXiv:1703.02645.
[Koller and Friedman 2009] Koller, D., and Friedman, N.
2009. Probabilistic graphical models: principles and techniques. MIT press.
[Leskovec et al. 2007] Leskovec, J.; Krause, A.; Guestrin,
C.; Faloutsos, C.; VanBriesen, J.; and Glance, N. 2007.
Cost-effective outbreak detection in networks. In Proceedings of the 13th ACM SIGKDD international conference on
Knowledge discovery and data mining, 420–429. ACM.
[Marbach et al. 2009] Marbach, D.; Schaffter, T.; Mattiussi,
C.; and Floreano, D. 2009. Generating realistic in silico gene
networks for performance assessment of reverse engineering
methods. Journal of computational biology 16(2):229–239.
[Meek 1997] Meek, C. 1997. Graphical models: Selecting
causal and statistical models.
[Meinshausen et al. 2016] Meinshausen, N.; Hauser, A.;
Mooij, J. M.; Peters, J.; Versteeg, P.; and Bühlmann, P. 2016.
Methods for causal inference from gene perturbation experiments and validation. Proceedings of the National Academy
of Sciences 113(27):7361–7368.
[Minoux 1978] Minoux, M. 1978. Accelerated greedy algo-

rithms for maximizing submodular set functions. Optimization Techniques 234–243.
[Nemhauser, Wolsey, and Fisher 1978] Nemhauser, G. L.;
Wolsey, L. A.; and Fisher, M. L. 1978. An analysis of
approximations for maximizing submodular set functions?i.
Mathematical Programming 14(1):265–294.
[Pearl and Bareinboim 2011] Pearl, J., and Bareinboim, E.
2011. External validity and transportability: A formal approach. JSM Proceedings, Miami Beach FL 157–171.
[Pearl 2009] Pearl, J. 2009. Causality. Cambridge university
press.
[Peters and Bühlmann 2012] Peters, J., and Bühlmann, P.
2012. Identifiability of gaussian structural equation models
with equal error variances. arXiv preprint arXiv:1205.2536.
[Peters, Bühlmann, and Meinshausen 2016] Peters,
J.;
Bühlmann, P.; and Meinshausen, N. 2016. Causal inference
by using invariant prediction: identification and confidence
intervals. Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 78(5):947–1012.
[Schölkopf et al. 2012] Schölkopf, B.; Janzing, D.; Peters, J.;
Sgouritsa, E.; Zhang, K.; and Mooij, J. 2012. On causal and
anticausal learning. arXiv preprint arXiv:1206.6471.
[Shanmugam et al. 2015] Shanmugam, K.; Kocaoglu, M.;
Dimakis, A. G.; and Vishwanath, S. 2015. Learning causal
graphs with small interventions. In Advances in Neural Information Processing Systems, 3195–3203.
[Shimizu et al. 2006] Shimizu, S.; Hoyer, P. O.; Hyvärinen,
A.; and Kerminen, A. 2006. A linear non-gaussian acyclic
model for causal discovery. Journal of Machine Learning
Research 7(Oct):2003–2030.
[Shpitser and Pearl 2006] Shpitser, I., and Pearl, J. 2006.
Identification of joint interventional distributions in recursive semi-markovian causal models. In Proceedings of the
National Conference on Artificial Intelligence, volume 21,
1219.
[Shpitser, VanderWeele, and Robins 2012] Shpitser, I.; VanderWeele, T.; and Robins, J. M. 2012. On the validity of
covariate adjustment for estimating causal effects. arXiv
preprint arXiv:1203.3515.
[Spirtes, Glymour, and Scheines 2000] Spirtes, P.; Glymour,
C. N.; and Scheines, R. 2000. Causation, prediction, and
search. MIT press.
[Verma and Pearl 1991] Verma, and Pearl, J. 1991. Equivalence and synthesis of causal models. In Proceedings of
UAI, 220–227.
[Verma and Pearl 1992] Verma, T., and Pearl, J. 1992. An
algorithm for deciding if a set of observed independencies
has a causal explanation. In Proceedings of UAI, 323–330.

