Sequentially learning the topological ordering of directed acyclic graphs with likelihood ratio scores
Authors: Gabriel Ruiz
OSCAR HERNAN MADRID PADILLA
Qing Zhou
Venue: Accepted by TMLR ()

Abstract: Causal discovery, the learning of causality in a data mining scenario, has been of strong scientific and theoretical interest as a starting point to identify "what causes what?'' Contingent on assumptions and a proper learning algorithm, it is sometimes possible to identify and accurately estimate an underlying directed acyclic graph (DAG), as opposed to a Markov equivalence class of graphs that gives ambiguity of causal directions. The focus of this paper is in highlighting the identifiability and estimation of DAGs through a sequential sorting procedure that orders variables one at a time, starting at root nodes, followed by children of the root nodes, and so on until completion. We demonstrate a novel application of this general sequential approach to estimate the topological ordering of the DAG corresponding to a linear structural equation model with a non-Gaussian error distribution family. At each step of the procedure, only simple likelihood ratio scores are calculated on regression residuals to decide the next node to append to the current partial ordering. The computational complexity of our algorithm on a $p$-node problem is $\mathcal{O}(pd)$, where $d$ is the maximum neighborhood size. Under mild assumptions, the population version of our procedure provably identifies a true ordering of the underlying DAG. We provide extensive numerical evidence to demonstrate that this sequential procedure scales to possibly thousands of nodes and works well for high-dimensional data. We accompany these numerical experiments with an application to a single-cell gene expression dataset. Our $\texttt{R}$ package with examples and installation instructions can be found at https://gabriel-ruiz.github.io/scorelingam/.

Published in Transactions on Machine Learning Research (12/2022)
Sequentially learning the topological ordering of directed acyclic graphs with likelihood ratio scores
Gabriel Ruiz*
rulag Sucla edu
Deportment of Statistics University of California, Los Angeles
Oscar Hernan Madrid Padilla
oscar.marrddOssaa.. ucla. edu
Department of Statistics
University of Colifornia, Los Angeles
Qing Zhout
show@stat. ucla edu
Department of Statistics
University of California, Los Angeles
Reviewed on OpenReview: https: //tps://openaees opw.revee...e.eeeeeo net/ forum? id- 4pCyIGI
Abstract
Causal discovery, the learning causality in a data mining scenario, has been st.oong scien- tific and theoretical interest a starting point identify "what causes what?" Contingent assumptions and proper learning algorithm, it is sometimes possible to identify and accurately estimate an underlying directed acyclic graph (DAG), as opposed to Markov equivalence class of graphs that gives ambiguity of causal directions. The focus of this paper is in highlighting the identifiability and estimation of DAGs through a sequential sorting procedure that orders variables one at a time, starting at root nodes, followed by children of the root nodes, and so on until completion. We demonst a novel application of this general sequential approach estimate the topological ordering of the DAG corresponding to a linear structural equation model with non-Gaussian error distribution family. At each step of the procedure, only simple likelihood ratio scores are calculated regression residuals to decide the next node to append to the current partial ordering. The computational complexity of our algorithm on a p-node problem is O(pd), where d is the maximum neighborhood size. Under mild assumptions, the population version our procedure provably identifies a true ordering of the underlying DAG. We provide extensive numerical evidence to demonstrate that this sequential procodure scales possibly thousands of nodes and works well for high-dimensiooal data. We accompany these numerical experiments with an application single-cell gene expression dataset. Our R package with examples and installation instructions can be found at htttp:://gabriel-riiz wtw.eeelllriiiigghghuuu wtallbiooccoooroo io/scorelingam/
1 Introduction
With observational data alone, causal inference using an accurate directed acyelic graph (DAG) has been shown to provide results that up to par with the quintessential randomized controlled experiment (Pearl
2009 Malinsky et al., 2019). However, it is difficult to imagine that this approach, with its reliance strong domain knowledge about the system of variables al hand, can be applied to cases with large numbers of variables and little background on how they are all related. For example, specifying a DAG may be difficult in bioinformatics and fields of science whore "big data" was previously unavailable and we are now trying to the DAG structure for Bayesian networks
get a grasp of it. On the other hand, causal
"Supported NSF grant DGE-1650604. Supported NSF grant 1952999..

---

Published in Transactions on Machine Learning Research (12/2022)
2 X
Z X
Figure 1: The original DAG (left) and its corresponding CPDAG (right): obtained by keeping the orient a ation of edges corresponding the v-structures - X + and Z + Y U, and removing the orientation from all other edges. Note the ambiguity about the causal dircction X - Y vs. x+ Y in the CPDAG.
from scratch- its own pitfalls, such as a super-exponentially increasing number of networks in the search space as the number of nodes grows and the fact that it is possible for multiple directed acyclic graphs and the Bayesian networks they encode to map to the same joint distribution-a phenomenon called Markov
equivalence (Frydenberg] 1990 Verma & Pearl 1990 Peters al. 2017 Verma & Pearl 2022 J. Nonetheless,
a large effort has been devoted to the structure learning of DAGs as a preliminary data mining step in the
scientific pipeline (Stekhoven al., 2012
Nandy et al. 2017 Malinsky & &Spirtes 2017: Henckel al. 2022:
Lee et al. 2022).
The present paper focuses on a sequential learning algorithm, which we show can be quite scaleable, along with novel identification theory for a specific application. The main task of the approach we advocate for is to sequentially estimate a topological ordering of the DAG, a permutation of node labels such that every parent must precode its children. To help with scalability in practice, we also make use of a priori known neighborhood sets, such as Markov blanket of a node. In order to demonstrate the thcoretical promise of this procedure, we discuss existing identification results that make linear structural equation model (SEM) first studied in (Shimizu et all (2006). The novelty of our application of the sequential sorting procedure to this SEM compared to the of for it is the scalability of our procedure to a large number of nodes in the underlying graphical model.
of it. We also provide new theory for
Representative methods for causal discovery under the assumption of no unobserved confounding (causal sufficiency) include the Peter-Clart (PC) algorithm (Spirtes & Glymour, 1991) and Greedy Equivalence Search (GES) (Chickering, 2002). ). The PC algorithm is a constraint- based method due to its use of conditional independence queries, while GES is considered score-based method for the objective function it sceks to optimize across the space of graphical models. Without additional structural assumptions, the best these methods can generally do in the limit of sample size (n + oc) is to obtain Markov equivalence class (MEC) of DAGs, visualized typically by a single Completed Paritially Directed Acyclic Graph (CPDAG) as in Figure ! Each DAG in the MEC, obtainable by orienting undirected edges in the CPDAG without introducing a cyclic path nor a "v-structure, encodes the same set of d-separalion relations that imply marginal and conditional independence relations between triplets of variable subsets in their underlying joint distribution (Spirtes& & Glymour 1991).
When additional assumptions are justified, such strict non-linearity of structural equations, or non- Gaussianity of noise terms in a linear structural equation model, a unique DAG can be identified (Bihlmann et al. 2014 Shimizu et al., 2006). When we are not willing to make the assumption of causal sufficiency, the Fast Causal Inference (FCI) algorithm provides an alternative at the cost of a potentially less precise, though possibly more accurate, graphical model compared to a DAG or CPDAG (Spirtes et al., 2000 Malinsky & Spirtes 2017). Beyond the case of independent and identically distributed (iid) copies from distribution
that our DAG of interest satisfies the Markov property with respect to, Glymour et al. (2019) and Peters et al. (2117) provide reviews on the trade-offs of different algorithms and what can and cannot be done when there is additional structure, such as the case that the system of variables varies in time. In the context of Earth system sciences, Runge el al. (2019) review causal discovery methods. Structure learning has also been explored for its possibility to explain the black-bor nature of state-of-the-ar deep learning architectures (Sani ct al. 2020). Morcover, Zhenng et al. (2018) and its extension to Zhong et al. (2020) provide an approach to optimize a non-convex score function in DAG space by using a smooth characterization of an adjacency
matrix's acyclicity constraint.
2

---

Published in Transactions on Machine Learning Research (12/2022)
1.1 Review of relevant work
Specilic to our lask of learning topological ordering for an underlying acyelic graphical model, we now review some relevant work. Let us first formally define a so-called topological ordering of a DAG, the target parameter we seek to estimate. Let [m]] = (m....,m} for integer m 2 1. Moreover, let PAg C [p] denote the nodes in the DAG with directed edges going into node k.
Definition 1.1. topological ordering for a DAG G given by a permutation π + [p] such that every parent node precedes its child in the ordering:
jE PAk = π (i) < (
Importantly, we note that the discrete search space across p! permutation functions in search of one that satisfies Definition 1.1 can be quite cumbersome (Raskutti & Uhler 2018 Solus et al. 2021 p. Several heuristic score-based methods have been developed to cope with the search space, however, remains the case that
score-basod approaches for ordering search are NP hard in general (Chickoring 1996 Ye ct al.] 2021). Along these lines, recent work by Ye et al. (2021) provides one approach under the case of linear Bayesian network
with Gaussian noise. The similar non-parametric approach of Solus et al. (2021) and Wang et al. (2017) requires consistent conditional independence testing procedure to decide the presence or absence of an edge in the DAG corresponding to a given # in the search space. The empirical results of these approaches are all promising. However, these works do not provide guarantees on whether the search can terminate at a point well before querying all permutations or DAGs in order to achieve optimal (statistically consistent) results.
Complementary to these advances, we here study a simpler approach to estimate a permutation for which the search has pre-determined number of steps: Oppd) in the number of least squares residual updates, where d S p is the maximum neighborhood size of node. The objective of our work is to estimate permulation from observed data X E Rexp of sample size n. We denote this estimate as â. To do so, we will apply Algorithm [ for t = 1,2,.... P until all nodes are sorted. At step t, given an input partial ordering A = (f(1), #(t 1)), the algorithm selects a node ftt & Ar to append to the estimated ordering by maximizing an appropriate score Ar; of interest.
such
Algorithm 1: Continue a Topological Ordering
The partial ordering A, and data matrix X E g2xp Result: The continued partial ordering A.+1
Data:
*tt) + arg m. 8k
A+11 + A U {f(t))
There exist structure learning methods that use the general approach in Algorithm to sequentially construct topological ordering. These approaches motivate our present work and include the following. Peters et al. (2014) apply Algorithm under an assumption of strictly nonlincar structural equations with additive noisc. Meanwhile, Ghoshal Honorio (2017), Park (2020). p Park & Kim (2020), and Chen el al. (2019) apply this sequential sorting procodure under a bounded conditional variance assumption: a V[X,)XXaa S for each e [p] and some unknown positive constants a S b restricted by the signal a parent sends its child node. Notably, Park (2020) contains extended discussion on the case of node's possibly non-linear relation with its parents, while Gao et al. (2020) further explore the computational scalability for the sequential application of Algorithm [ to estimate non-linear structural equation models under this bounded conditional variance assumption. With respect to linear SEMs with non-Gaussian noise, applications of Algorithm include Shimizu et al. (2011). Hyvärinen & Smith (2013). and Wang & Drton (2019), ; while Zeng et al. 2020 construct the topological ordering reverse starting with child-less nodes. We believe there exists potential to scale up the estimation of each of these models. We focus on the linear SEM with non-Gaussian noise here.
1.2 Paper Contribution and Outline
Our application of Algorithm 1 is to a linear SEM with non- non-Gaussian noise under causal sufficiency, which is known as the linear non-Gaussian acyclic model (LIGGAM). In terms of theoretical guarantees for the
3

---

Published in Transactions on Machine Learning Research (12/2022)
estimation of LiNGAM, Shimizu al. (2006) and Shimizu et al. (2011) provide identifiabilly results for the respective LiNGAM learning procedures that is, with knowledge of the true distribution defined by the LiNGAM and an oracle for conditional independence queries in the case of latter. Meanwhile Wang & Drton
2019 provide formal statistical consistency results for their LiNGAM-learning procedure.
Although the above methods have very nice theoretical guarantees, their practical application is limited as they do not presently scale well large graphs, say with thousands of nodes, as confirmed in the numerical results in this paper. Therefore, we develop fast sequential learning method that can estimate large graphs in practice. At each step of this method, node is selected to append to partial ordering, steps, where is the number of nodes in the underlying graph, full ordering of all the nodes will be produced. Compared to the existing works on LiNGAM, the main contributions of our work are:
that after
Based on a specified error distribution, we define a novel likelihood ratio score which used at each step in our sequential algorithm. The evaluation the likelihood ratio only involves linear regression and residual calculation. There are no tuning parameters.
2. We prove that at the population-level, this sequential algorithm will identify true ordering of the underlying DAG under proper assumptions on the LINGAM.
3. Our sequential method is computationally tractable with computational complexity O(p²) for the number of updates used in the entire algorithm. Tf prior knowledge on the Markov blankets of the nodes is provided, the computational complexity can be further redueed to O(pd), where d is the maximum size of the Markov blankets. This is in sharp contrast. to traditional score-based approaches
for ordering search, which are NP hard in general CCickering 1996 Ye et al. 2021
4. The computational complexity of the finite version of our procedure is O(pdmn) in the number of operations needed. The term m < is the maximum number of ancestors node can have (see Algortihm Meanwhile, the procedure in Shimizu et al. (2011) is stated to be O(pp M2 + p² in the number of operations, where M is the maximal rank found by the low-rank decomposition used in their kernel-based independence measure of choice. Moreover, Wann & Drton] (2019) address computational complexity in passing by stating that "the computational complexity of [their sorting algorithm] is exponential in [the maximum number of parents node can have], so in practice [the maximum mumber of parents a node can have] must remain relatively small."
The rest of the paper is organized as follows. Tn the rest of this section, we formally introduce the linear SEM of interest. Next, in Section 2 we will introduce our approach: 2.2 discusses the conditions for this approach to work; $2.3 provides a formal identifiability result; and $2.4 provides the finite sample version of the algorithm. Section 3 presents simulation results for our procedure for small and large-sized Bayesian networks, along with an application to single-cell gene expression data. Finally, we conclude with summary of our findings and discussion of future work.
1.3 Review of LiNGAM
We follow closely here the definition of a LiNGAM given by Shimizu et al. (2066)
Definition 1.2. (Lincar Non-Gaussian Acyelic Model) For p 2 2, let G be a DAG on p nodes and E Rexp be the weighted adjacency matrix of G such that Byt # means j PAk, the parent set of node k. Let e= (eI...., such that eş ~ g(=0%) independently, with g(: 0) density of non-Gaussian distribution parameterized by 0 E Ro, and E(c] = 0 (j = 1, ......... p)-
We say X E Rº follows a LiNGAM with DAG G if
X
Σ ByeX, + ck.
==
(1)
JEPAx
Denote X's probability density function (pdf) as f(x
4

---

Published in Transactions on Machine Learning Research (12/2022)
The scalar form of the linear SEM in Equation 1 can be rewritten in vector form as X B" + Put M (Tp - B)=² matrix with ones on its diagonal. Let ANg denote the ancestor set of node k: ANg means there exists a direct path starting at node a and ending at node k, a + -.... k. Then we arrive at X M in vector form and Xx = c₁ + ΣjEAN, MaJC, in scalar form for all k E [p]. Noting that M serves as a mixing matrix for the independent components in 6, we may think of the estimation of this linear SEM as an instance of Independent Component Analysis (ICA) (Hyvärinen & Oja 2000). Shimizu et al. 2006 discuss the connection between LiNGAM and ICA.
2 Methodology and algorithm
In this scction, we introduce both the population-level and finite-sample versions of our sorting procedure. We also show that our choice of summary score A;; X) in Algorithm [ will lead to the identification of a topological ordering of the true DAG G used to define the lincar SEM of Definition some useful notation.
1.2
We first go over
2.1 Summary of Notation
For a positive integer m, we write [m]= {1,2...., .m.. For any set S, [S] will denote its cardinality: the number of unique elements it contains. |||=0 = 0 means S = 0, the set with no elements. For sets T S [m], S [r], and matrix A E Rm A.s E RM is the sub-matrix defined by indexing columns S in A, while Ar. E g is the sub-matrix given by indexing rows T in A. Similarly, ATs E R|Tix|s| is the sub-matrix indexing rows T and columns S of A. Similarly, for vector FT we will write uT E RIT to denote the subset of entries indexed by T. We will also sometimes write [æy:j E T) to denote er. When the indexing sets are singletons, S= [k} or T =[}, we will simply write Ak, Aj., Ajk, and Uk, respectively. For two sets S and T, we will make use of set operations such as their intersection: SnT= = fa: E S and E T}; their union: SUT = (a: a E or a E T); and their difference: SYT = (a: a E S and a E T). For sets Si...., Sk, we denote their intersection as n Sj and their union as U Sz-
We will refer to a directed acyelic graph by G = (V,E). G is collection of singular vertices, V = [p], and vertex pairs, & C V x V. Specifically, (j.k) E E mcans that j has a directed edge into & in G. PAe denotes the parent set of node contains all nodes k such that kj,k) E E. Let NE S V denote the set of nodes that are thought to be most predictive of node k, which we will refer to as node k's neighborhood. Tet NEt = No n Ac, which is the subset of the neighborhood set that has been ordered at step of our procedure (Algorithm []. For the cases where |iaal 2 1, we will make use of least squares residuals for calculating the score S(k, As). At the population-level, the residual is
V: j E PAg ouly if (j) E The set of all of node j's children, CH₃,
Xx Xx
if N = otherwise
Rst =
(2)
X N.. x
where Ber is the least-squares regression coefficient vector,
BKe= (B X Ne.
x2 12.1)
E "B
Also,
Tik arg max Ex~ [ogg g( RRGi)) "
while p(rk; σke) is the density for N oZ. v[Rad),, i.. the normal density that matches the mean and variance of Ras. The corresponding sample analogues of Rai, Bie, lhtt and oke, are discussed in 2
2.2 Assumptions
Our main assumptions are on the distributions of the independent entries of the error vector c. We consider restricting our class of densities {g(: (ssss for the noise terms in Definition 1.2 to scale-location family
5

---

Published in Transactions on Machine Learning Research (12/2022)
in which the 0& 0 are the scale parameters, such the Laplace family of distributions, the Logistic family of distributions ions, Scaled-t distribution family (same degrees of freedom). This is summarized in Assumption
2.1
Assumption 2.1. Let U g(::8u with &o > and E[U] = For each ào = .P. the density of the error ej satisfics
Go
g(e; 0n) (00e//0zz00)
02
d (0,/0,)U, an equality in distribution.
That is, E
Our next assumption for the linear SEM interest is on linear combinations the noise terms. This condition is related Lemma 2.6 in $ 2.3.1 a key result about how to characterize the regression residuals of Equation 2 linear combinations of "independent components."
Assumption 2.2. For any 1, ... and any E IRP with at least two non-zero entries, the linear combination ae does not follow the same distribution as Ey.
Notable disagreements with Assumption 2.2 are when the Ej are all Gaussian distributed (not the case for LiNGAM), or when the Ej are all Poisson- Poisson-distribute distribbteedd Notable agreements with Assumption 2.2 (and Assumption 2.1 are the cases where the ey are all Taplace-distributedd, all Togistie-distributed, or all Scaled-t. distributed (same degrecs of freodom). More generally, we may start with any reference distribution, then specify the distribution for error terms in Definition 1.2 as re-scalings of the original reference distribution. So long as the characteristic function [Forbes al., 2010) for the sum of any two independent random variables following distribution in this scale-location family does not have the same form the characteristic function of the original distribution of summand, then Assumption 2.1 and Assumption
2.2 are satislied.
To allow for quicker sorting procedure in practice, we may make use of an a priori known support set for the neighborhood of each node in the DAG. We consider these neighborhood sets to arise basod on domain knowledge, previons studies, or a pre-processing step such as with neighborhood lasso regression (Meinshausen & Bühhlann, 2006). We highlight this usage Assumption 2.3
Assumption 2.3. For node k, denote its neighborhood estimate as No. Assume for each k = 1,2,.... that:
ÑA MB, = PA UCHHU U PAJM(
JECHz
where MB is known as the Markov Blanket of node k: the set of its parents PAs, its children CH2, and its co-parents UJECH, MA))
Remark 2.4. When we consider the population-level version of our algorithm in this section (i.e. we have infinite n); we can take N [p] for each i so that Assumption 2.3 holds trivially. For the finite sample version of our procedure discussed in Section 2.4, we will make use of Ordinary Least Squares (OLS) lincar regressions which require the design matrix to be of full column rank. So p & n, we may also take N₂ = (p) [k} for each k. In the case that p > n or p R n, sparse neighborhood sets, [& << n, will reduce the number of covariates in the required OLS regressions.
2.3 Our Choice of Likelihood Ratio Score
In Algorithm I we will select the next node continue our constructed topological ordering as:
log g( RRRs vike) o( RRi oke))
*()) arg max Ex) kg
(3)
Here, Ex= denotes expectation with respect to X's true density, f(x in Definition 1.2
Take a "valid" node continne the ordering A, node whose parent nodes are all contained in the set A And take "invalid" node to continue the ordering At a node whose parent. nodes not all in the set.

---

Published in Transactions on Machine Learning Research (12/2022)
Ar- The log-likelihodd ratio in Equation can be thought of score that allows to distinguish between valid and invalid nodes by telling "how non-Gaussian" the residual Rat is. If the residual is explained by a Gaussian distribution well relative to the non- Gaussian distribution in the assumed family, then we expect the log-likelihood ratio to be smaller. Otherwise, if the Gaussian density is not a good fit relative to g(ka Ta), then we have stronger evidence to believe that node & is a valid node to continue the ordering.
In Theorem 2.5 we claim that using Equation leads to the identification of a valid topological ordering our main result. The proof of Theorem [.5 is in Appendix AAI1 Key to the proof, note that Equation g can
also be written equivalently using the difference of two KL-divergence (Cover & Thomas 2005 terms:
(()= arg max (DK (Zse(raa))lo((rii orx)) D (fea(Tia)) (ru ma))). kga,
Here, fra(r is the probability density function of the residual RRt, which in general can be different from the specified density g( Mkr). The first term in Equation 4 deals with the divergence of the Rg 's true distribution from a Gaussian distribution. This term will be largest when k a valid node to append to the current ordering. The second term in Equation I concerns the deviation between Ry 's true distribution and the specified distribution given by the probability density function g(rni It will be non-zero when node & is an invalid node to append to the current ordering. Section 2.3.1 elaborates this explanation by introducing the key lemmas for the proof of Theorem 2.5
Theorem 2.5. Let X E Rp follow a LiNGAM with DAG g. If Assumplions 2. and 2.3 hold, then
applying Algorithm 1 at all steps = 1.2,.....7 with the score
S(k, A,) f(x)(s) log g (RRaine
@(Tri oke)
will identify permulotion म (#(1).... #pp) that is topological ordering of g.
Theorem 2.5 suggests that the maximization at each iteration in which we apply Algorithm 1 can be done easily. This differs from maximizing a score over whole ordering which may also lead to identification of the true MEC, but is in general NP hard (not tractable). Relatedly, Section [.3.2 gives additional motivation for the choice of At) in Theorem 2.5 as one that allows us to greedily optimize the mean log-likelihood when the full ordering is only partially discovered.
2.3.1 Key Lemmas for the proof of Theorem 2.5
The formal proof of Theorem 2.5 in $ A.1 below is an inductive application of the following reasoning when applying Algorithm [ at given step t. Lemma [.6] suggests that invalid nodes' residuals, Rat, are a linear combination of two or more entries in the vector e, while for valid nodes l we have Re = Eg. Under Assumption 2.2 this means that the term DEL (aa((r))llg)((ee) m)a)) in Equation 4 will be zero only if node k is valid to continue the ordering at step t.
Lemma 2.6 (Characterizing nodes' residuals as linear combinations of independent components).
Assume that A, is correct so far in the sense that for each a E Ae, we have PAg S At. Also assume Assumption 2.3 holds. We have that:
If S, is valid node to continue the ordering, i.e. PAd Aty then:
Rkc Xx S1.
X No. 6k-
Otherwise, if is not a valid node, then Rus a linear combination of more than one independent component in e.
The natural follow up question is what the behavior is for the term DEL (fra(Fa)) (Ra)) in Equation 4 when h is valid vs. invalid to continue the ordering. Lemma 2.7 provides this insight: for valid nodes to continue an ordering, this term's value is no less than the same term's value for invalid nodes. Given this and that the term Dxz (fuc(rx)) gg(rresk Tke)) Equation 4 will be zero only node & is valid continue the ordering at step t, Equation (equivalently Equation 3 helps us select a valid node continue the topological ordering at step L.
7

---

Published in Transactions on Machine Learning Research (12/2022)
Lemma 2.7 (KL Divergence from Gaussianity for valid and invalid nodes' residuals).
Let X E IR be LANGAM from Definition 1.2 that satisfies Assumplions 2.1 and 2.3 Assume that A is correct in the sense that PAa S Ae for all a E Ar. Let k E Ap]\ be an invalid node to continue the ordering in the sense that there crists j E PAk such that j E [p] As. And let & E [p] be a valid node to continue the ordering in the sense that PA()) S Ar- Then the least squares residual Re ~ See(re)) no closer to Gaussian than Rgt ~ [ka(TKa the sense that:
DEL (Fa((ra)lllori S DRL (falra)llou(ta))
where Tke and de are the respective densities for
RRi (( v(Rk)) and Ra ((B[Ra)) V(Rel).
In light of Lemma 2.6 Lemma 2.7 makes sense under Central Limit Theorem-like argument: a sum of two or more random variables is closer to Gaussian than each summand alone. Of particular note, key result that helps show why Lemma 2.7 holds Theorem 1 17.8.1 of Cover & Thomas (2005). a restatement of the enlropy-power inequality. This restalement says that the differential entropy for a sum any two independent. random variables, U and V, is no less than the differential entropy for the sum of two strategically defined Gaussian random variables, each having the same differential entropy as U and V (rather than the samne first. two moments), respectively.
2.3.2 Greedy Choice of a Factor to Optimize the Joint Likelihood Function
We now give further motivation to the choice of the mean log-likelihodd ratio score in Equation 3 Let vector
X ~ f(x), where f(x) is the density induced by LiNGAM in Definition 1.2 log-likelihood as function of the permutation m:
Consider X expected
[(T) = > Ex= logg (X, - [B"]" x;
(6)
where 07 is the corresponding scale parameter according to Assumption 2.1 Here, B is the acyelic weighted adjacency matrix that arises from a population-level least squares objective such that the x(j)-th column is given by:
arg min Ex=s(s) [(X*t) - - x
y =0 Vk s.t (A) 2
That is, the column BT =T) is comprised of the least squares coefficients when linearly regressing Xx()) onto its predecessors, if any, in the ordering given by n. Now let of be the density for the Gaussian distribution having the same first two moments as:
: -
Define
P [(T) = Σ Ex-s(a) [logo] (R7) and = Ex-s(a) (oo:/(((X0 V[X])]
Here, [xx;E[[]] V[x]) denotes the density for variaat Gaussian distribution with the same first and second order moments X. Due to the relation between B and the generalized Cholesky factorization of V[X]. Ye et al. (222] shows that we actually have the equality:
[())= = for all π.
(7)
The term x is notably a constant that does not depend on स. Thus, maximizing Equation 6 with respect to π is the same as maximizing the expected logelikcl ratio given by:
p
g (15:05) [(T)) K.
(C- ()( = Σ Ex Ex-
log
(8)
o7 (R.)
8

---

Published in Transactions on Machine Learning Research (12/2022)
With all this in mind, can think of our choice of node to append to the ordering A at step t as greedily choosing the largest summand,
E [x=)(x) log
( (R(()))
to add to the known log-likelihood ratio at step t:
t= log
(c (c- ())
E E Ex-s(a)
of af ( (75)
2 2
That is, our sequential application of Algorithm is attempting to greedily maximize Equation g one summand at lime.
2.4 Finite Sample Sorting Procedure
Assume that we have data matrix X E R"xp such that X,., the i-th row, is iid across = 1,2, .., from
distribution defined by LiNGAM satisfying Assumptions 2.1 and 2.2 Also let Assumption 2.3 hold.
Analogous to the population version of our sorting procedure in Section 2.3
consider:
Bxt
X X,N (x Net
X Nue
X
R
Nulx1
which exists so long as 1 1 [ikel n and X Ni is of full column rank almost surely. Further, we define RRte [xx
X. X
Nel = N
Rst
Bke
Nuc
the voctor of residuals which we will use to estimate the pertinent scale parameter of Equation B denoted as fike and Ōkt; respectively. Explicitly, we select the next node continue an ordering using the empirical analogue of the mean log- likkeiihood ratio in Equation 3
1
g( Mi fnk)) Σ log
*()) arg max
(9)
Φf R (R Ōzc)
where Re the i-th entry of the vector Ree, while o% = + REu 12 and the == arg max,, Σí log g( Ri m)- For example, if mu is the scale parameter for Laplace distribution, it can be seen that. fjat 1 R In this case, Equation is equivalent to
*(t) = arg max log Ōkt
arg max Rkellz kg.Ac [
(10)
Rg.Ar
iik
Remark 2.8. The Laplace update of Equation 10 excmplifies how simple the maximization of our likelihood ratio score is. After the regression of each unsorted node XN, & At, onto Net, we only need to compare the ratio between the two norms of the residual vector Rit across unsorted nodes to find (t).
Algorithm 2 shows the pseudo-code for the sorting procedure we use in practice with strategic update of regression residuals using partial regression that greatly reduces the computation cost. The computational complexity for Algorithm 2 is (pdmm) in the number of scalar to scalar multiplication, division, addition, or difference operations. Here, m < is the maximum number of ancestors a node can have, which is a bound on the iterations for the inner-mos most for loop in Algorithm
We have also provided the details on the estimation of the scale parameters for Logistic and Scaled-t distributions in Appendix C
9

---

Published in Transactions on Machine Learning Research (12/2022)
Algorithm 2: The sorting procedure in practice
Data: X E Rexp (slandardized), {N
Result: #(1), (2,....., (p)
# initialize mixing matrix
M + Ipxp
# initialize residual matrix
R +x
# initialize scores
8* S(k; R), k = 1,2, ,.
# sort the nodes
for 1= 1.,,,....,PP
do
((t) f- A. 8%
# update residuals for neighbors of selected node.
for & E No \e do
# update residuals with partial regression.
for a E {j: Matty M, 0} do -
MR
R. f R.
M R..
end
# update score
8 S(k, Ac: R)
end
end
3 Empirical Results
We now present empirical results for synthetic and real-data scenarios. For the boxplots of empirically observed distributions in Figures 2 3 and 7 the three horizontal lines on each box indicate the first, second, and third quartile, respectively, while the vertical lines indicate the range of the empirical distribution.
3.1 Simulations on Small Networks
We now present simulation results for networks that are on the smaller end with 35 < p S 223. These architectures are downloaded from the bnlearn. com Bayesian network repository. We compared our sorting procodure to other LiNGAM learning procodures. Due to their readily available code, the algorithms of interest are "DirectLiGGAM" (Shimizu et al., 2011). "HighDimLingam" (Wang Drton 2019). and "ScoreLiNGAM" (our procedure). For each simulation setting, we conduct replicates.
For each choice of G underlying LiNGAM, our synthetic data generation schema f.idd Uniform[-0.9, -0.4 U [0.4, 0.9] for cach kj, such that j E PAt. and otherwise set generated Bj
follows. We
Bi = O. We generated 0z i.idd ~ Uniform 00.4 0.7] across 1 S k S where 0₂ is the scale parameter for the error distributions as in Assumption 2.1 This specific range of simulation parameters for the weighted adjacency matrix entries and scale parameters chosen so that the marginal variance of X, does not explode when j is node with many ancestors. Finally, we varied sample size as n 0.5p. p. 2p, 10p, 50p. Note that 0.5p and n p represent the high-dinensional setting (p 2 n).
Next, dataset X gaxp of iid samples is drawn from the distribution given by the LINGAM parameterized (x) across & E [p]. Morcover, we varied the family of the densities in Assumption 2.1 to be the Laplace, the Logistic, or the Scaled-t distribution (10 degrees of freedom) scale-location families. Finally, ScoreLiNGAM and HighDimLiNGAM were run with knowledge of the true Markov blanket for cach node, while DirectLiNGAM was not as it does not have this option. Afterward, the data matrix X was standardized so that each column has sample standard deviation equal to and sample mean equal to 0.
by (B.0..... 8p) and having errors
10

---

Published in Transactions on Machine Learning Research (12/2022)
mildow (p=35)
alarm (p=37)
barley (p=48)
hailfinder (p=56)
hopar2 (p=70)
winStpts (p=76)
andes
(p=223)
0.75
Specifled: Laplace
0.25
0.00
proportion perents sorted after
0.75
Specifled:
0.50
0.25
0.00
Specilled: Scaledi4 (df-10)
0.75
0.50
0.25
(di-1
0.5
10
0.5
2
10
50
0.5
10
0.5
10
0.5
2
10
50
0.5
2
50
ratio: n/p
algorithm DiretLINGAM
HighDimLINGAM
ScoreL/NGAM
Figure 2: The simulation results comparing LiNGAM estimation procedures. Color indicates sorting procedure. The y axes each plot indicate sorting estimation error according to Equation []] The columns of the grid of plots corresponds to the bnlearn. com network architecture used. The rows indicate the true and specified scale-location family.
Figure 2 reports the results in terms of order estimation error (lower is better), which we define
Err(+) =
Σ (G,R)EE
1 > [k)]
(11)
Recall from 2.1 that & is the edge set of the true DAG G. The indicator 1 I # (k)) in Equation 11 states whether the ordering with á erroneously placed parent node j after child node k. The normalizing constant
E
ΣΣ 1[B #0, j=
the total number of directed edges in G, gives the order estimation error Err(f) value between and Given that the target paramcter of our work topological ordering as in Definition 1.1, the error metric in Equation [ concisely captures the extent of violation to this definition. Each edge j k in the true topological ordering, In this regard, Err(A) is the fraction of
DAG imposes an independent constraint these constraints not satisfied by the estimated ordering 4. If a node with many parents is placed at the start of an estimated ordering before many of those parents, then this error metric will penalize proportionally. Likewise, this penalty will occur proportionally if node with many children nodes is sorted after many of those children.
Our ScoMeLiNGAM achieved the highest accuracy for all high-dimensional settings (n S p). DirectLiNGAM became quite comparable until the sample size increased to 2p and did bit better than ScoreLiNGAM
11

---

Published in Transactions on Machine Learning Research (12/2022)
mildow (p=35)
alarm (p=37)
barley (p=:88
hailfiindr (p=56)
hopar2 (p=70)
winStpts (p=76)
andes
(p=:22)
0.75
Specified:
0.50
0.25
proportion of perents sorted
0.00
0.75
Specified: Scaled- (off-10)
0.50
0.25
0.00
0.5
10
0.5
2
10
50
50
0.5
10
50
0.5
10
0.5
2
50
0.5
2
50
ratio: n/p
algorithm DirectLINAM
HighDimLINGAM
SoaMeL/NGAM
Figure An extension of the simulation results in Figure Here, we mis-specify the noise distribution for ScoreLiNGAM.
2 comparing LiNGAM estimation procedures.
when n 2 10p (large sample size cases). Note that results are not presented for DirectLiNGAM when n 0.5p = p, because it is not applicable for p. For the Andes network, results for DirectLiNGAM are also not presented
this procedure takes about 118 minutes for single replicate, which adds up across 90 total replicates. On the other hand, HighDimLiNGAN is generally the least accurate algorithm across all networks and sample sizes. Recall that the data matrix X is re-scaled. The lack of accuracy change for Hig likely owed to the fact that this procodure is not invariant to a re-scaling the data, as ScoMeLiNGAM and NGAM are. We also note an apparent. good property of DirectLiGAM: does not require the MB input, yet it performs well terms of accuracy large values of
We also compared the three methods when the error distributions were mis-specified for ScoreLiNGAM in Figure 3 The true error distributions were Logistic or Scaled-t, but we still used the Laplace update of Equation 10 in ScoreLiNGAM. It is seen that its accuracy was comparable to the result when we correctly specified the error distributions (the other three columns), suggesting some robustness of our method to mis-specification in the error distribution family.
In terms of spoed, Figure summarizes this for the win95pts network. The advantage of our method is speed, with our method being less than 100 times faster the next fastest method. Note: HighDimLiNGAM's procedure is parallelized across 7 threads. For the setting with n/p 0.5 the win95pts network, ScoreLiNGAM is more than 1000 times faster than HighDimLiNGAM. Appendix B contains details about the implementation of each procedure, along with the machine used to run these experiments. Moreover, Figure in Appendix B contains sorting times for all the settings we considered; there is a similar pattern of ScoreLiNGAM's capability in terms of speed compared to the other methods.
In Appendix B.3 we present results on the barley, alarm, and mildew networks for the case that the data matrix is not re-scaled the results of this subsection. These results support our claim that HinADM
12

---

Published in Transactions on Machine Learning Research (12/2022)
Specified: Laplace
True: Laplace
log10 sorting time in seconds)
2
WId56uJM
-1
-2
-3
0.5
10
50
ratio: n/p
algorithm
DirectLiNGAM
HighDimLINGAM
ScoreLiNGAM
Figure 4: The log10 (avg. sorting time in seconds) scale for the various methods applied to the win95pts network.
is sensitive to re-scaling, while DirectLiNGAM and ScoreLiNGAM not. Next, Appendix B.4 obtains results on these same three networks for the case that neighborhood sets must be estimated. These results help us understand the sensitivity of our method to proper neighborhood specification.
proportion of parents sorted after a child
0.2
p=5000
0.1
0.0
0.2
p=10000
0.1
0.0
0.1
0.5
ratio: n/p
neighborhood
Most Correlated
Markov Blanket
Figure 5: Sorting errors for ScoreLiNGAM under p 5000, 10000 and 0.1p, 0.25p, 0.5p. Color indicates how the neighborhood sets are constructed.
13

---

Published in Transactions on Machine Learning Research (12/2022)
Mean Log-L Likelihood (Test Set)
-1250
-1275
-1300
-1325
-1350
Gaussian
Laplace
Residual Density Specification
Permutation Method
ScoreLiNGAM
Sparsebn
Figure 6: The mean log- liellikel on 1,000 genes for subsct of cells in the data of Yaao al. 2021 across 50 repetitions.
3.2 Larger Network Results
We also simulated large networks with 5000, 10000 and n = 0.1p, 0.25p, 0.5p to further demonstrate the scalability of ScoreLiGAM. We do not include results in these settings for DireetLiNGAM nor High- DimLiNGAM as they would take too long to run. The network generation is such that 5% of nodes are root nodes (no parents), and all other nodes have between 1 and 2 parents (with equi-probability which are selected at random from the set of predecessors in randomly generated permutation. Morcover, Lid. Býj ~ Uniform|-0.9, -0.4] U [0.4, 0.9] across (j,k) such that j E PAs, while 0₂ Lid. Uniform(0.255 0.9] across 1 1 k < is the scale parameter for the Laplace noise in the synthetic LiNGAM. A new LiNGAM is generated according to this schema for each data replicate.
Figure presents simulation results for ScoreLiNGAM with two different priori known neighborhood sets. "True Markov Blanket" means that set N&t M Bx for each with these oracle sets. The results for "10 Most Correlated" use 20% of the data to specify Nut the 10 most Pearson-correlated variables (in absolute value) to X for each 1 < k S and the other 80% of the data to estimate the topological ordering.
S p and run the sorting procedure
It encouraging to see in Figure 5 that the accuracy our method is high even for such a challenging high-dimensional setting. In fact, the average error rate quite comparable to that for the smaller networks reported in Figure 2 As expected, an accurate neighborhood method relatively quickly for large p, but its accuracy naturally is dictated by sample size. Figure in Appendix B contains the sorting times to go along with Figure 5
provides better sorting results. Further, our
3.3 Application: Single-Cel Cell Gene Expression Data
We present application of our method the data of Yao et al. 2021 We scek lincar SEM to model gene regulatory network, where each X₃ in Equation is the expression level of gene. We focused our attention their dataset for which isolated single cells were processed for RNA sequencing using SMART-Seq v4 (labeled "Mouse Cortex--lipppocampp (2019/2020)7 2020)"). Noting Yao et al. (2021) h'ss finding that cells' gene expressions cluster according to region and cell type, the focus is a subset of the data
Available http:////hh...........c..ccccaa.
compressed TSV format.
14

---

Published in Transactions on Machine Learning Research (12/2022)
follows. We focused on glutamatergic cells from the mice brains' primary visual cortex. We also focus cells for which injection materials are not specified (see Saleeba et al. (2019) for background on neuronal tracers). This takes us from 74,973 cells down to 7,159 the largest subset of all cell class, isocortex location, and injection material combinations. Moreover, a sizable amount of genes had expression mcasurements of exactly 0, so genes are also subset to those which were measured to be non-zero in 50% or more of these cells. This took us from 45,768 to 10,012 gene expression measurements taken across 7,159 cells.
3.3.1 Comparison to another scalable linear SEM estimation procedure
As for large simulated networks in Section 3.2 DirectLiNGAM and HighDimLiNGAM were too slow for this application. In order to compare ScoreLiNGAM another linear structural equation modeling procedure, we applicd the package sparsebn {Aragam et al., 2019) our data. This method is a scoro-basod method that maximizes a regularized Gaussian likelihood over the DAG space based on Aragam & Zhou (2015). To make comparisons across 50 repetitions computationally viable, we randomly select p 1000 of the original 10,012 genes. For each repetition, we randomly sampled 2,000 colls: half of the cells (n = 1000) were designated to be in the training set, and the other half in test set; cach data matrix was standardized such that columns have sample standard deviation 1 and sample mean 0.
In the training set, 20% of cells were randomly selected to estimate the Pearson correlation matrix. We specified the neighborhoods, Ne. for ScoreLiNGAM as the 50 genes j E [1000]] {[+}} with the largest Pearson correlation (in absolute value) with gene k. The remaining 80% of training data was used to estimate a topological ordering and the linear SEM's coefficients (via ordinary least squares). For Sparsebn, no priori neighborhood selection was used: parent sets for the linear SEM were learned with 100% of the training data using default options in the estimate. dag command, and the selection of the final DAG in the solution path was done by the recommended select. select.paramen parameter command. For Sparsebn, the linear SEM's model parameters were estimated according to the selected DAG, Moreover, the noise densities fit to the residuals in the training set are either Gaussian or Laplace.
As can be soen in Figure G the Laplace density specification for the additíve errors provides significantly higher mean log-likelihood on the test set compared to Gaussian density for both methods. This shows that the Laplace distribution, with its thicker tails than the Gaussian distribution, fits this data better. Furthermore, ScoreLiNGAM showed substantially higher test-data likelihood than Sparsebn under both error distributions for calculating the likelihood.
3.3.2 Application of ScoreLiNGAM to All 10,012 Genes
We now present the application of ScoreLiNGAM to all original 10, 012 genes discussed at the start of this section. The application is as follows:
1. We randomly split the gene expression measurement matrix with 7,159 cells (rows) into two folds having 358 (~ 5%) and 6801 (~ 95%) of the cells, respectively.
2. On the first fold, we ran neighborhood linear regressions in which we restrict coefficients to be non-negative via R's glmnet package (Fricdman et al. 2010) with no ridge or lasso regularization. We then selected the sets N& to correspond to genes N₂ S [p]{ such that coefficients are non- zero. The non-negative coefficient constraint for each linear regression, known as non-negative least squares, can itself be seen as a form of regularization that gives a sparse solution to the cocfficient vectors (Slawski & Hein 2013). Indeed, this constraint. the neighborhood regressions resulted neighborhood sets of size 45 to 277 genes (of 10,011 possible genes), with a median of 145 genes per neighborhood set. Our use of non-negative least squares regression is motivated by the use of non-negative linear regression to impute single-cell gene expression measurements (Li & Li 2018 It is also computationally faster compared to neighborhood Lasso regression.
3. We then randomly split the second fold into two folds, 2A and 2B, having 3401 (~ 47.5%) and 3400 (~~ 47.5%) of the original cells, respectively.
4. On Fold 2A, we estimated the permutation á using ScoMeL
15

---

Published in Transactions on Machine Learning Research (12/2022)
Adjusted
Unadjusted
0.190
0.185
50th Percentile
0.180
0.175
0.170
0.165
Coefficient of Determination's Quantiles
0.160
0.155
0.150
0.40
80th Percentile
Across all 10,012 Genes
0.39
0.38
0.37
0.36
0.35
0.34
0.54
90th Percentile
0.53
0.52
0.51
0.50
0.49
95th Percentile
Random
ScoreLiNGAM
Random
ScoreLiNGAM
Permutation Method
Figure 7: Across 30 replications, comparison of the estimated cocfficient of determination on Fold 2B for each gene. We summarize the coefficient of determination across genes by taking the median, 80th, 90th, and 95th percentiles.
5. On Fold 2B, the validation fold, estimated via linear least squares regression the 10,012 by 10,012 weighted adjacency matrix for the DAG corresponding to the topological ordering defined by A and such that the support of the k-th column is the set of indices N₂ n (#())oosict=00))
6. Using the weighted adjacency matrix estimate from the previous step, we then calculated linear least squares residuals from the predictions given by the estimated parent sot for each node k, on Fold 2B.
With the residuals the last step, we calculated the coefficient determination (RR): for each gene 1 E [10, 012), R3 an estimate of the proportion of variation explained by lincarly regressing gene k's measurement. the estimated parent. genes. It is minus the ratio of gene k's residual sum of squares and the total sum of squares.
The result of this application is summarized in Figure In order make a comparison, also calculated residuals from the linear SEM induced by randomly generated topological ordering, denoted in Figure 7 as "Random." The two linear SEMs, with topological ordering given by ScoreLiNGAM or randomly generated, select parent sets as the intersection of a node's neighborhood set, Ng, and the node's predecessors in the corresponding topological ordering if any. We repeat Steps 3-7 a total of 30 times.
Considcring that some genes may have more estimated parents than others, and that a coefficient of determination can be artificially large as the number of regressors increases, Figure also includes the adjusted coefficient of determination which incorporates a penalty for the number of regressors (Neter et. al.
16

---

Published in Transactions on Machine Learning Research (12/2022)
1983
The adjusted coefficient of determination is
n - 1
(1 R
RB - |PAR
where |PA&] is the number of estimated parents for node i and ng 3400 is the sample size in Fold 2B.
As we can see from Figure 7 ScoreLiNGAM gives higher coefficients of determination (adjusted and unadjusted) on the test datasets (validation fold 2B) compared to a randomly generated permutation across all random replications-at as summarized by the median, 80th, 90th, and 95th percentiles taken across the 10,012 coefficients of determination. Taking for granted the linearity assumption, the above higher-end percentiles of the R2 across all genes provide meaningful comparisons because it may very well be the case that a majority of the genes have quite random expression patterns. Based on the 90th percentile for the adjusted coefficient of determination in Figure 7 it appears that for 10% of all the genes, more than 51% of their expression variation is explainable by its estimated parents in the linear SEM given by ScoreLiNGAM. As shown by the 95th percentiles in Figure 7 for the top 5% of genes in terms of adjusted R², 79% or more of the genes' expression variation is explainable by its estimated parents in the linear SEM given by ScoreLiNGAM. These H² levels are significantly higher than those given by random permutations, with no overlap in the boxplots in Figure
Across the 30 replications, ScoreLiNGAM's sorting time for all 10,012 genes had median of 10.28 minutes, confirming its scalability for such large and high-dimensional datasets of p 10, 000 and 3, 000.
4 Discussion
In this paper, we claimed that sequentially applying Algorithm ! can give interesting structure learning results. We demonstrated this with a novel sequential procedure based on parametric specification that provides an alternative to the state of the art for the identifiability and estimation of a linear DAG model with non-Gaussian errors. We discussed the conditions, Assumptions 2.1 and 2.2 under which the proposed causal discovery procedure will identify the valid DAG. We also proposed a relatively simple procedure that can make strategic use of an a priori known neighborhood set for each node. Finally, we presented numerical evidence that our procedure scales to large dimensions, which is otherwise not the case for the state-of-the-art. of LINGAM. We accompanied these simulations with real-data application.
Our motivation for learning topological ordering is that it is relatively simpler problem to tackle compared to searching across the expeeeexponnttallly large space of DAG structures. Under Assumption 2.3 and regardless of whether the order estimate satisfics Definition 1.1, a graph output by our Algorithm 2 will have false positive edges corresponding to at least the number of times two non-adjacent nodes have a common child in the true graph. This is because Markov blanket contains all co-parents (or "spouses"). Moreover, the number of incorrectly oriented true edges will be the total number of true edges times the accuracy metric in Equation equation []] Thus, to produce a good graph estimate given our estimated ordering, one needs an additional step to apply sparse regression of each X, on the variables sorted in front of it to estimate the support (parent set). This is the well-studied support recovery problem in sparse linear regression.
With regard to the linearity, scale-location family, and a priori known Markov blanket assumptions for Theorem 2.5 we concede to the apparent limitation compared to assuming non-linear relations and making less parametric assumptions. However, to us, there appears to be a trade-off between more general modeling assumptions and scalability to a large number of variables, the latter of which our paper addresses. Recall the application to the large-scale gene expression data in $ 3.3.2 which demonstrates our method's improvement at lincarly explaining variation in large number of genes compared to a random permutation bascline. The concern about whether the entries of the error vector need to follow the same distribution family is an important one. It is not immediately clear how our current theoretical results, making use of information thcoretic quantitics, will come about without this assumption. Extensions of our work may like to explore whether this assumption can be relaxed. This relaxation is suggested by our sinulation results where the errors' distribution family is mis-specified.
17

---

Published in Transactions on Machine Learning Research (12/2022)
As practical manner, consider prospective applications to single cell gene expression data (scRNA-seq) in 3.3 Recent work suggests hierarchical structure between true (hidden) expressions and measured expressions with missing, and possibly zero, counts (Sarkar & Stephens 2021). Should the hierarchical nature be justified, further work on causal discovery for gene co-expression models may need to incorporate the fact that what we really would like a graphical model, possibly causal if a domain expert agrees, on the true (hidden) expressions. Along these lines, future causal discovery procedures for such data can build on the procedures of McDavi et al. (2019) and Yu et al. (2020), which themselves build on Gaussian graphical models, using heavier tail distribution for residuals we do here.
Further extensions of the work presented here include formal statistical guarantees along with extensions of the likelihood ratio approach to nonlinear SEMs.
5 Broader Impact Statement
Our scalable causal discovery method comes with extra modeling assumptions, not to mention the fact that we are making use of the causal sufficiency assumption (no unobserved confounders or selection bias). Therefore, it is imporlant for a responsible practitioner to not over-extend their conclusions with this or related data mining methods.
References
Aragam and Q Zhou. Concave Penalized Estimation Sparse Gaussian Bayesian Networks. Journal of Machine Learning Research, (6669))2273228)) 2015. URL http::/jmmr.... http://aww..oo/r/pppppr///////g/ html
B. Aragam, Jiaying Gu, and Qing Zhou. Learning Large-Scale Bayesian Networks with the sparsebn Package. Journal of Statistical Software, 91(11):1-38, 2019. doi: 10.18637/jss.0991.111 URL https
//ww. www.jstatsoft.com org/ index. phtp:jsss/artiill///ee/e////
P. Bühlmann, J. Peters, and J. Ernest. CAM: Causal Additive Models, High-Dimensional Order Search and Penalized Regression. Ann. Statist., 42(6))2526 2556, 12 2014. doi: 10.1214/14-AOSLI URL
https://doi org/ 10. 1214/14-11111 A0S1260
W. Chen, M. Drton, and Y.S. Wang. On Causal Discovery with an Equal-Variance Variance Assumption. Biometrika, 106(4) (06(4))9733-4) 990, 09 2019. ISSN 0006-3444. doi: 10.1093/bione /s22049. URL biomet/as2049
Chickering. Optimal Structure Identification With Greedy Search. Journal of Machine Learning Research, 15507-554, 01 2002. doi: 10.1162/15324433333111817
D. M. Chickering. Learning Bayesian Networks is NP- NP-Complete, pp. 121-130. Springer New York, New York, NY, 1996. ISBN 978-1-4612-240444 doi: 10. 111000//7----41122220000
12. URL https: //dgi.
12
T.M. Cover and J.A. Thomas. Inequalilies in Information Theory, chapter 17, pp. 657-687. John Wiley & Sons, Ltd, 2005. ISBN 9780471748823. doi: https:///oi.oo /ro.//////og////////o///o woo.ooooo10111.00000000.01.1 1000 //omlinelibrary, lalleyycom.. http://///ab/11111100//777888380
URL https
ch17
C. Forbes, M. Evans, N. Hastings, and B. Peacock. Statistical Distributions. Wiley-Blackwell, Hoboken, NJ, 4th edition, November 2010.
J. Friedman, Hastie, and R. Tibshirani. Regularization Paths for Generalized Linear Models via Coordinate Descent. Journal of Statistical Software, 33(1):1-22, 2010. doi: 10.18637/jss..033.. URL https: :
//www. jst.tssftttoorgg.o.oo.o. ortpe///t.org/33//111
M. Frydenberg. The chain graph Markov property. Scandinavion Journal of Statistics, 17:333-353, 1990.
18

---

Published in Transactions on Machine Learning Research (12/2022)
M. Gao, Y. Ding, and B. Aragam. A Polynomial-Time Algorithm for Learning Nonparametric Causal Graphs. Tn H. Tarochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Newral Information Processing Systems, volume 33, pp. 11599 11611. Curran Associates, Inc., 2020. URL https: Paper....p
A. Ghoshal and J. Honorio. Learning Identifiable Gaussian Bayesian Networks in Polynomial Time and Sample Complexity. In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS'17, pp. (460-6469, Red Hook, NY, USA, 2017. Curran Associates Inc. ISBN 9781510860964.
C. Glymour, K. Zhang, and P. Spirtes. Review of Causal Discovery Methods Based Graphical Models. Prontiers in Genetics, 10:524, 2019. ISSN 1664-8021. doi: 20389//gene.201990
URL https:
//www. www.rotiiessoooogrrrrtiiiiiioiii
3389/ffen..e9eeeeeeeeee 00524
L. Henckel, E. Perković, and M.H. Maathuis. Graphical criteria for efficient total effect estimation via adjustment in causal linear models. Journal of the Royal Statistical Sociely Series B, 84(2):579 599 April 2022. doi: 10.1111/rsb.12451. URL https://idees html
repec. wtw.:aas.rrpppc.rrreeeeoee http://wwwppppoooo//ooooo/oo/aaa
A. Hyvärinen and E. Oja. Independent Component Analysis: Algorithms and Applications. Neural Networks, 13(4):411-430, 2000. ISSN 0893-6080. doi:
URL https
//wwx www.sciencedirecc. http://cw.................oaoaao
A. Hyvärinen and S.M. Smith. Pairwise Likelihood Ratios for Estimation of Non-Gaussian Structural Equation Models. Journal of Machine Learning Research, 14(1):111 152, January 2013. ISSN 1532-4435.
J.R.J. Lec, R. Srinivasan, C.S. Ong, D. Alcjo, S. Schona, Shpitser, M. Sussman, G.J.R. Whitman, and D. Malinsky. Causal determinants of postoperative length of stay in cardiac surgery using causal graphical learning. The Journal of Thoracic and Cardiovascular Surgery, 2022. ISSN 0022-5223. doi: http://ps//www..eaccacaaaa.o.ooa
pi1/80022522322000000
W.V. Li and J.J. Li. An Accurate and Robust Imputation Method scImpute for Single-Cell RNA-Seq Data. Nature Communicationss, 9(1), March 2018. doi: 10.1188/4411667-01183330------ URL https:
adoi..org/0. (008//441477001-03200
D. Malinsky and P. Spirtes. Estimating bounds on causal effects in high-dimensional and possibly confounded systems. International Journal of Approrimate Reasoning, 88:371 384, 2017. ISSN 0888-613X. doi:
https///ooooor//00000///iiiaeii 200.000000.. URL http::/wwww..ccaeeeeeeeee...co.s
pJ1/S0888613117303822
D. Malinsky, I. Shpitser, and T. Richardson. A Potential Outcomes Calculus for Identifying Conditional Path- Specific Effects. In Kamalika Chaudhuri and Masashi Sugiyama (eds.), Proceedings of the Serenty-Second International Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings of Machine Learning Research, pp. 3080 3088. PMLR, 16 18 Apr 2019. URL https://proceedingss mlr www.mml...pp..eeeemeemme ptwss//8// malinsky1bb.hta
A. McDavid, R. Gottardo, N. Simon, and M. Drton. Graphical Models for Zero-Inflated Single Cell Gene Expression. The Annals of Applied Statistics, 13(2):848- 873, 2019. doi: 11224/18-AAAS1200 URL
https://doi.org org/ror//11 10. 1ew.1/21411111 1214/18-A0A81221 ADAS1213
N. Meinshausen and P. Bühlmann. High-Dimension Dimensional Graphs and Variable Selection with the Lasso. The Annals of Statistics, 34(3):1436 - 1462, 2006. doi: 10.1214/00053000000022 URL https
//doi //doi.org/10.1t /tte//0..22//0005500000000000000
P. Nandy, M.H. Maathuis, and T.S. Richardson. Estimating the Effect of Joint Interventions from Observational Data in Sparse High-Dinenssonall Settings. Ann. Slatist., 45(2):647766) 04 2017. doi: 10.1214/16AAOS1622 URL https://doi.orgg: 111.111-11111111111111
J. Neter, W. Wasserman, and M.H. Kutner. Applied Linear Regression Models. Trwin, Homewood III., 1983.
19

---

Published in Transactions on Machine Learning Research (12/2022)
G. Park. Identifiability of Additive Noise Models Using Conditional Variances. Journal of Machine Learning Research, 21 (75):1-334 2020. URL http://jnlr... httpp//rrrrgppppee//////119999 html
G. Park and Y. Kim. Identifiability of Gaussian Linear Structural Equation Models with Homogeneous and Heterogeneous Error Variances. Journal of the Korean Statistical Society, Jan 2020. URL https: //link. springer. com/ /ttp///articee/...00///222/2222/ 019-00019-7
J. Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press, USA, 2nd edition, 2009. ISBN 052189560X.
J. Peters, J.M. Mooij, D. Janzing, and B. Schölkopf. Causal Discovery with Contimuous Additive Noise Models. Journal of Machine Learning Research, 15(58):2009-200) 2014. URL http://jmlr. org/papers/ v15/peters14a htmml html
J. Peters, D. Janzing, and B. Schölkopf. Elements of Causal Inference: Foundations and Learning Algorithms. MIT Press, Cambridge, MA, USA, 2017.
G. Raskutti and C. Uhler. Learning directed acyelic graph models based sparsest permutations. Stat, 7(1): e183, 2018. doi: https:/ https://doi.owr ////i..org/1..100////taaaaa 1002/sta4..83. URL https://onnleeelbbaarry www.eeewwlarr..e.rr....ar.eeeeee wtw.ley.comddoil.. 10.1002/se sta4. 183 e183 sta4.183.
J. Runge, S. Bathiany, E. Bollt, G. Camps Valls, D. Coumou, E. Deyle, C. Glymour, M. Kretschmer, M.D. Mahecha, J. Muñoz-Mari, E.H. van Nes, J. Peters, R. Quax, M. Reichstein, M. Scheffer, B. Schölkopf, Spirtes, G. Sugihara, J. Sun, K. 'Zhang, and J. Zscceischherr Inferring Causation from Time Series in Earth System Sciences. Nature Communications, 10(1), June 2019. doi: 00.108//441467-00 10100--- URL
C. Saleeba, B. Dempsey, S. Le, A. Goodchild, and S. McMullan. A Student's Guide Neural Circuit Tracing. Frontiers in Neuroscience, 13:897, 2019. ISSN 1662-453X. doi: 10.3389/fmins.2019. 2019.00097. URL
bttps://ww/ frontiersin wttpeiinrgraaaraceee/e000 http:/aarttclee1..338/////innn 2019. 00897
N. Sani, D. Malinsky, and I. Shpitser. Explaining The Behavior of Black- Box Prediction Algorithms With Causal Learning, 2020. URL https://arxiv.org/abtt
02482
A. Sarkar and M. Stephens. Separating Measurement and Exprcssion Models Clarifies Confusion in Single- Cell RNA Sequencing Analysis. Nature genetics, 53(6))770 777, June 2021. ISSN 1061-4036. doi:
URL https: httpps://all.orr/1100/3///////// http:/10.033///1188-000111111111
S. Shimizu, P.O. Hoyer, A. Hyvärinen, and A. Kerminen. A Linear Non-Gaussian Acyclic Model for Causal Discovery. Journal of Machine Learning Research, ((72))2003-220) 2006. URL /ttp:///// org/papers/ 97/shimizu06a.hn html
S. Shimizu, T. Inazumi, Y. Sogawa, A. Hyvärinen, Y. Kawahara, T. Washio, P.O. Hoyer, and K. Bollen. DNGA A Direct Method for Learning Linear Non-Gaussian Structural Equation Model. J. Mach. Learn. Res., 12:1225 1248, July 2011. ISSN 1532-4435.
M. Slawski and M. Hein. Non-Negative Least Squares for High- Biigh-Dimessioall Linear Models: Consistency and Sparse Recovery Without Regularization. Electronic Journal of Statistics, 7:3004 3056, 2013. doi: 10.1214/13-E5 11-11/13--EJ886.. URL
L. Solus, Y. Wang, and C. Uhler. Consistency Guarantees for Greedy Permutation-Bant Based Causal Inference Algorithms. Biometrika, 108(4):795-81) 814, 01 2021. ISSN 0006-3444. doi: 10.1093 1.11099/bbomee,, biomet
asaa104. biomet/asaal04 URL
https://doi.org org/ 10. 1tt..00000//mmmm//aaaaa/00.
P. Spirtes and Glymour. An Algorithm for Fast Recovery of Sparse Causal Graphs. Social Science Computer Review, 91))622-7, 1991. doi: 10.1177/0994489910000000 URL https: http:///oo.oo//1.1177///
0894459331109900106
20

---

Published in Transactions on Machine Learning Research (12/2022)
P. Spirtes, C. Glymour, R. Scheines, S. Kauffman, V. Aimale, and F. Wimberly. Constructing Bayesian Network Models of Gene Expression Networks from Microarray Data. 1 2000. doi: 10.1184/
R1 R1/6491229.v1. URL https://kiittub. cmu wwwp///ticle///ooeeeeeeeeeeeee/e journal lawpuuaaalooorituttittaaaaaaeeee
Bayesian_ Network BaatarrKModell_eee of Gene Bearesssio__eewooeee Bacress_ffoooMMooM Microarray Datpp/aaaaaaaal///////
D.J. Stekhoven, I. Moraes, G. Sveinbjörnsson, L. Hennig, M.H. Maathuis, and P. Bühlmann. Causal Stability Ranking. Bioinformalies, 28(21):2819 2823, 09 2012. ISSN 1367-4803. doi: 10.1093/bioc www.n/oooaaaiis/eaaaeee. bbt5233.
URL https://doi http:///ooo.ogggg.....//////ooee
T.S. Verma and J. Pearl. Equivalence and Synthesis of Causal Models. Probabilislic and Causal Inference, 1990.
T.S. Verma and J. Pearl. Equivalence and Synthesis of Causal Models, pp. 221-236. Association for Computing Machinery, New York, NY, USA, 1 odition, 2022. ISBN 9781450395861. URL https :
//doi.org/11 10.g//00... 1145/3501714.350173
Y. Wang, L. Solus, K. Yang, and C. Uhler. Permutation-Based Causal Inference Algorithuss with Interventions.
I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.). Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL. https
//proceedings.newne wwwpiiingg..nuri......mm http::iigg...uurppppp...///////p hrtp:paww....iiiooooooooooooooan
Y.S. Wang and M. Drton. High-dimensional causal discovery under non-Gaussianityy 10 2019. ISSN 0006-3444. doi: 10.1093/biie biomet/as2055. URL
107(1):41 59
Z. Yao, C.T.J. van Velthoven, T.N. Nguyen, J. Goldy, A.E. Sedeno-Cortes, F. Baftizadeh, D. Bertagnolli, T. Casper, M. Chiang, K. Crichton, S.L. Ding, O. Fong, E. Garren, A. Glandon, N.W. Gouwens, J. Gray, L.T. Graybuck, M.J. Hawrylycz, D. Hirschstein, M. Kroll, K. Lathia, C. Lee, B. Levi, D. McMillen, S. Mok, T. Pham, Q. Ren, C. Rimorin, N. Shapovalova, J. Sulc, S.M. Sunkin, M. Ticu, A. Torkelson, H. Tung, K. Ward, N. Dee, K.A. Smith, B. Tasic, and H. Zeng. Taxonomy of Transcriptomic Cell Types Across the Isocortex and Hippocampal Formation. Cell, 184(12):3222 3241.e26, 2021. ISSN 0092-8674. doi: https://ddoiorrg/01110////// 2021................... URL pi1/8009926644210000188
Q. Ye, A.A. Amini, and Q. Zhou. Optimizing Regularized Cholesky Score for Order-Based Learning of Bayesian Networks. IEEE Tronsactions on Pattern Analysis and Machine Intelligence, (5150::555-337) 2021. doi: 10.1109/TPAMI. 200...20220222222222
S. Yu, M. Drton, and A. Shojaie. Directed graphical models and causal discovery for zero-inflated data, 2020. URL httpp:// org/ htgggbb///0...... 04150
Y. Zeng, Z. Hao, R. Cai, F. Xie, L. Ou, and R. Huang. A Causal Discovery Algorithm Based on the Prior Selection of Leaf Nodes. Neural Networks, 124:130-145, 2020. ISSN 0893-6080. doi: https:
(ttol.org/10.101/j..mmm www.ee.te.....2......... URL /tt.htt/s///w.e.eeeeeeeee www.sseedeeeeect... httpa//cciccccccoooocc/ccccaaaaa $0893608019304204
X. Zheng, B. Aragam, P. K. Ravikumar, and E.P. Xing. DAGs with NO TEARS: Continuous Optimization for Structure Learning. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems, vol- ume 31. Curran Associates, Inc., 2018. URL https://proceedinng neurips. /c/ppaper/211/ffilee pdf
X. Zheng, C. Dan, B. Aragam, P. Ravikumar, and E. Xing. Learning Sparse Nonparametric DAGs, In Silvia Chiappa and Roberto Calandra (eds.), Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pp. 3414-3425. PMLR, 26 28 Aug 2020. URL htttp:://proceeeeeeeee ottp:/cceddnnnm.m...rrree///////
mlr wwtppppwwsesseeeeepeppppe wwwp//////hhen.aa.mmm html
21

---

Published in Transactions on Machine Learning Research (12/2022)
A Proof of Theorem 2.5
A.1 Formal Proof of Theorem 2.5
Proof of Theorem 2.5
Our proof boils down to making the correct decision in Algorithm ! at step 1, then making the correct choice at step 2 assuming the choice in step was correct, and so on.
For the sake of induction, let us assume that Ar correct the sense that PAa S Ae for all a E A, This is true at the basc case when A, 0, since having made no ordering choices also means we have made no mistakes.
Let h St be invalid node to continue the ordering in the sense that PAg n D. And let l E S₁ be valid node to continue the ordering in the sense that PA(E) S At.
Lemma 2.7 tells us that the least squares residual Re ~ Sec(re) is no closer to Gaussian than REi ~ Ske(TK) in the sense that:
DKL (JRe(ria) ( ((((a))))))))) S DEL (Ja(ra)llén(ra))
Furthermore, regularity Assumption
2.2 ensures that:
DEL (Fa((ra)llaaanree ma)) > 0.
On the other hand, long properly specified the error density for node E, we have that:
DEL (Je(ra)lloerrssnna)) O.
Thus,
log 9x (K The) <(rni ax))
(Tas)
DeL (Fu((ra)llo(russo))) DEL (oe(rr)llonrrreinn
ge(Reg (Res ma)
< Bie(tr)
log
Dei (felra)llu((ra))
der (Rex)
Altogether, this implies that
max S(j,Ai) > S(k, Ar). JES:
and
i = arg max S(j, Ac);
jES,
since and l were arbitrary valid and invalid nodes, respectively.
So at step f, we will always make the correct choice for a node to continue the ordering.
A.2 Proofs of Lemma 2.6 and Lemma 2.7
In this section, we formally prove Lemin 2.6 and Lemma 2.7
A.2.1 Some Useful Shorthand Notation
Let define some new strategic sets which contain indices in [p], and review some we have been using already.
22

---

Published in Transactions on Machine Learning Research (12/2022)
The set
0
t=
A (f.,,....,. ((L-)),,((111)) 122
This the partial ordering at step 1, 1,,.... In our population-level identification results, we will typically assume it is correct at step which means that for all a As, PAa C As.
S = A, is the set of unordered nodes at step t.
MB, PA U & UjECH, PAj is the Markov Blanket of node k.
Ng is the Markov Blanket superset such that N₂ 2 MB, In finite data, we will typically estimate Ng by procedure such as neighborhood lasso regression, so this containment may not hold. For the sake of this section, because we are deriving quantities at the population-level we assume that N₂. is known and contains the true Markov blanket. Note that trivially, we may consider N₂ = [p]][}), and the results of this section would still hold.
Nu = A. n No is the intersection of the Markov blanket superset with the partial ordering.
Lu UU= U (t) AN,, which are either nodes of Nu ancestors of nodes in Nei- When A is correct, it is necessarily the case that L S A, for each l & AA.
Lie the complement of set Lus which either contains nodes in A. which are not in Lut, i.e. the nodes of As Lkt, which are unordered, i.e we have that S. C kt
Note that for each node &s St we write:
X MA.E= MRLLEEL + MaLg E
(12)
where the second equality holds since Lus U ULL = [p]. We can similarly write
X N M N.. = M NarLe, ela
(13)
We omit term with ezg. since by definition of Lit, the sub- sub-mixing matrix M Nung is matrix. Combining Equation 12 and Equation 13
REt
(MEtu
B Neckad
elm + MaLg $L
which will make of in the proof for Lemma 2.6 below.
A.2.2 Proof of Lemma 2.6
Proof of Lemma 2.0
Case Assume i is a valid node to continue the ordering the sense that PA₃ S. A.. We want show that Rxt ég. In this case, M&Lg has non-zero entry corresponding to only M = 1. This is because AN₁ support (M+) A(k) C Lat, which in turn holds because PAB S MB n A C Ngc. Thus we have
MRL ELE =
So have left to show that
[MALZe BL N Tua
ELLE
0.
Rocall that B is the weighted adjacency matrix for the underlying LiNGAM. We have that support(B..) = PA Let us index the entries of the column vector B.g by NEs and denote this as B Nock" One thing that could be helpful to prove is that if & is valid, then:
3ke B Nack
23

---

Published in Transactions on Machine Learning Research (12/2022)
Because support(B..) PAg and PAg S Nes, consider that
fXx = B + = X
B
Neck + ck.
Nut
with Ck X X Neck and E[ce] 0. Thus,
But E
NuL
X 1)
Nui
Nxi
XI
Nui
Bunk k
+ E
X Xix ex
E
X Nue
x Nin 1)
E
X
Nu.
X7
B Nack
(14)
Nk.
= B Nicch"
wanted.
It follows that X2 B Nick Nxx Ck = X No. + +eg. This then means that Rgt = Xx X Nue = Exs
we wanted to show.
Case Assume not valid node. All we need this case for our identifiability proof that Rs is linear combination of more than one independent component. This is the case because k is invalid to continue the ordering, then have that there exists at least one E PAg such that St (unordered) and therefore j E I Recall that:
Rgt (MELUs BE. Neckie
CLa + MALG, C
Note that it is necessarily the case that Maj # O, otherwise j g PA₁. Thus, Ra. includes the sum Majej + ee- That is, Rgt in this case a linear combination of more than one independent component in e. Note that Rat could be a linear combination of more entries in e, in addition to and Ek.
A.2.3 Some Information Theory Definitions and Results
We now present some straightforward information theoretic results. They are meant to help demonstrate that our surrogate optimization (now likelihood ratio) approach for Algorithm leads the identifiability of causal order. These lemmas are used later to prove Lemma 2.7 key result that says valid nodes j in a LINGAM closer to Gaussian compared to invalid nodes k, conditional on the nodes in N3. and Nas, respectively.
Definition A.1 (Differential Entropy).
For a continuous random variable X with density p(x), denote Ep(x) []] to be expectation with respect p(r). The differential entropy of X is given by:
1
h h = Ep(e)
log p(X)
Lemma A.2 (Restatement of the entropy power inequality).
Consider two independent random variables X p(z) and Y p(y), and let X' ~ NN v(X')) and Y' ~ N (EYY), V(Y')) be independent random variables such that h(X) = h(X') and h(Y) = h(Y"). Then:
n(X+Y) 2 h(X' h(X' +Y").
Proof. This is exactly Theorem 17.8.1 of Cover & Thomas (2055). so refer the reader their proof.
24

---

Published in Transactions on Machine Learning Research (12/2022)
Lemma A.3 (KL Divergence from Gaussianity). Let X ~ p(x) and ((x) the densily for X ~ N (E[X). Cool (X])-
DEL (p(x))la(x)) h(X) hX)).
(15)
As in Lemma A.2. let X' ~ (((I") such that h(X) h(X).. divergence from Gaussianity as:
We can equivalently write the KL
V[X]
DEL (p(x)]la(a)) = 2 v[X"]
log
Proof.
Because
1 h(X) 1xx xxx((t)) log {(X) | Exme(e) { log ((X) S
11
by properties of this normal distribution (namely, that Ellog ((X)) x v[x] = v[X]) we have that:
DEL (p(x)lla(x)) = h(X) h(X).
Noting that the differential entropy for any N(4, o²) is 1 log(2neo²) and our assumption that h(X) = h(X)), we arrive at the second equality:
[
DEL (p(x)]la(x))
log
2 2
Note that also V[X] = v[X] always non-negative.
2 V[X = DKL (p()))))la()))) 2 0, which is the case because KL-divergence is
This well known result also implies that the normal distribution is the maximum entropy distribution when we constrain the first and second order moments of each distribution to be the same.
Lemma A.4 (Same distance to Gaussianity).
Let & NO0 v(ck)) with density qu(-) for each k 1,2..... .p Also let d NOO, v(()) such that h(e()) h(ex). If in our LaNGAM salisfies Assumplion 2.1 then there erists constand 7 2 0 such that
DKL (x)e::00)lll
and
Vie V[4
= = = exp(27)
for all k = 1,2,..., =1,,,...
Proof.
From Lemma
A.3 we have that:
DKL (o(es; (x) h(Tc) - h()))
Noting Assumption 2.1 and properties of differential entropy under rescaling, it follows that for U ~ g(: 0o): h(e&): h(U) + log(O./Oo).
Let U' N(0, V(U')) such that h(U') = h(U). We have also that
h(()) h(U') log 0z/0o))
25

---

Published in Transactions on Machine Learning Research (12/2022)
based the construction of both 4 and U²
Similar to Ek, let U ~ NO, V(U]). Thus, regardless of k = 1, ..
., we have that:
[ = h((x) h(())) = h(U) h(U") =:
log
2 2TEV(ck))
A.2.4 Proof of Lemma 2.7
Proof of Lemma 2.7 7
For each j E [p), let ef be a normally distributed random variable such that h(c))) = h(ee); while ēÿ is distributed as (B[e,) (e,)). Here, for all j, k E (1, p), قل (unless j = k) and E (even if
j= =).
Recall also that for St
Rg = (MIL. - - Rgyking Ryytnn ELLo = MILgELG = Σ 6efEs
selpl
where the coelficients dij in the last equality are used for shorthand. Note that dgf always. And if j is invalid to continue the ordering, then also dig * 0 for at least other i E [pliii,, based on Lemma 2.6
The relation between the quantities of interest is as follows:
DEL (falrra)lll(rrs h(RRe) h (Rur)
by Lemma A.3
h
Σ Kske
-
Σ
deues
Notioe: Rec
KE(n)
#E(p)
$ h Σ
Σ Saae
by Lemma A.2
tE(p]
<ELp)
(16)
2ne Σ iElp) 53, Var(e;
log
by normality of the Es,
25ee Lielo)
iE(p) 82, Var(e()
2ne
EeE(r)
o7 Var(e))
log
by Lemma A.4
Eie(r) 87 Var(c)
DxL (falr))llu(tr)))
as wanted (Recall Rec Ee by Lemma 2.6
B More Figures
B.1 Sorting Time for Small Networks
The takeaway of Figure E is that ScoreLiNGAM is generally much faster. Consider the largest DAG, the Andes network (p 223), where the sorting time of ScoreLiNGAM is typically under 1 second across all sample sizes, while for High (parallelized across threads) the sorting procedure takes between 10-1000 seconds across sample sizes. We note that ScoreLiNGAM is written with C++ using the Armadillo linear algebra library and an wrapper via the Repp package, while DiroctLiNGAM is written in Python (https: //github. com/catt1//ingan with wrapper function in using the reticulate package that written by this paper's authors. HiighimLiNGGA is also written in C++ (https: https::///gthbbbb/ wttetttuu..oom/ssaaaaaaaa wwt.///aaaaaaa///gg/iiii..
26

---

Published in Transactions on Machine Learning Research (12/2022)
with wrapper, but it searches regressor subsets when computing low-dimensional linear regressions likely reason for its slower time despite 7 parallel threads. All simulations were run a Dell XPS 13 with Intel CM 17-8550U CPU @ 1.80GHz x 8, 8 GB RAM, and 64-bit Ubuntu 20.04.3 LTS OS.
Specified: Laplace
Specifled:
Laplace True: Logistic
Specified: Laplace True: Scaled-t (dt=1))
Specilled: Logistic True: Logistic
Specifled: Scalle-1 Scaled-((df=10)) True: Scaled-t (ct=1))
True: Laplace
Lg=d Wlre
12224 sapue
log10(sorting time in seconds)
2
2
995
P-70
Diopar2
-2
2
9 2.00000
0.5
10
50
0.5
10
50
0.5
10
50
0.5
10
0.5
2
10
ratio: n/p
algorithm
DirecILINGAM
HighDimLiNGAM
ScoreLINGAM
Figure 8: The simulation times for LiNGAM estimation procedures
B.2 Sorting Times for Large Networks
Figure
contains the sorting times to go along with Figure 5 in the main text.
B.3
Further simulations: Results without re-scaling the data matrix
See Figure 10 for results on the alarm, barley, and mildew networks.
B.4 Further simulations: Results without re-scaling the data matrix and neighborhood estimation
See Figure for results on the barley, alarm, and mildew networks.
C Obtaining the scale- parametrrr for Empirical Mean Log-likelihood in Equation
As discussed in the main text, our sequential algorithm at step t 2 1 in practice requires the estimation of the scale parameter, The in Equation g Here, we discuss the cstimator for the three parametric assumptions used this paper. We make use the respective definitions and properties Forbes et al. (011
Laplace Distribution: If e ~ Laplace(O. 0x), we have that 0% E[|c is the scale parameter. When m(5) is specified as the density for Laplace mhr), the maximum likelihood estimator we use practice is îu 1 Rige
27

---

Published in Transactions on Machine Learning Research (12/2022)
p=5000
P=10000
90
sorting time in seconds
60
30
0.1
0.25
0.5
0.1
0.25
0.5
ratio: n/p
neighborhcod
Figure 9: Sorting times for ScoreLiNGAM under p 5000, 10000 and 0.1p, 0.25p, 0.5p. Color indicates how the neighborhood sets are constructed.
Logistic Distribution: If ek Logistic(O,, 0k), then 0% the scale parameter. We have that Var[ex] 0z When (((m)) is specified as the density for Logistic(O, me), we find that the plug-in estimator îike = 3 Cox to work satisfactorily.
Scaled-t Distribution: If 6 Scaled-t(O, u, 0&), then we say is equal in distribution to the scale paramcter, tz. times U ~ t(0, Student's t-distributed random variable having mean 0 and degrees of freedom v > O. That is, Tw)) is specified as the density for Scaled-t(O, v, mh)) with assumed to known, we find that the plug-in estimator îlk = o to work satisfactorily.
d OgU. For v > 2, we have that Var[ex] = 0f (1²)) When
In Equation and in the plug-in estimators for the Logistic and Scaled-t specifications, we use
2
Rke
kk
7
28

---

Published in Transactions on Machine Learning Research (12/2022)
Results without standardizing data matrix, known Markov blankets
mildew (p=35)
alarm (p=37)
barley
0.6
Specilled: Laplace
proportion of perents sorted child
Specifled: Logisiic
0.4
0.2
Specilled: Scaled- (dfe 10) Scaled+ (dix10)
0.2
0.0
0.5
50
0.5
10
ratio: n/p
algorithm
DirctLINGAM
HihDMLLNNGA
ScoreLINGAM
Figure 10: An extension of the simulation results in Figure comparing LiNGAM estimation procedures. Here, we do not re-scale the data matrix.
29

---

Published in Transactions on Machine Learning Research (12/2022)
Results without standardizing data matrix, 10 Most Correlated Nodes as Neighbors
alarm (p=37)
barley
(p=35)
(p=45)
Specifled: Laplace
proportion parents sorted after child
0.0
Specilied: Logistte
Specifled: Scaled4 (dix10)
0.6
0.4
0.5
10
50
0.5
2
10
50
0.5
10
50
ratio: n/p
algorithm
HighDiLLNGAM
ScorLLINGAM
An extension of the simulation results in Figure comparing LiNGAM estimation procedures. Here, do re-scale the data matrix. For ScoreLiNGAM and HighDimLiNGAM, we use 20% of the observations to estimate the correlation matrix and neighborhood sets, and the remaining 80% of obscrvations to istte the topological orderings. For DirecLLIGGAM, all 100% of observations are used to estimate the ordering without neighborhood pre-specification
Figure 11:
30


Reviews:
[{"review_id": "fcvlIAzb68", "reviewer": "Reviewer_ZRRu", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "While I am not an expert in this area, I generally found the paper to be well-written. For example, I found it very helpful when the authors provided text descriptions and inutions on what the various equations represented (though see below for some specific questions). However, due to the dense notation required for this topic, I believe the readability would be significantly improved by adding an appendix, table, etc., which lists the notation used in the main text.\n\nThe reproducibility of the work is somewhat limited. While the authors do provide an R notebook for running some small examples, it does not follow typical R package structure. More importantly, the simulated and processed data files are not available. Thus, it is not possible to reproduce the results.\n\nAs a non-expert, I did not verify the proofs in detail. However, they seemed reasonable to me. However, I do have several technical questions.\n\n* Do the entries of the error vector need to follow the same distribution family? For example, could one term follow a Laplacian and another follow a logistic? or even Gaussian? Do the theoretical results depend on them all following distributions in the same family?\n* What is meant by a “valid” vs. an “invalid” node? From Lemma 2.6, this seems to refer to a node consistent with the true topological order (in the population setting). The term is used several times before that, though, so it should at least be defined clearly before use.\n\nThe empirical evaluation is modest. A strength is that the authors investigate the case when the distributional family is mis-specified. However, the comparison to HighDimLiNGAM can be improved by using an experimental design similar to that in Figure 4, in which both methods are given an empirically determined subset of candidates for the MB, rather than the ground truth. Also, since HighDimLiNGAM is not invariant to re-scaling, it would be meaningful to include results on it with the unscaled data.\n\nA major limitation of the empirical analysis is that it is not very deep. For example, the accuracy of HighDimLiNGAM does not seem to particularly depend on the size of the dataset; on the other hand, the proposed approach generally improves as more data is available. It would be insightful to the differences which lead to this behavior. Similarly, it is striking that DirectLiNGAM does not require the MB as input, yet it typically outperforms ScoreLiNGAM in terms of accuracy; I appreciate that the scalability is an important aspect of the proposed approach (and that DirectLiNGAM is not an option when n <= p), but it would still be useful to point to specific differences in the methods that lead to these results.\n\nI also have a few specific suggestions about presenting the results\n* I may have missed it, but I do not believe it is ever stated what the box plots are showing (quartiles, standard deviation, variance, standard error, etc.)\n* Sorting Figure 2 by number of variables rather than alphabetically\n* Including the density of the networks, in terms of number of edges compared to the number of possible edges (or similar)\n* Show standard deviation or similar in Figure 3 (though it may be very low)\n* Split Figure 6 into multiple plots or similar; especially since the y-axis values change so much from part to part, it is very difficult to follow. Also, use similar formatting (i.e., filled-in boxes) as the other box plots.\n\nThe gene expression analysis would be improved with some semi-qualitative analysis comparing the learned relationships (and strengths) with wet lab-validated interaction data. For example, the STRING database (https://string-db.org/) is one open resource which includes such interaction data.\n\nTypos, etc.\n* The references are not consistently formatted.\n* The second sentence of Section 1 is a run-on sentence.\n* The first sentence in Section 1.2 is missing a verb.\n* When referring to specific “Definitions”, “Assumptions”, etc., the first letter is always capitalized, but it is never capitalized when referring to specific “equations.”", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this work, the authors propose an approach for learning variable (or node) orderings, according to LiNGAM assumptions, by iteratively adding one new variable to the order. The variable added at each iteration is selected based on a likelihood ratio between a non-Gaussian and Gaussian term. Under a given set of assumptions, the procedure guarantees to find an ordering consistent with the true structure. A finite sample version of the approach is also given, which largely entails fitting (and efficiently updating) regression models. A set of empirical experiments suggest that the proposed approach is competitive with SOTA in terms of accuracy and has better computational scalability.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Critical changes\n\n* Add an appendix (table, etc.) that contains the major notation used in the main text\n* Add the processed data files and brief documentation for running the experiments in the text\n* Clarify earlier about “valid” vs. “invalid” nodes\n* Perform analysis similar to that in Figure 4 for both HighDimLiNGAM  and ScoreLiNGAM for a few settings and datasets, such as small-data settings for alarm and win95pts\n* Address the various points about the plots, etc., mentioned above\n* Properly format the references, including proper capitalization, etc.\n* Address the various typos mentioned above\n\nSuggested changes\n* Package the source code as a proper R package (“DESCRIPTION”, “R” folder, etc.)\n* Clarify about the requirement on the same distribution family, especially if it reveals something interesting about the approach\n* Run HighDimLiNGAM on the unscaled data, for at least a few settings to confirm that the observed trends are not due to scaling\n* Address the observation that HighDimLiNGAM does not improve with more data, but the proposed approach (and DirectLiNGAM) do improve as more data is available.\nAddress the observation that DirectLiNGAM is typically more accurate than ScoreLiNGAM despite not being given the MB\n* Show a few examples of learned causal relations which have wet lab support (from STRING, etc.)\n", "broader_impact_concerns": "No major ethical concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "z2ZNwy4S32", "reviewer": "Reviewer_BwK3", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The proposed method offers a tractable solution that can be improved with prior knowledge.\n\n- The main contribution is presented well. The proposed likelihood ratio score is natural. The progression of theorem 2.5, lemma 2.6, and lemma 2.7 is easy to follow. The CLT argument following lemma 2.7 provided useful intuition.\n\n- The experimental setup demonstrates the ability of the method to scale where existing methods do not. In general, the empirical evidence supports the main claims.\n\nWeaknesses\n\n- The presentation of the empirical results could be improved (see requests and suggestions). The empirical results shown are based on a limited number of repetitions (30, 50, 90).\n\n- Evaluation of misspecified prior knowledge (Markov blanket) not presented.", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work considers the problem of estimating a topological ordering of a DAG. A sequential procedure greedily adds nodes to an ordering based on a novel likelihood ratio score. The authors show at the population level, with assumptions of a LiNGAM (Linear Non-Gaussian Additive Model) holding, a true ordering will be identified. This framework allows prior knowledge of Markov blankets of nodes to improve computational efficiency. Synthetic and real data experiments show that this method is tractable for larger graphs where existing methods are not. Overall, the paper is well written, and care is taken by the authors to articulate the proposed method. The empirical evaluation answer most, but not all, of the questions I have regarding the performance of the method.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Requested changes:\n\n- The first sentence of the first paragraph of Section 1.2 is missing an “is”.\n\n- Figure 2 may be trying to summarize too many results at once. I would recommend separating the cases of correctly specified error distributions from the misspecified cases. The first figure will show that the method accomplishes the goal it is theoretically set up to accomplish, while the second figure will show that it is robust to a potential misspecification, which do not need to be presented simultaneously. It also may be worth sorting the dataset by p, rather than alphabetically. Currently, it is a bit difficult to process all the information being presented at once.\n\nSuggested changes to strengthen work:\n\n- The simulation studies utilize a portion of the data to create prior knowledge of the Markov blanket, and it is clearly stated in that section that a consistent method should be used for estimation. However, earlier in the paper, it is suggested that the prior information may be available through domain knowledge, which is susceptible to being incorrect. How robust is this method to incorrect specifications? A brief simulation study, perhaps better suited for an appendix, would be worthwhile to have. The presented evidence show that correct prior knowledge can improve computation ", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Kv8n7nilJz", "reviewer": "Reviewer_6gVS", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n* The exposition of the problem setting and the main ideas of the proposed approach are clear.\n* The algorithm is easy to digest and leads to discovering the true ordering of the underlying DAG in the population case.\n* The approach is flexible enough to incorporate prior knowledge about the structure of the DAG.\n* The method can scale to high dimensions with competitive runtimes.\n\n### Weaknesses\n* There are some details about the method that I found unclear. I have a few questions that I will elaborate on below.\n* The sample complexity of the method is unclear and not discussed at all.\n* The experimental setting does not include the population case to corroborate the theoretical results. Moreover, I found part of the experimental setting somewhat strange; see below.\n* The approach would work only for linear models, which is a significant limitation for real-world problems.", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper studies the problem of learning the topological ordering of the underlying DAG of a structural equation model (SEM) with non-Gaussian additive errors, also known as LiNGAM. The proposed approach learns the topological order by first identifying a root of the DAG and then sequentially appending nodes that would lead to a topological ordering compatible with the true DAG. The appending procedure is based on a likelihood-ratio test. The authors show that in the population setting, the algorithm identifies the true ordering of the DAG. Experiments are conducted to validate mainly the scalability aspect of the proposed algorithm.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "My following questions/comments are more or less in order. I would appreciate it if the authors could clarify them.\n\n### Major (please address if possible)\n\n1. On page 3, paragraph 2, it is stated that the method runs in $O(pd)$ number of LS updates. It would help to see the **total** computational complexity in terms of $p,d$ for the population case and $p,d,n$ for the finite sample regime, and how it compares to previous work. This is to make a fair comparison against existing work. For instance, it is clear from Algorithm 2, that the overall computational complexity is not $O(pd)$.\n\n1. There is another work by Ghoshal & Honorio, \"Learning linear structural equation models in polynomial time and sample complexity,\" that tackles the linear setting for general error distributions. I suggest also comparing against this work.\n\n1. Why focus on the ordering of the DAG instead of the whole structure? I am aware that if one recovers the ground-truth ordering then discovering the edges becomes a much simpler problem. However, given that the method does not include finite-sample guarantees, I would like to see how the method compares against existing approaches when learning the whole DAG.\n\n\n**About Experiments**\n1. Is there any particular reason why the entries $B_{jk}$ cannot be greater than 1 in absolute value?\n\n1. The scale parameter for each error is somewhat similar. Why not sample from a larger range than $[.4, .7]$? I ask this because it is known that other methods that assume equal variance are robust to small changes in the variances. Since the proposed method does not assume equal variances, I would suggest letting the scale parameter be sampled from a larger range.\n\n1. There should be simulations to corroborate the correctness of the approach in the population setting.\n\n1. I am not sure if the metric in equation 11 is the right one. In particular, given that the graphs are sparse, it is likely that the topological orderings are not unique. For instance, assume we have a naive method that always outputs an empty graph, then any permutation is a valid topological ordering according to Definition 1.1. Thus, what topological ordering is used to compare against the proposed approach? I believe this would make the comparison unfair against other methods that output the whole DAG and can potentially have many valid topological orderings. \n\n\n### Minor (might strengthen the work)\n\n1. Why not simply write $\\hat{\\pi}(t) \\gets \\arg \\max_{k \\notin \\mathcal{A}_t} \\mathcal{S}(k,\\mathcal{A};\\mathbf{X})$ in Algorithm 1 and avoid the unnecessary loop?\n\n1. Is there a typo in Assumption 2.1? In particular, should $g\\left(e ; \\theta_k\\right)=\\frac{\\theta_0}{\\theta_k} g\\left(\\theta_0 e / \\theta_k ; \\theta_0\\right)$ instead be $g\\left(e ; \\theta_k\\right)=\\frac{\\theta_k}{\\theta_0} g\\left(\\theta_0 e / \\theta_k ; \\theta_0\\right)$? Otherwise, I am not clear why $\\epsilon_k \\stackrel{d}{=}\\left(\\theta_k / \\theta_0\\right) U$ follows.\n\n1. While three distributions are mentioned that fulfill Assumption 2.1., are there more to add to the list? If so, I could not see any discussion on when it makes sense to assume a particular error distribution.\n\n1. The notation is a bit sloppy in some places. In the definition of LiNGAM, it is not assumed that the error means are zero. However, Assumption 2.1 gives the impression that the errors have mean zero. Also, at the beginning of Page 6, the density $\\phi$ is said to have mean zero, however, near the end of page 6, $\\kappa$ is defined with a normal distribution with mean $E[X]$, is that not zero?\n\n1. It would help to standardize the notation of the expectations of the likelihood ratios. In eq.(3) the expectation is w.r.t $f_{kt}(r_{kt})$, while in pages 6 and 7, the expectation is w.r.t. $f(x)$.\n\n1. In the definition of $(\\mathcal{L}-\\tilde{\\mathcal{L}})_t (\\hat{\\pi})$ in Page 7, should it be $ R\\_{\\hat{\\pi}(j)}^{\\hat{\\pi}} $ instead of $ R\\_{\\hat{\\pi}(j)}^{\\pi} $?\n\n1. In Algorithm 2, I think it is more clear to write $S(k, \\mathcal{A}_t;\\mathbf{R})$ instead of just $S(k;\\mathbf{R})$. This is to be consistent with the previous definition of $S$.\n\n\n", "broader_impact_concerns": "No broader impact statement is included in the paper. I believe the authors could try addressing how applying this method in practice could have negative ethical implications when wrongly applied or maliciously.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "extaSi87V8", "reviewer": "Reviewer_5PmC", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1) The paper is generally well written and technical matter is self contained \n\n2) The proposed method allows existing LiNGAM-style methods to scale to larger variable sets when neighborhoods are known or Markov blankets can be estimated\n\n3) The correctness results are rigorous \n\n4) Simulations are comprehensive and results indicate a slight gain in accuracy\n\n\nWeaknesses:\n\n1) While the main contribution of the paper is a causal order estimation procedure that is more computationally efficient than existing approaches, the theoretical complexities of these approaches are not provided\n\n2) The list of assumptions required for the procedure to work is somewhat extensive: linearity, non-Gaussianity, scale-location families, knowledge of neighborhoods\n\n3) In general, the significance of the neighborhood assumption is somewhat glossed over - it's not clear how easy it is to estimate Markov blankets accurately in general and how much errors made during this estimation affect the procedure\n\n4) The intuitive explanation of Theorem 2.5 could be better explained", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a new approach for learning a topological node ordering of a DAG (causal ordering) when causal relations are linear with non-Gaussian noise, e.g. as in the first step of the LiNGAM algorithm, based on likelihood ratios of linear regression residuals. \n\nThe approach is $O(p^2)$ in general and $O(pd)$ when neighborhoods of each node are known or Markov blankets can be estimated accurately and $d$ is the largest neighborhood size.\n\nThe procedure is proven to be correct when noise terms belong to a non-Gaussian scale-location family density.\n\nExperiments show the accuracy of the approach is favorable and more computationally efficient for smaller networks and it can scale to larger networks than can be handled by existing comparable approaches.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Most important requested changes:\n\n1) The paper would be greatly improved by providing a thorough discussion of the theoretical complexities of related methods for estimating causal orderings (and maybe adding a table listing each)\n\n2) The paragraph preceding Theorem 2.5 does not provide an intuitive explanation for why selecting the most non-Gaussian ratio is correct. This becomes clear in the following section, but it would be helpful to explain the general idea at this stage in the paper.\n\n\nOther possible changes that may improve the paper:\n\n1) Better understanding the robustness of the algorithm, compared to existing approaches. when estimated Markov blanket sets are slightly incorrect would improve the paper\n\n2) It's not clear whether the 'order estimation error' defined in (11) is new or has been used in existing approaches. It might improve the paper to also include a more well known ranking accuracy metric like Spearman correlation or at least provide more discussion for the justification of the proposed metric. Specifically, the proposed metric may overpenalize errors made at the beginning of the procedure when the ranking is otherwise correct compared to Spearman or a related metric.", "broader_impact_concerns": "No concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors have substantially revised the paper in response to reviewer comments. They have listed a few changes that are still pending. They should also review the claim of robustness to misspecification and identify ways in which it can be better aligned with their evidence if possible.
Decision: Accept (poster)