Multi-task Learning of Order-Consistent Causal Graphs
Authors: Xinshi Chen
Haoran Sun
Caleb Ellington
Eric Xing
Le Song
Venue: NeurIPS 2021 Poster (2021)
Keywords: Causal discovery
DAG estimation
multi-task learning
graphical models

Abstract: We consider the problem of discovering $K$ related Gaussian directed acyclic graphs (DAGs), where the involved graph structures share a consistent causal order and sparse unions of supports. Under the multi-task learning setting, we propose a $l_1/l_2$-regularized maximum likelihood estimator (MLE) for learning $K$ linear structural equation models. We theoretically show that the joint estimator, by leveraging data across related tasks, can achieve a better sample complexity for recovering the causal order (or topological order) than separate estimations. Moreover, the joint estimator is able to recover non-identifiable DAGs, by estimating them together with some identifiable DAGs.  Lastly, our analysis also shows the consistency of union support recovery of the structures. To allow practical implementation, we design a continuous optimization problem whose optimizer is the same as the joint estimator and can be approximated efficiently by an iterative algorithm. We validate the theoretical analysis and the effectiveness of the joint estimator in experiments.

Multi-task Learning of Order-Consistent Causal Graphs
Xinshi Chen"
Haoran Sun Georgia Institute Technology
Georgia Institute of Technology xinshi. cwwn@gatech. edu
haoransunggatech.. edu
Caleb Ellington Carnegie Mellon University cellingt@cs, cmu. edu
Eric Xing Carnegie Mellon University
Le Song BioMap MBZUAI
MBZUAI
eric. xwwg@iiznaii...... ae
1e. wwwg@mbzuai..c...omomoooo ae
Abstract
We consider the problem of discovering K related Gaussian directed acyclic graphs (DAGs), where the involved graph structures share a consistent causal order and sparse unions of supports. Under the multi-task learning setting, propose l //2-- regularized maximum likelihood estimator (MLE) for learning K linear structural equation models. We theoretically show that the joint estimator, by leveraging data across related tasks, can achieve a better sample complexity for recovering the causal order (or topological order) than separate estimations. Moreover, the joint estimator is able to recover non-identifiable DAGs, by estimating them together with some identifiable DAGs. Lastly, our analysis also shows the consistency of union support recovery of the structures. To allow practical implementation, we design continuous optimization problem whose optimizer is the same as the joint estimator and can be approximated efficiently by an iterative algorithm. We validate the theoretical analysis and the effectiveness of the joint estimator experiments.
Introduction
Estimating causal effects among a set of random variables is of fundamental importance in many disciplines such as genomics, epidemiology, health care and finance [1, 2, 3, 4, 5, 6]. Therefore, designing and understanding methods for causal discovery of great interests in machine learning.
Causal discovery from finite observable data is often formulated as a directed acyclic graph (DAG) estimation problem in graphical models. A major class of DAG estimation methods score-based, which search over the space of all DAGs for the best scoring one. However, DAG estimation remains a very challenging problem from both the computational and statistical aspects [7]. On the one hand, the number of possible DAG structures grows super-xxponentially in the number of random variables, whereas the number of observational sample size is normally small. On the other hand, some DAGs are non-ideenipee -on-identii from observational data even with infinitely many samples.
Fortunately, very often multiple related DAG structures need to be estimated from data, which allows us to leverage their similarity to improve the estimator. For instance, in bioinformatics, gene expression levels are often measured over patients with different subtypes [8, under various experimental conditions [10]. In neuroinformatics, fMRI signals are often recorded for multiple subjects for studying the brain connectivity network [11, 12]. In these scenarios, multiple datasets will be collected, and their associated DAGs are likely to share similar characteristics. Intuitively, it may be beneficial estimate these DAGs jointly.
Work done partially during the visit at MBZUAI (Mohamed hin Zayed University of Artificial Intelligence)
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

In this paper, we focus on the analysis of the multi-task DAG estimation problem where the DAG structures can be related by a consistent causal order and a (partially) shared sparsity pattern, but allowed to have different connection strengths and edges, and differently distributed variables. In this setting, we propose a joint estimator for recovering multiple DAGS based on a group norm regularization. We prove that the joint 4 /llla-pealizzd maximum likelihood estimator (MLE) can recover the causal order better than individual estimators.
Intuitively, it is not surprising that joint estimation is beneficial. However, our results provide quantitative characterization on the improvement in sample complexity and the conditions under which such improvement can hold. We show that:
For identifiable DAGs, if the shared sparsity pattern (union support) size is of order (1) in K where K is the number of tasks (DAGss), then the effective sample size for order recovery will be K where n is the sample size in each problem. Furthermore, as long as 8 is of order ((VK) in K, regularization leads to improvement in sample complexity. A non-identifiable DAG cannot be distinguished by single-task estimators even with indefinitely many observational data. However, non-identifiable DAGS can be recovered by our joint estimator if they are estimated together with other identifiable DAGS.
the joint estimator with group
Apart from the theoretical guarantee, design efficient algorithm for approximating the joint estimator through a formulation of the combinatorial search problem to a continuous programming This continuous formulation contains novel design of a learnable masking matrix, which plays an important role in ensuring the acyclicity and shared order for estimations in all tasks. An interesting aspect of our design is that can learn the masking matrix by differentiable search over continuous space, but the optimum must be contained in a discrete space of cardinality p! (reads p factorial, where p the number of random variables).
We conduct set of synthetic experiments to demonstrates the effectiveness of the algorithm and validates the theoretical results. Furthermore, we apply our algorithm to more realistic single-cell expression RNA sequencing data generated by SERGIO [13] based on real gene regulatory networks.
The remainder of the paper is organized as follows. In Section 2, we introduce the linear structural equation model (SEM) interpretation of Gaussian DAGs and its properties. Section is devoted to the statement of our main results, with some discussion on their consequences and implications. In Sec- tion 4, we present the efficient algorithm for approximating the joint estimator. Section 5 summarizes related theoretical and practical works. Experimental validations are provided in Section 6.
2 Background
substantial body of work has focused on the linear SEM interpretation of Gaussian DAGS [14, 15, 16, 17]. Let X = (X1, ;Xp) be d-mmennsional random variable. Then a linear SEM reads
x= G. X W, W ~ (0,
where Po is a x positive diagonal matrix which indicates the variances of the noise W. Go E RRxx is the connection strength matrix adjacency matrix. Each nonzero entry Goúj represents the direct causal effect of X; on Xj. This model implies that X is Gaussian, X ~ ~N(0, 2), where
: = (1= - -1
Causal order To. The nonzero entries of Go defines its causal order (also called topological order), which informs possible "parents" of each variable. A causal order can be represented by a permutation To over [p] = 21. .2.....)) We say Go is consistent with To if and only if
Goss + 0 + (i) To()).
There could exist more than one permutations that are consistent with DAG structure Go. so we denote the set of permutations that satisfy Eq. by llg. Once a causal order mg is identified, the con- nection strengths Go can be estimated by ordinary least squares regression which is a comparatively easier problem.
Identifiability and equivalent class. However, estimating the true causal order no very challeng- ing, largely due to the existence of the equivalent class described below.

---

Let Sp be the set of all permutations over [p]. For every E Sp, there exists a connection strength
matrix G(T),, which is consistent with æ, and diagonal matrix ((=) diag (T)( ,y (T)
such that the variance Σ in Eq. 2 equals to
= (I ((+))
(4)
and therefore the random variable X in Eq. can be equivalently described by the model [15]:
X= G T x + ( W(x) ~ ((T)
(5)
Without further assumption, the true underlying DAG, Go = G(ro). cannot be identified from the cquivalent class:
{G(T) E Sp}
(6)
based on the distribution of X, even if infinitely many samples are observed.
3 Joint estimation of multiple DAGs
In the multi-task setting, we consider K linear SEMs:
x Go (x) W((), for k == 1,... 1...., K.
The superscript notation (k) indicates the k-th model. As mentioned See 2, each model defines a random variable X ~ NO0 ((0) with g = (1 - ( -- (1- -()))=
3.1 Assumption
Similarity. In Condition 3.1 and Definition 3.1, we make few assumptions about the similarity of these models. Condition 3.1 ensures the causal orders the K DAGS do not have conflicts. Definition 3.1 measures the similarity of the sparsity patterns in the K DAGs. If the sparsity patterns are strongly overlapped, then the size of the support union s will be small. We do not enforce a strict constraint on the support union, but we will see in the later theorem that smaller can lead to better recovery performance.
Condition 3.1 (Consistent causal orders). There exists nonempty set of permutations TT S Sp such that VTO E II. it holds for all & [K] that Go.g) (k) * 0 To(i) < (i)
Definition 3.1 (Support union). Recall G(T) Eq. 6. The support union of the K DAGs associated with permutation π is denoted by S(x) ;= { (i 3Ą E [K] s.t. ((()) (T) A 0}. The support union of the K true DAGs (i(k) is So = (mr) We further denote == S and == sup rEs, |S(T))-
Identifiability. To ensure the consistency of the estimator based on least squared loss, we first assume that the DAGS to recovered are minimum-tran DAGs. Condition 3.2 (Minimum-tracce Recall the equivalent class defined in Eq. to Eq. Assume for all 1= ... K, ( = min trace((2(TT)
However, the minimum-trace DAG may not be unique without further assumptions, making the true DAG indistinguishable from other minimum-trace DAGs. Therefore, we consider the equal variance condition in Condition 3.3 which ensures the uniqueness of the minimum trace DAG. In this paper, assume the first K' < K models satisfy this condition, they identifiable. We do not make such an assumption on the other K - K' models, so the K K'' models may not be identifiable.
Condition 3.3 (Equal variance). For all ls 1,... K' with K' s K. the noise (Kk) (O ))) has equal variance with 22 =00 Ip.
3.2 4/llzpenalized joint estimator
Denote the sample matrix by (*) whose row vectors are i.i.d. samples from (0 I())). Based on the task similarity assumptions, we propose the following se/lz-penalized joint maximum likelihood estimator (MLE) for jointly estimating the connection strength matrices (G(A) )
A, (G(0)
arg min
x
x (()) (x) ()G(())(()))l) (7)
EE,, (D(M)
2n
3

---

Similar to the notation in Eq. 4. (T) indicates its consistency with ((T) denote the space of all DAGs that are consistent with æ. is notable that single π shared across all K tasks is optimized in Eq. 7, which respects Condition 3.1. The group norm over the set of K matrices is defined as
(G))))()))l,, :- E プラ |C] [|,,
where Gij G. [c]]] (T)....
((K) (T)).
will penalize the size of union support in a soft way. When K = 1, this joint estimator will be reduced to the 1 -peenalieed maximum likelihood estimation. Remark 3.1. The optimization in Eq. 7 is used for analysis only. A continuous program with the same optimizer will be discussed in Section 4 for practical implementation.
3.3 Main result: causal order recovery
We start with a few pieces of notations and definitions. Then the theorem statement follows. (k) Definition 3.2. Let 9j (T) denote the j-th column of (T)) Let
d = =u,, UkE[k] supp ( (T))), 9max == SuPrES,,
|| (((-K) (T) (T)
In Definition 3.2, d is the maximal number of parents (in union) in the DAGs (T)), which is also measure of the sparsity level. Imax is bounded by the maximal entry value in the matrices ( (T)). Condition 3.4 (Bounded spectrum). Assume for all k: = 1,.... K, the covariance matrix (k) is positive definite. There exists constants 0 < Amin S Amax co such that for all 1,... K,
(a) all eigenvalues of (k) are upper bounded by Amax (b) all eigenvalues of (A) are lower bounded by Amin-
Condition 3.5 (Omega-min). There exists constant The > 0 such that for any permutations и & Il. pk Σk_
(k) (x)2 (o
0 (x)2)2
Condition 3.5 with K' K 1 is called 'omega-min' condition in [14], so follows this terminology. In some sense, when Yhe is larger, (T) with π & IT is allowed to deviate less from the true variance (To) an with xo Πo, which will make more difficult separate the set To from its complement in a finite sample scenario. Ideally, we should allow The to be large, so that the recovery not only restricted to easy problems.
Now we are ready to present the recovery guarantee for causal order. Theorem 3.1 is a specific statement when the regularization parameter À follows the classic choice in (group) Lasso problems. more general statement which allows other is given in Appendix B along with the proof.
Theorem 3.1 (Causal order recovery). Suppose solve the joint optimization in Eq. 7 with specified regularization parameter V plogg for a set of K problems that satisfy Condition 3.1, 3.4 (a), 3.2, 3.3 and 3.5. the following conditions are satisfidd:
0(n, K, s) = P n
K
> K1 Ki
(8)
8 plogp K
c1 log K + cz(d + loggr K < koplogp.
(9)
(10)
then the following statements hold true:
(a) With probability at least - c3 exp Kg, + 1) log exp (-c,n), it holds that ll.-
(b) If in addition, n satisfies n. 2 Rzd((og K + log p) with à == max 10, 9 lla- and Condition 3.4 (b) holds, then with probability at least es exp (-Rald + 1)) 1)logp) exp (-cin) - exp( log log K).
c 1| = =0 (
1 Ek (1(() R
\
plogg nk

---

In this statement, C1, O2. C3, C1 are universal constants (i.e. independent of n. K, S, ( Ni, N2g, Rig, R4 are constants depending on Imax and Amax. and K5 a constant depending on Amin-
Discussion Eq. 8-Eq. 10. (i) In Eq. 8, Theorem 3.1 identifies sample complexity parameter ((n, K, K', ,,,)) s) Following the terminology in [18], our use of the term "sample complexity" for 0 reflects the dominant role plays in our analysis as the rate at which the sample size much grow in order to obtain a consistent causal order. More precisely, for scalings (n, K, K", P, s) such that (n, K, K', P, s) exceeds fixed critical threshold K1 The we show that the causal order can be correctly recovered with high probability.
(ii) The additional condition for the sample size n in Eq. 9 is the sample requirement for ordinary linear regression problem. It is in general much weaker than Eq. 8, unless K grows very large.
(iii) The last condition in Eq. 10 on K could be relaxed if tighter analysis on the distribution properties of Chi-squared distribution is available. However, is notable that this restriction on the size of K has already been weaker than many related works multi-task l sparse recovery, which either implicitly treat K as a constant [18, 19] or assume K o(log p) [20, 21].
Recovering non-identifiable DAGs. A direct consequence of Theorem 3.1 is that as long as the number of identifiable DAGs K' is non-zero, the joint estimator can recover the causal order of non-identifiable DAGs with high probability. This is not achievable in separate estimation even with infinitely many samples. Therefore, we show that how the information of identifiable DAGs helps recover the non-identifiable ones.
Effective sample size. As indicated by 0(n, K, K", K', 8) in Eq. 8, the effective sample size for recovering the correct causal order is nk'² if the support union size s is of order (11) in K. To show K the improvement in sample complexity, it is more fair to consider the scenario when the DAGs are identifiable, i.e., K' = K. In this case, it is clear that the parameter ((n, K, p) indicates a lower sample complexity relative to separate estimation as long as s is of order o(vK).
Separate estimation. Consider the special case of single task estimation with K= in which the joint estimator reduces to é ei-penalized MLE. We discuss how our result can recover previously known results for single DAG recovery. Unfortunately, existing analyses were conducted under different frameworks with different conditions. [14] and [22] are the most comparable ones since they were based the same omega-min condition (i.e., Condition 3.5), but they chose smaller regularization parameter A. In our proof, the sample complexity parameter is derived from (n, K, K', pp, 8) = pk'/ (xXvK) and Eq. 8 pk'/
K' = 1,
> K1 When K K" 1, this condition matches what is identified in [14] and [22] for recovering the order of single DAG.
(xVKK)
Error of ((R), To compare the estimation of (()) to the true DAG (. Theorem 3.1 (b) says the averaged error in F-norm goes to zero when nk + co. It decreases in K as long as s = o(K).
To summarize, Theorem 3.1 analyzes the causal order consistency for the joint estimator in Eq. 7. Order recovery is the most challenging component in DAG estimation. After â E IT has been identified, the DAG estimation becomes p linear regression problems that can be solved separately. Theorem 3.1 (b) only shows the estimation error of connection matrices in F-norm. To characterize the structure error, additional conditions are required.
3.4 Support union recovery
Theorem 3.1 has shown that á To E Th holds with high probability. Consequently, the support recovery analysis in this section is conditioned on this event. In fact, given the true order To, what remains is set of p separate / /2-penalized group Lasso problems, each of which the order mo plays a role of constraining the support set by the set ()()) := (i To(i) < j.. However, we need to solve p such problems simultaneously where p is large. A careful analysis required, and directly combining existing results will not give a high recovery probability.
In the following, we impose set of conditions and definitions, which are standard in many 1 sparse recovery analyses [18, 19], after which theorem statement follows. See Appendix C for the proof.
Definition 3.3. The union support of j-th columns of {Go (k)) ]AE[K] is denoted by RS, == {i E [p]: Bk E [K] s.t. G # 0}. The maximal cardinality is Tmax = supje(p] |RSy]-
5

---

Definition 3.4. Pa := SUpke[]]], es the conditional covariance matrices, where >0/ s
(>(k)
s is the maximal diagonal entry of
[ 5 (>(S) gge
Definition 3.5. Imin = (s [i /2/ /K represents the signal strength.
Condition 3.6 (Irrepresentable condition). There exists a fixed parameter E (0, such that
SuP, RS [[S S 1 7, where A(S).g == SuP)E[K]
(258))
S88)
(R)
Theorem 3.2 (Union support recovery). Assume on the subset of probability space where { E Hb} holds, and assume Condition 3.6. Assume the following conditions are satisfied
n 2 K67 logp,
(11)
K S e7 log P.
(12)
8A max
logp
2
p plogp Tmax
o(gmin),
(13)
Am T3
Amin
K
where fig is a constant depending on y,, Amin: Pu: muxx and o7 is some universal constant. Then w.p. at
least 1 maxx exp (-c,K log p) exp (-cg logp) - c.K exp (-cg Tmax - llgg))) the support (11:) and that ||@(UKK) k) ((I:K) l1 /12 / /VV o(gmin).
union of (1:K) is the same as that of GD
Discussion on Eq. 11 and Eq. 12. (i) Eq. 11 poses a sample size condition. The value Tmax is the sparsity overlap defined in Definition 3.3 (i). takes value in the interval [d, min P. dK) H. depending on the similarity in sparsity pattern.
(ii) The restriction on K in Eq. 12 plays a similar role as Eq. 10 in Theorem 3.1. This a stronger restriction, but also guarantees the stronger result of support recovery. Existing analyses on ///22 penalized group Lasso were not able to relax this constraint, neither, so some of them treated K as constant in the analysis [18, 19]. Recall that in Theorem 1. we were able allow K (p log p). Technically, this was achieved because in the proof of Theorem 3.1, we avoid analyzing the general recovery of group Lasso, but only its null-consistency (i.e., the special case of true structures having zero support), where tighter bound can be derived and it is sufficient for order recovery.
Benefit of joint estimation. Eq. 13 plays a similar role as Eq. 8 in Theorem 3.1. It specifies rate which the sample size must grow for successful union support recovery. As long as Tmax is of order o(vK), K will effectively reduce the second term in Eq. 8. Apart from that, the recover probability specified in Theorem 3.2 grows in K.
4 Algorithm
Solving the optimization in Eq. 7 by searching over all permutations т E Sp is intractable due to the large combinatorial space. Inspired by the smooth characterization of acyclic graph [16], propose a continuous optimization problem, whose optimizer is the same as the estimator in Eq. 7. Furthermore, we will design an efficient iterative algorithm approximate the solution.
4.1
Continuous program
We convert Eq. the following constrained continuous program
K
min
TERPX GO) ((0) PRPXP
x (k)
X (A)
(1))
1| + A||
(1.K) 11, + p1 - -T (14)
2n
k=1
subject to I(T)) == trace(e p=
(15)
where (k)) == ((k) of is element-wise multiplication between (k) and T. and px is axp matrix with entries equal to one. Eq. 15 is smooth 'DAG-ness' constraint proposed by NOTEARS [16], which ensures is acyelic. One can also use h(T) == trace((I T T) - proposed in [23].
We would like to highlight the novel and interesting design of the matrix T in Eq. 14. What makes Eq. 7 difficult to solve is the requirement that {G(K) must be DAGs and share the same order. A straightforward idea is to apply the smooth acyclic constraint to every ((+), but it is not clear how to enforce their consistent topological order. Our formulation realizes this by a single matrix T.

---

Algorithm 1: Joint Estimation Algorithm
Hyperparameters a, 1,
Initialize G(I:)) randomly;
for itr 1.... .M do
for itr' = 1, .. .M' do
[G T] - GradOptStep(/: ((: 3); > Gradient-baser update on
Vi,j E [p], G3j G(I:k) + |c G
max (0, [|] (| - b Proximal step
Vi,j E [p], Tig + sign(T.y) max {0, TTyl [xllo! l12]:
Proximal step
t- ++ TT();;
Dual ascent
a a (1+ +6;;
b. Typical rule [23]
To better understand the design rationale of T. recall in Eq. that a matrix G DAG of order if and only if its support set is in {(i.5) x(1) < x(j)}. The matrix T plays a role of restricting the support set of G(k) by masking its entries. Two examples are shown above. However, unlike the learning of masks in other papers which allows T to have any combi- nations of nonzero entries, here we need to exactly represent the support set for each æ. That is
21,2,3,4
= (4, 1,3)
$
$
0
0
0
0
1
0 0
0
0
0
0
1
0
0
T
=
0
0
1 1
0
1
0
T is from space Tp with p! elements: Tp= == (T E {0. 1]pxp Tij = Ti) < r(j)). Now
key question arises: How to perform a continuous and disferentiable search over Tp? ? The following finding motivates design:
T Tp = T E arg minreRxx> p11 TI? subject h(T) =]-
In other words, T' is continuous projection of 1pxp u to the space of DAGs. We then optimize the mask T in the continuous space R2*P but the optimal solution must be an element in the discrete space Tz- This observation also naturally leads the design of Eq. 14.
Finally, we want to emphasize that it is important for the optimal T have binary entries. Without this property, any nonzero value can scale the (G. T) pair to give equivalent masked DAG, i.e., Go (cG) (ET). This scaling equivalence will make the optimization hard to solve in practice.
Proofs for the above arguments and the equivalence between Eq. 7 and Eq. 14 are in Appendix F.
4.2 Iterative algorithm
We derive an efficient iterative algorithm using the Lagrangian method with quadratic penalty, which converts Eq. 14 an unconstrained problem:
max ((G():K) 3): = f(G(1:K), ,T; T) +)
(1:K) /lG//:))//
min
T,G C
ER
$20
where K( T; 1
G x
1
x
(2)
T1? f:2 + (T) + ah(T)?
2n
is dual variable, a is the coefficient for quadratic penalty, and f the smooth term in the objective.
We solve this min-max problem by alternating primal updates on (G(I:), T') and dual updates on B. Due to the non mn-smoothness of group norm, the primal update is based on proximal- gradicnnt method, where the proximal-operaton with respect to || 111//12 has closed form:
arg min g 1 E (k() - x (k) [| +
+ cllz(1:K))
11,/1a)3
(1:K)
0 la max(0, || (1:K) |1 c}.
(1:K)
which is group-wise soft-threshold. Since (k)) and T are multiplied together element-wisely inside the group norm, the proximal operator will be applied to both of them separately. Together with the
7

---

dual update for 3, 3 <- + Th(T) with 7 as the step size, we summarize the overall algorithm for solving Eq. 14 in Algorithm
Related work
Single DAG estimation. Unlike the large literature of research on undirected graphical models [24, 25, 26. 27], statistical guarantees for score-based DAG estimator have been available only in recent years. [14, 15, 17] have shown the DAG estimation consistency in high-dimensions, but they do not consider joint estimation. Nevertheless, some techniques in [14, 15] are useful for our derivations.
Multi-task learning. (i) Undirected graph estimation. There have been extensive studies the joint estimation of multiple undirected Gaussian graphical models [28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39]. (ii) DAG estimation. In contrast, not much theoretical work has been done for joint DAG learning. A few pieces of recent works addressed certain statistical aspects of multi-task DAG estimation [9, 20, 21], but [20, 21J tackle the fewer task regime with K o(log p). and [9] assumes the causal order is given. Another related work that we notice after the paper submission is [40], which also assumes the DAGs have a consistent causal order, but focuses on estimating the difference between two DAGs. (iii) Linear regression. Multi-task linear regression is also a related topic [41, 42, 43, 44, 45, 46, 47], because the techniques for analyzing group Lasso are used in our analysis [18, 19].
Practical Algorithms. Works on practical algorithm design for efficiently solving the score- based optimization are actively conducted [48, 49, 50, 51]. Our algorithm is most related to recent methods exploiting a smooth characterization of acyclicity, including NOTEARS [16] and several subsequent works [23, 52, 53, 54], but they only apply for single-task DAG estimation. Although algorithms for the multi-task counterpart were proposed a decade ago [11, 55, 56, 57], none of them leverage recent advances in characterizing DAGs and providing theoretically guarantees.
Experiments
6.1
Synthetic data
The set of experiments is designed to reveal the effective sample size predicted by Theorem 3.1, and demonstrate the effectiveness of the proposed algorithm. In the simulations, we randomly sample causal order ı and a union support set Sg. Then we randomly generate multiple DAGs that follow the order π and have edges contained in the set So. For each DAG, we construct linear SEM with standard Gaussian noise, and sample data points from it. On tasks with different combinations of (p, n. 8, K), we exam the behavior of the joint estimator, estimated by Algorithm on 64 tasks and report the statistics in the following for evaluation. In this experiment, take K' = K so that all the DAGs are identifiable. This simpler case will make it casier to verify the proposed algorithm and the rates in the theorem.
Success probability for order recovery. For each fixed tuple (p. n, K), we measure the sample complexity in terms of the parameter 6 specified by Theorem 1. Fig plots the success probability Pr(f E Πo), versus 0 = p/s h//s/nK//ppoo (plogy) at logarithmic scale. Theorem 3.1 predicts that the success probability should transition to once 0 exceeds a critical threshold. Curves in Fig 1 actually have sharp transitions, showing step-function behavior. The sharpness is moderated by the logarithmic scale in z-axis. Moreover, by scaling the sample size n using 8, the curves align well as predicted by the theory and have similar transition point, even though they correspond to very different model dimensions
thm:supportce-4mm
Recovery Probability
0.6
0.4
p=32
6=
0.2
-28
D 256
2
2
2
Figure 1: Success probability vs log 0_ Fig 2 shows the success probability in the form of Heat Maps, where the rows indicate an increase in per task sample size n, and the columns indicate an increase in the number of tasks K. The results show that the increases in these two quantities have similar effect to the success probability.
8

---

84
128
296
Number of Samples per Task
0:14
0:67
28
0.99
40
0.88
60
1067
0:23
0.91
0.7
0.81 0.85
0.37
80
0.31
0.52
0.88
0.97
0.078
6.42
0:67
0.84
0.016
0:34
0.54
0.82
320 163
0:58
0.8
0.94
160
0.5
0.58
0.88
0.97
0.98
0.17
0.3
82
0.94
0.076
0.39
07
0.85
$20
0.83
0:38
5.0160
0.34
16
32
8
16
52
2
16
32
16
52
Number of Tasks
Number of Tasks
Number of Tasks
Number Tasks
Figure 2: Heat map: Darker colors indicate lower success probability, and lighter colors higher.
0.5
P-
p- 123
2.6
2.4
0.6
k-4
k-4
cm3
2.4
k-8
km
0.4
32
03
k-332
k-16
0.3
16
k2:
0.3
K:32
K-16
6.
011
0.0
1.0
P5 1998)
13.9
6-1
K-1
K-1
0.6
6-2
k+2
1-2
kod
K-4
0.4
3.4
1-3
1-16
1-18
0.2
02
0.0
2
Inummber
[nummerr sampled)
(nummber
".juubeer samplee)
Figure 3: False discovery rate (FDR) and true positive rate (TPR) of the edges.
Effectiveness K structure recovery. In this experiment, we aim at verifying the effectiveness of the joint estimation algorithm for recovering the structures. For brevity. only report the numbers for false discovery rate (FDR) and true positive rate (TPR) of edges in Fig 3, but figures for additional metrics be found in Appendix G. In Fig 3, we observe consistent improvements when increasing the number tasks K. When per task sample size is small, this improvement reveals to be more obvious.
Comparison with other joint estimator. In this experiment, we compare our method MultiDAG with JointGES [21] models with p 32 and 64. Results in Table show that two algorithms have similar performance when = However, when increases, our method returns consistently better structures in terms of SHD.
Table 1: Comparison of MultiDAG (MD) and JointGES (JG) in SHD
p = 32
p == 64
n 10
20
80
320
10
20
80
πt 320
MD(k=1)
39 +
25 5
10 +
+ 3
104 7
77 + 8
29 +
19 +
MD(k=2)
37 +
2
8 8
4 3
103 +
68 +
22 22
13 + 6
MD(k=8)
5 5
1
4
2+1
85 +
42 4
13
6 3
MD(k=32)
23
1113
312
111
66 5
35
1013
3
JG(k=1)
31
19
8
6
100
53
18 +
JG(k=2)
19
945
7
99 10
51 10
20 -
21 10
IG(k=8)
30 5
19 5
12
10
82 10
10
20 -
27 -
IG(k=32)
26 + 4
1844
12
9
57
36 06400
19 1945
26 +
6.2 Recovery of gene regulatory network
We investigate how our joint estimator works on more realistic models, by conducting a set of experiments realistic gene expression data generated by SERGIO [13], which models the additive

---

effect of cooperative transcription factor binding across multiple gene regulators in parallel with protein degradation and noisy expression.
(a) True: Ecoli 100
(b) K 1,
1000
(c) K 10,n= 100
(d) K 1, 100
FDR: 0.11, TPR: 0.47
FDR:0.06, TPR:0.41
FDR: 0.18, TPR: 0.29
Figure 4: Visualization of the recovered DAG structures. Each light green colored pixel at position (i,j) indicates an edge from node i to j.
We conduct this experiment on the E. coli 100 network, which includes 100 known genes and eeieeeeeeeee
137 known regulatory interactions. To evaluate our algorithm, we generate multiple networks by rearranging and re-weighting edges at random in this network without violating the topological order. We simulate the gene expression from each network using SERGIO with a universal non-cooperative hill coefficient of 0.05, which works well with our linear recovery algorithm.
Fig 4 provides a visual comparison of the recovered structures. It can be seem from the true network that there a few key transcription factors that highlight several rows in the figure. These transcription factors are better identified by the two structures in (b) and (c), but not that clear in (d). Combining this observation with the more quantitative results in Table 2, we that the combination (K = 10, n = 100) achieves comparable performance to K = with the same total number of samples, and outperforms the single task estimation with n 100.
K
n
FDR
TPR 0.30
FPR
SHD
1
100
0.18 - 3.80-3
= 1.2c-3 3
0.001 + 7.3c- 7
104.61 + 3.2e1
10
100
0.07 + 1.2e
0.42 - 5e
0.001 + 2.8e-7
84.0 + 1.5e1
1000
0.09 + 2.9e-3
0.50
+ 1.4e-
0.002 + 9.2e-7
76.4 + 5.9e1
Table 2: Recovery across 25 independent initializations SERGIO for each experiment. FPR and SHD stand for false positive rate and structural hamming distance, respectively.
7 Conclusion and discussion
In this paper, we have analyzed the behavior of l 11/22-penaiized joint MLE for multiple DAG estimation tasks. Our main result is to show that its performance in recovering the causal order is governed by the sample complexity parameter (n,, K, K" p) s) in Eq. 8. Besides, have proposed an efficient algorithm for approximating the joint estimator via formulating a novel continuous programming, and demonstrated its effectiveness experimentally. The current work applies to DAGs that have certain similarity in sparsity pattern. It will be interesting to consider whether the joint estimation without the group-norm (and without the union support assumption) can also lead to similar improvement in causal order recovery.
Acknowledgement
Xinshi Chen is supported by the Google PhD Fellowship. This works is done partially during visit at MBZUAL We are grateful for the computing resources provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology, Atlanta. We thankful for PACE Research Scientist Fang (Cherry) Liu's excellent HPC consulting.
References
[1] Nir Friedman, Michal Linial, Iftach Nachman, and Dana Pe'er. Using bayesian networks to analyze expression data. Journal of computational biology, 7(3-))601-620,, 2000.
10

---

[2] James M Robins, Miguel Angel Hernan, and Babette Brumback. Marginal structural models and causal inference in epidemiology, 2000.
[3] Pedro Aguilera Aguilera, Antonio Feráández, Rosa Fernández. Rafael Rumí, and Antonio Salmerón. Bayesian networks in environmental modelling. Environmental Modelling & Software, 26(12) 2011.
[4] A Nicholson, F Cozman, M Velikova, JT Van Scheltinga, PJ Lucas, and M Spaanderman. Applications of bayesian networks exploiting causal functional relationships in bayesian net- work modelling for personalised healthcare. International Journal of Approximate Reasoning. 55(1):59-73, 2014.
[5] Bin Zhang, Chris Gaiteri, Liviu-Gabrie Bodea, Zhi Wang, Joshua McElwee, Alexci A Podtelezh- nikov, Chunsheng Zhang, Tao Xie, Linh Tran, Radu Dobrin, et al. Integrated systems approach identifies genetic nodes and networks in late-onset alzheimer's disease. Cell, 153(3):707-720, 2013.
[6] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: founda- tions and learning algorithms. The MIT Press, 2017.
[7] David Maxwell Chickering. Learning bayesian networks is np-complete. In Learning from data, pages 121-130. Springer, 1996.
[8] T Tony Cai, Hongzhe Li, Weidong Liu, and Jichun Xie. Joint estimation of multiple high- dimensional precision matrices. Statistica Sinica, 26(2):445, 2016.
[9] Jianyu Liu, Wei Sun, and Yufeng Liu. Joint skeleton estimation of multiple directed acyclic graphs for heterogeneous population. Biometrics, 75(1):36-44 2019.
[10] Teppei Shimamura, Seiya Imoto, Rui Yamaguchi, Masao Nagasaki, and Satoru Miyano. Infer- ring dynamic gene networks under varying conditions for transcriptomic network comparison. Bioinformatics,
2010.
[11] Shohei Shimizu. Joint estimation of linear non-gaussian acyclic models. Neurocomputing, 81:104-107, 2012.
[12] Aiying Zhang, Gemeng Zhang, Vince D Calhoun, and Yu-Ping Wang. Causal brain network in schizophrenia by a two-step bayesian network analysis. In Medical Imaging 2020: Imag- ing Informatics for Healthcare, Research, and Applications, volume 11318, page 1131817. International Society for Optics and Photonics, 2020.
[13] Payam Dibaeinia and Saurabh Sinha. Sergio: single-cell expression simulator guided by gene
regulatory networks. Cell Systems, 1 11(3):252-22 271211,, 2020.
[14] Sara Van de Geer, Peter Bühlmann, al. fo-penalized maximum likelihood for sparse directed acyclic graphs. Annals of Statistics, 41 (1(2)5336-567 2013.
[15] Bryon Aragam, Arash A. Amini, and Qing Zhou Learning directed acyclic graphs with penalized neighbourhood regression, 2017.
[16] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning. Advances in Neural Information Processing Systems, 31:9472-9483, 2018
[17] Bryon Aragam, Arash Amini, and Qing Zhou. Globally optimal score basee learning of directed acyclic graphs in high-dimensions. In Advances in Neural Information Processing Systems, volume 32, 2019.
[18] Guillaume Obozinski, Martin Wainwright, Michael Jordan, et al. Support union recovery in high-dimensional multivariate regression. The Annals of Statistics, 39(1):1-47, 2011.
[19] Weiguang Wang, Yingbin Liang, and Eric P Xing. Collective support recovery for multi-design multi-response linear regression. IEEE Transactions Information Theory, 61(1):513-534 2014.
[20] Yuhao Wang, Santiago Segarra, Caroline Uhler, et al. High-dimensional joint estimation of multiple directed gaussian graphical models. Electronic Journal of Statisties, 14(1):2439-244 2483,, 2020.
[21] Kyoungjae Lee and Xuan Cao. Bayesian joint inference for multiple directed acyelic graphs. arXiv preprint arXiv:2008.061900 2020.

---

[22] Magali Champion, Victor Picheny, and Matthieu Vignes. Inferring large graphs using 1 (1)- penalized likelihood (vol 28, pg 905, 2018). STATISTICS AND COMPUTING, 28(6): 2018.
[23] Yue Yu, Jie Chen, Tian Gao, and Mo Yu. Dag-gnn: Dag structure learning with graph neural networks. International Conference on Machine Learning, pages 7154-7163. PMLR, 2019.
[24] DR Cox and Nanny Wermuth, Multivariate Dependencies: Models, Analysis and Interpretation, volume 67. CRC Press, 1996.
[25] Nicolai Meinshausen, Peter Bühlmann, et al. High-dimensional graphs and variable selection with the lasso. Annals of statistics, 34(3):1436-148 2006.
[26] Ming Yuan and Yi Lin. Model selection and estimation in the gaussian graphical model. Biometrika, 94(1):19 35, 2007.
[27] Pradeep
Ravikumar, Garvesh Raskutti, Martin J Wainwright, and Bin Yu. Model selection in gaussian graphical models: High-dimensional consistency of 11-regularried mle. In NIPS, pages 1329-1336, 1336, 2008.
[28] Le Song, Mladen Kolar, and Eric P Xing. Keller: estimating time-varying interactions between genes. Bioinformaties, 25(12)i128-i136 2009.
[29] Jean Honorio and Dimitris Samaras. Multi-task learning of gaussian graphical models. In ICML, 2010.
[30] Jian Guo, Elizaveta Levina, George Michailidis, and Ji Zhu. Joint estimation of multiple graphical models. Biometrika, 98(1):1-15, 2011.
[31] Tony Cai, Weidong Liu, and Xi Luo. A constrained & minimization approach to sparse precision matrix estimation. Journal of the American Statistical Association, 106(494):594 607, 2011.
[32] Diane Oyen and Terran Lane. Leveraging domain knowledge in multitask bayesian network structure learning. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 26, 2012.
[33] Patrick Danaher, Pei Wang, and Daniela 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, 2014.
[34] Mladen Kolar, Le Song, Amr Ahmed, Eric P Xing, et al. Estimating time-varying networks. Annals of Applied Statistics, 4(1):94-123, 2010.
[35] Karthik Mohan, Palma London, Maryam Fazel, Daniela Witten, and Su-In Lee. Node-based learning of multiple gaussian graphical models. The Journal of Machine Learning Research, 15(1):445-48) -4804888888.... 2014.
[36] Christine Peterson, Francesco C Stingo, and Marina Vannucci. Bayesian inference of multiple gaussian graphical models. Journal of the American Statistical Association, 1106509):159-144 2015.
[37] Sen Yang, Zhaosong Lu, Xiaotong Shen, Peter Wonka, and Jieping Ye. Fused multiple graphical lasso. SIAM Journal on Optimization, 25(2):91 9-344,, 2015.
[38] André R Gonçalves, Fernando J Von Zuben, and Arindam Banerjee. Multi-task sparse structure learning with gaussian copula models. The Journal of Machine Learning Research, 17(1): 1234, 2016.
[39] Burak Varici, Saurabh Sihag, and Ali Tajer. Learning shared subgraphs in ising model pairs. In International Conference on Artificial Intelligence and Statistics, pages 3952-330 3960. PMLR, 2021.
[40] Asish Ghoshal, Kevin Bello, and Jean Honorio. Direct learning with guarantees of the difference dag between structural equation models. arXiv preprint arXiv:1906.1200 2019.
[41] Francis R Bach. Consistency of the group lasso and multiple kernel learning. Journal of Machine Learning Research, 96)), 2008.
[42] Sahand Negahban and Martin J Wainwright. Joint support recovery under high-dimensional scaling: Benefits and perils of f rrgalarrgalaiiiitii Advances Neural Information Processing Systems, 21:1 1t::161-1166 2008
12

---

[43] Peng Zhao, Guilherme Rocha, Bin Yu, et al. The composite absolute penalties family for grouped and hierarchical variable selection. The Annals of Statistics, 37(6A):346833499 2009.
[44] Karim Lounici, Massimiliano Pontil, Alexandre B Tsybakov, and Sara Van De Geer. Taking advantage of sparsity in multi-task learning. arXiv preprint arXiv:0903. 1468, 2009.
[45] Junzhou Huang, Tong Zhang, et al. The benefit of group sparsity. The Annals of Statistics, 38(4):1978-25 2004, 2010.
[46] Jianhui Chen, Jiayu Zhou, and Jieping Ye. Integrating low-rank and group-sparse structures for robust multi-task learning. Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 42-50, 2011.
[47] Lu Tian, Pan Xu, and Quanquan Gu. Forward backward greedy algorithms for multi-task learning with faster rates. In UAI, 2016.
[48] Mauro Scanagatta, Antonio Salmerón, and Fabio Stella. A survey bayesian network structure learning from data. Progress in Artificial Intelligence, 8(4):425-44 439, 2019.
[49] Ignavier Ng. Zhuangyan Fang. Shengyu Zhu, Zhitang Chen, and Jun Wang. Masked gradient- based causal structure learning. arXiv preprint vrrXiv/990.000
2019.
[50] Shengyu Zhu, Ignavier Ng, and Zhitang Chen. Causal discovery with reinforcement learning.
In International Conference on Learning Representations, 2020.
[51] Hasan Manzour, Simge Küçükkyavzz, Hao-Hsiam Wu, and Ali Shojaie. Integer programming for learning directed acyclic graphs from continuous data. Informs Journal on Optimization, 3(1):46-73, 2021.
[52] Sébastien I Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient- based neural dag learning. International Conference on Learning Representations, 2020.
[53] Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning sparse nonparametric dags. In International Conference on Artificial Intelligence and Statistics, pages 3414-3425. PMLR, 2020.
[54] Roxana Pamfil, Nisara Sriwattanaworachai,, Shaan Desai, Philip Pilgerstorfer, Konstantinos Georgatzis, Paul Beaumont, and Bryon Aragam. Dynotears: Structure learning from time-series data. In International Conference Artificial Intelligence and Statistics, pages 1595-1605. PMLR, 2020.
[55] Robert Tillman and Peter Spirtes. Learning equivalence classes of acyclic models with latent and selection variables from multiple datasets with overlapping variables. In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, pages 3-15. JMLR Workshop and Conference Proceedings, 2011.
[56] Robert E Tillman. Structure learning with independent non-identically distributed data. In Proceedings of the 26th Annual International Conference on Machine Learning, pages 1041- 1048, 2009.
[57] Lele Xu, Tingting Fan, Xia Wu, WeWei Chen, Xiaojuan Guo, Jiacai Zhang, and Li Yao. pooling-lingant algorithm for effective connectivity analysis of fmri data. Frontiers in computational neuroscience, 8:125, 2014.
13


Reviews:
[{"review_id": "cncp70Fm6wi", "reviewer": "Reviewer_x4HL", "summary": "The paper considers structure discovery of multiple but related DAGs in the sense that they share a common causal order. Authors propose an l1/l2-regularized MLE for multiple linear Gaussian models and design a continuous optimization method to obtain the estimates. Sample complexity analysis are also provided regarding causal order and union support recovery, which is interesting and novel. Experiments also validate the theoretical finding and the effectiveness of joint estimation. The paper is generally well-written, though some details seem missing in the experiment part.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "### Pros:\n\n- A quantitative theoretic analysis of sample complexity regarding causal order and union support recovery, which is interesting and novel.\n- An efficient algorithm to estimate the graph.\n- Paper is generally very well-written.\n\n### Cons:\n\nPerhaps the most major issue to me is the setting of the problem. \n-\tAuthors consider linear-Gaussian models where all the noise variances are identical, in particular, are set to be 1 among all datasets. The consideration for this choice is from identifiability. However, it is believed (at least I believe) that the condition of identical noise variances is quite restricted in practice, and it is not testable from observations. Consequently, the proposed theory and algorithm may have little use in practice.\n-\tAuthors are encouraged to also consider non-identifiable case. In this case, it might be beneficial to reduce some indeterminacies of the edges in Markov equivalence class, since you have some different but related DAGs?\n-\tWhile it is claimed the analysis could easily extend to a more general case where the noise variances are different across datasets, it is not clear how they would affect the structure recovery result. In particular, if a dataset has a small very noise variance while another has a large one, then is the result affected?\n-\tWhat if the DAGs do not share the same causal order but can have a very few violations? Can the proposed algorithm return any guaranteed results?\n\nOther issues in the experiment part:\n-\tThe paper is well-written in other parts, but this section seem to omit some import details: 1) running time is not reported. 2) since you have multiple DAGs, it is not clear what DAGs are counted in calculating the metrics (FDR, SHD); are they the average? or simply the result regarding the causal order? 3) standard deviations are not provided.\n-\tImport baselines are missing: for linear Gaussian models with equal noise variances, a well known method is to check the marginal variances of the variables; the largest one would be the source node; for other variables, apply linear regression with the previous variable and obtain residuals; again check the marginal variances of these residuals and follow the same steps; one would then obtain a causal order. How does this procedure compare with the proposed method, when, e.g., in the second experiment that a single dataset can have a number of samples?\n-\tIt would be also beneficial to consider some more settings: different graph structures and densities; different sample numbers across datasets.\n\n\n**update after reading authors' response** \n\nI appreciate authors' response which clearly extends the practical use of this work. I have decided to increase my score to 6 (tentatively), and I am also waiting for other reviewers' opinions on the pratical usage to see if it is possible to further increase the evaluation.\n\nTwo followed questions: \n1. Even without identifiability of the unique graph, how about thinking about the MEC? i.e., how close the estimated graph could be to the MEC?\n2. Please check the followed comment on the mentioned baseline. Sorry for any inconvience.\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I did not find any limitations and negative societal impact of the submission.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5"}, "extra_text": {}}, {"review_id": "PUtNUTbMZW4", "reviewer": "Reviewer_Masb", "summary": "Authors study the problem of learning the causal ordering and union support of DAGs from K linear SEMs. Authors propose a regularized joint estimator and show that this joint estimation enjoys a better sample complexity than single separate estimations. Moreover, guarantees for causal order and union support recovery are provided. For experiments, authors propose a continuous optimization problem that outputs the same causal order and union support than the joint estimator.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Here I present my review along the axes from the NeurIPS reviewing guidelines:\n\n**Originality**\n\n- Are the tasks or methods new? \n  - Yes. Authors study the problem of multitask learning of causal ordering which has been previously unexplored to the best of my knowledge.\n\n- Is the work a novel combination of well-known techniques?\n  - It is motivated by existing methods for multitask learning in the context of Markov random fields.\n\n- Is it clear how this work differs from previous contributions? \n  - Yes. The setting itself is different from existing ones.\n\n- Is related work adequately cited?\n  - Yes. I think the work can benefit from referencing additional works. See below.\n\n\n**Quality** \n\n- Is the submission technically sound? \n  - Yes. Proofs seem sound although I have not checked in full detail.\n\n- Are claims well supported (e.g., by theoretical analysis or experimental results)? \n  - Yes. Proos are provided in the appendix. Empirical validation is provided for synthetic and real-world data.\n\n- Are the methods used appropriate? \n  - Yes. The theoretical analyses use tools from optimization, linear algebra and concentration inequalities.\n\n- Is this a complete piece of work or work in progress? \n  - Complete work.\n\n- Are the authors careful and honest about evaluating both the strengths and weaknesses of their work?\n  - Discussions on the limitations of the presented approach can be beneficial.\n\n\n**Clarity** \n\n- Is the submission clearly written? \n  - Yes.  \n\n- Is it well organized? \n  - Yes. I was able to follow the main contributions without issues.\n\n- Does it adequately inform the reader? \n  - Yes. Experimental details are provided in the appendix. The work is reproducible.\n\n\n**Significance** \n\n- Are the results important? \n  - Yes. I think the topic is of interest to the graphical models community.\n\n- Are others (researchers or practitioners) likely to use the ideas or build on them? \n  - Probably yes. I think the method proposed is an interesting direction to the setting under study.\n\n- Does the submission address a difficult task in a better way than previous work? \n  - The work can be considered as a first set of results for multitask DAG learning.\n\n- Does it advance the state of the art in a demonstrable way? \n  - Yes.\n\n# Questions / Comments\n\n1. I believe the work might benefit from a discussion section on the limitations of the results. For instance, Condition 3.3 (Bounded spectrum) might be an assumption that will not work for the high-dimensional case. In such scenario, due to n << p, there will be zero eigenvalues and the condition will not hold. Can authors comment on this?\n\n2. The statement of Theorem 3.2. has been updated in the appendix, where authors claim that the appendix version is better than the one in the main paper. I am not sure if one is allowed to update this much a relevant part of the contributions, which is a major concern for me.\n\n3. I think the paper can benefit from referencing works that also have a \"direct\" or \"joint\" goal instead of separate estimations. For instance:\n\n    * In undirected networks: \"Learning Shared Subgraphs in Ising Model Pairs\" by Varici et al.\n    * In linear SEMs: \"Direct Learning with Guarantees of the Difference DAG Between Structural Equation Models\" by Ghoshal et al.\n\n  In both works above, there is also evidence of the advantages of a joint objective instead of performing separate estimations.\n\nTypos:\n* L51: p was not introduced.\n* L245: Appendix F not G.\n\n---\nI thank the authors for their response. After reading other reviews and responses, I feel more positive about the paper and will increase my score to 7. I think that the paper could benefit by adding discussions on some points brought by the reviewers.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I think the limitations can be stated more explicitly in a separate section or paragraph.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "ANmpX1oZgGw", "reviewer": "Reviewer_T7T2", "summary": "The paper proposes a multi-task DAG learning under equal variance and Gaussian assumptions. Consistency results are provided and new algorithm is proposed to solve the formulation. Empirically, results are shown to be consistent with the theoretical analysis. ", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strength:\n\nRecovery main results are serious work and appreciated.\n\nFormulation of objective is novel\n\nWeakness:\n\nStrong assumptions on equal variance and Gaussian, although I understand the necessity for the theoretical analysis. \n \nEmpirical evaluation is very weak. \n \n\nComments:\n- L51: it is not clear on what cardinality p is. \n- The model and method only work with Gaussian data, so the current title and introduction is somewhat misleading. Authors should make it clear that the model under consideration is only Gaussian DAG. \n- Equal variance and Gaussian assumptions: they are pretty strong assumption in the context of DAG identification and learning. Such a restricted space may limit its practical values. \n- causal order recovery theoretical analysis is thorough (although I did not check the proof details)\n- Typical Gaussian network structures are learned via precision matrices. It seems possible that these precision matrices can help learning the DAG, e.g., serve as an skeleton of the DAG learning. Have authors considered this aspect?\n- \"the optimal solution must be an element in the discrete space Tp\": would be good to have a formal statement on this.\n- Algorithm 1 uses a min-max formulation and solved slightly different from original NOTEAR. Can authors discuss the advantages of this method?\n- There seem no baseline methods for comparison.  This is a major concern. Why are [27][28] not compared?\n- Gene regulatory network: it seems K = 1, n = 1000 has a better TPR: 0.47 and better SHD.  The results are therefore not very convincing for the proposed approach. What about cases with K > 1 and n = 1k? ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Did not find such a section. ", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3.5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This was a very borderline paper. After the author response and discussion, the one reviewer recommending rejection raised their score to a weak accept. Please make sure to read all of the reviewer feedback carefully and make suggested/promised changes for the camera ready. In particular, reviewer x4HL had some additional followup concerns that should be addressed.
Decision: Accept (Poster)