On Low Rank Directed Acyclic Graphs and Causal Structure Learning
Authors: Zhuangyan Fang
Shengyu Zhu
Jiji Zhang
Yue Liu
Zhitang Chen
Yangbo He
Venue:  (2021)
Keywords: causal discovery
structure learning
low rank graphs
directed acyclic graphs

Abstract: Despite several important advances in recent years, learning causal structures represented by directed acyclic graphs (DAGs) remains a challenging task in high dimensional settings when the graphs to be learned are not sparse.  In this paper, we propose to exploit a low rank assumption regarding the (weighted) adjacency matrix of a DAG causal model to mitigate this problem. We demonstrate how to adapt existing methods for causal structure learning to take advantage of this assumption and establish several useful results relating interpretable graphical conditions to the low rank assumption. In particular, we show that the maximum rank is highly related to hubs, suggesting that scale-free networks which are frequently encountered in real applications tend to be low rank. We also provide empirical evidence for the utility of our low rank adaptations, especially on relatively large and dense graphs. Not only do they outperform existing algorithms when the low rank condition is satisfied, the performance is also  competitive even though the rank of the underlying  DAG may not be as low as is assumed.

Under review as a conference paper at ICLR 2021
ON Low RANK DIRECTED ACYCLIC GRAPHS AND CAUSAL STRUCTURE LEARNING
Anonymous authors Paper under double-blind review
ABSTRACT
Despite several important advances in recent years, learning causal structures represented by directed acyclic graphs (DAGs) remains challenging task in high dimensional settings when the graphs to be learned are not sparse. In this paper, we propose to exploit a low rank assumption regarding the (weighted) adjacency matrix of DAG causal model to mitigate this problem. We demonstrate how adapt existing methods for causal structure learning to take advantage of this assumption and establish several useful results relating interpretable graphical conditions to the low rank assumption. In particular, we show that the maximum rank is highly related to hubs, suggesting that scale-free networks which are frequently encountered in real applications tend to be low rank. We also provide empirical evidence for the utility of our low rank adaptations, especially on relatively large and dense graphs. Not only do they outperform existing algorithms when the low rank condition is satisfied, the performance is also competitive even though the rank of the underlying DAG may not be as low as is assumed.
INTRODUCTION
An important goal in many sciences is to discover the underlying causal structures in various domains, both for the purpose of explaining and understanding phenomena, and for the purpose of predicting effects of interventions (Pearl 2009). Due to the relative abundance of passively observed data as opposed to experimental data, how learn causal structures from purely observational data has been vigorously investigated (Peters et al.) 2017, Spirtes et al., 2000). In this context, causal structures are usually represented by directed acyelic graphs (DAGs) over a set of random variables.
For this task, existing methods can be roughly categorized into two classes: constraint- and score- based. The former use statistical tests extract from data a number of constraints in the form of conditional (in)dependence and seek to identify the class of causal structures compatible with those constraints (Meek) [995] Spirtes et al., 2000 Zhang 2008). The latter employ a score function to evaluate candidate causal structures relative to data and seek to locate the causal structure (or class of causal structures) with the optimal score. Due to the combinatorial nature of the acyclicity constraint (Chickering 1996; He et al. 2015). most coore-based methods rely on local heuristics perform the search. A particular example is the greedy equivalence search (GES) algorithm (Chickering 2002) that can find an optimal solution with infinite data and proper model assumptions.
Recently, Zheng et al. (2018) introduced a smooth acyclicity constraint w.r.t. graph adjacency matrix, and the task linear data models was then formulated as a continuous optimization problem with least-squares loss. This change of perspective allows using deep learning techniques to model causal mechanisms and has already given rise to several new algorithms for causal structure learning with non-linear data, e.g., Yu et al. (2019); Ng et al. (2019b a); Ke et al. 22019;; Lachapelle et al. (2020); Zheng et al. (2000), among others. While these new algorithms represent the current state of the art in many settings, their performance generally degrades when the target DAG becomes large and relatively dense, as seen from the empirical results reported in the referred works and also this paper. This issue is of course challenge to other approaches. [Ramsey et al. (2017) proposed fast GES for impressively large problems, but it works reasonably well only when the large structure is very sparse. The max-min hill-climbing (MMHC) (Tsamardinos al. 2006) relies local learning methods that often do not perform well when the target node has a large neighborhood. How to improve the performance relatively large and dense DAGs is therefore an important question.

---

Under review as a conference paper at ICLR 2021
In this work, we study the potential of exploiting kind of low rank assumption on the DAG structure to help address this problem. The rank of a graph that concerns us is the algebraic rank of its associated weighted adjacency matrix. Similar to the role of sparsity assumption on graph structures, treat the low rank assumption as methodological and it is not restricted to a particular DAG learning method. However, unlike sparsity assumption, it is much less apparent when DAGs tend to be low rank and how low rank DAGs behave. Thus, besides demonstrating the utility of exploiting a low rank assumption in causal structure learning, another important goal is to improve our understanding of the low rank assumption by relating the rank of graph to its graphical structure. Such a result also enables us to characterize the rank of graph from several structural priors and helps to choose rank related hyperparameters for the learning algorithm. Our contributions are summarized as follows:
We show how to adapt existing causal structure learning methods to take advantage of the low rank assumption, and provide strategy to select rank related hyperparameters utilizing the lower and upper bounds on the true rank, if they are available.
To improve our understanding of low rank DAGs, we establish some lower bounds on the rank of a DAG in terms of simple graphical conditions, which imply necessary conditions for DAGS to be low rank.
We also show that the maximum possible rank of weighted adjacency matrices associated with a directed graph is highly related to hubs in the graph, which suggests that scale-free free networks tend to be low rank. From this result, we derive several graphical conditions to bound the rank of DAG from above, providing simple sufficient conditions for low rank.
Empirically, we demonstrate that the low rank adaptations are indeed useful. Not only do they outperform the original algorithms when the low rank condition is satisfied, the performance is also very competitive even when the true rank is not as low as is assumed.
Related Work
The low rank assumption is frequently adopted in graph- based applications (Smith
et al. 2012: Zhou et al., 2013; Yao & Kwok] 2016 Frot et al. (201)). matrix completion and factorization (Recht (201] 1 Koltchinski et al. | |011 Cao et al. 2015 Davenport & Romberg 2016). network sciences (Hsieh et 2012 Huang et al. 2013 [ Zhang et al. 2017) and so on, but to our best knowledge, has not been used on the DAG structures in the context of learning causal DAGs. We notice two works Barik & Honorio (2019); Tichavský & Vomlel (2018) that assume low rank conditional probability tables in learning Bayesian networks, which are different from ours.
Also related are existing works that studied the rank of real weighted matrices described by a given simple unddrced/udirreeted graph. However, most works only considered the zero monzero pattern of off-diagonal entries (see, e.g.. Fallat & Hogben (2007); (Hogben (200);; Mitchell et al. (2110), whereas we also take into account the diagonal entries. This difference is crucial: if one only considers the off-diagonal entries, then the maximum rank over all possible weighted matrices is trivial and is always equal to the number of vertices. Consequently, many works focus on the minimum rank of a given graph, but to characterize exactly the minimum rank remains open, except for some special graph structures like trees (Hogben 2010). Apart from these works, [dmonds 1967) studied algebraically the maximum rank for matrices with a common zero-nonzent pattern. In Section [4] we use this result to relate the maximum possible rank to more interpretable graphical condition, which further implies several structural conditions of DAGs that may be easier to obtain in practice.
PRELIMINARIES
2.1 GRAPH TERMINOLOGG
A graph G is defined as a pair (V, E), where V = (X,, .x} is the vertex set and E C V2 denotes the edge set. We are particularly interested in directed (acyclic) graphs in the context of causal structure learning. For any S C V, we use pa(S, g). ch 9). and adj(S, 9) to denote the union of all parents, children, and adjacent vertices of the nodes of S in G, respectively. A graph is called weighted if every edge in the graph is associated with non-zero value. We will work with weighted graphs and treat unweighted graphs as a special case where the edge weights are set to 1. Weighted graphs can be treated algebraically via weighted adjaceney matrices. Specifically, the weighted adjacency matrix of a weighted graph 9 is a matrix W e &d where (i. is the weight of edge X; + X, and W (i..j) )) # ( if and only if X, - X, exists in g. The binary adjacency
2

---

Under review as a conference paper at ICLR 2021
matrix A E {0, 1] is such that A(i.j) = 1f - Xy in 9 and ((i.j) 0 otherwise. The rank of a weighted graph is defined as the rank of the associated weighted adjacency matrix.
2.2 CAUSAL STRUCTURE LEARNING AND RECENT GRADIENT-BASED METHOSS
A commonly used model in causal structure learning is the structural equation model (SEM) that describes data generating procedure. In a slight abuse of notation, we also use X; 's to denote random variables associated with the nodes in a graph G. Assuming 9 being DAG, then the SEM is given by
X; fi (pa( ec), i- 1.2.....dd
where f; is a deterministic function and e;'s are jointly independent noises. The SEM induces a marginal distribution P(X) over
[X1, X2. .... "Xal", and and P(X) are said to form a causal Bayesian network (Pearl 2009: Spirtes et al. 20000. The problem of causal structure learning
to infer the underlying causal DAG & based on the marginal distribution P(X), or more practically, an empirical version consisting of number of i.i.d. observations from P(X).
We next briefty review recently developed gradient- based methods that rely on a smooth characteriza- tion of acyclicity of directed graphs. These methods aim to find DAG that optimizes a score function and can be categorized into two classes. The first class of methods explicitly associates the target causal model with weighted adjacency matrix W and then estimate W by solving optimization problems in the following form:
min E Exx S(X, (,))) d)), subject to trace (e (ewo))= - d = 0, Wio
(1)
where h : .dd -> Rd is model function parameterized by W (and other possible parameter Φ) that aims to reconstruct X, (.) denotes score function between the true and reconstructed variables, notation denotes the element-wise product, and e.M is the matrix exponential of a square matrix M. The constraint was proposed by Zheng et al. (2018), which is smooth and holds and only W indicates a DAG. Methods in this class include: NOTEARS (Zheng et al., 2018), which targets linear models, with h(X; W. Φ) = WTX X and S(,..) being the Frobenius norm or equivalently the least-squares loss; and DAG-GNN (Yu et al., 2019) and the graph autoencoder approach (Ng) et al. 2019b), where neural networks are used for the function h with & being the weights of neural networks, and the score function be chosen as the evidence lower bound (Kingma & &Welling 2013). A sparsity inducing term may be further added when the causal graph is assumed to be sparse. These objectives are equivalent or are variants of some well studied score functions like the penalized maximum likelihood (Chickering) 2002 Van de Geer et al. 2013 Loh & Bühlmann, 2014).
The second class uses certain functions, with parameter 9, to construct weighted adjacency matrix (6)) (or a binary one A(f)) to represent the causal structure. These methods can be summarized as
min Ex~P(x) ((,,(( W (0), o)), subject to trace (e (((0))W(())) 0.6
d 0.
For example, GraN-DAG [Lacaaellee et al. 2020 and NOTEARS-MLP Zheng et al. 2020 respec- tively neural network path products and partial derivatives between variables to construct W(O).
The binary matrix A(0) can be obtained by sampling according to some distributions with learnable parameters, as used by Kalainathan et al. (2018) Ke al. 2019); Ng et al. (2019a); Zhu et al. (2020).
Before ending this section, we remark that while the gradient- based methods intend to learn a causal DAG, the learned DAG may not be identical to the underlying one for general SEMs due to the Markov equivalence (Spirtes et al. 22000: Peters et al. 2017). For such cases, one may convert the obtained DAG its corresponding Completed Partially Directed Acyclic Graph (CPDAG) as the estimate. Nevertheless, if the SEM is identifiable and proper score function is used, then the exact solution to the optimization problem is consistent, i.e., same as the true graph with probability 1; see, e.g., Shimizu et al. ) PPteers & Bühlmann (Ptee):Peters (2009). For further details and other technical issues like parameter optimization of the gradient- based methods, we refer the reader to the cited works and references therein.
|Peters et al. (2014); Zhang & Hyvärinen
3 EXPLOITING LOW RANK ASSUMPTION IN CAUSAL STRUCTURE LEARNING
This section shows how to adapt existing gradient- based methods to take advantage of the low rank assumption, by providing way for each class utilize this assumption using techniques from the
3

---

Under review as a conference paper at ICLR 2021
matrix completion literature. We remark that our adaptations with the low rank assumption not restricted to a particular learning algorithm; other DAG learning methods may potentially combine one of the proposed modifications for learning low rank causal graphs. too.
Matrix Factorization Since the weighted adjacency matrix W is explicitly optimized in the first class of methods, we can then apply the matrix factorization technique. Specifically, with estimate A for the graph rank, we can factorize W as W = UVT with U, V E R Problem I is then to optimize U and V that minimizes the score function under the DAG constraint, and has the same solution W (obtained from the product UV as the original one if f is greater than or equal to the true rank. Furthermore, if f < d. we have a much reduced number of parameters to optimize.
Nuclear Norm For the second class of methods, the adjacency matrix (0) is not an explicit parameter to be optimized. In such a case, we can adopt a commonly used technique add nuclear norm term (f)l with X > O being a tuning parameter, to the objective to induce low-rankness.
The optimization procedures these recent structure learning methods can directly incorporate the two adaptations as they are all gradient- gradient-based, though some extra care needs to be taken. Appendix C provides a detailed description of the optimization procedure and our implementation. The second approach is also feasible for the first class of methods, but we find that it does not work as well as the matrix factorization approach, possibly due the singular value decomposition to compute the Suub-)gradient w.r.t. W at each optimization step.
An acute reader may have noticed that we assumed a proper rank estimate f or a proper penalty parameter A. Yet knowing exactly the rank of the graph to be learned can be difficult in practice. Similar to the sparsity assumption, one may determine the hyperparameters 3 and À assisted by a validation dataset (or by cross- validation if the observed dataset is not sufficiently large). Alternatively, we can try different choices of the hyperparameters and then apply traditional score-based method where the search space is restricted to the resulting DAGs. However, since we are more concerned with relatively large and dense problems, the possible ranks may be too many to choose. As such, lower bound r and an upper bound ru on the graph rank would be beneficial-we need only consider ranks in [rt in the matrix factorization method, while the bounds are still useful by providing qualitative information for the nuclear norm approach: the lower an upper bound, the higher the tuning parameter À should be chosen. Moreover, lower bound can also justify the low rank assumption, i.e., if the lower bound is high, then the low rank assumption is likely to fail to hold.
4 GRAPHICAL BOUNDS ON RANKS
Obtaining exact algebraic information of a DAG such as its rank and eigenvalues may be infeasible in practice, because it may require a full knowledge of the graph to be learned. On the other hand, structural information, such as graph connectivity, distributions of in-degrees and out-degrees, and an estimate of number of hubs, is sometimes more accessible. As such, this section is devoted to relating the rank of graph to more easily interpretable graphical conditions, for the sake of a better understanding of what kinds of DAGs tend to satisfy the low rank assumption and for lower and upper bounds on the graph rank from certain structural priors.
4.1 PROBLEM SETTING
Consider a DAG G = (V, E) with weighted adjacency matrix W and binary adjacency matrix A. We aim to seek upper and lower bounds on (an using only the graphical structure. Specifically, we focus on the weighted adjacency matrices with the same binary adjacency matrix A, i.c., WA = {W E dd sign( WT) A). where sign(-) and |- are point-wise sign and absolute value functions, respectively. Notice that there exist trivial upper bound d - 1 and lower bound 0 for any DAG, but they are generally too loose for our purpose. the following, we investigate the maximum rank max feaxxrandie rank(W)), (W); ): W E }A} and minimum rank min (rank (W); W E WA} to find tighter upper and lower bounds for any W E WA. Before introducing two useful graph concepts, we comment that low rank DAGS are not necessarily sparse and vice versa; see discussion in Appendix A
Definition 1 (Height). Given DAG G (V, E) and vertex X; E V, the height of X', denoted by I(X.), is defined as the length of the longest directed path starting from X,, The height of G, denoted by IG)), is the length of the longest path in G.
4

---

Under review as a conference paper at ICLR 2021
10
10
(a) A DAG g
(b) Binary adjacency matrix
(c) A head-tail vertex cover
Figure 1: A DAG g with 12 vertices, 12 edges and height 3, where Vo X2, X1]
V1 = (Xs, Xo, X7). V2 = Xo). and V3 = [X10, X11, X12).
Definition 2 (Head-tail vertex cover). Lei 9 = (V, E) be directed graph and H, T be two subsets Vf (H, T) is called a head-tail vertex cover of G if every edge in 9 has its head vertex in H or its tail vertex in T. The size of head-tail vertex cover (H, T) is defined as [H] + (T).
As example, Figure Ic is head-tail vertex cover of G in Figure la where H = (X2, X4, Xs} (red nodes) and (Xx, Xg, X10) (blue nodes). The size of this vertex cover is 6.
4.2 LOWER BOUNDS
We first study lower bounds on the rank of weighted DAG. Define V-1 Ø and V. (X.) } for = 0, ...1,...44 ((),, Denote by 9 the induced subgraph of G over V,UV..1 Let 1() be the set of non-singleton connected components of G,,--1 and [C(G the cardi- nality. We have the following lower bounds.
Theorem 1. Let 9 be a DAG with binary adjacency matrix A. Then
Σ ((G) (C(G 2 (G)).
min(rank(WW); W E WA)
(3)
All the proof's in this paper are provided in Appendix B Theorem shows that rank (W) is greater than or equal the sum of the number of non-singleton connected components each Ga As Ga.a-1 has at least one non-singleton connected component, we obtain the second inequality. In other words, the rank of a weighted DAG is at least as high the length of the longest directed path. As an example, consider the graph shown in Figure [ One can verify that min (rank(W); W E WA} = 6, [C(91,0) 2, [C(G2
W WA) = 6 2 + +1 = 4 3. We remark that the bounds in Theorem [ may be loose in some cases. To characterize the minimum rank exactly is an on-going research problem (Hogben 2010).
1, [C(G3.2)
1, and ((G) = 3. Thus, we have min frinkran(W))
4.3 UPPER BOUNDS
We turn the more important issue for purpose, regarding upper bounds rank The next theorem shows that max max(rank (rank(W); W E WA} can be characterized exactly in graphical terms.
Theorem 2. Let be a directed graph with binary adjacency matrix A. Then max max(rank (W); W WA} is equal to the minimum size of the head-tail vertex cover of G. that is.
max(rank(W); W E WA} min[[H] + [T]; (H, T) is a head-tail vertex cover of g).
We comment that Theorem 2 holds for all directed graphs (not only DAGs), which may be of independent interest to other applications. A head-tail vertex cover of minimum size is called minimum head-tail vertex cover, which in general not unique. For a head-tail vertex cover (H, T). the vertices in H cover all the edges pointing towards these vertices while the vertices in T cover the edges pointing away. A head-tail cover relatively small size then indicates the presence of hubs, that is, vertices with relatively high in-degrees maximum rank of a weighted DAG is highly related to the presence hubs: a DAG with many hubs tends to have low rank. Intuitively, hub of high in in-degree (out-degree) is common effect (cause) of a number of direct causes (effect variables), comprising many V-structures (inverted V-structures). For example, in Figure Ta X is hub of V-structures and Xg is hub of inverted V-structures.
out-degrees. Therefore, Theorem 2 suggests that the
5

---

Under review as a conference paper at ICLR 2021
Such features are fairly common in real graph structures. Appendix A presents real network, called pathtfinee, which describes the causal relations among 109 variables Heckerman et al. 1992) with the center node being the parent of large number of other nodes. The tamous scale-free (SF) graphs also tend to have hubs. scale-free graph is one whose distribution of degree k follows a power law: P(k) k->, where r is the power parameter typically within [2,3] and P(E) denotes the fraction of nodes with degree k (Nikolova & Aluru 2012). It is observed that many real-world networks are scale-free, and some of them, such as gene regulatory networks, protein networks, and financial system network, may be viewed as causal networks (Guelzim et al., 2002, Barabasi & Oltvai 2004
Hartemink 2005 Eguíluz et al. 2005 Gao en Ren 2013 Ramsey et al., 2017). In particular, Barabasi & Oltvai (2004) claimed that most protein networks, some of which are directed and acyclic due to irreversible reactions, are the results of growth processes and preferential attachments, probably due to the gene duplication.
Empirically, the ranks of frele-free graphs are relatively low, especially in comparison to Erdös-Rényi (ER) random graphs (Mihail & Papadimitriou 2002). Figure provides simulated example where 7 is chosen from (2.3) and each reported value is over 100 random runs. As graph becomes denser, the graph rank also increases. However, for scale-free graphs with a relatively large y, the increase of their ranks much slower than that of Erdös-Rényi graphs; indeed, their ranks tend to stay fairly low even when the graph degree is large.
100
ER
+
SF-2
SF-F-3
80
60
40-
20
1
5
0
9 10
average degrea
Figure 2: 100-node graphs.
Theorem 2 can also be used to generate a low rank graph, or more precisely, a random DAG with a given rank and properly specified graph degree. Here we briefly describe the idea and leave the detailed algorithm to Appendix C.T first generate graph with r edges and rank r;; random edge is sampled without replacement and would added to the graph, if adding this edge does not increase the size of the minimum head-tail vertex cover; repeat the previous step until the pre-specified degree is reached or no edge could be added to the graph: finally, assign the edge weights randomly according to continuous distribution and the weighted graph will have rank with high probability.
The next two theorems report some looser but simpler upper bounds on rank(W).
Theorem 3. Let be a DAG with binary adjacency matrix A, and denote the set of vertices with at least one parent by Veh and those with at least one child by V pa. Then we have
((S)
min ( s/, |ch (V. 1) S Ppa.
man( W WEWWWWW WA)
199
(4)
min V 11 pa V+) < (Val,
max [ || 0 8 (()))
Since Vch and V pa are the non-root and the non-leaf vertices, respectively, the first two inequalities of (4) indicate that the maximum rank is bounded from above by the number of non-root vertices and also by the number of non-leaf vertices. The last inequality of (4) is generalization of the first two, which implies that the rank is likely to be low if most vertices have the same height.
Theorem 4. Let G be DAG with binary adjacency matrix A. Denote by skeleton (A) and moral(A)
the binary adjacency matrices of the skeleton and moral graph of G, respectively. Then we have
max(rank(W)): W E W.} I max((rnk(WW)) (W); sign([W]) = skeleton(A)}
< max)(rank(W)); sign([W]) = moral(A)).
The skeleton of DAG is the undirected graph obtained by removing all the arrowheads, and the moral graph is the undirected graph where two vertices are adjacent if they are adjacent or if they share common child in the DAG. This result is useful when the skeleton or the moral graph can be accurately estimated and the corresponding rank is low. In practice, we may use all available structural priors to obtain upper bounds on the underlying rank and choose the lowest as estimate.

---

Under review as a conference paper at ICLR 2021
5 EXPERIMENTS
This section reports empirical results of the low rank adaptations of existing methods, compared with their original versions. We choose NOTEARS (Zheng et al., 2018) for linear SEMs by adopting the matrix factorization approach, denoted as NOTEARS-low-rank, and use the nuclear norm approach in combination with GraN-DAG (Lachapelle et al., 2020) for a non-linear data model. Again we remark that the two methods are only demonstrations of the utility of low rank assumption, which can be potentially combined with other methods as well. For more information, we also include several benchmark methods: fast GES (Ramsey et al. 2017), PC (Spirtes et al. (2000), MMHC (Tsamardinos)
et al.) 2006), ICA-LiNGAM (Shimizu et al. 2006) specifically designed with non-Gaussian noises, for linear SEMs T and DAG-GNN (Yu et al., 2019). NOTEARS-ME MMLP (Zheng et al. 2020), and CAM (Bühlmann et al. 2014) for the non- non-linear case. Their implementations are described in Appendix C
We consider randomly sampled DAGs with specified ranks (the generating procedure was described in Section 4.3 and is given as Algorithm in Appendix [C.T). scale-free graphs, and a real network structure. For linear SEMs, the weights uniformly sampled from [-2, - 0.5] U [0.5, 2] and the noises are either standard Gaussian or standard exponential. For non- non-linear SEMs, we use additive Gaussian noise model with functions sampled from Gaussian processes with RBF kernel of bandwidth one. These data models are known to be identifiable (Shimizu et al.) 2006] Peters & Bühlmann] 2013) Peters et al. 2014). From each SEM, we then generate n = 3, 000 observations. We repeat ten times over different seeds for each experiment setting. Detailed information about the setup can be found in Appendix C.3 Below we mainly report structural Hamming distance (SHD) which takes into account both false positives and false negatives, and a smaller SHD indicates a better estimate.
5.1 LINEAR SEMS WITH
GRAPHS
We first consider linear SEMs on rank-specified graphs, with number of nodes d E {100, 300). rank r [0.1d], and average degree & E (2, 4, 6, 8). The true rank is assumed to be known and is used as the rank parameter in
For a better visualization, Figure 3 only reports the average SHDs, while the true positive rate, false discovery rate, and running time left to Appendix D We also show the results after using the interquartile range rule to remove outlier SHDs. We observe that the low rank assumption can greatly improve the performance of NOTEARS, reducing the SHDs by at least a half. For this data model, the fast GES has much higher SHDs (see also Appendix D). PC is too slow (for example, did not finish in 16 hours for a dataset with 100 nodes and degree 6). because some nodes may have a high in-degree. For the same reason, the skeleton may not be well estimated by MMHC; its performance is slightly worse than the fast GES and is not reported.
-dh
NOTEARS-LLI
NOTEARS -
-
NOTEARS- low--aak-nnkk
-X- ICA-LINAAM
-A-
NOTEARS (WO outliors)
NOTEARS (w/b outlers)
-9- NOTEARS loW rank (mi outliers)
-X-
LINAA (w/o outliers)
120
240
100
200
90
360
SHD
60
120
SHD
60
40
30
20
40
0
everage degree
average degree
everage degree
(a) 100, Gaussian
(b) d 300, Gaussian
(c) d = 100, Exponential
Figure 3: Average SHDs on rank-specified graphs. The models are linear SEMs with (a)-(b) Gaussian noises, and (c) exponential noises. The true rank is assumed to be known.
For more information regarding the role of sparsity, we include NOTEARS with an & penalty, named NOTEARS-LI. Here the & penalty weight chosen from (0.01, 0.02, 0.05, 0.1, 0.2, 0.53. Instead
Here we choose ICA-LMNGAM, other than alternative LINGAM methods like DirectLiNGAM (Shimizu et al., 2011). based on our empirical observation. Specifically, an implementation of ICA LiNGAM has a noticeably hetter performance than Directl.NGGA for relatively dense graphs. Please find detailed discussion and empirical comparison in Appendix D..
7

---

Under review as a conference paper at ICLR 2021
of relying on an additional validation dataset, we treat NOTEARS-LI favorably by picking the lowest SHD obtained from different weights for each dataset. As seen from Figure 3a NOTEARS-LA L1 is slightly better than NOTEARS when the average degree is 2. but is largely outperformed with relatively dense graphs. This observation was also reported in ZZeng et al. (2018). We conjecture that it is because our experiments consider relatively sufficient data and dense graphs. Moreover, the thresholding procedure controls false discoveries and may have similar effect to the L penalty,
Appendix D.I studies graphs with higher ranks, where it is observed that the advantage of NOTEARS- low-rank over NOTEARS decreases when the rank of the underlying DAG increases. Nevertheless, NOTEARS-IOw-rank is still competitive when the true rank is [d/2] and the factorized matrix has the same number of parameters as NOTEARS. We also conduct an empirical analysis with different sample sizes in Appendix [D.2 which shows that NOTEARS- l-w--kkk performs reasonably well when the sample size is small and tends to have better performance with larger number of samples.
Due to space limit, please find further details in the appendix.
5.2 LINEAR SEMs WITH SCALE-FREE GRAPHS
We next consider scale-free graphs with d = 100 nodes, average degree
6, and power 7 2.5. For this experiment, the minimum, maximum, and mcan ranks of generated graphs are 14, 24, and 18.7. respectively. Here we choose the rank parameter i from {20, 30, 40) for NOTEARS- PORSS-lo-rran low-ran.. As seen from Figure NOTEARS- low-rank with rank parameter f 20 performs the best, even though there are graphs with ranks greater than 20.
150
mean
120-
SHD
90-
60-
30
NOTEARS
20
20
40
rank
Figure 4: Scale-free graphs.
5.3 SENSITIVITY OF RANK PARAMETERS AND VALIDATION
So far we have assumed that the true rank or an accurate estimate is known. In this experiment, we conduct an empirical analysis with different rank parameters for linear Gaussian data model on rank-specified graphs with 100 nodes, degree and rank 10. We also include the validation based approach where 2, 000 samples are chosen as training dataset and the rest as validation dataset. We use the derived lower and upper bounds in Theorems T and 3 to obtain a range of possible rank parameters, assuming that the corresponding structural priors are available. Within this range, we then select evenly distributed rank parameters used with NOTEARS--oww-aan Finally, evaluate each learned DAG using the validation dataset and choose the DAG with the best score as our estimate.
to learn causal graphs.
As seen from Figure 5 NOTEARS-low-rand performs the best when the rank parameter is identical to the true rank, while the rank parameter chosen by validation has almost the same performance. Compared with NOTEARS on the same datasets, the low rank version performs well across range of rank parameters. Although this validation approach increases the total running time that depends on the number of candidate rank parameters, we believe that it is acceptable given the gained accuracy and also the fact that this strategy has been frequently adopted for tuning hyperparameters in practice.
350
C
3 150
250
120
90
+ HOTEARS-MIL
CAM
GraN-DAG, rank.
GraN-DAG
100
60-
GraN low rank.
50
GraN DAG, low rank. 1.0
GraN DAG, lowira Wk, val dation
NOTEARS 5
validution
average degree
Figure Different rank parameters.
Figure 6: Non-linear SEMs.
5.4 NON-LINEAR SEMS
For non-lin non-linear data models, we pick rank-specified graphs with 50 nodes, rank 5, and average degree k E (2, 4, 6, 8). To our knowledge, the selected benchmark methods CAM, NOTEARS-MIJ. and
8

---

Under review as a conference paper at ICLR 2021
GraN-DAG are state-of-the-ard methods on this data model. As demonstration of the low rank assumption. we apply the nuclear norm approach to GraN-DAG and choose from {0.3, 0.5, 1.0) as penalty weights. For validation,
use the same splitting ratio in Section 5.3 and consider more
penalty weights from (0.1, 0.2, 0.3, 0.5, 1, 2, 5}. Similarly, the learned graph that achieves the best score the validation dataset is chosen final estimate. Figure (and Appendix with more detailed result) shows that adding nuclear norm can improve the performance of GraN-DAG across a large range of weights when the graph is relatively dense. For degree the low rank version with validation achieves average SHD 77.4, while the SHDs of CAM, NOTEARS-MIA and original GraN-DAG are 131.9, 119.4, and 109.4, respectively.
5.5
REAL NETWORK
We apply the proposed method to the 5rth150 gene net- work, which is a DAG containing 107 genes and 150 edges. Its maximum rank is 40. Since the real dataset has only 22 samples, we instead use simulated data from linear Gaussian SEMs. We pick f from {36, 40, 44} and also use validation to select the rank parameter. We apply NOTEARS-L.I where the & penalty weight is chosen from (0.05, 0.1,0.2), 0.2). and similarly treat this method favorably by picking the lowest SHD for each dataset. The mean and median SHDs are shown in Figure 7 Using Student's 1-test, we find that with significance level 0.1, the results obtained with f = 44 and the validation approach are significantly better than NOTEARS. This experiment demonstrates again the utility of the low rank assumption, even when the true rank the graph is not very low.
40
10
NOTEARS L1
NOTEARS
tw33
ind4
validation
Figure 7: Real network.
6 CONCLUDING REMARKS
This paper studies the potential of low rank assumption in causal structure learning. Empirically, we show that the low rank adaptations perform noticeably better than existing algorithms when the low rank condition is satisfied, and also deliver competitive performances when the rank is not as low as assumed. Theoretically, we provide an improved understanding of what kinds of graphs tend to be low rank and possibility to obtain bounds on the underlying rank from several structural priors.
We treat the present work as our first step to incorporate low-rankness into causal DAG learning. future direction is to approximate a high rank DAG with a low rank one (possibly adding additional DAG that is sparse). While there a rich literature on low rank approximations of matrices and combining low-rankness with sparsity, it is non-trivial to to conclude under what conditions such an approximation is guaranteed to be effective to learn causal DAGs, Another direction is to compare the low rank assumption to other structural or parametric priors affecting model selection through marginal likelihood (Eggeling et al. 2019 Silander et al. 2007). Finally, it is also interesting to investigate if low rank DAG model implies any useful behavior in the data.
REFERENCES
Albert-Laszlo Barabasi and Zoltan N Oltvai. Network biology: Understanding the cell's functional organization. Nature Reviews Genetics, 5(2):101-113, 2004.
Adarsh Barik and Jean Honorio. Learning Bayesian networks with low rank conditional probability tables. In Advances in Neural Information Processing Systems (NeurIPS), 2019.
Dimitri P Bertsekas. Nonlinear Programming. Athena Scientific, 1999.
Peter Bühlmann, Jonas Peters, and Jan Ernest. CAM: Causal additive models, high-dimensional order search and penalized regression. The Annals of Statistics, 42(6):2526-2566, 2014.
Xiangyong Cao, Yang Chen, Qian Zhao, Deyu Meng, Yao Wang. Dong Wang, and Zongben Xu. Low-rank matrix factorization under general mixture noise distributions, In The IEEE International Conference on Computer Vision (ICCV), 2015.

---

Under review as a conference paper at ICLR 2021
David Maxwell Chickering. Learning Bayesian networks is NP-complete. In Learning from Data: Artificial Intelligence and Statistics V. Springer, 1996.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3(Nov):507-55) 2002.
Mark A. Davenport and Justin Romberg. An overview of low-rank matrix recovery from incomplete observations. IEEE Journal of Selected Topics in Signal Processing, 10(4):608-22,, 2016.
König Dénes. Gráfok és mátrixok. Matematikai és Fizikai Lapok, 38:116-119, 1931.
Jack Edmonds. Systems of distinct representatives and linear algebra. Journal of Research of the Notional Bureau of Standards B. Mathematics and Mathematical Physics, 71B(4):241-245,
Oct-Dec 1967.
Ralf Eggeling, Jussi Viinikka, Aleksis Vuoksenmaa, and Mikko Koivisto. On structure priors for learning Bayesian networks. In International Conference on Artificial Intelligence and Statistics (AISTATS), volume 89, pp. 1687-1695, 2019.
Victor M. Eguíluz, Dante R. Chialvo, Guillermo A. Cecchi, Marwan Baliki, and A. Vania Apkarian. Scale-free brain functional networks. Physical Review Letters, 94:018102, Jan 2005.
Shaun M Fallat and Leslie Hogben. The minimum rank of symmetric matrices described by a graph: survey, Linear Algebra and its Applications, 426(2-3):558-58) 2007.
Benjamin Frot, Preetam Nandy, and Marloes H. Maathuis. Robust causal structure learning with some hidden variables. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 8133):5594877 2019.
Bo Gao and Ruo en Ren. The topology of a causal network for the Chinese financial system. Physica A: Statistical Mechanics and its Applications, 392(13):2965 -976, 2013.
Nabil Guelzim, Samuele Bottani, Paul Bourgine, and François Képès. Topological and causal structure of the yeast transcriptional regulttrry network. Nature Genetics, 31:60-63, 2002.
Alexander J Hartemink. Reverse engineering gene regulatory networks. Nature Biotechnology, 23(5): 554-555, 2005.
Yangbo He, Jinzhu Jia, and Bin Yu. Counting and exploring sizes of Markov equivalence classes of directed acyelic graphs. Journal of Machine Learning Research, 16:2589-2609, 2015.
David Earl Heckerman, Eric Horvitz, and Bharat N Nathwani. Toward normative expert systems: Part i the pathfinder project. Methods of Information Medicine, (100))90- 1992.
Leslie Hogben. Minimum rank problems. Linear Algebra and Applications, (32(8):1961-3) 2010.
Cho-Jui Hsich, Kai- Yang Chiang, and Inderjit S. Dhillon. Low rank modeling of signed networks. In the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2012.
Jin Huang, Feiping Nie, Heng Huang, Yu Lei, and Chris Ding. Social trust prediction using rank-k matrix recovery. In International Joint Conference on Artificial Intelligence (IJCAI), 2013.
Diviyan Kalainathan, Olivier Goudet, Isabelle Guyon, David Lopez-Paz, and Michèle Sebag. Struc- tural agnostic modeling: Adversarial learning of causal graphs. arXiv preprint arXiv: 1803.049299 2018.
Nan Rosemary Ke, Olexa Bilaniuk, Anirudh Goyal, Stefan Bauer, Hugo Larochelle, Bernhard Schölkopf, Michael C. Mozer, Chris Pal, and Yoshua Bengio. Learning neural causal models from unknown interventions. arXiv preprint 191.00055 2019.
Diederik P Kingma and Max Welling. Auto-encoding variational bayes. In International Conference Learning Representations (ICLR), 2013.
10

---

Under review as a conference paper at ICLR 2021
Vladimir Koltchinskii, Karim Lounici, and Alexandre B. Tsybakov. Nuclear- penalization and optimal rates for noisy low-rank matrix completion. The Annals of Statistics, 39(5):3302-233) 10 2011.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient-based neural DAG learning. In International Conference on Learning Representations (ICLR), 2020.
Po-Ling Loh and Peter Bühlmann. High- High-dimensional learning of linear causal networks via inverse covariance estimation. The Journal of Machine Learning Research, 15(1):3065-3105, 2014.
Christopher Meek. Causal inference and causal explanation with background knowledge. In UAI, 1995.
Milena Mihail and Christos Papadimitriou. On the eigenvalue power law. In Randomization and Approximation Techniques in Computer Science, pp. 254-262, Berlin, Heidelberg, 2002. Springer.
Lon H. Mitchell, Sivaram K. Narayan, and Andrew M. Zimmer. Lower bounds in minimum rank problems. Linear Algebra and its Applications, 432(()::400 440, 2010.
Ignavier Ng, Zhuangyan Fang, Shengyu Zhu, Zhitang Chen, and Jun Wang. Masked gradient- causal structure learning. arXiv preprint arXiv: 1910.0857, 2019a.
Ignavier Ng, Shengyu Zhu, Zhitang Chen, and Zhuangyan Fang. graph autoencoder approach to causal structure learning. arXiv preprint arXiv: 1et.2.9111.. 2019b.
Olga Nikolova and Srinivas Aluru. Parallel Bayesian network structure learning with application to gene networks. Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, 2012.
Judea Pearl. Causality. Cambridge University Press, 2009.
Jonas Peters and Peter Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101(1):219-228, 2013.
Jonas Peters, Joris M Mooij. Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. The Journal of Machine Learning Research, (5()):2009-2053 2014.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of Causal Inference Foundations and Learning Algorithms. MIT Press, Cambridge, MA, USA, 2017.
Joseph Ramsey, Madelyn Glymour, Ruben Sanchez-Romero, and Clark Glymour. A million variables and more: the fast greedy equivalence search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images. International Journal of Data Science and Analytics, 3(2):121-129, 2017.
Benjamin Recht. A simpler approach to matrix completion. Journal of Machine Learning Research, 12(104):343-3433 2011.
Shohei Shimizu, Patrik Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(0ct):2003-2030 2006.
Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvärinen, Yoshinobu Kawahara, Takashi Washio, Patrik O Hoyer, and Kenneth Bollen. Directlingam: A direct method for learning linear non- Gaussaa structural equation model. Journal of Machine Learning Research, 12(Apr): 1225-1248, 2011.
Tomi Silander, Petri Kontkanen, and Petri Myllymaki. On sensitivity of the map Bayesian network structure to the equivalent sample size parameter. In Uncertainty in Artificial Intelligence (UAI), 2007.
Carl Smith, Frank Wood, and Liam Paninski. Low rank continuous-spi space graphical models. In International Conference Artificial Intelligence and Statisties (AISTATS), 2012.
11

---

Under review as a conference paper at ICLR 2021
Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, Prediction, and Search. MIT Press, second edition, 2000.
Petr Tichavský and Jifí Vomlel. Representations of bayesian networks by low-rank models. In International Conference on Probabilistic Graphical Models, 2018.
loannis Tsamardinos, Laura E. Brown, and Constantin F Aliferis. The max-min hill-climbing Bayesian network structure learning algorithm. Machine learning, 65(1):31-78, 2006.
Sara Van de Geer, Peter Bühlmann, al. lo-penalized maximum likelihood for sparse directed acyclic graphs. The Annals of Statistics, 41(2):536-677, 2013.
Quanming Yao and James T. Kwok. Greedy learning of generalized low-rank models. In International Joint Conference Artificial Intelligence, 2016.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In International Conference on Machine Learning (ICML), 2019.
Jiawei Zhang, Jianhui Chen, Shi Zhi, Yi Chang, Philip S. Yu, and Jiawei Han. Link prediction across aligned networks with sparse and low rank matrix estimation. IEEE 33rd International Conference Data Engineering. 2017.
Jiji Zhang. On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias. Artificial Intelligence, 172(16):1873 1896, 2008.
Kun Zhang and Aapo Hyvärinen. On the identifiability of the post-nonlinear causal model. In Uncertainty in Artificial Intelligence (UAI), 2009.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. DAGS with NO TEARS: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems (NeurIPS), 2018.
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric P Xing. Learning sparse non- parametric DAGs. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2020.
Ke Zhou, Hongyuan Zha, and Le Song. Learning social infectivity in sparse low-rank networks using multi-dimennional Hawkes processes. In International Conference on Artificial Intelligence and Statistics (AISTATSS, 2013.
Shengyu Zhu, Ignavier Ng. and Zhitang Chen. Causal discovery with reinforcement learning. In International Conference on Learning Representations (ICLR), 2020.
12

---

Under review as a conference paper at ICLR 2021
Appendix
A EXAMPLES AND DISCUSSIONS
We provide more examples and discussions in this section.
Minimum rank of the graph in Figure [ We first show that the minimum rank of the DAG structure in Figure 6. It is clear that the 6-th to 10-th rows of A are always linearly independent, so it suffices to show that the 11-th row linearly independent of the 6-th to 10-th rows. To see this, notice that if the 11-th row is a linear combination of the 6-th to 10-th rows, then A(11, 1) would be non-zero, which is a contradiction.
The pathfinder and arth150 networks Figure visualizes the pathfinder and arth150 networks that are mentioned in Sections 4.3 and 5 respectively. Both networks can be found at http: As one can see, these two networks contain hubs: the center note in the pathyfiner network has a large number of children, while the arth network contains many 'small' hubs, each which has ~~ children. We also notice that nearly all the hubs in the two networks have high out-degrees.
//www. www.aan.e....................... www.///rrepsssttoooooooooeoeereo
Figure 8: The pathfinder (left) and arth150 (right) networks.
Sparse DAGs and low rank DAGs A sparse DAG does not necessarily indicate low rank DAG, and vice versa. For example, a directed linear graph with d vertices has only d edges, i.c. Xa. while the rank of its binary adjacency matrix d 1. According to Theorems and 2 the maximum and minimum ranks of a directed linear graph are equal to its number of edges. Thus, directed linear graphs are sparse but have high ranks. On the other hand, for some non-sparse graphs,
X1 X2
assign the edge weights so that the resulting graphs have low ranks. simple example would be a fully connected directed balanced bipartite graph, as shown Figure D The definition of bipartite graphs its two parts contain the same number vertices. The rank of fully connected balanced bipartite graph with d vertices is 1 if all the edge weights the same (e.g., the binary adjacency matrix), but the number of edges d²
be found in Appendix B.1 A bipartite graph is called balanced
6
8
10
3
1
(a)
(b) Binary adjacency matrix
Figure 9: A fully connected directed balanced bipartite graph G and its binary adjacency matrix.
We also notice that there exist some connections between the maximum rank and the graph degree, or more precisely, the total number of edges in the graph, according to Theorem 2 Intuitively, if the
13

---

Under review as a conference paper at ICLR 2021
graph is dense, then need more vertices to cover all the edges. Thus, the size of the minimum tail vertex cover should be large. Explicitly providing a formula to characterize these two graph parameters is an interesting problem, which will be explored in the future.
B PROOFS
In this section, we present proofs for the theorems given in the main content.
B.1 PRELIMINARIES
A bipartite graph is a graph whose vertex set V can be partitioned into two disjoint subsets Vo and VI, such that the vertices within each subset are not adjacent one another. Vo and V are called the parts of the graph. A matching of a graph is a subset of its edges where no two of them share a common endpoint.
vertex cover of a graph is subset of the vertex set where every edge in the graph has at least one endpoint in the subset. The size of matching (vertex cover) is the number of edges (vertices) in the matching (vertex cover). A maximum matching of a graph a matching of the largest possible size and a minimum vertex cover is vertex cover of the smallest possible size. An important result about bipartite graphs is König's 's theorem (Dénes 1931), which states that the size of minimum vertex cover is equal to the size of maximum matching in bipartite graph.
Based on the heights of vertices in V, we can define a weak ordering among the vertices: X X; if and only if (XX.) I(X,), and X; X, if and only if (XX,) = l(X,). Given this weak ordering, we can group the vertices by their heights, and the resulting graph shows a hierarchical structure; see Figure 1 in the main text for an example. This hierarchical representation has some simple and nice properties. Let V, = [X2; (X,, = 8), 8 = 0, 1......, ((9), and let V 1 = Ø. We have: (1) for any given E {0, 1...., (({)) and two distinct vertices X2 V X1 and X2 are not adjacent, and (2) for any given s E (1,,)...,((( [[G)] and X; E V.. there is at least one vertex in V which child of X,, If we denote the induced subgraph of G over V, U V.-1 by 92.8-1, then G is a bipartite graph with V, and V,-1 as parts, and singletons in - (i.e., vertices that are not endpoints of any edge) only appear in V.-1-
For ease of presentation, we occasionally use index i to represent variable X, in the following sections.
B.2 PROOF OF THEOREM I
Proof. Let 9 = EV, Consider an equivalence relation, denoted by re, among vertices in V defined as follows: for any X,, X2 E V. X, if and only if (X,)) I(X,) and X; and X, are connected. Here, connected means that there is path between X; and X,. Below we use C(X,) to denote the equivalence class containing X,. Next, we define a weak ordering equivalence classes induced by ~x, by letting C(X.) E C(X,) if and only if ((X)) 2 ((x)). Then, we extend >,, total ordering p on V/ ~ The ordering p also induces weak ordering (denoted by p) on V: X; 2 X, if and only if ((2,) >p C(X,). Finally, we extend 5 to total ordering on V. It can be verified that y a topological ordering of G. that is, if we relabel the vertices according to y, then X; E pa(X,. G) if and only if i > and X; and X, are adjacent, and the adjacency matrix of G becomes lower triangular.
V/ ~ i.e., the
Assume that the vertices of G are relabeled according to y and we will consider the binary adjacency matrix A of the resulting graph throughout the rest of this proof. Note that relabelling is equivalent to applying a permutation onto the adjacency matrix, which does not change the rank. Let Vo- (1,2...., An - 1) for some hij 2 2. Then the ki-th row of A, denoted by ),, is the first non-zero row vector of A. Letting S (A(ki,)) then S contains subset of linearly independent vector(s) of the first bj rows of A. Suppose that we have visited the first m rows of A and = (A(k1,) :), A(kz, ,)... A(ks, -}} contains subset of linearly independent vector(s) of the first m rows of A, where b S m < d. If Xm+1 ~e Xt., then we add A(m++mm)))))) +. to S: otherwise, we keep S unchanged. We claim that the vectors in S are still linearly independent after the above step. Clearly, if we do not add any new vector, then S contains only linearly independent vectors. To show the other case, note that if ((Xm+1)
I(Xx,) 2 2 I(Xx,), then there an index i E Vixx41)-11 such that A(m + 1, # O, by the definition of height. Since 1(Xm+t) > 1(Xx,), we have
14

---

Under review as a conference paper at ICLR 2021
((xx,) S ((Xm+1) and thus A(k,, = 0 for all 1.2....,,t Therefore, A(m + 1, cannot be linearly represented by {A(kk ) 2....... .t..
and the vectors are linearly independent. On the other hand, if ((xm+1) I( then the definition of the equivalence relation implies
that x and XR, are disconnected, which means that X and Xk, do not share a common
child in V ((X ((Xm+))-1 Consequently, there is an index i V (((mm4)) ) such that A(m + 1, / O show that A(kj, i) for all 1,,,...,t Thus, the vectors
but Ak,, i) = O. Similarly, in S are still linearly independent.
After visiting all the rows in A, the number of vectors in S is equal to 1((9) [C(Gs,a-1) based on the definition of ~. The second inequality can be shown by noting that (( has at least one elements. The proof is complete.
B.3 PROOF OF THEOREM 2
Proof. Denote the directed graph by G (V,E). Edmonds (1967 Theorem 1) showed that max (rank(W) W E Wa} is equal to the maximum number of nonzero entries of A, no two of which lie in a common row column. Therefore, it suffices to show that the latter quantity is equal to the size of the minimum head-tail vertex cover. Let V V UV, where V = x {0} {(X,, 0); X; E VJ and V = V x (1) = {(X,, 1); X, E VJ. Now define a bipartite graph B = (V', E') where E' = 0( - (X (X;,))(X, Xj) E E). Denote by M set of nonzero entries of so that no two entries lie in the same row or column. Notice that M can be viewed as an edge set and two edges in .M share common endpoint. Thus, M is matching of B. Conversely, it can be shown by similar arguments that any matching of corresponds to set nonzero entries of A, no two of which lie in common row column. Therefore, max (rank (W), W E Wa) equals the size of the maximum matching of B, and further the size of the minimum vertex cover of B according König's theorem. Note that any vertex cover of B be equivalently transformed to a head-tail vertex cover of G, by letting H and T be the subsets of the vertex cover containing all variables in Vo and of the vertex cover containing all variables in Vi, respectively. Thus. (raxkrank(W) W E Wa} is equal to the size of the minimum head-tail vertex cover.
B.4 PROOF OF THEOREM 3
Proof. We start with the first inequality in Equation 44p. Let h hi,.... hp denote the heights where lch(V.). Let H V and T U .,. It is straightforward to see that (H, T) is a head-tail vertex cover. Thus, Equation (4) holds according to Theorem _ The second inequality can be shown similarly and its proof is omitted. For the third inequality, let
IV ch(V.)),, and ti....,ty .t the height where
>
argmax{[V.l. 0 I ((g)), and define H Ulom V, and = Uicm Vz- Then (H, T) is also a head-tail vertex cover and the third inequality follows from
Theorem
B.5 PROOF OF THEOREM
Proof. Notice that Theorem 2 holds for all directed graphs. This theorem then follows by treating the skeleton and the moral graph as directed graphs with loops, i.e., an undirected edge X, X, is treated as two directed edges X; - X, and X, - X,.
C IMPLEMENTATION DETAILS
In this section, we present an algorithm to generate random DAG with given rank, a low rank version of NOTEARS and GraN-DAG, and also a description of our experimental settings.
C.I GENERATING RANDOM DAGS
In Section 4.3
we briefly discuss the idea of generating a random DAG with a given rank. We now describe the detailed procedure in Algorithm In particular, we aim generate random DAG with nodes, average degree k and rank r. The first part of Algorithm after initialization is to sample a number N, representing the total number of edges, from binomial distribution B(d(d 1)
15

---

Under review as a conference paper at ICLR 2021
Algorithm I Generating random DAGS
Require: Number of nodes d, average degree k. and rank r.
Ensure: A randomly sampled DAG with the number nodes d, average degree and rank r.
1: Set M empty graph. Mp Ø and R = {(i,]): i < j, i.j 1.2....,d)
2: Set = k//dd 1).
3: Sample number N S(d(d-1)/2/) p), where B(n, p) is binomial distribution with parameters
n and p.
4: if N < then
5: return FAIL
6: end if
7: Sample r indices from 1....,d d - and store them in Mp in descending order.
8: for each i Mp do
9: Sample an index j from + 1 to d.
Add edge (i to M and remove (i,j) from R.
11: end for
12: while * Ø and [M| < N do
13: Sample an edge (i.j) from R and remove it from R.
14: if adding (i,j) to M does not change the size of the minimum head-tail vertex cover of M
then
15:
Add (i,j) to M.
16: end if
17:
end while
if [M < N then
19: return FAIL
20: end if
21: return M
where p- k/(dd/dd)) If < Algorithm would return FAIL since graph with N < edges could never have rank r. Otherwise, Algorithm samples an initial graph with r edges and rank r, by choosing edges such that no two of them share the same head points the same tail points, i... each row and each column of the corresponding adjacency matrix have at most one non-zero entry. Then, Algorithm [ sequentially samples edge from R containing all possible edges and checks whether adding this edge to the graph changes the size of the minimum head-tail vertex cover. If not, the edge will be added to the graph; otherwise, it will be removed from R. This because if graph G is a super-graph of another graph H. then the size of the minimum head-tail cover of G is no less than that of H. We repeat the above sampling procedure until there is edge in R the number of edges in the resulting graph reaches N. If the latter happens, the algorithm will return the generated graph; otherwise, will return FAIL.
The theoretic basis of Algorithm T is Theorem Note that the algorithm may not return a valid graph if the desired number N of edges cannot reached. This could happen if the input rank is too low while the input average degree is too high. With our experiment settings, we find it rare for Algorithm to fail return a desired graph.
C.2 OPTIMIZATION
For this part, we consider a dataset consisting of n i.i.d. observations from P(X) and consequently the expectations in Problems (I) and (2) are replaced by empirical means. Denote the design matrix by X E xd where each row of X corresponds to an observation and each column represents variable. Here we use NOTEARS (Zheng et al. 2018) and Gran- DAG (Lachapelle et al. 2020) from each class of methods examples and will describe their low rank versions in the following. Other gradient-Be based methods and their optimization procedures can be similarly modified to incorporate the low rank assumption.
16

---

Under review as a conference paper at ICLR 2021
Algorithm 2 Optimization procedure for NOTEARS-low-tankk
Require: Design matrix X, starting point (Uo, Vo. au), rate E (0, 1), tolerance 0, and threshold
w > 0.
Ensure: Locally optimal parameter W*
1: for 1 =1 1, 1.2,...doo do
2:
(Solve primal) U241, Ve+1 t arg mingyv Lp, aa) with p such that g(Ub+1 Vi+1) < cg(U,)
3:
(Dual ascent) a++1 + ar po(Ue+1 V2-)).
4:
if g(Uu++1V24 < then
5:
Set U U441 and V V1+1-
6: break
7: end if
8: end for
9:
(Thresholding) Set V = UVV 1[1*V= > w).
10:
return W
C.2.1 NOTEARS WITH LOW RANK ASSUMPTION
Following Section the optimization problem in our work be written
|X XUVTIE subject to trace (c (eVVTOUVT) - -d= 2n
min
(5)
where U, V E Rdx and o is the point-wise product. The constraint in Problem (5) holds if and only if UVT is weighted adjacency matrix of a DAG. This problem then be solved by standard numeric optimization methods such as the augmented Lagrangian method (Bertsekas 1999). In
particular, the augmented Lagrangian is given by
1 | XUV og((VV) + Glo(UVT)) 2n
L,,U,
where (UVV)) = trace (SUVTOUVT d, is the Lagrange multiplier, and p is the penalty parameter. The optimization procedure is summarized in Algorithm 2 similar to Zheng et al. (2018 Algorithm 1). Notice that here we do not include the ε penalty term (except for the first and last experiments in Sections 5.1 and 5.5 respectively), for the following reasons: (1) the thresholding procedure can also control false discoveries; (2) we consider relatively sufficient data for the experiments and NOTEARS with thresholding has been shown in /Zhnng et al. (2018) to perform consistently well even when the graph is sparse; (3) we are more concerned with relatively large and dense graphs, so sparsity assumption may be harmful, as shown also by Zheng et al. (2018); (4) the & penalty term requires a tuning parameter, which itself is not easy to choose.
Zheng al. 2018) used L-BFGS solve the unconstrained subproblem in Step 2. We alternatively use the Newton conjugate gradient method that is written in C. Empirically, these two optimizers much faster thanks
behave similarly in terms of the estimate performance, while the latter to its C implementation. The DAG constraint may not be satisfied exactly using iterative numeric methods, is a common practice pick small tolerance, followed by thresholding procedure on the estimated entries to obtain exact DAGs. In our implementation, we choose Uo and Vo to be the first f columns of the d x d identity matrices. Other parameter choices are: 0g = 0, c 0.25, e 10 and w 0.3, similar to those used in related methods on the same datasets (e.g., Zheng al. (2018): Yu al. (2019); Zhu et al. (2220). The chosen threshold our experiments and in the experiments of related works that use the same data model. In case the thresholded matrix is not a DAG, one may further increase the threshold until the resulting matrix corresponds to DAG.
0.3 works well in
After obtaining W*, we add an additional pruning step: we use linear regression to refit the dataset based on the structure indicated by W and then apply another thresholding (with 0.3) to the refitted weighted adjacency matrix. Both the Newton conjugate gradient optimizer and the pruning technique also applied to NOTEARS, which not only accelerate the optimization but also improve its performance by obtaining a much lower SHD, particularly for large and dense graphs. See Appendix D33 for empirical comparison.
17

---

Under review as a conference paper at ICLR 2021
C.2.2 GRAN-DAG WITH LOW RANK ASSUMPTION
We next consider a low rank version of GraN-DAG, The optimization problem can be written
min - n
Σ
logp ( pa(X, ((0))()) (0)l
1=1
(6)
subject to trace (c 0,
where x! is the 1-th sample of variable X, and pa(X,, (0)) means the 1-th sample of X,'s parents indicated by the adjacency matrix W()). Here, denotes the parameters of neural networks and (()) with non-negative entries is obtained from the neural network path products.
Problem can be solved similarly using augmented Lagrangian. The procedure is similar to Algorithm 2 and is the same to that used by GraN-DAG, with slight modifications: (1) the subproblem in Step 2 approximately solved using first-order methods; (2) the thresholding at Step 9 is replaced by variable selection method proposed by [Bhhlmann et al. (2014). The same variable selection or pruning method adopted by two other benchmark methods CAM and NOTEARS-MI in our experiment. Please refer to Lachapelle et al. (2020) and Bühhlman et al. (0144) for further details.
C.3 EXPERIMENT SETUP
In our experiments, we consider three data models: linear Gaussian SEMs, linear non Ganssaasi SEMs (linear exponential SEMs), and non-linear SEMs (Gaussian processes). Given randomly generated DAG G, the associated SEM is generated as follows:
Linear Gaussian
A linear Gaussian SEM is given by
X; =
Σ W +4, i= 11 ()E))(X..9)
where va(X, denotes X,'s parents in G and ej's 's are jointly independent standard Gaussian noises. In our experiments, the weights Wii,j)'s uniformly sampled from [-2, -0.5] [0.5, 2].
Linear Exponential A linear exponential SEM is also generated according to Equation 7 where &; ss are replaced by jointly independent Exp(1) random variables. The weights (i,jj)'s are sampled from [-2, -05] U [0.5, uniformly, too.
Gaussian Processes
We consider the following additive noise model:
X; = ()) + i= 1.2.......
where éss are jointly independent standard Gaussian noises and fi's are functions sampled from Gaussian processes with RBF kernel of bandwidth one.
We sample 3. 000 observations according the SEM. The reported results of each setting summa- rized over 10 repetitions with different seeds. The experiments are run on a Linux workstation with 16-core Intel Xeon 3.20GHz CPU and 128GB RAM.
C.4 BENCIIMAR METHOOSS
Existing causal structure learning methods used in our experiments all have available implementations, as listed below:
GES and PC: implementation of both methods is available through the py causal package at /t tittpps:i //githuh..ttbb...// com wttpccoooo.ceeeeeccoeoeeoeecooee bd2
We note that, the imple-
mentation py-causa] package is based on the CMU TETRAD project. in which the version of GES is indeed the fast GES algorithm proposed by Ramsey et al. (2017).
MMHC (Tsamardinos et al. 2006); an implementation available in the inneean package at htttss:/cRAA.. R- R-pro ject or. wttp:/aaaaaaccagannnnnn
18

---

Under review as a conference paper at ICLR 2021
CAM (Peters et al. 2014): its codes available through the CRAN R package repository at atttpss: cran. - p-pr wttp://eee.ttogo/////aootc.ccee
CAMM
NOTEARS (Zeng et al., 2018) and NOTEARS-MLP (Zheng et al, 2020 codes are available at the first author's s github repository https: /ttpts::/gitae thub. com/ xunzheng/
notears
GraN-DAG (Lachapelle al. 2020 an implementation is available at the first author's github repository ht / /t/////. /tt//githbbb com/ Graaooaaaaan/aaaaa-aaa DAG Note that for graphs of 50 nodes more, GraN-DAG performs preliminary neighborhood selection step to avoid overfitting.
DAG-GNN (Yu et al. 2019): the codes are available at the first author's github repository https://githut. com himmmfishmmmn 2200fffffmooo112134-------
ICA-LINGAM (Shimizu et al. 2006 an implementation is available at uthttpss//sitee
google.com/sitee wetesiteeseee sihee/sshing http://sshhmm2u0/iigggeee
In the experiments, we mostly use default hyperparameters unless otherwise stated.
D ADDITIONAL EXPERIMENTAL RESULTS
D.1 LINEAR SEMS WITH HIGHER RANKS
This experiment considers graphs of higher ranks. We use rank-specified random graphs with d 100 nodes and rank {30, 35, 40, 45, 50} on linear Gaussian SEMs. The results are shown in Figures 10a and 10b with degrees and 8, respectively. We observe that when the rank of the underlying graph becomes higher, the advantage of NOTEARS-low-rant over NOTEARS decreases. with rank r = 50 is still comparable to NOTEARS, and has a
Nonetheless, NOTEARS-low lower average SHD after removing outlier SHDs using the interquartile range rule.
20
150
NOTEARS
NOTEARS Mor-s-lowwrnkk.
940
5HD
NOTEARS low rank (wio outliers)
30
35
40
50
30
35
50
rank
rank
(a) Degree 2
(b) Degree
Figure 10: Average SHDs rank-specified graphs with higher ranks. The true rank is assumed to be known.
D.2 NOTEARS-LOW-RAL WITH DIFFERENT SAMPLE SIZES
We next empirically study the consistency of NOTEARS-
low-rank. Again, we use rank-specified random graphs (sampled according to Algorithm [) with 100 nodes, degree h 8, rank = 10, and linear Gaussian SEMs. also assume that the true rank is known. We fix the rank parameter f 10 and use different sample sizes ranging from 200 to 5. 000. From Figure II NOTEARS- low-rank performs reasonably well when the sample size small and tends to have better performance with a larger number of samples.
mean (wio outliers)
40
9:30
200
500
1000
2000
3000
5000
sample size
Figure 11: Different sample sizes.
19

---

Under review as a conference paper at ICLR 2021
D.3 FURTHER PRUNING
We compare the empirical results before and after applying the additional pruning technique described in Appendix C.2 The graphs are rank-specified with E {100, 300} nodes, rank r = [0.1d]. and degree l {2, 4, 6, 8). We again use linear Gaussian data model with equal noise variances to generate the datasets. The average SHDs are reported Figure [2] We see that applying an additional pruning step indeed improves the final performance of both NOTEARS and NOTEARS-low-rand especially on relatively large and dense graphs.
60-
50
250
NOTEARS, w/o pruning
200
NOTEARS-IoW-rank w/o pruning
150
NOTEARS
20
100
3
NOTEARS-low-rank
average degree
average degree
(a) 100 nodes, rank 10
(b) 300 nodes, rank 30
Figure 12: Average SHDs before and after pruning.
D.4 AN EMPIRICAL COMPARISON BETWEEN ICA-LINGAM AND DIRECTLINGAM
To our best knowledge, there two Python implementations of ICA-LINGAM (Shimizu al. 2006 released by the authors, available at pottps://sitte gowgle.oooooeeeeooooooooooooooo com www.oom/sst////e/ /sshimssiii zu 6/ lingam and https:: / thub. com http://t11///aaa/ respectively, where the latter is Python package containing several LiNGAM related methods. In the following, we use ICA- LiNGAM-pre and ICA-LiNGAM-cdt to denote these two implementations, respectively. For Di- rectLiNGAM (Shimizu et al. 201T), we only find a Python implementation available at the previously mentioned Python package containing ICA-LiNGAM-cde
Here we run DirectLiNGAM, ICA-LiNGAM-cdt, and ICA-LiNGAMMprr on linear exponential data models with 100-node and rank-10 graphs. The mean SHDs are reported below Table Based on this experimental result as well as our past experience, LireetLiNGAM usually has (slightly) better performance than ICA-LiNGAM-cdt while ICA-LiNGAM pre has a noticeably (if not much) better performance for relatively dense and large graphs. We are more concerned with relatively large and dense graphs and hence report the results achieved by ICA-LiNGAM-pu in the main paper.
Table 1: An empirical comparison between ICA-LiNGAM and DirectLiNGAM.
Degree
2
8
DireetLiNGAM
1.4
30.1
114.2
225.0
ICA-LiNGAM-cd
2.2
37.1
128.0
241.0
ICA-LINGAM-pre
7.0
31.7
61.7
72.6
D.S DETAILED EMPIRICAL RESULTS FOR EXPERIMENT 1 WITH LINEAR GAUSSIAN SEMS
Table reports detailed results including true positive rates (TPRs), false discovery rates (FDRs), structural Hamming distances (SHDs), and running time on rank aneciiecdf graphs with linear Gaussian data model. Here the true rank is assumed to be known and is used as the rank parameter NOTEARS- low-rank. We also test (fast) GES, MMHC, and PC. However, PC is too slow since some nodes may have a high in-degree (i.e., hubs) in large, dense, and low rank graphs. For the same reason, the skeleton may not be correctly estimated by MMHC, which has similar performance to that of GES. Therefore, we only include the results of GES for comparison. We treat GES favorably by regarding undirected edges positives if the true graph has directed edge in place of the undirected ones.
20

---

Under review as a conference paper at ICLR 2021
Table 2: Detailed results for linear Gaussian data model with equal noise variances.
100 nodes, rank 10
Degree
2
6
TPR
0.99 1 0.02
0.99 0.03
0,99 0.01
0.99 0.09
NOTEARS-
FDR
0.001 +: 0.003
0.01 0.00
0.02 1.004
0.02
low-rank
SHD
0.7 +: 1.6
4.7 10.2
8.5 13.7
9.15
Time (mins.)
1.9 +
2.7
5.8 2.6
7.5
TPR
0.90 0.06
0.00
0.05
0.89 0.03
FDR
0.07 0.03
0.04
0.06
0.04
NOTEARS
SHD
8..6
33.2
55.0 20.9
56.3 15.2
Time (mins.)
3.6 1.6
6.6 1.6
9.6
7.3
TPR
0.3 0.3
0,13 0.02
0.08 + 0.01
0.06 0.06++0000
FDR
0,03
0.82 0.02:0003
0.85 0.01
0.87 -.01
fast GES
SHD
141.6 1.1.6199.99 19.3
292.1 26.0
412.5 - 26.25
$21.9 22.3
Time
seconds
300 nodes, rank 30
Degree
2
TPR
0.99
0,99 0.01
0.99 0.01
0.97
NOTEARS- low-rank
0.03
0.02-0.00
0.04 0.02
0.06 0.08
SHD
12.0 10.6
13.9 14.8
50.0 = 32,1
106.0 169.5
Time (mins.)
57.6
76.5
158.7 = 94.7
262.9 144.4
TPR
0.01
0
0.93 = 0.02
0.91 0.02
FDR
0.07
0.06
0.08 = 0.04
0.09 0.06
NOTEARS
SHD
35.4 10.5
29.9
138.3 64.22
216.9 102.6
Time (mins.)
5.8
42.0 10.5
74.2
104.2
TPR
0.31 0.04
0.12
0.07 0.01
0.05 + 0.01
0.68 0.05
0.83=0.00
0.01
0.89 0.01
fast GES
SHD
427.0 58.1
883.3 60.5
1260.8 1 59.8
1608.4 67.5
Time
< 30 seconds
D.6 DETAILED RESULTS FOR EXPERIMENT WITH NON-LINEAR SEMS
Table 3 reports the detailed SHDs for each method in Section 5.4 We also mark in bold the best results from methods with or without low rank modifications.
Table 3: Detailed SHDs for Experiment 4 with non-linear SEMs.
Degree
DAG-GNN
50.1 8.2
96.9
182.0 14..9
CAM
1.99
19.4
98.91121.7
131.9
NOTEARS-MLP
6.48
47.8+200
86.0
119.422399
GraN-DAG
17.9
50.3
82.6
109.4 102.4
GraN-DAG, low rank, 0.3
23.2
45.8
65.9
87.1
GraN-DAG, low rank, 0.5
20.7
50.0
61.8+66.7
83.9
GraN-DAG, low rank, 1.0
42.7
68.7+76.2 67.9
83.2 76.9
GraN-DAG, low rank, validation
4.5
24.4
63.3
77.4 28.8
21


Reviews:
[{"review_id": "qwe5_DgHSFJ", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "# Summary\n\nThe paper develops several useful lower and upper bounds on the rank of DAGs — specifically minimum and maximum rank of all weighted matrices that induce the same DAG — in terms of various graphical properties like head-tail vertex cover, number of non-root and non-leaf vertices. The paper also bounds the rank of DAG in terms of the rank of its skeleton and moral graph. The paper proposes learning low-rank linear or non-linear structural equation models (SEMs) by adding simple norm constraints or matrix factorization to existing SEM learning methods. Through experiments on synthetic and real world data the authors demonstrate that when the underlying SEM is low-rank, exploiting this low-rank assumption in the learning process can lead to better performance. The authors also demonstrate that the rank can be estimated using the obtained bounds from a validation set.\n\n# Strengths\n\n1. The main contribution of the paper is a strong justification for learning SEMs under a low-rank assumption by showing that graphs with many hubs are low-rank. Existing theoretical results for learning SEMs show a polynomial dependence of the sample complexity on the maximum degree of the true SEM. Therefore, learning SEMs subject to rank constraints rather than sparsity constraints can be useful for graphs with hubs.\n2. The bounds on the rank of DAGs are generally useful beyond learning SEMs.\n\n# Weakness\n\n1. The paper does not propose any novel algorithms for learning low-rank DAGs, other than merely augmenting existing methods with nuclear norm constraint or using matrix factorization.\n2. The method for estimating the rank from the validation set is crude and computationally expensive.\n\n# Questions to address in rebuttal\n\n1. In Figure 2, is degree (x-axis) the maximum degree of a node graphs ?\n2. Figure 2 shows that the rank increases with the degree and that the rank is always larger than the degree. Therefore, even for graphs with hubs learning SEMs subject to sparsity constraints might still give better results than learning SEMs subject to rank constraints?\n3. More details are needed on how the rank is estimated from the validation set with a complete algorithm.\n\n# Post-rebuttal comments\nHello everyone,\n\nI have read the author's response and I am leaning towards rejection. The paper can be divided into two halves. The first half where the authors obtain bounds on ranks of DAGs is the main contribution of the paper and is clearly interesting. The second half of the paper tries to shoehorn these bounds into an algorithm for learning causal DAGs from observational data which is disappointing and is clearly below standard for the following reasons:\n\n1. The bounds depend on the underlying DAG which is unknown and therefore cannot be estimated from samples. Therefore the authors propose using \"structural priors\" to obtain these bounds. The authors don't mention where they get these structural priors from. Furthermore the bounds are only useful to restrict the hyper-parameter search space in the matrix factorization approach which is applicable to linear SEMs. These bounds can only be used \"qualitatively\" to guide selection of regularization penalty in the nuclear norm approach which is necessary for non-linear SEM methods.\n\n2. The theoretical results would still be useful if the authors could adequately demonstrate that for certain family of graphs the maximum degree can be high while the rank can be low therefore learning DAGs subject to sparsity constraints (whose sample complexity depend on the maximum degree) can perform worse than learning DAGs with rank constraints. However, this is not clear since in experiments the authors only show the SHD as a function of \"average degree\" and not \"maximum degree\". Figure 2 again compares rank against average degree and not maximum degree.\n\n3. The experiments are only performed in the low-dimensional regime at a fixed sample size (3000 samples and 300 nodes).\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "l0jRnhN39Ty", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "##########################################################################\n\nSummary:\n \nThe paper provides a new approch for learning a (possibly densely-connced) low-rank DAG models in the high dimensional settings. In particular, this paper provides how to exploit the property of the low-rank for recovering a underlying causal structure. Futher shown is that under what circumstance the low-rank assumption holds and heuristically confirms that thgrough simulation settings. Lastly, the proposed approach is compared against the state-of-the-art DAG learning algorithms that requirs the assumption of a sparse graph.\n\n##########################################################################\n\nReasons for score: \n\n \nOverall, I vote for accepting. This paper is well-written and delivers its main contribution really well. Futhermore it well summarizes the prior works on learning a causal graph. In addition, the main idea of recovering a graph under low-rank is novel. \nHowever, my major concern is about the simulations of the paper although I acknolwedge that most of relevant papers exploit a similar settings. It would be better to emphasize that the proposed algorithm attempts to learn a complete partial DAG, not a DAG. Although some of related paper fasely asserts that their approches recovers a DAG using conditional independence relationships or score function, I hope this paper clarifies this point.  \n\n \n##########################################################################Pros: \nPros:\n \n1. The paper solve a very important problem of causal inference. It seems to be practical and novel. \n \n2. The paper is really clear and convincing. \n \n##########################################################################\n\nCons: \nOne important comment from my side is that the way in which you simulate your models is severely biased. What I always do when I simulate models (and I think others should do something similar), is rescale edge weights for each\nnode such that if all parents have values with a standard-normal distribution, then the value of the node itself will also have a standard-normal distribution (assuming Gaussian additive noise). In this way one avoids that the variance of the variables blows up (or converges to 0) as one adds more and more nodes to the graph. Therefore, assuming a standard-normal error distribution (or error variance is large) is impractical.  \n\nFurthermore, in the densely-connected graph settings, it must be really careful to determine the range of edge weights; otherwise, the variance of the variables are again blowing up. Hence, in some points, the targeted graph is unrealistic in large-scale settings (d is large). Nevertheless, as an emerging field of learning DAG models in polynomial time with complete search, it should be accepted. However, for a better representation and fair comparison, it would be better to change simulation settings.  \n\nLastly, this paper does not explain a complete partial DAG that the proposed method is actually finds. In principle, there might be plenty of solutions for the considered optimization problem. Hence this paper would be clearer for new researcher in DAG model learning if it emphasizes CPDAG or PDAG. \n\n##### Update ######\nAlthough it is responded that the simulation setting used in the paper does cause blowing up samples or marginal variance, it is in general impossible or the setting assumes too sparse case where considered graphs are almost empty. In addition, as you mentioned, I also ackowledge that it is a widely-used setting; however, there are a lot of papers that are rejected because of the unfair simulation setting. I like the main idea of the paper a lot, and hence, I hope the authors set the simulation setting more carefully. \n\nFurthermore, it is really frustaring answer that the authors consider the only case where graph is uniquely idenfiable from the pure observations. As you know that is really rare when the number of nodes is large (p > 50). ", "extra_scores": {}, "extra_text": {}}, {"review_id": "52hC-S8MdUx", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "# Review\n## Summary of the paper and the review\n- Paper summary: The paper explores the possibility of exploiting the potential low-rank nature of the underlying causal graph when conducting causal structure learning. Building on the recent gradient-based methods in (causal) DAG selection, they describe ways of integrating this assumption into the algorithms for existing gradient-based methods that assume linear or non-linear structural equation models (SEMs). They also provide some upper and lower bounds to the graph rank given the knowledge of certain structural priors. They test their proposal by augmenting two existing gradient-based methods and comparing its performance to the originals and other methods in simulated data.\n\n- Review summary: The authors examine an intriguing and timely phenomenon, and provide some engaging arguments and results. However the study has some problems in its argumentation and experimentation. I have listed the perceived strengths and weaknesses of the paper below, including my questions for further clarification. I am looking forward to authors' improvements and responses regarding these issues.\n\n## Strengths of the study\n- The field of causal discovery can be considered as much a discussion regarding the nature of causality and how it manifests itself in observed joint distributions as it is a discussion regarding specific methods. I believe that it is very important that this work investigates the properties arising from the nature of the underlying distributions and how this should relate to the models and algorithms.\n- The authors discuss the structure of causal-DAGs in relation to widely observed complex network structures such as scale-free networks. Although the connections they put forth are not always very clear (see below), I believe that discussing the nature of causal structures in the real world is an important direction, since causal discovery inevitably requires certain structural priors (e.g. independent causal mechanisms).\n- Gradient-based DAG discovery methods constitute a promising direction towards causal structure discovery, especially in the presence of a large number of variables. The authors propose a method that can potentially improve any such method with relative ease.\n- The modification the authors suggest comes with acceptable computational cost. Considering an extensive potential hyperparameter search, this cost is not negligible, yet it can still be said that it only adds a scaling factor to the complexity.\n- Their method can be applied to both linear and nonlinear SEM's. Although the conceptual interpretations in the two case do not always translate, this is more due to the differing nature of the optimization problem in the two cases.\n- The authors provide lower and upper bounds for the rank estimations of the causal systems given the existence of certain structural priors.\n- The authors conduct simulated experiments where they vary a number of experimental conditions.\n- The paper is well-written.\n\n## Weaknesses of the study\n### The claimed ubiquity of low-rank structure in causal systems: \n- The authors methods aim to exploit underlying low-rankness in the (adjacency matrix of the) underlying causal graphs. When arguing when real life causal systems are likely to be low rank, the authors invoke the prevalence of scale-free networks in real life. This argument is confusing since the observation of scale-free networks are frequently made in relational contexts where a rich-get-richer phenomena are frequently assumed to drive the scale-free property (e.g. a popular social media personality becoming more and more connected as opposed to a typical user). \n    - I fail to see why we would also see this in a causal-DAG setting as well, why would we have a causal hub? The authors write that: \"It is observed that many real-world networks are scale-free, and some of them may be viewed as causal networks\" - the authors could at least present some examples from the literature (in addition to the Pathfinder and arth150 data they cite) as well as more reasoning to make a more persuasive case. \n- This might have been less of an issue if the authors showed that low-rank assumption would benefit performance even in the case of high rank causal networks, but their results do not provide strong arguments to that effect (Appendix D.1).\n- However, though being a weak point in their argumentation, this does not rule out the potential benefits of their methods, since the fact that the low-rankness assumption does not hold can be established in the cross-validation phase of a study using their method. \n\n### Upper and lower bounds on the rank of the adjacency matrix\n- The paper devotes a whole section to extracting lower and upper bounds from prior knowledge regarding the structure of the graph, however provide little as to how we would realistically have access to such information. Almost none of the structural information required by the theorems seem to be easily accessible enough. The authors mention that \"we usually have access to some structural information\" but do not go into detail.\n- Even when we have some structural information, it is not always clear how this information can be used in constraining the search space when this information is anything less than the knowledge of the causal structure itself. For example, even if we knew that our causal network had a number of causal hubs and thus is likely to be scale-free (which we perhaps could in some cases), it still is not clear how this knowledge could be directly used to obtain upper bounds. The relation is even more unclear for the nonlinear case since there is no direct way of relating $\\hat{r}$ and $\\lambda$.\n- This might not have been as big an issue if this was not a section supporting a central claim in an application oriented paper - indeed, the theorems are presented as guidelines to constrain the search space for the estimated rank of the adjacency matrix, so more discussion and demonstration of their practical use is naturally sought.\n- I believe that the authors could do two things:\n    - They could show how the prior information required in Section 4 can actually plausibly be known beforehand in realistic scenarios, and maybe show example(s) of this in one or more of their applications (unlike in Section 5.3 where the prior information is assumed to be known).\n    - Or they could make it clear that these are, though potentially useful as is or after contribution by other research, mostly theoretical demonstrations, and as of yet, their method would likely require extensive hyperparameter search as to the rank of the decomposition or the nuclear norm coefficient. This could possibly require them reevaluate their emphasis on these theoretical findings when introducing their research.\n\n### Experiments\n- Given the above concerns, the experiments take on more importance in demonstrating the utility of the methods proposed. Though the experiments have some important findings, a stronger experiments section would allow the authors to make a stronger claim and base their argument more empirically. At least _some_ additional experiments should be feasible given the running times of the algorithms provided in the Appendix.\n- Given that it is not clear why causal networks would be both low rank and dense, and given that the prior information to constrain the search space for the rank are not always accessible, I believe that a comparison with sparsity-inducing regularizers such as $\\ell_1$ or $\\ell_2$-regularizers  would be beneficial, with a correspondingly extensive hyperparameter search so that empirical superiority of the authors' method over others could be empirically established (This has only been conducted in Section 5.5, see below regarding these findings). Given that both would require extensive parameter search, for a practitioner it might not matter whether their regularizer is sparsity or low-rank-decomposition inducing.\n- The authors conducted experiments by applying their methods to modify NOTEARS and GraN-DAG. In the absence of application of their modification to other methods, justification for their method choice and why similar improvements should be expected for others as well would be helpful if provided. \n- Additional experiments when the true rank goes beyond $d/2$ (for both linear and nonlinear SEMs) would be beneficial in understanding the behavior of the given method in such cases. \n    - Also, this would demonstrate a potential contrast between the two methods: in the nonlinear case the cross-validation could produce a very small $\\lambda$ to practically disregard the low-rank assumption. However in the linear case, e.g. setting $\\hat{r} = r$ would lead to full rank assumption but would lead to $2d^2$ parameters, potentially leading to worse results than the original algorithm both in terms of performance and computation time. \n    - In the cases where $d/2 \\ll r$, could score comparison e.g. with an unfactorized NOTEARS method lead to the correct choice (that is, full-rank NOTEARS)?\n- Is there a specific reason for the absence of results with alternative $k$'s and with ICA-LiNGAM in Linear SEM experiments in scale-free networks in Section 5.2? If not, the inclusion of these would also be illustrative, especially given the latter's comptetive performance in previous experiments.\n- Given the centrality of scale-free graphs in their argumentation, is there a specific reason for the lack of nonlinear-SEM experiments with scale-free graphs?\n- Section 5.3: As noted before, it might be hard in many cases to provide realistic lower and upper bounds for the rank values. This might require an extensive search for $\\hat{r}$ and $\\lambda$. Given that the method requires a potentially large parameter search, its effect to general algorithm time complexity must also be mentioned, e.g. its runtime being at least 10x that of the base algorithm. However, I believe that this extra computation cost could be acceptable for increase in structure discovery accuracy.\n- The results presented in Figure 7 seem to include very wide interquartile ranges: are the differences between different methods/hyperparameters significant? \n    - Also in the description of this experiment, the authors note that the results show the \"utility of the low-rank assumption\" when the rank is not as low as assumed, but the results show almost equal results, and given that the assumed rank is not much lower than the real rank, I cannot see how this result supports their conclusions, especially given that this is conclusion is presented in the abstract.\n\n## Other comments\n- In Figures 2, 3, and 6 the horizontal axis reads \"Degree\" and the corresponding text mentions \"graph degree\", in Sections 4.3 and 5.1 . The authors must have meant average (in + out) degree as they did in Section 5.2 or Appendix C.1. I believe the paper would benefit from this being explicit everywhere. In Appendix A the authors seem to equate \"graph degree\" with total edges in the graph, but they must mean the average number of edges per node in the graph, so that it is consistent with the rest of the paper.\n- In Figure 2a, is average NOTEARS (w/o outliers) performance non-existent in the figure or follow exactly the same trajectory as NOTEARS? \n- At Pg. 6 in Experiments section the authors might have wanted to cite (Zheng et al. 2018) for the algorithm NOTEARS.\n- Figure 5 would benefit from adding the NOTEARS result from the Figure 3.\n- Is there a specific reason for authors preferring ICA-LiNGAM over the following DirectLiNGAM?\n- Parameter search spaces: In this case $\\hat{r}$ would have to be searched for in $[1, d/2]$ (vs. the original algorithm with no low-rank assumption, see above), however it is unclear what the search space for $\\lambda$ should be.  Also the authors would benefit from suggesting a search strategy for such a case.\n- As a future work, the research could benefit from comparing the inductive bias of the current method to the literature on structure priors [1] or parameter priors affecting model selection through marginal likelihood [2] in DAGs. \n\n## References\n1. Eggeling, Ralf, Jussi Viinikka, Aleksis Vuoksenmaa, and Mikko Koivisto. 2019. “On Structure Priors for Learning Bayesian Networks.” In _The 22nd International Conference on Artificial Intelligence and Statistics_, 1687–95. [http://proceedings.mlr.press/v89/eggeling19a.html](http://proceedings.mlr.press/v89/eggeling19a.html).\n2. Silander, Tomi, Petri Kontkanen, and Petri Myllymaki. 2007. “On Sensitivity of the MAP Bayesian Network Structure to the Equivalent Sample Size Parameter.” in _Proceedings of the Twenty Third Conference on Uncertainty in Artificial Intelligence_.", "extra_scores": {}, "extra_text": {}}, {"review_id": "2G4gJTC1fmm", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper attempts to exploit the low-rankness of the adjacency matrix of the DAG in Bayesian network structure learning. The overall framework is similar to NOTEARS, except that the adjacency matrix W is decomposed into low rank components W = UV'. To justify the approach, the paper also includes lower and upper bounds of the rank of DAGs, albeit mostly theoretical and not applicable to real experiments. \n\nThe paper is very solid in presenting mathematical facts and detailed algorithms. However, my main concern is about the fact that the algorithm requires knowledge (or guess) about the rank. In fact, the experiments in Section 5 already uses the ground truth rank information in NOTEARS-low-rank. Algorithm 1 is a great resource to be shared in the community, however in principal it shouldn't be needed to perform the experiments in Section 5. If one can gain accuracy benefit even without knowing the true rank, paying extra runtime cost is acceptable (Table 1). \n\nThere are also many works on combining low-rankness with sparsity, which I suggest the authors to consider as future steps.\n\nUpdate:\nThe authors have explained the issue raised in the review. It's not ideal that the algorithm requires the knowledge of rank beforehand, but it's okay if this point is clearly communicated in the paper. I would keep my current score. ", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies the low-rank properties of DAG models, and illustrates through proof-of-concept how low-rank-ness can be exploited in structure learning of DAGs. After a lengthy discussion amongst the reviewers, it became clear that although there are some interesting ideas here, there is not enough enthusiasm for this work in its current form. The results in Section 4 connecting rank to structural properties are interesting, but the reviewers were concerned by the lack of precise statements connecting these results to known ensembles such as scale-free graphs (even though the authors discuss some heuristic connections). In the end, despite considerable enthusiasm regarding these ideas and the importance of the problem studied, there remained too many concerns that require a major revision before acceptance.
Decision: Reject