# Inference for a Large Directed Acyclic Graph with Unspecified Interventions

Chunlin Li<EMAIL_ADDRESS>
School of Statistics,
University of Minnesota,
Minneapolis, MN 55455, USA Xiaotong Shen<EMAIL_ADDRESS>
School of Statistics,
University of Minnesota,
Minneapolis, MN 55455, USA Wei Pan<EMAIL_ADDRESS>
Division of Biostatistics,
University of Minnesota,
Minneapolis, MN 55455, USA

###### Abstract

Inference of directed relations given some unspecified interventions, that is,
the target of each intervention is unknown, is challenging. In this article,
we test hypothesized directed relations with unspecified interventions. First,
we derive conditions to yield an identifiable model. Unlike classical
inference, testing directed relations requires identifying ancestors and
relevant interventions of hypothesis-specific primary variables. Towards this
end, we propose a peeling algorithm based on nodewise regressions to establish
a topological order of primary variables. Moreover, we prove that the peeling
algorithm yields a consistent estimator in low order polynomial time. Second,
we propose a likelihood ratio test integrated with a data perturbation scheme
to account for the uncertainty of identifying ancestors and interventions.
Also, we show that the distribution of a data perturbation test statistic
converges to the target distribution. Numerical examples demonstrate the
utility and effectiveness of the proposed methods, including an application to
infer gene regulatory networks. The implementation of the proposed methods is
available at https://github.com/chunlinli/intdag.

Keywords: High-dimensional inference, data perturbation, structure learning,
identifiability, peeling algorithm

##  1 Introduction

Directed relations are essential to explaining pairwise dependencies among
multiple interacting units. In gene network analysis, regulatory gene-to-gene
relations are a focus of biological investigation (Sachs et al., 2005), while
in a human brain network, scientists investigate causal influences among
regions of interest to understand how the brain functions through the regional
effects of neurons (Liu et al., 2017). In such a situation, inferring directed
effects without other information is generally impossible because of the lack
of identifiability in a Gaussian directed acyclic graph (DAG) model, and hence
external interventions are introduced to treat a non-identifiable situation
(Heinze-Deml et al., 2018). For instance, the genetic variants such as single-
nucleotide polymorphisms (SNPs) can be, and indeed are increasingly, treated
as external interventions to infer inter-trait causal relations in a
quantitative trait network (Brown and Knowles, 2020) and gene interactions in
a gene regulatory network (Teumer, 2018; Molstad et al., 2021). In
neuroimaging analysis, scientists use randomized experimental stimuli as
interventions to identify causal relations in a functional brain network
(Grosse-Wentrup et al., 2016; Bergmann and Hartwigsen, 2021). However, the
interventions in these studies often have unknown targets and off-target
effects (Jackson et al., 2003; Eaton and Murphy, 2007). Consequently,
inferring directed relations while identifying useful interventions for
inference is critical. This paper focuses on the simultaneous inference of
directed relations subject to unspecified interventions, that is, without
known targets.

In a DAG model, the research has been centered on the reconstruction of
directed relations in observational and interventional studies (van de Geer
and Bühlmann, 2013; Oates et al., 2016; Heinze-Deml et al., 2018; Zheng et
al., 2018; Yuan et al., 2019). For inference, Bayesian methods (Friedman and
Koller, 2003; Luo and Zhao, 2011; Viinikka et al., 2020) have been popular.
Yet, statistical inference remains under-studied, especially for
interventional models in high dimensions (Peters et al., 2016; Rothenhäusler
et al., 2019). Recently, for observational data, Janková and van de Geer
(2018) propose a debiased test of the strength of a single directed relation
and Li et al. (2019) derives a constrained likelihood ratio test for multiple
directed relations.

Despite progress, challenges remain. First, inferring directed relations
requires identifying a certain DAG topological order (van de Geer and
Bühlmann, 2013), while the identifiability in a Gaussian DAG with unspecified
interventions remains under-explored. Second, the inferential results should
agree with the acyclicity requirement for a DAG. As a result, degenerate and
intractable situations can occur, making inference greatly different from the
classical ones. Third, likelihood-based methods for learning the DAG
topological order often use permutation search (van de Geer and Bühlmann,
2013) or continuous optimization subject to the acyclicity constraint (Zheng
et al., 2018; Yuan et al., 2019), where a theoretical guarantee of the actual
estimate (instead of the global optimum) has not been established for these
approaches. Recently, an important line of work (Ghoshal and Honorio, 2018;
Rajendran et al., 2021; Rolland et al., 2022) has focused on order-based
algorithms with computational and statistical guarantees. However, in Gaussian
DAGs, existing methods often rely on some error scale assumptions (Peters and
Bühlmann, 2014), which is sensitive to variable scaling like the common
practice of standardizing variables. This drawback could limit their
applications, especially in causal inference, as causal relations are
typically invariant to scaling.

To address the above issues, we develop structure learning and inference
methods for a Gaussian DAG with unspecified additive interventions. Unlike the
existing approach where structure discovery and subsequent inference are
treated separately, our proposal integrates DAG structure learning and testing
directed relations (also known as the natural direct effect in causal
inference), accounting for the uncertainty of structure learning for
inference. With suitable interventions called instrumental variables (IVs),
the proposed approach removes the restrictive error scale assumptions and
delivers creditable outcomes with theoretical guarantees in low order
polynomial time. This result indicates IVs, a well-known tool in causal
inference, can play important roles in structure learning even if some
interventions do not meet the IV criterion. More specifically, our
contributions are summarized as the following aspects.

  * •

For modeling, we establish the identifiability conditions for a Gaussian DAG
with unspecified interventions. In particular, the conditions allow
interventions on more than one target, which is suitable for multivariate
causal analysis (Murray, 2006). Moreover, we introduce the concepts of
nondegeneracy and regularity for hypothesis testing under a DAG model to
characterize the behavior of a test.

  * •

For methodology, we develop likelihood ratio tests for directed edges and
pathways in a super-graph of the true DAG, where the super-graph is formed by
ancestral relations and candidate interventional relations, offering the
topological order for inference. Furthermore, we reconstruct the super-graph
by the peeling algorithm, which automatically meets the acyclicity
requirement. By integrating structure learning with inference, we account for
the uncertainty of the super-graph estimation for the proposed tests via a
novel data perturbation (DP) scheme, which effectively controls the type-I
error while enjoying high statistical power.

  * •

For theory, we prove that the proposed peeling algorithm based on nodewise
regressions yields consistent results in
$O(p\times\log\kappa_{\max}^{\circ}\times(q^{3}+nq^{2}))$ operations almost
surely, where $p,q$ are the numbers of primary and intervention variables, $n$
is the sample size, and $\kappa_{\max}^{\circ}$ is the sparsity. On this
basis, we justify the proposed DP inference method by establishing the
convergence of the DP likelihood ratio to the null distribution and desired
power properties.

  * •

The numerical studies and real data analysis demonstrate the utility and
effectiveness of the proposed methods. The implementation of the proposed
tests and structure learning method is available at
https://github.com/chunlinli/intdag.

The rest of the article is structured as follows. Section 2 establishes model
identifiability and introduces the concepts of nondegeneracy and regularity.
Section 3 develops the proposed methods for structure learning and statistical
inference. Section 4 presents statistical theory to justify the proposed
methods. Section 5 performs simulation studies, followed by an application to
infer gene pathways with gene expression and SNP data in Section 6. Section 7
concludes the article. The Appendix contains illustrative examples and
technical proofs.

##  2 Gaussian directed acyclic graph with additive interventions

To infer directed relations among primary variables (variables of primary
interest) $\bm{Y}=(Y_{1},\ldots,Y_{p})^{\top}$, consider a structural equation
model with additive interventions:

$\displaystyle\bm{Y}=\bm{U}^{\top}\bm{Y}+\bm{W}^{\top}\bm{X}+\bm{\varepsilon},\quad\bm{\varepsilon}\sim
N(\bm{0},\bm{\Sigma}),\quad\bm{\Sigma}=\text{Diag}(\sigma^{2}_{1},\ldots,\sigma^{2}_{p}),$
(1)

where $\bm{X}=(X_{1},\ldots,X_{q})^{\top}$ is a vector of additive
intervention variables, $\bm{U}\in\mathbb{R}^{p\times p}$ and
$\bm{W}\in\mathbb{R}^{q\times p}$ are unknown coefficient matrices, and
$\bm{\varepsilon}=(\varepsilon_{1},\ldots,\varepsilon_{p})^{\top}$ is a vector
of random errors with $\sigma_{j}^{2}>0$; $j=1,\ldots,p$. In (1),
$\bm{\varepsilon}$ is independent of $\bm{X}$ but components of $\bm{X}$ can
be dependent. The matrix $\bm{U}$ specifies the directed relations among
$\bm{Y}$, where $U_{kj}\neq 0$ if $Y_{k}$ is a direct cause of $Y_{j}$,
denoted by $Y_{k}\rightarrow Y_{j}$, and $Y_{k}$ is called a parent of $Y_{j}$
or $Y_{j}$ a child of $Y_{k}$. Thus, $\bm{U}$ represents a directed graph,
which is required to be acyclic to ensure the validity of the local Markov
property (Spirtes et al., 2000). The matrix $\bm{W}$ specifies the targets and
strengths of interventions, where $W_{lj}\neq 0$ indicates $X_{l}$ intervenes
on $Y_{j}$, denoted by $X_{l}\to Y_{j}$. In (1), no directed edge from a
primary variable $Y_{j}$ to an intervention variable $X_{l}$ is permissible.

In what follows, we will focus on the DAG $G$ with nodes
$\bm{Y}=(Y_{1},\ldots,Y_{p})$, primary variable edges
$\mathcal{E}=\\{(k,j):U_{kj}\neq 0\\}$, and intervention edges
$\mathcal{I}=\\{(l,j):W_{lj}\neq 0\\}$.

###  2.1 Identifiability

Model (1) is generally non-identifiable without interventions
($\bm{W}=\bm{0}$) when errors do not meet some requirements such as the equal-
variance assumption (Peters and Bühlmann, 2014) and its variants (Ghoshal and
Honorio, 2018; Rajendran et al., 2021). Moreover, a model can be identified
when $\bm{\varepsilon}$ in (1) is replaced by non-Gaussian errors (Shimizu et
al., 2006) or linear relations are replaced by nonlinear ones (Hoyer et al.,
2009). Regardless, suitable interventions can make (1) identifiable. When
intervention targets are known, the identifiability issue has been studied
(Oates et al., 2016; Chen et al., 2018). However, it is less so when the exact
targets and strengths of interventions are unknown, as in many biological
applications (Jackson et al., 2003; Kulkarni et al., 2006), referring to the
case of _unspecified_ or _uncertain_ interventions (Heinze-Deml et al., 2018;
Eaton and Murphy, 2007; Squires et al., 2020).

We now categorize interventions as instruments and invalid instruments.

######  Definition 1 (Instrument)

An intervention is an instrument in a DAG if it satisfies:

  1. (A)

(Relevance) it intervenes on at least one primary variable;

  2. (B)

(Exclusion) it does not intervene on more than one primary variable.

Otherwise, it is an invalid instrument.

Here, (A) requires an intervention to be active, while (B) prevents
simultaneous interventions of a single intervention variable on multiple
primary variables. This is critical to identifiability because in this
situation an instrument on a cause variable $Y_{1}$ helps reveal its directed
effect on an outcome variable $Y_{2}$, which would otherwise be impossible.
Thus, the instruments are introduced to break the symmetry that results in
non-identifiability of directed relations $Y_{1}\rightarrow Y_{2}$ and
$Y_{2}\rightarrow Y_{1}$.

Next, we make some assumptions on intervention variables to yield an
identifiable model, where dependencies among intervention variables are
permissible.

###### Assumption 1

Assume that model (1) satisfies the following conditions.

  1. (1A)

(Positive definiteness) $\operatorname{\mathbb{E}}\bm{X}\bm{X}^{\top}$ is
positive definite, where $\operatorname{\mathbb{E}}$ is the expectation.

  2. (1B)

(Local faithfulness)
$\operatorname{Cov}(Y_{j},X_{l}\mid\bm{X}_{\\{1,\ldots,q\\}\setminus\\{l\\}})\neq
0$ when $X_{l}$ intervenes on $Y_{j}$, where $\operatorname{Cov}$ denotes the
covariance.

  3. (1C)

(Instrumental sufficiency) Each primary variable is intervened by at least one
instrument.

Assumption 1A imposes mild distributional restrictions on $\bm{X}$, permitting
discrete variables such as SNPs. Assumption 1B requires interventional effects
not to cancel out as multiple targets from an invalid instrument are
permitted. Importantly, if either Assumption 1B or Assumption 1C fails, model
(1) is generally not identifiable, as shown in Example 2 of Appendix A.1. In
Section 5, we empirically examine the situation when Assumption 1C is not met.

######  Proposition 1 (Identifiability)

Under Assumption 1, the parameters $(\bm{U},\bm{W},\bm{\Sigma})$ in model (1)
are identifiable.

Proposition 1 (proved in Appendix B.1) is derived for a DAG model with
unspecified interventions. This is in contrast to Proposition 1 of Chen et al.
(2018), which proves the identifiability of the parameters in a directed graph
with target-known instruments on each primary variable. Moreover, it is worth
noting that the estimated graph in Chen et al. (2018) may be cyclic and lacks
the local Markov property for causal interpretation (Spirtes et al., 2000).

###  2.2 Likelihood-based inference for a DAG

Our primary goal is to perform statistical inference of directed edges and
pathways. Let $\mathcal{H}\subseteq\\{(k,j):k\neq j,1\leq k,j\leq p\\}$ be an
edge set among primary variables $\\{Y_{1},\ldots,Y_{p}\\}$, where
$(k,j)\in\mathcal{H}$ specifies a (hypothesized) directed edge $Y_{k}\to
Y_{j}$ in (1). We shall focus on two types of testing with null and
alternative hypotheses $H_{0}$ and $H_{a}$. For simultaneous testing of
directed edges,

$\displaystyle H_{0}:U_{kj}=0;\text{ for all }(k,j)\in\mathcal{H}\quad\mbox{
versus }\quad H_{a}:U_{kj}\neq 0\text{ for some }(k,j)\in\mathcal{H};$ (2)

for simultaneous testing of pathways,

$\displaystyle H_{0}:U_{kj}=0;\text{ for some }(k,j)\in\mathcal{H}\quad\mbox{
versus }\quad H_{a}:U_{kj}\neq 0\text{ for all }(k,j)\in\mathcal{H},$ (3)

where $(\bm{U}_{\mathcal{H}^{c}},\bm{W},\bm{\Sigma})$ are unspecified nuisance
parameters and c is the set complement.

Note that $H_{0}$ in (3) is a composite hypothesis that can be decomposed into
sub-hypotheses

$H_{0,\nu}:U_{k_{\nu},j_{\nu}}=0,\quad\mbox{ versus
}\quad{H}_{a,\nu}:U_{k_{\nu},j_{\nu}}\neq 0;\quad\nu=1,\ldots,|\mathcal{H}|,$

where
$\mathcal{H}=\\{(k_{1},j_{1}),\ldots,(k_{|\mathcal{H}|},j_{|\mathcal{H}|})\\}$
and each sub-hypothesis is a directed edge test.

Now consider likelihood inference for (2). Given an independent and
identically distributed sample $(\bm{Y}_{i},\bm{X}_{i})_{i=1}^{n}$ from (1),
the log-likelihood is

$L(\bm{\theta},\bm{\Sigma})=-\frac{1}{2}\sum_{i=1}^{n}\|\bm{\Sigma}^{-1/2}((\bm{I}-\bm{U}^{\top})\bm{Y}_{i}-\bm{W}^{\top}\bm{X}_{i})\|^{2}_{2}-n\log\sqrt{\det(\bm{\Sigma})},$
(4)

where $\bm{\theta}=(\bm{U},\bm{W})$,
$\bm{\Sigma}=\text{Diag}(\sigma_{1}^{2},\ldots,\sigma_{p}^{2})$, and $\bm{U}$
is subject to the acyclicity constraint (Zheng et al., 2018; Yuan et al.,
2019) that no directed cycle is permissible in a DAG. The likelihood ratio is
defined as
$\text{Lr}=L(\widehat{\bm{\theta}}^{(1)},\widehat{\bm{\Sigma}})-L(\widehat{\bm{\theta}}^{(0)},\widehat{\bm{\Sigma}})$,
where $\widehat{\bm{\theta}}^{(1)}$ is a maximum likelihood estimate (MLE) in
(4) subject to the acyclicity constraint, $\widehat{\bm{\theta}}^{(0)}$ is an
MLE subject to an additional requirement $\bm{U}_{\mathcal{H}}=\bm{0}$, and
$\widehat{\bm{\Sigma}}$ is an estimator of $\bm{\Sigma}$.

In many statistical models, a likelihood ratio has a nondegenerate and
tractable distribution, for instance, a chi-squared distribution with degrees
of freedom $|\mathcal{H}|$. However, in (1), because of the acyclicity
constraint, some hypothesized edges in $\mathcal{H}$ can be absent given the
edge subset $\mathcal{E}\cap\mathcal{H}^{c}$ defined by the nonzero elements
of nuisance parameter $\bm{U}_{\mathcal{H}^{c}}$, where $\mathcal{E}$ denotes
the edge set of the whole DAG. As a result, Lr may converge to a distribution
with degrees of freedom less than $|\mathcal{H}|$ and the distribution may be
even intractable, making inference for a DAG greatly different from the
classical ones, as illustrated by Example 1.

Figure 1: A true DAG structure of five primary variables $Y_{1},\ldots,Y_{5}$
and five intervention variables $X_{1},\ldots,X_{5}$, where directed edges are
represented by solid arrows while dependencies among $\bm{X}$ are not
displayed.

###### Example 1

Consider the likelihood ratio test under null $H_{0}$ and alternative $H_{a}$
for the DAG displayed in Figure 1. Let $\mathcal{E}$ be the edge set of the
DAG.

  * •

$H_{0}:U_{21}=0$ versus $H_{a}:U_{21}\neq 0$, where $\mathcal{H}=\\{(2,1)\\}$.
Here, $(2,1)$ forms a cycle together with the edges in
$\mathcal{E}\cap\mathcal{H}^{c}$ (namely the edges not considered by the
hypothesis), violating the acyclicity constraint. Note that the maximum
likelihood value $L(\widehat{\bm{\theta}},\widehat{\bm{\Sigma}})$ with
$U_{21}\neq 0$ tends to be smaller than that with $U_{21}=0$ when a random
sample is obtained under $H_{0}$, especially so when the asymptotics kicks in
as the sample size increases. Consequently, the likelihood ratio Lr becomes
zero, constituting a degenerate situation.

  * •

$H_{0}:U_{45}=U_{53}=0$ versus $H_{a}:U_{45}\neq 0$ or $U_{53}\neq 0$, where
$\mathcal{H}=\\{(4,5),(5,3)\\}$. In this case, $\\{(4,5),(5,3)\\}$ forms a
cycle with the edges in $\mathcal{E}\cap\mathcal{H}^{c}$, violating the
acyclicity constraint. Since the likelihood value
$L({\bm{\theta}},\bm{\Sigma})$ tends to be maximized under the true graph when
data is sampled under $H_{0}$, we have
$L(\widehat{\bm{\theta}}^{(1)},\widehat{\bm{\Sigma}})=\max(L(\widehat{U}_{45},\widehat{U}_{53}=0,\widehat{\bm{U}}_{\mathcal{H}^{c}},\widehat{\bm{W}},\widehat{\bm{\Sigma}}),L(\widehat{U}_{45}=0,\widehat{U}_{53},\widehat{\bm{U}}_{\mathcal{H}^{c}},\widehat{\bm{W}},\widehat{\bm{\Sigma}}))$.
As a result, the likelihood ratio distribution becomes complicated in this
situation due to the dependence between the two components in
$L(\widehat{\bm{\theta}}^{(1)},\widehat{\bm{\Sigma}})$.

Motivated by Example 1, we introduce the concepts of nondegeneracy and
regularity.

######  Definition 2 (Nondegeneracy and regularity with respect to DAG)

  1. (A)

An edge $(k,j)\in\mathcal{H}$ is nondegenerate with respect to DAG $G$ if
$\\{(k,j)\\}\cup\mathcal{E}$ contains no directed cycle, where $\mathcal{E}$
denotes the edge set of $G$. Otherwise, $(k,j)$ is degenerate. Let
$\mathcal{D}\subseteq\mathcal{H}$ be the set of all nondegenerate edges with
respect to $G$. A null hypothesis $H_{0}$ is nondegenerate with respect to DAG
$G$ if $\mathcal{D}\neq\emptyset$. Otherwise, $H_{0}$ is degenerate.

  2. (B)

A null hypothesis $H_{0}$ is said to be regular with respect to DAG $G$ if
$\mathcal{D}\cup\mathcal{E}$ contains no directed cycle, where $\mathcal{E}$
denotes the edge set of $G$. Otherwise, $H_{0}$ is called irregular.

###### Remark 1

In practice, $\mathcal{D}$ is unknown and needs to be estimated from data.

Nondegeneracy ensures nonnegativity of the likelihood ratio. In testing (2),
regularity excludes intractable situations for the null distribution. In
testing (3), if $H_{0}$ is irregular, then $\mathcal{D}\cup\mathcal{E}$ has a
directed cycle, which means the hypothesized directed pathway cannot exist due
to the acyclicity constraint. Thus, regularity excludes the degenerate
situations in testing (3).

In what follows, we mainly focus on nondegenerate and regular hypotheses. For
the degenerate case, the p-value is defined to be one. For the irregular case
of edge test (2), we decompose the hypothesis into regular sub-hypotheses and
conduct multiple testing. For the irregular case of pathway test (3), the
p-value is defined to be one. More discussions on the implementation in
irregular cases are provided in Appendix A.5.

Finally, we introduce some notations for a DAG $G$ with primary variables
$\\{Y_{1},\ldots,Y_{p}\\}$, intervention variables $\\{X_{1},\ldots,X_{q}\\}$,
a directed edge set $\mathcal{E}\subseteq\\{(k,j):k\neq j,1\leq k,j\leq p\\}$,
and an intervention edge set $\mathcal{I}\subseteq\\{(l,j):1\leq l\leq q,1\leq
j\leq p\\}$. If there is a directed path $Y_{k}\to\cdots\to Y_{j}$ in $G$,
$Y_{k}$ is called an ancestor of $Y_{j}$ or $Y_{j}$ a descendant of $Y_{k}$.
Let $\mathcal{A}=\\{(k,j):Y_{k}\to\cdots\to Y_{j}\\}$ be the set of all
ancestral relations. In a DAG $G$, let
$\textnormal{{pa}}_{G}(j)=\\{k:(k,j)\in\mathcal{E}\\}$,
$\textnormal{{an}}_{G}(j)=\\{k:(k,j)\in\mathcal{A}\\}$,
$\textnormal{{d}}_{G}(j)=\\{k:(k,j)\in\mathcal{D}\\}$, and
$\textnormal{{in}}_{G}(j)=\\{l:(l,j)\in\mathcal{I}\\}$ be the parent,
ancestor, nondegenerate hypothesized parent, and intervention sets of $Y_{j}$,
respectively.

##  3 Methodology

This section develops the main methodology, including the peeling algorithm
for structure learning and the data perturbation inference for simultaneous
testing of directed edges (2) and pathways (3).

One major challenge to the likelihood approach is optimization subject to the
acyclicity requirement. This constraint imposes difficulty on not only
computation but also asymptotic theory. As a result, there is a gap between
the asymptotic distribution based on a global maximum and that on the actual
estimate which could be a local maximum (Janková and van de Geer, 2018; Li et
al., 2019). Moreover, the actual estimate may give an imprecise topological
order, tending to impact adversely on inference.

To circumvent the acyclicity requirement, we propose using a super-graph $S$
of the true DAG $G$, where $S$ is constructed by a novel peeling method
without explicitly imposing the acyclicity constraint. As a result, our method
is scalable with a statistical guarantee of the actual estimate. We will focus
on a specific super-graph $S$ consisting of primary variable edges
$\mathcal{E}_{S}:=\mathcal{A}\supseteq\mathcal{E}$ and candidate intervention
edges $\mathcal{I}_{S}:=\mathcal{C}\supseteq\mathcal{I}$, where $\mathcal{E}$
and $\mathcal{I}$ are sets of primary variable edges and intervention edges in
the original graph $G$, and set $\mathcal{C}$ is to be defined shortly after
(6). Note that $S$, whose primary variable edges consist of all ancestral
relations of $G$, is also a DAG. Given the super-graph $S$, the true DAG $G$
can be reconstructed by nodewise constrained regressions.

Let $\bm{\theta}_{S}=\big{(}\bm{U}_{\mathcal{A}},\bm{W}_{\mathcal{C}}\big{)}$
be $\bm{\theta}$ restricted to $S$. Note that
$\bm{\theta}_{S^{c}}=(\bm{U}_{\mathcal{A}^{c}},\bm{W}_{\mathcal{C}^{c}})=\bm{0}$
in the original graph $G$. Given $S$, the log-likelihood
$L(\bm{\theta},\bm{\Sigma})$ in (4) can be rewritten as

$\begin{split}L_{S}(\bm{\theta}_{S},\bm{\Sigma})&=-\sum_{j=1}^{p}\Big{(}\text{RSS}_{j}(\bm{\theta}_{S})/2\sigma_{j}^{2}+\log\sqrt{2\pi\sigma_{j}^{2}}\Big{)},\\\
\text{RSS}_{j}(\bm{\theta}_{S})&=\sum_{i=1}^{n}\Big{(}Y_{ij}-\sum_{(k,j)\in\mathcal{A}}U_{kj}Y_{ik}-\sum_{(l,j)\in\mathcal{C}}W_{lj}X_{il}\Big{)}^{2},\\\
\end{split}$ (5)

which involves $|\mathcal{A}|+|\mathcal{C}|$ parameters for $\bm{\theta}_{S}$.

Our plan is as follows. In Section 3.1, we hierarchically construct ${S}$ via
nodewise regressions of $Y_{j}$ over $\bm{X}$; $j=1,\ldots,p$ without the
acyclicity constraint for $\bm{U}$. Then in Section 3.2, on the ground of
constructed super-graph $S$ and (5), we develop likelihood ratio tests for
hypotheses (2) and (3).

###  3.1 Structure learning via peeling

This section develops a novel structure learning method based on a peeling
algorithm and nodewise constrained regressions to construct
$S=(\mathcal{A},\mathcal{C})$ in a hierarchical manner.

First, we observe an important connection between primary variables and
intervention variables. Rewrite (1) as

$\bm{Y}=\bm{V}^{\top}\bm{X}+\bm{\varepsilon}_{V},\quad\bm{\varepsilon}_{V}=(\bm{I}-\bm{U}^{\top})^{-1}\bm{\varepsilon}\sim
N(\bm{0},\bm{\Omega}^{-1}),$ (6)

where $\bm{\Omega}=(\bm{I}-\bm{U})\bm{\Sigma}^{-1}(\bm{I}-\bm{U}^{\top})$ is a
precision matrix and $\bm{V}=\bm{W}(\bm{I}-\bm{U})^{-1}$. Now, define the set
of candidate interventional relations $\mathcal{C}:=\\{(l,j):V_{lj}\neq 0\\}$.

###### Proposition 2

Suppose Assumption 1 is satisfied. Then $\mathcal{C}\supseteq\mathcal{I}$.
Moreover,

  1. (A)

if $V_{lj}\neq 0$, then $X_{l}$ intervenes on $Y_{j}$ or an ancestor of
$Y_{j}$;

  2. (B)

$Y_{j}$ is a leaf variable (having no child) if and only if there is an
instrument $X_{l}$ such that $V_{lj}\neq 0$ and $V_{lj^{\prime}}=0$ for
$j^{\prime}\neq j$.

The proof of Proposition 2 is deferred to Appendix B.2. Intuitively,
$V_{lj}\neq 0$ implies the dependence of $Y_{j}$ on $X_{l}$ through a directed
path $X_{l}\to\cdots\to Y_{j}$, and hence that $X_{l}$ intervenes on $Y_{j}$
or an ancestor of $Y_{j}$. Thus, the instruments on a leaf variable are
independent of the other primary variables conditional on the rest of
interventions. This observation suggests a method to reconstruct the DAG
topological order by sequentially identifying the leaves in the DAG and
removing the identified leaf variables.

Next, we discuss the estimation of $\bm{V}$ and construction of $S$.

####  3.1.1 Nodewise constrained regressions

We propose nodewise $\ell_{0}$-constrained regressions to estimate nonzero
elements of $\bm{V}$:

$\widehat{\bm{V}}_{\cdot j}=\mathop{\rm arg\,min\,}_{\bm{V}_{\cdot j}}\
\sum_{i=1}^{n}(Y_{ij}-\bm{V}_{\cdot
j}^{\top}\bm{X}_{i})^{2}\quad\text{s.t.}\quad\sum_{l=1}^{q}\operatorname{I}(V_{lj}\neq
0)\leq\kappa_{j};\quad j=1,\ldots,p,$ (7)

where $1\leq\kappa_{j}\leq q$ is an integer-valued tuning parameter
controlling the sparsity of $\bm{V}_{\cdot j}$ and can be chosen by BIC or
cross-validation.

To solve (7), we use $J_{\tau_{j}}(z)=\min(|z|/\tau_{j},1)$ as a surrogate of
the $\ell_{0}$-function $\operatorname{I}(z\neq 0)$ (Shen et al., 2012) and
develop a difference-of-convex (DC) algorithm with the $\ell_{0}$-projection
to improve globality of the solution of (7). At $(t+1)$th iteration, we solve
a weighted Lasso problem,

$\displaystyle\widetilde{\bm{V}}_{\cdot j}^{[t+1]}=\mathop{\rm
arg\,min\,}_{\bm{V}_{\cdot j}}\ \sum_{i=1}^{n}(Y_{ij}-\bm{V}_{\cdot
j}^{\top}\bm{X}_{i})^{2}+2n\gamma_{j}\tau_{j}\sum_{l=1}^{q}\operatorname{I}\big{(}|\widetilde{V}^{[t]}_{lj}|\leq\tau_{j}\big{)}|V_{lj}|;\quad
j=1,\ldots,p,$ (8)

where $\gamma_{j}>0$ is a tuning parameter and
$\widetilde{\bm{V}}^{[t]}_{\cdot j}$ is the solution of (8) at the $t$-th
iteration. The DC algorithm terminates at $\widetilde{\bm{V}}_{\cdot
j}=\widetilde{\bm{V}}_{\cdot j}^{[t]}$ such that $\|\widetilde{\bm{V}}_{\cdot
j}^{[t+1]}-\widetilde{\bm{V}}_{\cdot
j}^{[t]}\|_{\infty}\leq\sqrt{\text{tol}}$, where tol is the machine precision.
Then, we obtain the solution $\widehat{\bm{V}}_{\cdot j}$ of (7) by projecting
$\widetilde{\bm{V}}_{\cdot j}$ onto the $\ell_{0}$-constrained set
$\\{\|\bm{V}_{\cdot j}\|_{0}\leq\kappa_{j}\\}$.

1

2Specify $\gamma_{j}>0$, $\tau_{j}>0$, and $1\leq\kappa_{j}\leq q$. Initialize
$\widetilde{\bm{V}}_{\cdot j}^{[0]}$ with
$|\\{l:\widetilde{V}^{[0]}_{lj}>\tau_{j}\\}|\leq\kappa_{j}$.

3Compute $\widetilde{\bm{V}}_{\cdot j}=\widetilde{\bm{V}}^{[t]}_{\cdot j}$
such that $\|\widetilde{\bm{V}}_{\cdot j}^{[t+1]}-\widetilde{\bm{V}}_{\cdot
j}^{[t]}\|_{\infty}\leq\sqrt{\text{tol}}$.

4($\ell_{0}$-projection) Let
$B_{j}=\\{l:\sum_{l^{\prime}=1}^{q}\operatorname{I}(|\widetilde{V}_{l^{\prime}j}|\geq|\widetilde{V}_{lj}|)\leq\kappa_{j}\\}$.
Set  $\widehat{\bm{V}}_{\cdot j}=\mathop{\rm arg\,min\,}_{\bm{V}_{\cdot j}}\
\|\bm{\mathsf{Y}}_{j}-\bm{V}_{\cdot
j}^{\top}\bm{\mathsf{X}}\|_{2}^{2}\quad\text{s.t.}\quad V_{lj}=0\text{ for
}l\notin B_{j}.$

Algorithm 1 Constrained minimization via DC programming + $\ell_{0}$
projection

1

2Let $\widehat{\bm{V}}^{[0]}=\widehat{\bm{V}}$. Begin iteration with
$h=0,\ldots$: at iteration $h$,

3In the current sub-graph identify all instrument-leaf pairs $X_{l}\to Y_{j}$
satisfying $l=\mathop{\rm
arg\,min\,}_{l^{\prime}:\|\widehat{\bm{V}}^{[h]}_{l^{\prime}\cdot}\|_{0}\neq
0}\|\widehat{\bm{V}}^{[h]}_{l^{\prime}\cdot}\|_{0}$ and $j=\mathop{\rm
arg\,max\,}_{j^{\prime}}|\widehat{V}^{[h]}_{lj^{\prime}}|$. Let
$h_{\widehat{G}}(j)=h$.

4For a leaf $Y_{k}$, identify all directed relations $Y_{k}\to Y_{j}$ with
$h_{\widehat{G}}(j)=h-1$ when $\widehat{V}_{lj}\neq 0$ for each instrument
$X_{l}$ of $Y_{k}$ in the current sub-graph.

5Peeling-off all the leaves $Y_{j}$ by deleting $j$th columns from
$\widehat{\bm{V}}^{[h]}$ to yield $\widehat{\bm{V}}^{[h+1]}$.

6Repeat Steps 2-4 until all primary variables are removed.

7Compute $\widehat{S}=(\widehat{\mathcal{A}},\widehat{\mathcal{C}})$ and
$\widehat{\mathcal{D}}$:
$\begin{split}\widehat{\mathcal{A}}&=\Big{\\{}(k,j):Y_{k}\to\cdots\to
Y_{j}\Big{\\}},\quad\widehat{\mathcal{D}}=\Big{\\{}(k,j)\in\mathcal{H}:(j,k)\notin\widehat{\mathcal{A}}\Big{\\}},\\\
\widehat{\mathcal{C}}&=\Big{\\{}(l,j):|\widehat{V}_{lj}|>\tau_{j}\Big{\\}},\text{
where $\tau_{j}$ is specified in Algorithm \ref{algorithm:nodewise-
regression}.}\end{split}$

Algorithm 2 Reconstruction of topological layers by peeling

####  3.1.2 Peeling

Now, we describe a _peeling_ algorithm to estimate
$S=(\mathcal{A},\mathcal{C})$ based on $\bm{V}$.

Before proceeding, we introduce the notion of _topological height_ to
characterize the DAG topological order. Given a DAG $G$, the height $h_{G}(j)$
of $Y_{j}$ is defined as the maximal length of a directed path from $Y_{j}$ to
a leaf variable, where $0\leq h_{G}(j)\leq p-1$ with a leaf having zero
height. Then the primary variables of the same height form a topological layer
and the layers constitute a partition of primary variables in $G$. For
instance, in Figure 1, $\\{Y_{4},Y_{5}\\},\\{Y_{3}\\},\\{Y_{2}\\},\\{Y_{1}\\}$
are topological layers with height $0,1,2,3$, respectively.

Observe that the topological layer of height $h$ are leaves in the sub-graph
with primary variables of height no less than $h$. Based on Proposition 2, the
peeling algorithm sequentially identifies and peels off the primary variables
of height $h$; $h=0,1,2,\ldots$. Note that instruments in the sub-graph may
not be necessarily the instruments in the graph $G$ after some primary
variables are removed. However, the next proposition indicates that the
instruments in the corresponding sub-graph remain useful to reconstruct
ancestral relations.

###### Proposition 3

Let $Y_{k}$ have height $h_{G}(k)$. Let $G^{\prime}$ be the sub-graph with
primary variables of height $\geq h_{G}(k)$. Under Assumption 1, for a pair
$(k,j)$ with $h_{G}(j)=h_{G}(k)-1$, we have $Y_{k}\to Y_{j}$ if and only if
$V_{lj}\neq 0$ for each instrument $X_{l}$ of $Y_{k}$ in $G^{\prime}$.

Proposition 3 (proved in Appendix B.3) permits the construction of directed
relations between the topological layers of heights $h$ and $h-1$; $h\geq 1$
and thus ancestral relations.

The peeling algorithm is summarized in Algorithm 2 and a detailed illustration
is presented in Example 3 of Appendix A.2.

######  Remark 2 (Recovery of true graph $G$)

$\widehat{\mathcal{A}}$ estimates a superset of directed relations of $\bm{U}$
and $\widehat{\mathcal{C}}$ a superset of interventional relations of
$\bm{W}$. Therefore, under some technical conditions, with suitable tuning
parameters, the true graph $G$ can be reconstructed by estimating
$(\bm{U},\bm{W})$ through $\ell_{0}$-constrained regressions given the
estimated super-graph
$\widehat{S}=(\widehat{\mathcal{A}},\widehat{\mathcal{C}})$,

$\begin{split}\min_{(\bm{U}_{\textnormal{{an}}_{\widehat{S}}(j),j},\bm{W}_{\textnormal{{in}}_{\widehat{S}}(j),j})}\sum_{i=1}^{n}\left(Y_{ij}-\bm{U}_{\textnormal{{an}}_{\widehat{S}}(j),j}^{\top}\bm{Y}_{i,\textnormal{{an}}_{\widehat{S}}(j)}-\bm{W}_{\textnormal{{in}}_{\widehat{S}}(j),j}^{\top}\bm{X}_{i,\textnormal{{in}}_{\widehat{S}}(j)}\right)^{2}\\\
\text{s.t.}\quad\sum_{k\in\textnormal{{an}}_{\widehat{S}}(j)}\operatorname{I}(U_{kj}\neq
0)+\sum_{l\in\textnormal{{in}}_{\widehat{S}}(j)}\operatorname{I}(W_{lj}\neq
0)\leq\kappa_{j}^{\prime}.\end{split}$ (9)

The final estimates $\widehat{\bm{U}}_{\cdot
j}=\Big{(}\widehat{\bm{U}}_{\textnormal{{an}}_{\widehat{S}}(j),j},\bm{0}_{\textnormal{{an}}_{\widehat{S}}(j)^{c}}\Big{)}$
and $\widehat{\bm{W}}_{\cdot
j}=\Big{(}\widehat{\bm{W}}_{{\textnormal{{in}}}_{\widehat{S}}(j),j},\bm{0}_{{\textnormal{{in}}}_{\widehat{S}}(j)^{c}}\Big{)}$;
$j=1,\ldots,p$.

###  3.2 Likelihood ratio test via data perturbation

This subsection proposes an inference method for testing (2) and (3). We first
construct a likelihood ratio $\text{Lr}(S,\bm{\Sigma})$ in given
$(S,\bm{\Sigma})$ and then plug-in an estimate
$(\widehat{S},\widehat{\bm{\Sigma}})$. Next we perform a test via data
perturbation, accounting for the uncertainty of constructing $S$.

####  3.2.1 Likelihood ratio given $S$

Consider the test in (2). From (5), the likelihood ratio given
$(S,\bm{\Sigma})$ is

$\text{Lr}(S,\bm{\Sigma})=L_{S}(\widehat{\bm{\theta}}_{S}^{(1)},\bm{\Sigma})-L_{S}(\widehat{\bm{\theta}}_{S}^{(0)},\bm{\Sigma})=\sum_{j=1}^{p}\frac{\text{RSS}_{j}(\widehat{\bm{\theta}}^{(0)}_{S})-\text{RSS}_{j}(\widehat{\bm{\theta}}^{(1)}_{S})}{2\sigma_{j}^{2}}$

where $\widehat{\bm{\theta}}_{S}^{(1)}$ is the MLE given $(S,\bm{\Sigma})$,
and $\widehat{\bm{\theta}}_{S}^{(0)}$ is the MLE given $(S,\bm{\Sigma})$
subject to $\bm{U}_{\mathcal{H}}=0$.

Now, recall that $\textnormal{{d}}_{G}(j)=\\{k:(k,j)\in\mathcal{D}\\}$, where
$\mathcal{D}$ is the set of nondegenerate edges of $H_{0}$ with respect to $G$
(Definition 2). Of note, $\textnormal{{d}}_{S}(j)=\textnormal{{d}}_{G}(j)$
since an edge is nondegenerate with respect to $S$ if and only if it is
nondegenerate with respect to $G$. To eliminate the nuisance parameters,
observe that if $\textnormal{{d}}_{S}(j)=\emptyset$, then
${\text{RSS}}_{j}(\widehat{\bm{\theta}}^{(0)}_{S})-{\text{RSS}}_{j}(\widehat{\bm{\theta}}^{(1)}_{S})=0$,
because no $H_{0}$ constraint is imposed for $Y_{j}$. Hence,
$\text{Lr}(S,\bm{\Sigma})$ only summarizes the contributions from the primary
variables with nondegenerate hypothesized edges,

$\text{Lr}(S,\bm{\Sigma})=\sum_{j:\textnormal{{d}}_{S}(j)\neq\emptyset}\frac{{\text{RSS}}_{j}(\widehat{\bm{\theta}}^{(0)}_{S})-{\text{RSS}}_{j}(\widehat{\bm{\theta}}^{(1)}_{S})}{2\sigma_{j}^{2}},$
(10)

where the nuisance parameters for $Y_{j}$ with
$\textnormal{{d}}_{S}(j)=\emptyset$ are eliminated.

The likelihood ratio test statistic becomes
$\text{Lr}=\text{Lr}(\widehat{S},\widehat{\bm{\Sigma}})$. Here, $\widehat{S}$
is constructed by Algorithm 2 in Section 3.1 and
$\widehat{\bm{\Sigma}}=\text{Diag}(\widehat{\sigma}^{2}_{1},\ldots,\widehat{\sigma}^{2}_{p})$,

$\displaystyle\widehat{\sigma}^{2}_{j}=\frac{\text{RSS}_{j}(\widehat{\bm{\theta}}^{(1)}_{\widehat{S}})}{n-|{\textnormal{{an}}}_{\widehat{S}}(j)|-|{\textnormal{{in}}}_{\widehat{S}}(j)|},\quad
j=1,\ldots,p,$ (11)

where ${\textnormal{{an}}}_{\widehat{S}}(j)$ and
${\textnormal{{in}}}_{\widehat{S}}(j)$ are the estimated ancestor and
candidate intervention sets of $Y_{j}$.

####  3.2.2 Testing directed edges (2) via data perturbation

The likelihood ratio requires an estimation of $S$, where we must account for
the uncertainty of $\widehat{S}$ for accurate finite-sample inference.

To proceed, we consider the test statistic Lr based on a “correct” super-graph
$S\supseteq S^{\circ}$, where $S=(\mathcal{A},\mathcal{C})\supseteq
S^{\circ}=(\mathcal{A}^{\circ},\mathcal{C}^{\circ})$ means that
$\mathcal{A}\supseteq\mathcal{A}^{\circ}$ and
$\mathcal{C}\supseteq\mathcal{C}^{\circ}$ and ∘ denotes the truth.
Intuitively, a “correct” super-graph distinguishes descendants and
nondescendants and thus can help to infer the true directed relations defined
by the local Markov property (Spirtes et al., 2000) without introducing model
errors, yet may lead to a less powerful test when $S$ is much larger than
$S^{\circ}$. By comparison, a “wrong” graph $S\not\supseteq S^{\circ}$
provides an incorrect topological order, and a test based on a “wrong” graph
may be biased, accompanied by an inflated type-I error.

Now, let $\bm{\mathsf{Z}}_{n\times(p+q)}=(\bm{\mathsf{X}},\bm{\mathsf{Y}})$ be
the data matrix, where $\bm{\mathsf{Y}}_{n\times
p}=(\bm{Y}_{1},\ldots,\bm{Y}_{n})^{\top}$ and $\bm{\mathsf{X}}_{n\times
q}=(\bm{X}_{1},\ldots,\bm{X}_{n})^{\top}$. Moreover, let
$\bm{\mathsf{e}}_{j}\in\mathbb{R}^{n}$ be the $j$th column of error matrix
$\bm{\mathsf{e}}_{n\times
p}=(\bm{\varepsilon}_{1},\ldots,\bm{\varepsilon}_{n})^{\top}$. From (10),
assuming $\widehat{S}\supseteq S^{\circ}$ is a “correct” super-graph, the
likelihood ratio becomes

$\begin{split}\text{Lr}&=\sum_{j:{\textnormal{{d}}}_{\widehat{S}}(j)\neq\emptyset}\frac{\bm{\mathsf{Y}}_{j}^{\top}(\bm{\mathsf{P}}_{\widehat{A}_{j}}-\bm{\mathsf{P}}_{\widehat{B}_{j}})\bm{\mathsf{Y}}_{j}}{2\widehat{\sigma}_{j}}\\\
&=\sum_{j:{\textnormal{{d}}}_{\widehat{S}}(j)\neq\emptyset}\frac{(\bm{\mathsf{Y}}_{{\textnormal{{d}}}_{\widehat{S}}(j)}\bm{U}_{{\textnormal{{d}}}_{\widehat{S}}(j),j}+\bm{\mathsf{e}}_{j})^{\top}(\bm{\mathsf{P}}_{\widehat{A}_{j}}-\bm{\mathsf{P}}_{\widehat{B}_{j}})(\bm{\mathsf{Y}}_{{\textnormal{{d}}}_{\widehat{S}}(j)}\bm{U}_{{\textnormal{{d}}}_{\widehat{S}}(j),j}+\bm{\mathsf{e}}_{j})}{2\widehat{\sigma}_{j}},\\\
\widehat{\sigma}_{j}&=\frac{\bm{\mathsf{Y}}_{j}^{\top}(\bm{I}-\bm{\mathsf{P}}_{\widehat{A}_{j}})\bm{\mathsf{Y}}_{j}}{n-|\widehat{A}_{j}|}=\frac{\bm{\mathsf{e}}_{j}^{\top}(\bm{I}-\bm{\mathsf{P}}_{\widehat{A}_{j}})\bm{\mathsf{e}}_{j}}{n-|\widehat{A}_{j}|};\quad
1\leq j\leq p,\end{split}$ (12)

where
$\bm{\mathsf{P}}_{A}=\bm{\mathsf{Z}}_{A}(\bm{\mathsf{Z}}_{A}^{\top}\bm{\mathsf{Z}}_{A})^{-1}\bm{\mathsf{Z}}_{A}^{\top}$
is the projection matrix onto the span of columns $A$ in $\bm{\mathsf{Z}}$,
$\widehat{A}_{j}={\textnormal{{an}}}_{\widehat{S}}(j)\cup{\textnormal{{in}}}_{\widehat{S}}(j)\cup{\textnormal{{d}}}_{\widehat{S}}(j)$,
and
$\widehat{B}_{j}=({\textnormal{{an}}}_{\widehat{S}}(j)\cup{\textnormal{{in}}}_{\widehat{S}}(j))\setminus{\textnormal{{d}}}_{\widehat{S}}(j)$.
Of note,
$\bm{\mathsf{Y}}_{{\textnormal{{d}}}_{S}(j)}\bm{U}_{{\textnormal{{d}}}_{S}(j),j}=\bm{0}$
for all $j$ under $H_{0}$, while
$\bm{\mathsf{Y}}_{{\textnormal{{d}}}_{S}(j)}\bm{U}_{{\textnormal{{d}}}_{S}(j),j}\neq\bm{0}$
for some $j$ under $H_{a}$.

We propose a data perturbation (DP) method (Shen and Ye, 2002; Breiman, 1992)
to approximate the null distribution of Lr in (12). The idea behind DP is to
assess the sensitivity of the estimates through perturbed data
$\bm{Y}^{*}_{i}=\bm{Y}_{i}+\bm{\varepsilon}^{*}_{i}$, where
$\bm{\varepsilon}_{i}^{*}\sim N(0,\widehat{\bm{\Sigma}})$; $i=1,\ldots,n$.

Let
$(\bm{\mathsf{Z}}^{*},\bm{\mathsf{e}}^{*})=(\bm{\mathsf{X}},\bm{\mathsf{Y}}^{*},\bm{\mathsf{e}}^{*})$
be perturbed data, where $\bm{\mathsf{e}}^{*}_{n\times
p}=(\bm{\mathsf{e}}^{*}_{1},\ldots,\bm{\mathsf{e}}^{*}_{p})=(\bm{\varepsilon}^{*}_{1},\ldots,\bm{\varepsilon}^{*}_{n})^{\top}$
denotes perturbation errors. Given
$(\bm{\mathsf{Z}}^{*},\bm{\mathsf{e}}^{*})$, we obtain a perturbation estimate
$\widehat{S}^{*}$. In (12), under $H_{0}$,
$\text{Lr}=f(\bm{\mathsf{X}},\bm{\mu},\bm{\mathsf{e}})$ is a function of
intervention data $\bm{\mathsf{X}}$, conditional mean
$\bm{\mu}=\operatorname{\mathbb{E}}(\bm{\mathsf{Y}}\mid\bm{\mathsf{X}})\in\mathbb{R}^{n\times
p}$ and unobserved error $\bm{\mathsf{e}}$. However, the perturbation error
$\bm{\mathsf{e}}^{*}$ is _known_ given the DP sample
$(\bm{\mathsf{Z}}^{*},\bm{\mathsf{e}}^{*})$, suggesting a form $\text{Lr}^{*}$
based on Lr:

$\begin{split}\text{Lr}^{*}&\equiv
f(\bm{\mathsf{X}},\bm{\mathsf{Y}},\bm{\mathsf{e}}^{*})=\sum_{j:{\textnormal{{d}}}_{\widehat{S}}(j)\neq\emptyset}\frac{(\bm{\mathsf{e}}_{j}^{*})^{\top}(\bm{\mathsf{P}}_{\widehat{A}^{*}_{j}}-\bm{\mathsf{P}}_{\widehat{B}^{*}_{j}})\bm{\mathsf{e}}^{*}_{j}}{2\widetilde{\sigma}_{j}^{*}},\quad\widetilde{\sigma}^{*}_{j}=\frac{(\bm{\mathsf{e}}_{j}^{*})^{\top}(\bm{I}-\bm{\mathsf{P}}_{\widehat{A}^{*}_{j}})\bm{\mathsf{e}}^{*}_{j}}{n-|\widehat{A}^{*}_{j}|},\quad
1\leq j\leq p,\end{split}$ (13)

which mimics (12) when
$\bm{\mathsf{Y}}_{{\textnormal{{d}}}_{S}(j)}\bm{U}_{{\textnormal{{d}}}_{S}(j)}=\bm{0}$
under $H_{0}$, with $\bm{\mathsf{e}}_{j}$ replaced by
$\bm{\mathsf{e}}^{*}_{j}$. As a result, when $\widehat{S}^{*}\supseteq
S^{\circ}$, the conditional distribution of $\text{Lr}^{*}$ given
$\bm{\mathsf{Z}}$ well approximates the null distribution of Lr, where the
model selection effect is accounted for by assessing the variability of
$\widehat{A}^{*}_{j}$ and $\widehat{B}^{*}_{j}$ over different realizations of
$(\bm{\mathsf{Z}}^{*},\bm{\mathsf{e}}^{*})$.

In practice, we use Monte-Carlo to approximate the distribution of
$\text{Lr}^{*}$ given $\bm{\mathsf{Z}}$ in (13). In particular, we generate
$M$ perturbed samples
$(\bm{\mathsf{Z}}^{*}_{m},\bm{\mathsf{e}}_{m}^{*})_{m=1}^{M}$ independently
and compute $\text{Lr}^{*}_{m}$; $m=1,\ldots,M$. Then, we examine the
condition $\widehat{S}^{*}_{m}\supseteq S^{\circ}$ by checking its empirical
counterpart $\widehat{S}^{*}_{m}\supseteq\widehat{S}$. The DP p-value of edge
test in (2) is defined as

$\text{Pval}=\left(\sum_{m=1}^{M}\operatorname{I}(\text{Lr}_{m}^{*}\geq\text{Lr},\widehat{S}^{*}_{m}\supseteq\widehat{S})\right)/\left(\sum_{m=1}^{M}\operatorname{I}(\widehat{S}^{*}_{m}\supseteq\widehat{S})\right),$
(14)

where $\operatorname{I}(\cdot)$ is the indicator function.

###### Remark 3

Assume $\widehat{S}^{*}\supseteq S^{\circ}$ is a “correct” super-graph. A
naive likelihood ratio based on perturbed data
$(\bm{\mathsf{Z}}^{*},\bm{\mathsf{e}}^{*})$ is

$\text{Lr}(\widehat{S}^{*},\widehat{\bm{\Sigma}}^{*})=\sum_{j:{\textnormal{{d}}}_{\widehat{S}}(j)\neq\emptyset}\frac{(\bm{\mathsf{Y}}_{{\textnormal{{d}}}_{\widehat{S}}(j)}\bm{U}_{{\textnormal{{d}}}_{\widehat{S}}(j),j}+\bm{\mathsf{e}}_{j}+\bm{\mathsf{e}}^{*}_{j})^{\top}(\bm{\mathsf{P}}_{\widehat{A}^{*}_{j}}-\bm{\mathsf{P}}_{\widehat{B}^{*}_{j}})(\bm{\mathsf{Y}}_{{\textnormal{{d}}}_{\widehat{S}}(j)}\bm{U}_{{\textnormal{{d}}}_{\widehat{S}}(j),j}+\bm{\mathsf{e}}_{j}+\bm{\mathsf{e}}^{*}_{j})}{2\widehat{\sigma}_{j}^{*}}.$

Note that the error $\bm{\mathsf{e}}_{j}$ given $\bm{\mathsf{Z}}$ is
deterministic and does not vanish under either $H_{0}$ or $H_{a}$;
$j=1,\ldots,p$. Thus, the conditional distribution of
$\text{Lr}(\widehat{S}^{*},\widehat{\bm{\Sigma}}^{*})$ given $\bm{\mathsf{Z}}$
does not approximate the null distribution of Lr in (12), in contrast to the
DP likelihood ratio $\text{Lr}^{*}$ in (13).

####  3.2.3 Extension to hypothesis testing for a pathway

Next, we extend the DP inference for (3). Denote
$\mathcal{H}=\\{(k_{1},j_{1}),\ldots,(k_{|\mathcal{H}|},j_{|\mathcal{H}|})\\}$.
Then the test of pathways in (3) can be reduced to testing sub-hypotheses

$H_{0,\nu}:U_{k_{\nu},j_{\nu}}=0,\quad\mbox{ versus
}\quad{H}_{a,\nu}:U_{k_{\nu},j_{\nu}}\neq 0;\quad\nu=1,\ldots,|\mathcal{H}|,$

where each sub-hypothesis is a directed edge test. Given $(S,\bm{\Sigma})$,
the likelihood ratio for $\text{H}_{0,\nu}$ is
$\text{Lr}_{\nu}(S,\bm{\Sigma})=L_{S}(\widehat{\bm{\theta}}_{S}^{(1)},\bm{\Sigma})-L_{S}(\widehat{\bm{\theta}}_{S}^{(0,\nu)},\bm{\Sigma})$,
where $\widehat{\bm{\theta}}_{S}^{(0,\nu)}$ is the MLE under the constraint
that $U_{k_{\nu},j_{\nu}}=0$. Then for $\widehat{S}\supseteq S^{\circ}$, we
have

$\begin{split}\text{Lr}_{\nu}&=\frac{1}{2}\frac{(\bm{\mathsf{Y}}_{k_{\nu}}U_{k_{\nu},j_{\nu}}+\bm{\mathsf{e}}_{j_{\nu}})^{\top}(\bm{\mathsf{P}}_{\widehat{A}_{j_{\nu}}}-\bm{\mathsf{P}}_{\widehat{B}_{j_{\nu}}})(\bm{\mathsf{Y}}_{k_{\nu}}U_{k_{\nu},j_{\nu}}+\bm{\mathsf{e}}_{j_{\nu}})}{(\bm{\mathsf{e}}_{j_{\nu}})^{\top}(\bm{I}-\bm{\mathsf{P}}_{\widehat{A}_{j_{\nu}}})\bm{\mathsf{e}}_{j_{\nu}}/(n-|\widehat{A}_{j_{\nu}}|)},\\\
\text{Lr}_{\nu}^{*}&=\frac{1}{2}\frac{(\bm{\mathsf{e}}^{*}_{j_{\nu}})^{\top}(\bm{\mathsf{P}}_{\widehat{A}^{*}_{j_{\nu}}}-\bm{\mathsf{P}}_{\widehat{B}^{*}_{j_{\nu}}})\bm{\mathsf{e}}^{*}_{j_{\nu}}}{(\bm{\mathsf{e}}^{*}_{j_{\nu}})^{\top}(\bm{I}-\bm{\mathsf{P}}_{\widehat{A}^{*}_{j_{\nu}}})\bm{\mathsf{e}}^{*}_{j_{\nu}}/(n-|\widehat{A}^{*}_{j_{\nu}}|)},\end{split}$
(15)

where the distributions of $\text{Lr}_{\nu}^{*}$ given $\bm{\mathsf{Z}}$
approximates the null distributions of $\text{Lr}_{\nu}$. Finally, define the
p-value of pathway test in (3) as

$\text{Pval}=\max_{1\leq\nu\leq|\mathcal{H}|}\left(\sum_{m=1}^{M}\operatorname{I}(\text{Lr}_{\nu,m}^{*}\geq\text{Lr}_{\nu},\widehat{S}^{*}_{m}\supseteq\widehat{S})\right)/\left(\sum_{i=1}^{M}\operatorname{I}(\widehat{S}^{*}_{m}\supseteq\widehat{S})\right).$
(16)

Note that if any $\text{H}_{0,\nu}$ is degenerate, then $\text{Pval}=1$.

Algorithm 3 summarizes the DP method for hypothesis testing.

1Specify the Monte Carlo size $M$.

2Compute Lr and
$\widehat{\bm{\Sigma}}=\text{Diag}(\widehat{\sigma}^{2}_{1},\ldots,\widehat{\sigma}^{2}_{p})$
via Algorithm 2 with original data $\bm{\mathsf{Z}}$.

3(Parallel DP) Generate perturbed data
$(\bm{\mathsf{Z}}^{*}_{m},\bm{\mathsf{e}}_{m}^{*})=(\bm{\mathsf{X}},\bm{\mathsf{Y}}_{m}^{*},\bm{\mathsf{e}}_{m}^{*})$
in parallel, where $\bm{\varepsilon}_{m,i}^{*}\sim
N(\bm{0},\widehat{\bm{\Sigma}})$. For (2), compute $\text{Lr}^{*}_{m}$ based
on $(\bm{\mathsf{Z}}^{*}_{m},\bm{\mathsf{e}}_{m}^{*})$, while, for (3),
compute $\text{Lr}^{*}_{\nu,m}$; $\nu=1,\ldots,|\mathcal{H}|$ based on
$(\bm{\mathsf{Z}}^{*}_{m},\bm{\mathsf{e}}_{m}^{*})$; $m=1,\ldots,M$.

4Compute the p-value of a test as (14) or (16) accordingly.

Algorithm 3 DP likelihood ratio test

######  Remark 4 (Computation)

To accelerate the computation, we consider a parallelizable procedure in
Algorithm 3. In addition, we use the original estimate $\widehat{\bm{\theta}}$
as a warm-start initialization for the DP estimates, which effectively reduces
the computing time.

######  Remark 5 (Connection with bootstrap)

One may consider parametric or nonparametric bootstrap for Lr. The parametric
bootstrap requires a good initial estimate of $(\bm{U},\bm{W})$. Yet, it is
rather challenging to correct the bias of this estimate because of the
acyclicity constraint. By comparison, DP does not rely on such an estimate. On
the other hand, nonparametric bootstrap resamples the original data with
replacement. In a bootstrap sample, only about 63% distinct observations in
the original data are used in model selection and fitting, leading to
deteriorating performance (Kleiner et al., 2012), especially in a small
sample. As a result, nonparametric bootstrap may not well-approximate the
distribution of Lr, while DP provides a better approximation of Lr, taking
advantage of a full sample.

##  4 Theory

This section provides the theoretical justification for the proposed methods.

###  4.1 Convergence and consistency of structure learning

First, we introduce some technical assumptions to derive statistical and
computational properties of Algorithms 1 and 2. Let
$(\bm{\mathsf{Y}},\bm{\mathsf{X}})\in\mathbb{R}^{n\times(p+q)}$ be the data
matrix. Denote by $\bm{\zeta}$ a generic vector and
$\bm{\zeta}_{A}\in\mathbb{R}^{|A|}$ be the projection of $\bm{\zeta}$ onto
coordinates in $A$. Let $\kappa^{\circ}_{j}=\|\bm{V}^{\circ}_{\cdot j}\|_{0}$
and $\kappa^{\circ}_{\max}=\max_{1\leq j\leq p}\kappa_{j}^{\circ}$.

######  Assumption 2 (Restricted eigenvalues)

For a constant $c_{1}>0$,

$\min_{A:|A|\leq
2\kappa^{\circ}_{\max}}\min_{\\{\bm{\zeta}:\|\bm{\zeta}_{A^{c}}\|_{1}\leq
3\|\bm{\zeta}_{A}\|_{1}\\}}\frac{\|\bm{\mathsf{X}}\bm{\zeta}\|_{2}^{2}}{n\|\bm{\zeta}\|_{2}^{2}}\geq
c_{1}.$

######  Assumption 3 (Interventions)

For a constant $c_{2}>0$,

$\max\Big{(}\max_{1\leq l\leq
q}n^{-1}(\bm{\mathsf{X}}^{\top}\bm{\mathsf{X}})_{ll},\max_{A:|A|\leq
2\kappa^{\circ}_{\max}}\max_{1\leq
l\leq|A|}n((\bm{\mathsf{X}}_{A}^{\top}\bm{\mathsf{X}}_{A})^{-1})_{ll}\Big{)}\leq
c_{2}^{2}.$

######  Assumption 4 (Nuisance signals)

$\min_{V^{0}_{lj}\neq 0}|V^{0}_{lj}|\geq 100\frac{c_{2}}{c_{1}}\max_{1\leq
j\leq
p}\sqrt{\Omega_{jj}^{-1}\Big{(}\frac{\log(q)}{n}+\frac{\log(n)}{n}\Big{)}}.$
(17)

Assumptions 2-3, as a replacement of Assumption 1A, are satisfied with a
probability tending to one for isotropic subgaussian or bounded $\bm{X}$
(Rudelson and Zhou, 2012). Assumption 2 is a common condition for proving the
convergence rate of Lasso regression (Bickel et al., 2009; Zhang et al.,
2014). Assumption 4, as an alternative to Assumption 1B, specifies the minimal
signal strength over candidate interventions. A signal strength requirement
like Assumption 4 is a condition for establishing selection consistency in
high dimensional variable selection (Fan et al., 2014; Loh et al., 2017; Zhao
et al., 2018). Moreover, Assumption 4 can be further relaxed to a less
intuitive condition, Assumption 6; see Appendix A.3 for details.

######  Theorem 4 (Finite termination and consistency)

Under Assumptions 1-4 with constants $c_{1}<6c_{2}$, if the machine precision
$\text{tol}\ll 1/n$ is negligible, and tuning parameters of Algorithm 1
satisfy:

  1. (1)

$\gamma_{j}\leq c_{1}/6$,

  2. (2)

$\sqrt{32c_{2}^{2}\Omega_{jj}^{-1}(\log(q)+\log(n))/n}/\gamma_{j}\leq\tau_{j}\leq
0.4\min_{V^{0}_{lj}\neq 0}|V^{0}_{lj}|$,

  3. (3)

$\kappa_{j}=\kappa_{j}^{\circ}$,

for $j=1,\ldots,p$, then Algorithm 1 yields the global minimizer
$\widehat{\bm{V}}_{\cdot j}$ of (7) in at most
$1+\lceil\log(\kappa_{\max}^{\circ})/\log 4\rceil$ DC iterations almost
surely, where $\lceil\cdot\rceil$ is the ceiling function. Moreover, Algorithm
2 recovers $\mathcal{A}^{\circ}$,
$\mathcal{C}^{\circ}\equiv\\{l:V^{\circ}_{lj}\neq
0\\}\supseteq\mathcal{I}^{\circ}$, and $\mathcal{D}^{\circ}$ almost surely as
$n\to\infty$.

By Theorem 4 (proved in Appendix B.4), the computational complexity of
Algorithm 1 is the number of iterations multiplied by that of solving a
weighted Lasso regression, $O(\log\kappa^{\circ}_{\max})\times
O(q^{3}+nq^{2})$ (Efron et al., 2004). Also note that in Steps 2-6 of
Algorithm 2, no heavy computation is involved. Thus, the overall time
complexity of Algorithm 2 is bounded by
$O(p\times\log\kappa^{\circ}_{\max}\times(q^{3}+nq^{2}))$ for computing
$\widehat{\bm{V}}$ matrix. Finally, the peeling algorithm does not apply to
observational data ($\bm{W}=\bm{0}$). In a sense, interventions play an
essential role.

Although Theorem 4 establishes the consistent reconstruction of the peeling
algorithm for the super-graph $S$, it does not provide any uncertainty measure
of the presence of some directed edges of the true DAG. In what follows, we
will develop a theory for hypothesis tests to make valid inference concerning
directed edges of interest.

###  4.2 Inferential theory

######  Assumption 5 (Hypothesis-specific dimension restriction)

For a constant $\rho$,

$\max_{j:\textnormal{{d}}_{G^{\circ}}(j)\neq\emptyset}\frac{|\textnormal{{an}}_{G^{\circ}}(j)|+|\textnormal{{d}}_{G^{\circ}}(j)|+{\kappa_{j}^{\circ}}}{n}\leq\rho<1,\quad\text{
as }n\rightarrow\infty.$

Assumption 5 is a hypothesis-specific condition restricting the underlying
dimension of the testing problem. Usually,
$|\textnormal{{an}}_{G^{\circ}}(j)|\asymp\kappa_{j}^{\circ}\ll p$;
$j=1,\ldots,p$, which dramatically relaxes the condition
$n\gg{|\mathcal{D}|^{1/2}p\log(p)}$ for the constrained likelihood ratio test
(Li et al., 2019).

######  Theorem 5 (Empirical p-values)

Suppose Assumptions 1-5 are met and $H_{0}$ is regular. If the tuning
parameters in Algorithm 1 satisfy (1)-(3) in Theorem 4, then,

  1. (A)

(Test of directed edges (2))

$\begin{split}\lim_{\begin{subarray}{c}n\to\infty\\\
\bm{\theta}^{\circ}\textnormal{ satisfies H}_{0}\textnormal{ in
}\eqref{equation:link-
test}\end{subarray}}\lim_{M\rightarrow\infty}\operatorname{\mathbb{P}}_{\bm{\theta}^{\circ}}(\textnormal{Pval}<\alpha)=\alpha,&\quad\text{
if }H_{0}\text{ is nondegenerate};\\\ \lim_{\begin{subarray}{c}n\to\infty\\\
\bm{\theta}^{\circ}\textnormal{ satisfies H}_{0}\textnormal{ in
}\eqref{equation:link-
test}\end{subarray}}\lim_{M\rightarrow\infty}\operatorname{\mathbb{P}}_{\bm{\theta}^{\circ}}(\textnormal{Pval}=1)=1,&\quad\text{
if }H_{0}\text{ is degenerate}.\end{split}$

  2. (B)

(Test of directed pathways (3))

$\begin{split}\limsup_{\begin{subarray}{c}n\to\infty\\\
\bm{\theta}^{\circ}\textnormal{ satisfies H}_{0}\textnormal{ in
}\eqref{equation:pathway-
test}\end{subarray}}\lim_{M\rightarrow\infty}\operatorname{\mathbb{P}}_{\bm{\theta}^{\circ}}(\textnormal{Pval}<\alpha)=\alpha,&\quad\text{
if }H_{0}\text{ is nondegenerate with }|\mathcal{D}|=|\mathcal{H}|;\\\
\limsup_{\begin{subarray}{c}n\to\infty\\\ \bm{\theta}^{\circ}\textnormal{
satisfies H}_{0}\textnormal{ in }\eqref{equation:pathway-
test}\end{subarray}}\lim_{M\rightarrow\infty}\operatorname{\mathbb{P}}_{\bm{\theta}^{\circ}}(\textnormal{Pval}=1)=1,&\quad\text{
if }|\mathcal{D}|<|\mathcal{H}|.\end{split}$

By Theorem 5, the DP likelihood ratio test yields a valid p-value for (2) and
(3). Note that $|\mathcal{D}|$ is permitted to depend on $n$. Moreover,
Proposition 6 summarizes the asymptotics for directed edge test (2).

######  Proposition 6 (Asymptotic distribution of Lr of (2))

Assume that the assumptions of Theorem 5 are met. Under a regular $H_{0}$, as
$n\to\infty$,

$\begin{split}2\textnormal{Lr}\overset{d}{\longrightarrow}\chi^{2}_{|\mathcal{D}|},&\text{
if }H_{0}\text{ is nondegenerate with }|\mathcal{D}|>0\text{ being fixed};\\\
\frac{2\textnormal{Lr}-|\mathcal{D}|}{\sqrt{2|\mathcal{D}|}}\overset{d}{\longrightarrow}N(0,1),&\text{
if }H_{0}\text{ is nondegenerate with }|\mathcal{D}|\rightarrow\infty,\
\frac{|\mathcal{D}|\log|\mathcal{D}|}{n}\rightarrow 0.\end{split}$

The proofs of Theorem 5 and Proposition 6 are given in Appendix B.5 and B.6.

######  Remark 6 (Local structure for hypothesis testing)

Theorem 5 or Proposition 6 only requires correct identification of the local
structures of the DAG such as,
$\textnormal{{an}}_{S}(j)=\textnormal{{an}}_{G}(j)$,
$\textnormal{{d}}_{S}(j)=\textnormal{{d}}_{G}(j)$, and
$\textnormal{{in}}_{S}(j)$, for variables $Y_{j}$ with
$\textnormal{{d}}_{S}(j)\neq\emptyset$, as opposed to those of entire $S$.

Next, we analyze the local limit power of the proposed tests for (2) and (3).
Assume $\bm{\theta}^{\circ}=(\bm{U}^{\circ},\bm{W}^{\circ})$ satisfies
$H_{0}$. Let $\bm{\Delta}\in\mathbb{R}^{p\times p}$ satisfy
$\bm{\Delta}_{\mathcal{D}^{c}}=\bm{0}$ so that $\bm{U}^{\circ}+\bm{\Delta}$
represents a DAG. For nondegenerate and regular $H_{0}$, consider an
alternative $H_{a}$:
$\bm{U}_{\mathcal{H}}=\bm{U}_{\mathcal{H}}^{\circ}+\bm{\Delta}_{\mathcal{H}}$,
and define the power function as

$\beta(\bm{\theta}^{\circ},\bm{\Delta})=\operatorname{\mathbb{P}}_{H_{a}}(\text{Pval}<\alpha).$
(18)

######  Proposition 7 (Local power of edge test)

Assume that $H_{0}$ is nondegenerate and regular. Let
$\|\bm{\Delta}\|_{F}=\|\bm{\Delta}_{\mathcal{H}}\|_{F}=n^{-1/2}\delta$ when
$|\mathcal{D}|>0$ is fixed and
$\|\bm{\Delta}\|_{F}=|\mathcal{D}|^{1/4}n^{-1/2}h$ when
$|\mathcal{D}|\to\infty$, where $\delta>0$ and $\|\cdot\|_{F}$ is the matrix
Frobenius norm. If the assumptions of Theorem 5 are met, then under $H_{a}$,
as $n,M\to\infty$,

$\beta(\bm{\theta}^{\circ},\bm{\Delta})\geq\begin{cases}\operatorname{\mathbb{P}}\Big{(}\|\bm{Z}+c_{l}\sqrt{n}\bm{\Delta}\|_{2}^{2}>\chi^{2}_{|\mathcal{D}|,1-\alpha}\Big{)}&\text{if
}|\mathcal{D}|>0\text{ is fixed};\\\
\operatorname{\mathbb{P}}\Big{(}Z>z_{1-\alpha}-{c_{l}\|\bm{\Delta}\|^{2}_{2}}/{\sqrt{2|\mathcal{D}|}}\Big{)}&\text{if
}|\mathcal{D}|\to\infty,\frac{|\mathcal{D}|\log|\mathcal{D}|}{n}\to
0,\end{cases}$

where $\bm{Z}\sim N(\bm{0},\bm{I}_{|\mathcal{D}|})$, $Z\sim N(0,1)$, and
$\chi^{2}_{|\mathcal{D}|,1-\alpha}$ and $z_{1-\alpha}$ are the $(1-\alpha)$th
quantile of distributions $\chi^{2}_{|\mathcal{D}|}$ and $N(0,1)$,
respectively. In particular,
$\lim_{\delta\to\infty}\lim_{n\to\infty}\beta(\bm{\theta}^{\circ},\bm{\Delta})=1$.

######  Proposition 8 (Local power of pathway test)

Assume that $H_{0}$ is nondegenerate and regular with
$|\mathcal{D}|=|\mathcal{H}|$. Let
$\min_{(k,j)\in\mathcal{H}}|U^{\circ}_{kj}+\Delta_{kj}|=n^{-1/2}\delta$ when
$|\mathcal{H}|>0$ is fixed and
$\min_{(k,j)\in\mathcal{H}}|U^{\circ}_{kj}+\Delta_{kj}|=n^{-1/2}\delta\sqrt{\log|\mathcal{H}|}$
when $|\mathcal{H}|\to\infty$. If the assumptions of Theorem 5 are met, then
under $H_{a}$, as $n,M\to\infty$,

$\beta(\bm{\theta}^{\circ},\bm{\Delta})\geq
1-\frac{|\mathcal{H}|}{\sqrt{2\pi}}\exp\Big{(}-\Big{(}\delta\sqrt{\log|\mathcal{H}|}/\max_{1\leq
j\leq p}\Omega_{jj}-\sqrt{\chi^{2}_{1,1-\alpha}}\Big{)}^{2}/2\Big{)},$

where $Z\sim N\Big{(}\delta^{2}/\max_{1\leq j\leq p}\Omega_{jj},1\Big{)}$ and
$\chi^{2}_{1,1-\alpha}$ is the $(1-\alpha)$th quantile of distribution
$\chi^{2}_{1}$. Then,
$\lim_{\delta\to\infty}\lim_{n\to\infty}\beta(\bm{\theta}^{\circ},\bm{\Delta})=1$.

The proofs of Propositions 7 and 8 are deferred to Appendix B.7 and B.8.

##  5 Simulations

This section investigates the operating characteristics of the proposed tests
and the peeling algorithm via simulations. In simulations, we consider two
setups for generating $\bm{U}\in\mathbb{R}^{p\times p}$, representing random
and hub DAGs, respectively.

  * •

Random graph. The upper off-diagonal entries $U_{kj}$; $k<j$ are sampled
independently from $\\{0,1\\}$ according to Bernoulli$(1/p)$, while other
entries are zero. This generates a random graph with a sparse neighborhood.

  * •

Hub graph. Set $U_{1,2j+1}=1$ and $U_{2,2j+2}=1$ for $j=1,\ldots,\lfloor
p/2\rfloor-2$, while other entries are zero. This generates a hub graph, where
nodes 1 and 2 are hub nodes with a dense neighborhood.

Moreover, we consider three setups for intervention matrix
$\bm{W}\in\mathbb{R}^{q\times p}$, representing different scenarios. Setups A
and B are designed for inference, whereas Setup C in Section 5.2 is designed
to compare with the method of Chen et al. (2018) for structure learning. Let
$\bm{W}=(\bm{A}^{\top},\bm{B}^{\top},\bm{0}^{\top})^{\top}$, where
$\bm{A},\bm{B}\in\mathbb{R}^{p\times p}$ and
$\bm{0}\in\mathbb{R}^{(q-2p)\times p}$.

  * •

Setup A. Set $A_{jj}=B_{jj}=B_{j,j+1}=1$; $j=1,\ldots,p-1$, $A_{pp}=1$, while
other entries of $\bm{A}$, $\bm{B}$ are zero. Then, $X_{1},\ldots,X_{p}$ are
instruments for $Y_{1},\ldots,Y_{p}$, respectively, $X_{p+1},\ldots,X_{2p-1}$
are invalid instruments with two targets, and $X_{2p},\ldots,X_{q}$ represent
inactive interventions.

  * •

Setup B. Set $A_{jj}=A_{j,j+1}=B_{jj}=B_{j,j+1}=1$; $j=1,\ldots,p-1$,
$A_{pp}=1$, while other entries of $\bm{A}$, $\bm{B}$ are zero. Here, the only
valid instrument is $X_{p}$ on $Y_{p}$, and the other intervention variables
either have two targets or are inactive. Importantly, the instrumental
sufficiency (Assumption 1C) is not met.

To generate $(\bm{X},\bm{Y})$ for each setup, we sample $\bm{X}\sim
N(\bm{0},\bm{\Sigma}_{X})$ with
$(\Sigma_{X})_{ll^{\prime}}=0.5^{|l-l^{\prime}|}$; $1\leq l,l^{\prime}\leq q$
and sample $\bm{Y}$ according to (1) with
$(\bm{U},\bm{W},\sigma_{1}^{2},\ldots,\sigma_{p}^{2})$, where
$\sigma_{1}^{2},\ldots,\sigma_{p}^{2}$ are set to be equally spaced from $0.5$
to $1$.

###  5.1 Inference

We compare three tests in empirical type-I errors and powers in simulated
examples, namely, the DP likelihood ratio test (DP-LR) in Algorithm 3, the
asymptotic likelihood ratio test (LR), and the oracle likelihood ratio test
(OLR). Here LR uses $\text{Lr}(\widehat{S},\widehat{\bm{\Sigma}})$, while OLR
uses $\text{Lr}(S^{\circ},\widehat{\bm{\Sigma}})$ assuming that the super-
graph $S^{\circ}$ were known in advance. The p-values of LR and OLR are
computed via Proposition 6. The implementation details of these tests are in
Appendix C.

For the empirical type-I error of a test, we compute the percentage of times
rejecting $H_{0}$ out of 500 simulations when $H_{0}$ is true. For the
empirical power of a test, we report the percentage of times rejecting $H_{0}$
out of 100 simulations when $H_{a}$ is true under alternative hypotheses
$H_{a}$.

  * •

Test of directed edges. For (2), we examine two different hypotheses:

(i) $H_{0}:U_{1,20}=0$ versus $H_{a}:U_{1,20}\neq 0$. In this case,
$|\mathcal{D}|=1$.

(ii) $H_{0}:\bm{U}_{\mathcal{H}}=\bm{0}$ versus
$H_{a}:\bm{U}_{\mathcal{H}}\neq\bm{0}$, where
$\mathcal{H}=\\{(k,20):k=1,\ldots,15\\}$. In this case, $|\mathcal{D}|=15$.

Moreover, five alternatives $H_{a}:U_{1,20}=0.1l$ and
$\bm{U}_{\mathcal{H}\setminus\\{(1,20)\\}}=\bm{0}$; $l=1,2,3,4,5$, are used
for the power analysis in (i) and (ii). The data are generated by modifying
$\bm{U}$ accordingly.

  * •

Test of directed pathways. We examine the test of a directed path $Y_{1}\to
Y_{5}\to Y_{10}\to Y_{15}\to Y_{20}$, namely
$\mathcal{H}=\\{(1,5),(5,10),(10,15),(15,20)\\}$ in (3). Since (3) is a test
of composite null hypothesis, the data are generated under a graph with
parameters $(\bm{U},\bm{W},\bm{\Sigma})$ satisfying $H_{0}$, where
$\bm{U}_{\mathcal{H}}=\bm{0}$. Five alternatives
$H_{a}:\bm{U}_{\mathcal{H}}=\bm{0.1}l$; $l=1,2,3,4,5$ are used for the power
analysis.

Figure 2: Empirical type-I errors and powers of tests of directed edges. The
black dotted line marks the nominal level of significance $\alpha=0.05$.
Figure 3: Empirical type-I errors and powers of tests of a directed pathway.
The black dotted line marks the nominal level of significance $\alpha=0.05$.

For testing directed edges, as displayed in Figure 2, DP-LR and LR perform
well compared to the ideal test OLR in Setup A with Assumption 1 satisfied. In
setup B with Assumption 1C not fulfilled, DP-LR appears to have control of
type-I error, whereas LR has an inflated empirical type-I error compared to
the nominal level $\alpha=0.05$. This discrepancy is attributed to the data
perturbation scheme accounting for the uncertainty of identifying $S$.
However, both suffer a loss of power compared to the oracle test OLR in this
setup. This observation suggests that without Assumption 1C, the peeling
algorithm tends to yield an estimate $\widehat{S}\supseteq S^{\circ}$, which
overestimates $S^{\circ}$, resulting in a power loss.

For testing directed pathways, as indicated in Figure 3, we observe similar
phenomena as in the previous directed edge tests. Of note, both LR and DP-LR
are capable in controlling type-I error of directed path tests.

In summary, DP-LR has a suitable control of type-I error when there are
invalid instruments and Assumption 1C is violated. Concerning the power, DP-LR
and LR are comparable in all scenarios and their powers tend to one as the
sample size $n$ or the signal strength of tested edges increases. Moreover,
DP-LR and LR perform nearly as well as the oracle test OLR when Assumption 1
is satisfied. These empirical findings agree with our theoretical results.

###  5.2 Structure learning

This subsection compares the peeling algorithm with the Two-Stage Penalized
Least Squares (2SPLS, Chen et al. (2018)) in terms of the structure learning
accuracy. For peeling, we consider Algorithm 2 with an additional step (9) for
structure learning of $\bm{U}$. For 2SPLS, we use the R package BigSEM.

Figure 4: Structural Hamming distances (SHDs) for the peeling algorithm and
2SPLS, where a smaller value of SHD indicates a better result.

2SPLS requires that all the intervention variables to be target-known
instruments in addition to Assumption 1C. Thus, we consider an additional
Setup C.

  * •

Setup C. Let $\bm{W}=(\bm{I}_{p\times p},\bm{0})^{\top}\in\mathbb{R}^{q\times
p}$. Then $X_{1},\ldots,X_{p}$ are valid instruments for $Y_{1},\ldots,Y_{p}$,
respectively, and other intervention variables are inactive.

For 2SPLS, we assign each active intervention variable to its most correlated
primary variable in Setups A-C. In Setup C, this assignment yields a correct
identification of valid instruments, meeting all the requirements of 2SPLS.

For each scenario, we compute the average structural Hamming distance
$\text{SHD}(\widehat{\bm{U}},\bm{U}^{\circ})=\sum_{k,j}|\operatorname{I}(\widehat{U}_{kj}\neq
0)-\operatorname{I}(U^{\circ}_{kj}\neq 0)|$ over 100 runs. As shown in Figure
4, the peeling algorithm outperforms 2SPLS, especially when there are invalid
instruments and Assumption 1C is violated.

###  5.3 Comparison of inference and structure learning

Figure 5: Rejection rates for $H_{0}:U_{1,20}=0$ versus
$H_{a}:U_{1,20}=1/\sqrt{n}$ by DP inference and structure learning. For
structure learning, $H_{0}$ is rejected if $\widehat{U}_{1,20}\neq 0$. The red
lines indicate the results of DP inference using the significance level
$\alpha=0.05$. The other colored lines display the results of structure
learning using different sparsity parameter values $\kappa=1,2,3,4,5$. The
simulation is repeated for 500 times and $\kappa=2$ is chosen by BIC in over
90% cases.

This subsection compares the proposed DP testing method against the proposed
structure learning method in (9) in terms of inferring the true graph
structure. To this end, we consider Setup A in Section 5.1 with $p=30,q=100$,
and the hypotheses

$H_{0}:U_{1,20}=0$ versus $H_{a}:U_{1,20}=1/\sqrt{n}$.

For DP inference, we use $\alpha=0.05$ and choose the tuning parameters by BIC
as in previous experiments; see Appendix C for details. For structure
learning, we reject the null hypothesis when $\widehat{U}_{1,20}\neq 0$.

As displayed in Figure 5, when the null hypothesis $H_{0}$ is true, the DP
testing method controls type-I error very close to the nominal level of 0.05,
whereas the type-I error of the structure learning varies greatly depending on
the tuning parameter selection. Under the alternative hypothesis $H_{a}$, the
DP inference enjoys high statistical power than structure learning methods
when $n\geq 200$. Interestingly, the power of structure learning diminishes as
$n$ increases. This observation is in agreement with our theoretical results
in Theorem 4, suggesting that consistent reconstruction requires the smallest
size of nonzero coefficients to be of order $\sqrt{\log(n)/n}$ with the tuning
parameter $\tau$ of the same order (fixing $p,q$). In this case, the edge
$U_{1,20}$ is of order $1/\sqrt{n}$, which is less likely to be reconstructed
as $n$ increases. In contrast, Proposition 7 indicates that a DP test has a
non-vanishing power when the hypothesized edges are of order $1/\sqrt{n}$.

Figure 5 demonstrates some important distinctions between inference and
structure learning. When different tuning parameters are used, the structure
learning results correspond to different points on an ROC curve. Although it
is asymptotically consistent when optimal tuning parameters are used,
structure learning lacks an uncertainty measure of graph structure
identification. As a result, it is nontrivial for structure learning methods
to trade-off the false discovery rate and detection power in practice. This
makes the interpretation of such results hard, especially when they heavily
rely on hyperparameters as in Figure 5. By comparsion, DP inference aims to
maximize statistical power while controlling type-I error at a given level,
offering a clear interpretation of its result. This observation agrees with
the discussions in the literature on variable selection and inference
(Wasserman and Roeder, 2009; Meinshausen and Bühlmann, 2010; Lockhart et al.,
2014; Candes et al., 2018) and it justifies the demand for inferential tools
for directed graphical models.

##  6 ADNI data analysis

This section applies the proposed tests to analyze an Alzheimer’s Disease
Neuroimaging Initiative (ADNI) dataset. In particular, we infer gene pathways
related to Alzheimer’s Disease (AD) to highlight some gene-gene interactions
differentiating patients with AD/cognitive impairments and healthy
individuals.

The raw data are available in the ADNI database (https://adni.loni.usc.edu),
including gene expression, whole-genome sequencing, and phenotypic data. After
cleaning and merging, we have a sample size of 712 subjects. From the KEGG
database (Kanehisa and Goto, 2000), we extract the AD reference pathway
(hsa05010, https://www.genome.jp/pathway/hsa05010), including 146 genes in the
ADNI data.

For data analysis, we first regress the gene expression levels on five
covariates – gender, handedness, education level, age, and intracranial
volume, and then use the residuals as gene expressions in the following
analysis. Next, we extract the genes with at least one SNP at a marginal
significance level below $10^{-3}$, yielding $p=63$ genes as primary
variables. For these genes, we further extract their marginally most
correlated two SNPs, resulting in $q=63\times 2=126$ SNPs as unspecified
intervention variables, in subsequent data analysis. All gene expression
levels are normalized.

The data contains individuals in four groups, namely, Alzheimer’s Disease
(AD), Early Mild Cognitive Impairment (EMCI), Late Mild Cognitive Impairment
(LMCI), and Cognitive Normal (CN). For our purpose, we treat 247 CN
individuals as controls while the remaining 465 individuals as cases (AD-MCI).
Then, we use the gene expressions and the SNPs to reconstruct the ancestral
relations and infer gene pathways for 465 AD-MCI and 247 CN control cases,
respectively.

(a) AD-MCI

(b) CN

(c) AD-MCI

(d) CN

Figure 6:  Display of the subnetworks associated with genes APP and CASP3. (a)
and (b): Solid/dashed arrows indicate significant/insignificant edges at
$\alpha=0.05$ after adjustment for multiplicity by the Bonferroni-Holm
correction. (c) and (d): Solid arrows indicate the reconstructed edges using
2SPLS (Chen et al., 2018).

(a) AD-MCI

(b) CN

Figure 7: The p-values of pathway tests (3) by the proposed tests for the AD-
MCI and CN groups, where p-values are adjusted for multiplicity by the
Bonferroni-Holm correction and solid/dashed arrows indicate
significant/insignificant pathways at $\alpha=0.05$.

In the literature, genes APP, CASP3, and PSEN1 are well-known to be associated
with AD, reported to play different roles in AD patients and healthy subjects
(Julia and Goate, 2017; Su et al., 2001; Kelleher and Shen, 2017). For this
dataset, we conduct hypothesis testing on edges and pathways related to genes
APP, CASP3, and PSEN1 in the KEGG AD reference (hsa05010) to evaluate the
proposed DP inference by checking if DP inference can discover the differences
that are reported in the biomedical literature. First, we consider testing
$H_{0}:U_{kj}=0$ versus $H_{a}:U_{kj}\neq 0$, for each edge $(k,j)$ as shown
in Figure 6 (a) and (b). Moreover, we consider two hypothesis tests of
pathways $H_{0}:U_{kj}=0$ for some $(k,j)\in\mathcal{P}_{\ell}$ versus
$H_{a}:U_{kj}\neq 0$ for all $(k,j)\in\mathcal{P}_{\ell}$; $\ell=1,2$, where
the two pathways are specified by
$\mathcal{P}_{1}=\\{\text{PSEN1}\to\text{CAPN1}\to\text{CDK5R1}\\}$, and
$\mathcal{P}_{2}=\\{\text{PSEN1}\to\text{CAPN2}\to\text{CDK5R1}\\}$. See
Figure 7. Of note, for clear visualization, Figure 6 (a), (b) and Figure 7
only display the edges related to hypothesis testing, whereas the ancestral
relations are reconstructed using $p=63$ genes and $q=126$ SNPs for AD-MCI and
CN groups separately.

In Figures 6-7, the significant results under the level $\alpha=0.05$ after
the Holm-Bonferroni adjustment for $2\times(9+2)=22$ tests are displayed. In
Figures 6, the edge test in (2) exhibits a strong evidence for the presence of
directed connectivity $\\{$APP $\to$ APBB1, APP $\to$ GSK3B, FADD $\to$
CASP3$\\}$ in the AD-MCI group, but no evidence in the CN group. Meanwhile,
this test suggests the presence of connections $\\{$TNFRSF1A $\to$ CASP3, FADD
$\to$ CASP8$\\}$ in the CN group but not so in the AD-MCI group. In both
groups, we identify directed connections $\\{$TNFRSF1A $\to$ FADD, TNFRSF1A
$\to$ CASP8$\\}$. In Figure 7, the pathway test (3) supports the presence of a
pathway PSEN1 $\to$ CAPN1 $\to$ CDK5R1 in the AD-MCI group with a p-value of
$0.044$ but not in the CN group with a p-value of $0.33$. The pathway PSEN1
$\to$ CAPN2 $\to$ CDK5R1 appears insignificant at $\alpha=0.05$ for both
groups. Also noted is that some of our discoveries agree with the literature
according to the AlzGene database (alzgene.org) and the AlzNet database
(https://mips.helmholtz-muenchen.de/AlzNet-DB). Specifically, GSK3B
differentiates AD patients from normal subjects; as shown in Figures 6, our
result indicates the presence of connection APP $\to$ GSK3B for the AD-MCI
group, but not for the CN group, the former of which is confirmed by Figure 1
of Vanleuven (2011). The connection APP $\to$ APBB1 also differs in AD-MCI and
CN groups, which appears consistent with Figure 3 of Bu (2009). Moreover, the
connection CAPN1 $\to$ CDK5R1, in the pathway PSEN1 $\to$ CAPN1 $\to$ CDK5R1
discovered in AD-MCI group, is found in the AlzNet database (interaction-ID
24614, https://mips.helmholtz-muenchen.de/AlzNet-DB/entry/show/1870). Finally,
as suggested by Figure 8, the normality assumption in (1) is adequate for both
groups.

By comparison, as shown in Figure 6 (c) and (d), gene APP in the reconstructed
networks by 2SPLS (Chen et al., 2018) is not connected with other genes,
indicating no regulatory relation of APP with other genes in the AD-MCI and CN
groups. However, as a well-known gene associated with AD, APP is reported to
play different roles in controlling the expressions of other genes for AD
patients and healthy people (Matsui et al., 2007; Julia and Goate, 2017). Our
results in Figure 6 (a) and (b) are congruous with the studies: the
connections of APP with other genes are different in our estimated networks
for AD-MCI and CN groups.

Figure 8: Normal quantile-quantile plots of studentized residuals of the AD-
MCI and CN groups.

In summary, our findings seem to agree with those in the literature (Julia and
Goate, 2017; Su et al., 2001; Kelleher and Shen, 2017), where the subnetworks
of genes APP, CASP3 in Figure 6 and PSEN1 in Figure 7 differentiate the AD-MCI
from the CN groups. Furthermore, the pathway PSEN1 $\to$ CAPN1 $\to$ CDK5R1 in
Figure 7 seems to differentiate these groups, which, however, requires
validation in biological experiments.

##  7 Summary

This article proposes structure learning and inference methods for a Gaussian
DAG with interventions, where the targets and strengths of interventions are
unknown. A likelihood ratio test is derived based on an estimated super-graph
formed by ancestral relations and candidate interventional relations. This
test accounts for the statistical uncertainty of the construction of the
super-graph based on a novel data perturbation scheme. Moreover, we develop a
peeling algorithm for the super-graph construction. The peeling algorithm
allows scalable computing and yields a consistent estimator. The numerical
studies justify our theory and demonstrate the utility of our methods.

The proposed methods can be extended to many practical situations beyond
biological applications with independent and identically distributed data. An
instance is to infer directed relations between multiple autoregressive time
series (Pamfil et al., 2020), where the lagged variables and covariates can
serve as interventions for each time series.

Acknowledgments

The authors would like to thank the action editor and anonymous referees for
helpful comments and suggestions. The research is supported by NSF grants
DMS-1712564, DMS-1952539 and NIH grants R01GM126002, R01HL116720, R01HL105397,
R01AG069895, R01AG065636, R01AG074858, and U01AG073079.

##  A Illustrative examples and discussions

###  A.1 Identifiability of model (1) and Assumption 1

The parameter space for model (1) is

$\\{(\bm{U},\bm{W},\bm{\Sigma}):\bm{U}\in\mathbb{R}^{p\times p}\text{
represents a DAG},\ \bm{W}\in\mathbb{R}^{p\times q},\
\bm{\Sigma}=\text{diag}(\sigma_{1}^{2},\ldots,\sigma^{2}_{p})\\}.$

As suggested by Proposition 1, Assumption 1 (1A-1C) suffices for
identification of every parameter value in the parameter space. Next, we show
by examples that if Assumption 1B or 1C is violated then model (1) is no
longer identifiable.

Before proceeding, we rewrite (1) as

$\bm{Y}=\bm{V}^{\top}\bm{X}+\bm{\varepsilon}_{V},\quad\bm{\varepsilon}_{V}=(\bm{I}-\bm{U}^{\top})^{-1}\bm{\varepsilon}\sim
N(\bm{0},\bm{\Omega}^{-1}),$

where $\bm{\Omega}=(\bm{I}-\bm{U})\bm{\Sigma}^{-1}(\bm{I}-\bm{U}^{\top})$ is a
precision matrix and $\bm{V}=\bm{W}(\bm{I}-\bm{U})^{-1}$.

######  Example 2 (Identifiability)

In model (1), consider two non-identifiable bivariate situations: (1) $p=2$
and $q=4$ and (2) $p=2$ and $q=2$.

  1. (1)

Model (1) is non-identifiable when Assumption 1B breaks down. Consider two
different models with different parameter values:

$\displaystyle\bm{\theta}:$ $\displaystyle
Y_{1}=X_{1}+X_{2}+X_{3}+\varepsilon_{1},\quad$ $\displaystyle
Y_{2}=Y_{1}-X_{2}+X_{3}+X_{4}+\varepsilon_{2},$ (19)
$\displaystyle\widetilde{\bm{\theta}}:$ $\displaystyle
Y_{1}=0.5Y_{2}+0.5X_{1}+X_{2}-0.5X_{4}+\widetilde{\varepsilon}_{1},\quad$
$\displaystyle Y_{2}=X_{1}+2X_{3}+X_{4}+\widetilde{\varepsilon}_{2},$ (20)

where $\varepsilon_{1},\varepsilon_{2}\sim N(0,1)$ are independent, and
$\widetilde{\varepsilon}_{1}\sim N(0,0.5)$, $\widetilde{\varepsilon}_{2}\sim
N(0,2)$ are independent. As depicted in Figure 9, (19) satisfies Assumption
1C. However, Assumption 1B is violated given that
$\operatorname{Cov}(Y_{2},X_{2}\mid\bm{X}_{\\{1,3,4\\}})=0$ and $X_{2}$ is an
intervention variable of $Y_{2}$. This is because the direct interventional
effect of $X_{2}$ on $Y_{2}$ are canceled out by its indirect interventional
effect through $Y_{1}$. Similarly, (20) satisfies Assumption 1C but
$\operatorname{Cov}(Y_{1},X_{4}\mid\bm{X}_{\\{1,2,3\\}})=0$ violating
Assumption 1B. In this case, it can be verified that $\bm{\theta}$ and
$\widetilde{\bm{\theta}}$ correspond to the same distribution
$\operatorname{\mathbb{P}}(\bm{Y}\mid\bm{X})$, because they share the same
$(\bm{V},\bm{\Omega})$ even with different values of
$(\bm{U},\bm{W},\bm{\Sigma})$. Hence, it is impossible to infer the directed
relation between $Y_{1}$ and $Y_{2}$.

  2. (2)

Model (1) is non-identifiable when Assumption 1C breaks down. Consider two
different models with different parameter values:

$\displaystyle\bm{\theta}:$ $\displaystyle
Y_{1}=X_{1}+X_{2}+\varepsilon_{1},\quad$ $\displaystyle
Y_{2}=Y_{1}+X_{2}+\varepsilon_{2},$ (21)
$\displaystyle\widetilde{\bm{\theta}}:$ $\displaystyle
Y_{1}=0.5Y_{2}+0.5X_{1}+\widetilde{\varepsilon}_{1},\quad$ $\displaystyle
Y_{2}=X_{1}+2X_{2}+\widetilde{\varepsilon}_{2},$ (22)

where $\varepsilon_{1},\varepsilon_{2}\sim N(0,1)$ are independent, and
$\widetilde{\varepsilon}_{1}\sim N(0,0.5)$, $\widetilde{\varepsilon}_{2}\sim
N(0,2)$ are independent. Note that (21) and (22) satisfy Assumption 1B. In
(21), $Y_{2}$ does not have any instrumental intervention although it has an
invalid instrument $X_{2}$. Similarly, in (22), neither does $Y_{1}$ have any
instrumental intervention while having an invalid instrument $X_{1}$. As in
the previous case, $\bm{\theta}$ and $\widetilde{\bm{\theta}}$ yield the same
distribution $\operatorname{\mathbb{P}}(\bm{Y}\mid\bm{X})$ because they share
the same $(\bm{V},\bm{\Omega})$ even with different values of
$(\bm{U},\bm{W},\bm{\Sigma})$. In this case, it is impossible to infer the
directed relation between $Y_{1}$ and $Y_{2}$.

Figure 9: (a) Display of DAG defined by (19). (b) Display of DAG defined by
(20). Figure 10: (a) Display of DAG defined by (21). (b) Display of DAG
defined by (22).

###  A.2 Illustration of Algorithm 2

We now illustrate Algorithm 2 by Example 3.

###### Example 3

Consider model (1) with $p=q=5$,

$\begin{aligned} Y_{1}&=X_{1}+\varepsilon_{1},\\\
Y_{3}&=0.5Y_{2}+X_{5}+\varepsilon_{3},\\\ \end{aligned}\qquad\begin{aligned}
Y_{2}&=0.5Y_{1}+X_{3}+\varepsilon_{2},\\\
Y_{4}&=0.5Y_{3}-0.1Y_{1}+X_{2}+\varepsilon_{4},\\\
\end{aligned}\qquad\begin{aligned} Y_{5}&=X_{4}+\varepsilon_{5},\\\
&\end{aligned}$ (23)

where $\varepsilon_{1},\ldots,\varepsilon_{5}\sim N(0,1)$ are independent.
Then (23) defines a DAG as displayed in Figure 1(a). For illustration, we
generate a random sample of size $n=40$ and compute $\widehat{\bm{V}}$ by
Algorithm 1. In particular,

$\bm{V}^{\circ}=\begin{pmatrix}1&0.5&0.25&0.025&0\\\ 0&0&0&1&0\\\
0&1&0.5&0.25&0\\\ 0&0&0&0&1\\\
0&0&1&0.5&0\end{pmatrix},\quad\widehat{\bm{V}}=\begin{pmatrix}0.92&0.48&0.27&0&0\\\
0&0&0&1.08&0\\\ 0&1.03&0.52&0.21&0\\\ 0&0&0&0&1.06\\\
0&0&0.98&0.55&0\end{pmatrix}.\quad$

Algorithm 2 proceeds as follows.

  * •

$h=0$: We have the set $\\{X_{2},X_{4}\\}$ of instruments on leaves.

    * –

$X_{2}$ is identified as an instrument of leaf node $Y_{4}$ ($X_{2}\to Y_{4}$)
as $\widehat{V}_{24}\neq 0$ is the only nonzero in the row 2 with the smallest
(positive) row $\ell_{0}$-norm.

    * –

$X_{4}$ is identified as an instrument of leaf node $Y_{5}$ ($X_{4}\to Y_{5}$)
as $\widehat{V}_{45}\neq 0$ is the only nonzero in row 4 with the smallest
(positive) row $\ell_{0}$-norm.

Then $h_{G}(4)=h_{G}(5)=0$ and $\\{Y_{4},Y_{5}\\}$ are removed.

  * •

$h=1$: We have the set $\\{X_{5}\\}$ of instrument on leaf.

    * –

$X_{5}$ is identified as an instrument of a leaf node $Y_{3}$ ($X_{5}\to
Y_{3}$) in the subgraph for $Y_{1},Y_{2},Y_{3}$ given that
$\widehat{V}_{53}\neq 0$ is the only nonzero element in the row with the
smallest (positive) row $\ell_{0}$-norm of the submatrix for
$Y_{1},Y_{2},Y_{3}$. Since $\widehat{V}_{54}\neq 0$ and $h_{G}(4)=0$, we have
$Y_{3}\to Y_{4}$ by Proposition 3.

Then $h_{G}(3)=1$ and $\\{Y_{3}\\}$ is removed.

  * •

$h=2$: We have the set $\\{X_{3}\\}$ of instrument on leaf.

    * –

$X_{3}$ is identified as an instrument of a leaf node $Y_{2}$ ($X_{3}\to
Y_{2}$) similarly in the subgraph for $Y_{1},Y_{2}$ given that
$\widehat{V}_{32}\neq 0$ is the largest nonzero element in its row of the
submatrix. Since $\widehat{V}_{33}\neq 0$ and $h_{G}(3)=1$, we have $Y_{2}\to
Y_{3}$.

Then $h_{G}(2)=2$ and $\\{Y_{2}\\}$ is removed.

  * •

$h=3$: We have the set $\\{X_{1}\\}$ of instrument on leaf.

    * –

$X_{1}$ is an instrument of $Y_{1}$ ($X_{1}\to Y_{1}$). Since
$\widehat{V}_{12}\neq 0$ and $h_{G}(2)=2$, we have $Y_{1}\to Y_{2}$.

Then $h_{G}(1)=3$, $\\{Y_{1}\\}$ is removed, and the peeling process is
terminated.

Finally, Step 6 of Algorithm 2 identifies the ancestral relations

$\widehat{\mathcal{A}}=\\{(1,2),(2,3),(3,4),(1,3),(1,4),(2,4)\\}$

and the candidate interventional relations

$\widehat{\mathcal{I}}=\\{(1,1),(1,2),(1,3),(2,4),(3,2),(3,3),(3,4),(4,5),(5,3),(5,4)\\}.$

In Example 3, $\widehat{\bm{V}}\neq\bm{V}^{\circ}$, suggesting that the
selection consistency of matrix $\bm{V}^{\circ}$ is unnecessary for Algorithm
2 to correctly reconstruct ancestral relations $\mathcal{A}^{\circ}$; see also
Section A.3.

###  A.3 Relaxation of Assumption 4

Assumption 4 in Theorem 4 leads to consistent identification for
$\bm{V}^{\circ}$. Now, we discuss when Algorithm 2 correctly reconstruct
ancestral relations $\mathcal{A}^{\circ}$ without requiring Assumption 4.

###### Assumption 6

There exists $(\kappa_{j}^{*},\tau^{*}_{j})$ such that
$\kappa_{j}^{*}=|\\{l:|V^{\circ}_{lj}|\geq\tau_{j}^{*}\\}|$ and

$\begin{split}&\min\Big{(}\min_{l:X_{l}\to
Y_{j}}|V^{\circ}_{lj}|,\min_{\begin{subarray}{c}l:X_{l}\to Y_{k}\to Y_{j}\\\
\text{ with }h(k)=h(j)+1\end{subarray}}|V^{\circ}_{lj}|\Big{)}\\\
&\geq\tau^{*}_{j}\geq
100\frac{c_{2}}{c_{1}}\sqrt{\Omega_{jj}^{-1}(\kappa^{\circ}_{j}-\kappa^{*}_{j}+1)\Big{(}\frac{\log(q)}{n}+\frac{\log(n)}{n}\Big{)}};\quad
1\leq j\leq p.\end{split}$ (24)

Assumption 6 requires a subset of nonzero entries in $\bm{V}^{\circ}$ to
exceed a certain signal strength. These signals enable us to reconstruct the
topological layers and $\mathcal{A}^{\circ}$. Note that Assumption 6 reduces
to Assumption 4 when $\kappa^{*}_{j}=\kappa^{\circ}_{j}$.

###### Theorem 9

Under Assumptions 1-3 and 6 with constants $c_{1}<6c_{2}$, if the machine
precision $\text{tol}\ll 1/n$ is negligible and tuning parameters of Algorithm
1 satisfy:

  1. (1)

$\gamma_{j}\leq c_{1}/6\sqrt{\kappa^{\circ}_{j}-\kappa^{*}_{j}+1}$,

  2. (2)

$\sqrt{32c_{2}^{2}\Omega_{jj}^{-1}(\log(q)+\log(n))/n}/\gamma_{j}\leq\tau_{j}\leq
0.4\tau^{*}_{j}$,

  3. (3)

$\kappa_{j}=\kappa_{j}^{\circ}$,

for $j=1,\ldots,p$, then Algorithm 1 terminates in at most
$1+\lceil\log(\kappa_{\max}^{\circ})/\log 4\rceil$ iterations almost surely.
Moreover, $\widehat{\mathcal{A}}=\mathcal{A}^{\circ}$,
$\widehat{\mathcal{C}}\supseteq\\{(l,j):|V_{lj}^{\circ}|\geq\tau^{*}_{j}\\}\supseteq\mathcal{I}^{\circ}$,
and $\widehat{\mathcal{D}}=\mathcal{D}^{\circ}$ almost surely as $n\to\infty$.

The proof of Theorem 9 is given in Appendix B.9.

###  A.4 Comparison of strong faithfulness and Assumption 4 (or 6)

In the literature, a faithfulness condition is usually assumed for
identifiability up to Markov equivalence classes (Spirtes et al., 2000). For
discussion, we formally introduce the concepts of faithfulness and strong
faithfulness.

Consider a DAG $G$ with node variables $(Z_{1},\ldots,Z_{p+q})^{\top}$. Nodes
$Z_{i}$ and $Z_{j}$ are adjacent if $Z_{i}\to Z_{j}$ or $Z_{j}\to Z_{i}$. A
_path_ (_undirected_) between $Z_{i}$ and $Z_{j}$ in $G$ is a sequence of
distinct nodes $(Z_{i},\ldots,Z_{j})$ such that all pairs of successive nodes
in the sequence are adjacent. A nonendpoint node $Z_{k}$ on a path
$(Z_{i},\ldots,Z_{k-1},Z_{k},Z_{k+1},\ldots,Z_{j})$ is a _collider_ if
$Z_{k-1}\to Z_{k}\leftarrow Z_{k+1}$. Otherwise it is a _noncollider_. Let
$A\subseteq\\{1,\ldots,p+q\\}$, where $A$ does not contain $i$ and $j$. Then
$\bm{Z}_{A}$ _blocks_ a path $(Z_{i},\ldots,Z_{j})$ if at least one of the
following holds: (i) the path contains a noncollider that is in $\bm{Z}_{A}$,
or (ii) the path contains a collider that is not in $\bm{Z}_{A}$ and has no
descendant in $\bm{Z}_{A}$. A node $Z_{i}$ is d-separated from $Z_{j}$ given
$\bm{Z}_{A}$ if $\bm{Z}_{A}$ block every path between $Z_{i}$ and $Z_{j}$;
$i\neq j$ (Pearl, 2000).

According to Uhler et al. (2013), a multivariate Gaussian distribution of
$(Z_{1},\ldots,Z_{p+q})^{\top}$ is said to be $\varsigma$-strong faithful to a
DAG with node set $\mathcal{V}=\\{1,\ldots,p+q\\}$ if

$\min_{A\subseteq\mathcal{V}\backslash\\{i,j\\}}\Big{\\{}|{\operatorname{Corr}}(Z_{i},Z_{j}\mid\bm{Z}_{A})|:Z_{i}\text{
is not d-separated from }Z_{j}\text{ given }\bm{Z}_{A}\Big{\\}}>\varsigma,$
(25)

for $1\leq i\neq j\leq p+q$, where $\varsigma\in[0,1)$,
${\operatorname{Corr}}$ denotes the correlation. When $\varsigma=0$, (25) is
equivalent to faithfulness. For consistent structure learning (up to Markov
equivalence classes), it often requires that
$\varsigma\gtrsim\sqrt{s_{0}\log(p+q)/n}$, where $s_{0}$ is a sparsity
measure; see Uhler et al. (2013) for a survey. For a pair $(i,j)$, the number
of possible sets for $A$ is $2^{(p+q-2)}$. If $Z_{i}\to Z_{j}$, then
$\text{Corr}(Z_{i},Z_{j}\mid\bm{Z}_{A})\neq 0$ for any $A$. Therefore, for
this $(i,j)$ pair alone, (25) could require exponentially many conditions.
Indeed, (25) is very restrictive in high-dimensional situation (Uhler et al.,
2013).

By comparison, Algorithm 2 yields consistent structure learning based on
Assumption 4 or 6 instead of strong faithfulness. In some sense, Assumption 4
or 6 requires the sufficient signal strength that is analogous to the
condition for consistent feature selection (Shen et al., 2012). This
assumption may be thought of as an alternative to strong faithfulness. As
illustrated in Example 4, Assumption 4 or 6 is less stringent than strong
faithfulness.

######  Example 4 (Faithfulness)

For simplicity, assume $\bm{X}\sim N(\bm{0},\bm{I})$. Consider model (1) with
$p=q=3$,

$\displaystyle Y_{1}=W_{11}X_{1}+\varepsilon_{1},\quad
Y_{2}=U_{12}Y_{1}+W_{22}X_{2}+\varepsilon_{2},\quad
Y_{3}=U_{13}Y_{1}+U_{23}Y_{2}+W_{33}X_{3}+\varepsilon_{3},$

where $\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}\sim N(0,1)$ are
independent and $U_{12},U_{13},U_{23},W_{11},W_{22},W_{33}\neq 0$. Let
$\bm{Z}=(Y_{1},Y_{2},Y_{3},X_{1},X_{2},X_{3})^{\top}$. Since directed
relations among $\bm{X}$ are not of interest, (25) becomes

$\begin{split}\min_{\begin{subarray}{c}\bm{Z}_{A}=(\bm{Y}_{A_{1}},\bm{X}_{A_{2}}):\\\
A_{1}\subseteq\\{i,j\\}^{c},A_{2}\end{subarray}}&\Big{\\{}|{\operatorname{Corr}}(Y_{i},Y_{j}\mid\bm{Z}_{A})|:Y_{i},Y_{j}\text{
are not d-separated given }\bm{Z}_{A}\Big{\\}}>\varsigma,\\\
\min_{\begin{subarray}{c}\bm{Z}_{A}=(\bm{Y}_{A_{1}},\bm{X}_{A_{2}}):\\\
A_{1}\subseteq\\{j\\}^{c},A_{2}\subseteq\\{l\\}^{c}\end{subarray}}&\Big{\\{}|{\operatorname{Corr}}(Y_{j},X_{l}\mid\bm{Z}_{A})|:Y_{j},X_{l}\text{
are not d-separated given }\bm{Z}_{A}\Big{\\}}>\varsigma,\end{split}$ (26)

for each pair $(i,j)$ with $i\neq j$ and each pair $(l,j)$. Then strong-
faithfulness in (26) assumes 152 conditions for the correlations. By
comparison, (17) requires the absolute values of
$V_{11},V_{12},V_{13},V_{22},V_{23},V_{33}\gtrsim\sqrt{\log(q)/n}$, which in
turn requires the minimal absolute value of six correlations
$\gtrsim\sqrt{\log(q)/n}$,
(i) ${\operatorname{Corr}}(Y_{1},X_{1}\mid X_{2},X_{3})$, (ii)
${\operatorname{Corr}}(Y_{2},X_{1}\mid X_{2},X_{3})$, (iii)
${\operatorname{Corr}}(Y_{3},X_{1}\mid X_{2},X_{3})$,
(iv) ${\operatorname{Corr}}(Y_{2},X_{2}\mid X_{1},X_{3})$, (v)
${\operatorname{Corr}}(Y_{3},X_{2}\mid X_{1},X_{3})$, (vi)
${\operatorname{Corr}}(Y_{3},X_{3}\mid X_{1},X_{2})$.
Importantly, (i)-(vi) are required in (26), suggesting that strong-
faithfulness is more stringent than (17).

###  A.5 Irregular hypothesis

Assume, without loss of generality, that
$\widehat{\mathcal{D}}=\mathcal{D}^{\circ}$ and
$\widehat{\mathcal{A}}=\mathcal{A}^{\circ}$ are correctly reconstructed in the
following discussion.

  * •

For testing of directed edges (2), suppose $H_{0}$ is irregular, namely,
$\mathcal{D}^{\circ}\cup\mathcal{E}^{\circ}$ contains a directed cycle. This
implies that a directed cycle exists in
$\widehat{\mathcal{D}}\cup\widehat{\mathcal{A}}$. In this situation, we
decompose $H_{0}$ into sub-hypotheses $H_{0}^{(1)},\ldots,H_{0}^{(\nu)}$, each
of which is regular. Then testing $H_{0}$ is equivalent to multiple testing
for $H_{0}^{(1)},\ldots,H_{0}^{(\nu)}$. For instance, in Example 1,
$H_{0}:U_{45}=U_{53}=0$ is irregular, and $H_{0}$ can be decomposed into
$H_{0}^{(1)}:U_{45}=0$ and $H_{0}^{(2)}:U_{53}=0$.

  * •

For testing of directed pathways (3), if $H_{0}$ is irregular, then
$\widehat{\mathcal{D}}\cup\widehat{\mathcal{A}}$ has a directed cycle. The
p-value is defined to be one in this situation since no evidence supports the
presence of the pathway.

##  B Technical proofs

###  B.1 Proof of Proposition 1

Suppose that $\bm{\theta}=(\bm{U},\bm{W},\bm{\Sigma})$ and
$\widetilde{\bm{\theta}}=(\widetilde{\bm{U}},\widetilde{\bm{W}},\widetilde{\bm{\Sigma}})$
render the same distribution of $(\bm{X},\bm{Y})$. We will prove that
$\bm{\theta}=\widetilde{\bm{\theta}}$.

Denote by $G(\bm{\theta})$ and $G(\widetilde{\bm{\theta}})$ the DAGs
corresponding to $\bm{\theta}$ and $\widetilde{\bm{\theta}}$, respectively.
First, consider $G(\bm{\theta})$. Without loss of generality, assume $Y_{1}$
is a leaf node. By Assumption 1C, there exists an instrumental intervention
with respect to $G(\bm{\theta})$, say $X_{1}$. Then,

$\displaystyle\operatorname{Cov}(Y_{j},X_{1}\mid\bm{X}_{\\{2,\ldots,q\\}})=0,$
$\displaystyle\quad j=2,\ldots,p,$ (27)
$\displaystyle\operatorname{Cov}(Y_{1},X_{1}\mid\bm{Y}_{A},\bm{X}_{\\{2,\ldots,q\\}})\neq
0,$ $\displaystyle\quad\text{for any }A\subseteq\\{2,\ldots,p\\}.$ (28)

By the local Markov property (Spirtes et al., 2000), (28) implies that
$X_{1}\to Y_{1}$ in $G(\widetilde{\bm{\theta}})$. Suppose $Y_{1}$ is not a
leaf node in $G(\widetilde{\bm{\theta}})$. Assume, without loss of generality,
that $Y_{1}\to Y_{2}$ with
$h_{G(\widetilde{\bm{\theta}})}(1)=h_{G(\widetilde{\bm{\theta}})}(2)+1$. Then
$\operatorname{Cov}(Y_{2},X_{1}\mid\bm{X}_{\\{2,\ldots,q\\}})=0$ but $X_{1}\to
Y_{1}$ and $Y_{1}\to Y_{2}$ with
$h_{G(\widetilde{\bm{\theta}})}(1)=h_{G(\widetilde{\bm{\theta}})}(2)+1$, which
contradicts to Assumption 1B. This implies that if $Y_{1}$ is a leaf node in
$G(\bm{\theta})$ then it must be a leaf node in $G(\widetilde{\bm{\theta}})$.
In both $G(\bm{\theta})$ and $G(\widetilde{\bm{\theta}})$, the parents and
interventions of $Y_{1}$ can be identified by

$\begin{split}\operatorname{\mathbb{E}}(Y_{1}\mid\bm{Y}_{\\{2,\ldots,p\\}},\bm{X})&=\operatorname{\mathbb{E}}(Y_{1}\mid\bm{Y}_{\textnormal{{pa}}_{G(\bm{\theta})}(1)},\bm{X})=\operatorname{\mathbb{E}}(Y_{1}\mid\bm{Y}_{\textnormal{{pa}}_{G(\bm{\theta})}(1)},\bm{X}_{\textnormal{{in}}_{G(\bm{\theta})}(1)}),\\\
\operatorname{\mathbb{E}}(Y_{1}\mid\bm{Y}_{\\{2,\ldots,p\\}},\bm{X})&=\operatorname{\mathbb{E}}(Y_{1}\mid\bm{Y}_{\textnormal{{pa}}_{G(\widetilde{\bm{\theta}})}(1)},\bm{X})=\operatorname{\mathbb{E}}(Y_{1}\mid\bm{Y}_{\textnormal{{pa}}_{G(\widetilde{\bm{\theta}})}(1)},\bm{X}_{\textnormal{{in}}_{G(\widetilde{\bm{\theta}})}(1)}).\end{split}$

Consequently, $Y_{1}$ has the same parents and interventions in
$G(\bm{\theta})$ and $G(\widetilde{\bm{\theta}})$.

The forgoing argument is applied to other nodes sequentially. First, we remove
$Y_{1}$ with any directed edges to $Y_{1}$, which does not alter the joint
distribution of $(\bm{Y}_{\\{2,\ldots,p\\}},\bm{X})$ and the subgraph of nodes
$Y_{2},\ldots,Y_{p}$. By induction, we remove the leafs in $G(\bm{\theta})$
until it is empty, leading to $G(\bm{\theta})=G(\widetilde{\bm{\theta}})$.
Finally, $\bm{\theta}=\widetilde{\bm{\theta}}$ because they have the same
locations for nonzero elements and these model parameters (or regression
coefficients) are uniquely determined under Assumption 1A (Shojaie and
Michailidis, 2010). This completes the proof.

###  B.2 Proof of Proposition 2

It follows from Assumption 1B that $\mathcal{C}=\\{(l,j):V_{lj}\neq
0\\}\supseteq\mathcal{I}$.

#### Proof of (A)

Note that the maximal length of a path in a DAG of $p$ nodes is at most $p-1$.
Then it can be verified that $\bm{U}$ is nilpotent in that
$\bm{U}^{p}=\bm{0}$. An application of the matrix series expansion yields that
$(\bm{I}-\bm{U})^{-1}=\bm{I}+\bm{U}+\ldots+\bm{U}^{p-1}$. Using the fact that
$\bm{V}=\bm{W}(\bm{I}-\bm{U})^{-1}$ from (6), we have, for any $1\leq l,j\leq
p$,

$V_{lj}=\sum_{k=1}^{p}W_{lk}\Big{(}I_{kj}+U_{kj}+\ldots+(\bm{U}^{p-1})_{kj}\Big{)},$

where $U_{kj}$ is the $(k,j)$th entry of $\bm{U}$. If $V_{lj}\neq 0$, then
there exists $k$ such that $W_{lk}\neq 0$ and $(\bm{U}^{r})_{kj}\neq 0$ for
some $0\leq r\leq p-1$. If $r=0$, then we must have $k=j$, and $X_{l}\to
Y_{j}$. If $r>0$, then $X_{l}\to Y_{k}$ and $Y_{k}$ is an ancestor of $Y_{j}$.

#### Proof of (B)

First, for any leaf node variable $Y_{j}$, by Assumption 1, there exists an
instrument $X_{l}\to Y_{j}$. If $V_{lj^{\prime}}\neq 0$ for some
$j^{\prime}\neq j$, then $Y_{j}$ must be an ancestor of $Y_{j^{\prime}}$,
which contradicts the fact that $Y_{j}$ is a leaf node variable.

Conversely, suppose that $V_{lj}\neq 0$ and $V_{lj^{\prime}}=0$ for
$j^{\prime}\neq j$. If $Y_{j}$ is not a leaf node variable, then there exists
a variable $Y_{j^{\prime}}$ such that $Y_{j}\to Y_{j^{\prime}}$ with
$h_{G}(j)=h_{G}(j^{\prime})+1$, that is $U_{jj^{\prime}}\neq 0$ and
$(\bm{U}^{r})_{jj^{\prime}}=0$ for $r>1$. Then
$V_{lj^{\prime}}=W_{lj}U_{jj^{\prime}}\neq 0$, a contradiction.

###  B.3 Proof of Proposition 3

Suppose $Y_{k}\to Y_{j}$ and $h_{G}(j)=h_{G}(k)-1$. Let $X_{l}$ be an
instrument of $Y_{k}$ in $G^{\prime}$. Then there are two cases: (1) $X_{l}$
intervenes on $Y_{k}$ but does not intervene on $Y_{j}$, namely $X_{l}\to
Y_{k}$ but $X_{l}\not\to Y_{j}$; (2) $X_{l}$ intervenes on $Y_{k}$ and $Y_{j}$
simultaneously, namely $X_{l}\to Y_{k}$ and $X_{l}\to Y_{j}$. For (1),
$V_{lj}=W_{lk}U_{kj}\neq 0$. For (2), Assumption 1B implies that $V_{lj}\neq
0$. This holds for every instrument of $Y_{k}$ in $G^{\prime}$, and the
desired result follows.

Conversely, suppose for each instrument $X_{l}$ of $Y_{k}$ in $G^{\prime}$, we
have $V_{lj}\neq 0$. Let $X_{l^{\prime}}$ be an instrument of $Y_{k}$ in $G$,
which is also an instrument in $G^{\prime}$. Then $0\neq
V_{l^{\prime}j}=W_{l^{\prime}k}U_{kj}$, which implies $U_{kj}\neq 0$. This
completes the proof.

###  B.4 Proof of Theorem 4

The proof proceeds in two steps: (A) Show that $\\{l:\widehat{V}_{lj}\neq
0\\}=\\{l:V_{lj}^{\circ}\neq 0\\}$ almost surely for $j=1,\ldots,p$; and (B)
Show that $\widehat{S}=S^{\circ}$ if $\widehat{\bm{V}}$ satisfies the property
in (A). As a result, $\widehat{\mathcal{D}}=\mathcal{D}^{\circ}$ follows
immediately.

#### Proof of (A)

Let $A_{j}^{\circ}=\\{l:V^{\circ}_{lj}\neq 0\\}$ and
$A^{[t]}_{j}=\\{l:|\widetilde{V}^{[t]}_{lj}|\geq\tau_{j}\\}$ be the estimated
support of penalized solution at the $t$-th iteration of Algorithm 1. For the
penalized solution, define the false negative set
$\text{FN}^{[t]}_{j}=A^{\circ}_{j}\setminus A_{j}^{[t]}$ and the false
positive set $\text{FP}^{[t]}_{j}=A_{j}^{[t]}\setminus A_{j}^{\circ}$; $t\geq
0$. Consider a “good” event
$\mathscr{E}_{j}=\\{\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}\leq
0.5\gamma_{j}\tau_{j}\\}\cap\\{\|\widehat{\bm{V}}^{\circ}_{\cdot
j}-\bm{V}^{\circ}_{\cdot j}\|_{\infty}\leq 0.5\tau_{j}\\}$, where
$\widehat{\bm{\xi}}_{j}=\bm{\mathsf{Y}}_{j}-\bm{\mathsf{X}}\widehat{\bm{V}}^{\circ}_{\cdot
j}$ is the residual of the least squares estimate
$\widehat{\bm{V}}^{\circ}_{\cdot j}$ such that
$A_{j}^{\circ}=\\{l:\widehat{V}^{\circ}_{lj}\neq 0\\}$; $j=1,\ldots,p$. We
shall show that $\text{FN}^{[t]}_{j}$ and $\text{FP}^{[t]}_{j}$ are eventually
empty sets on event $\mathscr{E}_{j}$ which has a probability tending to one.

First, we show that if $|A_{j}^{\circ}\cup A^{[t-1]}_{j}|\leq
2\kappa_{\max}^{\circ}$ on $\mathscr{E}_{j}$ for $t\geq 1$, then
$|A_{j}^{\circ}\cup A^{[t]}_{j}|\leq 2\kappa^{\circ}_{\max}$, to be used in
Assumption 2. Now suppose $|A_{j}^{\circ}\cup A^{[t-1]}_{j}|\leq
2\kappa^{\circ}_{\max}$ on $\mathscr{E}_{j}$ for $t\geq 1$. By the optimality
condition of (8),

$\begin{split}\Big{(}\widehat{\bm{V}}^{o}_{\cdot
j}-\widetilde{\bm{V}}^{[t]}_{\cdot
j}\Big{)}^{\top}\Big{(}-\bm{\mathsf{X}}^{\top}(\bm{\mathsf{Y}}_{j}-\bm{\mathsf{X}}\widetilde{\bm{V}}^{[t]}_{\cdot
j})/n+\gamma_{j}\tau_{j}\nabla\|\widetilde{\bm{V}}_{(A^{[t-1]}_{j})^{c},j}^{[t]}\|_{1}\Big{)}\geq
0,\end{split}$

where $\widetilde{\bm{V}}^{[t]}$ is defined in (8). Plugging
$\widehat{\bm{\xi}}_{j}=\bm{\mathsf{Y}}_{j}-\bm{\mathsf{X}}\widehat{\bm{V}}^{\circ}_{\cdot
j}$ into the inequality and rearranging it, we have
$\|\bm{\mathsf{X}}(\widetilde{\bm{V}}_{\cdot j}^{[t]}-\widehat{\bm{V}}_{\cdot
j}^{\circ})\|_{2}^{2}/n$ is no greater than

$\begin{split}&\Big{(}\widetilde{\bm{V}}_{\cdot
j}^{[t]}-\widehat{\bm{V}}_{\cdot
j}^{\circ}\Big{)}^{\top}\Big{(}\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n-\gamma_{j}\tau_{j}\nabla\|\widetilde{\bm{V}}^{[t]}_{(A_{j}^{[t-1]})^{c},j}\|_{1}\Big{)}\\\
&=\Big{(}\widetilde{\bm{V}}_{A_{j}^{\circ}\setminus
A_{j}^{[t-1]},j}^{[t]}-\widehat{\bm{V}}_{A_{j}^{\circ}\setminus
A_{j}^{[t-1]},j}^{\circ}\Big{)}^{\top}\Big{(}\bm{\mathsf{X}}^{\top}_{A_{j}^{\circ}\setminus
A^{[t-1]}_{j}}\widehat{\bm{\xi}}_{j}/n-\gamma_{j}\tau_{j}\nabla\|\widetilde{\bm{V}}^{[t]}_{A_{j}^{\circ}\setminus
A^{[t-1]}_{j},j}\|_{1}\Big{)}\\\
&\quad+\Big{(}\widetilde{\bm{V}}_{(A_{j}^{\circ}\cup
A^{[t-1]}_{j})^{c},j}^{[t]}-\widehat{\bm{V}}_{(A_{j}^{\circ}\cup
A^{[t-1]}_{j})^{c},j}^{\circ}\Big{)}^{\top}\Big{(}\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}/n-\gamma_{j}\tau_{j}\nabla\|\widetilde{\bm{V}}^{[t]}_{(A_{j}^{\circ}\cup
A^{[t-1]}_{j})^{c},j}\|_{1}\Big{)}\\\
&\quad+\Big{(}\widetilde{\bm{V}}_{A^{[t-1]}_{j}\setminus
A_{j}^{\circ},j}^{[t]}-\widehat{\bm{V}}_{A^{[t-1]}_{j}\setminus
A_{j}^{\circ},j}^{\circ}\Big{)}^{\top}\bm{\mathsf{X}}^{\top}_{A_{j}^{[t-1]}\setminus
A_{j}^{\circ}}\widehat{\bm{\xi}}/n.\end{split}$ (29)

Note that

$\Big{(}\widetilde{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{[t]}-\widehat{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{\circ}\Big{)}^{\top}\nabla\Big{\|}\widetilde{\bm{V}}^{[t]}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}\Big{\|}_{1}=\Big{\|}\widetilde{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{[t]}-\widehat{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{\circ}\Big{\|}_{1}.$

Then (29) is no greater than

$\begin{split}&\Big{\|}\widetilde{\bm{V}}_{A_{j}^{\circ}\triangle
A^{[t-1]}_{j},j}^{[t]}-\widehat{\bm{V}}_{A_{j}^{\circ}\triangle
A^{[t-1]}_{j},j}^{\circ}\Big{\|}_{1}\Big{(}\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}+\gamma_{j}\tau_{j}\Big{)}\\\
&+\Big{\|}\widetilde{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{[t]}-\widehat{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{\circ}\Big{\|}_{1}\Big{(}\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}-\gamma_{j}\tau_{j}\Big{)},\end{split}$
(30)

where $\triangle$ denotes the symmetric difference. Note that
$\|\bm{\mathsf{X}}(\widetilde{\bm{V}}_{\cdot j}^{[t]}-\widehat{\bm{V}}_{\cdot
j}^{\circ})\|_{2}^{2}/n\geq 0$. Rearranging the inequality yields that

$\begin{split}&\Big{(}\gamma_{j}\tau_{j}-\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}\Big{)}\Big{\|}\widetilde{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{[t]}-\widehat{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{\circ}\Big{\|}_{1}\\\
&\leq\Big{(}\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}+\gamma_{j}\tau_{j}\Big{)}\Big{\|}\widetilde{\bm{V}}_{A_{j}^{\circ}\triangle
A^{[t-1]}_{j},j}^{[t]}-\widehat{\bm{V}}_{A_{j}^{\circ}\triangle
A^{[t-1]}_{j},j}^{\circ}\Big{\|}_{1}.\end{split}$

On event $\mathscr{E}_{j}$,
$\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}\leq\gamma_{j}\tau_{j}/2$,
implying that $\|\widetilde{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{[t]}-\widehat{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{\circ}\|_{1}\leq
3\|\widetilde{\bm{V}}_{A_{j}^{\circ}\triangle
A^{[t-1]}_{j},j}^{[t]}-\widehat{\bm{V}}_{A_{j}^{\circ}\triangle
A^{[t-1]}_{j},j}^{\circ}\|_{1}\leq 3\|\widetilde{\bm{V}}_{A_{j}^{\circ}\cup
A^{[t-1]}_{j},j}^{[t]}-\widehat{\bm{V}}_{A_{j}^{\circ}\cup
A^{[t-1]}_{j},j}^{\circ}\|_{1}$. Note that $|A_{j}^{\circ}\cup
A^{[t-1]}_{j}|\leq 2\kappa^{\circ}_{\max}$. By Assumption 2 and (30),

$\begin{split}c_{1}\|\widetilde{\bm{V}}_{\cdot
j}^{[t]}-\widehat{\bm{V}}_{\cdot
j}^{\circ}\|_{2}^{2}\leq&\Big{(}\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}+\gamma_{j}\tau_{j}\Big{)}\Big{\|}\widetilde{\bm{V}}_{A_{j}^{\circ}\triangle
A_{j}^{[t-1]},j}^{[t]}-\widehat{\bm{V}}_{A_{j}^{\circ}\triangle
A_{j}^{[t-1]},j}^{\circ}\Big{\|}_{1}\\\
&+\Big{(}\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}-\gamma_{j}\tau_{j}\Big{)}\Big{\|}\widetilde{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{[t]}-\widehat{\bm{V}}_{(A_{j}^{\circ}\cup
A_{j}^{[t-1]})^{c},j}^{\circ}\Big{\|}_{1}\\\
\leq&\Big{(}\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}+\gamma_{j}\tau_{j}\Big{)}\Big{\|}\widetilde{\bm{V}}_{A_{j}^{\circ}\triangle
A_{j}^{[t-1]},j}^{[t]}-\widehat{\bm{V}}_{A_{j}^{\circ}\triangle
A_{j}^{[t-1]},j}^{\circ}\Big{\|}_{1}.\end{split}$ (31)

By the Cauchy-Schwarz inequality,
$\Big{\|}\widetilde{\bm{V}}_{A_{j}^{\circ}\triangle
A_{j}^{[t-1]},j}^{[t]}-\widehat{\bm{V}}_{A_{j}^{\circ}\triangle
A_{j}^{[t-1]},j}^{\circ}\Big{\|}_{1}\leq\sqrt{A_{j}^{\circ}\triangle
A^{[t-1]}_{j}}\Big{\|}\widetilde{\bm{V}}_{\cdot
j}^{[t]}-\widehat{\bm{V}}^{\circ}_{\cdot j}\Big{\|}_{2}$. Thus,
$c_{1}\|\widetilde{\bm{V}}_{\cdot j}^{[t]}-\widehat{\bm{V}}_{\cdot
j}^{\circ}\|_{2}\leq 1.5\gamma_{j}\tau_{j}\sqrt{2\kappa^{\circ}_{\max}}$,
since $|A_{j}^{\circ}\triangle A_{j}^{[t-1]}|\leq|A_{j}^{\circ}\cup
A_{j}^{[t-1]}|\leq 2\kappa^{\circ}_{\max}$ and
$\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}\leq
0.5\gamma_{j}\tau_{j}$. By the condition (1) of Theorem 4, $\gamma_{j}\leq
c_{1}/6$ and $\|\widetilde{\bm{V}}_{\cdot j}^{[t]}-\widehat{\bm{V}}_{\cdot
j}^{\circ}\|_{2}/\tau_{j}\leq\sqrt{\kappa^{\circ}_{\max}}$. On the other hand,
for any $l\in\text{FP}^{[t]}_{j}=A^{[t]}_{j}\setminus A_{j}^{\circ}$, we have
$|\widetilde{V}_{lj}^{[t]}-\widehat{V}_{lj}^{\circ}|=|\widetilde{V}_{lj}^{[t]}|>\tau_{j}$.
Thus, $\sqrt{|\text{FP}^{[t]}_{j}|}\leq\|\widetilde{\bm{V}}_{\cdot
j}^{[t]}-\widehat{\bm{V}}_{\cdot
j}^{\circ}\|_{2}/\tau_{j}\leq\sqrt{\kappa^{\circ}_{\max}}$. Consequently,
$|A_{j}^{\circ}\cup A^{[t]}_{j}|=|A_{j}^{\circ}|+|\text{FP}^{[t]}_{j}|\leq
2\kappa^{\circ}_{\max}$ on $\mathscr{E}_{j}$ for $t\geq 1$.

Next, we estimate the number of iterations required for termination. Note that
the machine precision is negligible. The termination criterion is met when
$A^{[t]}_{j}=A^{[t-1]}_{j}$ since the weighted Lasso problem (8) remains same
at the $t$th and $(t-1)$th iterations. To show that
$A^{[t]}_{j}=A^{\circ}_{j}$ eventually, we prove that
$|\text{FN}^{[t]}_{j}|+|\text{FP}^{[t]}_{j}|<1$ eventually.

Now, suppose $|\text{FN}^{[t]}_{j}|+|\text{FP}_{j}^{[t]}|\geq 1$. For any
$l\in\text{FN}^{[t]}_{j}\cup\text{FP}_{j}^{[t]}$, by Assumption 4,
$|\widetilde{V}^{[t]}_{lj}-\widehat{V}^{\circ}_{lj}|\geq|\widetilde{V}^{[t]}_{lj}-V^{\circ}_{lj}|-|\widehat{V}^{\circ}_{lj}-V^{\circ}_{lj}|\geq\tau_{j}-0.5\tau_{j}$,
so
$\sqrt{|\text{FN}^{[t]}_{j}|+|\text{FP}^{[t]}_{j}|}\leq\|\widetilde{\bm{V}}_{\cdot
j}^{[t]}-\widehat{\bm{V}}_{\cdot j}^{\circ}\|_{2}/0.5\tau_{j}$. By (31) and
the Cauchy-Schwarz inequality, $c_{1}\|\widetilde{\bm{V}}_{\cdot
j}^{[t]}-\widehat{\bm{V}}_{\cdot j}^{\circ}\|_{2}\leq
1.5\gamma_{j}\tau_{j}\sqrt{|\text{FN}_{j}^{[t-1]}|+|\text{FP}^{[t-1]}_{j}|}.$
By conditions (1) and (2) for $(\tau_{j},\gamma_{j})$ in Theorem 4, we have

$\begin{split}\sqrt{|\text{FN}^{[t]}_{j}|+|\text{FP}^{[t]}_{j}|}\leq\frac{\|\widetilde{\bm{V}}_{\cdot
j}^{[t]}-\widehat{\bm{V}}_{\cdot
j}^{\circ}\|_{2}}{0.5\tau_{j}}&\leq\frac{3\gamma_{j}}{c_{1}}\sqrt{|\text{FN}_{j}^{[t-1]}|+|\text{FP}_{j}^{[t-1]}|}\\\
&\leq 0.5\sqrt{|\text{FN}_{j}^{[t-1]}|+|\text{FP}_{j}^{[t-1]}|}.\end{split}$

Hence,
$\sqrt{|\text{FN}^{[t]}_{j}|+|\text{FP}^{[t]}_{j}|}\leq(1/2)^{t}\sqrt{|A^{\circ}_{j}|+|A^{[0]}_{j}|}$.
In particular, for $t\geq 1+\lceil\log\kappa^{\circ}_{j}/\log 4\rceil$,
$|\text{FN}_{j}^{[t]}|+|\text{FP}^{[t]}_{j}|<1$ implying that
$\text{FN}_{j}^{[t]}=\emptyset$ and $\text{FP}_{j}^{[t]}=\emptyset$.

Let $t_{\max}=1+\lceil\log\kappa^{\circ}_{j}/\log 4\rceil$. Then
$\text{FN}_{j}^{[t_{\max}]}=\text{FP}^{[t_{\max}]}_{j}=\emptyset$ on event
$\mathscr{E}_{j}$. Under condition (3), we have
$\\{l:\widetilde{V}^{[t_{\max}]}_{lj}\neq 0\\}=\\{l:V_{lj}^{\circ}\neq 0\\}$.
To bound $\operatorname{\mathbb{P}}(\bigcup_{j=1}^{p}\mathscr{E}_{j}^{c})$,
let
$\bm{\eta}=\bm{\mathsf{X}}^{\top}(\bm{I}-\bm{\mathsf{P}}_{A^{\circ}_{j}}){\bm{\xi}}_{j}$
and
$\bm{\eta}^{\prime}=n(\bm{\mathsf{X}}^{\top}_{A^{\circ}_{j}}\bm{\mathsf{X}}_{A^{\circ}_{j}})^{-1}\bm{\mathsf{X}}^{\top}_{A_{j}^{\circ}}{\bm{\xi}}_{j}$,
where
$\bm{\xi}_{j}=((\varepsilon_{V})_{1j},\ldots,(\varepsilon_{V})_{nj})^{\top}\in\mathbb{R}^{n}$.
Then $\bm{\eta}\in\mathbb{R}^{q}$ and
$\bm{\eta}^{\prime}\in\mathbb{R}^{|\kappa^{\circ}_{j}|}$ are Gaussian random
vectors, where $\mathop{\rm
Var}(\eta_{l})\leq\Omega_{jj}^{-1}(\bm{\mathsf{X}}^{\top}\bm{\mathsf{X}})_{ll}\leq
n\Omega_{jj}^{-1}c_{2}^{2}$ and $\mathop{\rm
Var}(\eta^{\prime}_{l})\leq\Omega_{jj}^{-1}((\bm{\mathsf{X}}^{\top}_{A^{\circ}_{j}}\bm{\mathsf{X}}_{A^{\circ}_{j}})^{-1})_{ll}\leq
n\Omega^{-1}_{jj}c_{2}^{2}$. Then

$\begin{split}\operatorname{\mathbb{P}}\Big{(}\|\bm{\mathsf{X}}^{\top}\widehat{\bm{\xi}}_{j}/n\|_{\infty}>0.5\gamma_{j}\tau_{j}\Big{)}&=\operatorname{\mathbb{P}}\Big{(}\|\bm{\eta}/n\|_{\infty}>0.5\gamma_{j}\tau_{j}\Big{)}\\\
&\leq\sum_{l=1}^{q}\operatorname{\mathbb{P}}\Big{(}|\eta_{l}|\geq
0.5n\gamma_{j}\tau_{j}\Big{)}\\\ &\leq
q\int_{\frac{0.5\sqrt{n\Omega_{jj}}\gamma_{j}\tau_{j}}{c_{2}}}^{\infty}\frac{e^{-t^{2}/2}}{\sqrt{2\pi}}dt\leq
q\sqrt{\frac{2}{\pi}}\exp\Big{(}-\frac{n\Omega_{jj}\gamma^{2}_{j}\tau^{2}_{j}}{8c^{2}_{2}}\Big{)},\end{split}$

and similarly,

$\begin{split}\operatorname{\mathbb{P}}\Big{(}\|\widehat{\bm{V}}_{\cdot
j}^{\circ}-\bm{V}_{\cdot
j}^{\circ}\|_{\infty}>0.5\tau_{j}\Big{)}\leq\sum_{l\in
A^{\circ}_{j}}\operatorname{\mathbb{P}}\Big{(}|\eta^{\prime}_{l}|>0.5\tau_{j}\Big{)}\leq\kappa^{\circ}_{\max}\sqrt{\frac{2}{\pi}}\exp\Big{(}-\frac{n\Omega_{jj}\tau^{2}_{j}}{8c^{2}_{2}}\Big{)}.\end{split}$

Under conditions (1) and (2) for $(\tau_{j},\gamma_{j})$ in Theorem 4,

$\begin{split}&\operatorname{\mathbb{P}}\Big{(}\Big{\\{}l:\widetilde{V}^{[t_{\max}]}_{lj}\neq
0\Big{\\}}=\Big{\\{}l:V_{lj}^{\circ}\neq 0\Big{\\}};\ 1\leq j\leq p\Big{)}\\\
&\geq
1-\operatorname{\mathbb{P}}\Big{(}\bigcup_{j=1}^{p}\mathscr{E}_{j}^{c}\Big{)}\\\
&\geq
1-p\sqrt{\frac{2}{\pi}}\Big{(}e^{-4(\log(q)+\log(n))+\log(q)}+e^{-144c_{1}^{-2}c_{2}^{2}(\log(q)+\log(n))+\log(\kappa^{\circ}_{j})}\Big{)}\geq
1-\sqrt{\frac{2}{\pi}}pq^{-3}n^{-4}.\end{split}$

By Borel-Cantelli lemma, $\\{l:\widetilde{V}^{[t_{\max}]}_{lj}\neq
0\\}=\\{l:V_{lj}^{\circ}\neq 0\\}$; $j=1,\ldots,p$ almost surely as
$n\to\infty$.

So far, we have $\widetilde{\bm{V}}_{\cdot
j}=\widetilde{\bm{V}}^{[t_{\max}]}_{\cdot j}=\widehat{\bm{V}}^{\circ}_{\cdot
j}$; $j=1,\ldots,p$, almost surely. It remains to show that
$\widehat{\bm{V}}^{\circ}_{\cdot j}$ is a global minimizer of (7);
$j=1,\ldots,p$, with a high probability. Note that Assumptions 2 and 4 imply
the degree of separation condition (Shen et al., 2013). By Theorem 2 of Shen
et al. (2013),

$\operatorname{\mathbb{P}}\Big{(}\widehat{\bm{V}}_{\cdot j}^{\circ}\text{ is
not a global minimizer of }\eqref{pseudo-likelihood};1\leq j\leq p\Big{)}\leq
3p\exp\Big{(}-2(\log(q)+\log(n))\Big{)}.$

implying that $\widehat{\bm{V}}_{\cdot
j}=\widetilde{\bm{V}}^{[t_{\max}]}_{\cdot j}=\widehat{\bm{V}}^{\circ}_{\cdot
j}$ is a global minimizer of (7) almost surely, as $n\to\infty$. This
completes the proof.

#### Proof of (B)

Assume $\widehat{\bm{V}}$ satisfies the properties in (A). Then Propositions 2
and 3 holds true for $\widehat{\bm{V}}$. Clearly, we have
$\widehat{\mathcal{C}}=\mathcal{C}^{\circ}$. Thus, we need to show
$\widehat{\mathcal{A}}=\mathcal{A}^{\circ}$ and
$\widehat{\mathcal{D}}=\mathcal{D}^{\circ}$. We shall prove this by induction.

For $h\geq 0$, the set of instruments on leaves in the graph $G^{[h]}$ (with
primary variables $\bm{Y}_{\\{j:h_{G}(j)\geq h\\}}$ and intervention variables
$\bm{X}$) is

$\mathscr{I}^{[h]}=\Big{\\{}l:l=\mathop{\rm
arg\,min\,}_{l^{\prime}:\|\widehat{\bm{V}}^{[h]}_{l^{\prime}\cdot}\|_{0}\neq
0}\|\widehat{\bm{V}}^{[h]}_{l^{\prime}\cdot}\|_{0}\Big{\\}}=\Big{\\{}l:\|\bm{V}^{[h]}_{l\cdot}\|_{0}=1\Big{\\}}.$

Hence, $X_{l}$ is an instrument of leaf variable $Y_{j}$ in $G^{[h]}$, when
$l\in\mathscr{I}^{[h]}$ and $j=\mathop{\rm
arg\,max\,}_{j^{\prime}}|V_{lj^{\prime}}|$. Moreover, $h_{G}(j)=h$. By
Proposition 3 it also recovers the edges $Y_{k}\to Y_{j}$ such that
$h_{G}(k)=h_{G}(j)+1$. In model (1), after removing $\\{Y_{j}:h_{G}(j)\leq
h\\}$, the local Markov property among $(\bm{Y}_{\\{j:h_{G}(j)\geq
h+1\\}},\bm{X})$ remain intact provided that $\\{j:h_{G}(j)\geq
h+1\\}\neq\emptyset$. Then we increment $h$ and repeat the procedure until all
primary variables are removed, namely, $\\{j:h_{G}(j)\geq h+1\\}=\emptyset$.

As a result, Algorithm 2 correctly identifies the topological height
$h_{G}(j)$ for each primary variable $Y_{j}$ and the edges $Y_{k}\to Y_{j}$
such that $h_{G}(k)=h_{G}(j)+1$. Consequently, the ancestral relations can be
reconstructed by

$\mathcal{A}=\Big{\\{}(k_{0},k_{d}):Y_{k_{0}}\to\cdots\to Y_{k_{d}},\
h_{G}(k_{\nu})=h_{G}(k_{\nu+1})+1,\ \nu=0,\ldots,d-1\Big{\\}}.$

Therefore, Step 6 of Algorithm 2 recovers
$\widehat{\mathcal{A}}=\mathcal{A}^{\circ}$ and
$\widehat{\mathcal{D}}=\mathcal{D}^{\circ}$ almost surely as $n\to\infty$.
This completes the proof of (B).

###### Lemma 10

Let

$T_{j}=\frac{\bm{\mathsf{e}}_{j}^{\top}(\bm{\mathsf{P}}_{A^{\circ}_{j}}-\bm{\mathsf{P}}_{B^{\circ}_{j}})\bm{\mathsf{e}}_{j}}{\bm{\mathsf{e}}_{j}^{\top}(\bm{I}-\bm{\mathsf{P}}_{A^{\circ}_{j}})\bm{\mathsf{e}}_{j}/(n-|A^{\circ}_{j}|)},\quad
T_{j}^{*}=\frac{(\bm{\mathsf{e}}^{*}_{j})^{\top}(\bm{\mathsf{P}}_{A^{\circ}_{j}}-\bm{\mathsf{P}}_{B^{\circ}_{j}})\bm{\mathsf{e}}^{*}_{j}}{(\bm{\mathsf{e}}^{*}_{j})^{\top}(\bm{I}-\bm{\mathsf{P}}_{A^{\circ}_{j}})\bm{\mathsf{e}}^{*}_{j}/(n-|A^{\circ}_{j}|)},$
(32)

where
$A^{\circ}_{j}=\textnormal{{an}}_{S^{\circ}}(j)\cup\textnormal{{in}}_{S^{\circ}}(j)\cup\textnormal{{d}}_{S^{\circ}}(j)$,
$B^{\circ}_{j}=(\textnormal{{an}}_{S^{\circ}}(j)\cup\textnormal{{in}}_{S^{\circ}}(j))\setminus\textnormal{{d}}_{S^{\circ}}(j)$.
Then $(T_{1},\ldots,T_{p})$ are independent and $(T_{1}^{*},\ldots,T_{p}^{*})$
are conditionally independent given $\bm{\mathsf{Z}}$. Moreover,

$T^{*}_{j}/(|A^{\circ}_{j}|-|B^{\circ}_{j}|)\mid\bm{\mathsf{Z}}\overset{d}{=}T_{j}/(|A^{\circ}_{j}|-|B^{\circ}_{j}|)\sim
F_{|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|};\quad 1\leq j\leq p,$

where $F_{d_{1},d_{2}}$ denotes the F-distribution with $d_{1}$ and $d_{2}$
degrees of freedom.

### Proof of Lemma 10

Let $\bm{\mathsf{Z}}=(\bm{\mathsf{Y}},\bm{\mathsf{X}})$ as in (12). Given data
submatrix $\bm{\mathsf{Z}}_{A^{\circ}_{j}}$,

$\frac{\bm{\mathsf{e}}_{j}^{\top}(\bm{\mathsf{P}}_{A^{\circ}_{j}}-\bm{\mathsf{P}}_{B^{\circ}_{j}})\bm{\mathsf{e}}_{j}}{\sigma^{2}_{j}}\mid\bm{\mathsf{Z}}_{A^{\circ}_{j}}\sim\chi^{2}_{|A^{\circ}_{j}|-|B^{\circ}_{j}|},\quad\frac{\bm{\mathsf{e}}_{j}^{\top}(\bm{I}-\bm{\mathsf{P}}_{A^{\circ}_{j}})\bm{\mathsf{e}}_{j}}{\sigma^{2}_{j}}\mid\bm{\mathsf{Z}}_{A^{\circ}_{j}}\sim\chi^{2}_{n-|A^{\circ}_{j}|},$

and they are independent, because
$\bm{\mathsf{P}}_{A^{\circ}_{j}}-\bm{\mathsf{P}}_{B^{\circ}_{j}}$ and
$\bm{I}-\bm{\mathsf{P}}_{A^{\circ}_{j}}$ are orthonormal projection matrices,
$\bm{\mathsf{e}}_{j}\mid\bm{\mathsf{Z}}_{A^{\circ}_{j}}\sim
N(\bm{0},\sigma_{j}^{2}\bm{I}_{n})$, and
$(\bm{\mathsf{P}}_{A^{\circ}_{j}}-\bm{\mathsf{P}}_{B^{\circ}_{j}})(\bm{I}-\bm{\mathsf{P}}_{A^{\circ}_{j}})=\bm{0}$.
Then, for any $t\in\mathbb{R}$, the characteristic function
$\operatorname{\mathbb{E}}\exp(\iota tT_{j}/(|A^{0}_{j}|-|B^{0}_{j}|))$ can be
written as ($\iota$ is the imaginary unit)

$\operatorname{\mathbb{E}}\left(\operatorname{\mathbb{E}}\left(\exp(\iota
tT_{j}/(|A^{\circ}_{j}|-|B^{\circ}_{j}|))\mid\bm{\mathsf{Z}}_{A^{\circ}_{j}}\right)\right)=\operatorname{\mathbb{E}}\psi_{|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|}(t)=\psi_{|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|}(t),$

where $\psi_{|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|}$ is the
characteristic function of F-distribution with degrees of freedom
$|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|$. Hence,
$T_{j}/(|A^{\circ}_{j}|-|B^{\circ}_{j}|)\sim
F_{|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|}$. Similarly,
$T_{j}^{*}/(|A^{\circ}_{j}|-|B^{\circ}_{j}|)\sim
F_{|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|}$; $j=1,\ldots,p$.

Next, we prove independence for $\bm{T}$ and $\bm{T}^{*}$ via a peeling
argument. Let $\bm{t}=(t_{1},\ldots,t_{p})$. Let $j$ be a leaf node of the
graph $G^{\circ}$. Then $\bm{T}_{-j}\mid\bm{\mathsf{Y}}_{-j},\bm{\mathsf{X}}$
is deterministic, where $\bm{T}_{-j}$ is the subvector of $\bm{T}$ with the
$j$th component removed. The characteristic function of $\bm{t}^{\top}\bm{T}$
is

$\begin{split}\operatorname{\mathbb{E}}\exp(\iota\bm{t}^{\top}\bm{T})&=\operatorname{\mathbb{E}}\Big{(}\operatorname{\mathbb{E}}\Big{(}\exp(\iota
t_{j}T_{j})\mid\bm{\mathsf{Y}}_{-j},\bm{\mathsf{X}}\Big{)}\exp(\iota\bm{t}^{\top}_{-j}\bm{T}_{-j})\Big{)}\\\
&=\psi_{|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|}(t_{j})\operatorname{\mathbb{E}}\exp(\iota\bm{t}_{-j}^{\top}\bm{T}_{-j}),\end{split}$
(33)

where $\psi_{|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|}$ is the
characteristic function of the F-distribution with degrees of freedom
$(|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|)$. Next, let $j^{\prime}$
be a leaf node of the graph $G^{\prime}$ and apply the law of iterated
expectation again, where $G^{\prime}$ is the subgraph of $G^{\circ}$ without
node $j$. Repeat this procedure and after $p$ steps
$\operatorname{\mathbb{E}}\exp(\iota\bm{t}^{\top}\bm{T})=\prod_{j=1}^{p}\psi_{|A^{\circ}_{j}|-|B^{\circ}_{j}|,n-|A^{\circ}_{j}|}(t_{j})$,
which implies $(T_{1},\ldots,T_{p})$ are independent. Similarly, $\bm{T}^{*}$
also has independent components given $\bm{\mathsf{Z}}$ and has the same
distribution as $\bm{T}$. This completes the proof.

###  B.5 Proof of Theorem 5

#### Proof of (A)

Without loss of generality, assume $M$ is sufficiently large. For any
$x\in\mathbb{R}$,

$\begin{split}&|\operatorname{\mathbb{P}}(\textnormal{Lr}\leq
x)-\operatorname{\mathbb{P}}(\textnormal{Lr}(S^{\circ},\widehat{\bm{\Sigma}})\leq
x)|\\\ &\leq\operatorname{\mathbb{E}}|\operatorname{I}(\textnormal{Lr}\leq
x)-\operatorname{I}(\textnormal{Lr}(S^{\circ},\widehat{\bm{\Sigma}})\leq
x)|\\\ &=\operatorname{\mathbb{E}}|\operatorname{I}(\textnormal{Lr}\leq
x,\widehat{S}\neq
S^{\circ})+\operatorname{I}(\textnormal{Lr}(S^{\circ},\widehat{\bm{\Sigma}})\leq
x)(\operatorname{I}(\widehat{S}=S^{\circ})-1)|\\\ &\leq
2\operatorname{\mathbb{P}}(\widehat{S}\neq S^{\circ}).\end{split}$ (34)

From (13),
$\textnormal{Lr}^{*}(S^{\circ},\widehat{\bm{\Sigma}}^{*})=\sum_{j:\textnormal{{d}}_{S^{\circ}}(j)\neq\emptyset}T_{j}^{*}$.
By Lemma 10,
$\operatorname{\mathbb{P}}(\textnormal{Lr}(S^{\circ},\widehat{\bm{\Sigma}})\leq
x)=\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}(S^{\circ},\widehat{\bm{\Sigma}}^{*})\leq
x\mid\bm{\mathsf{Z}})$. Note that
$\textnormal{Lr}^{*}=\textnormal{Lr}^{*}(\widehat{S}^{*},\widehat{\bm{\Sigma}}^{*})$.
Then for any $x\in\mathbb{R}$,

$\begin{split}&|\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}\leq
x\mid\bm{\mathsf{Z}})-\operatorname{\mathbb{P}}(\textnormal{Lr}(S^{\circ},\widehat{\bm{\Sigma}})\leq
x)|\\\ &=|\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}\leq
x,\widehat{S}^{*}\neq
S^{\circ}\mid\bm{\mathsf{Z}})+\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}\leq
x,\widehat{S}^{*}=S^{\circ}\mid\bm{\mathsf{Z}})-\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}(S^{\circ},\widehat{\bm{\Sigma}}^{*})\leq
x\mid\bm{\mathsf{Z}})|\\\
&\leq|\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}\leq x,\widehat{S}^{*}\neq
S^{\circ}\mid\bm{\mathsf{Z}})|+|\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}(S^{\circ},\widehat{\bm{\Sigma}}^{*})\leq
x,\widehat{S}^{*}\neq S^{\circ}\mid\bm{\mathsf{Z}})|\\\ &\leq
2\operatorname{\mathbb{P}}(\widehat{S}^{*}\neq
S^{\circ}\mid\bm{\mathsf{Z}}).\end{split}$ (35)

Since (34) and (35) hold uniformly in $x$, we have

$\sup_{x\in\mathbb{R}}|\operatorname{\mathbb{P}}(\textnormal{Lr}\leq
x)-\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}\leq
x\mid\bm{\mathsf{Z}})|\leq 2\operatorname{\mathbb{P}}(\widehat{S}\neq
S^{\circ})+2\operatorname{\mathbb{P}}(\widehat{S}^{*}\neq
S^{\circ}\mid\bm{\mathsf{Z}}).$

By Theorem 4, we have $\operatorname{\mathbb{P}}(\widehat{S}\neq S^{\circ})\to
0$, which holds uniformly for all $\bm{\theta}$ which satisfy the Assumptions
1C-5. For $\operatorname{\mathbb{P}}(\widehat{S}^{*}\neq
S^{\circ}\mid\bm{\mathsf{Z}})$, the error terms of (7) in the perturbed data
$\bm{\mathsf{Z}}^{*}$ are rescaled with
$\Omega_{jj}^{-1}+\widehat{\sigma}^{2}_{j}\leq 2\Omega_{jj}^{-1}$. By Theorem
4, $\operatorname{\mathbb{P}}(\widehat{S}^{*}\neq
S^{\circ})=\operatorname{\mathbb{E}}\operatorname{\mathbb{P}}(\widehat{S}^{*}\neq
S^{\circ}\mid\bm{\mathsf{Z}})\to 0$, which implies
$\operatorname{\mathbb{P}}(\widehat{S}^{*}\neq
S^{\circ}\mid\bm{\mathsf{Z}})\overset{p}{\longrightarrow}0$ as
$n\rightarrow\infty$ by the Markov inequality. Consequently,
$\sup_{x\in\mathbb{R}}|\operatorname{\mathbb{P}}(\textnormal{Lr}\leq
x)-\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}\leq x\mid\bm{\mathsf{Z}})|\to
0$. For $|\mathcal{D}|=0$, $\operatorname{\mathbb{P}}(\textnormal{Lr}=0)\to
1$, $\operatorname{\mathbb{P}}(\textnormal{Lr}^{*}=0\mid\bm{\mathsf{Z}})\to
1$, and $\operatorname{\mathbb{P}}(\textnormal{Pval}=1)\to 1$. For
$|\mathcal{D}|>0$,
$P(\textnormal{Lr}^{*}\geq\textnormal{Lr}\mid\bm{\mathsf{Z}})\to\textnormal{Unif}(0,1)$
and $\operatorname{\mathbb{P}}(\textnormal{Pval}<\alpha)\to\alpha$. This
completes the proof of (A).

#### Proof of (B)

Let
$\textnormal{Pval}_{k}=M^{-1}\sum_{m=1}^{M}\operatorname{I}(\text{Lr}^{*}_{k,m}\geq\text{Lr})$,
the p-value of sub-hypothesis $\text{H}_{0,k}$. For
$|\mathcal{D}|<|\mathcal{H}|$, there exists an edge
$(i_{k},j_{k})\in\mathcal{H}$ but $(i_{k},j_{k})\notin\mathcal{D}$. Then by
(A), $\operatorname{\mathbb{P}}(\textnormal{Pval}=\textnormal{Pval}_{k}=1)\to
1$. For $|\mathcal{D}|=|\mathcal{H}|$, note that as $n,M\to\infty$,

$\operatorname{\mathbb{P}}\Big{(}\textnormal{Pval}<\alpha\Big{)}=\operatorname{\mathbb{P}}\Big{(}\textnormal{Pval}_{1}<\alpha,\ldots,\textnormal{Pval}_{|\mathcal{H}|}<\alpha\Big{)}\leq\operatorname{\mathbb{P}}\Big{(}\textnormal{Pval}_{1}<\alpha\Big{)}\to\alpha.$

Now, define a sequence $\bm{U}_{\mathcal{H}}^{(r)}$; $r=1,2,\ldots$ such that
$\bm{U}_{(i_{1},j_{1})}^{(r)}=0$ and
$\min_{k=2}^{|\mathcal{H}|}|\bm{U}_{(i_{k},j_{k})}^{(r)}|\geq c>0$. Thus,
$\bm{U}_{\mathcal{H}}^{(r)}$; $r=1,2,\ldots$ satisfies $H_{0}$. By Proposition
7, $\text{Pval}_{k}\overset{p}{\longrightarrow}0$ for $k\geq 2$ as
$r\to\infty$ . Hence,

$\limsup_{\begin{subarray}{c}n\to\infty\\\ \bm{U}_{\mathcal{H}}^{(r)}\text{
satisfies
H}_{0}\end{subarray}}\operatorname{\mathbb{P}}_{\bm{\theta}^{(r)}}\Big{(}\text{Pval}<\alpha\Big{)}\geq\sup_{r}\lim_{n,p,q\to\infty}\operatorname{\mathbb{P}}_{\bm{\theta}^{(r)}}\Big{(}\text{Pval}<\alpha\Big{)}=\alpha.$

This completes the proof.

###  B.6 Proof of Proposition 6

Since $\operatorname{\mathbb{P}}(\widehat{S}=S^{\circ})\to 1$, it suffices to
consider $\text{Lr}(S^{\circ},\widehat{\bm{\Sigma}})$. For $|\mathcal{D}|=0$,
we have $\operatorname{\mathbb{P}}(\text{Lr}=0)\to 1$. Now, assume
$|\mathcal{D}|>0$. Then

$\begin{split}2\text{Lr}(S^{\circ},\widehat{\bm{\Sigma}})=\underbrace{\sum_{j:\textnormal{{d}}_{S^{\circ}}(j)\neq\emptyset}\frac{\bm{\mathsf{e}}_{j}^{\top}(\bm{\mathsf{P}}_{A^{\circ}_{j}}-\bm{\mathsf{P}}_{B^{\circ}_{j}})\bm{\mathsf{e}}_{j}}{\sigma^{2}_{j}}}_{R_{1}}+\underbrace{\sum_{j:\textnormal{{d}}_{S^{\circ}}(j)\neq\emptyset}\left(\frac{\sigma^{2}_{j}}{\widehat{\sigma}^{2}_{j}}-1\right)\frac{\bm{\mathsf{e}}_{j}^{\top}(\bm{\mathsf{P}}_{A^{\circ}_{j}}-\bm{\mathsf{P}}_{B^{\circ}_{j}})\bm{\mathsf{e}}_{j}}{\sigma_{j}^{2}}}_{R_{2}}.\end{split}$

To derive the asymptotic distribution of $R_{1}$, we apply the peeling
strategy with law of iterated expectation as in the proof of Lemma 10. Then we
have
$\bm{\mathsf{e}}_{j}^{\top}(\bm{\mathsf{P}}_{A^{0}_{j}}-\bm{\mathsf{P}}_{B^{0}_{j}})\bm{\mathsf{e}}_{j}/\sigma^{2}_{j}$;
$\textnormal{{d}}_{S^{\circ}}(j)\neq\emptyset$ are independent. Therefore,
$R_{1}\sim\chi^{2}_{|\mathcal{D}|}$.

To bound $R_{2}$, we apply Lemma 1 of Laurent and Massart (2000). By
Assumption 5,

$\begin{split}\operatorname{\mathbb{P}}\left(\max_{j:\textnormal{{d}}_{S^{\circ}}(j)\neq\emptyset}\left|\frac{\widehat{\sigma}^{2}_{j}}{\sigma^{2}_{j}}-1\right|\geq
4\sqrt{\frac{\log|\mathcal{D}|}{(1-\rho)n}}+8\frac{\log|\mathcal{D}|}{(1-\rho)n}\right)&\leq
2\exp(-\log|\mathcal{D}|).\end{split}$

Hence,
$\max_{j:\textnormal{{d}}_{S^{\circ}}(j)\neq\emptyset}\left|\frac{\sigma^{2}_{j}}{\widehat{\sigma}_{j}^{2}}-1\right|\leq
8\sqrt{\frac{\log|\mathcal{D}|}{(1-\rho)n}}$ with probability tending one as
$n,p\rightarrow\infty$. Thus

$|R_{2}|\leq|R_{1}|\max_{j:\textnormal{{d}}_{S^{\circ}}(j)\neq\emptyset}\left|\frac{\sigma^{2}_{j}}{\widehat{\sigma}_{j}^{2}}-1\right|\leq
O_{\operatorname{\mathbb{P}}}\left(|\mathcal{D}|\sqrt{\frac{\log|\mathcal{D}|}{n}}\right).$

Consequently, as $n,p,q\rightarrow\infty$, if $|\mathcal{D}|>0$ is fixed, then
$2\text{Lr}\overset{d}{\longrightarrow}\chi^{2}_{|\mathcal{D}|}$; if
$|\mathcal{D}|\to\infty$, then
$(2|\mathcal{D}|)^{-1/2}(2\text{Lr}-|\mathcal{D}|)\overset{d}{\longrightarrow}N(0,1)$.
This completes the proof.

###  B.7 Proof of Proposition 7

Let $\bm{\theta}^{(n)}=(\bm{U}^{\circ}+\bm{\Delta},\bm{W}^{\circ})$. Then

$L(\bm{\theta}^{(n)},\bm{\Sigma}^{\circ})-L(\bm{\theta}^{\circ},\bm{\Sigma}^{\circ})=\sum_{j:\textnormal{{d}}^{\circ}(j)\neq\emptyset}\Big{(}\sqrt{n}\bm{\eta}_{j}^{\top}\bm{\Delta}_{\cdot
j}-\frac{1}{2}\sqrt{n}\bm{\Delta}_{\cdot
j}\Big{(}n^{-1}\bm{\mathsf{Y}}^{\top}_{\textnormal{{d}}^{\circ}(j)}\bm{\mathsf{Y}}_{\textnormal{{d}}^{\circ}(j)}\Big{)}\sqrt{n}\bm{\Delta}_{\cdot
j}\Big{)},\\\ $

where
$\bm{\eta}_{j}=(n^{-1/2}\bm{\mathsf{Y}}_{\textnormal{{d}}^{\circ}(j)}^{\top}\bm{\mathsf{e}}_{j},\bm{0}_{|\textnormal{{d}}^{\circ}(j)^{c}|})\in\mathbb{R}^{p}$.
For the likelihood ratio, it suffices to consider
$\textnormal{Lr}(S^{\circ},\widehat{\bm{\Sigma}})$, since
$P(\widehat{S}=S^{\circ})\to 1$. Under $H_{0}$, we have
$2\text{Lr}=\sum_{j:\textnormal{{d}}^{\circ}(j)\neq\emptyset}T_{j}$. Note that
$T_{j}=\sum_{r=1}^{|\textnormal{{d}}^{\circ}(j)|}(\bm{\mathsf{q}}_{j,r}^{\top}\bm{\mathsf{e}}_{j})^{2}+o_{\operatorname{\mathbb{P}}}(1)$
by Proposition 6, where
$\bm{\mathsf{P}}_{A^{\circ}_{j}}-\bm{\mathsf{P}}_{B^{\circ}_{j}}=\sum_{r=1}^{|\textnormal{{d}}^{\circ}(j)|}\bm{\mathsf{q}}_{j,r}\bm{\mathsf{q}}_{j,r}^{\top}$.
Let
$\bm{\mathsf{Q}}_{j}=(\bm{\mathsf{q}}_{j,1},\ldots,\bm{\mathsf{q}}_{j,|\textnormal{{d}}^{\circ}(j)|})\in\mathbb{R}^{n\times|\textnormal{{d}}^{\circ}(j)|}$.
Then

$\begin{pmatrix}\bm{\mathsf{Q}}_{j}^{\top}\bm{\mathsf{e}}_{j}\\\
\bm{\eta}_{j}\end{pmatrix}\ \vline\
\bm{\mathsf{Y}}_{\textnormal{{an}}(j)},\bm{\mathsf{X}}\sim
N\left(\bm{0},\begin{pmatrix}\bm{I}_{|\textnormal{{d}}^{\circ}(j)|}&n^{-1/2}\bm{\mathsf{Q}}_{j}^{\top}\bm{\mathsf{Y}}_{\textnormal{{d}}^{\circ}(j)}\\\
n^{-1/2}\bm{\mathsf{Y}}_{\textnormal{{d}}^{\circ}(j)}^{\top}\bm{\mathsf{Q}}_{j}&n^{-1}\bm{\mathsf{Y}}^{\top}_{\textnormal{{d}}^{\circ}(j)}\bm{\mathsf{Y}}_{\textnormal{{d}}^{\circ}(j)}\end{pmatrix}\right),$
$\bm{\mathsf{Q}}^{\top}_{j}\bm{\mathsf{e}}_{j}\mid\bm{\eta}_{j},\bm{\mathsf{Y}}_{\textnormal{{an}}(j)},\bm{\mathsf{X}}\sim
N\Big{(}n^{-1/2}\bm{\mathsf{Y}}_{\textnormal{{d}}^{\circ}(j)}^{\top}\bm{\mathsf{Q}}_{j}\bm{\eta}_{j},\bm{\mathsf{Y}}_{\textnormal{{d}}^{\circ}(j)}^{\top}(\bm{I}_{n}-\bm{\mathsf{Q}}_{j}\bm{\mathsf{Q}}_{j}^{\top})\bm{\mathsf{Y}}_{\textnormal{{d}}^{\circ}(j)}\Big{)}.$

Next, let $k$ be a leaf node of the graph $G^{\circ}$. For fixed
$|\mathcal{D}|>0$, after change of measure,

$\begin{split}&\beta(\bm{\theta}^{\circ},\bm{\Delta})\\\
&\geq\liminf_{n\to\infty}\operatorname{\mathbb{E}}_{\bm{\theta}^{\circ}}\left(\operatorname{I}\left(\sum_{j=1}^{p}\|\bm{\mathsf{Q}}^{\top}\bm{\mathsf{e}}_{j}\|^{2}_{2}>\chi^{2}_{|\mathcal{D}|,1-\alpha}\right)\exp(L(\bm{\theta}^{(n)},\bm{\Sigma}^{\circ})-L(\bm{\theta}^{\circ},\bm{\Sigma}^{\circ}))\right)\\\