Recursive Causal Structure Learning in the Presence of Latent Variables and Selection Bias
Authors: Sina Akbari
Ehsan Mokhtarian
AmirEmad Ghassami
Negar Kiyavash
Venue: NeurIPS 2021 Poster (2021)
Keywords: Causal Discovery
Latent Discovery
Bayesian Network
Structure Learning
Latent Variables
Selection Bias
Removable Variable
Recursive Structure Learning

Abstract: We consider the problem of learning the causal MAG of a system from observational data in the presence of latent variables and selection bias. Constraint-based methods are one of the main approaches for solving this problem, but the existing methods are either computationally impractical when dealing with large graphs or lacking completeness guarantees. We propose a novel computationally efficient recursive constraint-based method that is sound and complete. The key idea of our approach is that at each iteration a specific type of variable is identified and removed. This allows us to learn the structure efficiently and recursively, as this technique reduces both the number of required conditional independence (CI) tests and the size of the conditioning sets. The former substantially reduces the computational complexity, while the latter results in more reliable CI tests. We provide an upper bound on the number of required CI tests in the worst case. To the best of our knowledge, this is the tightest bound in the literature. We further provide a lower bound on the number of CI tests required by any constraint-based method.  The upper bound of our proposed approach and the lower bound at most differ by a factor equal to the number of variables in the worst case. We provide experimental results to compare the proposed approach with the state of the art on both synthetic and real-world structures.

Recursive Causal Structure Learning in the Presence of Latent Variables and Selection Bias
Sina Akbari
Ehsan Mokhtarian
Department of Computer and Communication Sciences EPFL., Lausanne, Switzerland sina. akbari@epfl.
Department of Computer and Communication Sciences EPFL., Lausanne, Switzerland ehsan www.mmktaaaaanepll..cm
AmirEmad Ghassami
Negar Kiyavash
Department of Computer Science Johns Hopkins University, Baltimore, USA aghassa1 juu. edu
College of Management of Technology EPFL, Lausanne, Switzerland negar kiyavashhepfll ch
Abstract
We consider the problem of learning the causal MAG of system from observa- tional data in the presence of latent variables and selection bias. Constraint based methods of the main approaches for solving this problem, but the existing methods are either computationally impractical when dealing with large graphs or lacking completeness guarantees. We propose novel computationally efficient recursive constraint-bar approach is that each iteration specific type of variable is identified and removed. This allows us to learn the structure efficiently and recursively, as this technique reduces both the number of required conditional independence (CI) tests and the size of the conditioning sets. The former substantially reduces the computational complexity, while the latter results in more reliable CI tests. We provide an upper bound on the number of required CI tests in the worst case. To the best our knowledge, this the tightest bound in the literature. We further provide lower bound on the number of CI tests required by any constraint based method. The upper bound of our proposed approach and the lower bound most differ by a factor equal to the number of variables in the worst case. We provide experimental results to compare the proposed approach with the state of the art on both synthetic and real-world structures.
method that is sound and complete. The key idea
Introduction
Learning the causal structure among the of variables in the system the initial step for performing statistical inference tasks such as estimating the reward of policy in off-policy evaluation [23] [2]. etc. the literature, structure learning for the most part done under the assumption that all the
variables in the system are observed [21 I0 [7 I1l 23 However, in many applications in real-life
systems, this assumption violated. Moreover, the accessible data may contain selection bias, i.e..
some of the variables may have been conditioned on.
The problem of causal structure learning is significantly more challenging when unmeasured (latent) confounders and selection variables exist in the system. This is because the set of directed acyclic graphs (DAGs) independence models, which is the predominant modeling approach in the absence of unobserved variables, is not closed under marginalization and conditioning [19] That is, there does not necessarily exist a DAG over the observed variables that demonstrate one-to-one map with
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

U
W
T
Y
W
Y
W
Z
Z
z
(a) Gv
(b) Gv|O
(c) O (2) s]
(d) G (W]Ss)
Figure 1: Effect of removing variable the MAG over the remaining variables.
the conditional independence relationships in the observational distribution Pois, where O and S denote the observed variables and selection variables, respectively. To address this problem, several extensions of the DAG models, such acyelic directed mixed graphs (ADMGs) [18] induced path
graphs (IPGs) 21 and maximal ancestral graphs (MAGs) 19] are introduced in the literature.
The main approaches for structure learning include constraint-based and score-based methods [21 27
13 4
There also methods that require specific assumptions on the data generating modules, such as requiring linearity 30], linearity and Gon-Gaussianity of the noises [20] or additivity of the noise with specific types of non-linearity X (See
2 for a summary of structure learning approaches.)
Constraint based methods are the most commonly used methods for structure learning in the presence of latent variables and selection bias 21 24
[T The main idea in these methods is to find the structure which is most consistent with the conditional independence (CT) relationships in the data 21 However, the sheer number of tests required by these methods prohibits applying them to systems with large number of variables.
Several methods are proposed in the literature reduce the number of CI tests needed in constraint- based methods, specifically when there are no latent and selection variables in the system. For instance, 2I proposed the seminal PC algorithm for graphs with bounded degree, which has polynomial complexity in the number of vertices. ITO 17 I1 proposed using Markov boundary information to reduce the number of required CI tests. the size of the Markov boundaries or the in-degree of the variables is bounded, these methods achieve quadratic complexity in the number of the variables. However, the majority of the work on causal structure learning in the presence of latent and selection variables do not provide any analysis for the required number of CI tests. As an exception, for sparse graphs, and given the exact value of the maximum degree of the MAG as side information, [ proposed an algorithm that requires polynomial number of CI tests in the number of variables. Additionally, 4 proposed modification of the FCI algorithm, called RFCI, with specific attention to its time complexity. However, RFCI is not complete; that is, the output of this algorithm does not capture all the relationships the data.
In this paper, we propose novel recursive constraint- method for causal structure learning in the presence of latent confounders and selection bias. We use MAGs the graphical representation of the system. The main idea of our recursive approach is that in each iteration, we choose a particular variable of the system, say X, and locally identify its adjacent variables. Then, we recursively learn the structure over the rest of the variables using the marginal distribution Po Posssss. Note that the choice of X cannot be arbitrary. For instance, consider the DAG Gy in Figure Ta where the variable U is latent. The causal MAG over O = {T, W, Y, Z} is shown in Figure Ib as MAG M. As seen in Figure Tc| if we start with the choice of = T. we can correctly learn the subgraph of M over (Y, W, 2 J. whereas if we start with X = W, we will end up learning the graph in Figure Id which has two extra edges (highlighted in red) between Y, Z and Y, T that do not exist in M (we will revisit this example in Section
Our main contributions are as follows.
We introduce the notion of removable variable in MAGs, which is variable that can be removed from the causal graph without changing the m-separation relations (Definition F.. We further represent method to test the removability of a variable given the observational data (Theorem 2
We propose algorithm called L-MARVEL for causal structure learning in the presence of latent and selection variables. We show that our method sound and complete (Theorem B and performs On +n (M))24) 20))222040 ((M) (MV) CI tests in the worst case (Proposition 30, where denotes the number

---

of variables and (M is the maximum size of the union of parents, district and the parents of district of a vertex in the MAG (Equation (5)
We show that any constraint-based algorithm requires 22n2 + n (MM) CI tests in the worst case (Theorem D. Comparing this lower bound with our upper bound demonstrates the efficiency of our proposed method.
To sum up, the purpose and desirability of the proposed recursive algorithm for structure learning is two fold. First, since choose specific (removable) variables in each iteration (with the property of having small Markov boundary), we ensure that the number of required CI tests in each iteration, and hence in total, remains small. Therefore, we can significantly reduce the time complexity compared to non-recursive approaches. Second, by virtue of the gradual reduction of the order of the graph over the iterations, the size of the conditioning sets used in each CI test is reduced, which results in more reliable CI tests with smaller errors and more accurate results.
This paper is organized as follows. In Section we review the preliminaries, present the terminology, and formally describe the problem. In Section 3 we present the L-MARVEL method along with its analysis. In Section 3.3 we also provide the universal lower bound on the complexity of every constraint-basee based method. Finally, Section presents a comprehensive set of experiments to compare L-MARVEL with various algorithms synthetic and real-world structures.
2 Preliminaries and problem description
2.1 Terminology
mixed graph over the of vertices V is graph containing three types of edges - and ++. The two ends of an edge are called marks. There are two kinds of marks: arrowhead (>) and tail (-). If there exists a directed edge X + Y in the graph, we say X is parent of Y and Y a child of X. For bi-directed edge X + Y, we say X and Y are spouses. For an undirected edge X - Y. X and Y are called neighbors. In all of the aforementioned cases, we say X and Y are adjacent. The skeleton of G is an undirected graph with the same set of vertices V where there is an edge between X and Y if they are adjacent in G. A path from X to Y where every vertex on the path a child of its preceding vertex is called directed path. If a directed path exists from X to Y. X is called an ancestor of Y. We assume every vertex an ancestor of itself. We denote by Pa(X), (X), Ch(X), (XX.. X))(X), and Anc(X), ), the set of parents, children, neighbors, adjacent vertices, and ancestors of X, respectively. The district set of a variable X, denoted by (X), is the set of variables that have a path to X comprised of only bidirectional edges. By Pa*(X' (X) we denote the union of parents, district set, parents of district set, and the neighbors of a variable i.e.,
Pat (X)= = Pa(X) U Xi) U Pa(
UN(())
Uppercase capitals indicate single vertices, whereas bold letters denote sets of vertices. For set of vertices X. Anc(X) = Jx exAnc(X). ). A non- mon-endpoint vertex X on a path is called a collider, if both of the edges incident to X on the path have an arrowhead at X. A path P is collider path if every non-endpoint vertex on PP is a collider on P. A path P between the vertices X and Y is called an m-connecting active path relative to a set Z C (X Y), if (i) every non-collider on P is not Z).
a member of Z, and (ii) every collider on P belongs to Anc([XXY
Definition 1 (m-separation).. Suppose G is a mixed graph. A set Z m-separates X and Y in G. denoted by (X 1 g if there is m-connecting path between X and Y relative to Z in N We call Z a separating set for X and Y. We drop the subscript G whenever is clear from context.
A directed cycle exists in mixed graph if X Y and Y E Anc(X). An almost directed cycle exists in a mixed graph when X <> Y and Y E Anc(X). A mixed graph is said to be ancestral, if it does not contain directed cycles or almost-directed cycles, and for any undirected edge X Y, X and Y have no parents or spouses. An ancestral graph is called maximal if for any pair of non-adjacent
'The motivation behind this definition is that the local Markov property does not necessarily hold when causal sufficiency is violated, but if X Anc(Dis(X))) ())) then Pa tX) separates X from its non-descendants. See the supplementary material for proots.
'DAGs are subclass of mixed graphs. Note that for DAGs, this definition reduces to d-separation. See [4 for the definition of d d-separation.
3

---

vertices, there exists a set of vertices that m-separates them. A mixed graph is called a Maximal Ancestral Graph (MAG) if is ancestral and maximal. A MAG is called directed acyelic graph (DAG) if it has only directed edges.
A DAG 9 can be projected into unique MAG over subset of its vertices with the following projection, referred to as an embedded pattern in (28]
Definition 2 (Latent projection). Suppose a DAG over OUUUSS The projection of g
over O conditioned on S, denoted by ((O is a MAG over verfices O constructed as follows:
(i) Skeleton: X.Y E O are adjacent in [[[))] if there exists an inducing path $ in 9 between X and Y relative to S).
(ii) Orientation: For each pair of adjacent variables X, Y in G[O[S), the edge between X and Y is oriented as X - if E Anc([Y] US) and & Anc( [X] US); as X Y if X É Anc( ( US) and Y & Anc((X) S); and as X - Y if X E An) U U) and
Y E Anc( [xcc(X]uux US).
The above projection is the unique projection which satisfies the following property
19
(X g[Z) + (X YIZ U g)
Two MAGs are called Markov equivalent if they impose the same m-separations. A class of Markov equivalent MAGs can be represented as a (maximally informative) partially-oriented oriented ancestral graph (PAG), where the PAG contains the skeleton and all the invariant edge marks in the class.
Let P be the joint distribution over a of variables V. For X,Y V, S (XX,, Y), a conditional independence (CI) test in P on the triplet (X, Z, Y") yields independence, denoted by (X IL Y p. if P(X Y. Z) PZ) We drop the subscript P when it clear from context. Suppose 9 is DAG over V, i.e., each vertex of G corresponds to a variable of V. We say P is faithful with respect to G. if (X IL r + (X 1 Y [Z)g, g. i.e., the conditional independence in distribution / is equivalent to m-separation in the DAG g.
2.2 Problem description
We consider system with the set of variables V OULUS and the joint distribution Pv, where O, L. and S denote the set of observed, latent, and selection variables, respectively. Each variable X V is generated X fx (Pa( (X), where fx is deterministic function, Pa XX is the set of parents of X, i.e., the set of variables that have a direct causal effect on X. and ex is the exogenous noise corresponding to X. We assume all noise variables are jointly independent. This model is referred to as structural equations model (SEM) 15]. The causal graph of the system, which represents the causal relations among the variables, is denoted by Gy. Gy is directed graph over V. i.e., each vertex is associated with a variable and a directed edge exists from each variable in PX( to X, for all X E V. We assume that Gv is a DAG, and its latent projection over O conditioned on S is denoted by M == (G]())), We will call M the ground truth MAG. Further, we assume that Py is faithful with respect to Gv. which along with Equation ? implies that for each X,Y
Y E O
and Z S Ox(X
(X T Y = (X Z/Z)ro
(3)
Given the observational data from Pois. i.c., the marginal distribution the observed variables, conditioned on the selection variables, we consider the problem of learning the PAG that represents the Markov equivalence class (MEC) of M.
3 L-MARVEL Algorithm
In this section, we present Latent MARVEL (L-MARVEL.) algorithm to learn the PAG over O that represents the system. This algorithm relies notion similar to the MARVEL algorithm proposed by [II for DAG learning when all the variables are observable. Our approach relies on the Markoy boundary information as input.
3 An inducing path between X and Y relative to (L. S). where and S are disjoint sets not containing X and Y, is path on which every non-collider is member of I, and every collider belongs to Anc Y) US. We will use vertex and variable interchangeably throughout the paper.

---

Algorithm 1: L-MARVEL.
1: Input: O Pois
2: Output: PAG M
3: Mbo <- ComputeMb((., Pois)
4: A < Initialization O, Mbo)
5: A < L-MARVEL(O, Pois: Mbo. A)
6: Create Mt according adjacencies in A and orient it maximally using rules 0-10 of 31
1: Function L-MARVEL(V, sy Mby, A)
2: if [V] then
3: return A
4: else
5:
(X1, X2,.. (9,...X(v)) Xivi) t- Sort V based on the Markov boundary size in ascending order.
6:
for = --1n to V do
7:
(Adj(X;:, A) + FindAdjacent( X,, Mbv(X,). Pvis. A)
8:
isR <- -sRemovable(X; Mby (b,) Pvis, Ad)(X;)) % Main step of the algorithm.
9:
issR is true then
(Mbvx., A) UpdateMb(X, (X,, Adj(X;), Mbv, Pvis, A) return L-MARVEL(V (X,). Pr(x.))ss Mbr(xx.). A)
11:
Definition (Markov boundary). Suppose V S O. Markov boundary of X V with respect V is minimal set of variables Z C (XX), such that X is independent of the rest of the variables of conditioned on Z U.
Under faithfulness, Markov boundary of X E V O with respect to V. denoted by Mby (X), is
unique and it consists of all the variables that have collider path to X in Gv S (29 16 We indicate by Mby, the Markov boundaries of all of the variables in V with respect to V.
Throughout the algorithm, the data structure A stores the pairs of vertices that have been identified to be adjacent, and the separating sets found for non-adjacent vertices far. As the first step, the Markov boundary information with respect to O identified using one of the standard methods in the literature, as discussed in Section 3.2 Then A is initialized with separating sets implied by the Markov boundary information, i.e., for any X and Y & Mbo (X). Mbo (X) is a separating set for X, Y. is updated when a new separating set is discovered for a pair of vertices, two vertices determined to be adjacent. After initializing A in line 4, we call the L-MARVEL function over O, which recursively identifies all the adjacent pairs of vertices, i.e., the skeleton of M, and discovers separating set for all non-adjacent pairs of vertices. This information suffices to maximally orient the edge marks at the end of the algorithm using the complete set of orientation rules in 3
Our learning procedure is outlined in Algorithm T
L-MARVEL works as follows. It chooses variable X, identifies Adj(X) (i.e., the set of variables adjacent to X). and then recursively learns the structure over O (X). discarding X. This is desirable as the problem size decreases at each iteration, which results in substantial reduction in the computational complexity. Morcover, performing CI tests of high order is avoided. If the learned graph, i.c., Gv(O) \ {X) [S] is equal to the induced subgraph of Gy [ over O (X we can add X to this graph and connect it to its adjacent variables with edge. As discussed in the example in Figure []] this is not true for an arbitrary X. We show that for certain vertices, called removable, indeed apply such a recursive learning procedure. Next, define what makes variable
removable.
Definition 4 (Removable). Suppose is a MAG over V, X E V, and H is the induced subgraph of over V {X}. X is a removable vertex in i and H impose the same m-separation relations over V {X}. That is, for any vertices Y, W E V (XX and Z S V VXX,Y,WW,
(Y 1 Z
(Y L
(4)
In the case that G is DAG, Definition 4 reduces to what [ proposed for DAGs. However, their tests for identifying removability fail when causal sufficiency violated. Next, we provide graphical characterization of removable variables.
5

---

V2
X
m
Z
(a) Condition 1
(b) Condition 2
Figure 2: Graphical characterization of a removable variable. The edge marks indicated by a star can be either tail or am arrowhead.
Theorem 1. Vertex X is removable in a MAG M over the variables V, if and only if
1. for any Y Ad)(X) and Z Ch(X) \ (Y), Y and Z are adjacent, and
2. for any collider path u (X, Vi..., Vm, Y) and Z V VX Vm3 such that (x. (X., V _ Pa(Z), Y and Z are adjacent.
Figure represents the graphical constraints of Theorem Figure 2a depicts the first condition, where 7 is either child or a neighbor of X, and Y E X while Figure 2b depicts a collider path where X and Vis are parents of Z. Theorem [ states that is removable if and only the edges highlighted in red are present in both cases. See Appendix A formal proof and further discussion on Theorem [ The next proposition clarifies why removable variables are exactly those that can be removed at each iteration in our recursive approach.
Proposition 1. Suppose S and X E V. ( is equal to the induced subgraph of GV [V over ] if and only if X is removable in Gv (V)))
Appendix B includes the proofs of our results. Identifying removable variable at each iteration is the core of L-MARVEL. We will discuss an efficient algorithm to determine whether a variable is removable in Section T At each iteration, given the set of remaining variables V, these variables are sorted in ascending order of their Markov boundary size. Starting with the variable with the smallest Markov boundary, we search for its adjacent vertices within its Markov boundary. If Y E (by(X;) is not adjacent to X,, then X; and have a separating set in (b,))x)(()) (T6] Hence, identifying Ad)(X,) can be performed using a brute-force search in the Markov boundary, using most [Mb (x:v(( (2) (X.)) CI tests. In Section 3.3 we show that the loop in line 6 of Algorithm
D never reaches variables with large Markov boundaries, and this guarantees that both the number of CI tests and their order remains small throughout the algorithm. We then determine whether X; is removable given Mby (X,) and AX), using the efficient algorithm we shall discuss Xx,
in Section 3.1 We continue this procedure until identify the first removable variable Then, we remove X from the set of remaining variables, and update the Markov boundaries with respect V (X). which is the input to the next iteration. The latter does not require the discovery of Markov boundaries from scratch, and is implemented as we shall see in Section 3.2
The rest of the section is dedicated to showing how efficiently identify removable variable (Section 3.I. how to update the Markov boundary information (Section 3. 3.2h. and the analysis of the algorithm (Section
3.1 Testing removability in MAGs
The following theorem presents the conditions of removability of a variable using CI tests within the Markov boundary. This theorem excludes a particular structure of the MAG M. where M has a cycle of the length of at least four that contains only undirected edges, and this cycle has no chords. We shall discuss Appendix [ which specific structure of the DAG Gy this MAG represents, and why it is required to exclude this specific structure. As we shall see in Appendix D such MAGs imply very restrictive structure over the selection variables.
Theorem 2. Suppose the edge-induced subgraph of M over the undirected edges (i.e., the edges due to selection bias) is chordal. Let G = Gv V S) for some V O. X E V is removable in G if and only if for every Y E XX) and Z E Mby (X), at least one of the following holds.
Condition 1: 3W C Mby (Y, ZJ: Y Z(W.

---

Condition 2: MA C Mbb Y + 2 (X).
Furthermore, the set of removable vertices in G is non-empry.
Using Theorem and given Xdj(X) and (X)), Algorithm tests the removability of X by perform- ing ([AAjj( X) Mbv(X)) 2(Mbv(X)) CI tests. Note that the removability test only performed for variables with small Markov boundaries, which keeps both the number of CI tests and the size of the conditioning sets small, as we shall see in Section 3.3
Algorithm 2: IsRemovable Determine whether X is removable.
1: Input: (X, Mbv (X), Pvis, Adj(X))
2: for Y Ad)(X), Z E (X) do
3: if' Condition of Theorem 2 does not holds then
4: if Condition of Theorem 2 does not holds then
5: Return False
% X is not removable. % X is removable.
6: Return True
Conditions of Theorem 2 can be checked in different orders, although we have witnessed in our experiments that checking these conditions in the order of Algorithm increases the accuracy.
3.2 Markoy boundary discovery and updating Markov boundaries
L-MARVEL requires Markov boundary information for initialization. Several algorithms have been proposed the literature for discovering the Markov boundaries [0] 17 26 29]. For instance, TC 17 algorithm states that
(X M. (V),
X E MYv(Y) and Y E Mbv(X),
where V C O. Grow-Shrink (GS) algorithm ITC and its modifications, including IAMB and its variants (26) address Markov boundary discovery by performing more CI tests with smaller conditioning sets. These algorithms require linear number of CI tests in the number of variables to determine the Markov boundary of certain variable, i.e., quadratic number of CI tests to discover the entire Markov boundaries. However, given the challenging nature Markov boundary discovery. these algorithms might fail to accurately discover this information in some settings. We need to utilize one of these methods to initially discover the Markov boundaries, but the subsequent update of the boundaries throughout the later iterations is performed within L-MARVEL as we shall next discuss.
Updating Markoy boundaries: Let Mby be the input to an iteration of L-MARVEL where X is identified as removable and we need learn Mbv By definition of removability, the latent projection of Gy over V (X{) } is the induced subgraph of Gy V [V As result, removing X has two effects: 1. X removed from all Markov boundaries, and for Y, 2 E (X). if all of the collider paths between Y and Z in GV S pass through X, then Y and Z must be excluded from each others Markov boundary. Note that in the latter case, Y, Z E (Xb. The latter update is performed using Z)), or equivalently, (Y J.. Z (bb Y, 2)).
single CI test, i.e., (Y J.. Z| Mbbx(2)((X, We choose the CI test with the smaller conditioning set among the two. If the dependency does not hold, we remove Y, Z from each other's Markov boundary.
3.3 Analysis
First, we state the soundness and completeness L-MARVEL the following theorem.
Theorem 3. Suppose the distribution Pv over V O S is faithful to the DAG Gy. If the conditional independence relations among all variables in given S is provided to L-MARVEL, the output of L-MARVEL is the PAG representing the Markov equivalence class of [O]S))
Let + (H) denote the maximum size of Pat (-) (defined in I in MAG H, i.c.,
A (H)
max XEM
Pat (X)]-
(5)
In experiments, used TC.

---

Next, we provide an upper bound on the size of the Markov boundary of a removable variable. Proposition 2. If X is a removable variable in MAG H with vertices V, then (bv(()) - A (H).
L-MARVEL processes variables in the ascending order of their Markov boundary size at each iteration and stops when the first removable variable is identified. Therefore, Proposition 2 guaranteess that all the processed vertices each iteration have Markov boundaries smaller than A in (V))V) where V is the set of remaining variables. This number gets smaller during the algorithm. as L-MARVEL keeps only subgraph of the imput. Note that this bound applies to the size of the conditioning sets of CI tests performed in functions Aijcenj and IsRemovable, since the conditioning sets are subset of the Markov boundary. Furthermore, it results in the following upper bound on the number of CI tests.
Proposition 3. The number of conditional independence tests Algorithm performs on a MAG M of order n, in the worst case, upper bounded by
On + n (M))
(6)
The quadratic term in the upper bound of Equation 6 that algorithms such as GS, TC, IAMB, etc. discover the Markov boundary of each variable requiring only linear number of CI tests in n2
is for initial Markov boundary discovery. Note
To the best of our knowledge, this is the tightest bound in the literature. The following lower bound on all of the constraint based structure learning algorithm demonstrates the efficiency of L-MARVEL.
Theorem 4. The number of conditional independence tests of the form (X 11 Z required by any constraint- algorithm on a MAG M of order n, in the worst case, is lower bounded by
Din² + n (M))
Comparing this lower bound with our achievable upper bound, we can see that the complexity of L-MARVEL the worst case is merely different by a factor which is at most the number of observed variables in the worst case.
4 Experiments
We report empirical results on both synthetic (random graphs) and real-world structures available in the Bayesian network repository the benchmark for structure learning in the literature. We evaluate and compare L-MARVEI to various algorithms, namely the constraint-
methods FCI (21], 27 in terms of both computational
RFCI 4 and MBCS* [LG and the hybrid method M3HC
complexity and accuracy. Following the convention in [ 27 16 22 the data is generated according linear SEM with additive Gaussian noise, where all the variables of the system (including latent and selection) are generated as linear combinations of their parents plus a Gaussian noise. For each system, we simulate data from Po,s, the data available to all the algorithms. We use TC [I7] algorithm to learn the initial Markov boundaries, To make a fair comparison among the algorithms, we feed the Markov boundary information to all the algorithmss. that is, algorithms start from graph where the edges between vertices that are not in each other's ' Markov boundary are already deleted. This is similar to the ideas in 17 For CI tests, we use Fisher Z-transformation 6 with significance level a = 0.01 for all the algorithms, and a = 2/n² for TC [7]. In all the experiments, each point on the plots and each entry of the table represents an average of 50 datasets, where the latent and selection variables were chosen uniformly at random.
Random Structures: We used two different generating processes to obtain our random graphs. 1. MAGs corresponding to DAGs generated by Erdos-Renyl model G(ñ, p)
[5]. where ň denotes the total number of the variables, and 2. MAGs corresponding to random DAGs where each vertex
has maximum of 3 or 4 parents, similar to the setting in [27 I 2 Figures Ba and 3b illustrate the performance of the algorithms on Erdos Renyi graphs, whereas Figures 3c 3d
and Be represent
the performance of these methods on the latter generative model. The coefficients of the linear
model and the standard deviation of the exogenous noises are chosen uniformly at random from
buleemm. cow.aam.coommbbppsttoorrr///
Thh implementation of L-MARVEL is available in github Mokhaaarnnnntttaiaaaaaaaaaaaaa MARVEL
8

---

200
(a) MAGs corresponding G( 60.3 and latent 10%.
(b) MAGs corresponding Cū, 7..
and latent rate 20%
Grsph it
(c) MAGs corresponding to random DAGs with maximum of parents for each variable, latent rate 10%.
Otack Drsermo
Glack Ovder
Grech Order (n
(d) MAGs corresponding to random DAGs with a maximum of parents for each variable, latent rate = 20%
Dpni Spt
1.0 2011
Grach
(e) MAGS corresponding to random DAGs with a maximum of parents each variable, latent rate 10%.
L-MARVEL
M3HC
Figure 3: Performance of various algorithms random graphs with significance level 0.01 and 50]O] samples are available. Figures (a) and (b) demonstrate the evaluation over MAGs corresponding Erdos-Reny graphs, while (c), (d) and (e) represent the MAGG corresponding random DAGs with bounded number of parents for each variable (sample size 50]OO).
(0.5,, 2) and (1, V3), respectively. We did not continue running algorithms that not capable of keeping up with the cohort as the order the graphs grew. Moreover, the runtime of M3HC is not reported in the plots as it does not fit into the scale of the plots. As seen in the plots, L-MARVEL demonstrates substantially lower computational complexity in terms of number of the CI tests and runtime compared to the other algorithms, while maintaining high accuracy (the highest among the cohort in most of the cases). We also observed a low size for the conditioning sets in our CI tests for L-MARVEL. Only RFCI performs better than L-MARVEL in this metri
Benchmark Structures: Algorithms are evaluated on benchmark structures, where 5% to 10% of the variables are assumed to be latent, and 5% of them are selection variables. Latent and selection variables are chosen uniformly at random for each dataset. The coefficients of the linear model and
Note that RFCI avoids performing too many CI tests but with the caveat of lacking completeness.

---

the standard deviation of the noises are chosen uniformly at random from (0.5, 1) and (v0.5, 1). respectively. Our experiments, summarized in Table demonstrate that L-MARVEL outperforms the other algorithms both in terms of computational complexity and the accuracy of the learned structure. NA entries for FCI demonstrate that the runtime exceeds a certain threshold.
Table 1:
Performance of various algorithms on the benchmark structures, when 5% to 10% of the variables are latent and ~ 5% of them are selection variables (sample size 50(O)).
Structure
Insurance
Alarm
Ecoli
Barley
Hailfinder
Carpo
Arth1.0
(10), |C (S))
(22,3,2)
(31,4,2)
(40.3.3)
(40,5.3)
(50.3.3)
(53,4,4)
(95,6,6)
RFCI L-MARVEL
#CI tests
272
235
227
894
333
569
1185
Runtime
0.03
0.04
0.05
0.16
0.07
0.12
0.36
Fi-score
0.85
0.88
0.82
0.92
0.97
0.89
Precision
0.97
0.98
0.97
0.98
0.98
0,99
0.99
Recall
0.76
0.87
0.81
0.72
0.87
0.96
0.82
#CI tests
947
981
4314
2158
256754
11670
2644794
Runtime
0.14
0.20
0.86
0.44
62.22
2.59
1047 44
F1-score
0.76
0.89
0.85
0.73
0.88
0.94
0.87
Precision
0.99
1.00
1.00
1.00
1.00
1.00
1.00
Recall
0.63
0.81
0.74
0.58
0.79
0.89
0.77
#CI tests
7117
6899
56781
117566
NA
123198
NA
Runtime
1.13
1.25
13.22
25.78
NA
31.41
NA
FCI
Fi-score
0.75
0.88
0.83
0.70
NA
0.45
NA
Precision
0.99
1.00
1.00
1.00
NA
0.48
NA
Recall
0.61
0.80
0.72
0.54
NA
0.42
NA
#CT tests
640
335
499
2649
502
1221
3225
MBCS*
Runtime
0.12
0.11
0.17
0.77
0.19
0.46
1.94
F1-score
0.80
0.90
0.86
0.76
0.89
0.96
0.86
Precision
0.98
0.98
0.98
0.99
0,99
0.99
0.99
Recall
0.68
0.84
0.77
0.62
0.82
0.94
0.76
#CT tests
896
674
3033
1731
139788
8354
793754
M3HC
Runtime
13.66
4.19
6.64
12.53
47.72
7.42
322.33
F1-score
0.75
0.87
0.84
0.71
0.86
0.92
0.84
Precision
0.99
1.00
1.00
0.99
1.00
1.00
0.99
Recall
0.62
0.78
0.73
0.56
0.77
0.85
0.74
More comprehensive experimental results including the effect of the sample size, wider range of latent and selection rates, and assessments different settings of parameters and structures, along with alternative metrics are reported in Appendix
Concluding Remarks
We proposed a recursive structure learning approach capable of handling latent and selection variables. The recursive technique significantly reduced the number of required CI tests (and hence the time complexity). Also, since the order of the graph becomes smaller over the iterations, the recursive approach reduces the size of the conditioning performance of the tests. We provided an upper bound on the complexity of the proposed method as well as a lower bound for any constraint-based method. The upper bound of our proposed approach and the lower bound at most differ by a factor equal to the number of variables in the worst case, which demonstrates the efficiency of the proposal. We compared the performance the proposed method with several oftte-of-the-an approaches on both synthetic and real-world structures. The results showed improvement both performance and complexity on almost all the setups. We note that the performance of the proposed method is reliant on the accuracy of the Markov boundary information that used in the algorithm, Devising efficient and high accuracy approaches for learning the Markov boundary of the variables is left as important direction for future work.
in each CI test, which leads an improved
10

---

Acknowledgments and Disclosure of Funding
The work presented in this paper was in part supported by Office of Naval Research (ONR) under grant number W91 1911NF-15-1-0479
References
[1] Jean RS Blair and Barry Peyton. An introduction to chordal graphs and clique trees. In Graph theory and sparse matrix computation, pages 1-29. Springer, 1993.
[2] Kiattikun Chobtham and Anthony Constantinou. Bayesian network structure learning with causal effects in the presence of latent variables. arXiv preprint wwXi:2205./4338. 2020.
[3] Tom Claassen, Joris Mooij, and Tom Heskes. Learning sparse causal models is not np-hard. arXiv preprint
2013.
[4] Diego Colombo, Marloes H Maathuis, Markus Kalisch, and Thomas S Richardson. Learning high-dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics, pages 294-321, 2012.
[5] Paul Erdös and Alfréd Rényi. On the evolution of random graphs. Publications of the Mathe- maticul Institute of the Hungarian Academy of Sciences, 5:17-61, 1960.
[6] Ronald A Fisher. Frequency distribution of the values of the correlation coefficient in samples from an indefinitely large population. Biometrika, 10(4):507-521, 1915.
[7] Delbert Fulkerson and Oliver Gross. Incidence matrices and interval graphs. Pacific journal of mathematics, 15(3):835-855, 1965.
[8] Patrik O Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in neural information processing systems, pages 689-696 2009.
[9] Nan Jiang and Lihong Li. Doubly robust off-policy value evaluation for reinforcement learning. In International Conference on Machine Learning, pages 652-661. PMLR, 2016.
[10] Dimitris Margaritis and Sebastian Thrun. Bayesian network induction via local neighborhoods. Advances Neural Information Processing Systems, 12:505-511, 1999.
[11] Ehsan Mokhtarian, Sina Akbari, AmirEmad Ghassami, and Negar Kiyavash. A recursive markov boundary-based approach to causal structure learning. In The KDD''1 Workshop Causal Discovery, pages 26-54. PMLR, 2021.
[12] Susan A Murphy. Optimal dynamic treatment regimes. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65(2):331-355 2003.
[13] Juan Miguel Ogarrio, Peter Spirtes, and Joe Ramsey. A hybrid causal search algorithm for latent variable models. In Conference Probabilistic Graphical Models, pages 368-379. PMLR, 2016.
[14] Judea Pearl. Probabilistic reasoning in intelligent systems: Networks of plausible inference. Morgan Kaufmann, 1988.
[15] Judea Pearl. Causality. Cambridge university press, 2000.
[16] Jean-Philippe Pellet and André Elisseeff. Finding latent causes in causal networks: an efficient approach based on markov blankets. Neural Information Processing Systems Foundation, 2008
[17] Jean-Philippe Pellet and André Flisseeff. Using markov blankets for causal structure learning. Journal of Machine Learning Research, 90ul):1295-1344 2008
[18] Thomas Richardson. Markov properties for acyelic directed mixed graphs. Scandinavian Journal of Statistics, 355 2003.

---

[19] Thomas Richardson, Peter Spirtes, et al. Ancestral graph markov models. The Annals of Statistics, 30(4):92--1030, 2002.
[20] Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, and Antti Kerminen. linear non gaussian acyelic model for causal discovery. Journal of Machine Learning Research, 7(0ct):2003-2030, 2006.
[21] Peter Spirtes, Clark Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
[22] Eric V Strobl. A manste algorithm for causal discovery with cycles, latent variables and selection bias. International Journal of Data Science and Analytics, 8(1):33-56, 2019.
[23] Philip Thomas and Emma Brunskill. Data-efficient off- off-policy policy evaluation for reinforce- ment learning. In International Conference on Machine Learning, pages 2139-2148. PMIR, 2016.
[24] Robert E Tillman, David Danks, and Clark Glymour. Integrating locally learned causal structures with overlapping variables. In NIPS, pages 1665-1672. Citeseer, 2008.
[25] loannis Tsamardinos, Constantin Aliferis, and Alexander Statnikov. Time and sample efficient discovery of markov blankets and direct causal relations. In Proceedings of the ninth ACM SIGKDD international conference on Knowledge discovery and data mining, pages 673-678, 2003.
[26] Ioannis Tsamardinos, Constantin F Aliferis, Alexander R Statnikov, and Er Statnikov. Al- gorithms for large scale markov blanket discovery. In FLAIRS conference, volume 2, pages 376-380, 2003.
[27] Konstantinos Tsirlis, Vincenzo Lagani, Sofia Triantafillou, and Ioannis Tsamardinos. On scoring maximal ancestral graphs with the max-min hill climbing algorithm. International Journal of Approximate Reasoning. 102:74-85, 2018.
[28] Thomas Verma and Judea Pearl. Equivalence and synthesis of causal models. UCLA, Computer Science Department, 1991.
[29] Kui Yu, Lin Liu, Jiuyong Li, and Huanhuan Chen. Mining markov blankets without causal sufficiency. IEEE transactions on neural networks and learning systems, 29(12):6333-6447 2018.
[30] Chi Zhang, Bryant Chen, and Judea Pearl. A simultaneous discover-identify approach causal inference in linear models. In Proceedings of the AAAI Conference Artificial Intelligence, volume 34, pages 10318-10325, 2020.
[31] Jiji Zhang. Causal reasoning with ancestral graphs. Journal of Machine Learning Research, 1437-1474, 2008.
[32] Kun Zhang, Bernhard Schölkopf, Peter Spirtes, and Clark Glymour. Learning causality and causality-related learning: Some recent progress. National Science Review, 5(1):26-29, 2017.
12


Reviews:
[{"review_id": "j8Qv9QMnynF", "reviewer": "Reviewer_kohp", "summary": "The authors propose a method for learning a MAG from observational data in the presence of latent variable and selection, a problem addressed previously by FCI, RFCI, etc. The procedure works by identifying a removable variable, removing it, and so on, which renders the problem simpler and simpler as the algorithm proceeds. One benefit of the procedure is that it reduces the number of conditional independence tests needed for the search.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I found the approach nice, of identifying removable variables and then removing them one by one, as a way of optimizing the problem. The theory was clear.\n\nI found the paper to be original and useful. I've never seen this approach before. The quality of the theory was good. The writing was fine.\n\nI did have an issue with the experimental section. Table 1. It seems, though I can’t find text to support it, that F1, precision, and recall in this table are for adjacency. Is this right? If so, one wonders about orientation--could those stats be included?\n\nAlso, for Table 1, what packages were used for the various algorithms? I couldn’t find this information, though maybe it was there.\n\nAgain, Table 1, what graph were the results being compared to, and was it the same in each case? The true DAG? The true PAG?\n\nIn line 298, an especially wide interval about zero was excluded for random selection of coefficients. Is this necessary?\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I don't see how this could have negative societal impact.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "d3JHIsji9vc", "reviewer": "Reviewer_TRky", "summary": "The authors consider the problem of causal structure learning when there are latent variables and selection bias (i.e., estimating a MAG from observational data). While this is a well studied problem, most existing procedures can be computationally expensive when the number of variables is large and perform a large number of conditional independence tests which can suffer from statistical inefficiency. The authors propose a procedure which only considers subsets of the observed variables by checking whether certain variables are \"removable.\" This results in less computation and also increases statistical efficiency by generally requiring smaller conditioning sets when checking conditional independencies. The authors show that this improves estimation accuracy (as measured by F1 scores) and uses less computation.     ", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I think the problem that the article addresses is interesting and the approach is novel. Overall, I think this is a nice contribution as the procedure retains the same theorettical guarantees as previous methods, but results in improved empirical performance. However, since the authors give the same theoretical guarantees as FCI+ (i.e., sound and complete), I think it is quite important to compare RLCD to FCI+.\n\nMajor comments:\n- It's not immediately obvious to me that every $$V \\subseteq \\mathcal{O}$$ contains a removable node, but it seems that Line 257 requires this. Prop 2 requires X to be removable in order for the upper bound to hold. Therefore, if there are no removable nodes, line 257 may not hold.   \n-Although they present simulations which show that the RLCD procedure is drastically reduces computation / the number of CI tests, I am surprised they did not consider FCI+ as a competitor in the simulations, despite citing the paper. FCI+ is also sound and complete and is polynomial in |V| when an upper bound on the sparsity is given (an exact number is not required as the current authors state). This is not too different from RLCD which requires an initial estimate of the Markov boundary (who's size could be used as an upper bound on the degree). Also, comparing the theoretical results could be helpful. In particular, FCI+ states the computational cost in terms of maximum degree k, which in some cases could be smaller than the size of Pa+\n\nMinor comments:\n- Perhaps I'm missing something in Figure 1 (c), but when T is removed, shouldn't this result in a bidirected edge between W and Z? For instance, in (a), conditioning on W, the path Z, T, W, U, Y d-connects Z and Y and implies that Y \\not \\perp Z | W. However in C, Y \\perp Z | W. \n- Pg 1, line 74: it would be helpful to define n here since it is not defined until later\n- Pg 5, line 181: \"substantially decrease\" should be \"a substantial decrease\"\n- Pg 5, line 201: Should all the X's be X_1?    \n- Pg 6, line 230: It seems that the LHS of the equation should be \\not \\perp", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "NA", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "6cz48g48diA", "reviewer": "Reviewer_r1m9", "summary": "A causal discovery method is presented for MAG learning based on removing an appropriately chosen variable, then recursively learning the structure over the remaining variables. Also, a lower bound is given for the number of CI tests that any constraint-based method for this problem needs, and the number of CI tests by the proposed method compares favorably to this bound.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The algorithm is based on interesting ideas that are put together well (how to approach the structure learning problem recursively, and how to use Markov boundaries to pick the right node to recurse on). Additionally, the paper is well-written. I am not sufficiently familiar with constraint-based algorithms for MAG learning to evaluate the theoretical results directly, or to confirm their novelty in comparison with other recent work in this field (but as far as I can tell, this approach to MAG learning is novel). But the empirical validation looks properly performed, and its results are very positive. The reported improvements in running time as well as accuracy make this a significant result.\n\nUPDATE: I have read the authors' response, and my score is unchanged.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors have proposed a recursive algorithm for learning the structure of a Bayesian network in the presence of latent variables, based on the concept of a removable variable. All of the reviewers agreed that the approach is novel and interesting, and there is a consensus to recommend acceptance.
Decision: Accept (Poster)