Learning Directed Acyclic Graphs From Partial Orderings

Consider a DAG $G=(V,E)$ with the node set $V$, and the edge set $E\subset V\times V$. For the general problem, we assume that $V$ is partitioned into $L$ sets, $V_{1},\ldots,V_{L}$ such that for any $\ell\in\{1,\ldots,L\}$, the nodes in $V_{\ell}$ cannot be parents of nodes in any set $V_{\ell^{\prime}},\ell^{\prime}<\ell$. Such a partition defines a layering of $G$ (Manzour et al., 2021), denoted $V_{1}\prec V_{2}\prec\cdots\prec V_{L}$. In fact, a valid layering can be found for any DAG, though some layers may contain a single node. As such, the notion of layering is general: We make no assumption on the size of each layer, the causal ordering of variables in each set $V_{\ell}$, or interactions among them, except that $G$ is a DAG.

To simplify the presentation, we primarily focus on two-layer, or bipartite, DAGs, and defer the discussion of more general cases to Section 5. Let $V_{1}={{\mathcal{X}}}\equiv\{X_{1},\ldots,X_{p}\}$ be the nodes in the first layer and $V_{2}={{\mathcal{Y}}}\equiv\{Y_{1},\ldots,Y_{q}\}$ be those in the second layer. In the case of causal inference for multiple outcomes discussed in Section 1, ${{\mathcal{X}}}$ represents the exposure variables, and ${{\mathcal{Y}}}$ represents the outcome variables.

Throughout the paper, we denote individual variables by regular upper case letters, e.g., $X_{k}$ and $Y_{j}$, and sets of variables as calligraphic upper case letters, e.g., ${{\mathcal{X}}}_{-k}$ and ${{\mathcal{X}}}_{S}$, where $S$ is a set containing more than one index. We also refer to nodes of the network by using both their indices and the corresponding random variables; for instance, $k\in{{\mathcal{X}}}$ and $X_{k}$. We denote the parents of a node $j$ in $G$ by $\textnormal{pa}_{j}$, its ancestors and children by $\mathrm{an}_{j}$ and $\textnormal{ch}_{j}$, respectively, and its adjacent nodes by $\textnormal{adj}_{j}$, where $\textnormal{adj}_{j}=\textnormal{pa}_{j}\cup\textnormal{ch}_{j}$. We may also represent the edges of $G$ using its adjacency matrix ${\Theta}$, which satisfies ${\Theta}_{jj}=0$ and ${\Theta}_{jj^{\prime}}\neq 0$ if and only if $j^{\prime}\in\textnormal{pa}_{j}$. For any bipartite DAG with layers ${{\mathcal{X}}}$ and ${{\mathcal{Y}}}$ (${{\mathcal{X}}}\prec{{\mathcal{Y}}}$), we can partition ${\Theta}$ into the following block matrix

where the zero constraint on the upper right block of ${\Theta}$ follows from the partial ordering. Here $A$ and $C$ contain the information on the edges amongst $X_{k}$ $(k=1,\ldots,p)$ and $Y_{j}$ $(j=1,\ldots,q)$, respectively. Both of these matrices can be written as lower-triangular matrices (Shojaie and Michailidis, 2010). However, there are generally no constraints on $B$. We denote by $H$ the subgraph of $G$ containing edges from ${{\mathcal{X}}}$ to ${{\mathcal{Y}}}$ only, i.e., those corresponding to entries in $B$. Our goal is to estimate $H$ using the fact that ${{\mathcal{X}}}\prec{{\mathcal{Y}}}$.

Figure 2a shows a simple example of a two-layer DAG $G$ with layers consisting of $p=q=2$ nodes. This example illustrates the difference between full causal orderings of variables in a DAG and partial causal orderings: In this case, the graph admits two full causal orderings, namely $\mathcal{O}_{1}=\{X_{1},X_{2},Y_{1},Y_{2}\}$ and $\mathcal{O}_{2}=\{X_{1},Y_{1},X_{2},Y_{2}\}$; either of these orderings can be used to correctly discover the structure of the graph. Here, the partial ordering of the variables defined by the layering of the graph to sets ${{\mathcal{X}}}$ and ${{\mathcal{Y}}}$, i.e., ${{\mathcal{X}}}\prec{{\mathcal{Y}}}$, can be written as $\mathcal{O}^{\prime}=\{\{X_{1},X_{2}\},\{Y_{1},Y_{2}\}\}$. This ordering determines that $X_{k}$s should appear before $Y_{j}$s in the causal ordering but does not restrict the ordering of $\{X_{1},X_{2}\}$ and $\{Y_{1},Y_{2}\}$. In this example, only $\mathcal{O}_{1}$ is consistent with the partial ordering $\mathcal{O}^{\prime}$. Moreover, $\mathcal{O}^{\prime}$ is also consistent with $\{X_{2},X_{1},Y_{2},Y_{1}\}$, which is not a valid causal ordering of $G$, indicating that a partial causal order is not sufficient for learning the full DAG.

In this section, we discuss the challenges of estimating the graph $H$ and why this problem cannot be solved using simple approaches. Since the partial ordering of nodes, $V_{1}\prec V_{2}$, provides information about direction of causality between the two layers of the network, we focus on simple constraint-based methods (Spirtes et al., 2001), which learn the network edges based on conditional independence relationships among nodes. This requires conditional independence relations among variables to be compatible with the edges in $G$; formally, the joint probability distribution $\mathcal{P}$ needs to be faithful to $G$ (Spirtes et al., 2001). (As discussed in Section 4, our algorithm requires a weaker notion of faithfulness; however, for simplicity, we consider the classical notion of faithfulness in this section.)

Given the partial ordering of variables—which means that $Y_{j}$s cannot be parents of $X_{k}$s—one approach for estimating the bipartite graph $H$ is to draw an edge from $X_{k}$ to $Y_{j}$ whenever $Y_{j}$ is dependent on $X_{k}$ given all other nodes in the first layer, ${{\mathcal{X}}}_{-k}\equiv\{X_{k^{\prime}},k^{\prime}\neq k\}$. Formally, denoting by $Y_{j}{\bot\negthickspace\negthickspace\bot}X_{k}$ the conditional independence of two variables $Y_{j}$ and $X_{k}$, we define

to emphasize that the estimate is obtained without conditioning on any $Y_{j^{\prime}}\neq Y_{j}$.

Figure 2b shows the estimated $H^{(0)}$ in the setting where true causal relationships are linear. In this setting, edges in $H^{(0)}$ represent nonzero coefficients in linear regressions of each $Y_{j}$ on its parents in $G$, $\textnormal{pa}_{j}$. It can be seen that, in this example, the true causal effects from $X_{k}$s to $Y_{j}$s are in fact captured in $H^{(0)}$; these are shown in solid blue lines in Figure 2b. However, this simple example suggests that $H^{(0)}$ may include spurious edges, shown by orange dashed edges in the figure. The next lemma formalizes and generalizes this finding. The proof of this and other results in the paper are gathered in the Appendix A.

Lemma 2.1 suggests that if $G$ includes edges among $Y_{j}$s, failing to condition on $Y_{j}$s can result in falsely detected causal effects from $X_{k}$s to $Y_{j}$s. Thus, without knowledge of interactions among $Y_{j}$s, one may consider an estimator that corrects for both $X_{-k}$ and $Y_{-j}$ when trying to detect the casual effects of $X_{k}$ on $Y_{j}$; in other words, we may declare an edge from $X_{k}$ to $Y_{j}$ whenever $Y_{j}{\medspace/\negmedspace\negmedspace\negthickspace\bot\negthickspace\negthickspace\bot}X_{k}\mid\{{{\mathcal{X}}}_{-k},{{\mathcal{Y}}}_{-j}\}$. We denote the resulting estimate of $H$ as $H^{(-j)}$:

Unfortunately, as Figure 2c shows, estimation based on this model may also include false positive edges. As the next lemma clarifies, the false positive edges in this case are due to the conditioning on common descendants of a pair of $X_{k}$ and $Y_{j}$ that are not connected in $G$, which is sometimes referred to as Berkson’s Paradox (Pearl, 2009).

While false positive edges in $H^{(0)}$ are caused by failing to condition on necessary variables, the false positive edges in $H^{(-j)}$ are caused by conditioning on extra variables that are not part of the correct causal order of variables. More generally, the partial ordering information does not lead to a simple estimator that correctly identifies direct causal effects of exposures, $X_{k}\,(k=1,\ldots,p)$, on outcomes, $Y_{j}\,(j=1,\ldots,q)$.

Building on the above findings, in the next section we first present a modification of the PC algorithm that leverages the partial ordering information. We then present a general framework for leveraging the partial ordering of variables more effectively, especially when the goal is to learn the graph $H$.

Given an ideal test of conditional independence, we can estimate $H$ by first learning the skeleton of $G$, using the PC Algorithm (Kalisch and Bühlmann, 2007) and then orienting the edges in $H$ according to the partial ordering of nodes. However, such an algorithm uses the ordering information in a post hoc way—the information is not utilized to learn the edges in $H$.

Alternatively, we can incorporate the partial ordering directly into the PC algorithm. Recall that PC starts from a fully-connected graph and iteratively removes edges by searching for conditional independence relations that corresponds to graphical d-separations, to this end, it uses separating sets of increasing sizes, starting with sets of size 0 (i.e., no conditioning) and incrementally increasing the size as the algorithm progresses. To modify the PC algorithm, we use the fact that if two nodes $j$ and $k$ are non-adjacent in a DAG, then for any of their d-separators, say $S$, the set $S\cap\mathrm{an}(\{j,k\})$ is also a d-separator; see Chapter 6 of Spirtes et al. (2001). Given a partial ordering of variables, this rule allows excluding from the PC search step any nodes that are in lower layers than nodes $j$ and $k$ under consideration. Then, under the assumptions of the PC algorithm, the above modification, which we refer to as PC+, enjoys the same population and sample guarantees.

Given its construction, the partial ordering becomes more informative in PC+ as the number of layers increases: in two-layer settings, for learning the edges in $H$, PC+ is exactly the same as PC, and cannot utilize the partial ordering. However, as the results in Section 6.2 show, PC+ remains ineffective in graphs with larger number of layers, motivating the development of the framework proposed in the next section.

In this section, we propose a new framework for learning DAGs from partial orderings. The proposed approach is motivated by two key observations in Lemmas 2.1 and 2.2: First, the graphs $H^{(0)}$ and $H^{(-j)}$ include all true causal relationships from nodes $X_{k}$ $(k=1,\ldots,p)$ in the first layer to nodes $Y_{j}$ $(j=1,\ldots,q)$ in the second layer. Second, both graphs may also include additional edges; $H^{0}$ due to not conditioning on parents of $Y_{j}$ in ${{\mathcal{Y}}}$, and $H^{-j}$ due to conditioning on common children of $X_{k}$ and $Y_{j}$.

Let $S^{(0)}_{j}:=\{k:(k\to j)\in H^{(0)}\}$ and $S^{(-j)}_{j}:=\{k:(k\to j)\in H^{(-j)}\}$. The next lemma, which is a direct consequence of Lemmas 2.1 and 2.2, characterizes the intersection of these two sets.

Lemma 3.1 implies that even though $H^{(0)}_{j}\cap H^{(-j)}_{j}$ may contain more edges than the true edges in $H$, these additional edges must be in some special configuration. Based on this insight, our approach examines the edges in $S^{(0)}_{j}\cap S^{(-j)}_{j}$, for each $j\in{{\mathcal{Y}}}$, to remove the spurious edges efficiently. Let $k\in S^{(0)}_{j}\cap S^{(-j)}_{j}$ and $(k\to j)\notin H$. Suppose, for simplicity, that conditional independencies are faithful to the graph $G$. Then, there must be some set of variables $\mathcal{Z}\subset{{\mathcal{X}}}\cup{{\mathcal{Y}}}$ such that $X_{l}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y_{j}|\mathcal{Z}$. In general, to find such a set of variables, we will need to search among subsets of ${{\mathcal{X}}}\cup{{\mathcal{Y}}}$. However, utilizing the partial ordering, we can significantly reduce the complexity of this search. Specifically, we can restrict the search to the conditional Markov blanket, introduced next for general DAGs.

Our new framework builds on the idea that by limiting the search to the conditional Markov blankets, given the nodes in the previous layer, we can significantly reduce the search space and the size of conditioning sets. Moreover, the conditional Markov blanket can be easily inferred together with $H^{(0)}$ and $H^{(-j)}$. Coupled with the observations in Lemmas 2.1 and 2.2, especially the fact that conditioning on the nodes in the previous layer does not remove true causal edges, these reductions lead to improvements in both computational and statistical efficiency.

The new framework is summarized in Algorithm 1. The algorithm has two main steps: a screening loop, where supersets of relevant edges are identified by leveraging the partial ordering information, and a searching loop, similar to the one in the PC algorithm, but tailored to searching over the conditional Markov blankets. To this end, we will next show that in order to learn the edges between any node in ${{\mathcal{X}}}$ and a node in ${{\mathcal{Y}}}$, it suffices to search over subsets of the conditional Markov blanket of nodes in ${{\mathcal{Y}}}$.

The next lemma characterizes key properties of the conditional Markov blanket, which facilitate efficient learning of the graph $H$. Specifically, we show that if a distribution satisfies the intersection property of conditional independence—i.e., if $X\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y\mid Z\cup S$ and $X\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Z\mid Y\cup S$ then $X\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y\cup Z\mid S$—then all Markov blankets and conditional Markov blankets are unique.

We next define the faithfulness assumption needed for correct causal discovery from partial orderings. This assumption is trivially weaker than the general notion of strong faithfulness (Zhang and Spirtes, 2002).

Our main result, given below, describes how conditional Markov blankets can be used to effectively incorporate the knowledge of partial ordering into DAG learning and reduce the computational cost of learning causal effects of ${{\mathcal{X}}}_{k}$s on ${{\mathcal{Y}}}_{j}$s.

It follows directly from Theorem 3.5 that the proposed framework in Algorithm 1—i.e., taking the intersection of $H^{(0)}$ and $H^{(-j)}$ and then searching within conditional Markov blankets to remove additional edges—recovers the correct bipartite DAG $H$.

Next, we show that instead of inferring $H^{(0)}$ and $H^{(-j)}$ separately and taking their interception, we can use $H^{(0)}$ to infer $H^{(0)}\cap H^{(-j)}$ directly. More specifically, noting that Algorithm 1 only relies on $S^{(0)}_{j}\cap S^{(-j)}_{j}$ and ${\mathrm{cmb}}_{{{\mathcal{X}}}}(j)$ for each $j\in{{\mathcal{Y}}}$, the next lemma shows that given $S^{(0)}_{j}$, we can learn $S^{(0)}_{j}\cap S^{(-j)}_{j}$ and ${\mathrm{cmb}}_{{{\mathcal{X}}}}(j)$ without having to learn $S^{(-j)}_{j}$ separately. This implies that Algorithm 1 need not learn the unconditional Markov blankets, but only the conditional Markov blankets.

Lemma 3.7 characterises the conditional dependence sets used in Algorithm 1. Using these characterizations, we can cast the set inference problems in Algorithm 1 as variable selection problems in general regression settings. Let

Then, $S^{(0)}_{j}$ and $S^{(1)}_{j}$ can be obtained by selecting the relevant variables when regressing $Y_{j}$ onto ${{\mathcal{X}}}$ and ${{\mathcal{X}}}_{S^{(0)}_{j}}\cup{{\mathcal{Y}}}_{-j}$, respectively. We can then obtain the two sets used in Algorithm 1 from $S^{(1)}_{j}$:

There are multiple benefits to directly inferring $H^{(0)}\cap H^{(-j)}$—i.e., by estimating $S^{(0)}_{j}$ and $S^{(1)}_{j}$—instead of separately inferring $H^{(0)}$ and $H^{(-j)}$. First, the graph $H^{(-j)}$ could be hard to estimate in high-dimensional settings, as learning $H^{(-j)}$ requires performing tests conditioned on $p+q-2$ variables. In contrast, $S_{j}^{(1)}$ in (5) only requires performing tests conditioning on at most $\max\left\{p-1,\max_{j}|S^{(0)}_{j}|+q-2\right\}$ variables. Moreover, the target $H^{(0)}\cap H^{(-j)}$ is often sparse (see Lemma 3.1) and can thus be efficiently learned in high-dimensional settings. In an extreme example, suppose each node in ${{\mathcal{X}}}$ has exactly one outgoing edge into a node in ${{\mathcal{Y}}}$, and the first $q-1$ nodes in ${{\mathcal{Y}}}$ have as their common child $Y_{q}$. In this case, $H^{(-j)}$ is fully connected and hard to learn, whereas $H^{(0)}\cap H^{(-j)}$ only has $p$ edges.

The second advantage of directly inferring $H^{(0)}\cap H^{(-j)}$ is that, by using the conditional dependence formulation of sets as in Lemma 3.7, we can show that even if the sets $S^{(0)}$ and $S^{(1)}_{j}$—consequently, $S^{(0)}_{j}\cap S^{(-j)}_{j}$ and ${\mathrm{cmb}}_{{{\mathcal{X}}}}(j)$—are not inferred exactly, the algorithm is still correct as long as no false negative errors are made.

We refer to the framework proposed in Algorithm 1 for learning DAGs from partial ordering the PODAG framework. In the next section, we discuss specific algorithms for learning direct casual effects of $X_{k}$s on $Y_{j}$s. These algorithms utilize the fact that given the partial ordering of nodes in $G$, the conditional Markov blanket of each node $j\in{{\mathcal{Y}}}$ can be found efficiently by testing for conditional dependence of $Y_{j}$ and $Y_{j^{\prime}},j^{\prime}\neq j$ after adjusting for the effect of the nodes in the first layer.

Given multivariate Gaussian observations, the most direct way to identify the sets satisfying the conditional independence relations in $S^{(0)}_{j}$ and $S^{(1)}_{j}$ is to use the sample partial correlation $\hat{\rho}$ and reject the hypothesis of conditional independence when $|\hat{\rho}|>\xi$ for some suitable threshold $\xi$. The next results establishes the screening guarantee for such an approach in low-dimensional settings.

The conditioning sets in $S^{(0)}_{j}$ and $S^{(1)}_{j}$—namely, ${{\mathcal{X}}}_{-k}$ and ${{\mathcal{X}}}_{S^{(0)}_{j}}\cup{{\mathcal{Y}}}_{-j}\setminus Z$—involve large number of variables. Therefore, in high dimensions, when $p,q\gg n$, it is not feasible to learn $S^{(0)}_{j}$ and $S^{(1)}_{j}$ using partial correlations. To overcome this challenge, and using the insight from Lemma 3.7, we treat the screening problem as selection of non-zero coefficients in linear regressions. This allows us to use tools from high-dimensional estimation, including various regularization approaches, for the screening step.

The success of the screening step requires the event $\mathcal{A}\left(\widehat{S}^{(0)}_{j},\widehat{S}^{(1)}_{j}\right)$—defined in Theorem 4.4—to hold with high probability. That is, we need to identify supersets of $S^{(0)}_{j}$ and $S^{(1)}_{j}$ such that $\left|\widehat{S}^{(0)}_{j}\cap\widehat{S}^{(1)}_{j}\right|\leq n$. A simple approach to screen relevant variables is sure independence screening (SIS, Fan and Lv, 2008), which selects the variables with largest marginal association at a given threshold $t\in(0,1)$. Defining the SIS estimators

we can obtain the following result following Fan and Lv (2008).

By using simple marginal associations, SIS offers a simple and efficient strategy for screening; however, it may result in large sets for the searching step. An alternative to SIS screening is to select the sets $\widehat{S}_{j}^{(0)}$ and $\widehat{S}_{j}^{(1)}$ using two penalized regressions. For instance, using the lasso penalty, $\widehat{S}_{j}^{(0)}$ for $j\in{{\mathcal{Y}}}$ can be found as

Similarly, $S_{j}^{(1)}$ can be learned using a lasso regression on a different set of variables; specifically, given any set $S\subseteq{{\mathcal{X}}}$, we define

The next result establishes the success of the screening step when using lasso estimators in high dimensions.

The above two high-dimensional screening approaches—SIS and lasso—are two extremes: SIS obtains the estimates of $S^{(0)}_{j}$ and $S^{(1)}_{j}$ without adjusting for any other variables, whereas lasso adjusts for all other variables (either from the first layer or the two layers combined). These two extremes represent examples of broader classes of methods based on marginal or joint associations. They differ with respect to computational and sample complexities of screening and searching steps: Compared with SIS, screening with lasso may result in smaller sets for the searching step but may be less efficient for screening. Nevertheless, both approaches are valid choices for screening under milder conditions than those needed for consistent variable selection. This is because the additional assumption required for the screening step is implied by the layering-adjacency-faithfulness assumption, which, as discussed before, is weaker than the strong faithfulness assumption needed for the PC algorithm (see also Section 6.2). While the above two choices seem natural, other intermediates, including conditioning on small sets (similar to Sondhi and Shojaie, 2019) can also be used for screening.

Together with Proposition 4.5, 4.6, or 4.7, Theorem 4.4 establishes the consistency of Algorithm 1 for linear SEMs, in low- and high-dimensional settings. While we only considered the simple case of Gaussian noise, Algorithm 1 and the results in this section can be extended to linear SEMs with sub-Gaussian errors (Harris and Drton, 2013). Moreover, many other regularization methods can be used for screening instead of the lasso. The only requirement is that the method is screening-consistent. For instance, inference-based procedures, such as debiased lasso (van de Geer et al., 2014; Zhang and Zhang, 2014; Javanmard and Montanari, 2014) can also fulfill the requirement of our screening step.

Similar ideas can also be used for learning DAGs from other distributions and/or under more general SEM structures. In Appendix B, we outline the generalization of the proposed algorithms to causal additive models (Bühlmann et al., 2014). The Gaussian copula transformation of Liu et al. (2009, 2012) and Xue and Zou (2012) can be similarly used to extend the proposed ideas to non-Gaussian distribution. We leave the detailed exploration of these generalizations to future research.

In the previous sections, we focused on the problem of learning edges between two layers, ${{\mathcal{X}}}$ and ${{\mathcal{Y}}}$. In particular, Algorithm 1 provides a framework that first finds a superset of the true edges, and then uses a searching loop to remove the spurious edges. The same idea can be applied to learning edges within the second layer, ${{\mathcal{Y}}}$: As suggested by Lemma 3.3, for each $j\in{{\mathcal{Y}}}$, all of its adjacent nodes are contained in $S^{(1)}_{j}$. This suggests that we can learn edges in ${{\mathcal{Y}}}$ by simply modifying Algorithm 1 to run the searching loop on $\left\{(k,j):k\in S^{(1)}_{j},j\in{{\mathcal{Y}}}\right\}$ instead of only on those between ${{\mathcal{X}}}$ and ${{\mathcal{Y}}}$. After recovering the skeleton, edges among ${{\mathcal{Y}}}$ can be oriented using Meek’s orientation rules (Meek, 1995). Given correct d-separators and skeleton, the rules can orient as many edges as possible without making mistakes.

In order to successfully recover within layer edges from observational data, we need to assume faithfulness among these layers. Following Ramsey et al. (2006), we characterize the faithfulness condition required to learn the DAG via constraint-based algorithms as two parts: (a) adjacency faithfulness, which means that neighboring nodes are not associated with conditional independence; and (b) orientation faithfulness, which means that parents in any v-structure are not conditionally independent given any set containing their common child. The orientation faithfulness can be defined in our context as follows.

Given the additional information from partial ordering, to describe the graphical object learned by the new framework, we need to introduce a new notion of equivalence. This is because the Markov equivalent class of DAGs can be reduced by using the known partial ordering. For example, if the only conditional independent relation among three variables $X_{i},X_{j},X_{k}$ is $X_{i}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptstyle\perp$\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X_{k}|X_{j}$, but we know $X_{i}\prec\{X_{j},X_{k}\}$, then the edge $X_{j}\to X_{k}$ is identifiable. We define partial-ordering-Markov-equivalence simply as Markov equivalence restricted to partial ordering. This equivalence class can be represented by a maximally oriented partial DAG, or maximal PDAG (Perkovic et al., 2018). With this notion, we can describe the target of learning of within- and between-layer edges as learning the maximal PDAG of the true $G$ given the background information ${{\mathcal{X}}}\prec{{\mathcal{Y}}}$.

The theory and algorithm developed in Section 4 can be extended to scenarios with multiple layers. To facilitate this extension, we introduce a general representation of the problem. Suppose $V$ is a random vector following some distribution Markov to a graph $G=(V,E)$. Suppose $G$ admits a partial-ordering $\mathcal{O}=\{V_{1}\prec\cdots\prec V_{L}\}$ where $V=\cup_{\ell=1}^{L}V_{\ell}$. Parallel to the notation in 2-layer case, we define, for each $j\in V_{\ell}$, $1\leq\ell\leq L$,

The next lemma generalizes Lemma 3.1.