Optimal Experiment Design for Causal Discovery from Fixed Number of Experiments
AmirEmad Ghassami Saber Salehkaleybar Negar Kiyavash


[
Optimal Experiment Design for Causal Discovery from
Fixed Number of Experiments





AmirEmad Ghassamighassam2@illinois.edu
Department of ECE, Coordinated Science Laboratory,

			University of Illinois at Urbana-Champaign Urbana, IL 61801 USA
Saber Salehkaleybarsabersk@illinois.edu
Coordinated Science Laboratory,

			University of Illinois at Urbana-Champaign Urbana, IL 61801 USA
Negar Kiyavashkiyavash@illinois.edu
Department of ECE and ISE, Coordinated Science Laboratory,

			University of Illinois at Urbana-Champaign Urbana, IL 61801 USA





boring formatting information, machine learning, ICML

0.3in
]




 
 
 

We study the problem of causal structure learning over a set of random variables when the experimenter is allowed to perform at most M experiments in a non-adaptive manner. We consider the optimal learning strategy in terms of minimizing the portions of the structure that remains unknown given the limited number of experiments in both Bayesian and minimax setting. We characterize the theoretical optimal solution and propose an algorithm, which designs the experiments efficiently in terms of time complexity. We show that for bounded degree graphs, in the minimax case and in the Bayesian case with uniform priors, our proposed algorithm is a  ρ-approximation algorithm, where ρ is independent of the order of the underlying graph. Simulations on both synthetic and real data show that the performance of our algorithm is very close to the optimal solution.






§ INTRODUCTION






Causal structures are commonly represented by directed acyclic graphs (DAGs), where the vertices of the graph are random variables and a directed edge from X to Y indicates that variable X is a direct cause of variable Y <cit.>.
Given a set of variables, there are two main approaches for uncovering the causal relationships among them.
First to perform conditional independence tests on the set of variables based on observational measurements <cit.>.
The second involves intervening on some of variables to recover their causal effect on the rest of the variables.
Unlike observational only tests, sufficient experiments involving interventions can uniquely identify the underlying causal graph completely.

In the study of intervention-based inference approach,
often a setup in which the experimenter performs M experiments on the set of variables is considered. In each experiment a set of at most k variables are intervened on.
In this setting, two natural questions arise:




  * What is the smallest required number of experiments in order to learn all the causal relations?

  * For a fixed number of experiments, what portion of the causal relationships can be learned?




The first problem has been addressed in the literature under different assumptions (see the Related Work subsection).
To the best of our knowledge, the second question has not been studied in the literature, and it is this question that we address herein.
Specifically, we consider a setup with M experiments, each containing exactly one intervention.
The reason we consider single-intervention experiments is that in many applications, such as some experiments in biology, simultaneous intervention in multiple variables may not be feasible.
For the underlying structure on the variables, we assume that the number of cycles of length three in the structure is negligible compared to the order of the graph. Structures satisfying this assumption arise in many applications.
For instance, causal structure in gene regulatory network (GRN) for some bacteria such as the Escherichia coli (E-coli) and Saccharomyces cerevisiae (S. cerevisiae) has a tree like structure (see Figure <ref>), and hence, satisfies our assumption <cit.>.

Contributions.
Unlike most of the previous work, we utilize a hybrid inference scheme instead of an adaptive approach (Subsection <ref>).
In this approach, first an observational test, such as the IC algorithm <cit.>, is performed on the set of variables. This test reveals the skeleton as well as the orientation of some of the edges of the causal graph. Next, based on the result of the initial test, the complete set of M experiments is designed in a non-adaptive manner. 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 Subsection <ref>.
We study the problem of structure learning for both Bayesian and minimax settings.
In Section <ref>, we present the optimal solution for both settings. This solution is optimal in the sense of recovering the structure that minimizes the loss for a given number of interventional experiments. Finding this optimal solution is in general computationally intense. In Section <ref>, we propose ProBal algorithm, which finds the set of experiments in a computationally efficient manner. We show that for bounded degree graphs, in the minimax setting and in the Bayesian settings with uniform prior, our proposed algorithm is a ρ-approximation algorithm, where ρ is independent of the order of the underlying graph.
In Section <ref>, using synthetic and real data, we show that the performance of ProBal is very close to the optimal solution.

Related Work.

The best known algorithms for general purely observational recovery approaches are IC <cit.> and PC <cit.>.
Such purely observational approaches reconstruct the causal graph up to Markov equivalence classes, and hence, the direction of some of the edges may remain unresolved.
Of course under some conditions, full causal structure learning using merely observational data is feasible <cit.>.

There is a large body of research on learning causal structures using interventional data <cit.>. Specifically,
Pearl <cit.>, considers the SEM model and defines so-called do-intervention to infer the causal relations among a set of variables.

A similar approach for representing interventions is adopted in <cit.>, but it allows interventions to have non-degenerate distributions.
Woodward <cit.> proposed another type of intervention which unlike Pearl's, does not depend on a specified model of the causal relations among the random variables.
Peters et al. introduced invariant causal prediction <cit.>, which is a causal discovery method that uses different experimental settings to predict the set of ancestors of a variable.
In that work, data comes from different unknown experimental settings (which could results from interventions).
See <cit.> for some validations of this method.

Regarding the first question discussed earlier, <cit.> consider the complete graph as the underlying causal structure to obtain the worst case bounds on the number of required experiments.
In that work, both cases of experiments containing bounded and unbounded number of interventions are studied.
In <cit.>, it has been shown that there is  a connection between the problem of finding a separating system in a graph and designing a proper set of experiments for causal inference, and hence, results from combinatorics help in finding the fundamental bounds.

In <cit.>, two algorithms that minimize the number of experiments in the worst case are developed. The proposed 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 <cit.>, 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 allows to learn chordal graphs completely.





§ MODEL DESCRIPTION





 §.§ Preliminaries


In this subsection we introduce some definitions and concepts that we require later.



Consider a directed graph D=(V,E) with vertex set V and set of directed edges E.
D is a DAG if it is a finite graph with no directed cycles.



A DAG D is called causal if its vertices represent random variables V={X_1, ...,X_n} and a directed edges (X_i,X_j) indicates that variable X_i is a direct cause of variable X_j.

We consider a structural equation model <cit.>, which is a collection of n equations X_i=f_i(PA_X_i,N_i), i=1,...,n, where PA_X_i denotes the parents of X_i in D, and N_i's are jointly independent noise variables. We assume here that in our network, all variables are observable.
Also, throughout the rest of the paper, we assume the faithfulness assumption on the probability distribution.












Two causal DAGs D_1 and D_2 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 <cit.>.



A v-structure is a structure containing two converging directed edges whose tails are not connected by an edge. v-structures are also known as immorality and a graph with no immorality is called a moral graph.





Using purely observational data  (referred to which as the null experiment by <cit.>), one can utilize  a “complete" conditional independence based algorithm to learn the causal structure as much as possible. By complete we mean that the algorithm is capable of distinguishing all the orientations up to the Markov equivalence class. Such an algorithm includes performing a conditional independence test followed by applying the Meek rules <cit.>.

On the other hand, interventions can enable us to differentiate among the different causal structures inside a Markov equivalence class. Define an intervention I on variable X∈ V(D) as removing the influence of all the variables on X and randomizing the value of this variable. We denote this intervention by I=X. An inference algorithm consists of a set of M experiments[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 randomization.] ℰ_total={ℰ_1, ℰ_2, ..., ℰ_M}, where each experiment contains k interventions, i.e., ℰ_i={I_1^(i), I_2^(i), ..., I_k^(i)} for 1≤ i ≤ M. As shown in <cit.>, observing the result of the null experiment,  one can find the orientation of the edge between any two variables X_i and X_j, if there exists ℰ_k∈ℰ_total such that
(X_i∈ℰ_k , X_j∉ℰ_k) or 
(X_j∈ℰ_k , X_i∉ℰ_k).


An inference algorithm may be (1) adaptive, in which case the experiments are performed sequentially and the information obtained from the previous experiments is used to design the next one; (2) passive, in which all the experiments are designed beforehand;
(3) hybrid, in which the experimenter first performs a pure observational study to obtain the skeleton and some of the orientations in the causal DAG, and then designs the rest of the experiments in a passive manner.
The third approach is referred to as the passive setup by <cit.>, while <cit.> use the term passive for a setting in which the interventions are selected without performing the null experiment.





 §.§ Problem Definition



Since in many practical applications, it is resource-intensive to perform the experiments sequentially, we will investigate the hybrid approach for the design of experiments. 
A nice feature of the hybrid approach is that it allows us to parallelize the performance of the experiments.

For example, in the study of GRNs introduced in Section <ref>, the GRN of all E-coli cells are the same and experiments can be performed on different cells simultaneously. 



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.





Let Z be the set of edges whose orientations are identified after the obseravtional test. Consider the moral graph D\ Z. As noted in <cit.>, D\ Z consists of a set of disjoint chordal moral graphs {G_1, G_2, ..., G_K} (see Figure <ref> as an example). To learn the structure, it suffices to learn {G_1, G_2, ..., G_K}. We focus on one such graph, say G. Hence, in the remainder of the paper, G is assumed to be a moral chordal DAG.


Let S be the set of all edges in a chordal graph which do not belong to any triangle.
We refer to a connected component of G\ S of order larger than one as “cyst”.




We assume that the number of triangles, compared to the order of the chordal graph is negligible. Formally, we assume that the proportion of the number of triangles to the order of G goes to zero as n=|V(G)| tends to infinity.


Assumption <ref> implies that both the number of cysts and their orders are negligible.

As mentioned in Section <ref>, structures with this property occur frequently in genomic and other applications.


As mentioned earlier, in some scenarios the experimenter may be restricted to perform a limited number of experiments.
Hence, we focus on the following problem: If the experimenter is allowed to perform M experiments, each of size k=1. What portion of the graph could be reconstructed on average and in the worst case?
 We shall distribute our total budget of M experiments over the components {G_1, G_2, ..., G_K} proportionally to their size. Thus, we assume that we are capable of performing m experiments ℰ={ℰ_1, ℰ_2, ..., ℰ_m}⊆ℰ_total on component G.



Consider the graph obtained by contracting all the vertices of a cyst into a single vertex. Clearly, this merging graph would be a tree,

and hence, we work on a moral tree structure after the contractions. If we are required to intervene on a vertex in the resulting tree that corresponds to a cyst, we will intervene on all the vertices of that cyst. This is in some cases necessary as shown in Figure <ref>. Note that after intervening in all the nodes of a cyst, we recover all the orientations inside the cyst. We emphasize that some edge orientations remain unidentified in the cysts that we have not intervened on; however, due to Assumption <ref>, the number of such edges will not scale with n. Also, due to Assumption <ref>, we can count intervening on a contracted cyst as one intervention without impacting the scaling results.


Root variable is a variable for which the number of edges entering that variable is zero.




A moral chordal graph G has only one root unless they all belong to the same cyst.


It suffices to prove that we cannot have two roots in two separate cysts. We prove this by contradiction.
Suppose roots r_1 and r_2 exist in two separate cysts. Since the component is assumed to be connected, there exists a vertex v where a directed path from r_1 to v meets a directed path from r_2 in v. Let the parents of v on these paths  be p_1 and p_2, respectively. If p_1 and p_2 are not adjacent, we will have a v-structure on the set {p_1,v,p_2} which is a contradiction. Otherwise, p_1 and p_2 are adjacent and {p_1,v,p_2} are in a cyst. Without loss of generality, assume p_1 is a parent of p_2. Let p_3 be the parent of p_2 on the path from r_2. In this case, we will have the same argument for the set  {p_1,p_2,p_3}, and we either have a v-structure, or p_3 is also in the cyst containing v. Continuing this argument, we either have a v-structure or r_1 and r_2 are in the same cyst as v, which is a contradiction.




In a moral chordal graph, randomizing the root suffices to fully learn the orientations of all the edges not belonging to any cysts.


Let r be the root vertex. Consider an edge e=(u,v).
Since e does not belong to any directed cycle, and there is a directed path from r to u and also a directed path from r to v, e should be directed from u to v, otherwise it contradicts Lemma <ref>.


Therefore, in the sequel, we focus on a moral tree structure and from Lemma <ref>, we know that randomizing the root variable will give us all the orientations.

We model this problem as follows. Let 𝒫={P^θ:θ∈Θ} be the set of probability distributions over the location of the root of the tree. We assume that the probability distributions of interest are all positive.

The following remark is a consequence of Lemma <ref>.


For each location of the root, only one moral tree is consistent with the skeleton of the tree. Therefore, 𝒫 is the set of probably distributions over the realizations of the moral tree.



Define U(ℰ) as the set of edges whose orientation is not found after performing the experiment set ℰ. For the given skeleton obtained from the observational test, let T_v be the moral tree of order n with vertex v as its root.
We define the loss of an experiment set ℰ on T_v as
l(ℰ,T_v)=|U(ℰ)|,
and the average loss of the experiment set under distribution P^θ as

    L_θ(ℰ)=∑_vP^θ(v)l(ℰ,T_v).


The problem of finding the best experiment set for the worst case, could be stated as the following minimax problem:

    min_ℰmax_θL_θ(ℰ),
    
    s.t. |ℰ|=m,
       |ℰ_i|=1   ∀ i: 1≤ i≤ m

In some real life applications, the experimenter may have prior knowledge about the possible location of the root in the tree. That is, the probability distribution P^θ∈𝒫 over the location of the root in the tree is known.  

In such a setting, we can investigate the following Bayesian version of the problem:




    min_ℰL_θ(ℰ),
    
    s.t. |ℰ|=m,
       |ℰ_i|=1   ∀ i: 1≤ i≤ m




§ OPTIMAL SOLUTION



Consider a moral tree T of order n.

In this structure, every non-root vertex has incoming degree d^-(v)=1,
the root vertex, r has incoming degree d^-(r)=0, and recall that intervening on the root identifies the whole tree.

Assume the experiment set ℰ={{I^(1)_1},...,{I^(m)_1}} was performed on the moral tree T. Let ℐ={I^(1)_1,...,I^(m)_1} be the set of variables on which we intervened. In this case, the subgraph T\ℐ will be a forest containing J components {C_1,...,C_J}.



Performing experiment set ℰ, on T_v rooted at v,

    l(ℰ,T_v) = {[           0        v∈ℐ,; |C_j|-|B_j|     v∈ C_j, ].

where C_j∈{C_1,...,C_J}, and B_j=C_j∩ N(ℐ), where, N(ℐ) denotes the set of neighbors of variables in ℐ.





After performing a single intervention I=X, the orientation of all the edges entering the descendants of X and the edges entering X itself will be recovered. Thus
l(X,T_v)=|ND_v(X)|-1_{X≠ v}.
As a result, intervening on variables in ℐ, if v∈ℐ, since |ND_v(v)|=0 the loss would be equal to zero; otherwise, if v∈ C_j we will have:

l(ℰ,T_v)=|⋂_X∈ℐND_v(X)
|-|⋂_X∈ℐND_v(X)∩ N(ℐ)|
=|C_j|-|B_j|.










The average loss could be bounded as follows



    ∑_j=1^JP^θ(C_j)|C_j|-m≤ L_θ(ℰ)≤∑_j=1^JP^θ(C_j)|C_j|-1,






Using Lemma <ref>, we have



    L_θ   (ℰ)=∑_v P^θ(v)l(ℰ,T_v)
       =∑_v∈ℐP^θ(v)×0
    +∑_j=1^J∑_v∈ C_jP^θ(v)(|C_j|-|B_j|)
       =∑_j=1^J(|C_j|-|B_j|)∑_v∈ C_jP^θ(v)
    
    
       =∑_j=1^JP^θ(C_j)|C_j|
    -∑_j=1^JP^θ(C_j)|B_j|.

Note that for all j, 1≤|B_j|≤ m, and the result is immediate.




Since we have assumed that m≪ n, in the sequel we will focus on minimizing

    L̂_θ(ℰ)=∑_j=1^JP^θ(C_j)|C_j|.



In the special case of uniform P^θ, the function L̂ could be obtained as
L̂_θ(ℰ)=1/n∑_j=1^J|C_j|^2.





Bayesian Setting. In the Bayesian setting, we seek the set of experiments that minimizes L̂_θ. This can be done by checking all nm possible vertex sets of size m. This brute-force solution is computationally intensive.

In Section <ref> we will introduce an efficient approximation algorithm instead.




Minimax Setting. As mentioned in Section <ref>, in the minimax setting we are interested in finding optimal ℰ in min_ℰmax_θL̂_θ(ℰ). Note that the loss is maximized if all the probability mass of root is put on the largest component. That is



    min_ℰmax_θL̂_θ(ℰ)
       =min_ℰmax_θ∑_j=1^JP^θ(C_j)|C_j|
       ≤min_ℰmax_j|C_j|.


Therefore, it suffices to choose a variable set ℐ={X_1,...,X_m}, which minimizes max_j |C_j|.

We can again use a brute-force solution which suffers from the same computational complexity issues as in the Bayesian case. We will propose an approximation algorithm for the minimax approach in Section <ref>.



§ EFFICIENT LEARNING ALGORITHM






In this section, we propose an algorithm for finding experiment sets efficiently for both Bayesian and minimax settings. First, we define the concept of a separator vertex, which plays a key role in the proposed algorithm.
We shall see that in our method, we require the prior on the location of the root variable.
Using separators, we propose the probability balancer (ProBal) algorithm, which allows for experiment design. The main idea is to iteratively decompose the tree into two subtrees, referred to as segments, sharing a separator vertex. In Subsection <ref>, we consider the minimax extension in which no prior is not available to the experimenter.







Let lobes of a vertex v in a tree be the remaining components in the graph after removing v.
A vertex v  in a tree T is a separator if the probability of the root being in each of its lobes is less than 1/2.



There exists a separator vertex in any tree.




Consider any arbitrary vertex v_1 in the tree. If all its lobes have probability less than 1/2, then it is a separator; otherwise, only one of its lobes, say lobe b, has probability larger than or equal to 1/2. Consider v_2, which is the neighbor of v_1 in b. Since the probabilities should add up to 1, the lobe connected to v_2 through v_1 should have probability less than 1/2 and hence we continue with checking the probability of the other lobes of v_2. This process will result in finding a separator because there are no cycles in a tree and we have assumed that the tree is finite.




The lobes of a separator vertex can be partitioned into two wings such that the probability of the root being in each of them is less than 2/3.




Suppose v is a separator vertex with lobes b_1,⋯,b_l sorted in ascending order of their probabilities  P^θ(b_i)=p_i. We also add lobe b_0 with p_0=0.
Let j be the largest index such that ∑_i=0^jp_i≤2/3.
If j=l, any arbitrary partitioning of the lobes is acceptable. Assume j<l.
 If ∑_i=j+1^lp_i≤2/3, then {{b_0,⋯,b_j},{b_j+1,⋯,b_l}} is the desired partitioning. Otherwise, we have ∑_i=0^j+1p_i>2/3 and ∑_i=j+1^lp_i>2/3, which implies that p_j+1>1/3 and since v is a separator, p_j+1<1/2. Therefore, P^θ({b_0,⋯ b_l}\{b_j+1})<2/3, and {{b_j+1},{{b_0,⋯ b_l}\{b_j+1}}} is the desired partitioning.







ProBal algorithm searches for m variables in a given tree T in an iterative manner.
It starts with the original tree as the initial segment and in each round it breaks it into smaller segments in the following manner. Let 𝒢 be the set containing all the segments.
At each round the algorithm picks the segment G_m with largest P^θ(G) in the set 𝒢, and finds the most suitable separator (described below) and adds it to the intervention set ℐ (if it is not already in ℐ).
This is done using the function FindSep.
 Then using the function Div, the algorithm divides G_m into two new segments G_1 and G_2, and replaces G_m with {G_1,G_2} in the set 𝒢 unless they have a star structure with the used separator as the center vertex. The reason that we ignore the star structures is that since the center is already chosen to be intervened on, the orientation of all the edges of the structure are discovered and there is no need for further interventions on the other variables in the star. The process of choosing separators continues until m variables are collected or all the graph is resolved. The set ℐ will be returned as the set of intervention variables.
The functions FindSep(·) and Div(·) are described below.


 

  * In the function FindSep(·), first we normalize the values of probability in the input segment to obtain distribution P^θ. For any variable X∈ G_m, we compute the probability of root being in the lobes of X, and partition the lobes into two wings W_1^*(X) and W_2^*(X) such that the probability of the root being in the wings is as balanced as possible. Define the unbalancedness of X as s(X) |P^θ(W_1^*(X))-1/2P^θ(W_1^*(X)∪ W_2^*(X))|. The function returns variables X^* with minimum s(X).

  * The function Div(·) outputs segments G_1=G_m\ W_2^*(X^*) and G_2=G_m\ W_1^*(X^*). In both segments, it sets the probability of X^* to zero.







(a) A leaf will not be chosen as the separator more than once.

(b) If the chosen separator X^* in G_m is not a leaf, and the segments G_1 and G_2 are produced, then
max{|V(G_1)|,|V(G_2)|}<|V(G_m)|.




(a) Suppose a leaf variable X is chosen as the separator in one of the rounds of the algorithm. Consequently, its probability will be set to zero in the segment containing it. Since X is a leaf, one of its wings will be empty, and hence after normalization of the probability in function FindSep(·), the wings of X will have probabilities 0 and 1, while any other variable X' in the segment with non-zero probability, will have wings that are balancing the measure 1-P^θ(X') (Note that all the variables in this segment cannot have zero probability, because the distribution is assumed to be positive and also all the other variables could not have been picked as separators before, otherwise the algorithm would not have kept this segment). Therefore, the function FindSep(·) will not choose X.

(b) Since X^* is not a leaf, |V(W_2^*(X^*))|≠ 0, and since |V(G_1)|=|V(G_m)|-|V(W_2^*(X^*))|, we have |V(G_1)|<|V(G_m)|. Similarly |V(G_2)|<|V(G_m)|.




ProBal algorithm runs in time O(n^3).




One can find the probability of each lobe of a vertex in linear time by running Depth-first search (DFS) algorithm. Therefore, FindSep(·) runs in time O(n^2). Additionally, Div(·) runs in o(n^2). From Lemma <ref>, ProBal algorithm will end in at most n rounds. Thus, the time complexity of the algorithm is O(n^3).






 §.§ Analysis


In this subsection we find bounds on the performance of ProBal algorithm. We will show that in the case of a uniform prior, the proposed algorithm is a ρ-approximation algorithm, where ρ is independent of the order of the graph.





After running ProBal algorithm for r rounds and obtaining the experiment set ℰ, the loss L̂_θ(ℰ) is upper bounded as



    L̂_θ(ℰ)≤(2/3)^⌊log_2(r+1) ⌋n.







First we claim that with r=2^k-1, the set 𝒢 defined in the algorithm which contain at most 2^k segments, has the property
max_G∈𝒢P^θ(G)≤(2/3)^k.
We use induction to prove this claim. The base of the induction is clear from Proposition <ref> and the fact that we set the probability of the separator itself to zero in function Div(·). For the induction step we need to show that after r=2^k+1-1 rounds max_G∈𝒢P^θ(G)≤(2/3)^k+1.
By the induction hypothesis with r=2^k-1 rounds max_G∈𝒢P^θ(G)≤(2/3)^k. Now, after the extra 2^k rounds, if a different segment were divided in each of those rounds, by Proposition <ref>, the desired result could be concluded. Otherwise, at least one of the segments, say G', was not divided while there exists another segment, say G”, which was divided more than once.
This implies that in the second dividing of G”, at least one of its sub-segments, say G_1”, obtained from the first division, had a larger probability than P^θ(G'). But by Proposition <ref>, we have P^θ(G_1”)≤2/3(2/3)^k. Therefore, P^θ(G')≤(2/3)^k+1.

Each component C_j belongs to a segment G∈𝒢. Hence, by the claim above, for all j, P^θ(C_j)≤(2/3)^k. This concludes

L̂_θ(ℰ)=∑_j=1^JP^θ(C_j)|C_j|

≤(2/3)^kn
≤(2/3)^⌊log_2(r+1) ⌋n.




In the following we prove that in the case of uniform prior on bounded-degree graphs,  ProBal is a ρ-approximation algorithm, where ρ is independent of the order of the graph. To this end, we first obtain a lower bound on the loss of the optimum algorithm introduced in Section <ref>.


Consider tree T of order n with maximum degree Δ(T) and uniform probability distribution over the location of the root of the tree. 
Then for any experiment set ℰ, we have

(n-m)^2/n(Δ(T)m-1)≤L̂_θ(ℰ).





For J components, from Corollary <ref>, the loss function may be lower bounded as follows

L̂_θ(ℰ)=1/n∑_j=1^J|C_j|^2≥1/nJ(n-m/J)^2.

Since the tree is connected, the maximum number of components created by an experiments set of size m is Δ(T)m-1, which implies the result.




Consider tree T of order n with maximum degree Δ(T).
In the case of uniform prior distribution, if m≤ϵ n, ProBal is a ρ-approximation algorithm, where
ρ=3/2(m∧ r)^log_22/3(Δ(T)m-1)/(1-ϵ)^2,
where m∧ r min{r,m}.

ρ is constant in n, polynomial of degree less than 0.42 in m, and linear on Δ(T).




Proof is immediate from Theorem <ref> and Lemma <ref>.






 §.§ Minimax Setting



As shown in (<ref>) in the minimax setting, we need to minimize the size if the largest component generated from the set of intervention variables.
If in ProBal, instead of balancing the probability of existence of the root variable in wings, we chose the separators to balance the number of vertices in the wings, then regardless of the probability distribution P^θ, it would guarantee that in each round of the algorithm, the largest segment G_m  would be divided into two segments, where each of these segments would contain at most 2/3|V(G_m)|+1 vertices. Therefore, at most the orientation of 2/3|V(G_m)| of the edges in G_m may not be found.
This is equivalent to running ProBal algorithm on a uniform distribution P^θ.



Consider tree T of order n with maximum degree Δ(T). Let L̂_W(ℰ) be the minimax loss of the experiment set ℰ.

(a)
After running ProBal for r rounds and obtaining the experiment set ℰ, the loss L̂_W(ℰ) is upper bounded as follows:
L̂_W(ℰ)≤(2/3)^⌊log_2(r+1) ⌋n.

(b)
For any experiment set ℰ, we have
n-m/Δ(T)m-1≤L̂_W(ℰ).

(c)
If m≤ϵ n, ProBal is a ρ-approximation algorithm, where
ρ=3/2(m∧ r)^log_22/3(Δ(T)m-1)/1-ϵ.




Since the performance of ProBal for the minimax case is the same as the Bayesian case with uniform distribution, the proof of part (a) follows from Theorem <ref>.
For part (b), since the tree is connected, the maximum number of components created by an experiments set of size m is Δ(T)m-1, and we can minimize the order of the largest, by making the orders equal. That is, the order of the largest component is at least n-m/Δ(T)m-1.
The proof of part (c) is immediate from parts (a) and (b).






§ EXPERIMENTAL RESULTS



In this section we evaluate the performance of our proposed algorithm on both synthetic and real data.



 §.§ Synthetic Data


We generated 1000 instances of trees based on Barabási-Albert (BA) model <cit.>, and bounded degree (BD) model created according to Galton-Watson branching process <cit.>. For both models we considered uniform and degree based distributions for the location of the root of the tree. In the degree based distribution, the probability of vertex v being the root is proportional to its degree.

Figure <ref> depicts the loss of ProBal as well as the loss for the optimal solution with respect to the order of the tree and the number of the interventions. As shown in this figure, in all cases, the performance of ProBal algorithm is very close to the optimal solution.
There are worst case scenarios where special graphs with specifically designed distributions can reach the upper-bounds, but as seen in Figure <ref>, for many distributions, ProBal algorithm works much better than what is predicted by our theoretic upper bound.



 §.§ Real Data

We examined the performance of ProBal on real data. The graph that we work on is the GRN of E-coli bacteria which we sourced from the RegulonDB database <cit.>.
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.


Figure <ref> depicts the normalized average loss of ProBal with respect to the total budget for the number of interventions. As seen in this figure, seven interventions are enough to reconstruct more than 95 percent of the network.




















§ 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 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 assumed that the underlying structure on the variables contains negligible number of triangles compared to the size of the graph. This assumption is satisfied in many applications such as the structure of GRN for some bacteria. 
We addressed the following question: “How much of the causal structure can be learned when only a limited number of experiments are available?” We characterized the optimal solution to this problem and then proposed an algorithm, which designs the experiments in a time efficient manner. We showed that for bounded degree graphs, in both the minimax setting and the Bayesian settings with uniform prior, our proposed algorithm is a  ρ-approximation algorithm, where ρ is independent of the order of the underlying graph. 
We examined our proposed algorithm on synthetic as well as real datasets. The results show that the performance of our proposed algorithm is very close to the optimal theoretical solution. One direction of future work is to extend this experiment design problem to even more general causal structures.












icml2016











§ PROOF OF THEOREM <REF>


First we claim that with r=2^k-1, the set 𝒢 defined in the algorithm which contain at most 2^k segments, has the property
max_G∈𝒢P^θ(G)≤(2/3)^k.
We use induction to prove this claim. The base of the induction is clear from Proposition <ref> and the fact that we set the probability of the separator itself to zero in function Div(·). For the induction step we need to show that after r=2^k+1-1 rounds max_G∈𝒢P^θ(G)≤(2/3)^k+1.
By the induction hypothesis with r=2^k-1 rounds max_G∈𝒢P^θ(G)≤(2/3)^k. Now, after the extra 2^k rounds, if a different segment is divided in each of those rounds, by Proposition <ref>, the desired result is concluded; otherwise, at least one of the segment, say G', is not divided while there exists another segment, say G”, which is divided more than once.
This implies that in the second dividing of G”, at least one of its sub-segments, say G_1”, obtained from the first division, had larger probability than P^θ(G'). But by Proposition <ref>, we have P^θ(G_1”)≤2/3(2/3)^k. Therefore, P^θ(G')≤(2/3)^k+1.

Each component C_j belongs to a segment G∈𝒢. Hence, by the claim above, for all j, P^θ(C_j)≤(2/3)^k. This concludes

L̂_θ(ℰ)=∑_j=1^JP^θ(C_j)|C_j|

≤(2/3)^kn
≤(2/3)^⌊log_2(r+1) ⌋n.

















