High-Dimensional Joint Estimation
of Multiple Directed Gaussian Graphical Models
Yuhao Wang1,2
1

Santiago Segarra2

Caroline Uhler1,2

Laboratory for Information and Decision Systems, Massachusetts Institute of Technology
2 Institute for Data, Systems, and Society, Massachusetts Institute of Technology

arXiv:1804.00778v1 [] 3 Apr 2018

Abstract
We consider the problem of jointly estimating multiple related directed acyclic graph (DAG) models
based on high-dimensional data from each graph. This problem is motivated by the task of learning
gene regulatory networks based on gene expression data from different tissues, developmental stages or
disease states. We prove that under certain regularity conditions, the proposed `0 -penalized maximum
likelihood estimator converges in Frobenius norm to the adjacency matrices consistent with the datagenerating distributions and has the correct sparsity. In particular, we show that this joint estimation
procedure leads to a faster convergence rate than estimating each DAG model separately. As a corollary
we also obtain high-dimensional consistency results for causal inference from a mix of observational
and interventional data. For practical purposes, we propose jointGES consisting of Greedy Equivalence
Search (GES) to estimate the union of all DAG models followed by variable selection using lasso to
obtain the different DAGs, and we analyze its consistency guarantees. The proposed method is illustrated
through an analysis of simulated data as well as epithelial ovarian cancer gene expression data.

1

Introduction

Directed acyclic graph (DAG) models, also known as Bayesian networks, are widely used to model causal
relationships in complex systems across various fields such as computational biology, epidemiology, sociology, and environmental management [1, 11, 30, 35, 39]. In these applications we often encounter highdimensional datasets where the number of variables or nodes greatly exceeds the number of observations.
While the problem of structure identification for undirected graphical models in the high-dimensional setting is quite well understood [34, 25, 10, 4, 44], such results are just starting to become available for directed
graphical models. The difficulty in identifying DAG models can be attributed to the fact that searching over
the space of DAGs is NP-complete in general [5].
Methods for structure identification in directed graphical models can be divided into two categories
and hybrids of these categories. Constraint-based methods, such as the prominent PC algorithm, first
learn an undirected graph from conditional independence relations and in a second step orient some of
the edges [14, 39]. Score-based methods, on the other hand, posit a scoring criterion for each DAG model,
usually a penalized likelihood score, and then search for the network with the highest score given the observations. An example is the celebrated Greedy Equivalence Search (GES) algorithm, which can be used to
greedily optimize the `0 -penalized likelihood such as the Bayesian Information Criterion (BIC) [6]. Highdimensional consistency guarantees were recently obtained for the PC algorithm [19] and for score-based
methods [23, 28, 43].
Existing methods have focused on estimating a single directed graphical model. However, in many applications we have access to data from related classes, such as gene expression data from different tissues,
cell types or states [24, 37], different developmental stages [2], different disease states [41], or from different
perturbations such as knock-out experiments [8]. In all these applications one would expect that the underlying regulatory networks are similar to each other, since they stem from the same species, individual or cell
1

type, but also have important differences that drive differentiation, development or a certain disease. This
raises an important statistical question, namely how to jointly estimate related directed graphical models in
order to effectively make use of the available data.
Various methods have been proposed for jointly estimating undirected Gaussian graphical models. To
preserve the common structure, Guo et al. [15] suggested to use a hierarchical penalty and Danaher et
al. [7] suggested the use of a generalized fused lasso or group lasso penalty. While both approaches achieve
the same convergence rate as the individual estimators, Cai et al. [3] were able to improve the asymptotic
convergence rate of joint estimation using a weighted constrained `∞ /`1 minimization approach. Bayesian
methods have been proposed for this problem [32] as well. Related works also include [27], where it is
assumed that the networks differ only locally in a few nodes and [21, 38], where the assumption is that the
networks are ordered and related by continuously changing edge weights.
In this paper, we propose a framework based on `0 -penalized maximum likelihood estimation for jointly
estimating related directed Gaussian graphical models. We show that the joint `0 -penalized maximum likelihood estimator (MLE) achieves a faster asymptotic convergence rate as compared to the individual estimators. In addition, by viewing interventional data as data coming from a related network, we show that
the interventional BIC scoring function proposed in [16] can be obtained as a special case of the joint
`0 -penalized maximum likelihood approach presented here. Our theoretical consistency guarantees also explain the empirical findings of [16], namely that estimating a DAG model from interventional data usually
leads to better recovery rates as compared to estimating a DAG model from the same amount of purely
observational data. These theoretical results are based on the global optimum of `0 -penalized maximum
likelihood estimation. To overcome the computational bottleneck of this optimization problem we propose
a greedy approach (jointGES) for solving this problem by extending GES [6] to the joint estimation setting.
We analyze its properties from a theoretical point of view and test its performance on synthetic data and
gene expression data from epithelial ovarian cancer.
The remainder of this paper is structured as follows. In Section 2, we review some relevant background
related to DAG models and introduce notation for the joint DAG estimation problem studied in this paper. In
Section 3, we present the joint `0 -penalized maximum likelihood estimator and jointGES, an adaptation of
GES for solving this optimization problem. Section 4 establishes results regarding the statistical consistency
of the `0 -penalized MLE and jointGES whereas Section 5 presents the implications for learning DAG models
from a mix of observational and interventional data. In Section 6, we illustrate the performance of our
proposal in a simulation study and an application to the analysis of gene expression data. We conclude with
a short discussion in Section 7. The proofs of supporting results are contained in the Appendix.

2

Preliminaries

In Section 2.1 we introduce DAGs, their relation to linear structural equation models, and discuss statistical
features enjoyed by random vectors following these models. In Section 2.2 we briefly review existing
approaches for estimating a single directed graphical model from observational data. Finally, Section 2.3
describes a setting where multiple related directed graphical models exist.

2.1

Directed acyclic graphs and linear structural equation models

Let G = (V, E) denote a DAG with vertices V = [p] = {1, · · · , p} and directed edges E ⊆ V × V , where
|G| denotes the cardinality of E. Let A ∈ Rp×p be the adjacency matrix specifying the edge weights of the
underlying DAG G, i.e., Aij 6= 0 if and only if (i, j) ∈ E. Also, let  ∼ N (0, Ω) denote a p-dimensional
2

multivariate Gaussian random variable with zero mean and diagonal covariance matrix Ω. In this work, we
assume that the observed random vector X = (X1 , · · · , Xp ) ∈ Rp is generated according to the following
linear structural equation model (SEM).
X = AT X + .

(1)

Hence X follows a multivariate Gaussian distribution with zero mean and covariance matrix Σ, where the
inverse covariance (or precision) matrix Θ = Σ−1 is given by
Θ = (I − A)Ω−1 (I − A)T .

(2)

Let Paj (G) denote the parents of node j in G, then it follows from (1) that the distribution of X can be
factorized as
P(X) =

p
Y

P(Xj |XPaj (G) ).

j=1

In addition, such a factorization of P according to G is equivalent to the Markov assumption with respect to
G [22, Theorem 3.27]. Formally, given i, j ∈ V and an arbitrary subset of nodes S ⊂ V \ {i, j}, then
j is d-separated from k | S in G

⇒

Xj ⊥
⊥ Xk |XS in P.

(3)

If the implication (3) holds bidirectionally, then P is said to be faithful [39] with respect to G. Note that there
exist DAGs G1 and G2 that encode the same d-separations and hence the same conditional independence
relations. Such DAGs are said to belong to the same Markov equivalence class.
A consequence of the acyclicity of G is that there exists at least one permutation π of [p] such that
Aij = 0 for all π(i) ≥ π(j). Putting it differently, if the rows and columns of A are reordered according
to π, then the resulting matrix is strictly upper triangular. Hence, if such a permutation π is known a priori,
one can obtain the SEM parameters (A, Ω) from Θ according to the following steps [cf. (2)]. First, we
reorder Θ according to π. Then, we perform on the reordered Θ an upper-triangular-plus-diagonal Cholesky
decomposition to obtain (A0 , Ω0 ). Finally, we revert the ordering by permuting the rows and columns of
A0 and Ω0 according to π −1 and obtain the sought (A, Ω). For an arbitrary permutation π and a given Θ,
we denote by (Aπ , Ωπ ) the Cholesky decomposition parameters obtained from the procedure just described.
Alternatively, one can obtain (Aπ , Ωπ ) by solving p linear regressions [cf. (1)]. More precisely, we can
obtain each column of Aπ by regressing Xj only on those Xi such that π(i) < π(j) for all j. Once Aπ is
obtained, one can estimate the variance of  in (1) to get Ωπ . In the rest of the paper, we denote by (A0 , Ω0 )
and Θ0 the true parameters of the data-generating SEM and the associated precision matrix, respectively.
Moreover, we denote by (A0π , Ω0π ) the SEM parameters obtained from the described procedure when the
true precision matrix Θ0 is used. Notice that (A0 , Ω0 ) = (A0π , Ω0π ) if π is any permutation consistent
with the true underlying DAG G0 . The DAG Gπ corresponding to the non-zero entries of Aπ is known as
the minimal I-MAP (independence map) with respect to π. The minimal I-MAP with the fewest number of
edges is called minimal-edge I-MAP [43]. If P is faithful with respect to a DAG G, then G is a minimal-edge
I-MAP of P [33, 43].
Furthermore, it has been shown in [31] that all DAGs in a Markov equivalence class share the same
skeleton – i.e., the set of edges when directions are ignored – and v-structures. A v-structure is a triplet
(, k, `) ⊆ V such that (j, k), (`, k) ∈ E but j and ` are not connected in either direction. This motivates
the representation of a Markov equivalence class as a completely partially directed acyclic graph (CPDAG),
3

which is a graph containing both directed and undirected edges. A directed edge means that all DAGs in
the Markov equivalence class share the same direction for this edge whereas an undirected edge means that
both directions for that specific edge are present within the class. In the same way, one can represent a
subset of a Markov equivalence class via a partially directed acyclic graph (PDAG), where the directions
of the edges are only determined by the graphs within the subset. In particular, some undirected edges in
a CPDAG would become directed edges in a PDAG representing a subset of the class. Notice that both
DAGs and CPDAGs are special cases of PDAGs, where the former represents a single graph and the latter
represents the whole equivalence class.

2.2 `0 -penalized maximum likelihood estimation for a single DAG model
We denote by X̂ ∈ Rn×p our observed data, where each row of X̂ represents a realization of the random
vector X. We say that we are in the low-dimensional setting if asymptotically p remains a constant as n →
∞. By contrast, whenever p → ∞ as n → ∞, we say that we are in the high-dimensional setting. Assuming
faithfulness, Chickering [6] shows that GES achieves a consistent estimator in the low-dimensional setting
by optimizing the following objective – also known as the Bayesian information criterion (BIC) –
(Â, Ω̂) := argmax `n (X̂; A, Ω) − λ2 kAk0 ,

(4)

A∈A, Ω∈D+

where λ2 = 21 logn n , A denotes the set of all valid adjacency matrices associated with DAGs, D+ is the set
of non-negative diagonal matrices, and `n is the likelihood function
!

X̂ T X̂
(5)
· (I − A)Ω−1 (I − A)T + log det (I − A)Ω−1 (I − A)T .
`n (X̂; A, Ω) := −trace
n
In the high-dimensional setting, van de Geer and Bühlmann [43] give consistency guarantees for the global
optimum of (4) when the collection A is further constrained to contain only adjacency matrices with at most
d incoming edges for each node, where d = O(n/ log p). More precisely, they show that there exists some
parameter λ2  logn p such that the optimum (Â, Ω̂) in (4) converges in Frobenius norm to (A0π̂ , Ω0π̂ ) for
increasing n and p, where π̂ is a permutation consistent with Â, i.e.,

kÂ − A0π̂ k2F + kΩ̂ − Ω0π̂ k2F = O λ2 |G0 | .

(6)

Notice, however, that (6) does not guarantee statistical consistency since π̂ need not be a permutation consistent with the true underlying DAG. Moreover, (6) does not hold for any permutation π̂ consistent with Â,
but [43] shows the existence of at least one such permutation. In addition, it is shown in [43] that the number
of non-zero elements in Â, A0π̂ , and A0 are all of the same order of magnitude, i.e., |Ĝ|  |G0π̂ |  |G0 |.

2.3

Collection of DAGs

Consider the setting where not all the observed data comes from the same DAG, but rather from a collection
of DAGs {G (k) = (V, E (k) )}K
k=1 that share the same node set V = [p]. In addition, we assume that all
DAGs in a collection are consistent with some permutation π. This precludes a scenario where (i, j) ∈ E (k)
0
and (j, i) ∈ E (k ) for some k 6= k 0 . This is a reasonable assumption in, e.g., the analysis of gene expression
data, where regulatory links may appear or disappear, but they in general do not change direction.
4

(k) }K
(k) }K
Denote by {(A(k) , Ω(k) )}K
k=1 a set of SEMs on the K DAGs {G
k=1 and by {X̂
k=1 the data
(k)
generated from each SEM, where we observe nk realizations for each DAG G . In this way, each row of
the data matrix X̂ (k) ∈ Rnk ×p corresponds to a realization of the random vector X (k) defined as
T

X (k) = A(k) X (k) + (k)

with (k) ∼ N (0, Ω(k) ).

Collections of DAGs arise for example naturally when considering data from perfect (also known as
hard) interventions [9]. Consider a non-intervened DAG G with SEM parameters (A, Ω) [cf. (1)]. Then a
perfect intervention on a subset of nodes Ik ⊂ V gives rise to the interventional distribution
T

with Ik ∼ N (0, ΩIk ),

X Ik = AIk X Ik + Ik

where AIijk = 0 if j ∈ Ik and AIijk = Aij otherwise, and the diagonal matrix ΩIk satisfies ΩIiik = Ωii if
i 6∈ Ik [16, 17]. We denote the DAG given by the non-zero entries of AIk by G Ik .
(k)
(k)
(k)
In accordance with the notation introduced in Section 2.1, we denote by G0 and (A0 , Ω0 ) the true
data-generating DAG and SEM parameters for class k, respectively, and by π0 a permutation that consistent
(k)
(k)
(k)
(k)
with A0 for all classes k ∈ [K]. Moreover, we denote by G0π and (A0π , Ω0π ) the DAG and SEM
(k)
parameters obtained from the Cholesky decomposition of the true precision matrix Θ0 when permuted by
(k)
(k)
π. We denote by Σ0 the true covariance matrix of the SEM for class k, i.e., the inverse of Θ0 . Finally,
(k)
(k)
union
we define G0union as the union of all G0 – i.e., an edge appears in G0union if it appears in any G0 – and G0π
(k)
Ik
Ik
as the union of G0π . For interventional data, we use (A0 , Ω0 ) to denote the true SEM parameters after
intervening on targets Ik .

3

Joint estimation of multiple DAGs

We first present a penalized maximum likelihood estimator that is the natural extension of (4) for the case
where a collection of DAGs is being estimated. Since this involves minimizing k·k0 , we then discuss a
greedy approach that alleviates the computational complexity of this estimator.

3.1

Joint `0 -penalized maximum likelihood estimator

With d denoting a pre-specified sparsity level and wk = nk /n indicating the proportion of observed data
from DAG k, we propose the following estimator:
n
o
π̂, {(Â(k) , Ω̂(k) )}K
k=1 :=

argmax
π,{(A(k) ,Ω(k) )}K
k=1

K
X
k=1

wk `nk (X̂

(k)

;A

(k)

(k)

,Ω

)−λ

2

K
X
k=1

|A(k) |

(7)
0

subject to A(k) ∈ Aπ , kA(k) k∞,0 ≤ d, Ω(k) ∈ D+ ∀k,
where Aπ is the set of all adjacency matrices consistent with permutation π and the matrix norm k · k∞,0
computes the maximum `0 -norm across the rows of the argument matrix. The optimization problem in (7)
seeks to maximize a weighted log-likelihood of the observations (where more weight is given to SEMs with
more realizations) penalized by the support of the union of all estimated DAGs. To see why this is true,
P
(k)
(k)
notice that k K
k=1 |A |k0 counts the number of (i, j) entries for which Aij 6= 0 for at least one graph k.
This penalization on the union of estimated DAGs promotes overlap in the supports of the different A(k) .
5

Algorithm 1 JointGES for joint `0 -penalized maximum likelihood estimation of multiple DAGs.
Input: Collection of observations X̂ (1) ∈ Rn1 ×p , · · · , X̂ (K) ∈ RnK ×p , sparsity constraint d, penalization
parameters λ1 and λ2
Output: Collection of weighted adjacency matrices Â(1) , · · · , Â(K)
1: Apply GES to find Ĝ union , an approximate solution to the following optimization problem
#
!
"
p
K


X
X
(k)
(k)
argmin
wk
min log kX̂j − X̂Paj (G) ak22 + λ21 |Paj (G)|
|Paj (G)|
G
a∈R
(8)
j=1 k=1
subject to
2:

max |Paj (G)| ≤ d
j

union by solving Kp sparse reEstimate the weighted adjacency matrices {Â(k) }K
k=1 consistent with Ĝ
gressions of the form
(k)

âj

=

1
(k)
kX̂j − X̂ (k) ak22 + λ22 kak1 .
n
union
k
a | supp(a)⊆Paj (Ĝ
)
argmin

Regarding the constraints in (7), the first constraint imposes that all estimated DAGs are consistent with
the same permutation π, which is itself an optimization variable. This constraint is in accordance with our
assumption in Section 2.3 and drastically reduces the search space of DAGs. The second constraint ensures
that the maximum in-degree in all graphs is at most d, and the last constraint imposes the natural requirement
that all noise covariances are diagonal and non-negative.
Notice that (7) is a natural extension of (4). Indeed, for the case K = 1 the objective in (7) immediately
boils down to that in (4). Moreover, when there is only one graph and π can be selected freely, the constraint
A(1) ∈ Aπ is effectively identical to A(1) ∈ A, i.e., the constraint in (4). Finally, observe that in (7)
we have included the additional maximum in-degree constraint required in the high-dimensional setting
[cf. discussion after (5)].

3.2

JointGES: Joint greedy equivalence search

The `0 norms as well as the optimization over all permutations π render the problem of (7) non-convex, thus,
hard to solve efficiently. In this section, we present a greedy approach to find a computationally tractable
approximation to a solution to (7). The algorithm, which we term JointGES, is succinctly presented in
Algorithm 1 and consists of two steps.
In the first step of Algorithm 1 we recover Ĝ union , our estimate of the union of all the DAGs to be
inferred. We do this by finding an approximate solution to (8) via the implementation of GES [6]. The
objective (scoring function) in (8) consists of two terms. The first term is given by the sum of the log(k)
likelihoods of the achievable residues when regressing the jth column of X (k) , denominated as Xj , on
(k)

XPaj (G) for each node j and DAG k. In [43], van de Geer and Bühlmann further show that if we keep
the underlying DAG G fixed,
estimator proposed in (4) is equivalent to optimizing
 the maximum likelihood

Pp 
2
. Thus, the first term in (8) corresponds to the first term
j=1 mina∈R|Paj (G)| log kX̂j − X̂Paj (G) ak2
in the objective of (7). The second term penalizes the size of the parent set of each node in the graph to
be recovered, effectively penalizing the number of edges in the graph. In this way, the scoring function in
6

(8) promotes a sparse G in the same way that the objective of (7) promotes the union of all K recovered
graphs to have a sparse support. Additionally, it is immediate to see that the scoring function in (8) is
decomposable [6], a key feature that enables the implementation of GES to find an approximate solution.
Once we have obtained the union of all sought DAGs Ĝ union from step 1, in the second step of our algorithm
we estimate the DAGs Ĝ (1) , · · · , Ĝ (K) by searching over the subDAGs of Ĝ union . More precisely, for each
(k)
(k)
node j we estimate its parents in Ĝ (k) by regressing Xj on X
using lasso, where the support of
union
Paj (Ĝ

(k)
âj

)

Ĝ (k) .

corresponds to the set of parents of j in
To summarize, Algorithm 1 recovers K DAGs by first estimating the union of all these DAGs Ĝ union
using standard GES and then inferring the specific weight adjacency matrices Â(k) via a lasso regression,
while ensuring consistency with the previously estimated Ĝ union .

4

Consistency guarantees

The main goal of this section is to provide theoretical guarantees on the consistency of the solution to
Problem (7) in the high-dimensional setting. Our main result is presented in Theorem 4.9. In Section 4.3
we present a laxer statement of consistency based on milder conditions.

4.1

Statistical consistency of the joint `0 -penalized MLE

A series of conditions must hold for our main result to be valid. We begin by stating these conditions
followed by the formal consistency result in Theorem 4.9. The rationale behind these conditions and their
implications are discussed after the theorem in Section 4.2.
(1)

(K)

Condition 4.1. All DAGs G0 , · · · , G0

are minimal-edge I-MAPs.

Condition 4.2. There exists a constant σ02 that bounds the variance of all the observed processes, i.e.,
(k)
maxk,i [Σ0 ]ii ≤ σ02 .
(k)

(k)

Condition 4.3. The smallest eigenvalues of all Σ0 are non-zero, i.e. mink Λmin (Σ0 ) = Λmin > 0.
Condition 4.4. There exists some constant α such that, for all k, the maximum allowable in-degree d is
bounded as d ≤ αnk / log p.
Condition 4.5. For all π and j there exist some constants α̃ and cs > 2 such that
)!
(
1
3
n
n
union
, 7
.
|Paj (G0π )| + cs ≤ α̃ min
K 7 (log p)3
K (log n)2 log p
Condition 4.6. The number of DAGs K satisfies K = o(log p) and the amount of data associated with each
DAG is comparable in the sense that n1  n2  · · ·  nK .
PK
(k)
union | ≤ c
Condition 4.7. There exists some constant ct > 0 such that |G0π
t
k=1 wk |G0π | for any permutation π.
Condition 4.8. There exist constants η1 and η2 such that 0 ≤ η1 < 1, 0 < η02 < (1 − η1 )/ct , and
(
p
q
)
X
log
p/n
(k)
(k)
1 [A0π ]i,j >
p/|G0union | ∨ 1
≥ (1 − η1 )|G0π |,
η0
i,j

7

(9)

for all permutations π and graphs k ∈ [K], where 1{·} denotes the indicator function and ct is as in
Condition 4.7.
With the above conditions in place, the following result can be shown.
Theorem 4.9. If Conditions 4.1-4.8 hold and λ is chosen such that λ2 

log p
n



p
|G0union |


∨ 1 , then there exists

a constant c > 0 , that depends on cs , such that with probability 1 − exp(−cp) the solution to (7) satisfies
K
X

(k)

wk kÂ

−

(k)
A0π̂ k2F

+

k=1

K
X


(k)
wk kΩ̂(k) − Ω0π̂ k2F = O λ2 |G0union | .

(10)

k=1

Furthermore, denoting by Ĝ the union of the graphs Ĝ (k) associated with the K recovered adjacency matrices Â(k) , we have that
union
|Ĝ|  |G0π̂
|  |G0union |.

(11)

The proof of Theorem 4.9 is given in Appendix A.2. To intuitively grasp the result in the above theorem,
assume that the number of edges in G0union is proportional to the number of nodes p so that λ2 |G0union | → 0 for
increasing n as long as n > p log p. Hence, under these conditions, (10) guarantees that for the recovered
(k)
permutation π̂, the estimated adjacency matrix Â(k) converges to A0π̂ in Frobenius norm for all k. This not
only implies that both adjacency matrices have similar structure, but also that the edge weights are similar.
Moreover, from (11) it follows that the number of edges in the estimated graph Ĝ, i.e., |Ĝ| is similar to the
union |. More importantly,
number of edges in the union of all minimal I-MAPs with permutation π̂, i.e., |G0π̂
union
|Ĝ| is also similar to the number of edges in the true union graph |G0 |. Despite these guarantees, it should
be noted that the permutation π̂ need not coincide with the permutation π of the true graphs to be recovered.
We now assess the benefits of performing joint estimation of the K DAGs as opposed to estimating
them separately. To do so, we compare the guarantees in Theorem 4.9 to those developed in [43] for
separate estimation. More precisely, the application of the consistency bound reviewed in (6) yields that for
the separate estimation of K DAGs, when we are in the
 setting where
 all K DAGs are highly overlapping
log p
p
2
(cf. Condition 4.7), by choosing λ such that λ  n
∨ 1 , one can guarantee that
|G union |
0

K
X

wk kÂ

(k)

−

(k)
A0π̂(k) kF

k=1

+

K
X


(k)
wk kΩ̂(k) − Ω0π̂(k) k2F = O Kλ2 |G0union | ,

(12)

k=1

where it should be noted that in the separate estimation the recovered permutation π̂ (k) can vary with k.
A direct comparison of (10) and (12) reveals that performing joint estimation improves the accuracy by
a factor of K from Ω(K logn p ) to Ω( logn p ). Hence, for joint estimation the accuracy scales with the total
number of samples n, showing that our procedure yields maximal gain from each observation, even if the
data is generated from K different DAGs. Moreover, the result in (10) holds under slightly milder conditions
than those needed for (12) to hold since Condition 4.8 is a relaxed version of the beta-min condition in [43].
A more detailed discussion about the conditions of Theorem 4.9 is given next.

4.2

Conditions for Theorem 4.9

It has been shown in [33] that if a data-generating distribution is faithful with respect to G, then G must
be a minimal-edge I-MAP. By enforcing the latter for every true graph, Condition 4.1 imposes a milder
8

requirement compared to the well-established faithfulness assumption [39]. Conditions 4.2-4.4 ensure that
we avoid overfitting and provide bounds for the noise variances. These are direct adaptations from Conditions 3.1-3.3 in [43]. Condition 4.5 is required to bound the difference between the sample variances of our
observations and the true variances, and is related to Condition 3.4 in [43] butadapted to our joint inference
n1/3
setting. Notice that Condition 4.5 is trivially satisfied when p = O K 7/3
. Condition 4.6 follows from
log n
the bounds for sample variances shown in [3]. Intuitively, we are imposing the natural restriction that the
number of DAGs is small compared to the number of nodes p in each DAG and the total number of observations n. Moreover, given that our objective is to draw estimation power from the joint inference of multiple
graphs, we require that each DAG is associated with a non-vanishing fraction of the total observations.
Condition 4.7 enforces that, for every permutation π, the number of edges in the union of all recovered
graphs is proportional to the weighted sum of the edges in every graph as K → ∞. In particular, this requires
P
(k)
(k)
the individual graphs G0π to be highly overlapping. To see why this is the case, notice that K
k=1 wk |G0π | is
(k)
upper bounded by the maximum number of edges across graphs G0π . Consequently, Condition 4.7 enforces
the number of edges in the union of graphs to be proportional to the number of edges in the single graph
with most edges, thus requiring a high level of overlap. Imposing high overlap for all permutations π might
seem too restrictive in some settings. Nonetheless, Condition 4.7 can sometimes be derived from apparently
less restrictive conditions as the following example illustrates.
P
(k)
Consider the more relaxed bound |G0union | ≤ ct K
k=1 wk |G0 |, which is equivalent to requiring Condition 4.7 to hold but only for permutations consistent with the true graph G0union . In the following example, we
show that this might be sufficient for Condition 4.7 to hold. Suppose that G0union consists of two connected
respectively defined on the subsets of nodes V1 and V2 . Moreover, assume
and G 00 union
components G 0 union
0
0
(k)
(k)
that the subDAGs of G0 over V1 (denoted by G 0 0 ) are identical for all k. Putting it differently, the differ(k)
ences between the DAGs G0 are limited to the second connected component. In addition, assume that for
0 union |. Then, for any permutation π,
all possible permutations π2 of nodes V2 we have that |G 00 union
0π2 | ≤ |G 0
where we denote by π1 (respectively π2 ) the restriction of π to the node set V1 (respectively V2 ), we have
that
union

union

union
|G0π
| = |G 0 0π1 | + |G 00 0π2 | ≤

K
X

(k)

wk |G 0 0π1 | +

K
X
k=1

k=1

(k)

wk |G 0 0 | ≤ 2

K
X

(k)

wk |G0π |,

k=1

from where Condition 4.7 is satisfied for ct = 2. This example shows that learning the structure of large
components that are common across the different DAGs is not affected by the changes in the smaller components of these DAGs. Despite the previous example, Condition 4.7 might still be too restrictive in some
settings. In this respect, a more relaxed requirement is stated in Section 4.3 under Condition 4.7’.
Condition 4.8 requires that, for every permutation π and every graph k, the value of at least a fixed
(k)
proportion (1 − η1 ) of the edges in G0π is above the ‘noise level’, i.e., the lower bound within the indicator
function in (9). Intuitively, if the true weight of many edges is close to zero then correct inference of the
graphs would be impossible since the true edges would be mistaken with spurious ones. Thus, it is expected
that the weights of a sufficiently large fraction of the edges have to be sufficiently large. Condition 4.8 is
the right formalization of this intuition. Moreover, notice that a straightforward replication
of the beta-min
p
condition introduced
in
[43]
would
have
required
the
‘noise
level’
to
scale
with
log
p/n
k , instead of the
p
smaller scaling of log p/n required in (9). In this sense, Condition 4.8 is a relaxed version of the extension
of the beta-min condition to the setting of joint graph estimation.
9

4.3

Consistency under milder conditions

As previously discussed, in some settings Condition 4.7 might be too restrictive. Hence, in this section we
present a consistency statement akin to Theorem 4.9 that holds for a milder version of Condition 4.7. More
precisely, consider the alternative condition stated next.
Condition 4.7’. Let ct (π) be some function of π that scales as a constant for permutations consistent with
PK
(k)
union | ≤ c (π)
G0union and scales as o(K) for all other permutations such that |G0π
t
k=1 wk |G0π | for all π.
Observe that for permutations π consistent with the true union graph G0union , Condition 4.7’ boils down
to the previously discussed Condition 4.7. However, for all other permutations, ct (π) need not be a constant
and is allowed to grow with K. Intuitively, for all permutations not consistent with G0union we are not
(k)
requiring a high level of overlap among all the graphs G0π . Nonetheless, since ct (π) = o(K) we do require
union to be ‘sparser’ than the extreme case in which all graphs G (k) are disjoint.
G0π
0π
In order to account for the fact that ct depends on the permutation π in Condition 4.7’, we have to modify
Condition 4.8 accordingly, resulting in the following alternative statement.
Condition 4.8’. Let Cmax := max ct (π), then there exist constants η1 and η2 such that 0 ≤ η1 < 1,
π

0 < η02 < (1 − η1 ), and
(
p
)
q
X
Cmax log p/n
(k)
(k)
≥ (1 − η1 )|G0π |,
1 [A0π ]i,j >
p/|G0union | ∨ 1
η0
i,j

for all permutations π and graphs k, where 1{·} denotes the indicator function.
The following consistency result holds for the alternative set of conditions.
Theorem 4.10. Under Conditions 4.1-4.6, 4.7’, and 4.8’ and with λ such that λ2  Cmax logn p



p
|G0union |


∨1 ,

then there exists a constant c > 0 that depends on cs such that with probability 1 − exp(−cp), the solution
to (7) satisfies that, at least for one k ∈ [K],


(k)
(k)
(k)
kÂ(k) − A0π̂ k2F + kΩ̂(k) − Ω0π̂ k2F = O λ2 |G0 | .
(13)
Furthermore, denoting by Ĝ (k) the graph associated with Â(k) for the k ∈ [K] satisfying (13), we have that
(k)

(k)

|Ĝ (k) |  |G0π̂ |  |G0 |.

(14)

The proof is given in Appendix A.3. Condition 4.7’ is milder than Condition 4.7 and this relaxation
entails a corresponding loss in the guarantees of recovery: Comparing (13) and (14) with (10) and (11)
immediately reveals that what could be guaranteed for the ensemble of graphs in Theorem 4.9 can only be
guaranteed for a single graph in Theorem 4.10, thereby explaining the trade-off in relaxing the conditions.
On the other hand, the result in Theorem 4.10 still draws inference power from the joint estimation of
multiple graphs since neither (13) nor (14) can be shown using existing results for separate estimation. To
be more precise, as discussed in Section 4.2, when performing separate estimation, theoretical guarantees
are based on thep
assumption that at least a fixed proportion of the edge weights are above the
p ‘noise level’,
which scales as log p/nk . However, Condition 4.8’ requires the noise level to scale with Cmax log p/n
which, given the fact that Cmax = o(K), is not large enough to achieve the guarantee needed for separate
estimation. In addition, the convergence rate of Ω(Cmax logn p ) in (13) is still faster than the corresponding
10

convergence rate of Ω(K logn p ) associated with separate estimation [cf. discussion after (12)]. We close
Section 4 with the following remark on the consistency guarantees of jointGES.
P
log nk
Remark 4.11 (Consistency of jointGES). In the low-dimensional setting, by choosing λ21 = K
k=1 wk 2nk
and assuming faithfulness, it follows from [6] that the first step in Algorithm 1 is guaranteed to recover the
true Markov equivalence class of G0union , denoted by M(G0union ), in the limit of large data. This allows us
(k)
union ) as inputs to the second step of
to recover {G0 }K
k=1 by successively considering all DAGs in M(G0
union
union
Algorithm 1, and then selecting the DAG G
∈ M(G0 ) whose output from step 2 is the sparsest.
For the high-dimensional setting, as even the global optimum of (7) is not guaranteed to recover the true
G0union (cf. Theorems 4.9 and 4.10), jointGES is not consistent in general. Recently, Maathuis et al. [28]
showed consistency of GES for single DAG estimation under more restrictive assumptions than the ones
here considered. Although of potential interest, further strengthening the presented conditions to guarantee
consistency of jointGES in the high-dimensional setting is not pursued in the current paper.

5

Extension to interventions

We show how our proposed method for joint estimation can be extended to learn DAGs from interventional
data. It is natural to consider learning from interventional data as a special case of joint estimation since
the DAGs associated with interventions are different but closely related. In this section we mimic some of
the developments of Sections 3 and 4 but specialized for the case of interventional data. More precisely,
we first propose an optimization problem akin to (7) and then state the consistency guarantees in the highdimensional setting of the associated optimal solution.
Recall from Section 2.3 that the true adjacency matrix AI0k of the SEM associated with an intervention
in the nodes Ik is identical to the true adjacency A0 of the non-intervened model except that [AI0k ]ij = 0
for all j ∈ Ik . In this way, our assumption that there exists a common permutation π consistent with all
DAGs under consideration (cf. Section 2.3) is automatically satisfied for interventional data. Additionally,
assuming that we observe samples X̂ Ik from K different models corresponding to the respective intervention
of the nodes in {Ik }K
k=1 , the knowledge of the intervened nodes can be incorporated into our optimization
problem as follows [cf. (7)].
o
n
π̂, Â, Ω̂, {(ÂIk , Ω̂Ik )}K
k=1 =

K
X

argmax

π,A,Ω,{(AIk ,ΩIk )}K
k=1

subject to

wk `nk (X̂ Ik ; AIk , ΩIk ) − λ2 kAk0

(15a)

k=1

A ∈ Aπ , kAk∞,0 ≤ d, Ω ∈ D+ ,

(15b)

AIijk
ΩIjjk

= 0 ∀j ∈ Ik ,

(15c)

∈ D+ .

(15d)

= Aij ∀j 6∈ Ik ,

AIijk
Ik

= Ωjj ∀j 6∈ Ik , Ω

From the solution of (15) we obtain an estimate for the non-intervened SEM (Â, Ω̂) as well as K estimates
for the corresponding intervened models (ÂIk , Ω̂Ik ). The objective in (15a) is equivalent to that in (7) where
we leverage the fact that the union of all intervened graphs results in the non-intervened one under the
implicit assumption that no single node has been intervened in every experiment. Alternatively, if some
nodes were intervened in all experiments, objective (15a) would still be valid since enforcing zeros in the
unobservable portions of A does not affect the recovery of the intervened adjacency matrices AIk . The
constraints in (15b) impose that A has to be consistent with permutation π and with bounded in-degree, and
Ω has to be a valid covariance matrix for uncorrelated noise. Putting it differently, (15b) enforces for the
11

non-intervened SEM what we impose separately for all SEMs in (7). The constraints in (15c) impose the
known relations between the intervened and the non-intervened adjacency matrices. Finally, (15d) constrains
the matrices ΩIk to be consistent with the base model on the non-intervened nodes while still being a valid
covariance on the intervened ones.
Even though it might seem that in (15) we are estimating K + 1 SEMs (the base case plus the K
intervened ones), from the previous reasoning it follows that the effective number of optimization variables
is significantly smaller. To be more specific, for a given π, once A is fixed then all the adjacency matrices AIk
are completely determined. Moreover, for a fixed Ω, the only freedom in ΩIk corresponds to the diagonal
entries associated with intervened nodes in Ik . In this way, it is expected that for a given number of samples,
the joint estimation of K SEMs obtained from interventional data [cf. (15)] outperforms the corresponding
estimation from purely observational data [cf. (7)].
Recalling that we denote by (AI0π̂k , ΩI0π̂k ) the parameters recovered from the Cholesky decomposition of
the true precision matrix ΘI0k under the assumption of consistency with permutation π̂, the following result
holds.


Corollary 5.1. If Conditions 4.1-4.8 hold and λ is chosen as λ2  logn p |Gp0 | ∨ 1 , then there exist constants
c1 , c2 > 0 such that with probability 1 − c1 exp(−c2 p), the solution to (15) satisfies
K
X

wk kÂIk − AI0π̂k k2F = O(λ2 |G0 |).

(16)

k=1

Furthermore, denoting by Ĝ the graph associated with the recovered adjacency matrix Â for the nonintervened model, we have that
|Ĝ|  |G0π̂ |  |G0 |.

(17)

The proof is given in Appendix A.4. A quick comparison of Theorem 4.9 and Corollary 5.1 seems to
indicate that the consistency guarantees of observational and interventional data are very similar. Nonetheless, recovery from interventional data is strictly better as we argue next. As discussed after Theorem 4.9,
the results presented do not guarantee that the permutation π̂ recovered coincides with the true permutation
of the nodes. In principle, one could recover a spurious permutation π̂ (different from the true permutation
π) that correctly explains the observed data [cf. (10) and (16)] and leads to sparse graphs [cf. (11) and (17)].
However, the more interventions we have, the smaller the set of spurious permutations π̂ that can be recovered, as we illustrate in the following example. Figure 1 portrays the existence of a spurious permutation
that could be recovered from observational data but cannot be recovered from interventional data. More
precisely, Figure 1(a) presents the two true DAGs we aim to recover, where the second one is obtained
by intervening on node 2. By contrast, Figure 1(b) shows the DAGs that are obtained when performing
Cholesky decompositions on the true precision matrices under the spurious permutation π1 . Notice that the
sparsity levels of the DAGs in both figures are the same. In general, one could recover π1 instead of π0 from
observational data, but one would never recover π1 from interventional data. To see this, simply notice from
the figure that [AI0π2 1 ]32 6= 0 whereas for the interventional estimate [ÂI2 ]32 = 0 [cf. (15c)], thus, the error
terms in (16) cannot vanish for π̂ = π1 . This example also indicates that it is preferable to intervene on
multiple targets in the same experiment instead of doing interventions one at a time. This observation is in
accordance with new genetic perturbation techniques, such as Perturb-seq [8].
From a practical perspective, the objective in (15) corresponds to the same scoring function as GIES [16].
Therefore, GIES can be used to obtain an approximate solution to (15). A simulation study using GIES was
12

π0 = (1, 2, 3)

I1 = ;

I 2 = {2}

1

I

A 01 :

π1 = (3, 2, 1)

I1 = ;

3

1

I

A 02 :

2
3

I 2 = {2}

I

1
A 0π
:
1

2
3

3

I

2
A 0π
:
1

2
1

(a)

2
1

(b)

Figure 1: Interventional data can avoid the recovery of spurious permutations. (a) True DAGs to be recovered. (b) DAGs obtained from the Cholesky decomposition consistent with π1 . The spurious permutation
π1 does not satisfy (16) for cases where node 2 is intervened.
performed in [16, Section 5.2] showing that in line with the theoretical results obtained in this section, not
only identifiability increases, but also the estimates obtained using interventional data are better than with
the same amount of purely observational data.

6

Experiments

In this section, we present numerical experiments on both synthetic (Section 6.1) and real (Section 6.2) data
that support our theoretical findings.

6.1

Performance evaluation of joint causal inference

We analyze the performance of the joint recovery of K different DAGs where we vary K ∈ {3, 5}. For
all experiments, we set the number of nodes p = 100 and the total number of observations n = 600. In
addition, we selected the number of samples from each DAG to be the same, i.e., n1 = . . . = nK = n/K.
For each experiment, the true DAGs were constructed in two steps. First, we generated a core graph that is
shared among the K DAGs under consideration. We did this by generating a random graph from an ErdősRényi model with 100 edges in expectation, and then oriented the edges according to a random permutation
of the nodes. In the second step, we sampled 30 additional edges specific to each DAG uniformly at random.
(1)
(K)
This procedure results in the generation of a collection of true underlying DAGs G0 , . . . , G0 . Associated
(k)
(k)
with each DAG, we generated a true adjacency matrix A0 and a true diagonal error covariance Ω0 . For
the latter, we drew each error variance independently and uniformly from the interval [1, 2.25]. Regarding
the adjacency matrices, we drew the edge weights independently and uniformly from [−1, −0.1] ∪ [0.1, 1] to
ensure that they are bounded away from zero. Notice that we did not put any constraints on the edge weights
that are in the shared core structure for different DAGs: the same edge can change its weight in different
DAGs, or even flip sign.
For each setting, we randomly generated 100 collections of DAGs and data associated with them. We
then estimated the DAGs from the data via two different methods: a joint estimation procedure using the
jointGES algorithm presented in Section 3.2 and a separate estimation procedure using the well-established
GES method [6].
To assess performance of the two algorithms, we considered two standard measures, namely the structural Hamming distance (SHD) [42] and the receiver operating characteristic (ROC) curve. SHD is a com13

(a) Average SHD

(b) Average ROC curve

Figure 2: (a): Average SHD as a function of the scaling constant c for joint and separate GES; (b): Average
ROC curve obtained by varying the tuning parameters. JointGES consistently achieves a better performance
for both values of K.
monly used metric based on the number of operations needed to transform the estimated DAG into the true
one [19, 42]. Hence, a smaller SHD value indicates better performance. The ROC curve plots the true
positive rate against the false positive rate for different choices of tuning parameters.
In Figure 2 (a), we selected the `0 -penalization parameter λ21 = c logn p with scaling constant c ∈
{1, 2, 3, 4, 5} in both joint and separate estimation and then plotted average SHD as a function of the scaling
constant c averaged over the K DAGs to be recovered and the 100 realizations. The penalization parameter
λ2 in the second step of the joint estimation procedure was chosen based on 10-fold cross validation. In
Figure 2 (b) we plotted the average ROC curve where for each choice of tuning parameter, we computed the
true positive and false positive rates by averaging over the K DAGs to be recovered and the 100 realizations.
It is clear from the two figures that in general joint inference achieves better performance, which matches
our theoretical results in Section 4. However, Figure 2 (a) shows also that jointGES performs worse than
separate estimation for small scaling constants (c = 1). Note that this is in line with our theoretical findings
in Theorem 4.10, which imply that whenever Condition 4.7 – which sometimes is a restrictive assumption –
is violated, we need to choose a larger penalization parameter.

6.2

Gene regulatory networks of ovarian cancer subtypes

To assess the performance of the proposed joint `0 -penalized maximum likelihood method on real data, we
analyzed gene expression microarray data of patients with ovarian cancer [41]. Patients were divided into six
subtypes of ovarian cancer, labeled as C1-C6, where C1 is characterized by significant differential expression
of genes associated with stromal and immune cell types and with a lower survival rate as compared to the
other 5 subtypes. We hence grouped the subtypes C2-C6 together and our goal was to infer the differences
in terms of gene regulatory networks in ovarian cancer that could explain the different survival rates. The
gene expression data in [41] includes the expression profile of n = 83 patients with C1 subtype and n = 168
patients with other subtypes. We implemented our jointGES algorithm (Algorithm 1) to jointly learn two
gene regulatory networks: one corresponding to the C1 subtype GC1 and another corresponding to the
other five subtypes together GC2−6 . In addition, as in [3], we focused on a particular pathway, namely the
apoptosis pathway. Using the KEGG database [20, 29] we selected the genes in this pathway that were
associated with at most two microarray probes, resulting in a total of p = 76 genes.
14

Method
JointGES
GES
PC

|GC1 |
50
68
14

|GC2−C6 |
73
101
30

|GC1 ∩ GC2−C6 |
48
32
9

Table 1: Number of edges in the DAGs estimated by different methods for the gene regulatory network of
subtype C1 (|GC1 |) and all other subtypes (|GC2−6 |). The last column shows the number of edges shared
between both inferred graphs.
Table 1 lists the number of edges discovered by jointGES as well as for two separate estimation methods, namely using the GES [6] and PC [39] algorithms. All three methods were combined with stability
selection [26] in order to increase robustness of the output and provide a fair comparison. As expected, the
two graphs inferred using jointGES share a significant proportion of edges, whereas the overlap is markedly
smaller for the two separate estimation methods. Interestingly, under all estimation methods the network for
the C1 subtype contains fewer edges than the network of the other subtypes, thereby suggesting that GC1
could lack some important links that are associated with patient survival.
To obtain more insights into the relevance of the obtained networks, we analyzed the inferred hub nodes
in the three networks. For our analysis we defined as hub nodes those nodes, where the sum of the in- and
out-degree was larger than 5 in the union of the two DAGs. For jointGES, this union is given by Ĝ union , the
graph identified in the first step of Algorithm 1. The hub nodes identified by jointGES are CAPN1, CTSD,
LMNB1, CSF2RB, BIRC3. Among these, CAPN1 [12], CTSD [40], LMNB1 [36], and BIRC3 [18] have
been reported as being central to ovarian cancer in the existing literature. In addition, CSF2RB was also
discovered by joint estimation of undirected graphical models on this data set [3]. The hub nodes discovered
by GES are ATF4, BIRC2, CSF2RB, TUBA1C, MAPK3, while PC only discovered the hub node CSF2RB.
While we were not able to validate the relevance of any of these genes for ovarian cancer in the literature,
interestingly, CSF2RB has been identified as a hub node by all methods, thereby suggesting this gene as an
interesting candidate for future genetic intervention experiments.

7

Discussion

In this paper we presented jointGES, an algorithm for the joint estimation of multiple related DAG models
from independent realizations. Joint estimation is of particular interest in applications where data is collected
not from a single DAG, but rather multiple related DAGs, such as gene expression data from different tissues,
cell types or from different interventional experiments. JointGES first estimates the union of DAGs Ĝ union
by applying a greedy search to approximate the joint `0 -penalized maximum likelihood estimator and then
it uses variable selection to discover each DAG as a subDAG of Ĝ union . From an algorithmic perspective,
jointGES is to the best of our knowledge the first method to jointly estimate related DAG models in the
high-dimensional setting. From a theoretical perspective, we presented consistency guarantees on the joint
`0 -penalized maximum likelihood estimator. and showed that the accuracy bound scales with the total
number of samples, rather than the number of samples from each DAG. As a corollary to this result, we
obtained consistency guarantees for `0 -penalized maximum likelihood estimation of a causal graph from a
mix of observational and interventional data. Finally, we validated our results via numerical experiments on
simulated and real-world data, showing that the proposed jointGES algorithm for joint inference outperforms
separate-inference approaches using well-established algorithms such as PC and GES.
The present work serves as a platform for the potential development of multiple future directions. These
15

directions include: i) relaxing the assumption that all DAGs must be consistent with the same underlying
permutation; ii) extending jointGES to the setting where the samples come from K related DAGs but it
is unknown a priori which particular DAG each sample comes from; this is for example of interest in the
analysis of gene expression data from tumors or tissues that consist of a mix of cell types; iii) extending
jointGES to allow for latent confounders.

A

Appendix: Theoretical analysis of statistical consistency results

Throughout the appendix we develop the proofs of Theorems 4.9 and 4.10. To facilitate understanding, we
first provide a high-level explanation of the rationale behind the proof. If we have data generated from a
single DAG and we are given a permutation π consistent with the true DAG a priori, then we can estimate
(Â, Ω̂) by performing p regressions as explained in Section 2.1. By contrast, when the permutation is
unknown and we need to consider all the possible permutations, the total number of regressions to run
increases to p · p!. However, these regressions are not independent and, intuitively, by bounding the noise
level of a subset of these regressions, we can derive bounds for the noise on the other ones. We characterize
the ‘noise level’ of these regressions by analyzing the asymptotic properties of three random events. More
precisely, whenever these events hold – and we show that they hold with high probability –, the noise is small
enough so that the error of the `0 -penalized maximum likelihood estimator can also be bounded. Finally, we
use this upper bound in the error to show that the recovered graph converges to a minimal I-MAP that is as
sparse as the true DAG.
The remainder of this appendix is organized as follows. In Section A.1 we define the three random
events previously mentioned and show that each of them holds with high probability. Section A.2 then
leverages the definition of these events to show Theorem 4.9, our main result. In Section A.3 we prove
Theorem 4.10, which is based on relaxed conditions but similar proof techniques as in the previous section. Finally, Section A.4 fleshes out the proof of Corollary 5.1, our result applicable to the setting for
interventional data.
Throughout the appendix, we use the following notation. Let âj denote the j-th column of Â and ω̂j
denote the j-th diagonal entry of Ω̂. Also, denote by a0jπ and ω0jπ the j-th column of A0π and the j-th
diagonal entry of Ω0π , respectively.

A.1

Random events

Our proofs of Theorems 4.9 and 4.10 are inspired by the main proof in [43] and, thus, also based on a
set of random events. However, the events considered here differ from those in [43] since, as explained in
Section 4.1, a naive application of the guarantees in [43] to the joint estimation scenario does not achieve
optimal learning rates [cf. discussion after (12)]. Intuitively, the rate gain achieved here comes from the
assumption that all the DAGs share the same permutation, allowing us to effectively reduce the size of the
search space.
In our proofs we consider three random events E1 , E2 , and E3 that are respectively stated – along with
proofs showing that they hold with high probability – in Sections A.1.1, A.1.2, and A.1.3.
Random event E1

A.1.1
(k)

(k)

Let jπ ∈ Rn denote the residual when regressing Xj
(k)

Xj

T (k)

(k)

(k)

(k)

on XS with S = {i | π(i) < π(j)}, i.e., jπ :=

− X (k) a0jπ . Similarly, let ˆjπ ∈ Rn denote the regression residual from the sampled data, i.e.,
16

(k)

(k)

ˆjπ := X̂j

(k)

− X̂ (k) a0jπ . Consider a generic set {A(k) }K
k=1 of adjacency matrices consistent with a given
(k)

(k)

permutation π, where the columns of A(k) are denoted by A(k) := (a1 , . . . , ap ), and let G union denote the
union of the support of A(1) , . . . , A(K) . Then, event E1 is defined as
( p K
p X
K
X X wk (k)T
X
2
wk
(k)
(k)
(k)
(k) (k)
E1 := 2
ˆjπ X̂ (aj − a0jπ ) ≤ δ1
X̂ (k) (aj − a0jπ )
n
n
2
(18)
j=1 k=1 k
j=1 k=1 k
o

union
+ λ21 |G union | + |G0π
| /δ1 , ∀ permutations π , and ∀ {A(k) }K
k=1 consistent with π ,
p
for some constant δ1 > 0 and some λ1  (log p)/n. On random event E1 a uniform inequality holds
(k)
across all K DAGs for the sample correlation between the regression residual jπ and any random variable
(k)

T

spanned by the random vector XS , i.e., X (k) v for any v ∈ Rp with vi = 0 for all i 6∈ S. Notice that
(k)
(k)
for convenience for further steps of the analysis, this generic vector v is written as aj − a0jπ in (18).
(k)

Furthermore, for simplicity in the rest of this appendix, we denominate the space spanned by XS as the
(k)
projection space of jπ . Intuitively, one could foresee E1 holding since the expected correlation between the
(k)

regression residual jπ and XS is equal to zero, and therefore the sample correlation can be upper bounded
by a sum of terms that converge to zero as n → ∞ as in (18).
We now show that random event E1 holds with high probability, a result stated in Theorem A.2. Essential
T (k)
(k)
towards proving this result is the observation that, since the random variable X (k) (aj − a0jπ ) lies within
(k)

the projection space of jπ , these two random variables are independent. We can therefore deal with the
(k)

(k)

randomness in ˆjπ and X̂S separately, one at a time. To formally leverage this observation, we rely on
Lemma 7.4 of [43], stated next for completeness.
Lemma A.1. [43, Lemma 7.4] Let Z be a fixed n × m matrix and e1 , · · · , en be independent N (0, σe2 )distributed random variables. Then for all t > 0
!
p
p
P
sup |eT Za|/n ≥ σe ( 2m/n + 2t/n) ≤ exp(−t).
kZak22 /n≤1

Based on the above lemma and recalling that Aπ denotes the set of adjacency matrices consistent with a
given permutation π, we can show the following result.
Theorem A.2. Assume that Conditions 4.2 and 4.6 hold, then for all t > 0 and all δ1 > 0,
P

max
π

sup
{A(k) }K
k=1 ∈Aπ

2

p X
K
X
wk
j=1 k=1

nk

(k)T
(k)
ˆjπ X̂ (k) (aj

≥

−

(k)
a0jπ )

16σ02 (t +

− δ1

p X
K
X
wk

(k)

n
j=1 k=1 k

union |)
2 log p)(|G union | + |G0π
nδ1

(1)T

(K)T

(k)

X̂ (k) (aj − a0jπ )

2
2

≤ exp(−t).
(1)T

(K)T

Proof. Let ˆjπ and a0jπ be the concatenated vectors ˆjπ := (ˆ
jπ , . . . , ˆjπ )T and a0jπ := (a0jπ , . . . , a0jπ )T .
Also, define the block diagonal matrix X̂ := diag(X̂ (1) , · · · , X̂ (K) ). We denote by Ajπ ⊂ RpK the set that
(k)
contains all vectors that can be formed by vertically concatenating the jth columns aj for all k and that
satisfy that
n
o
(k)
(k)
(k)
Ajπ := aj ∈ RpK | ∀i, if ∃k such that ai,j 6= 0, then Xi ⊥
⊥ jπ for all k .
17

Based on this, and recalling that Paj (·) denotes the set of parent nodes of j in the argument graph, we define
the random event Bjπ as
(
Bjπ := ∃aj ∈ Ajπ :

ˆTjπ X̂(aj − a0jπ ) /n

sup

(19)

kX̂(aj −a0jπ )k22 /n≤1

s
2K(|Paj

≥ σ0 

(G union )|

+

union )|)
|Paj (G0π

n

s
2(t +

+

union )| log p
|Paj (G0π

n


+ 2 log p) 
.


union )|) non-zero entries,
Combining the facts that: i) aj − a0jπ may have at most K(|Paj (G union )| + |Paj (G0π
2
and ii) the variance of each element of ˆjπ is upper bounded by σ0 (cf. Condition 4.2), we may apply
Lemma A.1 to show that

union
P(Bjπ ) ≤ exp −t − |Paj (G0π
)| log p − 2 log p .
(20)
union . Putting
As can be seen from (19), event Bjπ depends exclusively on the set of parents of node j in G0π
union
union , then
it differently, if for two permutations π1 and π2 node j has the same set of parents in G0π1 and G0π
2
the random events Bjπ1 and Bjπ2 coincide, since Ajπ , a0jπ and ˆjπ would all be the same for π ∈ {π1 , π2 }.
union
If we denote
 by Πj (m) the set of permutations where node j has exactly m parents in G0π , then there are
p
at most m unique events Bjπ for all π ∈ Πj (m). We therefore obtain that


 
[
p
P
Bjπ  ≤
exp (−t − m log p − 2 log p) ≤ exp (−t − 2 log p) .
m
π∈Πj (m)

Applying a union bound on the events Bjπ across all nodes j and permutations π yields that




p X
p−1
[[
X
[
P
Bjπ  ≤
P
Bjπ  ≤ p2 exp (−t − 2 log p) = exp(−t).
j

π

j=1 m=1

(21)

π∈Πj (m)

Combining (21) and (19) it follows that with probability at least 1 − exp(−t), for all j, π and all aj ∈ Ajπ
we have that
q

q
|ˆ
Tjπ X̂(aj − a0jπ )|
union
union
union
≤ σ0
2K(|Paj (G
)| + |Paj (G0π )|) + 2(t + |Paj (G0π )| log p + 2 log p) .
kX̂(aj − a0jπ )k2
0(k) }K
Based on the collection of adjacency matrices {A(k) }K
k=1 ∈ Aπ we define another collection {A
k=1
(k)
where each column a0 j is given by
( (k)
(k)T
(k)
(k)
aj
if ˆjπ X̂ (k) (aj − a0jπ ) ≥ 0,
0 (k)
aj =
(k)
(k)
2a0jπ − aj
otherwise.
(k)

(k)

(k)

(k)

Notice that the positions of the non-zero entries in a0 j − a0jπ coincide with those in aj − a0jπ . By also
using the fact that
kX̂(aj − a0jπ )k22 =

K
X
k=1

(k)

(k)

kX̂ (k) (aj − a0jπ )k22 =

K
X
k=1

18

(k)

(k)

kX̂ (k) (a0 j − a0jπ )k22 = kX̂(a0j − a0jπ )k22 ,

we have that for all j and π with probability at least 1 − exp(−t), it holds that
(k)
(k)T
jπ X̂ (k) (aj
k=1 |ˆ

PK

(k)

− a0jπ )|

(k)T
(k)
jπ X (k) (a0 j −
k=1 |ˆ
kX̂(a0j − a0jπ )k2

PK
=

(k)

a0jπ )|

=

|ˆ
Tjπ X̂(a0j − a0jπ )|

(22)
kX̂(aj − a0jπ )k2
kX̂(a0j − a0jπ )k2
q

q
union
union
union
2K(|Paj (G
)| + |Paj (G0π )|) + 2(t + |Paj (G0π )| log p + 2 log p) .
≤ σ0

In the above expression we may use that ∀ δ1 , a, b > 0, it follows ab ≤ δ1 a2 + b2 /δ1 in order to obtain
2|ˆ
Tjπ X̂(a0j − a0jπ )| ≤ δ1 kX̂(a0j − a0jπ )k22 +
(23)
2
q
q
union )|) +
union )| log p + 2 log p)
/δ1 .
4σ02
2K(|Paj (G union )| + |Paj (G0π
2(t + |Paj (G0π
By combining this with the fact that ∀a, b > 0, (a + b)2 ≤ 2(a2 + b2 ) and the fact that K = o(log p) from
Condition 4.6, we get that
q
2
q
union
union
union
2K(|Paj (G
)| + |Paj (G0π )|) + 2(t + |Paj (G0π )| log p + 2 log p)

(24)

union
≤ 4(t + 2 log p)(|Paj (G union )| + |Paj (G0π
)|).

By replacing (24) back in (23), we obtain that
2|ˆ
Tjπ X̂(a0j − a0jπ )|/n ≤ δ1 kX̂(a0j − a0jπ )k22 /n +

union )|)
16σ02 (t + 2 log p)(|Paj (G union )| + |Paj (G0π
.
nδ1

Rewriting the absolute value in the left hand side as the sum of the corresponding K absolute values and
adding the above inequality for j = 1, . . . , p we get that, with probability at least 1 − exp(−t),
2

p X
K
X

(k)T

(k)

(k)

wk |ˆ
jπ X̂ (k) (aj − a0 jπ )|/nk

j=1 k=1

≤

p 
X

δ1 kX̂(aj −

a0jπ )k22 /n

j=1

union )|)
16σ02 (t + 2 log p)(|Paj (G union )| + |Paj (G0π
+
nδ1

Finally, by noticing that kX̂(aj − a0jπ )k22 /n =
of the theorem, concluding the proof.
A.1.2

PK

k=1 wk kX̂

(k) (a(k) −a(k) )k2 /n
k
j
0jπ 2


.

we recover the statement

Random event E2
(k)

(k)

(k)

Let ω0jπ denote the j-th diagonal entry of Ω0π , E2 holds whenever the empirical variances of all jπ , i.e.,
(k)

(k)

kˆ
jπ k22 /nk are close to the true variances ω0jπ . More formally,
E2 :=


p X
K
X


j=1 k=1

wk

(k)
kˆ
jπ k22 /nk −
(k)
ω0jπ

19

(k) !2
ω0jπ

≤ 4 λ22



union
p + |G0π
|
,


(25)

p
for some λ2  (log p)/n. Mimicking the development for event E1 , we now show that E2 also holds with
high probability. This result is stated in Theorem A.4. Similar to the proof of Theorem A.2, we first consider
the asymptotic property for a particular node j and permutation π, and then leverage this to get a uniform
bound across all permutations and nodes. For this proof, we use the following lemma stated in [3], which
also follows from [45]. After the lemma, we formally state our result.
Lemma A.3. [3, Lemma 2] Suppose X1 , · · · , Xn are K-dimensional random vectors satisfying EXi = 0
and kXi k2 ≤ M for 1 ≤ i ≤ n. We have for any β > 0 and x > β


P
1/2
P(k ni=1 Xi k2 ≥ x) ≤ P kN k2 ≥ (x − β)/λmax + c1 K 5/2 exp(−c2 K −5/2 β/M ),
P
where λmax is the largest eigenvalue of Cov( ni=1 Xi ), N is a K-dimensional standard normal random
vector and c1 , c2 are positive constants.
Theorem A.4. Assume Conditions 4.5 and 4.6 hold, then there exist constants c1 , c2 , c3 > 0 such that
!
(k)
(k) !2
p X
K
union | log p
X
kˆ
jπ k22 /nk − ω0jπ
cs p log p + |G0π
wk
P ∃π :
≥ c1
≤ c2 exp(−c3 log p),
(k)
n
ω
j=1 k=1

0jπ

where cs is the constant defined in Condition 4.5.
(k)

Proof. We begin by analyzing the asymptotic properties of jπ for all k given a fixed permutation π and
node j. More specifically, consider the following random event


(k) 2
(k) !2
K
X

union
kˆ
jπ k2 /nk − ω0jπ
2 log p (|Paj (G0π )| + cs )
Cjπ :=
wk
≥
c
,
(26)
1
(k)


n
ω
k=1

0jπ

(k)

for some positive constant c1 . Following the proof of Theorem 1 in [3], we define ut
!

(k) 2
(k)
√
[ˆ

]
−ω

jπ t
0jπ
 wk
if t ≤ nk ,
(k)
(k)
nk
ω0jπ
ut :=


0
otherwise.
(1)

(K) T
)

as follows

(k)

Denoting by ut = (ut , · · · , ut

the random vector collecting all ut , by definition it follows that
2
(k)
(k) !2
K
n
X
X
kˆ
jπ k22 /nk − ω0jπ
wk
=
ut .
(k)
ω
t=1
k=1
0jπ
2

A straightforward substitution in (26) allows us to rewrite the event Cjπ as
)
( n
union )| + c )
X
log p (|Paj (G0π
s
ut ≥ c1 λn
with λ2n :=
Cjπ :=
.
n
t=1

(27)

2

Notice that in order to apply Lemma A.3 to bound the probability of occurrence of Cjπ , we would need
kut k2 to be smaller than a constant M , which is not true in general. Hence, we split Cjπ into two subevents,
enabling the utilization of Lemma A.3. More precisely, whenever the following random event holds
n
o
(k)
1/2−a
Fa := |ut | ≤ λ−1
K
/n,
∀t,
k
,
n
20

then the `2 norm of ut is smaller than
1−a
kut k2 ≤ Ma := λ−1
/n
n K

where a is a free parameter that will be fixed later in the proof. In detail, we bound the probability of Cjπ
according to
P(Cjπ ) ≤ P(Cjπ |Fa )P(Fa ) + P(¬Fa ),
(28)
where we use Lemma A.3 to bound P(Cjπ |Fa ) and where P(¬Fa ) can be estimated from the chi-squared
tail bound [13].
We first focus on bounding P(Cjπ |Fa ). Let us introduce a new variable ũt obtained by truncating the
(1)
(K)
tail of ut . Formally, recalling that 1{·} denotes the indicator function, we have that ũt := (ũt , · · · , ũt )
where
n
o
h
n
oi
(k)
(k)
(k)
(k)
(k)
1/2−a
1/2−a
ũt := ut 1 |ut | ≤ λ−1
/n − E ut 1 |ut | ≤ λ−1
/n .
n K
n K
Notice that whenever
the
h
n random event Fa holds,
oi then ut and ũt follow the same distribution except for a
(k)
(k)
(k)
(k)
(k)
1/2−a
−1
shift vt := E ut 1 |ut | ≤ λn K
/n . Putting it differently, the distribution of ut − ũt is a
P
(k)
(k)
constant
vt whenever we are on the random event Fa . This implies that k nt=1 ut k2 ≤ n maxt,k |vt | +
Pn
k t=1 ũt k2 .
(k)

Therefore, if we can guarantee that n|vt | = o(1)λn P
for all t and k, then there must exist
P a constant
0 < δ < 1 such that if we are on the random event Fa , k nt=1 ut k2 ≥ c1 λn implies that k nt=1 ũt k2 ≥
(1 − δ)c1 λn . Equivalently, we may write
!
!
n
n
X
X
P(Cjπ |Fa )P(Fa ) ≤ P
ũt ≥ (1 − δ)c1 λn Fa P(Fa ) ≤ P
ũt ≥ (1 − δ)c1 λn , (29)
t=1

t=1

2

2

where the second inequality follows from Bayes’ theorem, and we can bound the last term by applying
(k)
Lemma A.3 since ũt is bounded by definition. We now show that, indeed, n|vt | = o(1)λn for all t and k.
From the cumulative tail bound of a chi-squared random variable with one degree of freedom we have that
√
(k)
(k)
P(ut ≥ l) ≤ exp(−ηnl/ K) for some constant η > 0. Based on this, we can estimate the scale of vt
with respect to p and n as
h
n
oi
h
n
oi

(k)
(k)
(k)
(k)
(k)
1/2−a
−1 1/2−a
−a
|vt | = E ut 1 |ut | ≤ λ−1
K
/n
=
E
u
1
|u
|
>
λ
K
/n
≤ exp −η 0 λ−1
n
n
n K
t
t
(k)

for some 0 < η 0 < η, where the second equality follows from the fact that ut has zero mean. Notice
(k)
(k)
that in the last inequality, ut has been absorbed into the exponential term. As |vt | decays exponentially
(k)
with respect to λ−1
n , we have that n|vt | = o(1)λn . Having justified this, we may now apply Lemma A.3
to the rightmost term in (29) in order to bound P(Cjπ |Fa )P(Fa ). From the definition of ũt it follows that
P
P
(k)
λmax {Cov( nt=1 ũt )} ≤ λmax {Cov( nt=1 ut )}. Furthermore, since the variables ut are independently
distributed for all t, we have that


!
(k) 2
(k)
n
X
[ˆ

]
−
ω
wk
jπ t
0jπ
(k)
 ≤ c2 /n
var 
var
ut
=
(k)
nk
ω
t=1

0jπ

21

P
for some constant c2 > 0. This also implies that λmax {Cov( nt=1 ũt )} ≤ c2 /n. Applying Lemma A.3,
where we select x = (1 − δ)c1 λn , β = δ 0 x for some arbitrary positive constant 0 < δ 0 < 1 and, M =
1−a /n, it follows that
λ−1
n K
!


n
X
2
2
a− 72
2
P
ũt ≥ (1 − δ)c1 λn ≤ exp(−c3 (nλn − K)) + exp
log K − c4 K
nλn ,
(30)
5
t=1

2

for some constants c3 , c4 > 0 that increase if constant c1 is increased. In addition, by choosing a = 7/2,
there must exist a large enough c1 such that c3 , c4 > 1 and therefore
!
n
X
P
(31)
ũt ≥ (1 − δ)c1 λn ≤ exp(−nλ2n ) = exp(−|Paj (G union )| log p − cs log p).
t=1

2

Replacing (31) into (29) gives us the sought exponential bound for the first summand in (28).
We are now left with the task of finding a bound for P(¬Fa ). By relying on the fact that n(1)  · · · 
n(K) (cf. Condition 4.6), we get that (maxk nk ) K ≤ c5 n for some constant c5 and therefore
n
o
(k)
1/2−a
P(¬Fa ) ≤ c5 n
max
P |ut | ≥ λ−1
K
/n
.
(32)
n
1≤k≤K,1≤l≤n

Following this, in order to bound the probability that ¬Fa holds we further rely on the tail bound of the
(k)
chi-squared random variable ut to obtain


−a
−a
P(¬Fa ) ≤ c5 n exp −ηλ−1
= c5 exp log n − ηλ−1
.
n K
n K
Recalling the definition of λn from (27), Condition 4.5 implies that
3
λ−1
n
≥ log p(|Paj (G union )| + cs )/α̃ 2
7/2
K

√
and

α̃

λ−1
n
≥ log n.
K 7/2

(33)

By recalling that we have fixed a = 27 , it follows that there exists a constant η 0 such that

√


−a
−a
−1 −a
−1 −a
≤ c5 exp −η 0 λ−1
,
P(¬Fa ) ≤ c5 exp log n − ηλ−1
K
≤
c
exp
α̃λ
K
−
ηλ
K
5
n K
n
n
n
where we have used the second inequality in (33). Furthermore, by leveraging the first inequality in (33) we
obtain that



3
−a
0
union
2
≤
c
exp
−η
log
p
(|Pa
(G
)|
+
c
)/α̃
P(¬Fa ) ≤ c5 exp −η 0 λ−1
K
,
(34)
5
j
s
n
thus obtaining an exponential bound for the second summand in (28).
Having found exponential bounds for both summands in (28), it follows that for c1 > 0 sufficiently large
and α̃ sufficiently small we have that P(Cjπ ) ≤ (1 + c5 ) exp(−|Paj (G union )| log p − cs log p). Following an
argument based on union bounds similar to the one presented in the proof of Theorem A.2, we have that


(k) 2
(k) !2
K
union
X
kˆ
jπ k2 /nk − ω0jπ
log p (|Paj (G
)| + cs )
P
wk
, ∀ j, π 
≤ c1
(k)
n
ω0jπ
k=1


p [
p
[
[
≥ 1 − P
Cjπ  ≥ 1 − (1 + c5 ) exp(−(cs − 2) log p)
j=1 m=1 π∈Πj (m)

22

for some constant c1 > 0. It is immediately implied from the previous expression that


(k)
(k) !2
p X
K
union
X
kˆ
jπ k22 /nk − ω0 jπ
cs p log p + |G0π | log p 
≥ c1
P∃π :
wk
≤ (1 + c5 ) exp(−(cs − 2) log p),
(k)
n
ω0 jπ
j=1 k=1
thus recovering the statement of the theorem (since cs > 2 by Assumption 4.5) after accordingly renaming
the constants on the right hand side.
A.1.3

Random event E3
(k)

(k)

K
Event E3 is defined as the intersection of 2K events that we denote by {E3a }K
k=1 and {E3b }k=1 , where
(k)
(k)
(k)
events E3a and E3b are specific to the kth SEM. More specifically, event E3a ensures that all the estimated
(k)
noise variances ω̂j associated with the kth SEM are finite and bounded away from zero. Formally, we
define the following events for k = 1, . . . , K,
n


o
(k)
(k)
(k)
E3a := min ω̂j , 1/ω̂j
≥ 1/β 2 , for j = 1, . . . , p ,
(35)
(k)

for some β > 0. Event E3b imposes a universal lower bound on the norm achievable by any linear combinations of the data associated with the kth DAG. Mathematically, we consider the ensuing events for
k = 1, . . . , K,
n


o
p
√
(k)
(k)
E3b := kX̂ (k) vk2 / nk ≥ δ3 − λ3
kvk0 kvk2 , ∀v ∈ Rp ,
(36)
(k)

for some δ3 > 0 and λ3 

p
(k)
(k)
(k)
(log p)/nk . Based on (35) and (36) we define events E3 := E3a ∩ E3b , and
E3 :=

K
\

(k)

E3 .

(37)

k=1
(k)

Leveraging the fact that Condition 4.4 enforces the maximum in-degree of each G0π to be at most αnk / log p
for some positive constant α, we can generalize Lemmas 7.5 and 7.7 from [43] into the following lemma.
Lemma A.5. [43, Lemmas 7.5 and 7.7] Assume Conditions 4.2, 4.3, 4.4 and 4.5 hold and that
r
p
√
2(t + log p)
3 Λmin /4 −
− 3σ0 α + α̃ ≥ 1/β > 0,
n
for some t > 0. Based on this, define
s
s
p
log p
2(t + log p)
(k)
λ3 := 3σ0
, and δ3 := 3 Λmin /4 −
nk
nk
(k)

(k)

Then P(E3 ) ≥ 1 − 5 exp(−t) and on E3

it holds that
√
(k)
(k)
(k)
(k)
kX̂ (k) (aj − a0j π̂ )k2 / nk ≥ kaj − a0j π̂ k2 /β 2 .
(k)

(38)

Lemma A.5 shows that under certain conditions the events E3 hold with high probability, thus playing
a role analogous to that of Theorem A.2 for event E1 and Theorem A.4 for event E2 .
With the events E1 , E2 , and E3 defined and having shown under which conditions these hold with high
probability, in the next section we leverage these events to prove Theorem 4.9.
23

A.2

Proof of Theorem 4.9

A.2.1

Bounds on new probability space

Through direct manipulation of the likelihood function, in Lemma A.6 we show that the global optimum
(k)
(k) K
(Â(k) , Ω̂(k) )K
k=1 converges to the SEMs (A0π̂ , Ω0π̂ )k=1 , where π̂ is some permutation consistent with the
estimated adjacency matrices Â(1) , · · · , Â(K) . Lemma A.6 is inspired by Lemma 7.1 in [43].
Lemma A.6. Assume we are on E1 ∩E2 ∩E3 and Condition 4.2 holds. Consider anregularizer in (7) satisfying
o
λ2 > λ21 /δ1 + λ22 /δ2 with 0 < δ1 < 1/β 2 and 0 < δ2 < 1/(2β 2 σ02 ). Then, π̂, {(Â(k) , Ω̂(k) )}K
k=1 the
global optimum of (7), satisfies
(k)

X
X

p X
p X
K
K
1
1
(k)
(k) (k)
2
− δ1
wk kX (âj − a0j π̂ )k2 /nk +
− δ2
wk
β2
2β 4 σ04
j=1 k=1
j=1 k=1


union |)
λ2 |G union |
λ2 (p + |G0π̂
λ21 λ22
2
+ λ −
−
+ 1 0π̂ .
|Ĝ| ≤ λ2 |G0union | + 2
δ1
δ2
δ2
δ1

ω̂j



(k) !2

− ω0 j
(k)

ω̂j

(39)

Proof. By definition, the global optimum must satisfy
K
X

K
X

wk `nk (X̂ (k) ; Â(k) , Ω̂(k) ) − λ2 |Ĝ| ≥

k=1

(k)

(k)

wk `nk (X̂ (k) ; A0 , Ω0 ) − λ2 |G0union |.

(40)

k=1
(k)

(k)

Let π̂ denote any permutation consistent with all Â(k) . Since the value of the likelihood `nk (X̂ (k) ; A0 , Ω0 )
(k)
is completely determined by the precision matrices {Θ0 }K
k=1 , it then follows that the likelihood functions
(k)
(k)
(k)
(k)
(k)
(k)
`nk (X̂ ; A0 , Ω0 ) and `nk (X̂ ; A0π̂ , Ω0π̂ ) achieve the same value. We therefore replace the former by
the latter in (40) and expand the definition of the likelihood function in (5) to obtain
p+

p X
K
X

(k)
wk log ω̂j

2

+ λ |Ĝ| ≤

j=1 k=1

p X
K
X

(k)

wk

kˆ
j π̂ k22 /nk

j=1 k=1

(k)
ω0j π̂

+

p X
K
X

(k)

wk log ω0j π̂ + λ2 |G0union |.

j=1 k=1

Basic manipulations transform the above expression into the following inequality




(k) 2
(k)
p X
p X
K
K
X
X
kˆ

k
/n
ω̂j
k
j π̂ 2
wk 
− 1 + λ2 |G0union |.
wk log  (k)  + λ2 |Ĝ| ≤
(k)
ω0j π̂
ω0j π̂
j=1 k=1
j=1 k=1
(k)

Since we are on E3 , we have that 1/ω̂j
(k)

(41)

≤ β 2 [cf. (35)]. By combining this with Condition 4.2 we can

(k)

≤ β 2 σ02 . Then, using the Taylor expansion log(1 + x) ≤ x − x2 /(2(1 + t)2 ), for
 (k) 
ω̂j
−1 < x ≤ t, we can further replace log
in (41) to obtain
(k)
further bound ω0j π̂ /ω̂j

ω0j π̂

p X
K
X
j=1 k=1


wk 

(k)

ω̂j

(k)

− ω0j π̂
(k)

ω̂j


2


(k)
(k) 2
p X
K
X
ω0j π̂
kˆ

k
/n
k
1
j π̂ 2
+

− 1 + λ2 |Ĝ| ≤
wk 
− 1 + λ2 |G0union |.
(k)
2β 4 σ04 ω̂ (k)
ω0j π̂
j=1 k=1
j


(42)
24

(k)

Finally, using the fact that X̂j
(k)

ω̂j

(k)

= kX̂j

(k)

(k)

(k)

= ˆj π̂ + X̂ (k) a0j π̂ , we also rewrite ω̂j

(k)

(k)

(k)

as

(k)T

(k)

(k)

(k)

− X̂ (k) âj k22 /nk = kX̂ (k) (âj − a0j π̂ )k22 /nk − 2ˆ
j π̂ X̂ (k) (âj − a0j π̂ )/nk + kˆ
j k22 /nk .

By replacing the above into (42), we get that
p X
K
X

(k)

wk

(k)

kX̂ (k) (âj − a0j π̂ )k22 /nk
(k)
ω̂j

j=1 k=1

≤2

+

p X
K
X

1
2β 4 σ04

p X
K
X

(k)

(k)

− ω0j π̂

(k)
ω̂j

(k)

ˆj π̂ X̂ (k) (âj − a0j π̂ )/nk
(k)

+

ω̂j

j=1 k=1



−

wk 

(k)

ω̂j

j=1 k=1

(k)T

wk



(k)
p X
K
X
kˆ
j π̂ k22 /nk −
wk 
(k)
ω̂j
j=1 k=1

2
 + λ2 |Ĝ|

p X
K
X


wk 

(43)
(k)

kˆ
j π̂ k22 /nk

j=1 k=1

(k)

ω0j π̂


− 1



(k)
ω0j π̂


+ λ2 |G0union |.

In order to further bound the expression in (43), notice that the first summand in the right hand side of the
inequality corresponds to the sum of all empirical correlation coefficients. Leveraging that we are under the
assumption that E1 holds [cf. (18)], we have that
2

p X
K
X
j=1 k=1

(k)T

wk

(k)

(k)

ˆj π̂ X̂ (k) (âj − a0j π̂ )/nk
(k)

ω̂j

≤ δ1

p X
K
X

(k)

(k)

wk kX̂ (k) (âj − a0j π̂ )k22 /nk +

j=1 k=1

λ21
|Ĝ|.
δ1

(44)

In order to bound the second and third terms, we first restate their difference as follows







(k)
(k)
(k)
(k)
(k)
(k)
(k)
p X
p X
K
K
X
X
kˆ
j π̂ k22 /nk − ω0j π̂
kˆ
j π̂ k22 /nk − ω0j π̂ ω0j π̂ − ω̂j
kˆ
j π̂ k22 /nk
 −wk 

.
wk 
wk 
− 1 =
(k)
(k)
(k)
(k)
ω̂
ω
ω
ω̂
j=1 k=1
j=1 k=1
j
j
0j π̂
0j π̂
(45)
P
n
2
2
k=1 |vk | , we further bound (45)
j=1 |uj |

Pn

Pn

Next, by using Cauchy-Schwarz inequality, i.e. | i=1 ui vi |2 ≤
as




(k)
(k)
(k)
p X
p X
K
K
X
X
kˆ
j π̂ k22 /nk − ω0j π̂
kˆ
j π̂ k22 /nk
−
wk 
− 1
wk 
(k)
(k)
ω̂j
ω0j π̂
j=1 k=1
j=1 k=1



2 1/2 
2 1/2
(k)
(k)
(k)
(k)
p X
p X
K
K
X
X
ω0j π̂ − ω̂j
kˆ
j π̂ k22 /nk − ω0j π̂
  
 
wk 
≤
wk 
(k)
(k)
ω0j π̂
ω̂j
j=1 k=1
j=1 k=1

(46)

From the fact that event
q E2 holds [cf. (25)], we can upper bound the first of the two factors in the right
hand side of (46) by 2 λ22 (p + |G0union (π̂)|). Further, relying on the inequality 2ab ≤ a2 /δ2 + δ2 b2 for any
δ2 > 0, it follows that




(k)
(k)
(k)
p X
p X
K
K
X
X
kˆ
j π̂ k22 /nk − ω0j π̂
kˆ
j π̂ k22 /nk
−
wk 
wk 
− 1
(47)
(k)
(k)
ω̂
ω
j=1 k=1
j=1 k=1
j
0j π̂

2
(k)
(k)
p X
K
2
union
X
ω
−
ω̂
λ (p + |G0 (π̂)|)
j 
0j π̂
≤ 2
+ δ2
wk 
.
(k)
δ2
ω̂
j=1 k=1

25

j

By replacing (44) and (47) into (43), we recover (39), as we wanted to show.
From Lemma A.6 it follows that the global optimum of (7) corresponds to a minimal I-MAP, but no
claim is made about the sparsity level of this I-MAP. In order to show that the solution is indeed sparse, we
union | cannot be much
must rely in Conditions 4.7 and 4.8. Consequently, in Lemma A.7 we show that |G0π̂
union
larger than |Ĝ|. Then, in Thm. A.8 we further show how to cancel out |G0π̂ | with |Ĝ| in (39) to obtain our
main result.
Lemma A.7. Assume Condition 4.7 holds and let λ̃ > 0 be such that
K
X

(k)

union
|.
wk kÂ(k) − A0π̂ k2F ≤ λ̃2 |G0π̂

(48)

k=1

Consider also constants η1 , η2 with 0 ≤ η1 < 1 and 0 < η22 ct < 1−η1 such that

P

i,j 1

n

o
(k)
[A0π̂ ]i,j ≥ λ̃/η2 ≥

(k)

(1 − η1 )|G0π̂ |. Then, it follows that
|Ĝ| ≥

1 − η1 − η22 ct union
|G0π̂ |.
ct

(49)
(k)

Proof. Let N (k) and M(k) be the sets of entries satisfying N (k) := {(i, j) : |[A0π̂ ]i,j | ≥ λ̃/η2 } and
(k)
M(k) := {(i, j) : |[Â(k) ]i,j − [A0π̂ ]i,j | ≥ λ̃/η2 }. From these definitions it follows that
K
X

K

wk |N (k) ∩ M(k) |

k=1

λ̃2 X
wk
≤
η22
k=1

(k)

X

|[Â(k) ]i,j − [A0π̂ ]i,j |2 ≤

K
X

(k)

wk kÂ(k) − A0π̂ k2F .

k=1

(i,j) ∈ N (k) ∩M(k)

Leveraging inequality (48) and Condition 4.7 we further have that
K
X

wk |N

(k)

∩M

(k)

|≤

union
η22 |G0π̂
|

k=1

≤

η22 ct

K
X

(k)

wk |G0π̂ |.

(50)

k=1
C

Notice that for all (i, j)-th entries in the set N (k) ∩ M(k) it must be that |[Â(k) ]i,j | > 0. Hence, N (k) ∩
C

M(k) corresponds to a subset of non-zero entries of Â(k) , which in turn corresponds to a subset of edges
in Ĝ. From this we can infer that
|Ĝ| =

K
X
k=1

wk |Ĝ| ≥

K
X

C

wk |N (k) ∩M(k) | =

k=1

K
X

wk (|N (k) | − |N (k) ∩ M(k) |) ≥ (1 − η1 − η22 ct )

k=1

K
X

(k)

wk |G0π̂ |,

k=1

where the last inequality follows by combining (50) with the definition of η1 in the statement of the lemma.
The proof concludes by replacing Condition 4.7 in the above inequality.
Theorem A.8. Assume Conditions 4.1, 4.7 and 4.8 hold, and suppose that there exist some constants δB
and 0 < δs < 1 as well as λ and λ0 that scale as λ2  λ20  logn p (p/|G0union | ∨ 1) such that
δB

K
X

(k)

union
wk kÂ(k) − A0π̂ k2F + λ2 δs |Ĝ| ≤ λ2 |G0union | + λ20 |G0π̂
|.

k=1

26

(51)

If the constant η0 in Condition 4.8 is sufficiently small, then there exist some constants δs0 , cg , c0g > 0 such
that
δB

K
X

(k)

wk kÂ(k) − A0π̂ k2F + (λ2 δs − λ20 δs0 )|Ĝ| ≤ λ2 |G0union |

(52)

k=1

and
union
| ≥ c0g |G0union |.
|Ĝ| ≥ cg |G0π̂

(53)
(k)

union |. Suppose that for
Proof. Using Conditions 4.1 and 4.7, we have that |G0union | ≤ ct max |G0 | ≤ ct |G0π̂
k

some λ̃ > 0 one has that λ̃2  λ2  λ20 and λ̃2 δB ≥ λ2 ct + λ20 , then it follows from (51) that
δB

K
X

(k)

union
union
|.
| ≤ λ̃2 δB |G0π̂
wk kÂ(k) − A0π̂ k2F ≤ λ2 |G0union | + λ20 |G0π̂

k=1

q
Let η2 be some constant defined as η2 := η0 λ̃/ logn p (p/|G0union | ∨ 1), then we can rewrite Condition 4.8 as
that
o
X n (k)
(k)
1 [A0π̂ ]i,j ≥ λ̃/η2 ≥ (1 − η1 )|G0π̂ |.
i,j

Moreover, for η0 sufficiently small, η2 is also guaranteed that 0 < η22 ct < 1 − η1 . We could therefore
apply Lemma A.7 and get that |Ĝ| ≥
1−η1 −η22 ct
ct

1−η1 −η22 ct union
|G0π̂ |,
ct

which completes the proof of (53) by choosing

c0g

cg =
and = cg · ct . Notice that in order to apply Lemma A.7, it is required that η22 ct < 1 − η1
which is guaranteed by the assumption that η0 is sufficiently small. Leveraging the first inequality in (53),
λ20 ct
union | in (51) by
|Ĝ| in order to obtain
we can replace λ20 |G0π̂
1−η −η 2 c
1

δB

K
X

wk kÂ

(k)

−

(k)
A0π̂ k2F

k=1

2 t

λ20 ct
+ λ δs −
1 − η1 − η22 ct


2



|Ĝ| ≤ λ2 |G0union |.

(54)

Notice that (54) coincides with the sought expression (52) upon substituting δs0 = ct /(1 − η1 − ct η22 ).
A.2.2

Proof of Theorem 4.9

It follows from Theorem A.2, Theorem A.4, and Lemma A.5 that there exist some constants λ1 , λ2 with
(k) 2

p
 log
λ21  λ22  logn p as well as some λ3
nk for all k such that with probability 1 − exp(−c log p) for
some constant c > 0, the random event E1 ∩ E2 ∩ E3 occurs.
We may then apply Lemma A.6 to show that there exist some constants δB , δW such that with high
probability, for any λ > 0 satisfying λ2 > λ21 /δ1 + λ22 /δ2 ,

δB β

2

p X
K
X
j=1 k=1

wk kX

(k)

(k)
(âj

−

(k)
a0j π̂ )k22 /nk

+ δW

K
X

(k)

wk kΩ̂(k) − Ω0π̂ k2F + λ2 δs |Ĝ|

(55)

k=1

≤ λ2 |G0union | +

27

union |)
λ22 (p + |G0π̂
λ2 |G union |
+ 1 0π̂ .
δ2
δ1

Pp



(k)

(k)

ω̂j −ω0 j

2

(k)

Note that compared with Lemma A.6, we have replaced j=1
by kΩ̂(k) − Ω0π̂ k2F . This
(k)
ω̂j
P
(k)
(k)
(k)
1
is
follows from combining the facts that pj=1 (ω̂j − ω0 j )2 is equal to kΩ̂(k) − Ω0π̂ k2F and that (k)
ω̂j

1/β 2

bigger than
on the random event E3 . In addition, the constants β, σ0 , δ1 and δ2 in Lemma A.6 have
λ2
λ2
been absorbed into the new constants δB and δW . We also replaced λ2 − δ11 − δ22 by λ2 δs for some
0 < δs < 1.
By applying Lemma A.5 [cf. (38)] we may bound the first summand on the left hand side of (55)
P
(k) − A(k) k2 . Furthermore, replacing λ and λ by some λ that scales as λ2 
by δB K
1
2
0
0
k=1 wk kÂ
0π̂ F
log p
union | ∨ 1), we obtain that
(p/|G
0
n
δB

K
X

wk kÂ

(k)

−

(k)
A0π̂ k2F

+ δW

k=1

K
X

(k)

union
wk kΩ̂(k) − Ω0π̂ k2F + λ2 δs |Ĝ| ≤ λ2 |G0union | + λ20 |G0π̂
|.

(56)

k=1

By applying (56), we have that for a constant η0 small enough, (52) and (53) in Theorem A.8 holds by
choosing λ as that λ2  logn p (p/|G0union | ∨ 1) and λ2 δs > λ20 δs0 . Moreover, from (52) we further infer that
δB

K
X

(k)

wk kÂ(k) − A0π̂ k2F + λ2 δs00 |Ĝ| ≤ λ2 |G0union |,

(57)

k=1

where the constant δs00 is chosen such that λ2 δs00 = λ2 δs − λ20 δs0 in (52). From (57) it can thus be inferred
union | ≥ c0 /c |G union |, we recover the
that |Ĝ| ≤ |G0union |/δs00 . Combining this with (53) and the fact that |G0π̂
g g 0
union |. For the relation between |G union | and
first part of (11) in the statement of the theorem, i.e., |Ĝ|  |G0π̂
0π̂
union | ≤ |Ĝ|/c ≤ |G union |/(δ 00 · c ) and |G union | ≥ c0 /c |G union |. Finally, to recover
|G0union |, we use that |G0π̂
g
g
g
s
g
0
0
0π̂
(10) we combine (56) with (11).

A.3

Proof of Theorem 4.10

We first introduce a lemma that will be instrumental in proving Theorem 4.10 and that can be obtained by
generalizing Lemmas 7.2 and 7.3 in [43].
(k)

Lemma A.9. [43, Lemmas 7.2 and 7.3] Suppose for some δB , δs , λ0 , λ > 0 one has n
that δB kÂ(k) −A0π̂ ko2F +
P
(k)
(k)
(k)
λ2 δs |Ĝ (k) | ≤ λ2 |G0 | + λ20 |G0π̂ |. Let λ̃2 δB ≥ λ2 + λ20 and assume that i,j 1 [A0π̂ ]i,j ≥ λ̃/η2 ≥
(k)

(1 − η1 )|G0π̂ |. Then

(k)
δB kÂ(k) − A0π̂ k2F + λ2 δs −

λ20
1 − η1 − η22



(k)

|Ĝ (k) | ≤ λ2 |G0 |

and
(k)

(k)

|Ĝ (k) | ≥ (1 − η1 − η22 )|G0π̂ | ≥ (1 − η1 − η22 )|G0 |.
In order to show Theorem 4.10, we begin the proof just like for Theorem
4.9 in Section A.2.2 until we
P
get to expression (56). It then follows from Condition 4.7’ and |Ĝ| ≥ K
w
|Ĝ (k) | that there exists some
k
k=1
28

log p
union | ∨ 1), where C
λ02
max is defined in Condition 4.8, such that for any λ > 0,
0  Cmax n (p/|G0

δB

K
X

wk kÂ

(k)

−

(k)
A0π̂ k2F

+ δW

k=1
K
X

≤ λ2 ct (π0 )

(k)

wk |G0 | +

k=1

λ02

λ2

· ct (π0 ) and

Let
:=
0
any λ > 0,

δs0

K
X

wk kΩ̂

(k)

−

(k)
Ω0π̂ k2F

2

+ λ δs

k=1
K
X
(k)
λ02
wk |G0π̂ |.
0
k=1

K
X

wk |Ĝ (k) |

k=1

:= δs /ct (π0 ), it then follows that there must exist at least one k such that for

(k)

(k)

(k)

(k)

δB kÂ(k) − A0π̂ k2F + δW kΩ̂(k) − Ω0π̂ k2F + λ02 δs0 |Ĝ (k) | ≤ λ02 |G0 | + λ02
0 |G0π̂ |
Since according to Condition 4.7’, ct (π) scales as a constant for permutations consistent with G0union , we have
that δs0 is still a constant and λ0  λ. In this case, it follows from Lemma A.9 and Condition 4.8’ that there
exists some constant 0 < δs < 1 and δs0 > 0 such that by choosing λ0 as that λ0  Cmax logn p (p/|G0union | ∨ 1)
0
and λ02 δs > λ02
0 δs ,
 (k)
(k)
(k)
0
δB kÂ(k) − A0π̂ k2F + λ02 δs − λ02
| ≤ λ02 |G0 |
0 δs |Ĝ
(k)

(k)

It also follows from Lemma A.9 that |Ĝ (k) | ≥ cg |G0π̂ | ≥ cg |G0 | for some positive constant cg . Mimicking
the arguments employed in the proof of Theorem 4.9 from (57) until the end of the proof, one can show that
expressions (13) and (14) in the statement of Theorem 4.10 hold true.

A.4

Proof of Corollary 5.1

The following lemma is instrumental in showing the corollary.
Lemma A.10. [17, Lemma 4] Given fixed G, the maximum likelihood estimator in (15) can be written as




p


X
X
X
n
n
−j
(k)
(k)
(k)
k
 min
log 
kX̂j − X̂Paj (G) · ak22 /nk  +
wk log kX̂j k22 /nk 
p+
|Paj (G)| n
n
−j
a∈R
k:j6∈Ik

j=1

k:j∈Ik

where n−j is the total number of data from the interventions where node j is not intervened, i.e., n−j =
P
k:j6∈Ik nk .
Recall that in the interventional setting G0union is given by the true graph G0 of the non-intervened model,
(k)
(k)
and that the K models to be inferred (A0 , Ω0 ) correspond to the interventional models (AI0k , ΩI0k ). De(k)
noting by (π̂, Â, Ω̂) the (non-intervened) global optimum of (15), let ω̂j denote the empirical variance of
(k)

(k)

the random variable Xj − X (k) âj if j ∈ Ik and the empirical variance of Xj
Lemma A.10 that the global optimum satisfies that




p
X
X
X
nk (k) 
(k)
 n−j log 
p+
ω̂
+
wk log ω̂j  + λ2 |Ĝ|
n
n−j j
j=1

k:j6∈Ik

otherwise. It follows from

k:j∈Ik

≤

p X
K
X

(k)

wk

kˆ
j π̂ k22 /nk

j=1 k=1

29

(k)
ω0j π̂

+

p X
K
X
j=1 k=1

(k)

wk log ω0j π̂ + λ2 |G0union |.

P
P
wk a k ) ≥ K
Then, applying the inequality log( K
k=1 wk log ak for any choices of a1 , . . . , aK > 0
k=1
PK
and w1 , . . . , wK > 0 with k=1 wk = 1, we have that
p+

p X
K
X
j=1 k=1

(k)
wk log ω̂j

2

+ λ |Ĝ| ≤

p X
K
X

(k)

wk

kˆ
j π̂ k22 /nk

j=1 k=1

(k)
ω0j π̂

+

p X
K
X

(k)

wk log ω0j π̂ + λ2 |G0union |.

j=1 k=1

In this case, Corollary 5.1 immediately follows from the proof of Theorem 4.9.

References
[1] P. A. Aguilera, A. Fernández, R. Fernández, R. Rumı́, and A. Salmerón. Bayesian networks in environmental modelling. Environmental Modelling & Software, 26(12):1376–1388, 2011.
[2] M. N. Arbeitman, E. EM. Furlong, F. Imam, E. Johnson, B. H. Null, B. S. Baker, M. A. Krasnow, M. P.
Scott, R. W. Davis, and K. P. White. Gene expression during the life cycle of drosophila melanogaster.
Science, 297(5590):2270–2275, 2002.
[3] T. T. Cai, H. Li, W. Liu, and J. Xie. Joint estimation of multiple high-dimensional precision matrices.
The Annals of Statistics, 38:2118–2144, 2015.
[4] T. T. Cai, W. Liu, and X. Luo. A constrained `1 minimization approach to sparse precision matrix
estimation. Journal of the American Statistical Association, 106(494):594–607, 2011.
[5] D. M. Chickering. Learning Bayesian networks is NP-complete. In Proceedings of the Fifth International Workshop on Artificial Intelligence and Statistics, 1995.
[6] D. M. Chickering. Optimal structure identification with greedy search. Journal of Machine Learning
Research, 3(Nov):507–554, 2002.
[7] P. Danaher, P. Wang, and D. M. Witten. The joint graphical lasso for inverse covariance estimation
across multiple classes. Journal of the Royal Statistical Society: Series B (Statistical Methodology),
76(2):373–397, 2014.
[8] A. Dixit, O. Parnas, B. Li, J. Chen, C. P. Fulco, L. Jerby-Arnon, N. D. Marjanovic, D. Dionne, T. Burks,
R. Raychowdhury, B. Adamson, T. M. Norman, E. S. Lander, J. S. Weissman, N. Friedman, and
A. Regev. Perturb-seq: Dissecting molecular circuits with scalable single-cell RNA profiling of pooled
genetic screens. Cell, 167(7):1853–1866, 2016.
[9] F. Eberhardt, C. Glymour, and R. Scheines. On the number of experiments sufficient and in the worst
case necessary to identify all causal relations among n variables. In Proceedings of the Twenty-First
Conference on Uncertainty in Artificial Intelligence, pages 178–184. AUAI Press, 2005.
[10] J. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covariance estimation with the graphical
lasso. Biostatistics, 9(3):432–441, 2008.
[11] N. Friedman, M. Linial, I. Nachman, and D. Peter. Using Bayesian networks to analyze expression
data. Journal of Computational Biology, 7:601–620, 2000.
30

[12] D. M. Gau, J. L. Lesnock, B. L. Hood, R. Bhargava, M. Sun, K. Darcy, S. Luthra, U. Chandran, T. P.
Conrads, R. P. Edwards, J. L. Kelley, T. C. Krivak, and P. Roy. BRCA1 deficiency in ovarian cancer is
associated with alteration in expression of several key regulators of cell motility–a proteomics study.
Cell Cycle, 14(12):1884–1892, 2015.
[13] B. George. Probability inequalities for the sum of independent random variables. Journal of the
American Statistical Association, 57(297):33–45, 1962.
[14] C. Glymour, R. Scheines, P. Spirtes, and K. Kelly. Discovering Causal Strucure. Academic Press,
1987.
[15] J. Guo, E. Levina, G. Michailidis, and J. Zhu.
Biometrika, 98(1):1–15, 2011.

Joint estimation of multiple graphical models.

[16] A. Hauser and P. Bühlmann. Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs. Journal of Machine Learning Research, 13(Aug):2409–2464,
2012.
[17] A. Hauser and P. Bühlmann. Jointly interventional and observational data: Estimation of interventional
Markov equivalence classes of directed acyclic graphs. Journal of the Royal Statistical Society: Series
B (Statistical Methodology), 77(1):291–318, 2015.
[18] J. Jönsson, K. Bartuma, M. Dominguez-Valentin, K. Harbst, Z. Ketabi, S. Malander, M. Jönsson,
A. Carneiro, A. Måsbäck, G. Jönsson, and M. Nilbert. Distinct gene expression profiles in ovarian
cancer linked to Lynch syndrome. Familial Cancer, 13:537–545, 2014.
[19] M. Kalisch and P. Bühlmann. Estimating high-dimensional directed acyclic graphs with the PCalgorithm. Journal of Machine Learning Research, 8(Mar):613–636, 2007.
[20] M. Kanehisa, S. Goto, Y. Sato, M. Furumichi, and T. Mao. KEGG for integration and interpretation of
large-scale molecular data sets. Nucleic Acids Research, 40(D1):D109–D114, 2011.
[21] M. Kolar, L. Song, A. Ahmed, and E. P. Xing. Estimating time-varying networks. The Annals of
Applied Statistics, pages 94–123, 2010.
[22] S. L. Lauritzen. Graphical Models, volume 17. Clarendon Press, 1996.
[23] P. Loh and P. Bühlmann. High-dimensional learning of linear causal networks via inverse covariance
estimation. Journal of Machine Learning Research, 15(1):3065–3105, 2014.
[24] E. Z. Macosko, A. Basu, R. Satija, J. Nemesh, K. Shekhar, M. Goldman, I. Tirosh, A. R. Bialas,
N. Kamitaki, E. M. Martersteck, J. J. Trombetta, D. A. Weitz, J. R. Sanes, A. K. Shalek, A. Regev, and
S. A. McCarroll. Highly parallel genome-wide expression profiling of individual cells using nanoliter
droplets. Cell, 161(5):1202–1214, 2015.
[25] N. Meinshausen and P. Bühlmann. High-dimensional graphs and variable selection with the lasso. The
Annals of Statistics, pages 1436–1462, 2006.
[26] N. Meinshausen and P. Bühlmann. Stability selection. Journal of the Royal Statistical Society. Series
B: Statistical Methodology, 72(4):417–473, 2010.
31

[27] K. Mohan, P. London, M. Fazel, D. Witten, and S. Lee. Node-based learning of multiple Gaussian
graphical models. The Journal of Machine Learning Research, 15(1):445–488, 2014.
[28] P. Nandy, A. Hauser, and M. H. Maathuis. High-dimensional consistency in score-based and hybrid
structure learning. The Annals of Statistics (to appear), 2018.
[29] H. Ogata, S. Goto, K. Sato, W. Fujibuchi, H. Bono, and M. Kanehisa. KEGG: Kyoto encyclopedia of
genes and genomes. Nucleic Acids Research, 27(1):29–34, 1999.
[30] J. Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, 2000.
[31] J. Pearl and T. S. Vermal. Equivalence and synthesis of causal models. In Proceedings of Sixth
Conference on Uncertainty in Artijicial Intelligence, pages 220–227, 1991.
[32] C. Peterson, F. C. Stingo, and M. Vannucci. Bayesian inference of multiple Gaussian graphical models.
Journal of the American Statistical Association, 110(509):159–174, 2015.
[33] G. Raskutti and C. Uhler. Learning directed acyclic graphs based on sparsest permutations. Stat (to
appear), 2018.
[34] P. Ravikumar, M. J. Wainwright, G. Raskutti, and B. Yu. High-dimensional covariance estimation
by minimizing `1 -penalized log-determinant divergence. Electronic Journal of Statistics, 5:935–980,
2011.
[35] J. M. Robins, M. A. Hernan, and B. Brumback. Marginal structural models and causal inference in
epidemiology. Epidemiology, 11(5):550–560, 2000.
[36] A. D. Santin, F. Zhan, S. Bellone, M. Palmieri, S. Cane, E. Bignotti, S. Anfossi, M. Gokden, D. Dunn,
J. J. Roman, T. J. O’Brien, E. Tian, M. J. Cannon, J. Shaughnessy, and S. Pecorelli. Gene expression profiles in primary ovarian serous papillary tumors and normal ovarian epithelium: identification
of candidate molecular markers for ovarian cancer diagnosis and therapy. International Journal of
Cancer, 112(1):14–25, 2004.
[37] A. K. Shalek, R. Satija, J. Shuga, J. J. Trombetta, D. Gennert, D. Lu, P. Chen, R. S. Gertner, J. T.
Gaublomme, N. Yosef, S. Schwartz, B. Fowler, S. Weaver, J. Wang, X. Wang, R. Ding, R. Raychowdhury, N. Friedman, N. Hacohen, H. Park, A. P. May, and A. Regev. Single cell RNA Seq reveals
dynamic paracrine control of cellular variation. Nature, 510(7505):363, 2014.
[38] L. Song, M. Kolar, and E. P. Xing. Keller: estimating time-varying interactions between genes. Bioinformatics, 25(12):i128–i136, 2009.
[39] P. Spirtes, C. Glymour, and R. Scheines. Causation, Prediction and Search. MIT Press, 2000.
[40] E. A. Stronach, G. C. Sellar, C. Blenkiron, G. J. Rabiasz, K. J. Taylor, E. P. Miller, C. E. Massie,
A. Al-Nafussi, J. F. Smyth, D. J. Porteous, and H. Gabra. Identification of clinically relevant genes
on chromosome 11 in a functional model of ovarian cancer tumor suppression. Cancer Research,
63(24):8648–8655, 2003.
[41] R. W. Tothill, A. V. Tinker, J. George, R. Brown, S. B. Fox, S. Lade, D. S. Johnson, M. K. Trivett, D. Etemadmoghadam, B. Locandro, N. Traficante, S. Fereday, J. A. Hung, Y. Chiew, I. Haviv,
32

Australian Ovarian Cancer Study Group, D. Gertig, A. deFazio, and D. D.L. Bowtell. Novel molecular subtypes of serous and endometrioid ovarian cancer linked to clinical outcome. Clinical Cancer
Research, 14(16):5198–5208, 2008.
[42] I. Tsamardinos, L. E. Brown, and C. F. Aliferis. The max-min hill-climbing Bayesian network structure
learning algorithm. Machine learning, 65(1):31–78, 2006.
[43] S. van de Geer and P. Bühlmann. `0 -penalized maximum likelihood for sparse directed acyclic graphs.
The Annals of Statistics, 41(2):536–567, 04 2013.
[44] M. Yuan and Y. Lin. Model selection and estimation in regression with grouped variables. Journal of
the Royal Statistical Society: Series B (Statistical Methodology), 68(1):49–67, 2006.
[45] A. Y. Zaitsev. On the Gaussian approximation of convolutions under multidimensional analogues of
S.N. Bernstein’s inequality conditions. Probability Theory and Related Fields, 74(4):535–566, 1987.

33

