summaries: 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.  The 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. The 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.
abstractText: 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.
references.authors: Albert-Laszlo Barabasi Zoltan N Oltvai
references.title: Network biology: Understanding the cell’s functional organization
references.venue: Nature Reviews Genetics,
references.authors: Adarsh Barik Jean Honorio
references.title: Learning Bayesian networks with low rank conditional probability tables
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Dimitri P Bertsekas
references.title: Nonlinear Programming
references.venue: Athena Scientific,
references.authors: Peter Bühlmann Jonas Peters Jan Ernest
references.title: CAM: Causal additive models, high-dimensional order search and penalized regression
references.venue: The Annals of Statistics,
references.authors: Xiangyong Cao Yang Chen Qian Zhao Deyu Meng Yao Wang Dong Wang Zongben Xu
references.title: Low-rank matrix factorization under general mixture noise distributions
references.venue: In The IEEE International Conference on Computer Vision (ICCV),
references.authors: David Maxwell Chickering
references.title: Learning Bayesian networks is NP-complete
references.venue: In Learning from Data: Artificial Intelligence and Statistics V. Springer,
references.authors: David Maxwell Chickering
references.title: Optimal structure identification with greedy search
references.venue: Journal of Machine Learning Research,
references.authors: Mark A. Davenport Justin Romberg
references.title: An overview of low-rank matrix recovery from incomplete observations
references.venue: IEEE Journal of Selected Topics in Signal Processing,
references.authors: Kőnig Dénes
references.title: Gráfok és mátrixok
references.venue: Matematikai és Fizikai Lapok,
references.authors: Jack Edmonds
references.title: Systems of distinct representatives and linear algebra
references.venue: Journal of Research of the Notional Bureau of Standards - B. Mathematics and Mathematical Physics,
references.authors: Ralf Eggeling Jussi Viinikka Aleksis Vuoksenmaa Mikko Koivisto
references.title: On structure priors for learning Bayesian networks
references.venue: In International Conference on Artificial Intelligence and Statistics (AISTATS),
references.authors: Victor M. Eguı́luz Dante R. Chialvo Guillermo A. Cecchi Marwan Baliki A. Vania Apkarian
references.title: Scale-free brain functional networks
references.venue: Physical Review Letters,
references.authors: Shaun M Fallat Leslie Hogben
references.title: The minimum rank of symmetric matrices described by a graph: a survey
references.venue: Linear Algebra and its Applications,
references.authors: Benjamin Frot Preetam Nandy Marloes H. Maathuis
references.title: Robust causal structure learning with some hidden variables
references.venue: Journal of the Royal Statistical Society: Series B (Statistical Methodology),
references.authors: Bo Gao Ruo en Ren
references.title: The topology of a causal network for the Chinese financial system
references.venue: Physica A: Statistical Mechanics and its Applications,
references.authors: Nabil Guelzim Samuele Bottani Paul Bourgine François Képès
references.title: Topological and causal structure of the yeast transcriptional regulatory network
references.venue: Nature Genetics,
references.authors: Alexander J Hartemink
references.title: Reverse engineering gene regulatory networks
references.venue: Nature Biotechnology,
references.authors: Yangbo He Jinzhu Jia Bin Yu
references.title: Counting and exploring sizes of Markov equivalence classes of directed acyclic graphs
references.venue: Journal of Machine Learning Research,
references.authors: David Earl Heckerman Eric J Horvitz Bharat N Nathwani
references.title: Toward normative expert systems: Part i the pathfinder project
references.venue: Methods of Information in Medicine,
references.authors: Leslie Hogben
references.title: Minimum rank problems
references.venue: Linear Algebra and its Applications,
references.authors: Cho-Jui Hsieh Kai-Yang Chiang Inderjit S. Dhillon
references.title: Low rank modeling of signed networks
references.venue: In the 18th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining,
references.authors: Jin Huang Feiping Nie Heng Huang Yu Lei Chris Ding
references.title: Social trust prediction using rank-k matrix recovery
references.venue: In International Joint Conference on Artificial Intelligence (IJCAI),
references.authors: Diviyan Kalainathan Olivier Goudet Isabelle Guyon David Lopez-Paz Michèle Sebag
references.title: Structural agnostic modeling: Adversarial learning of causal graphs
references.venue: arXiv preprint arXiv:1803.04929,
references.authors: Nan Rosemary Ke Olexa Bilaniuk Anirudh Goyal Stefan Bauer Hugo Larochelle Bernhard Schölkopf Michael C. Mozer Chris Pal Yoshua Bengio
references.title: Learning neural causal models from unknown interventions
references.authors: Diederik P Kingma Max Welling
references.title: Auto-encoding variational bayes
references.venue: In International Conference on Learning Representations (ICLR),
references.authors: Vladimir Koltchinskii Karim Lounici Alexandre B. Tsybakov
references.title: Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion
references.venue: The Annals of Statistics, 39(5):2302–2329,
references.authors: Sébastien Lachapelle Philippe Brouillard Tristan Deleu Simon Lacoste-Julien
references.title: Gradient-based neural DAG learning
references.venue: In International Conference on Learning Representations (ICLR),
references.authors: Po-Ling Loh Peter Bühlmann
references.title: High-dimensional learning of linear causal networks via inverse covariance estimation
references.venue: The Journal of Machine Learning Research,
references.authors: Christopher Meek
references.title: Causal inference and causal explanation with background knowledge
references.venue: In UAI,
references.authors: Milena Mihail Christos Papadimitriou
references.title: On the eigenvalue power law
references.venue: In Randomization and Approximation Techniques in Computer Science,
references.authors: Lon H. Mitchell Sivaram K. Narayan Andrew M. Zimmer
references.title: Lower bounds in minimum rank problems
references.venue: Linear Algebra and its Applications,
references.authors: Ignavier Ng Zhuangyan Fang Shengyu Zhu Zhitang Chen Jun Wang
references.title: Masked gradient-based causal structure learning
references.venue: arXiv preprint arXiv:1910.08527,
references.authors: Ignavier Ng Shengyu Zhu Zhitang Chen Zhuangyan Fang
references.title: A graph autoencoder approach to causal structure learning
references.venue: arXiv preprint arXiv:1911.07420,
references.authors: Olga Nikolova Srinivas Aluru
references.title: Parallel Bayesian network structure learning with application to gene networks
references.venue: In Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis,
references.authors: Jonas Peters Peter Bühlmann
references.title: Identifiability of Gaussian structural equation models with equal error variances
references.authors: Jonas Peters Joris M Mooij Dominik Janzing Bernhard Schölkopf
references.title: Causal discovery with continuous additive noise models
references.venue: The Journal of Machine Learning Research,
references.authors: Jonas Peters Dominik Janzing Bernhard Schölkopf
references.title: Elements of Causal Inference - Foundations and Learning Algorithms
references.authors: Joseph Ramsey Madelyn Glymour Ruben Sanchez-Romero Clark Glymour
references.title: 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
references.venue: images. International Journal of Data Science and Analytics,
references.authors: Benjamin Recht
references.title: A simpler approach to matrix completion
references.venue: Journal of Machine Learning Research,
references.authors: Shohei Shimizu Patrik O Hoyer Aapo Hyvärinen Antti Kerminen
references.title: A linear non-Gaussian acyclic model for causal discovery
references.venue: Journal of Machine Learning Research,
references.authors: Shohei Shimizu Takanori Inazumi Yasuhiro Sogawa Aapo Hyvärinen Yoshinobu Kawahara Takashi Washio Patrik O Hoyer Kenneth Bollen
references.title: Directlingam: A direct method for learning a linear non-Gaussian structural equation model
references.venue: Journal of Machine Learning Research,
references.authors: Tomi Silander Petri Kontkanen Petri Myllymaki
references.title: On sensitivity of the map Bayesian network structure to the equivalent sample size parameter
references.venue: In Uncertainty in Artificial Intelligence (UAI),
references.authors: Carl Smith Frank Wood Liam Paninski
references.title: Low rank continuous-space graphical models
references.venue: In International Conference on Artificial Intelligence and Statistics (AISTATS),
references.authors: Peter Spirtes Clark N Glymour Richard Scheines
references.title: Causation, Prediction, and Search
references.venue: MIT Press, second edition,
references.authors: Petr Tichavskỳ Jiřı́ Vomlel
references.title: Representations of bayesian networks by low-rank models
references.venue: In International Conference on Probabilistic Graphical Models,
references.authors: Ioannis Tsamardinos Laura E Brown Constantin F Aliferis
references.title: The max-min hill-climbing Bayesian network structure learning algorithm
references.venue: Machine learning,
references.authors: Sara Van de Geer Peter Bühlmann
references.title: `0-penalized maximum likelihood for sparse directed acyclic graphs
references.venue: The Annals of Statistics,
references.authors: Quanming Yao James T. Kwok
references.title: Greedy learning of generalized low-rank models
references.venue: In International Joint Conference on Artificial Intelligence,
references.authors: Yue Yu Jie Chen Tian Gao Mo Yu
references.title: DAG-GNN: DAG structure learning with graph neural networks
references.venue: In International Conference on Machine Learning (ICML),
references.authors: Jiawei Zhang Jianhui Chen Shi Zhi Yi Chang Philip S. Yu Jiawei Han
references.title: Link prediction across aligned networks with sparse and low rank matrix estimation
references.venue: In IEEE 33rd International Conference on Data Engineering,
references.authors: Jiji Zhang
references.title: On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias
references.venue: Artificial Intelligence,
references.authors: Kun Zhang Aapo Hyvärinen
references.title: On the identifiability of the post-nonlinear causal model
references.venue: In Uncertainty in Artificial Intelligence (UAI),
references.authors: Xun Zheng Bryon Aragam Pradeep Ravikumar Eric P. Xing
references.title: DAGs with NO TEARS: Continuous optimization for structure learning
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Xun Zheng Chen Dan Bryon Aragam Pradeep Ravikumar Eric P Xing
references.title: Learning sparse nonparametric DAGs
references.venue: In International Conference on Artificial Intelligence and Statistics (AISTATS),
references.authors: Ke Zhou Hongyuan Zha Le Song
references.title: Learning social infectivity in sparse low-rank networks using multi-dimensional Hawkes processes
references.venue: In International Conference on Artificial Intelligence and Statistics (AISTATS),
references.authors: Shengyu Zhu Ignavier Ng Zhitang Chen
references.title: Causal discovery with reinforcement learning
references.venue: In International Conference on Learning Representations (ICLR),
references.title: Rd×r̂ and ◦ is the point-wise product. The constraint in Problem (5) holds if and only if UV T is a weighted adjacency matrix of a DAG. This problem can then be solved by standard numeric optimization methods such as the augmented Lagrangian method (Bertsekas
references.authors: − d
references.title: α is the Lagrange multiplier, and ρ > 0 is the penalty parameter. The optimization procedure is summarized in Algorithm 2, similar to Zheng et al
references.authors: Zheng
references.title: 2018) to perform consistently well even when the graph is sparse; (3) we are more concerned with relatively large and dense graphs, so a sparsity assumption may be harmful
references.authors: Yu
references.title: r̂ columns of the d× d identity matrices. Other parameter choices are: α0 = 0, c = 0.25, = 10−6, and w = 0.3, similar to those used in related methods on the same datasets (e.g., Zheng et al
sections.heading: 1 INTRODUCTION
sections.text: 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 to 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 acyclic graphs (DAGs) over a set of random variables.
  For this task, existing methods can be roughly categorized into two classes: constraint- and scorebased. The former use statistical tests to 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, 1995; 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 a class of causal structures) with the optimal score. Due to the combinatorial nature of the acyclicity constraint (Chickering, 1996; He et al., 2015), most score-based methods rely on local heuristics to 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 on 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. (2019); Lachapelle et al. (2020); Zheng et al. (2020), 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 in this paper. This issue is of course a 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 et al., 2006) relies on local learning methods that often do not perform well when the target node has a large neighborhood. How to improve the performance on relatively large and dense DAGs is therefore an important question.
  In this work, we study the potential of exploiting a 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 a sparsity assumption on graph structures, we 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 a graph to its graphical structure. Such a result also enables us to characterize the rank of a 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 a 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 networks tend to be low rank. From this result, we derive several graphical conditions to bound the rank of a 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., 2019), matrix completion and factorization (Recht, 2011; Koltchinskii et al., 2011; Cao et al., 2015; Davenport & Romberg, 2016), network sciences (Hsieh et al., 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); Tichavskỳ & 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 directed/undirected graph. However, most works only considered the zero-nonzero pattern of off-diagonal entries (see, e.g., Fallat & Hogben (2007); Hogben (2010); Mitchell et al. (2010)), 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, Edmonds (1967) studied algebraically the maximum rank for matrices with a common zero-nonzero pattern. In Section 4, we use this result to relate the maximum possible rank to a more interpretable graphical condition, which further implies several structural conditions of DAGs that may be easier to obtain in practice.
sections.heading: 2 PRELIMINARIES
sections.heading: 2.1 GRAPH TERMINOLOGY
sections.text: A graph G is defined as a pair (V,E), where V = {X1, X2, · · · , Xd} is the vertex set and E ⊂ V2 denotes the edge set. We are particularly interested in directed (acyclic) graphs in the context of causal structure learning. For any S ⊂ V, we use pa(S,G), ch(S,G), and adj(S,G) 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 a 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 adjacency matrices. Specifically, the weighted adjacency matrix of a weighted graph G is a matrix W ∈ Rd×d, where W (i, j) is the weight of edge Xi → Xj and W (i, j) 6= 0 if and only if Xi → Xj exists in G. The binary adjacency
  matrix A ∈ {0, 1}d×d is such that A(i, j) = 1 if Xi → Xj in G and A(i, j) = 0 otherwise. The rank of a weighted graph is defined as the rank of the associated weighted adjacency matrix.
sections.heading: 2.2 CAUSAL STRUCTURE LEARNING AND RECENT GRADIENT-BASED METHODS
sections.text: 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 Xi’s to denote random variables associated with the nodes in a graph G. Assuming G being a DAG, then the SEM is given by
  Xi = fi (pa(Xi,G), i) , i = 1, 2, . . . , d, where fi is a deterministic function and i’s are jointly independent noises. The SEM induces a marginal distribution P (X) over X = [X1, X2, · · · , Xd]T , and G and P (X) are said to form a causal Bayesian network (Pearl, 2009; Spirtes et al., 2000). The problem of causal structure learning is to infer the underlying causal DAG G based on the marginal distribution P (X), or more practically, an empirical version consisting of a number of i.i.d. observations from P (X).
  We next briefly review recently developed gradient-based methods that rely on a smooth characterization of acyclicity of directed graphs. These methods aim to find a 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 a weighted adjacency matrix W and then estimate W by solving optimization problems in the following form:
  min W,φ
  EX∼P (X) S ( X,h(X;W,φ) ) , subject to trace ( eW◦W ) − d = 0, (1)
  where h : Rd → Rd is a model function parameterized by W (and other possible parameter φ) that aims to reconstruct X , S(·, ·) denotes a score function between the true and reconstructed variables, notation ◦ denotes the element-wise product, and eM is the matrix exponential of a square matrix M . The constraint was proposed by Zheng et al. (2018), which is smooth and holds if and only if W indicates a DAG. Methods in this class include: NOTEARS (Zheng et al., 2018), which targets linear models, with h(X;W,φ) = WTX 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 can 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 to or are variants of some well studied score functions like the penalized maximum likelihood (Chickering, 2002; Van de Geer et al., 2013; Loh & Bühlmann, 2014).
  The second class uses certain functions, with parameter θ, to construct a weighted adjacency matrix W (θ) (or a binary one A(θ)) to represent the causal structure. These methods can be summarized as
  min θ, φ
  EX∼P (X) S ( X,h(X;W (θ), φ) ) , subject to trace ( eW (θ) ◦W (θ) ) − d = 0. (2)
  For example, GraN-DAG (Lachapelle et al., 2020) and NOTEARS-MLP (Zheng et al., 2020) respectively use neural network path products and partial derivatives between variables to construct W (θ). The binary matrix A(θ) can be obtained by sampling according to some distributions with learnable parameters, as used by Kalainathan et al. (2018); Ke et 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., 2000; Peters et al., 2017). For such cases, one may convert the obtained DAG to its corresponding Completed Partially Directed Acyclic Graph (CPDAG) as the estimate. Nevertheless, if the SEM is identifiable and a 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. (2006); Peters & Bühlmann (2013); Peters et al. (2014); Zhang & Hyvärinen (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.
sections.heading: 3 EXPLOITING LOW RANK ASSUMPTION IN CAUSAL STRUCTURE LEARNING
sections.text: This section shows how to adapt existing gradient-based methods to take advantage of the low rank assumption, by providing a way for each class to utilize this assumption using techniques from the
  matrix completion literature. We remark that our adaptations with the low rank assumption are 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 an estimate r̂ for the graph rank, we can factorize W as W = UV T with U, V ∈ Rd×r̂. Problem (1) 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 T ) as the original one if r̂ is greater than or equal to the true rank. Furthermore, if r̂ d, we have a much reduced number of parameters to optimize. Nuclear Norm For the second class of methods, the adjacency matrix W (θ) is not an explicit parameter to be optimized. In such a case, we can adopt a commonly used technique to add a nuclear norm term λ‖W (θ)‖∗, with λ > 0 being a tuning parameter, to the objective to induce low-rankness. The optimization procedures in these recent structure learning methods can directly incorporate the two adaptations as they are all 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 to the singular value decomposition to compute the (sub-)gradient w.r.t. W at each optimization step.
  An acute reader may have noticed that we assumed a proper rank estimate r̂ or a proper penalty parameter λ. 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 r̂ 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, a lower bound rl and an upper bound ru on the graph rank would be beneficial—we need only consider ranks in [rl, ru] 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, a 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.
sections.heading: 4 GRAPHICAL BOUNDS ON RANKS
sections.text: 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 a 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.
sections.heading: 4.1 PROBLEM SETTING
sections.text: 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 rank(W ) using only the graphical structure. Specifically, we focus on the weighted adjacency matrices with the same binary adjacency matrix A, i.e.,WA = {W ∈ Rd×d ; sign(|W |) = 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. In the following, we investigate the maximum rank max{rank(W );W ∈ WA} and minimum rank min{rank(W );W ∈ WA} to find tighter upper and lower bounds for any W ∈ WA. Before introducing two useful graph concepts, we comment that low rank DAGs are not necessarily sparse and vice versa; see a discussion in Appendix A.
  Definition 1 (Height). Given a DAG G = (V,E) and a vertex Xi ∈ V, the height of Xi, denoted by l(Xi), is defined as the length of the longest directed path starting from Xi. The height of G, denoted by l(G), is the length of the longest path in G.
  Definition 2 (Head-tail vertex cover). Let G = (V,E) be a directed graph and H,T be two subsets of V. (H,T) is called a head-tail vertex cover of G if every edge in G has its head vertex in H or its tail vertex in T. The size of a head-tail vertex cover (H,T) is defined as |H|+ |T|.
  As an example, Figure 1c is a head-tail vertex cover of G in Figure 1a, where H = {X2, X4, X8} (red nodes) and T = {X8, X9, X10} (blue nodes). The size of this vertex cover is 6.
sections.heading: 4.2 LOWER BOUNDS
sections.text: We first study lower bounds on the rank of a weighted DAG. Define V−1 = ∅ and Vs = {Xi; l(Xi) = s} for s = 0, 1, . . . , l(G). Denote by Gs,s−1 the induced subgraph of G over Vs∪Vs−1. Let C(Gs,s−1) be the set of non-singleton connected components of Gs,s−1 and |C(Gs,s−1)| the cardinality. We have the following lower bounds. Theorem 1. Let G be a DAG with binary adjacency matrix A. Then
  min{rank(W ) ; W ∈ WA} ≥ ∑l(G)
  s=1 |C(Gs,s−1)| ≥ l(G). (3)
  All the proofs in this paper are provided in Appendix B. Theorem 1 shows that rank(W ) is greater than or equal to the sum of the number of non-singleton connected components in each Gs,s−1. As Gs,s−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 as the length of the longest directed path. As an example, consider the graph shown in Figure 1. One can verify that min{rank(W );W ∈ WA} = 6, |C(G1,0)| = 2, |C(G2,1)| = 1, |C(G3,2)| = 1, and l(G) = 3. Thus, we have min{rank(W );W ∈ WA} = 6 > 2 + 1 + 1 = 4 > 3. We remark that the bounds in Theorem 1 may be loose in some cases. To characterize the minimum rank exactly is an on-going research problem (Hogben, 2010).
sections.heading: 4.3 UPPER BOUNDS
sections.text: We turn to the more important issue for our purpose, regarding upper bounds on rank(W ). The next theorem shows that max{rank(W );W ∈ WA} can be characterized exactly in graphical terms. Theorem 2. Let G be a directed graph with binary adjacency matrix A. Then 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 ∈ 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 a minimum head-tail vertex cover, which in general is 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 of a relatively small size then indicates the presence of hubs, that is, vertices with relatively high in-degrees or out-degrees. Therefore, Theorem 2 suggests that the maximum rank of a weighted DAG is highly related to the presence of hubs: a DAG with many hubs tends to have low rank. Intuitively, a hub of high in-degree (out-degree) is a common effect (cause) of a number of direct causes (effect variables), comprising many V-structures (inverted V-structures). For example, in Figure 1a, X8 is a hub of V-structures and X9 is a hub of inverted V-structures.
  Such features are fairly common in real graph structures. Appendix A presents a real network, called pathfinder, which describes the causal relations among 109 variables (Heckerman et al., 1992) with the center node being the parent of a large number of other nodes. The famous scale-free (SF) graphs also tend to have hubs. A scale-free graph is one whose distribution of degree k follows a power law: P (k) ∼ k−γ , where γ is the power parameter typically within [2, 3] and P (k) 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 scale-free graphs are relatively low, especially in comparison to Erdös-Rényi (ER) random graphs (Mihail & Papadimitriou, 2002). Figure 2 provides a simulated example where γ 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 γ, the increase of their ranks is much slower than that of Erdös-Rényi graphs; indeed, their ranks tend to stay fairly low even when the graph degree is large.
  Theorem 2 can also be used to generate a low rank graph, or more precisely, a random DAG with a given rank r and a properly specified graph degree. Here we briefly describe the idea and leave the detailed algorithm to Appendix C.1: first generate a graph with r edges and rank r; a random edge is sampled without replacement and would be 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 a continuous distribution and the weighted graph will have rank r with high probability.
  The next two theorems report some looser but simpler upper bounds on rank(W ).
  Theorem 3. Let G be a DAG with binary adjacency matrix A, and denote the set of vertices with at least one parent by Vch and those with at least one child by Vpa. Then we have
  max{rank(W ) ; W ∈ WA} ≤  ∑l(G) s=1 min (|Vs|, |ch(Vs)|) ≤ |Vpa|,∑l(G)−1 s=0 min (|Vs|, |pa(Vs)|) ≤ |Vch|,
  |V| −max{|Vs| ; 0 ≤ s ≤ l(G)}.
  (4)
  Since Vch and Vpa 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 a 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 a 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 ∈ WA} ≤max{rank(W ) ; sign(|W |) = skeleton(A)} ≤max{rank(W ) ; sign(|W |) = moral(A)}.
  The skeleton of a 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 a 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 one as our estimate.
sections.heading: 5 EXPERIMENTS
sections.text: 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;1 and DAG-GNN (Yu et al., 2019), NOTEARS-MLP (Zheng et al., 2020), and CAM (Bühlmann et al., 2014) for the 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 1 in Appendix C.1), scale-free graphs, and a real network structure. For linear SEMs, the weights are uniformly sampled from [−2,−0.5] ∪ [0.5, 2] and the noises are either standard Gaussian or standard exponential. For 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 & Bü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.
sections.heading: 5.1 LINEAR SEMS WITH RANK-SPECIFIED GRAPHS
sections.text: We first consider linear SEMs on rank-specified graphs, with number of nodes d ∈ {100, 300}, rank r = d0.1de, and average degree k ∈ {2, 4, 6, 8}. The true rank is assumed to be known and is used as the rank parameter r̂ in NOTEARS-low-rank. For a better visualization, Figure 3 only reports the average SHDs, while the true positive rate, false discovery rate, and running time are 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, it 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.
  For more information regarding the role of sparsity, we include NOTEARS with an `1 penalty, named NOTEARS-L1. Here the `1 penalty weight is chosen from {0.01, 0.02, 0.05, 0.1, 0.2, 0.5}. Instead
  1Here we choose ICA-LiNGAM, 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 better performance than DirectLiNGAM for relatively dense graphs. Please find a detailed discussion and an empirical comparison in Appendix D.4.
  of relying on an additional validation dataset, we treat NOTEARS-L1 favorably by picking the lowest SHD obtained from different weights for each dataset. As seen from Figure 3a, NOTEARS-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 Zheng 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 a similar effect to the `1 penalty.
  Appendix D.1 studies graphs with higher ranks, where it is observed that the advantage of NOTEARSlow-rank over NOTEARS decreases when the rank of the underlying DAG increases. Nevertheless, NOTEARS-low-rank is still competitive when the true rank is dd/2e 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-low-rank performs reasonably well when the sample size is small and tends to have a better performance with a 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 k = 6, and power γ = 2.5. For this experiment, the minimum, maximum, and mean ranks of generated graphs are 14, 24, and 18.7, respectively. Here we choose the rank parameter r̂ from {20, 30, 40} for NOTEARS-low-rank. As seen from Figure 4, NOTEARSlow-rank with rank parameter r̂ = 20 performs the best, even though there are graphs with ranks greater than 20.
sections.heading: 5.3 SENSITIVITY OF RANK PARAMETERS AND VALIDATION
sections.text: 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 8, 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 1 and 3 to obtain a range of possible rank parameters, assuming that the corresponding structural priors are available. Within this range, we then select 7 evenly distributed rank parameters used with NOTEARS-low-rank to learn causal graphs. Finally, we evaluate each learned DAG using the validation dataset and choose the DAG with the best score as our estimate.
  As seen from Figure 5, NOTEARS-low-rank 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 a 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.
sections.heading: 5.4 NON-LINEAR SEMS
sections.text: For non-linear data models, we pick rank-specified graphs with 50 nodes, rank 5, and average degree k ∈ {2, 4, 6, 8}. To our knowledge, the selected benchmark methods CAM, NOTEARS-MLP, and
  GraN-DAG are state-of-the-art methods on this data model. As a 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, we use the same splitting ratio as 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 on the validation dataset is chosen as final estimate. Figure 6 (and Appendix D.6 with a more detailed result) shows that adding a nuclear norm can improve the performance of GraN-DAG across a large range of weights when the graph is relatively dense. For degree 8, the low rank version with validation achieves average SHD 77.4, while the SHDs of CAM, NOTEARS-MLP, and original GraN-DAG are 131.9, 119.4, and 109.4, respectively.
  5.5 REAL NETWORK
  We apply the proposed method to the arth150 gene network, 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 r̂ from {36, 40, 44} and also use validation to select the rank parameter. We apply NOTEARS-L1 where the `1 penalty weight is chosen from {0.05, 0.1, 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 t-test, we find that with significance level 0.1, the results
  obtained with r̂ = 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 of the graph is not very low.
sections.heading: 6 CONCLUDING REMARKS
sections.text: 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 is assumed. Theoretically, we provide an improved understanding of what kinds of graphs tend to be low rank and a 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. A future direction is to approximate a high rank DAG with a low rank one (possibly adding an additional DAG that is sparse). While there is a rich literature on low rank approximations of matrices and combining low-rankness with sparsity, it is non-trivial to us 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 a low rank DAG model implies any useful behavior in the data.
sections.heading: Appendix
sections.heading: A EXAMPLES AND DISCUSSIONS
sections.text: We provide more examples and discussions in this section.
  Minimum rank of the graph in Figure 1 We first show that the minimum rank of the DAG structure in Figure 1 is 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 is 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 8 visualizes the pathfinder and arth150 networks that are mentioned in Sections 4.3 and 5, respectively. Both networks can be found at http: //www.bnlearn.com/bnrepository. As one can see, these two networks contain hubs: the center note in the pathfinder network has a large number of children, while the arth150 network contains many ‘small’ hubs, each of which has 5 ∼ 10 children. We also notice that nearly all the hubs in the two networks have high out-degrees.
  Sparse DAGs and low rank DAGs A sparse DAG does not necessarily indicate a low rank DAG, and vice versa. For example, a directed linear graph with d vertices has only d − 1 edges, i.e. X1 → X2 → · · · → Xd, while the rank of its binary adjacency matrix is d − 1. According to Theorems 1 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, we can assign the edge weights so that the resulting graphs have low ranks. A simple example would be a fully connected directed balanced bipartite graph, as shown in Figure 9. The definition of bipartite graphs can be found in Appendix B.1. A bipartite graph is called balanced if its two parts contain the same number of vertices. The rank of a fully connected balanced bipartite graph with d vertices is 1 if all the edge weights are the same (e.g., the binary adjacency matrix), but the number of edges is d2/4.
  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
  graph is dense, then we need more vertices to cover all the edges. Thus, the size of the minimum head-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.
sections.heading: B PROOFS
sections.text: In this section, we present proofs for the theorems given in the main content.
sections.heading: B.1 PRELIMINARIES
sections.text: A bipartite graph is a graph whose vertex set V can be partitioned into two disjoint subsets V0 and V1, such that the vertices within each subset are not adjacent to one another. V0 and V1 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. A vertex cover of a graph is a subset of the vertex set where every edge in the graph has at least one endpoint in the subset. The size of a matching (vertex cover) is the number of edges (vertices) in the matching (vertex cover). A maximum matching of a graph is a matching of the largest possible size and a minimum vertex cover is a vertex cover of the smallest possible size. An important result about bipartite graphs is König’s theorem (Dénes, 1931), which states that the size of a minimum vertex cover is equal to the size of a maximum matching in a bipartite graph.
  Based on the heights of vertices in V, we can define a weak ordering among the vertices: Xi Xj if and only if l(Xi) > l(Xj), and Xi ∼ Xj if and only if l(Xi) = l(Xj). 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 Vs = {Xi; l(Xi) = s}, s = 0, 1, . . . , l(G), and let V−1 = ∅. We have: (1) for any given s ∈ {0, 1, . . . , l(G)} and two distinct vertices X1, X2 ∈ Vs, X1 and X2 are not adjacent, and (2) for any given s ∈ {1, 2, . . . , l(G)} and Xi ∈ Vs, there is at least one vertex in Vs−1 which is a child of Xi. If we denote the induced subgraph of G over Vs ∪Vs−1 by Gs,s−1, then Gs,s−1 is a bipartite graph with Vs and Vs−1 as parts, and singletons in Gs,s−1 (i.e., vertices that are not endpoints of any edge) only appear in Vs−1.
  For ease of presentation, we occasionally use index i to represent variable Xi in the following sections.
sections.heading: B.2 PROOF OF THEOREM 1
sections.text: Proof. Let G = (V,E). Consider an equivalence relation, denoted by ∼, among vertices in V defined as follows: for any Xi, Xj ∈ V, Xi ∼ Xj if and only if l(Xi) = l(Xj) and Xi and Xj are connected. Here, connected means that there is a path between Xi and Xj . Below we use C(Xi) to denote the equivalence class containing Xi. Next, we define a weak ordering π on V/ ∼, i.e., the equivalence classes induced by ∼, by letting C(Xi) π C(Xj) if and only if l(Xi) ≥ l(Xj). Then, we extend π to a total ordering ρ on V/ ∼. The ordering ρ also induces a weak ordering (denoted by ρ̄) on V: Xi ρ̄ Xj if and only if C(Xi) ρ C(Xj). Finally, we extend ρ̄ to a total ordering γ on V. It can be verified that γ is a topological ordering of G, that is, if we relabel the vertices according to γ, then Xi ∈ pa(Xj ,G) if and only if i > j and Xi and Xj are adjacent, and the adjacency matrix of G becomes lower triangular. Assume that the vertices of G are relabeled according to γ 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 V0 = {1, 2, . . . , k1 − 1} for some k1 ≥ 2. Then the k1-th row of A, denoted by A(k1, ·), is the first non-zero row vector of A. Letting S = {A(k1, ·)}, then S contains a subset of linearly independent vector(s) of the first k1 rows of A. Suppose that we have visited the first m rows of A and S = {A(k1, ·), A(k2, ·), . . . , A(kt, ·)} contains a subset of linearly independent vector(s) of the first m rows ofA, where k1 ≤ m < d. IfXm+1 Xkt , then we addA(m+1, ·) 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 l(Xm+1) > l(Xkt) ≥ · · · ≥ l(Xk1), then there is an index i ∈ Vl(Xm+1)−1 such that A(m + 1, i) 6= 0, by the definition of height. Since l(Xm+1) > l(Xkt), we have
  l(Xkt) ≤ l(Xm+1)− 1 and thus A(kj , i) = 0 for all j = 1, 2, . . . , t. Therefore, A(m+ 1, ·) cannot be linearly represented by {A(kj , ·); j = 1, 2, . . . , t} and the vectors in S are linearly independent. On the other hand, if l(Xm+1) = l(Xkt), then the definition of the equivalence relation ∼ implies that Xm+1 and Xkt are disconnected, which means that Xm+1 and Xkt do not share a common child in Vl(Xm+1)−1. Consequently, there is an index i ∈ Vl(Xm+1)−1 such that A(m + 1, i) 6= 0 but A(kt, i) = 0. Similarly, we can show that A(kj , i) = 0 for all j = 1, 2, . . . , t. Thus, the vectors in S are still linearly independent.
  After visiting all the rows in A, the number of vectors in S is equal to ∑l(G) s=1 |C(Gs,s−1)| based on the definition of ∼. The second inequality can be shown by noting that C(Gs,s−1) has at least one elements. The proof is complete.
sections.heading: B.3 PROOF OF THEOREM 2
sections.text: Proof. Denote the directed graph by G = (V,E). Edmonds (1967, Theorem 1) showed that max{rank(W );W ∈ WA} is equal to the maximum number of nonzero entries of A, no two of which lie in a common row or 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′0 ∪V′1, where V′0 = V × {0} = {(Xi, 0);Xi ∈ V} and V′1 = V × {1} = {(Xi, 1);Xi ∈ V}. Now define a bipartite graph B = (V′ ,E′) where E′ = {(Xi, 0) → (Xj , 1); (Xi, Xj) ∈ E}. Denote byM a set of nonzero entries of A so that no two entries lie in the same row or column. Notice thatM can be viewed as an edge set and no two edges inM share a common endpoint. Thus,M is a matching of B. Conversely, it can be shown by similar arguments that any matching of B corresponds to a set of nonzero entries of A, no two of which lie in a common row or column. Therefore, max{rank(W ),W ∈ WA} equals the size of the maximum matching of B, and further the size of the minimum vertex cover of B according to König’s theorem. Note that any vertex cover of B can 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 V′0 and of the vertex cover containing all variables in V ′ 1, respectively. Thus, max{rank(W ),W ∈ WA} is equal to the size of the minimum head-tail vertex cover.
sections.heading: B.4 PROOF OF THEOREM 3
sections.text: Proof. We start with the first inequality in Equation (4). Let h1, . . . , hp denote the heights where |Vs| < |ch(Vs)|, and t1, . . . , tq the height where |Vs| > |ch(Vs)|. Let H = ∪pi=1Vhi and T = ∪qi=1Vti . It is straightforward to see that (H,T) is a head-tail vertex cover. Thus, Equation (4) holds according to Theorem 2. The second inequality can be shown similarly and its proof is omitted. For the third inequality, let m = argmax{|Vs| : 0 ≤ s ≤ l(G)}, and define H = ∪i>mVi and T = ∪i<mVi. Then (H,T) is also a head-tail vertex cover and the third inequality follows from Theorem 2, too.
sections.heading: B.5 PROOF OF THEOREM 4
sections.text: 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 Xi −Xj is treated as two directed edges Xi → Xj and Xj → Xi.
  C IMPLEMENTATION DETAILS
  In this section, we present an algorithm to generate a random DAG with a given rank, a low rank version of NOTEARS and GraN-DAG, and also a description of our experimental settings.
sections.heading: C.1 GENERATING RANDOM DAGS
sections.text: 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 1. In particular, we aim to generate a random DAG with d nodes, average degree k, and rank r. The first part of Algorithm 1 after initialization is to sample a number N , representing the total number of edges, from a binomial distribution B(d(d − 1)/2, p)
  Algorithm 1 Generating random DAGs Require: Number of nodes d, average degree k, and rank r. Ensure: A randomly sampled DAG with the number of nodes d, average degree k, and rank r.
  1: Set M = empty graph, Mp = ∅, and R = {(i, j); i < j, i, j = 1, 2, ..., d}. 2: Set p = k/(d− 1). 3: Sample a numberN ∼ B(d(d−1)/2, p), whereB(n, p) is a binomial distribution with parameters n and p. 4: if N < r then 5: return FAIL 6: end if 7: Sample r indices from 1, . . . , d− 1 and store them in Mp in descending order. 8: for each i in Mp do 9: Sample an index j from i+ 1 to d. 10: Add edge (i, j) to M and remove (i, j) from R. 11: end for 12: while R 6= ∅ 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 18: if |M | < N then 19: return FAIL 20: end if 21: return M
  where p = k/(d− 1). If N < r, Algorithm 1 would return FAIL since a graph with N < r edges could never have rank r. Otherwise, Algorithm 1 samples an initial graph with r edges and rank r, by choosing r edges such that no two of them share the same head points or the same tail points, i.e., each row and each column of the corresponding adjacency matrix have at most one non-zero entry. Then, Algorithm 1 sequentially samples an 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 is because if a graph G is a super-graph of another graphH, then the size of the minimum head-tail cover of G is no less than that ofH. We repeat the above sampling procedure until there is no edge in R or the number of edges in the resulting graph reaches N . If the latter happens, the algorithm will return the generated graph; otherwise, it will return FAIL.
  The theoretic basis of Algorithm 1 is Theorem 2. Note that the algorithm may not return a valid graph if the desired number N of edges cannot be 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 1 to fail to return a desired graph.
sections.heading: C.2 OPTIMIZATION
sections.text: For this part, we consider a dataset consisting of n i.i.d. observations from P (X) and consequently the expectations in Problems (1) and (2) are replaced by empirical means. Denote the design matrix by X ∈ Rn×d, where each row of X corresponds to an observation and each column represents a variable. Here we use NOTEARS (Zheng et al., 2018) and Gran-DAG (Lachapelle et al., 2020) from each class of methods as examples and will describe their low rank versions in the following. Other gradient-based methods and their optimization procedures can be similarly modified to incorporate the low rank assumption.
  Algorithm 2 Optimization procedure for NOTEARS-low-rank Require: Design matrix X, starting point (U0, V0, α0), rate c ∈ (0, 1), tolerance > 0, and threshold
  w > 0. Ensure: Locally optimal parameter W ∗.
  1: for t = 1, 2, . . . do 2: (Solve primal) Ut+1, Vt+1 ← arg minU,V Lρ(U, V, αt) with ρ such that g(Ut+1V Tt+1) <
  cg(UtV T t ).
  3: (Dual ascent) αt+1 ← αt + ρg(Ut+1V Tt+1). 4: if g(Ut+1V Tt+1) < then 5: Set U∗ = Ut+1 and V ∗ = Vt+1. 6: break 7: end if 8: end for 9: (Thresholding) Set W ∗ = U∗V ∗T ◦ 1(|U∗V ∗T | > w).
  10: return W ∗
sections.heading: C.2.1 NOTEARS WITH LOW RANK ASSUMPTION
sections.text: Following Section 3, the optimization problem in our work can be written as
  min W
  1
  2n
  ∥∥X−XUV T∥∥2 F , subject to trace ( eUV T ◦UV T ) − d = 0, (5)
  where U, V ∈ Rd×r̂ and ◦ is the point-wise product. The constraint in Problem (5) holds if and only if UV T is a weighted adjacency matrix of a DAG. This problem can then be solved by standard numeric optimization methods such as the augmented Lagrangian method (Bertsekas, 1999). In particular, the augmented Lagrangian is given by
  Lρ(U, V, α) = 1
  2n
  ∥∥X−XUV T∥∥2 F + αg(UV T ) + ρ 2 |g(UV T )|2,
  where g(UV T ) := trace ( eUV T ◦UV T ) − d, α is the Lagrange multiplier, and ρ > 0 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 `1 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 Zheng 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 a sparsity assumption may be harmful, as shown also by Zheng et al. (2018); (4) the `1 penalty term requires a tuning parameter, which itself is not easy to choose.
  Zheng et al. (2018) used L-BFGS to solve the unconstrained subproblem in Step 2. We alternatively use the Newton conjugate gradient method that is written in C. Empirically, these two optimizers behave similarly in terms of the estimate performance, while the latter can run much faster thanks to its C implementation. The DAG constraint may not be satisfied exactly using iterative numeric methods, so it is a common practice to pick a small tolerance, followed by a thresholding procedure on the estimated entries to obtain exact DAGs. In our implementation, we choose U0 and V0 to be the first r̂ columns of the d× d identity matrices. Other parameter choices are: α0 = 0, c = 0.25, = 10−6, and w = 0.3, similar to those used in related methods on the same datasets (e.g., Zheng et al. (2018); Yu et al. (2019); Zhu et al. (2020)). The chosen threshold w = 0.3 works well in 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 a DAG.
  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 w = 0.3) to the refitted weighted adjacency matrix. Both the Newton conjugate gradient optimizer and the pruning technique are 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 D.3 for an empirical comparison.
sections.heading: C.2.2 GRAN-DAG WITH LOW RANK ASSUMPTION
sections.text: We next consider a low rank version of GraN-DAG. The optimization problem can be written as
  min θ − 1 n n∑ l=1 d∑ i=1 log p ( X (l) i | pa(Xi,W (θ)) (l); θ ) + λ‖W (θ)‖∗
  subject to trace ( eW (θ) ) − d = 0,
  (6)
  where X(l)i is the l-th sample of variable Xi and pa(Xi,W (θ)) (l) means the l-th sample of Xi’s parents indicated by the adjacency matrix W (θ). Here, θ denotes the parameters of neural networks and W (θ) with non-negative entries is obtained from the neural network path products.
  Problem (6) 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 is approximately solved using first-order methods; (2) the thresholding at Step 9 is replaced by a variable selection method proposed by Bühlmann et al. (2014). The same variable selection or pruning method is adopted by two other benchmark methods CAM and NOTEARS-MLP in our experiment. Please refer to Lachapelle et al. (2020) and Bühlmann et al. (2014) for further details.
sections.heading: C.3 EXPERIMENT SETUP
sections.text: In our experiments, we consider three data models: linear Gaussian SEMs, linear non-Gaussian SEMs (linear exponential SEMs), and non-linear SEMs (Gaussian processes). Given a randomly generated DAG G, the associated SEM is generated as follows:
  Linear Gaussian A linear Gaussian SEM is given by Xi = ∑
  Xj∈pa(Xi,G)
  W (j, i)Xj + i, i = 1, 2, . . . , d, (7)
  where pa(Xi,G) denotes Xi’s parents in G and i’s are jointly independent standard Gaussian noises. In our experiments, the weights W (i, j)’s are uniformly sampled from [−2,−0.5] ∪ [0.5, 2].
  Linear Exponential A linear exponential SEM is also generated according to Equation (7), where i’s are replaced by jointly independent Exp(1) random variables. The weightsW (i, j)’s are sampled from [−2,−0.5] ∪ [0.5, 2] uniformly, too.
  Gaussian Processes We consider the following additive noise model:
  Xi = fi(pa(Xi,G)) + i, i = 1, 2, . . . , d, (8)
  where i’s 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 are summarized 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.
sections.heading: C.4 BENCHMARK METHODS
sections.text: Existing causal structure learning methods used in our experiments all have available implementations, as listed below:
  • GES and PC: an implementation of both methods is available through the py-causal package at https://github.com/bd2kccd/py-causal. We note that, the implementation of py-causal 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 is available in the bnlearn package
  at https://CRAN.R-project.org/package=bnlearn.
  • CAM (Peters et al., 2014): its codes are available through the CRAN R package repository at https://cran.r-project.org/web/packages/CAM.
  • NOTEARS (Zheng et al., 2018) and NOTEARS-MLP (Zheng et al., 2020): codes are available at the first author’s github repository https://github.com/xunzheng/ notears.
  • GraN-DAG (Lachapelle et al., 2020): an implementation is available at the first author’s github repository https://github.com/kurowasan/GraN-DAG. Note that for graphs of 50 nodes or more, GraN-DAG performs a 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://github.com/fishmoon1234/DAG-GNN.
  • ICA-LiNGAM (Shimizu et al., 2006): an implementation is available at https://sites. google.com/site/sshimizu06/lingam.
  In the experiments, we mostly use default hyperparameters unless otherwise stated.
sections.heading: D ADDITIONAL EXPERIMENTAL RESULTS
sections.heading: D.1 LINEAR SEMS WITH HIGHER RANKS
sections.text: This experiment considers graphs of higher ranks. We use rank-specified random graphs with d = 100 nodes and rank r ∈ {30, 35, 40, 45, 50} on linear Gaussian SEMs. The results are shown in Figures 10a and 10b with degrees 2 and 8, respectively. We observe that when the rank of the underlying graph becomes higher, the advantage of NOTEARS-low-rank over NOTEARS decreases. Nonetheless, NOTEARS-low-rank with rank r = 50 is still comparable to NOTEARS, and has a lower average SHD after removing outlier SHDs using the interquartile range rule.
sections.heading: D.2 NOTEARS-LOW-RANK WITH DIFFERENT SAMPLE SIZES
sections.text: We next empirically study the consistency of NOTEARSlow-rank. Again, we use rank-specified random graphs (sampled according to Algorithm 1) with d = 100 nodes, degree k = 8, rank r = 10, and linear Gaussian SEMs. We also assume that the true rank is known. We fix the rank parameter r̂ = 10 and use different sample sizes ranging from 200 to 5, 000. From Figure 11, NOTEARSlow-rank performs reasonably well when the sample size is small and tends to have a better performance with a larger number of samples.
sections.heading: D.3 FURTHER PRUNING
sections.text: We compare the empirical results before and after applying the additional pruning technique described in Appendix C.2. The graphs are rank-specified with d ∈ {100, 300} nodes, rank r = d0.1de, and degree k ∈ {2, 4, 6, 8}. We again use linear Gaussian data model with equal noise variances to generate the datasets. The average SHDs are reported in Figure 12. We see that applying an additional pruning step indeed improves the final performance of both NOTEARS and NOTEARS-low-rank, especially on relatively large and dense graphs.
sections.heading: D.4 AN EMPIRICAL COMPARISON BETWEEN ICA-LINGAM AND DIRECTLINGAM
sections.text: To our best knowledge, there are two Python implementations of ICA-LiNGAM (Shimizu et al., 2006) released by the authors, available at https://sites.google.com/site/sshimizu06/ lingam and https://github.com/cdt15/lingam, respectively, where the latter is a Python package containing several LiNGAM related methods. In the following, we use ICALiNGAM-pre and ICA-LiNGAM-cdt to denote these two implementations, respectively. For DirectLiNGAM (Shimizu et al., 2011), we only find a Python implementation available at the previously mentioned Python package containing ICA-LiNGAM-cdt.
  Here we run DirectLiNGAM, ICA-LiNGAM-cdt, and ICA-LiNGAM-pre on linear exponential data models with 100-node and rank-10 graphs. The mean SHDs are reported below in Table 1. Based on this experimental result as well as our past experience, DirectLiNGAM usually has a (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-pre in the main paper.
sections.heading: D.5 DETAILED EMPIRICAL RESULTS FOR EXPERIMENT 1 WITH LINEAR GAUSSIAN SEMS
sections.text: Table 2 reports detailed results including true positive rates (TPRs), false discovery rates (FDRs), structural Hamming distances (SHDs), and running time on rank-specified graphs with linear Gaussian data model. Here the true rank is assumed to be known and is used as the rank parameter in NOTEARSlow-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 a similar performance to that of GES. Therefore, we only include the results of GES for comparison. We treat GES favorably by regarding undirected edges as true positives if the true graph has a directed edge in place of the undirected ones.
sections.heading: D.6 DETAILED RESULTS FOR EXPERIMENT 4 WITH NON-LINEAR SEMS
sections.text: 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.
title: ON LOW RANK DIRECTED ACYCLIC GRAPHS