Truncated Matrix Power Iteration for Differentiable DAG Learning
Authors: Zhen Zhang
Ignavier Ng
Dong Gong
Yuhang Liu
Ehsan M Abbasnejad
Mingming Gong
Kun Zhang
Javen Qinfeng Shi
Venue: NeurIPS 2022 Accept (2022)
Keywords: DAG
Structure Learning
Causal Discovery

Abstract: Recovering underlying Directed Acyclic Graph (DAG) structures from observational data is highly challenging due to the combinatorial nature of the DAG-constrained optimization problem. Recently, DAG learning has been cast as a continuous optimization problem by characterizing the DAG constraint as a smooth equality one, generally based on polynomials over adjacency matrices. Existing methods place very small coefficients on high-order polynomial terms for stabilization, since they argue that large coefficients on the higher-order terms are harmful due to numeric exploding. On the contrary, we discover that large coefficients on higher-order terms are beneficial for DAG learning, when the spectral radiuses of the adjacency matrices are small, and that larger coefficients for higher-order terms can approximate the DAG constraints much better than the small counterparts. Based on this, we propose a novel DAG learning method with efficient truncated matrix power iteration to approximate geometric series based DAG constraints. Empirically, our DAG learning method outperforms the previous state-of-the-arts in various settings, often by a factor of $3$ or more in terms of structural Hamming distance. 

Truncated Matrix Power Iteration for Differentiable DAG Learning
Zhen Zhang" '1 Ignavier Ng 22 Dong Gong,, Yuhang Liu', Ehsan M Abbasnejad!, Mingming
Gong Kun Zhang 2,5 and Javen Qinfeng Shi!
The University of Adelaide Carnegie Mellon University The University of New South Wales *The University of Melbourne
Mohamed bin Zayed University of Artificial Intelligence
yarhunngg2.uuuunggg...g Iiu01, ,ehsan jaatismeeadjjaaen... @ adelaide edu.au
fignavierng @ emu.edu dong. gong @ unsw. eu mingming trrgmmnng-gong gong @ unimell edu.
Abstract
Recovering underlying Directed Acyclic Graph (DAG) structures from observa- tional data is highly challenging due to the combinatorial nature of the DAG- constrained optimization problem. Recently. DAG learning has been cast as a continuous optimization problem by characterizing the DAG constraint as a smooth equality one, generally based on polynomials over adjacency matrices. Existing methods place very small coefficients on high-order polynomial terms for stabiliza- tion, since they argue that large coefficients on the higher-order terms are harmful due to numeric exploding. On the contrary, higher-order terms are beneficial for DAG learning, when the spectral radiuses of the adjacency matrices small, and that larger coefficients for higher-order terms can approximate the DAG constraints much better than the small counterparts. Based on this, we propose novel DAG learning method with efficient truncated ma- trix power iteration to approximate geometric series based DAG constraints. Empir- ically, DAG learning method outperforms the previous state-of-the-aae in vari- ous settings, often by a factor of 3 or more in terms of structural Hamming distance.
discover that large coefficients
Introduction
Recovery of Directed Acyclic Graphs (DAGs) from observational data is classical problem in many fields, including bioinformatics 31 49
machine learning 19], and causal inference [39 The
graphical model produced by DAG learning algorithm allows one to decompose the joint distribution over the variables of interest in a compact and flexible manner, and under certain assumptions these graphical models can have causal interpretations.
DAG learning methods can be roughly categorized into constraint-based and score-based methods. [37
Most constraint-based approaches, e.g., PC
FCI
rely on conditional independence The score-based approaches, including exact
[1]
38
test and thus may require large sample size
33
methods based on dynamic programming [18
A search [47 48 and integer programming D model the validity of a graph according to some score function and are often formulated and solved discrete optimization problem. A key challenge for
35
36
133, as well as greedy methods like GES
score-based methods is the combinatorial search space of DAGs T0.
There is a recent interest in developing continuous optimization methods for DAG learning in the past few years. In particular, Zheng al. [50] develops an algebraic characterization of DAG constraint
"Equal Contribution
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

based on continuous function, specifically matrix exponential, of the weighted adjacency matrix of a graph, and applies to estimate linear DAGs with equal noise variances using continuous constrained optimization 4 S The resulting algorithm is called NOTEARS. Yu et al. [5 proposes
an alternative DAG constraint based on powers of binomial to improve practical convenience. Wei et al. [44] unifies these two DAG constraints by writing them in terms of a general polynomial one. These DAG constraints have been applied to more complex settings, e.g., with nonlinear relationships 20 25 15 [5], time series [26]. confounders [6], and flexible score functions [53] and have been shown to produce competitive performance. Apart from the polynomial based constraints, Lee et al.
[21], Zhu et al. [52] present alternative formulations of the DAG constraints based on the spectral radius of the weighted adjacency matrix. Instead of solving a constrained optimization problem, 4d proposes NOCURL which employs algebraic representation of DAGs such that continuous unconstrained optimization can be carried out in the DAG space directly. However, the focus of these works [211 46 52
Yu et al.
is to improve the scalability, and the performance of the resulting
methods may degrade.
To avoid potential numeric exploding that may be caused by high-order terms in the volynomial-basee DAG constraints, previous works use very small coefficients for high-order terms in DAG constraints More precisely, the coefficients decrease quickly to very small values for the relatively higher order terms that may contain useful higher-order information. This may cause difficulties for the DAG constraints to capture the higher-order information for enforeing acyclicity. With such DAG constraints, a possible solution to maintain higher-order information for acyclicity is to apply an even larger penalty weight 45 50), which, however, causes ill ill-conditioning issues demonstrated by Ng et al. (24). Moreover, these small coefficients will be multiplied on the gradients of higher-order terms during the optimization process, leading to gradient vanishing. In this paper, show that, on the contrary, larger coefficients may be safely used to construct more informative higher-order terms without introducing numeric exploding. The reason is that the weighted adjacency matrix of a DAG must be nilpotent; thus the candidate adjacency matrix often has a very small spectral radius. As result, the higher-order power of the matrix may be very close to zero matrix, and larger coefficients may be safely used without causing numeric exploding. We thus propose to use geometric series-based DAG constraint without small coefficients for more accurate and robust DAG learning. To relieve the computational burden from the geometric series, we propose an efficient Truncated Matrix Power Iteration (TMPI) algorithm with theoretically guaranteed error bound. We summarize our contributions as the following:
45 500
We demonstrate that the key challenge for DAG learning is the gradient vanishing issue instead of gradient and numeric exploding. For sparse graphs, the gradient vanishing issue can be severe.
We design DAG constraint based on finite geometric series which is an order-d polynomial over adjacency matrices to escape from gradient vanishing. We show that the relatively large coefficients on the higher-order terms would not result in gradient exploding in practice.
We show that there exists some constant & & d, s.f. our order-d polynomial constraint can be reduced to order-k polynomial without expanding its feasible set. Though finding such value of k requires solving the NP-hard longest simple path problem. we develop a simple heuristic to find such i with bounded error to the exact, order-d, DAG constraint. Based on this result, we propose a DAG constraint based on efficient Truncated Matrix Power Iteration (TMPI), and conduct experiments demonstrate that the resulting DAG learning method outperforms previous methods by a large margin.
We provide a systematical empirical comparison of the existing DAG constraints, includ- ing the polynomial and spectral radius based constraints, well as the algebraic DAG representation of NOCURL.
2
Preliminaries
DAG Model and Linear SEM
Given a DAG model (a.k.a. Bayesian network) defined over random
vector x [21, .xa]
and DAG G, the distribution P(x) and DAG G are assumed to satisfy We say that x follows a linear Structural Equation Model (SEM) if
the Markov assumption 27 39
x = B T x cc
(1)

---

where B is the weighted adjacency matrix representing the DAG G, and = [e1. e2.... ................................ eal" denotes the exogenous noise vector consisting of independent random variables. With slight abuse of notation, we denote by ((B) the graph induced by weighted adjacency matrix B. Moreover, we do not distinguish between random variables and vertices or nodes, and use these terms interchangeably.
Our goal is to estimate the DAG & from nt i.id. samples of x, indicated by the matrix X E IRaxd In general, the DAG 9 can only be identified up its Markov equivalence class under the faithfulness [39 or sparsest Markov representation [30] assumption. It has been shown that for linear SEMs with homoscedastic errors, from which the noise terms are specified up to a constant
and for linear
non-Gaussian SEMs, from which no more than one of the noise term is Gaussian 34]] the true DAG can be fully identified. In this work, we focus on the linear SEMs with equal noise variances [8] which can be viewed as an instance of the former.
DAG Constraints for Differentiable DAG Learning Recently, Zheng al. 50 reformulated the DAG learning problem as the following continuous optimization problem with convex least-squares objective but non-convex constraint,
min BERUXA
|
X B + +
Ill, subject to hesp [B @ B 0,
(2)
where -II and || || denote the Frobenius norm and elementwise & norm, respectively, 1j is the regularization coefficient that controls the sparsity of B, and @ denotes the Hadamard product, which is used to map to a positive weighted adjacency matrix with the same structure. Given any B e Baxd the function heep is defined as 30
heep (B) - (r(e -
d=
Σ
and can be efficiently computed using various approaches 21. The key here is that hesp
B
and only G(B) is a DAG. The formulation involves hard constraint, and Zheng et 500
adopted the augmented Lagrangian method 4
5 to solve the constrained optimization problem.
Ng al. [24] demonstrated that augmented Lagrangian method here behaves similarly to the
classical quadratic penalty method 15 29
in practice, and showed that it is guaranteed to return
DAG under certain conditions.
To improve practical convenience, Yu et al. 45 proposed similar DAG constraint using an order-d polynomial based on powers binomial, given by
Abin [B) = tr (1+ (I+aB) =d tr Σ a
(4)
where is hyperparameter and is often set to 1/d, e.g., in the implementation of DAG-GNN A5]. This function can be evaluated using Olog d) matrix multiplications using the exponentiation by squaring algorithm. Both hesp (B) and hoin (B) are formulated by Wei et al. as unified one
hpoly (B) = tr Σ cB
4 > 0, i= =11,2.... ,d.
(5)
DAG Learning with Truncated Matrix Power Iteration
In this section, we first show that gradient vanishing is one main issue for previous DAG constraints and propose truncated geometric series based DAG constraints to escape from the issue in Section Then in Section 3.2 the theoretical properties of our DAG constraints are analyzed. In Section we provide an efficient algorithm for our DAG constraints. In Section 3.4 between our algorithm and previous works.
3.3
we discuss the relation
3.1
Gradient Vanishing Issue and Truncated Geometric DAG Constraint
Theoretically, any positive values of e suffice to make the DAG constraint term /pply (B) in Eq. 5 valid one, i.e., it equals zero if and only if the graph ((B) is acyclic. Thus previously small
3

---

GES
MMHC
NOCURL
MSE Single_ter term
MSE Exponential
MSE Binomial
MSE NOBEARS
MSE TMPI(OUrs)
10
10a
Norm oradiinn
10
SHD
10
10"10
Genmetric
104
Exponewiaal
10
Binomial
50
75
100
30
100
Order
#Nodes
Figure 1: Previous polynomial based DAG constraints, including exponential (Eq. B and binomial (Eq. with a = 1/d) constraints, suffer from gradient vanishing. Left: The norm gradient on higher-order terms w.r.t. order on 100-node graph generated using the ER2 setting in Section [ Note that the gradients on different order terms refer V tr(B)), where, for different polynomial based DAG constraints, the coefficients c; different. Right: We propose a geometric series based DAG constraint, named Truncated Matrix Power Iteration (TMPI), to escape from vamished gradients, which leads to a significant performance improvement in terms of Structural Hamming Distance (SHD) on ER2 graphs (see Section 4 for full experimental details and results).
i
coefficients for higher-order terms are preferred since they are known to have better resistance to numeric exploding higher-order polynomials 45 [50]. These small coefficients forbid informative higher-order terms and cause severe gradient vanishing, as nilpotent properties of DAG result in candidate adjacency matrix with very small spectral radius (see empirical analysis in Figure left). In this case, the weight updates during the optimization process can hardly utilize this higher-order information from the vanishingly small gradients. As a consequence, only lower-order information can be effectively captured, which may be detrimental to the quality of final solution.
To avoid the severe vanishing gradient issue in existing DAG constraints (e.g., Eq. propose to discard the coefficients and use finite geometric series based function:
hgoo B tr Σ
B
(6a)
V hgo (B)
Σ (i 1i+1)B
(6b)
The gradient of hgco w.r.t. B denoted by Viligoo (B) is an order-(d 1) polynomial of B.
Compared and 4. Eq. 6a) is more precise DAG constraint escaped from from severe gradient vanishing. Also is notable that due to the nilpotent properties of weighted adjacency matrices of DAGs, the higher-order terms in (6a) and (bb) often comes with small norms; thus, DAG constraint may not suffer from numeric exploding in practice, based on the observations in Section
4
In Figure M we compared the norm of gradients on higher-order terms of different DAG constraints, where we observe that for previous exponential-baan -ased 3 and binomial-based 4 DAG constraints, the gradients of higher-order terms converges to 0 very quickly as the order increases. Meanwhile, the gradient of higher-order terms in (6a) converges to 0 with much slower speed and thus informative higher-order terms are preserved.
We further show that those non-informative higher-order terms that are close to zero in (6a) can be safely dropped form more efficient DAG constraint with some i & d, where we often have I < d in practice:
1 k
Σ
B
Naively, the DAG constraint requires ((k) matrix multiplication which is far less than O(d) for 6a and under mild conditions the error of is tightly bounded (see Proposition B for details). In

---

Figure D right we show that our DAG constraint 7 (named TMPI) attains far better performance than previous ones escapes from gradient vanishing.
3.2 Theoretical Properties of Truncated Geometric DAG Constraint
We analyze the properties of the power of adjacency matrices, which play an important role in constructing DAG constraints. First notice that, for the 1thh power of an adjacency matrix E 2
its entry (B" ); is nonzero if and only if there is a length-k path (not necessarily simple) from node X; node X, in graph (BB). This straightforwardly leads the following property, which is corollary of 44 Lemma 1].
Proposition 1. Ler B E IR be the weighted adjacency matrix of a graph G with d vertices. G is DAG if and only if 0.
Proposition ! can be used to construct a single-term DAG constraint term, given by
Rsingle ( ||
(8)
where -| denotes an elementwise bp norm. However, as the candidate adjacency matrix gets closer acyclic, the spectral radius of B will become very small. As a result, the value h as well as its gradient V singte (B) may become too small to be represented in limited machine precision. This suggests the use of the polynomial based constraints for better performance.
Without any prior information of the graph, one has to compute the dh power of matrix B, either using the DAG constraint terms in Eq. ( or (5). However, if the length of the longest simple path in the graph is known, we can formulate lower-order polynomial based DAG constraint.
Proposition 2. Let B E Rdxd be the weighted adjacency matrix of a graph 9 with d vertices, and k be the length of the longest simple pattl in graph G. The following three conditions are equivalent: 0, and (3) byranc (B) O.
(I) G is DAG, B k+1
The intuition of the above proposition is follows. For the graph whose largest simple cycle has length of k, it suffices to compute the :t power of its adjacency matrix, since computing any power larger than simply traveling in the cycles repeatedly. Therefore, Proposition suggests that bjuee (B) is valid DAG constraint.
Unfortunately, the longest simple path problem for general graphs is known to be NP-hard 2 and thus we have to find an alternative way to seek the proper k. For a candidate adjacency matrix B, a simple heuristic to zero matrix. Recall that Proposition T Iiplies that the adjacency matrix of a DAG is nilpotent; therefore, when a candidate adjacency matrix B is close to being acyclic, its spectral radius would be close to zero, which implies that the B would converge quickly to zero as increases. Thus, if we the above heuristic to find k; it is highly possible that find some < d.
iterate over the power of matrices to find some s.t. matrix B is close enough
By using an approximate & instead of the exact one, the DAG constraint Eq. becomes an approximate one with the following bounds. Denoting by llso the elementwise infinity norm (i.e., maximum norm), we have the following proposition.
Proposition 3. Given B Rdddd there exists l < d such that || 1l00 < (k+1)a then we have
1 (de)d/k-1 2+1/ke141/R
< lgoo (B) 1 (B) < 1- --dd)1//
0
V h geo B
V h M [B) < (k + 1)dell V h (B)lle-
Remark 1. The condition < [k+1jd in the proposition above such that the inequalities are
non-vacuous, specifically since we must have V Bg (B) - V (B)ll / h (B)lFF < 1.
We use simple path to refer to path that repeats no vertex, except that the first and last may be the same vertex 17 363].
5

---

Proposition 3 provides an efficient way to bound the enror of truncated constraints h (B) to the geometric series based DAG constraints (B). However, it also demonstrates possible limitation for all polynomial based DAG constraints. As we know, when a solution is close to DAG, its spectral radius will be close to zero. In this situation, if the longest simple path in a graph is long, it is unavoidable to use high-order terms to capture this information while such high-order terms might be close to O. In the worst case, the higher-order terms may underflow due to limited machine precision. Thus some specific numerically stable method may be required, and we leave this for future work.
3.3 Dynamic Truncated Matrix Power Iteration
Proposition B indicates that it may be safe to discard those matrix powers with index larger than k. In
this case, the DAG constraint and its gradient correspond to order-k and order-(k: 1) polynomials of B, respectively. straightforward way to compute the DAG constraint as well as the gradient is provided in Algorithm[[ which call TMPI. The general polynomial based DAG constraint in 5 requires (dd) matrix multiplication, while Algorithm Saatumm@reuuce reduces the number of multiplication to ( In practice, & can be much smaller than d; for example, when recovering 100-node Erdös-Réón DAG with expected degree 4, k is often smaller than 20.
Algorithm Fast TMPI
Algorithm Truncated Matrix Power Iteration (TMPI)
Input: B Rexd 2:0:6
Input: B E Em
Output: k(B). BBh(B)
1: itte
Output: ((B). VB)(B)
By < B, t- B, <- B (By stores S.BB), B,
1: i
2,b(B) tr(B). Vah(BB 4-1
stores V (gtr(f((B)). and B, stores B')
while & do
3: while & do
3:
Vgh(B) + V//(BBB (B) + +i (T
B, + By, B old + Bg B + Bp
4:
4:
- BB 1-1
5:
B, B Bp B
old old
B, B, B
old old
5:
h(B) t h(B) + tr
6:
if [ ||xB then break
6:
B, + B +B Bg +i(B, B, B - B,)
old
old
i+i+1
7:
if |Bplls & then break
end while
8:
i i 2
9: Output ((B). (B)
9:
end while
10: Output ir LB, B,
Fast Truncated Matrix Power Iteration We exploit the specific structure of geometric series to further accelerate the TMPI algorithm from (kk Ollog k)) matrix multiplications. This is acheived by leveraging the following recurrence relations, with a proof provided in Appendix A.4
Proposition 4. Given any d d real matrix B, let [.(B) B+ +....++++.............. +B' and h (B) = tr (r()(()). Then we have the following recurrence relations:
fiz = f.B (B),
V (B) h, T ( (B
This straightfowardly implies the following corollary, with proof given in Appendix A.5
Corollary 1. Given any matrix d x d real matrix B, let fi = B +
B and h (B)
tr():()). Then we have have the following recurrence relations:
faa(B) (I) (B)+
V (B) -V BBB + TVBT +i ( (BB) +B² B B²)
Instead of increasing the index i (i.e., power of B) arithmetically by 1, in line 9 of Algorithm I Corollary I suggests using its recurrence relations to increase the index geometrically by factor of 2, until the condition | |ox < is satisfied. The full procedure named Fast TMPI is described in Algorithm which requires Ollog k) matrix multiplications and O(d2) additional storage. An

---

Similar strategy can be applied to the original geometric series based constraint (baj) form a fast algorithm that requires Olog d) matrix multiplications. A comparison of the time complexity between our proposed constraint and the existing ones is provided in Appendix B.2
example implementation is given in Listing T in Appendix B.1
The Full Optimization Framework For NOTEARS, we apply the same augmented Lagrangian to convert the constrained optimization problem into
framework 4 5 adopted by Zheng al.
500
the following unconstrained problem
min BER
X - X
P
+ B
f 2 trunc (B B)
aah
Irune BBB),
where the parameters p and a are iteratively updated using the same strategy as Zheng et al. 50) We also apply the proposed DAG constraint to GOLEM [3 that adopts likelihood-basse objective with soft constraints, leading to the unconstrained optimization problem
d 2
min Bg
log X X-XBI log det(I B) + m B + 12 &B
(10)
where 71, 1/2 are the regularization coefficients. In both formulations above, we use Algorithm | or dynamically compute Mrme (B @ B) (and the value of k)) during each optimization iteration based on predefined tolerance e.
3.4 Relation to Previous Works
Our algorithm closely related the following DAG constraint considered by Zheng al. (50
Lemma (Zheng et al. [50 Proposition 1J). Ler B iaxd with spectral radius smaller than one be the weighted adjacency matrix of graph 9 with d vertices. a DAG if and only if
tr(1 B =
(11)
When the spectral radius is smaller than one, the inverse of l1 B is simply the following infinite geometric series:
8
(1-B) -
Σ B
(12)
Zheng et al. 150 argued that 12 is not well defined for B with spectral radius larger than and in practice the series 12 may explode very quickly. Based on this argument, they adopted the matrix exponential for constructing the DAG constraint. For the same reason, Yu et al. small coefficients for higher-order terms in polynomial based DAG constraints. However, as we illustrate in Section the DAG constraints proposed by Yu al. 45 Zheng al. (50) may suffer from gradient vanishing owing to the small coefficients for higher-order terms. In fact, our constraint 7 serves good approximation for (IT since we are actually looking for an approximate inverse of 1 B with bounded error.
[55] also used very
Proposition 5. Given B E Rdxd there exists such that |
1lo0 < < 1/d, then we have
(I-
B)
< de
[I B B)-
1=
Previous works 24 pointed out that for DAG constraints with form h(B B), the Hadamard product would lead to vanishing gradients in the sense that the gradient VBh(B B) 0 if and only if (BB @ B) = 0. Ng et al. [24] further showed that this may lead to ill-conditioning issue when solving the constrained optimization problem. One may replace the Hadamard product with absolute value, as in Wei et al. 44] but the non- monsmmottheess of absolute value may make the optimization problem difficult solve and lead to poorer performance 24 As pointed out by Ng et al. 24]. Quasi- Newtown methods, e.g., L-BFGS 77]. still perform well even when faced with the ill ill-condtiioning issue caused by vanishing gradients of constraints with form A(B c B). this work, we identify another type of gradient vanishing issue caused by (1) the higher-order terms and small coefficients function itself and (2) the small spectral radius of candidate matrix and further provide an effective solution by constructing a Truncated Matrix Power Iteration based DAG constraint.
7

---

4
Experiments
We compare the performance of different DAG learning algorithms and DAG constraints to demon- strate the effectiveness of the proposed TMPI based DAG constraint. We provide the main experi- mental results in this section, and further empirical studies can be found in Appendix
4.1
DAG Learning
We first conduct experiments on synthetic datasets using similar settings as previous works [23 46 50
In summary, random Erdös-Réng Rényi graphs [14] are generated with d nodes and kd expected edges (denoted by ERk). Edge weights generated from uniform distribution over the union of two intervals [-2, -0.5 U [0.5, 2.0] are assigned to each edge to form a weighted adjacency matrix B. Then n samples are generated from the linear SEM x B + to form an n x d data matrix X, where the noise e are i.i.d. sampled from Gaussian, Exponential, or Gumbel distribution. We set the sample size n 1000 and consider 4 different number of nodes d = 10, 30, 50, 100 unless otherwise stated. For each setting, we conduct 100 random simulations to obtain an average performance.
Methods and Baselines We consider continuous DAG learning methods based on both Mean Squared Error (MSE) loss (see Eq. (9) and likelihood-based (see Eq. (TO) loss from GOLEM [ For both losses, we include the following DAG constraints: (1) MSE/Likelihood loss+Exponential DAG constraint G i.e. NOTEARS [50/GOLEM 23 (2) MSE/Likelihood loss+Binomial DAG constraint 4 from DAG-GNN :55::
: (3) MSE/LLkellhhodd Spsctrl radius based DAG constraint from NOBEARS e1: (4) MSE/Likelihood Single Term DAG constraint R (5) MSE/Likeli hood loss+ Truncated Matrix Power Iteration (TMPI) based DAG constraint 7 For the NOCURL method
46 since it is a two stage methods where the initialization stage rely MSE loss and Binomial DAG constraint, we include it as baseline, but do not combine it with likelihood loss. The other baselines include GES M
and MMHC 41 42
GES
MMHC
NOCURL
MSE Single_tem
-1- MSE Exponential
MSE Binomial
MSE NOBEARS
MSE TMPI(Ours)
20
10a
SHD
SHD
10
10a
10a
30
50
100
100
10
50
100
Nodess
#Nodes
#Nodes
1.0
10a
100
SHD0
SHD
SHD
101
10
10
30
50
100 10
30
100
10
30
50
100
*Nodes
*Nodes
*Nodes
Gaussian
Exponential
Gumbel
Figure 2: DAG learning results of MSE loss (e.g. Eq. ( based algorithms and baselines in terms of SHD (lower is better) on ER2 (top) and ER3 (bottom) graphs, where our algorithm consistently outperforms others almost by an order of magnitude. Error bars represent standard errors over 100 simulations.
For MSE loss based DAG learning methods 9. we set 1) = 0.1 and use the same strategy as Zheng et al. (50) to update the parameters P and a. For likelihood loss based DAG learning methods IO), we set 7)1 0.02 and 7/2 = 5.0 Ng et al. [ For the Binomial DAG constraints 4 the parameter a is set to 1/d as Yu et al. [5 For our TMPI DAG constraint, we set the parameter
e 10-6 6 For the NOCURL method, we use their public available implementation and default
Thh source code is avaiahle at here
'The implementation is provided here
8

---

parameter setting
For GES and MMPC, we use the implementation and default parameter setting in
the CalsalDiscovvrr/oooooooo package [6
GES
MMHC
NOCURL
GOLEM Single_term
GOLEM Exponential
GOLEM Binomial
GOLEM NOBEARS
GOLEM TMPI(OurS)
10
10
10a
SHD
10
SID
101
SHD
101
10a
10"
100
10"1
10
30
50
100
50
100
10
50
100
#Nodes
10a
10a
10a
SHD
SHD
5HD
101
10
109
10a
10a
10
30
50
100
10
30
50
100
10
30
50
100
#Nodes
#Nodes
#Nodes
Gaussian
Exponential
Gumbel
Figure 3: DAG learning results of likelihood loss (i.e. Eq. ITC based algorithms and baselines terms of SHD (lower is better) on ER2 (top) and ER3 (bottom) graphs, where our algorithu consistently outperforms others in a large margin in almost all cases. Error bars represent standard errors over 100 simulations.
Results The results for MSE/likelihood based DAG learning algorithm are shown in Figure where we compare the Structural Hamming Distance (SHD) of different algorithms. The complete re- sults for other metrics reported in Appendix CA The original NOTEARS (i.e. MSE+-Exponental) / GOLEM (i.e. DOLEM+EExoneettall has similar performance as MSE/GOLEM+Binon may not be surprising because the coefficients of each order term in the Exponential and Binomial constraints are actually very close (see Figure | for an example). The performance of algorithms using single-term DAG constraint often not comparable the other polynomial-baas based constraints since it suffers more from gradient vanishing. The performance of algorithms using the spectral radius based constraints (BEARS) are also worse than the ones using Exponential and Binomial constraints. The algorithms using our TMPI constraints outperforms all others in almost all cases, which demonstrate the effectiveness of our DAG constraint.
This
Additional Experiments We conduct additional experiments with scale-free graphs, smaller sample size (i.e,, 200 samples), larger graphs (i.e., 500-node graphs), and denser graphs (i.e., ER6 graphs). Due to the space limit, the results of these experiments are provided in Appendices C.2 and C.6 respectively, which further demonstrate the effectiveness of our proposed DAG constraint. We also compare our method with NOFEARS 4A in Appendix
4.2 Extension to Nonlinear Case
Although the description of our algorithm in Section 3 focuses on the linear case, it can be straight- forwardly extended to the nonlinear case. In this section, we extend our DAG constraint to nonlinear case by replacing the binomial DAG constraint in DAG-GNN 45 NOTEARS-MIL [5] and GRAN-DAG 20 with our TMPI constraint.
Synthetic Data
We compare our algorithm and DAG-GNN on 50-node ER graph with nonlinear SEM x B cos(x + e over 5 simulations. Our algorithm achieves SHD of 22.2 + 4.2 while DAG-GNN achieves SHD of 25.2 + 4.5. We simulate nonlinear dataset with 50-node ER I graphs and each function being an MLP as Zheng et al. [51) The original NOTEARS-MLP has SHD 16.9 + 1.5, while COT with our TMPI constraint achieves an SHD of 14.9 + 1.3.
https: thub. com/ www.sssmoon200mm 123///43--4-A00 mit 4e2alf4 on 14 Jun 2021. Better performance may be attained by tuning the hyper-parameters, but the current best reported performance slightly worse than NOTEARS in most settings.
http////ommmaaeeeereeeee http://maaaaonffffaaaaaeemeemmee efficient py from com-

---

Real Data We also run an experiment on the Protein Signaling Networks dataset B1 consisting of n 7466 samples with d = 11 nodes. The dataset is widely used to evaluate the performance of DAG learning algorithms: however, since the identifiability of the DAG is not guaranteed, we compare both SHD and SHDC (SHD of CPDAG), following Lachapelle et al. [20]. this experiment, we replace the original DAG constraints in DAG-GNN ES] and Gran-DAG 20 with our TMPI DAG
constraints. The results are shown in Table II showing that by replacing the DAG constraint with our TMPI constraint, there is a performance improvement on SHD and SHDC for almost all cases.
4.3 Ablation Study
We use the ER3 dataset in previous section to conduct an ablation study to compare the performance of the Geometric constraint (6a) and our TMPI constraint 7 with both the naive implementation
(i.e., linear complexity) and fast implementation (i.e., logarithmic complexity) from Algorithms I and 2 respectively. The results are shown in Figure + The SHD between Geometric constraint and TMPI constraint in both the naive and fast implementations are similar, and our Fast TMPI algorithm runs significantly faster than the other three, especially when the graph size is large. We also conduct further empirical study to illustrate the numeric difference between different k for our DAG constraint; due to limited space, detailed discussion is provided in Appendix Ci
+ 100
EMP
DOLLMM Clesnatric
-+- DOLEX IMP
GOLEM Fast
Table I:
Experimental results on Sachs dataset. Our TMPI DAG constraint consis- tently outperforms other DAG constraints.
60-
10
SHD
SHDC
60
10"
DAG-GNN
16
21
20-
DAG-GNN/IMPI (Ours)
17
Gran DAG
13
11
10
30
104ada
50
1.0
WWodd
160
Gran: DAG/TMPI (Ouns)
Gran-DAG++
Figure 4: Comparison of Geometric and TMPI constraints, all algorithm achieves similar performance but fast algorithm runs significantly faster. Left: Running time; Right: SHD, Error bars represent standard errors over 100 simulations.
Gram DTGH+ (Ours)
12
5 Conclusion and Future Work
While polynomial constraints- DAG learning algorithms achieve promising performance, they often suffer from gradient vanishing in practice. We identified one important source of gradient vanishing: improper small coefficients for higher-order terms, and proposed to use the geometric series based DAG constraints to escape from gradient vanishing. We further proposed efficient algorithm that can evaluate the constraint with bounded error in Ollog k)) complexity, where is often far smaller than the number of nodes. By replacing the previous DAG constraints with ours, the performance of DAG learning can often be significantly improved in various settings.
A clear next step would be finding strategies to adaptively adjust the coefficients for higher-order terms in polynomial based constraints. In a nutshell, polynomial- based constraints may explode on candidate graphs that far from acyclic, and also suffer from gradient vanishing problem when the graph is nearly acyclic. Fixed coefficients are not able to solve both problems at the same time. Another possible direction is use the sparse properties in DAGS make DAG learning more efficient. Since the adjacency matrix of DAG must be nilpotent, there might exist more compact representation of the adjacency matrix, such low-rank CANDECOMPIPARAFAA (CP) decomposition, which can potentially be used to derive an efficient algorithm for very large graphs.
Acknowledgements
We thank the anonymous reviewers, especially Reviewer iGwL, for their helpful comments. ZZ, DG, YL, EA and JQS were partially supported by Centre for Augmented Reasoning at the Australian Institute for Machine Learning. IN and KZ were partially supported by the National Institutes of Health (NIH) under Contract R01 by the NSF-Convergence Accelerator Track-D award #2134901, by grant from Apple Inc., and by grant from KDDI Research Inc.. MG was supported by ARC DE2 10101624.
10

---

References
[1] A. Al-Mohy and N. Higham. A new scaling and squaring algorithm for the matrix exponential. SIAM Journal on Matrix Analysis and Applications, 31, 2009.
[2] P. Bader, S. Blanes, and F. Casas. Computing the matrix exponential with an optimized taylor polynomial approximation. Mathematics, 7(12):1174, 2019.
[3] A.-L. Barabási and R. Albert. Emergence of scaling in random networks. Science, 286(5439): 509-512, 1999.
[4] D. P. Bertsekas. Constrained Optimization and Lagrange Multiplier Methods. Academic Press, 1982.
[5] D. Bertsekas. Nonlinear Programming. Athena Scientific, 2nd edition, 1999.
[6] R. Bhattacharya, T. Nagarajan, D. Malinsky, and I. Shpitser. Differentiable causal discovery under unmeasured confounding. In International Conference on Artificial Intelligence and Statistics, 2021.
[7] R. H. Byrd, P. Lu, Nocedal, and C. Zhu. A limited memory algorithm for bound constrained optimization. SIAM Journal on scientific computing, 16(5):1190-1208, 1995.
[8] D. M. Chickering. Learning Bayesian networks NP-complete. In Learning from Data: Artificial Intelligence and Statistics V. Springer, 1996.
[9] D. M. Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3(Nov):507-554, 2002.
[10] M. Chickering, D. Heckerman, and C. Meek. Large-sample learning of bayesian networks is np-hard. Journal of Machine Learning Research, 5, 2004.
[11] D. Colombo, M. H. Maathuis, M. Kalisch, and T. S. Richardson. Learning high-dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics, pages 2012.
[12] D. Coppersmith and S. Winograd. Matrix multiplication via arithmetic progressions. Journal of Symbolic Computation, 93):251-2200 1990. ISSN 0747-7171. Computational algebraic complexity editorial.
[13] J. Cussens. Bayesian network learning with cutting planes. In Conference on Uncertainty in Artificial Intelligence, 2011.
[14] P. Erdös and A. Rényi. On random graphs I. Publicationes Mathematicae, 6:290-297 297, 1959.
[15] R. Fletcher. Practical Methods of Optimization.. Wiley-Interscience 1987.
[16] D. Kalainathan and O. Goudet. Causal discovery toolbox: Uncover causal relationships in
python. arXiv preprint (oXXv:9903.022777 2019.
[17] D. E. Knuth. The Art of Computer Programming, Volume I: (3rd Ed.) Sorting and Searching. 1997.
[18) M. Koivisto and K. Sood. Exact Bayesian structure discovery in Bayesian networks. Journal of Machine Learning Research, 5(Dec):549-51) 2004.
[19] D. Koller and N. Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
[20] Lachapelle, Brouillard, Deleu, and S. Lacoste-Julien. Gradient- neural DAG learning. In International Conference on Learning Representations, 2020.
[21] H.-C Lee, M. Danieletto, R. Miotto, S. T. Cherng, and J. T. Dudley. Scaling structural learning with no-bears infer causal transcriptome networks. In PACIFIC SYMPOSIUM ON BIOCOMPUTING 2020, pages 391-402. World Scientific, 2019.
[22] I.. Loh and P. Bühlmann. High-dimensional learning of linear causal networks via inverse covariance estimation. Journal of Machine Learning Research, 2013.
[23] I. Ng, A. Ghassami, and K. Zhang. On the role of sparsity and DAG constraints for learning linear DAGs. Advances in Neural Information Processing Systems, 33, 2020.
[24] I. Ng. S. Lachapelle, N. R. Ke, S. Lacoste-Julien, and K. Zhang. On the convergence of continuous constrained optimization for structure learning. International Conference on Artificial Intelligence and Statistics, 2022.

---

[25] I. Ng, S. Zhu, Z. Fang, H. Li, Z. Chen, and J. Wang. Masked gradient- gradient-based causal structure learning. In SIAM International Conference on Data Mining, 2022.
[26] R. Pamfil, N. Sriwattanaworachii, S. Desai, P. Pilgerstorfer, P. Beaumont, K. Georgatzis, and B. Aragam. DYNOTEARS: Structure learning from time-series data. In International Conference on Artificial Intelligence and Statistics, 2020.
[27] J. Pearl. Models, reasoning and inference. Cambridge, UK: Cambridge University Press, 19, 2000.
[28] Peters and P. Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, sequuniia unconstrained minimization techniques.
228, 2013.
[29] M. D. Powell. Nonlinear
The Computer Journal, 12(3). 1969.
[30] G. Raskutti and C. Uhler. Learning directed acyclic graph models based on sparsest permutations. Stat, 7(1)= 2018.
[31] K. Sachs, O. Perez, D. Pe'er, D. A. Lauffenburger, and G. Nolan. Causal -ignaaingg networks derived from multiparameter single-cell data. Science, 308(5721):523-5299 2005.
[32] A. Schrijver. Combinatorial optimization: polyhedra and efficiency, volume 24. Springer Science & Business Media, 2003.
[33] R. D. Shah and J. Peters. The hardness of conditional independence testing and the generalised covariance measure. The Annals of Statistics, 48(3):1514-1538 2020.
[34] S. Shimizu, P. O. Hoyer, A. Hyvärinen, and A. Kerminen. A linear non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(Oct):2003 2030,, 2006.
[35] T. Silander and P. Myllymäki. A simple approach for finding the globally optimal Bayesian network structure. In Conference on Uncertainty in Artificial Intelligence, 2006.
[36] A. P. Singh and A. W. Moore. Finding optimal Bayesian networks by dynamic programming. Technical report, Carnegie Mellon University, 2005.
[37] Spirtes and C. Glymour. An algorithm for fast recovery of sparse causal graphs. Social science computer review, 91):622-7,, 1991.
[38] P. Spirtes, C. Meek, and T. Richardson. Causal inference in the presence of latent variables and selection bias. In Conference on Uncertainty in Artificial Intelligence, 1995.
[39] P. Spirtes, C. N. Glymour, R. Scheines, and D. Heckerman. Causation, prediction, and search. MIT press, 2000.
[40] V. Strassen. Gaussian elimination is not optimal. Numerische Mathematik, 13:354-356, 1969.
[41]
Tsamardinos, C. F. Aliferis, and A. Statnikov. Time and sample efficient discovery of markoy blankets and direct causal relations. In Proceedings of the ninth ACM SIGKDD international conference on Knowledge discovery and data mining, pages 673-678, 2003.
[42] I. Tsamardinos, L. E. Brown, and C. F. Aliferis. The max hill-climbing bayesian network structure learning algorithm. Machine learning, 65(1):31-78, 2006.
[43] M. J. Vowels, N. C. Camgoz, and R. Bowden. D'ya like dags? survey on structure learning and causal discovery. arXiv preprint N.Niv:2103.02.:: 2021.
[44] D. Wei, T. Gao, and Y. Yu. DAGS with no fears: A closer look at continuous optimization for learning bayesian networks. arXiv preprint arXiv:2010.09/33 133, 2020.
[45] Y. Yu, J. Chen, T. Gao, and M. Yu. DAG-GNN: Dag structure learning with graph neural networks. In International Conference Machine Learning, pages 7154-7163, PMLR, 2019.
[46] Y. Yu, T. Gao, N. Yin, and Q. Ji. DAGs with no curl: An efficient dag structure learning approach. In Proceedings of the 38th International Conference on Machine Learning, 2021.
[47] C. Yuan and B. Malone. Learning optimal Bayesian networks: A shortest path perspective. Journal of Artificial Intelliencee Research, 48(1):23 65, 2013.
[48] C. Yuan, B. Malone, and X. Wu. Learning optimal Bayesian networks using A* search. In International Joint Conference Artificial Intelligence, 2011.
12

---

[49] B. Zhang, C. Gaiteri, L.-G. Bodea, Z. Wang, J. McElwee, A. A. Podtelezhnikov, C. Zhang, Xie, L. Tran, R. Dobrin, al. Integrated systems approach identifies genetic nodes and networks in late-onset alzheimer's disease. Cell, 153(3):707-720, 2013.
[50] X. Zheng, B. Aragam, P. Ravikumar, and E. P. Xing. DAGs with No Tears: Continuous optimization for structure learning. arXiv preprint AKM3000203 2018.
[51] X. Zheng, C. Dan, B. Aragam, P. Ravikumar, and E. Xing. Learning sparse nonparametric dags. In International Conference on Artificial Intelligence and Statistics, pages 3414-3425. PMLR, 2020.
[52) R. 7hu, A. Pfadler, 7. Wu, Y. Han, X. Yang, F. Ye, Z. Qian, J. Zhou, and B. Cui. Efficient and scalable structure learning for bayesian networks: Algorithms and applications. In 202/ IEEE 37th International Conference on Data Engineering (ICDE), pages 2613-2624. IEEE, 2021.
[53] S. Zhu, I. Ng. and Z. Chen. Causal discovery with reinforcement learning, In International Conference on Learning Representations, 2020.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [No]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes] See Appendices
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL)? [Yes]
(b) Did you specify all the training details (e.g.. data splits, hyperparameters, how they were chosen)? [Yes] See Section | and Appendices
(c) Did you report error bars (c.g., with respect to the random seed after running experi-
ments multiple times)? [Yes] See Section
4
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Appendices
4. If you are using existing assets (e.g., code, data, models) or curating/releusing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license of the assets? [Yes] See Appendices
(c) Did you include any new assets either in the supplemental material or as URL? [No]
(d) Did you discuss whether and how consent was obtained from people whose data you're 're're
using/curating? [Yes] We use publicly available datasets
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "pSU7t0VGTB_", "reviewer": "Reviewer_XcDH", "summary": "This is an interesting paper that seeks to provide new non-cyclic constraint formulations for DAG learning. The key novelty here is the employment of finite geometric series DAG constraint to escape from gradient vanishing . The idea of reducing the polynomial constraint to order-k polynomial also provides some very interesting benefits. The experimental comparison with SoTA methods is comprehensive, although limited to the comparison on accuracy. Overall, the work is well done and provides some new methodology to this field.", "questions": "1. For the experimental results, besides SHD and SHDC, the number of learnt edges should also be provided.\n2. In numerical tests, the methods were tested only on relatively sparse graphs (with <=3 averaged number of edges). However, it is often more challenging to learn the graph structure on denser graphs. Can the authors provide comparison on at least one denser graph case, say, ER6?\n3. Besides SHD, can the authors provide a comparison on computational time, to provide evidence on the claimed efficiency?\n4. NOFEARS has also claimed to resolve the vanishing gradient issue. Can the authors compare their method with NOFEARS?\n5. According to what I read, it seems DAG no curls aim to accelerate the learning, not to improve the accuracy. Moreover, its accuracy depends on the tunable penalty parameter, which is possibly why the authors observed a deteriorated performance here if they didn't optimize that penalty parameter. This fact should be at least pointed out and discussed in the paper.", "limitations": "The limitations were discussed adequately.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: \nA new type of DAG constraint is proposed, and analysis is provided. Empirical results show that this new constraint helps to improve the learning accuracy, in terms of SHD.\n\nWeakness:\nGiven that the possibility of using general polynomial constraints was already proposed and proved in NOFEARS, and they have also proposed a new formulation to resolve the vanishing gradient issue, the analysis in this paper seems incremental. One of the main contribution is an efficient algorithm to evaluate the constraint. However, only the comparison on SHD was provided. Therefore, there is not presented evidence that the proposed method is more efficient than the SoTA methods.", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "SXwTFqbOc-1", "reviewer": "Reviewer_SUap", "summary": "This paper proposes a new regularization term for DAG estimation, whose role is to enforce that the connectivity matrix is nilpotent.\n\nThe contribution basically amounts to depart from state of the art that considers infinite sums of power terms as a regularization to replace them by a finite sum of un-weighted power terms. \nThis simple change is motivated in a thorough theoretical way and a fast algorithm is proposed. Performance is good.", "questions": "* I am not super knowledgeable about DAG, so that I may be missing something trivial, but I don't understand proposition 2 (3). I do understand that \\tilde{B}^{k+1}=0, but I don't really understand why we should have tr(\\sum_i \\tilde{B}^i)=0. Your explanation \"Therefore, Proposition 2 suggests that h_{trunc}^k(\\tilde{B})=0 is a valid DAG constraint\" is not illuminating in this respect. Please elaborate a bit more for non-knowledgeable readers like me.\n* In \"The full optimization framework\", I don't understand why you are picking B\\odot B. It's probably obvious but once again, I need a sentence on that.\n* I need more details on the nonlinear real case. Could you please detail that better ? Notably the model and how the algorithm is modified.\n* L99: \"give\"-> given\n* L138: \"attains far more better\" far better\n* L225: \"we are actually find an\": looking for\n* L306: \"while [...] moderate\" reads weird because \"moderate\" is a negative adjective. Replacing \"moderate\"-> \"promising\" makes the sentence closer to what you probably mean (basically \"although the previous stuff is good, it suffers from some weakness and fix that\")\n* L314: \"nut-shell\"->nutshell", "limitations": "models able to handle bigger data are mentioned at the end, but as I wrote it would be good to show that the proposed method works fine on moderate but bigger graphs.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 4, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "* The paper focuses on the very focused problem of estimating a DAG connectivity matrix. This highly focused contribution is a strength in my view: it stays focused and everything about the paper looks good in my view, from the state of the art to the theoretical developments.\n* As another strength, I must mention that the whole paper is very factual and written in a neutral tone. For once, people are not pretending they are \"reinventing\" or \"rethinking\" something...\n\n* As a weakness, I could mention the limited scale of the experiments. I understand that the paper is mostly theoretical, but I guess that DAG estimation would be an important topic for much larger networks. Maybe going to millions is out of scope, but doing experiments with intermediate scales like thousands would be a desirable feature.", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "38Eqj9BVbdZ", "reviewer": "Reviewer_CoNX", "summary": "This paper consider continuous score-based DAG learning optimization problem, which is optimizing a least square objective with continuous constraint that characterizes acyclicity for directed graphs. The paper recognizes the gradient vanishing is a important issue in this optimization problem, and most of literature have this issue since they put small weights on the higher order terms in the constraint. The paper proposes a new constraint to avoid gradient vanishing with a trucation to accelerate and ensure stability. Some experiments are conducted to show the effectiveness. ", "questions": "- The paper claims that the gradient vanishing is a key challenge for continuous score-based DAG learning, while current literature all put small weights on higher order terms in the constraint, thus share this problem. Based on the paper, Fig 1 only shows TMPI constraint can control the gradient of higher order terms, but it is still not clear why and how gradient vanishing is important. Does it hinder the optimization process? convergence of algorithm? or statistical property of the estimator?\n\n- If the improment of TMPI is on optimization side, it would be interesting to see the runtime comparison among different methods. And it would be better to provide discussion or intuition why TMPI performs much better in terms of graph recovery (SHD) than other methods while they are essentially optimizing the same least square objective under DAG constraint.\n\n-  How is the proposed constraint related to the ones other than NOTEARS in existing work? Does it also resolve the issues recognized by others (not gradient vanishing)? If not, is it obvious how they can be combined?", "limitations": "One limitation is discussed right after *Remark 1* and also in *Conclusion*, about all polynomial based DAG constraints, which is when the longest path in the graph is long, the higher order term in constraint matters. It is interesting to see discussion on some particular examples, like Markov chain, especially when the coefficients $B_{jk} > 1$.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strength\n\n- Solving continuous score-based DAG learning optimization problem is a topic with significance;\n- Lots of experiments are conducted to show the effectiveness of proposed new constraint compared to existing work;\n- Smooth presentation and discussion about this thread of continuous DAG learning literature.\n\n# Weakness\n\n- It is unclear why gradient vanishing is important, and comparison with existing work is not detailed enough, see Questions;\n- Most experiments are about graph recovery, since this is an optimization paper, runtime comparison with other methods is more interesting.", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "0_qw1gPqkW", "reviewer": "Reviewer_iGwL", "summary": "This paper is on learning directed acyclic graphical models using a continuous optimization formulation, where acyclicity is enforced through a matrix polynomial constraint. It is shown that previously proposed polynomial constraints (NOTEARS and similar) suffer from vanishing higher-order terms. A polynomial with unit coefficients is proposed to address this problem, as well as a truncated, lower-order version of this polynomial. The approximation error of the truncated polynomial is bounded. Furthermore, a fast algorithm is presented to compute the new polynomial constraint as well as its gradient. Experiments are performed on synthetic data generated from linear and nonlinear structural equation models, as well as one real dataset. The new polynomial constraint can reduce structural Hamming distance by large factors (i.e. 2-4X) over existing polynomial constraints.", "questions": "- Propositions 3 and 5: Which norms are referred to by $\\lVert \\cdot \\rVert_\\infty$ and $\\lVert \\cdot \\rVert_2$?\n- Proposition 1: Is the reference to [37] in the proof actually to [35]? Is this proposition a corollary of [35, Lemma 1]?\n- Figure 1 and related discussion: Clarify what \"gradient\" and \"gradient on higher-order terms\" mean.\n\nMore for curiosity:\n- The SHD of all algorithms, and in particular the strong baselines Exponential and Binomial, are worse by a large factor for the MSE loss in Figure 2 compared to the likelihood loss in Figure 3. The SHDs in Figure 2 are also worse by a similarly large factor compared to previous work with a similar experimental setting. Compare the top row (ER2) in Figure 2 to e.g. the leftmost columns in Figures 1, 5, 6, 7 in [35], which are also for ER2, $n = 1000$, the same values of $d$, and the same noise types. Any idea why?", "limitations": "Satisfactory", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**\n+ Furthers our knowledge about continuous constraints for enforcing acyclicity, by pointing out a shortcoming of previously proposed polynomial constraints and proposing a different polynomial to remedy the problem.\n+ The structural Hamming distance improvements in the experiments are large (factors of 2-5).\n+ The degree $k$ of the polynomial is adapted to each candidate adjacency matrix $\\tilde{B}$.\n\n**Weaknesses:**\nI like the paper especially for the first reason above. However, the current version suffers from several technical issues in my view.\n- Propositions 3 and 5: There appear to be several problems.\n    1. First, which norms are meant by the notation $\\lVert \\cdot \\rVert_\\infty$ and $\\lVert \\cdot \\rVert_2$? Previously, $\\lVert \\cdot \\rVert_p$ was defined as the element-wise $p$-norm but I suspect this is no longer the case here.\n    2. Based on the proof in Appendix A.3, I will assume that $\\lVert \\cdot \\rVert_p$ now means the induced operator norm, i.e., $\\lVert A \\rVert_p = \\sup_{\\lVert x \\rVert_p = 1} \\lVert Ax \\rVert_p$. I am assuming this because in the unnumbered equation above (14), I believe $\\lVert A \\rVert_\\infty$ does bound the spectral radius $\\rho(A) = \\max_i \\lvert \\lambda_i(A) \\rvert$ from above. However, on the left-hand side of this inequality and in (14), the induced 2-norm $\\lVert A \\rVert_2$ does not equal $\\rho(A)$ for a general non-symmetric matrix (the Wikipedia article on spectral radius has a simple counterexample). The adjacency matrix $\\tilde{B}$ and its powers are almost certainly non-symmetric.\n    3. Continuing with the assumption that $\\lVert A \\rVert_2$ is the induced 2-norm, the same error occurs in the second-last line of Appendix A.3. On the other hand, I don't think we can replace $\\lVert \\cdot \\rVert_2$ by $\\rho(\\cdot)$ as we need the sub-additivity property of norms to go from the third-last line to the second-last line, and the spectral radius is not a norm.\n    4. Proposition 5: Again assuming that $\\lVert A \\rVert_2$ is the induced 2-norm, the same error occurs in the proof in Appendix A.5.\n    5. Remark 1: The big-O notation seems to be ignoring the $1 - \\epsilon^{1/k}$ in the denominators. I don't think this can be done as $k$ can be fairly large so $\\epsilon^{1/k}$ is close to 1. Moreover, $k$ is not a constant as it depends on $\\epsilon$.\n- Proposition 2: I think the proof of (3) $\\to$ (1) is incomplete. Since it is still to be proven that $\\mathcal{G}$ is a DAG, one must still address directed paths of length greater than $k$, and probably argue that they cannot exist. Also, the phrases \"shorter than $k$\" and \"less than $k$\" should be \"less than or equal to.\"\n- Proposition 1: The proof of the \"if\" part in Appendix A.1 refers to [37], but I could not find any mention of \"nilpotent\" in the paper by Yuan and Malone (2013). Perhaps the reference is to [35] (Wei et al., 2020), wherein Lemma 1 states that $\\mathcal{G}$ is acyclic if and only if $\\tilde{B}$ is nilpotent. Indeed, Proposition 1 is a corollary of [35, Lemma 1] since the $\\tilde{B}^d = 0$/nilpotent $\\to$ DAG direction is exactly the same. For the other direction, [35, Lemma 1] implies that $\\tilde{B}$ must be nilpotent, and if it is, then $\\tilde{B}^k = 0$ for a $k \\leq d$, hence $\\tilde{B}^d = 0$ also. If the authors agree with this, then the implication from [35, Lemma 1] needs to be acknowledged.\n- Figure 1 and corresponding discussion: \n    - It is not clear what gradient is referred to here and what \"gradient on higher-order terms\" means. I am guessing it is the gradient of the constraint function with respect to $\\tilde{B}$ and the \"higher-order terms in the gradient\" refer to the terms in (7) and similar expressions for the gradient of (3), (4). Perhaps it would be better to show the general form of the gradient and then specialize it to the different cases.\n     - The statement \"nilpotent properties of DAG result in candidate adjacency matrix with very small spectral radius\" needs to be justified mathematically. How this contributes to the vanishing gradient terms seen in Figure 1 should also be better explained.\n\n**Minor comments:**\n- Line 43, \"these works\": Which works does this refer to?\n- Eq. (4): Is this missing a $- d$ on the LHS?\n- \"path\": I believe this should be \"directed path\" to be more precise.\n- Eq. (9): Is $p$ the intended exponent rather than $2$?\n- Algorithm 2, line 5: I wonder whether the $\\tilde{B}$ should be $\\tilde{B}_p^{old}$.\n- Section 4.1, Results: While the improvements seen here are large, it would be nice to mention for context that factor of $> 2$ improvements upon NOTEARS have been achieved before, e.g. with the KKT search algorithm of [35].\n- Figure 4 caption, \"our fast\": Words missing here. In general, the paper could use careful proofreading for typos, missing/extra words, etc.\n- Line 310, $k$ is some constant: As mentioned above, I don't think this is true as it depends on $\\epsilon$ and $\\tilde{B}$.", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept