Learning causal graphs using variable grouping according to ancestral relationship

Let $X=(x_{1},\ldots,x_{p})^{\top}$ be a $p$-dimensional random vector. In the following, we identify $X$ with the variable set. Assume that the causal relationship among the variables is acyclic and linear. LiNGAM [3] is expressed by

where the error term $e=(e_{1},\ldots,e_{p})^{\top}$ is assumed to be independently distributed as continuous non-Gaussian distributions. $B=\{b_{ij}\}$ is a $p\times p$ coefficient matrix. $b_{ij}$ represents the direct causal effect from $x_{j}$ to $x_{i}$. $b_{ij}=0$ indicates the absence of the direct causal effect from $x_{j}$ to $x_{i}$. We note that it is possible to transform the matrix $B$ into a strictly lower triangular matrix by permuting the rows and columns when a DAG defines the causal relationship between $X$.

The model (1) is rewritten by

where $I$ denotes the $p\times p$ identity matrix. Noting that (2) is equivalent to the independent component analysis (ICA) model, Shimizu et al. [3] showed that $B$ is identifiable and proposed an algorithm for estimating $B$. The algorithm is known as ICA-LiNGAM.

However, as the number of variables increases, ICA-LiNGAM is more likely to converge to a locally optimal solution, which reduces the accuracy of causal DAG estimation. To overcome this problem, Shimizu et al. [4] proposed another algorithm for estimating $B$ named DirectLiNGAM. When linear regression analyses are conducted in two ways for each pair of variables, with one of them as the dependent variable and the other as the independent variable, if the independent variable and the residuals are mutually independent in only one of the models, the independent variable in that model will precede the dependent variable in causal order. After determining the causal order of all variables, DirectLiNGAM estimates the set of parents of each variable by sparse estimating a linear regression model with each variable as the dependent variable and all variables preceding it in the causal order as independent variables. This approach ensures no edges violate the estimated causal order. DirectLiNGAM is also less accurate when the number of variables is large relative to the sample size. DirectLiNGAM tends to output redundant edges even when the sample size is large.

The computational cost of DirectLiNGAM is higher than that of ICA-LiNGAM. The time complexity of DirectLiNGAM is $O(np^{3}M^{2}+p^{4}M^{3})$, where $n$ is the sample size, $p$ is the number of variables, and $M$ is the maximal rank found by the low-rank decomposition used in the kernel-based independence measure [4]. Since DirectLiNGAM requires estimating a linear regression model, DirectLiNGAM is not feasible when $n$ is smaller than $p$.

Cai et al. [1] proposed a divide-and-conquer approach called the scalable causation discovery algorithm (SADA) to enhance the scalability of causal DAG learning algorithms. SADA first splits the variable set $X$ into two or more subsets. Then, a causal DAG learning algorithm such as DirectLiNGAM is applied to each group of variables, and finally, all the results are merged to estimate the entire causal DAG. Because SADA applies the causal DAG learning algorithm to each group with a smaller number of variables, if the variables are grouped correctly, SADA is expected to improve the estimation accuracy when the sample size is small relative to the number of the entire variable set. SADA is feasible even if the sample size is smaller than the number of variables as long as it is larger than the number of variables in each group.

The procedure for grouping $X$ into two subsets by SADA is as follows. SADA first randomly selects two variables $x_{i},x_{j}\in X$ satisfying $x_{i}\mathop{\perp\!\!\!\!\perp}x_{j}\mid(X\setminus\{x_{i},x_{j}\})$ and finds the smallest $\hat{X}\subseteq X\backslash\{x_{i},x_{j}\}$ with respect to the inclusion relation satisfying $x_{i}\mathop{\perp\!\!\!\!\perp}x_{j}\mid\hat{X}$. The disjoint subsets $X_{1}$, $X_{2}$, and $C$ are computed according to the following procedure.

1. 1.

   Initialize with $X_{1}=\{x_{i}\}$, $X_{2}=\{x_{j}\}$, $C=\{\hat{X}\}$.

2. 2.

   For all $w\in X\setminus(X_{1}\cup X_{2}\cup C)$

   • (a)

     if $w\mathop{\perp\!\!\!\!\perp}X_{2}\mid\hat{C}$ for $\exists\hat{C}\subseteq C$,
     then $X_{1}\leftarrow X_{1}\cup\{w\}$

   • (b)

     if $w\mathop{\perp\!\!\!\!\perp}X_{1}\mid\hat{C}$ for $\exists\hat{C}\subseteq C$,
     then $X_{2}\leftarrow X_{2}\cup\{w\}$

   • (c)

     else $C\leftarrow C\cup\{w\}$

3. 3.

   For all $s\in C$

   • (a)

     if $s\mathop{\perp\!\!\!\!\perp}X_{2}\mid\hat{C}$ for $\exists\hat{C}\subseteq C\setminus\{s\}$,
     then $X_{1}\leftarrow X_{1}\cup\{s\}$ and $C\leftarrow C\setminus\{s\}$

   • (b)

     if $s\mathop{\perp\!\!\!\!\perp}X_{1}\mid\hat{C}$ for $\exists\hat{C}\subseteq C\setminus\{s\}$,
     then $X_{2}\leftarrow X_{2}\cup\{s\}$ and $C\leftarrow C\setminus\{s\}$

4. 4.

   Return $X_{1}$, $X_{2}$ and $C$

In step 2, the intermediate set $C$ tends to be large, and step 3 downsizes $C$ as much as possible.

The output of SADA is two subsets $V_{1}=X_{1}\cup C$ and $V_{2}=X_{2}\cup C$. This algorithm can be repeated recursively to group variables into smaller subsets. In the implementation, the lower bound of the number of variables in each group is set to $\theta$, and SADA is applied to $V_{1}$ or $V_{2}$ to create even smaller groups only when $|V_{1}|>\theta$ or $|V_{2}|>\theta$. If the number of variables is less than or equal to $\theta$, or if we cannot find $x_{i},x_{j}\in V_{k}$ satisfying $x_{i}\mathop{\perp\!\!\!\!\perp}x_{j}\mid(V_{k}\setminus\{x_{i},x_{j}\})$, $k=1,2$, no further grouping is made.

We illustrate the SADA procedure to split $X$ into two subsets where the true causal DAG is the one in Fig 1. Table I shows the process of grouping variables $x_{1},\ldots,x_{9}$ into two subsets. Since both $x_{3}\mathop{\perp\!\!\!\!\perp}x_{9}\mid\{x_{1},x_{2},x_{4}\ldots x_{8}\}$ and $x_{3}\mathop{\perp\!\!\!\!\perp}x_{9}$ hold in the model defined by the DAG in Fig. 1, we can initialize with $X_{1}=\{x_{3}\}$, $X_{2}=\{x_{9}\}$ and $C=\emptyset$. We then use the CI tests to determine one by one to which of $X_{1}$, $X_{2}$, and $C$ the remaining variables belong. This procedure returns $X_{1}=\{x_{1},x_{2},x_{3},x_{5},x_{6}\}$, $X_{2}=\{x_{7},x_{8},x_{9}\}$, $C=\{x_{4}\}$ and hence $V_{1}=\{x_{1},x_{2},x_{3},x_{4},x_{5},x_{6}\}$ and $V_{2}=\{x_{4},x_{7},x_{8},x_{9}\}$. In this case, we note that the marginal models for $V_{1}$ and $V_{2}$ are defined by the sub-DAGs induced by $V_{1}$ and $V_{2}$, respectively. If the initial values of $X_{1}$, $X_{2}$, and $C$ or the scanning order of $w\in X\setminus(X_{1}\cup X_{2}\cup C)$ changes, $V_{1}$ and $V_{2}$ may also change.

Suppose that $X$ is grouped into $V_{1},\ldots,V_{K}$ by the SADA’s variable grouping procedure. SADA applies a causal structure learning algorithm such as DirectLiNGAM to each of the $K$ groups. Let $\hat{G}_{1}=(V_{1},E_{1}),\ldots,\hat{G}_{K}=(V_{K},E_{K})$ be the estimated $K$ DAGs for each group, where $E_{k}\subset V_{k}\times V_{k}$, $k=1,\ldots,K$ are the set of directed edges in $\hat{G}_{k}$. Then, the entire causal DAG is estimated by $\hat{G}=(X,E_{1}\cup\cdots\cup E_{K})$.

However, SADA has several practical problems. In SADA, the marginal model $\hat{G}_{k}$, $k=1,\ldots,K$ is not necessarily the causal model induced by $V_{k}$ even if we knew the correct CI relationships between variables. Assume that the true causal DAG is Fig 2 and that $\theta\leq 3$. Since $x_{1}\mathop{\perp\!\!\!\!\perp}x_{4}\mid(x_{2},x_{3})$, we can initialize with $X_{1}=\{x_{1}\}$, $X_{2}=\{x_{4}\}$, and $C=\{x_{2},x_{3}\}$, which does not need further variable checking. As a result, we group the original variable set into $V_{1}=\{x_{1},x_{2},x_{3}\}$ and $V_{2}=\{x_{2},x_{3},x_{4}\}$. For $V_{1}$, the marginal model is defined by the sub-DAG induced by $x_{1}$, $x_{2}$ and $x_{3}$. However, the marginal model with respect to $V_{2}$ is not the model defined by the sub-DAG induced by $(x_{2},x_{3},x_{4})$ because $x_{2}\mathop{\not\perp\!\!\!\!\perp}x_{3}$.

SADA does not guarantee that it always returns a DAG. Cycles may appear when DAG estimates for each variable group are merged. Suppose $V_{1}$ and $V_{2}$ share $x_{1}$ and $x_{2}$. If $x_{1}\to x_{2}\in E_{1}$, $x_{2}\to x_{1}\in E_{2}$, the cycle $x_{1}\to x_{2}\to x_{1}$ appears after merging $\hat{G}_{1}$ and $\hat{G}_{2}$. Therefore, SADA needs to remove redundant edges in cycles.

In addition, SADA requires too many CI tests for grouping variables. Consider the simple model where $x_{1},\ldots,x_{p}$ are mutually independent and no edge exists in the true DAG. In this case, any $x_{i}$ and $x_{j}$ satisfy $x_{i}\mathop{\perp\!\!\!\!\perp}x_{j}\mid(X\setminus\{x_{i},x_{j}\})$. In the worst case, $2^{p-2}$ CI tests are required to ensure $\hat{C}=\emptyset$. Thus, SADA is not feasible for high-dimensional data.

Zhang et al. [2] proposed a refined grouping scheme called the causality partitioning (CAPA). Define the order of the CI test by the number of variables in the conditioning set. CAPA reduces time complexity by setting an upper limit on the order of CI tests and using only CI tests of orders below that limit.

When the maximum order of the CI test is $\sigma_{\max}$, the time complexity of CAPA is shown to be $O(p^{\sigma_{\max}+2})$. In implementation, $\sigma_{\max}$ is set to small numbers, like two or three.

The procedure for grouping $X$ into two subsets by CAPA is as follows. Let $\sigma$ be the current order of the conditioned set. Let $M^{\sigma}=\{M^{\sigma}_{ij}\}$ be a $p\times p$ matrix such that

Let $G^{\sigma}=(X,E^{\sigma})$ be an undirected graph with $M^{\sigma}$ as its adjacency matrix. Let $d_{i}=\sum_{r=1,r\neq i}^{p}M^{\sigma}_{ir}$, $i=1,\ldots,p$ be the degrees of $x_{i}$ in $G^{\sigma}$.

1. 1.

   Initialize with $\sigma=0$

2. 2.

   Initialize with $A=B=C=D=\emptyset$

3. 3.

   For all $x_{i}\in X$ in ascending order by $d_{i}$, $i=1,\ldots,p$ do at most one of the followings so that $A\neq\emptyset$, $B\neq\emptyset$ at output

   • (a)

     $A\leftarrow A\cup\{x_{i}\}$ if $\forall x_{j}\in B$, $(x_{i},x_{j})\notin E_{\sigma}$

   • (b)

     $B\leftarrow B\cup\{x_{i}\}$ if $\forall x_{j}\in A$, $(x_{i},x_{j})\notin E_{\sigma}$

4. 4.

   $C\leftarrow X\setminus(A\cup B)$

5. 5.

   For all $x_{i}\in A\cup B$
   $D\leftarrow D\cup\{x_{i}\}$ if $\exists x_{j}\in C$, $(x_{i},x_{j})\in E_{\sigma}$

6. 6.

   $X_{1}=A\cup C\cup D$, $X_{2}=B\cup C\cup D$

7. 7.

   Return $\{X_{1},X_{2}\}$ if $\max(|X_{1}|,|X_{2}|)<|X|$
   Else if $\sigma=\sigma_{\max}$, exit
   Else $\sigma\leftarrow\sigma+1$ and go to 2

Let $S=X_{1}\cap X_{2}$. When we know the true CI relationship between $X$, $S$ always d-separates $X_{1}\setminus S$ and $X_{2}\setminus S$ in the true causal DAG. Therefore, the marginal model for $X_{1}$ and $X_{2}$ is guaranteed to be defined by the sub-DAG induced by $X_{1}$ and $X_{2}$, respectively. In general, the ascending order of $d_{i}$, $i=1,\ldots,n$ is not unique, and the output $\{X_{1},X_{2}\}$ may change depending on the choice of the ascending order. In Step 3, both 3(a) and 3(b) may be possible depending on $x_{i}$, in which case one of them is randomly selected, which may also change the output.

We illustrate the CAPA procedure with $\sigma=1$ to create two subsets where the true causal DAG is the one in Fig 1. When $\sigma=1$, $G_{\sigma}$ is an undirected graph in Fig. 3. Table II shows an example of the results of Step 3 and Step 4 when $x_{3},x_{5},x_{6},x_{9},x_{2},x_{8},x_{1},x_{7},x_{4}$ is used as the order of $x_{i}$ by $d_{i}$. Since $C=\{x_{4}\}$, we have $D=\{x_{1},x_{5},x_{6},x_{7}\}$. Therefore $X_{1}$ and $X_{2}$ are

respectively.

In this example, $|X_{1}|=7$, $|X_{2}|=7$, so the size of each group is not so small. CAPA can also be repeated recursively to split variables into smaller subsets. CAPA tends to fail to group finely in causal DAGs that contain vertices with high outdegree and low indegree.

RCD [19] is a novel causal structure learning algorithm that can be applied even when the model contains latent confounders. RCD also assumes LiNGAM, i.e., the causal relationships are linear, and the exogenous variables are independently distributed as continuous non-Gaussian distributions. The processes in RCD are divided into the following three parts: ancestral relationship finding, parental relationship finding, and confounder determining.

The first step in the RCD is to determine the ancestral relationship between variables by repeatedly conducting simple linear regressions and independence tests on each variable pair, in a similar way as when determining causal order in DirectLiNGAM, and to create a list of ancestor sets for each variable. Let $\widehat{Anc}_{i}$ be the estimated ancestor set of $x_{i}$.

Once ancestral relationships are estimated, RCD extracts parent-child relationships between each pair of variables using CI tests. Assume $x_{j}\in\widehat{Anc}_{i}$. For $x_{i}$ and $x_{j}$, if $x_{i}\mathop{\not\perp\!\!\!\!\perp}x_{j}\mid\widehat{Anc}_{i}\setminus\{x_{j}\}$, then $x_{j}$ can be determined as a parent of $x_{i}$. If a variable pair in the RCD’s output for which the direction of causality cannot be identified exists, we conclude that there are latent confounders between them.

Our study assumes that the causal model does not contain latent confounders. RCD can also be applied to the model without latent confounders. RCD is interpreted as another divide-and-conquer approach to improve scalability because it estimates the causal DAG by estimating the parent-child relationship between variables within each variable’s ancestor set. Our proposed method also uses the algorithm to find ancestral relationships in RCD to group variables.

RCD may return a graph that contains cycles even when the sample size is large. This is because errors in estimating ancestral relationships lead to an incorrect estimation of the parent-child relationships.

The proposed method also assumes LiNGAM (1). Like SADA and CAPA, the proposed method begins by grouping variables into multiple subsets based on the ancestral relationships between variables. The ancestral relationships are estimated in the same way as in RCD. When the true causal DAG is connected, the proposed method defines the variable grouping by the maximal elements in the family of the union of each variable and its ancestors. Thus, if the ancestral relationships are correctly estimated, variables are partitioned into as many groups as the number of sink nodes in the true causal DAG. The proposed method groups finely for DAGs with high outdegree and low indegree. The proposed method’s time complexity for variable grouping is $O(p^{3})$, which is the same order as CAPA with $\sigma_{\max}=1$.

When the ancestral relationships among the variables are correctly estimated, the marginal model for variables in each group obtained by the proposed method is the LiNGAM defined by the sub-DAG induced by the variables in each group.

In RCD, the estimated ancestral relationships are used to estimate parent-child relationships among variables. In contrast, the estimated ancestral relationships are only used to group variables in the proposed method. For each group, DirectLiNGAM is applied to estimate sub-DAGs, which are then merged to estimate the entire causal DAG.

Section III provides a more comprehensive description of our proposed method. Section IV compares the performance of the proposed method, the original DirectLiNGAM, CAPA, and RCD in the absence of latent confounders by computer experiments.

This subsection summarizes the algorithm for finding the ancestral relationships between variables in $X$ in RCD [19]. To determine the ancestral relationship between $x_{i}$ and $x_{j}$, $i\neq j$, Maeda and Shimizu [19] considered the simple regression models

where $u_{i}$ and $u_{j}$ are error terms. Maeda and Shimizu [19] focused on the independence relationship between the independent variable and the error term in each model. The following proposition holds for the ancestral relationship between $x_{i}$ and $x_{j}$.

This algorithm first checks which conditions 1 through 4 in Proposition 1 holds for each variable pair $(x_{i},x_{j})$. For implementation, the error terms $u_{i}$ and $u_{j}$ are replaced with the OLS residuals. If $(x_{i},x_{j})$ satisfies condition 4, the determination of the ancestral relationship between $(x_{i},x_{j})$ is withheld. Let $R$ be the set of variable pairs $(x_{i},x_{j})$ that satisfies condition 4. Let $CA^{*}_{ij}$ be the set of common ancestors of $(x_{i},x_{j})$ found during checking Proposition 1 for all pairs $(x_{i},x_{j})$.

Assume that $(x_{i},x_{j})\in R$. Then, $CA_{ij}^{*}$ forms part of the backdoor path for $x_{i}$ and $x_{j}$. To remove the influence of $CA^{*}_{ij}$ from $x_{i}$ and $x_{j}$, consider the regression models

where we identify $CA^{*}_{ij}$ with the vector of common ancestors of $(x_{i},x_{j})$. Furthermore, consider the following regression models for the error terms $v_{i}$ and $v_{j}$,

Then Proposition 1 is generalized as follows.

Proposition 1 is the case where $CA_{ij}=\emptyset$. For implementation, the error terms $v$ and $u$ are replaced with OLS residuals. If $(x_{i},x_{j})$ satisfies condition 4, the determination of the ancestral relationship between $x_{i}$ and $x_{j}$ is withheld. After checking Proposition 2 for all $(x_{i},x_{j})\in R$, update $CA^{*}_{ij}$ and $R$. If $R\neq\emptyset$, recheck Proposition 2. Theoretically, if the model contains no latent confounders, the procedure in Proposition 2 can be repeated until $R=\emptyset$ to completely determine the ancestral relationship of all $(x_{i},x_{j})$.

The pseudo-code of this algorithm is described in Algorithm V-A in Appendix A. Since the linearity is assumed between variables, Pearson’s correlation test [20] can be used to test the marginally independence, $x_{i}\mathop{\perp\!\!\!\!\perp}x_{j}$ and $v_{i}\mathop{\perp\!\!\!\!\perp}v_{j}$. For testing the independence of independent variables and residuals, Shimizu and Maeda [19, 21] used the Hilbert-Schmidt independence criterion (HSIC) [22]. This paper uses a Kernel-based conditional independence (KCI) test [23] because the accuracy of the test remains the same, and the computation time is faster.

Simply checking the conditions of Proposition 1 or Proposition 2 may yield ancestor relationships such that the estimated graph contains directed cycles due to errors in testing. The proposed method adds a heuristic procedure to avoid directed cycles. See Algorithm V-A in the Appendix for details.

This subsection introduces a procedure for grouping $X$ into several subsets according to the estimated ancestral relationships among variables and merging the estimated causal DAGs for each group.

Let $\widehat{Anc}_{i}$ be the estimated ancestor set of $x_{i}$. Define ${\cal V}$ by the output of the following procedures.

1. 1.

   Let ${\cal V}_{0}$ be the family of maximal sets in

   $$\left\{\{x_{i}\}\cup\widehat{Anc}_{i},\;i=1,\ldots,p\right\}$$

2. 2a)

   If $|\bm{v}|>1$ for all $\bm{v}\in{\cal V}_{0}$, then ${\cal V}\leftarrow{\cal V}_{0}$

3. 2b)

   Otherwise, ${\cal V}\leftarrow{\cal V}_{0}$
   For all $\bm{v}\in{\cal V}_{0}$ with $|\bm{v}|=1$

   $${\cal V}\leftarrow({\cal V}\setminus\bm{v})\cup\left\{\bigcup_{\bm{v}^{\prime}\in{\cal V}_{0},\bm{v}^{\prime}\neq\bm{v}}\{\bm{v}\cup\bm{v}^{\prime}\}\right\}$$

4. 3)

   Return ${\cal V}$

We call the series of procedures for obtaining ${\cal V}$ the causal ancestral-relationship-based grouping (CAG). If the ancestral relationships among $X$ are correctly estimated and ${\cal V}={\cal V}_{0}$, each element of $\mathcal{V}$ is the union of a sink node in $G$ and its ancestors. We also note that $\bigcup_{\bm{v}\in{\cal V}}\bm{v}=X$. Now, we have the following theorem.

Theorem 1 claims that we can consistently estimate the sub-DAGs $G_{\bm{v}}$ for $\bm{v}\in{\cal V}$ by applying a causal structure learning algorithm such as DirectLiNGAM to $\bm{v}$. Once we obtain the estimates $\hat{G}_{\bm{v}}=(\bm{v},E_{\bm{v}})$, where $E_{\bm{v}}=(\bm{v}\times\bm{v})\cap E$, the entire causal graph $G$ can be estimated by $\hat{G}=(X,\bigcup_{\bm{v}\in{\cal V}}E_{\bm{v}})$.

When the sample size is small, many ancestral relationships cannot be detected due to the type II error of the CI test in the CAG, so there may be many variables for which the ancestor set is empty, resulting in many groups consisting of a single variable in ${\cal V}_{0}$. Using such grouping would make the causal DAG estimates too sparse. So here, the group consisting of a single variable in ${\cal V}_{0}$ is merged with the other groups according to the above procedure 2b. Note that Theorem 1 holds even if ${\cal V}_{0}$ contains a group consisting of one variable.

Table III presents the correct $Anc_{i}$ and $\bm{v}_{i}$ for $i=1,\ldots,9$ of the DAG in Figure 1. In this case, we can easily see that ${\cal V}_{0}={\cal V}=\{\bm{v}_{3},\bm{v}_{5},\bm{v}_{6},\bm{v}_{9}\}$. By applying the causal structure learning algorithm to each element of ${\cal V}$, we can consistently estimate the sub-DAGs as shown in Fig. 4. By merging the sub-DAGs as $(X,E_{3}\cup E_{5}\cup E_{6}\cup E_{9})$, we can obtain the true causal DAG in Figure 1.

In the same way as SADA, the estimated entire DAG may contain cycles due to errors in learning DAGs for each group, even if the correct ancestral relationships are estimated. The proposed method implements topological sorting on the merged graph to check if a directed cycle exists in the merged graph. In SADA, the significance of each coefficient $b_{ij}$ for all edges in a cycle is assessed by the Wald test, and the least significant edge in the cycle is removed from the cycle. In the proposed method, since the ancestral relationships between variables are estimated, directed cycles can also be eliminated by removing edges in the cycle that contradict the estimated ancestral relationships. Another possible approach would be removing the edge in the cycle with the smallest absolute value of coefficient $|b_{ij}|$. The computer experiment in Section IV compares the performance of these methods for removing cycles.

It suffices to focus on the number of CI tests to evaluate the time complexity of variable grouping by CAG. The case that requires the most CI tests is when the true causal DAG is fully connected, i.e., $Anc_{i}=Pa_{i}$ for all $x_{i}\in X$. In that case, since the proposed algorithm performs CI tests while removing a source node one by one from a DAG, the number of required CI tests is

As mentioned in Section II-A, the time complexity of DirectLiNGAM is $O(np^{3}M^{2}+p^{4}M^{3})$. Even allowing for the time complexity of grouping variables, if $n$ is small relative to $p$, CAG is expected to reduce overall computation time compared to the original DirectLiNGAM without variable grouping.

Let $n_{\max}$ be the largest cardinality of $\bm{v}\in{\cal V}$. Then, the required sample size is reduced to $n_{\max}+1$.

The number of variables $p$ and the number of edges $|E|$ were set to $(p,|E|)=(10,5)$, $(10,9)$, $(20,19)$, $(30,29)$ and $(40,39)$. The sample sizes $n$ were set to $11,25,50,100$ and $500$.

The error terms $e_{i}$ were generated from the uniform distribution $U(-1,1)$. The true causal DAGs were randomly generated with $(p,|E|)$ fixed at each iteration. In these experiments, the true causal DAGs were assumed to be sparse, and the indegrees of all vertices were set to at most one. The nonzero elements of the coefficient matrix $B$ were randomly generated from the uniform distribution $U([-1,-0.5]\cup[0.5,1])$ at each iteration. The number of iterations was set to 100.

Let $|E_{c}|$ be the number of correctly detected edges in the 100 estimated DAGs. Denote the number of redundant and missing edges in the 100 estimated DAGs by $|E_{r}|$ and $|E_{m}|$, respectively. Note that if an edge in an estimated DAG is in the wrong direction, both $|E_{r}|$ and $|E_{m}|$ will add one. We evaluated the performance of CAG-LiNGAM and the other methods using the following indices (e.g., Zhu et al. [24]).

1. 1.

   The precision

   $$Pre:=\frac{|E_{c}|}{|E_{c}|+|E_{r}|}$$

   is the ratio of correctly estimated edges among the edges in the estimated DAG;

2. 2.

   The recall

   $$Rec:=\frac{|E_{c}|}{|E_{c}|+|E_{m}|}=\frac{|E_{c}|}{|E|}$$

   is the ratio of correctly estimated edges among the edges in the true DAG;

3. 3.

   $F$-measure

   $$F:=\frac{2\times Pre\times Rec}{Pre+Rec}$$

   is the harmonic mean of the precision and the recall;

4. 4.

   $Time$ is the running time in seconds for estimating 100 DAGs.

Let $\alpha_{P}$ be the significance level of Pearson’s correlation test in CAG and RCD. Let $\alpha_{CI}$ be the significance level of KCI in CAG, RCD, and CAPA. $\alpha_{P}$ and $\alpha_{CI}$ are set to

respectively.

Experiments for $p=40$ were conducted on a machine with a 3.0 GHz Intel Core i9 processor and 256 GB memory, and experiments for $p=30$ were performed on a machine with a 3.0 GHz Intel Core i9 processor and 128 GB memory. The other experiments were conducted on a machine with a 2.1 GHz Dual Xeon processor and 96 GB memory.