Year: 2021

Abstract:
Discovering the causal relationship via recovering the directed acyclic graph (DAG) structure from the observed data is a challenging combinatorial problem. This paper proposes an extremely fast, easy to implement, and high-performance DAG structure recovering algorithm. The algorithm is based on the Cholesky factorization of the covariance/precision matrix. The time complexity of the algorithm is O(pn + p), where p and n are the numbers of nodes and samples, respectively. Under proper assumptions, we show that our algorithm takes O(log(p)) or O(p) samples to exactly recover the DAG structure under proper assumptions. In both time and sample complexities, our algorithm is better than previous algorithms. On synthetic and real-world data sets, our algorithm is significantly faster than previous methods and achieves state-of-the-art performance.

heading: 1 INTRODUCTION
text: As Schelling had said: “The whole world is thoroughly to caught in reason, but the question is: how did it get caught in the network of reason in the first place?” (Kuhn, 1942; Žižek & von Schelling, 1997), people found that learning the causal inferences between the variables is a fundamental problem and has many applications in biology, machine learning, medicine, and economics. The problem usually is considered as finding a directed acyclic graph (DAG) from an observational joint distribution. Unfortunately, learning the DAG structure from the observations is proved to be an NP-hard problem (Chickering, 1995; Chickering et al., 2004).
The problem is generally formulated as the structural equation model (SEM), where the variable of a child node is a function of its parents with additional noises. Depending on the types of functions (linear or non-linear) and noises (Gaussian, Gumbel, etc.), there are several SEM families, e.g., Spirtes et al. (2000); Geiger & Heckerman (1994); Shimizu et al. (2006). In general, the graph can be identified from the joint distribution only up to Markov equivalence classes. Zhang & Hyvarinen (2012); Peters et al. (2014); Peters & Bühlmann (2014); Gao et al. (2020) propose several SEM forms that make the graph fully identifiable from the observed data.
Various algorithms had been proposed to deal with the problem. Search-based algorithms (Chickering, 2002; Friedman & Koller, 2003; Ramsey et al., 2017; Tsamardinos et al., 2006; Aragam & Zhou, 2015; Teyssier & Koller, 2005; Ye et al., 2019; Lv et al., 2021) generally adopt a score (e.g., BIC (Peters et al., 2014) score, Cholesky score (Ye et al., 2019), remove-fill score (Squires et al., 2020)) to measure the fitness of different graphs over data and then search over the legal DAG space to find the structure that achieves the highest score. However, exhaustive search over the legal DAG space is infeasible when p is large (e.g., there are 4.1e18 DAGs for p = 10 (Sloane et al., 2003)). Those algorithms go in quest of a trade-off between the performance and the time complexity.
Since Zheng et al. (2018) proposed an approach that converts the traditional combinatorial optimization problem into a continuous program, many methods (Yu et al., 2019; Lee et al., 2019; Ng et al., 2019a;b; Zheng et al., 2020; Lachapelle et al., 2020; Squires et al., 2020; Zhu et al., 2021) have been proposed. Those algorithms formalize the problem as a data reconstruction task with various differentiable constraints on the DAG adjacent matrix and solve it via the augmented Lagrangian method. These algorithms are able to utilize neural networks to approximate the complicated relations between the features in the observed data and achieve good performances. Recently, reinforcement learning based algorithms (Zhu et al., 2020; Wang et al., 2021) also improved the performance by exploring the possible DAG structure candidates. The algorithms update the parameters of the model
via policy gradient as long as it explored a better DAG structure with a higher reward which measures how well an explored structure meets the requirement of DAG and the observed data.
Topology order search algorithms (TOSA) (Ghoshal & Honorio, 2017; 2018; Chen et al., 2019; Gao et al., 2020; Park, 2020) decompose the DAG learning problem into two phases: (i) Topology order learning via conditional variance of the observed data; (ii) Graph estimation depends on the learned topology order. Those algorithms reduce the computation complexity into polynomial time and are guaranteed to recover the DAG structure under some identifiable assumptions. Our method in this paper is also a topology order search algorithm and it merges the two phases in TOSA into one. In each iteration, it attempts to find a child or a contemporary of the current node. Meanwhile, it also determines the corresponding column vector of the adjacent matrix. The mergence brings three main differences: First, the topology order in TOSA is recovered purely based on the conditional variance of the observed data, whereas our method may also take the sparsity of the adjacent matrix into account; Second, the graph LASSO methods, which are commonly adopted to estimate the graph in the second phase in TOSA, encourage the sparsity of the precision matrix, whereas our method is able to encourage the sparsity of the adjacent matrix; Third, the time complexity is reduced significantly. To be specific, the time complexity of our algorithm is O(p2n + p3), while the fastest algorithm before is O(p5n) (Park, 2020; Gao et al., 2020). Here p and n are the numbers of nodes and samples, respectively. In addition, under proper assumptions, we show that our algorithm takes O(log(p)) or O(p) samples to exactly recover the DAG structure. Compared with previous TOSA algorithms, the sample complexity of our method is much better. Experimental results on synthetic data sets, proteins data sets, and knowledge base data set demonstrate the efficiency and effectiveness of our algorithm. For synthetic data sets, compared with previous baselines, our algorithm improves the performance with a significant margin and at least tens or hundreds of times faster. For the proteins data set, we achieve state-of-the-art performance. For the knowledge base data set, we can observe many reasonable structures of the discovered DAG. Our code is uploaded as supplementary material and will be open-sourced upon the acceptance of this paper.
The rest of this paper is organized as follows. In Section 2, we present our algorithm together with the theoretical analysis. In Section 3, numerical results on synthetic data sets, proteins data set, and knowledge base data set are given. Finally, the paper is concluded in Section 4.
Notations. The symbol ‖ · ‖ stands for the Euclid norm of a vector or the spectral norm of a matrix. For a vector x = [x1,x2, . . . ,xp] ∈ Rp, ‖ · ‖1 stands for the `1-norm, i.e., ‖x1‖ = ∑p i=1 |xi|. For a matrix X = [Xij ] ∈ Rm×n, ‖ · ‖2,∞ stands for the two-to-infinity norm, i.e., ‖X‖2,∞ = max1≤i≤m ‖Xi,:‖; ‖ · ‖max stands for the max norm, ‖X‖max = maxi,j |Xij |.

heading: 2 CAUSAL DISCOVERY VIA CHOLESKY FACTORIZATION (CDCF)
text: In this section, we first present some preliminaries on DAG, then motivating our algorithm. Next, the detailed algorithm and theoretical guarantees for the exact recovery of the algorithm are given.

heading: 2.1 PRELIMINARIES
text: We assume the observed data is entailed by a DAG G = (p, V,E), where p is the number of nodes, V = {v1, ..., vp} and E = {(vi, vj)|i, j ∈ {1, ...p}} represent the set of nodes and edges, respectively. Each node vi is corresponding to a random variable Xi. The observed data matrix X = [x1, ...,xp] ∈ Rn×p where xi is consisting of n i.i.d observations of the random variable Xi. The joint distribution of X is P (X) = ∏p i=1 P (Xi|PaG(Xi)), where PaG(Xi) := {Xj |(vi, vj) ∈ E} is the parents of node Xi.
Given X , we seek to recover the latent DAG topology structure for the joint probability distribution (Hoyer et al., 2008; Peters et al., 2017). Generally, X is modeled via a structural equation model (SEM) with the form
Xi = fi(PaG(Xi)) +Ni, (i = 1, ..., p),
where fi is an arbitrary function representing the relation between Xi and its parents, Ni is the jointly independent noise variable.
In this paper, we focus on the linear SEM defined by Xi = Xwi +Ni, (i = 1, ..., p),
where wi ∈ Rp is a weighted column vector. Let W = [w1, . . . ,wp] ∈ Rp×p be the weighted adjacency matrix, N = [n1, . . . ,np] ∈ Rn×p be an additive independent noise matrix, where ni is n i.i.d observations following the noise variable Ni. Then the linear SEM model can be formulated as
X = XW +N . (1)
We assume the noise deviation of the child variable is approximately larger than that of its parents (see Theorem 2.1 for details). Following this assumption, a classical identifiable form of SEM is the linear-Gaussian SEM, where all Ni are i.i.d. and homoscedastic (Peters & Bühlmann, 2014).

heading: 2.2 ALGORITHM MOTIVATION
text: As proposed in McKay et al. (2003); Nicholson (1975), a graph is DAG if and only if the corresponding weighted adjacent matrix W can be decomposed into
W = PTPT, (2)
where P is a permutation matrix, T is a strict upper triangular matrix, i.e., Tij = 0 for all i ≤ j.
We denote the scaled permuted data matrix as X̂ = 1√ n XP , the scaled permuted noise matrix as N̂ = 1√ n NP , and the permutation order [i∗1, i ∗ 2 . . . , i ∗ p] = [1, 2, . . . , p]P . We can rewrite (1) as
X̂ = X̂T + N̂ .
Then it follows that X̂ = N̂(I − T )−1. (3)
Let E(N̂TN̂) = Σ̂2∗ = Σ̂TΣ̂, (4)
where Σ̂2∗ is the covariance matrix of the noise variables, Σ̂ is upper triangular – the Cholesky factor of Σ̂2∗. Let the diagonal entries of Σ̂ be σ
2 i∗1 , σ2i∗2 , . . . , σ 2 i∗p . We know that σ2i∗k is the conditional variance of Ni∗k .
Now using (3) and (4), we have the covariance matrix of the permuted data:
Ĉ∗ = E(X̂TX̂) = (I − T )−TE(N̂TN̂)(I − T )−1 = (I − T )−TΣ̂TΣ̂(I − T )−1. (5)
Let L = (I − T )−TΣ̂T, then Ĉ∗ = LLT , which is the Cholesky factorization of the covariance matrix Ĉ∗ since L is lower triangular. Furthermore, we can see that the diagonal entries of L are the same as that of Σ̂, i.e., Lkk = σi∗k , the conditional variances of Xi∗k and Ni∗k are the same.
The task becomes to find the permutation i∗ = [i∗1, i ∗ 2, . . . , i ∗ p] and an upper triangular matrix U such that U−TU−1 is a good approximation of the empirical estimation of the permuted covariance matrix Ĉ = 1nX T :,i∗X:,i∗ , and U satisfies some additional constraints, such as the sparsity, etc.

heading: 2.3 ALGORITHM
text: We iteratively find the permutation i and calculate U via the Cholesky factorization. Assume that ik−1 = [i1, . . . , ik−1] and Uk−1 = U1:k−1,1:k−1 are settled, and we have
C1:k−1,1:k−1 = 1
n XT:,ik−1X:,ik−1 + λI = U −T k−1U −1 k−1, (6)
where λ > 0 is a diagonal augmentation parameter which we will give detailed discussion latter. Next, we show how to find ik and the last column of Uk.
For the time being, let us assume ik is known, we show how to compute the last column of Uk. Let U−1k =
[ U−1k−1 yk
0 αk ] , then[
U−1k−1 yk 0 αk
]T [ U−1k−1 yk
0 αk
] = [ U−Tk−1U −1 k−1 U −T k−1yk
yTkU −1 k−1 α 2 k+‖yk‖ 2
] = 1
n
[ XT:,ik−1 X:,ik−1+λI X T :,ik−1 X:,ik
XT:,ik X:,ik−1 ‖X:,ik‖ 2+λ
] ,
Algorithm 1 Causal Discovery via Cholesky Factorization (CDCF) 1: input: Data matrix X ∈ Rn×p, Truncate Threshold ω > 0, and tuning parameter γ. 2: output: Adjacent Matrix A. 3: Set i = [1, 2, . . . , p], R = ‖X‖22,∞ and λ = γ log p n R;
4: Set ` = argmin {‖X:,i1‖, ‖X:,i2‖, . . . , ‖X:,ip‖}; 5: Exchange i1 and i` in i; Set U1 =
√ n
‖X:,i`‖2+λ ;
6: for k = 2, 3, . . . , p do 7: for j = k, k + 1, . . . , p do 8: yj = 1 nU T k−1X T :,ik−1 X:,ij ;
9: αj = √ 1 n‖X:,ij‖2 + λ− ‖yj‖2;
10: end for 11: (V) ` = argmink≤j≤p α2j ;
(S) ` = argmink≤j≤p ‖Uk−1yj‖1; (VS) ` = argmink≤j≤p ‖Uk−1yj‖1 √∣∣α2j − 1k−1 ∑k−1h=1 1[Uk−1]2hh ∣∣; 12: Exchange ik and i` in i;
13: Set Uk = [ Uk−1 − 1α`Uk−1y`
0 1α`
] ;
14: end for 15: return A = [TRIU(TRUNCATE(Up, ω))]REVERSE(i),REVERSE(i).
where the last equality dues to (6). It follows that
yk = 1
n UTk−1X T :,ik−1 X:,ik , αk =
√ 1
n ‖X:,ik‖2 + λ− ‖yk‖2. (7)
And direct calculation gives rise to Uk = [ U−1k−1 yk
0 αk
]−1 = [ Uk−1 − 1αkUk−1yk
0 1αk
] . (8)
By (8), once ik is settled, we can obtain the last column of Uk. Our task remains to select ik from {1, . . . , p} \ {i1, . . . , ik−1}. There are several ways to accomplish this task. We propose three criteria to select ik. First, we need to compute αj and yj by (7) for all possible j (ij ∈ {1, . . . , p} \ {i1, . . . , ik−1}). Then we select ik according to one of the following criteria:
(V) ik = argmink≤j≤p α2j . Under the assumption that the noise variance of the child variable is approximately larger than that of its parents, it is reasonable/natural to select the index that has the lowest estimation of the noise variance. This criterion is guaranteed to find the correct permutation i∗ with high probability, which is shown in Section 2.4.
(S) ik = argmink≤j≤p ‖Uk−1yj‖1. Using (3) and (6), we know that Up intends to estimate (I − T )Σ̂−1. When the adjacent matrix T is sparse and the noise variables are independent (i.e., Σ̂ is diagonal), we would like to select the index that leading to the most sparse column of Uk. This criterion is especially useful when the number of samples is small, see Tables B.1, B.2 and B.3 in appendix.
(VS) ik = argmink≤j≤p ‖Uk−1yj‖1 √∣∣α2j − 1k−1 ∑k−1h=1 1[Uk−1]2hh ∣∣. We empirically combine
criterion (V) and criterion (S) together to take both aspects (variance and sparsity) into account. Numerically, we found that this criterion achieves the best performance in real-world data.
The diagonal augmentation trick in (6) is commonly used for an invertible and good conditioned estimation of the covariance matrix (see e.g., (Ledoit & Wolf, 2004)). Such a trick not only ensures that our algorithm does not break down due to the singularity of the sample covariance matrix, but also stabilizes the Cholesky factorization, especially when the sample is insufficient. In addition, by setting λ = O( log pn ), the error bound between the population covariance matrix and the augmented sample covariance matrix does not become worse (see Lemma ?? in the appendix). This trick
significantly improves the ability to recover the DAG, especially when the samples are insufficient, see Tables B.4, B.5 and B.6 in appendix.
The detailed algorithm is summarized in Algorithm 1. Some comments and implementation details follow. Line 4, we select the very initial value ` = argmin {‖X:,i1‖, ‖X:,i2‖, . . . , ‖X:,ip‖}. Line 5, we exchange i1 and i` in i and calculate U1 = √ n
‖X:,i`‖2+λ . Lines 6 to 14, we iteratively calculate
Uk and update permutation order i until all the indices are settled. Line 15, we truncate U , take its strict upper triangular part (denoted by “TRIU”) and re-permute the predicted adjacent matrix back to the original order according to the permutation order i. Specifically, the truncation is done column-wisely. By (8), the value of [Up]:,k is inversely proportional to αk. So, for column k, we set ωk =
ω αk , and do the truncation: [Up]ik is set to zero if |[Up]ik| < ωk. On output, node i connects to node j in G if |Aij | > 0.
Time Complexity Note that we do not have to re-calculate the matrix multiplication of XT:,ik−1X:,ij in line 8 since we can calculate C at the cost of O(p
2n) at first. Besides, at step k, we have already calculate UTk−2X T :,ik−1
X:,ij at previous step, we only need to calculate the last entry of yj , which is the inner product between two k dimensional vectors, at the cost ofO(p) in worst case. Overall, the time complexity of CDCF is O(p3 + p2n). When n > p, the complexity becomes O(p2n), which is equivalent to the complexity of calculating the covariance matrix. Additionally, the inner loop (lines 7 to 10) of CDCF can be done in parallel, which makes the algorithm friendly to run on GPU and suitable for large scale calculations.

heading: 2.4 EXACT DAG STRUCTURE RECOVERY
text: The following theorem tells that our algorithm is able to recover the DAG exactly with high probability under proper assumptions.
Theorem 2.1 Let x ∈ Rp be a zero-mean random vector, C = E(xxT) ∈ Rp×p be the covariance matrix. Let x1, . . . ,xn be n independent samples, Ĉ = 1n ∑n k=1 xkx T k be the sample covariance estimator. Assume ‖C − Ĉ‖ ≤ for some > 0. Denote Ĉλ = Ĉ + λI , where λ = O( ) ≥ 0 is a parameter. Let the Cholesky factorizations of C = ExxT and Ĉλ be C = LLT and Ĉλ = L̂L̂T, respectively, where L and L̂ are both lower triangular. For the linear SEM model (1), assume (2) and (4), and for k ∈ PaG(j), δ = infk∈PaG(j) δjk > 0, where
δjk = σ 2 i∗j + ‖Σ̂n[(I − T )−1]k:j−1,k‖2 − σ2i∗k .
If δ ≥ 4( + λ) and ‖L−1‖2( + λ) < 34 , then CDCF-V is able to recover P exactly. In addition, it holds that
‖TRIU(Up)− T ‖max ≤ 4‖Σ̂−1∗ (I − T )T‖22,∞‖(I − T )Σ̂−T∗ ‖2,∞( + λ), where TRIU(Up) stands for the strictly upper triangular part of Up, Up is the output of outer loop of Algorithm 1 with criterion (V).
we know that when T is sparse, we may recover its topology structure by truncating Up.
Proposition 1 Let Ni,: be independent bounded, or sub-Gaussian, or regular polynomial-tail, then for n > N( ), it holds ‖Ĉxx −Cxx‖ ≤ , w.h.p. Specifically,
N( ) ≥ C1 log p (‖(I − T )−1‖2‖Cnn‖ )2 , for bounded class;
N( ) ≥ C2 p (‖(I − T )−1‖2‖Cnn‖ )2 , for the sub-Gaussian class;
N( ) ≥ C3 p (‖(I − T )−1‖2‖Cnn‖ )2(1+r−1) , for the regular polynomial tail class.
The proofs of are provided in the Appendix A. This theorem and proposition also indicates the sample complexity of our algorithm is O(p). This sample complexity is better than the sample complexities of previous methods, see Table 2.1 for a detailed comparison.

heading: 3 EXPERIMENTS
text: In this section, we apply our algorithm to synthetic data sets, proteins data set and knowledge base data set, respectively, to illustrate the efficiency and effectiveness of our algorithm.

heading: 3.1 LINEAR SEM
text: We evaluate the proposed methods on simulated graphs from two well-known ensembles of random graph types: Erdös–Rényi (ER) (Gilbert, 1959) and Scale-free (SF) (Barabási & Albert, 1999). The average edge number per node is denoted after the graph type. For example, ER2 represents two edges per node on average. After the graph structure is settled, we assign uniformly random edge weights to obtain a weight matrix W . We generate the observation data X from the linear SEM with three noise distributions: Gaussian, Gumbel, Exponential.
We chose our baseline methods as NOTEARS (Zheng et al., 2018), DAG-GNN (Yu et al., 2019), CORL (Wang et al., 2021), NPVAR (Gao et al., 2020), and EQVAR (Chen et al., 2019). Other methods such as PC algorithm (Spirtes et al., 2000), LiNGAM (Shimizu et al., 2006), FGS (Ramsey et al., 2017), MMHC (Tsamardinos et al., 2006), L1OBS (Schmidt et al., 2007), CAM (Bühlmann et al., 2013), RL-BIC2 (Zhu et al., 2020), A*LASSO (Xiang & Kim, 2013), LISTEN (Ghoshal & Honorio, 2018), US (Park, 2020) perform worse than or approximately equal to the selected baselines, and the results can be found in the corresponding papers.
Table 3.1 presents the structural Hamming distance (SHD) of baseline methods and our method on 3000 samples (n = 3000). Nodes number p is noted in the first column. Graph type and edge level are noted in the second column. We only report the SHD of different algorithms due to page limitation, and we find that other metrics such as true positive rate (TPR), false discovery rate (FDR), false positive rate (FPR), and F1 score have the similar comparative performance with SHD. We also test bottom-up EQVAR which is equivalent to LISTEN, the result is worse than top-down EQVAR (EV-TD) in this synthesis experiment, so we do not include the result in the table. For p = 1000 graphs, we only report the result of EV-TD and CDCF since other algorithms spend too much time (longer than a week) to recover a DAG. We test our algorithms with different variations according to criteria (V, S, VS) introduced in Section 2.3, and with diagonal augmentation trick noted by a “+” as postfix. For example, "CDCF-V" means CDCF with V criterion and λ = 0, and "CDCF-V+" means CDCF with V criterion and λ = O( log pn ). The implementation details are in the Appendix B. We report the result of CDCF-V+ here, and the results of other CDCF variations can be found in Appendix Table B.4. We run our methods on ten randomly generated graphs and report the mean and variance in the table. Figure 3.1 plots the SHD results tested on 100 nodes graph recovering from different sample sizes. We choose EV-TD and high dimension top down (EV-HTD) as baselines when p > n and p ≤ n, respectively. We can see from the results, CDCF-V+ achieves significantly better performance comparing with previous baselines.
Table 3.2 shows the running time which is tested on a 2.3 GHz single Intel Core i5 CPU. Besides, parallel calculation of the matrix multiplication on GPU makes the algorithm even faster. Recovering 5000 and 10000 nodes graph from 3000 samples on an A100 Nvidia GPU is approximately 400 and 2400 seconds, respectively. For comparison, EV-TD costs approximately 100 hours to recover a 1000 nodes DAG from 3000 samples. As illustrated in the table, CDCF is approximately dozens or hundreds of times faster than EV-TD and LISTEN, and tens of thousands times faster than NOTEARS as CDCF does not have to update the parameters with gradients.
Due to the page limitation, further experiments and discussions of the ablation study (Figures B.3 to B.14, Tables B.1 to B.6), choice of λ (Tables B.7 to B.10), and performances on different noise distribution (Figures B.1, B.2) and deviation (Tables B.11, B.12, B.13) are given in Appendix B.

heading: 3.2 PROTEINS DATA SET
text: We consider a bioinformatics data set (Sachs et al., 2005) consisting of continuous measurements of expression levels of proteins and phospholipids in the human immune system cells. This is a widely used data set for research on graphical models, with experimental annotations accepted by the biological research community. Following the previous algorithms setting, we noticed that different previous papers adopted different observations. To included them all, we considered the observational 853 samples from the "CD3, CD28" simulation tested by Teyssier & Koller (2005); Lachapelle et al. (2020); Zhu et al. (2020) and all 7466 samples from nine different simulations tested by Zheng et al. (2018; 2020); Yu et al. (2019).
We report the experimental results on both settings in Table 3.3. The implementation codes of the baselines are introduced in the appendix, and we use the default settings of the hyper-parameters provided in their codes. The evaluate metric is FDR, TPR, FPR, SHD, predicted nodes number (N), precision (P), F1 score. As the recall score is equal to TPR, we do not include it in the table. In both settings, CDCF-VS+ achieves state-of-the-art performance. 1 Several reasons make the recovered graph not exactly the same as the expected one. The ground truth graph suggested by the paper is mixed with directed and indirect edges. Under the settings of SEM, the node "PKA" is quite similar to the leaf nodes since most of its edges are indirect while the ground truth graph notes it as the out edges. Non-linear would not be an impact issue here since NOTEARS and our algorithm both achieve decent results. In the meantime, we do not deny that further extension of our algorithm to non-linear representation would witness an improvement on this data set.

heading: 3.3 KNOWLEDGE BASE DATA SET
text: We test our algorithm on FB15K-237 data set (Toutanova et al., 2015) in which the knowledge is organized as {Subject, Predicate,Object} triplets. The data set has 15K triplets and 237 types of predicates. In this experiment, we only consider the single jump predicate between the entities, which
1For NOTEARS-MLP, Table 3.3 reported the results reproduced by the code provided in Zheng et al. (2020).
have 97 predicates remained. We want to discover the causal relationships between the predicates. We organize the observation data as each sample corresponds to an entity with awareness of the position (Subject or Object), and each variable corresponds to a predicate in this knowledge base.
In Figure 3.2, we give the adjacent weighted matrix of the generated graph and several examples with high confidence (larger than 0.5). In the left figure, the label of the axis notes the first capital letter of the domain of the relations. Some of them are replaced with a dot to save space. The exact domain name and the picture with the full predicate name are provided in the appendix. The domain clusters are denoted in black boxes at the diagonal of the adjacent matrix. The red boxes denoted the cross-domain relations that are worth paying attention to. Consistent with the innateness of human sense, the recovered relationships inside a domain are denser than those across domains. In the cross-domain relations, we found that the predicate in domain "TV" ("T") has many relations with the domain "Film" ("F"), the domain "Broadcast" (last row) have many relations with the domain "Music" ("M"). Several cases of the predicted causal relationships are listed on the right side of Figure 3.2, we can see that the discovered indication relations between predicates are quite reasonable.

heading: 4 CONCLUSION AND FUTURE WORK
text: In this paper, we proposed a topology search algorithm for the DAG structure recovery problem. Our algorithm is better than the existing methods in both time and sample complexities. To be specific, the time complexity of our algorithm isO(p2n+ p3), while the fastest algorithm before isO(p5n) (Park, 2020; Gao et al., 2020), where p and n are the numbers of nodes and samples, respectively. Under different assumptions, our algorithm takes O(log(p)) or O(p) samples to exactly recover the DAG structure. Experimental results on synthetic data sets, proteins data sets, and knowledge base data set demonstrate the efficiency and effectiveness of our algorithm. For synthetic data sets, compared with previous baselines, our algorithm improves the performance with a significant margin and at least tens or hundreds of times faster. For the proteins data set, we achieve state-of-the-art performance. For the knowledge base data set, we can observe many reasonable structures of the discovered DAG.
The proposed algorithm is under the assumption of linear SEM. Generalization of CDCF to nonlinear SEM would be a valuable and important research topic. Learning the representation of the observed data for better structure reconstruction via the CDCF algorithm, which requires the algorithm differentiable, is also an attractive problem. To deal with the extremely large-scale problems, such as millions of nodes, implementing CDCF via sparse matrix storage and calculation on the GPU is a promising way to further improve computational performance.

heading: A PROOF OF THEOREM 2.1
text: In this section, we first give several lemmas, then prove Theorem 2.1.
Lemma A.1 Let x ∈ Rp be a zero-mean random vector, C = E(xxT) ∈ Rp×p be the covariance matrix. Let x1, . . . ,xn be n independent samples, Ĉ = 1n ∑n k=1 xkx T k be the sample covariance estimator. Assume ‖C − Ĉ‖ ≤ for some > 0. Denote Ĉλ = Ĉ + λI , where λ = O( ) ≥ 0 is a parameter. Let the Cholesky factorizations of C = ExxT and Ĉλ be C = LLT and Ĉλ = L̂L̂T, respectively, where L and L̂ are both lower triangular. If ‖L−1‖2( + λ) < 34 , then
|‖Li,:‖2 − ‖[L̂λ]i,:‖2| ≤ + λ = O( ), for 1 ≤ i ≤ p; (9)
|[L−1]ij − [L̂−1]ij | ≤ 4‖L−1‖22,∞‖L−T‖2,∞( + λ) = O( ), for i > j. (10)
Proof. For all 1 ≤ i ≤ p, we have
|‖Li,:‖2 − ‖L̂i,:‖2| = |Cii − [Ĉλ]ii| ≤ ‖C − Ĉλ‖ ≤ ‖C − Ĉ‖+ λ ≤ + λ, (11)
which completes the proof for (9).
Next, we show (10). Let
L−1L̂ = I + F , (I + F )(I + F )T = I +E. (12)
We know that
L̂−1 −L−1 = [(I + F )−1 − I]L−1 = −F (I + F )−1L−1, (13)
E = L−1L̂L̂TL−T − I = L−1(Ĉλ −C)L−T. (14)
Then it follows from (13) that for i > j
|[L−1]ij−[L̂−1]ij | ≤ ‖Fi,1:i−1‖‖[(I+F )−1L−1]:,j‖ ≤ ‖Fi,1:i−1‖‖(I+F )−1‖‖L−T‖2,∞. (15)
First, we give an upper bound for ‖(I + F )−1‖. Using (12), we have (I + F )−T(I + F )−1 = (I +E)−1. It follows
‖(I + F )−1‖ = ‖(I + F )−T(I + F )−1‖ 12 = ‖(I +E)−1‖ 12
≤ 1√ 1− ‖E‖ ≤ 1√ 1− ‖L−1‖2‖Ĉλ −C‖ , (16)
where the last inequality uses (14).
Second, we give upper bound for ‖Fi,1:i−1‖. It follows from the second equality of (12) that
(1 + Fii) 2 + ‖Fi,1:i−1‖2 = 1 +Eii. (17)
Therefore,
‖Fi,1:i−1‖2 ≤ |(1 + Fii)2 − 1|+Eii (a) ≤ L̂ 2 ii −L2ii L2ii +Eii
(b) ≤ + λ L2ii + ‖L−1‖22,∞‖Ĉλ −C‖ (c) ≤ 2‖L−1‖22,∞( + λ), (18)
where (a) uses (12), (b) uses (9) and (14), (c) uses ‖C − Ĉ‖ ≤ . Substituting (18) and (16) into (15), we get
|[L−1]ij − [L̂−1]ij | ≤ 2‖L−1‖22,∞‖L−T‖2,∞ + λ√
1− ‖L−1‖2( + λ) . (19)
The conclusion follows since ‖L−1‖2( + λ) < 34 .
Theorem A.2 Let x ∈ Rp be a zero-mean random vector, C = E(xxT) ∈ Rp×p be the covariance matrix. Let x1, . . . ,xn be n independent samples, Ĉ = 1n ∑n k=1 xkx T k be the sample covariance estimator. Assume ‖C − Ĉ‖ ≤ for some > 0. Denote Ĉλ = Ĉ + λI , where λ = O( ) ≥ 0 is a parameter. Let the Cholesky factorizations of C = ExxT and Ĉλ be C = LLT and Ĉλ = L̂L̂T, respectively, where L and L̂ are both lower triangular. For the linear SEM model (1), assume (2) and (4), and for k ∈ PaG(j), δ = infk∈PaG(j) δjk > 0, where
δjk = σ 2 i∗j + ‖Σ̂n[(I − T )−1]k:j−1,k‖2 − σ2i∗k .
If δ ≥ 4( + λ) and ‖L−1‖2( + λ) < 34 , then CDCF-V is able to recover P exactly. In addition, it holds that
‖TRIU(Up)− T ‖max ≤ 4‖Σ̂−1∗ (I − T )T‖22,∞‖(I − T )Σ̂−T∗ ‖2,∞( + λ),
where TRIU(Up) stands for the strictly upper triangular part of Up, Up is the output of outer loop of Algorithm 1 with criterion (V).
Proof. For SEM model (1), denote Ĉ∗ = E( 1nX̂ TX̂), Σ̂2∗ = E( 1nN̂ TN̂) = Σ̂Tn Σ̂n, we have (5), i.e.,
Ĉ∗ = (I − T )−TΣ̂2∗(I − T )−1 = (I − T )−TΣ̂Tn Σ̂n(I − T )−1. (20)
When the permutation i∗ = [i∗1, . . . , i ∗ p] is exactly recovered, then Up in CDCF-V satisfies
Ĉλ = 1
n XT:,i∗X:,i∗ + λI = U −T p U −1 p . (21)
Denote i∗j = [i ∗ 1, . . . , i ∗ j ] for all j = 1, . . . , p. Consider the kth diagonal entries of (20) and (21). By calculations, we get
[Ĉ∗]kk = [(I − T )−1]T:,kΣ̂Tn Σ̂n[(I − T )−1]:,k = σ2i∗k + ‖uk‖ 2, (22) [Ĉλ]kk = 1
n ‖Xi∗k‖
2 + λ = 1
U2kk + ‖ûk‖2, (23)
where
uk = [Σ̂n]1:k−1,1:k−1(Ik−1 − T1:k−1,1:k−1)−1T1:k−1,k, ûk = 1
n UTk−1X T :,i∗k−1 X:,i∗k . (24)
Using ‖C − Ĉ‖ ≤ , we have
|[Ĉ∗]kk − [Ĉλ]kk| ≤ ‖C − Ĉλ‖ ≤ ‖C − Ĉ‖+ λ ≤ + λ. (25)
By Lemma A.1, we have
|‖uk‖2 − ‖ûk‖2| ≤ + λ. (26)
Using (22), (23), (25) and (26), we get
|σ2i∗k − 1
U2kk | ≤ 2( + λ). (27)
Assume that i∗1, . . . , i ∗ k−1 (k ≥ 1) are all correctly recovered. And without loss of generality, for k ∈ PaG(j), we also assume Tk:j−1,j 6= 0 (otherwise, jth and kth columns are exchangeable, and i forms another equivalence topology order to the same DAG (Sedgewick & Wayne, 2011)). Then we
have for k ∈ PaG(j) that 1
n ‖Xi∗j ‖
2 + λ− ‖[ûj ]1:k−1‖2 (a) = [Ĉ∗]jj + [Ĉλ]jj − [Ĉ∗]jj − ‖[ûj ]1:k−1‖2
(b) ≥ [Ĉ∗]jj − ( + λ)− ‖[uj ]1:k−1‖2 − ( + λ) (c) = σi∗j + ‖[uj ]k:j−1‖
2 − 2( + λ) (d) ≥ σi∗k + δ − 2( + λ) (e) = [Ĉ∗]kk − ‖uk‖2 + δ − 2( + λ) (f) ≥ [Ĉλ]kk − ‖ûk‖2 + δ − 4( + λ) (g) = 1
n ‖Xi∗k‖ 2 + λ− ‖ûk‖2 + δ − 4( + λ),
where (a) uses (23), (b) and (f) uses (25) and Lemma A.1, (c) uses (22), (d) dues to the assumption σi∗j ≥ σi∗k for k ∈ PaG(j), (e) uses (22), (g) uses (23). Therefore, using δ > 4( + λ), we have
1 n ‖Xi∗j ‖ 2 + λ− ‖[ûj ]1:k−1‖2 > 1 n ‖Xi∗k‖ 2 + λ− ‖ûk‖2,
which implies that i∗k can be correctly recovered. So, overall speaking, CDCF-V is able to recover the permutation P .
The upper bound for ‖TRIU(Up)− T ‖max follows from Lemma A.1. The proof is completed.
Proposition 2 Let Ni,: be independent bounded, or sub-Gaussian, 2 or regular polynomial-tail, 3 then for n > N( ), it holds ‖Ĉxx −Cxx‖ ≤ , w.h.p. Specifically,
N( ) ≥ C1 log p (‖(I − T )−1‖2‖Cnn‖ )2 , for bounded class;
N( ) ≥ C2 p (‖(I − T )−1‖2‖Cnn‖ )2 , for the sub-Gaussian class;
N( ) ≥ C3 p (‖(I − T )−1‖2‖Cnn‖ )2(1+r−1) , for the regular polynomial tail class.
Proof. For SEM model (1), we have
‖Ĉxx−Cxx‖ ≤ ‖(I−T )−1‖2‖Ĉnn−Cnn‖ ≤ ‖(I−T )−1‖2‖Cnn‖‖C − 12 nn ĈnnC − 12 nn −I‖, (28)
where Cxx = ExxT, Cnn = EnnT are the covariance matrices for x and n, respectively, Ĉxx, Ĉnn are the sample covariance matrices for x and n, respectively. The three results listed above follow from Corollary 5.52, Theorem 5.39 in Vershynin (2010), Theorem 1.1 in Srivastava & Vershynin (2013), respectively.
2A random vector z is isotropic and sub-Gaussian if EzzT = I and and there exists constant C > 0 such that P(|vTz| > t) ≤ exp(−Ct2) for any unit vector v. Here by “Ni,: is sub-Gaussian” we mean that C − 1 2 nn N T i,: is an isotropic and sub-Gaussian random vector. 3A random vector z is isotropic and regular polynomial-tail if EzzT = I and there exist constants r > 1, C > 0 such that P(‖V z‖2 > t) ≤ Ct−1−r for any orthogonal projection V and any t > C · rank(V ). Here by “Ni,: is regular polynomial-tail” we mean that C − 1 2 nn N T i,: is an isotropic and regular polynomial-tail random vector.

heading: B ADDITIONAL EXPERIMENTS
text: Here we provide implementation details and additional experiment results.
Figures B.1, B.2 provide the results of Gumbel and Exponential noises, respectively. As we can see from the result, our algorithm still performs better than Eqvar method in different noise types.
Tables B.1, B.2 , B.3, B.4, B.5, B.6 give results on 100 nodes over different sample sizes and variances of our CDCF methods. As noted in Algorithm 1, we have V, S, VS as different criteria to select the current column, "+" representing the sample covariance matrix augmented with the scalar matrix log p n I . The truncation threshold on column i is ωi = 3.5/αi, where αi is the diagonal value of the Cholesky factor. According to the results, the algorithm "V+" achieves the best performance as the sample size is relatively large. When the sample size is small, the criterion according to sparsity shows very effective performance improvement. We also test different choices over λ = β log pn , β ∈ {0.0, 1.0..., 9.0}, the result is given in Table B.7, B.8, B.9, B.10. Empirically, β ∈ {1.0, 2.0} achieves better results. In practice, one can sample a relatively small and labeled sub-graph of the DAG to test the hyper-parameter setting then apply to large unlabeled the DAG graph.
To test the performance limitation of our methods, we provide the results of SHD on different sample number and node number in Figures B.3 to B.14 where the x-axis represents the sample number (in thousand), the y-axis denotes the node number, the color represents the value of log2(SHD + 1) (the brighter the better). We provide the figures for CDCF-V+, CDCF-S+, and CDCF-VS+ on variances graph and noise types. The figures are drawn on the mean results over ten random seeds. The figures show that the graph can be exactly recovered on 800 nodes at approximately 6000 samples. Comparing CDCF-V+ with CDCF-S+, we find that criterion (S) damages the performance when the sample number is relatively large. When sample number ∈ {1500, 3000} and node number ∈ {400, 800}, CDCF-S+ achieves better performance. Such trend can also be demonstrated in Tables B.1, B.2, B.3. CDCF-VS+ alleviates the poor performance of CDCF-S when the data is sufficient and achieves good performance on real-world data set.
We also test the performance on linear SEM with monotonously increased noise variance. Concretely, assume the topology order is i = {i1, ..., ip}, we set the noise variance of node k as σk = 1 + ik/p. We test the results on Gaussian, Gumbel, and Exponential noise with monotonous noise variance. The results are reported in Tables B.11, B.12 and B.13. As the results indicated, even with different noise levels, our algorithms achieve good performance and are able to exactly recover the DAG structure when the data is sufficient.
In the result for knowledge base data set, the axis labels of Figure 3.2 are ‘Film’, ‘People’, ‘Location’, ‘Music’, ‘Education’, ‘Tv’, ‘Medicine’, ‘Sports’, ‘Olympics’, ‘Award’, ‘Time’, ‘Organization’, ‘Language’, ‘MediaCommon’, ‘Influence’, ‘Dataworld’, ‘Business’, ‘Broadcast’ from left to right for x-axis and top to bottom for y-axis, respectively. The adjacent matrix plotted here is re-permuted to make the relations in the same domain close to each other. We keep the adjacent matrix inside a domain an upper triangular matrix. Such typology is equivalent to the generated matrix with the original order.
Baseline Implementations The baselines are implemented via the codes provided from the following links:
• NOTEARS, NOTEARS-MLP: https://github.com/xunzheng/notears • NPVAR: https://github.com/MingGao97/NPVAR • EQVAR, LISTEN: https://github.com/WY-Chen/EqVarDAG • CORL: https://github.com/huawei-noah/trustworthyAI/tree/master/gcastle • DAG-GNN: https://github.com/fishmoon1234/DAG-GNN

Review summaries:
This paper proposes a method for recovering the causal graph of additive linear models from purely observational data, under some an identifiability assumption, that seems to be related to the forward step-size assumption of [1]. Their algorithm is based on iteratively identifying a root of the causal graph based on its conditional variance. Once a topological order is learned, the graph is constructed by thresholding the Cholesky factor of the permuted precision matrix. The proposed algorithm is then tested and compared on both synthetic and real-world data.

The paper works on causal discovery in the linear Gaussian case, on which the identifiability is based on (Peters & Bühlmann, 2014). The proposed method is based on Cholesky factorization and has better efficiency/time-complexity performance than the related state-of-art methods. Moreover, it also provides a theoretical analysis of the resulted graph, which is appreciated. The experiments can support the claims.

This paper develops a new algorithm for learning linear structural equation models using cholesky factorization. This paper explains that the proposed algorithm is consistent in high dimensional settings and computational feasible. This paper thoroughly discusses the recent studies of learning linear SEMs and provide a clear motivation. Furthermore, this paper provides a lot of numerical experiments to support its theoretical findings.