DAGs with NO TEARS: Continuous Optimization for Structure Learning

Introduction

Learning directed acyclic graphs (DAGs) from data is an NP-hard problem (Chickering, 1996; Chickering et al., 2004), owing mainly to the combinatorial acyclicity constraint that is difficult to enforce efficiently. At the same time, DAGs are popular models in practice, with applications in biology (Sachs et al., 2005), genetics (Zhang et al., 2013), machine learning (Koller and Friedman, 2009), and causal inference (Spirtes et al., 2000). For this reason, the development of new methods for learning DAGs remains a central challenge in machine learning and statistics.

In this paper, we propose a new approach for score-based learning of DAGs by converting the traditional combinatorial optimization problem (left) into a continuous program (right):

The main thrust of this work is to re-formulate score-based learning of DAGs so that standard smooth optimization schemes such as L-BFGS (Nocedal and Wright, 2006) can be leveraged. To accomplish this, we make the following specific contributions:

We develop an equality-constrained program for simultaneously estimating the structure and parameters of a sparse DAG from possibly high-dimensional data, and show how standard numerical solvers can be used to find stationary points.

We demonstrate the effectiveness of the resulting method in empirical evaluations against existing state-of-the-arts. See Figure 1 for a quick illustration and Section 5 for details.

We compare our ouput to the exact global minimizer (Cussens, 2012), and show that our method attains scores that are comparable to the globally optimal score in practice, although our methods are only guaranteed to find stationary points.

Most interestingly, our approach is very simple and can be implemented in about 50 lines of Python code. As a result of its simplicity and effortlessness in its implementation, we call the resulting method NOTEARS: Non-combinatorial Optimization via Trace Exponential and Augmented lagRangian for Structure learning. The implementation is publicly available at https://github.com/xunzheng/notears.

Background

2 Previous work

Popular score functions include BDe(u) (Heckerman et al., 1995), BGe (Kuipers et al., 2014), BIC (Chickering and Heckerman, 1997), and MDL (Bouckaert, 1993). Unfortunately, (4) is NP-hard to solve Chickering (1996); Chickering et al. (2004) owing mainly to the nonconvex, combinatorial nature of the optimization problem. This is the main drawback of existing approaches for solving (4): The acyclicity constraint is a combinatorial constraint with the number of acyclic structures increasing superexponentially in $d$ (Robinson, 1977). Notwithstanding, there are algorithms for solving (4) to global optimality for small problems (Ott and Miyano, 2003; Singh and Moore, 2005; Silander and Myllymaki, 2006; Xiang and Kim, 2013; Cussens, 2012; Cussens et al., 2017). There is also a wide literature on approximate algorithms based on order search (Teyssier and Koller, 2005; Schmidt et al., 2007; Scanagatta et al., 2015, 2016), greedy search (Heckerman et al., 1995; Chickering, 2003; Ramsey et al., 2016), and coordinate descent (Fu and Zhou, 2013; Aragam and Zhou, 2015; Gu et al., 2018). By searching over the space of topological orderings, the former order-based methods trade-off the difficult problem of enforcing acyclicity with a search over $d!$ orderings, whereas the latter methods enforce acyclicity one edge at a time, explicitly checking for acyclicity violations each time an edge is added. Other approaches that avoid optimizing (4) directly include constraint-based methods (Spirtes and Glymour, 1991; Spirtes et al., 2000), hybrid methods (Tsamardinos et al., 2006; Gámez et al., 2011), and Bayesian methods (Ellis and Wong, 2008; Zhou, 2011; Niinimäki et al., 2016).

The intractable form of the program (4) has led to a host of heuristic methods, often borrowing tools from the optimization literature, but always resorting to clever heuristics to accelerate algorithms. Here we briefly discuss some of the pros and cons of existing methods. While not all methods suffer from all of the problems highlighted below, we are not aware of any methods that simultaneously avoid all of them.

Broadly speaking, there are two camps: Approximate algorithms and exact algorithms, the latter of which are guaranteed to return a globally optimal solution. Exact algorithms form an intriguing class of methods, but as they are based around an NP-hard combinatorial optimization problem, these methods remain computationally intractable in general. For example, recent state-of-the-art work (Cussens et al., 2017; Chen et al., 2016) only scale to problems with a few dozen nodes (Van Beek and Hoffmann, 2015).Cussens (2012) reports experiments with $d>60$ under a constraint on the maximum parent size. Older methods based on dynamic programming methods (Ott and Miyano, 2003; Singh and Moore, 2005; Silander and Myllymaki, 2006; Xiang and Kim, 2013; Loh and Bühlmann, 2014) also scale to roughly a few dozen nodes. By contrast, state-of-the-art approximate methods can scale to thousands of nodes (Ramsey et al., 2016; Aragam and Zhou, 2015; Scanagatta et al., 2015, 2016).

Local vs. global search.

Arguably the most popular approaches to optimizing (4) involve local search, wherein edges and parent sets are added sequentially, one node at a time. This is efficient as long as each node has only a few parents, but as the number of possible parents grows, local search rapidly becomes intractable. Furthermore, such strategies typically rely on severe structural assumptions such as bounded in-degree, bounded treewidth, or edge constraints. Since real-world networks often exhibit scale-free and small-world topologies (Watts and Strogatz, 1998; Barabási and Albert, 1999) with highly connected hub nodes, these kinds of structural assumptions are not only difficult to satisfy, but impossible to check. We note here promising work towards relaxing this assumption for discrete data (Scanagatta et al., 2015). By contrast, our method uses global search wherein the entire matrix $W$ is updated in each step.

Model assumptions.

The literature on DAG learning tends to be split between methods that operate on discrete data vs. methods that operate on continuous data. When viewed from the lens of (3), the reasons for this are not clear since both discrete and continuous data can be considered as special cases of the general score-based learning framework. Nonetheless, many (but not all) of the methods cited already only work under very specific assumptions on the data, the most common of which are categorical (discrete) and Gaussian (continuous). Since (3) is agnostic to the form of the data and loss function, there is significant interest in finding general methods that are not tied to specific model assumptions.

Conceptual clarity.

Finally, on a higher level, a significant drawback of existing methods is their conceptual complexity: They are not straightforward to implement, require deep knowledge of concepts from the graphical modeling literature, and accelerating them involves many clever tricks. By contrast, the method we propose in this paper is conceptually very simple, requires no background on graphical models, and can be implemented in just a few lines of code using existing black-box solvers.

3 Comparison

It is instructive to compare existing methods for learning DAGs against other methods in the machine learning literature. We focus here on two popular models: Undirected graphical models and deep neural networks. Undirected graphical models, also known as Markov networks, is recognized as a convex problem (Yuan and Lin, 2007; Banerjee et al., 2008) nowadays, and hence can be solved using black-box convex optimizers such as CVX (Grant and Boyd, 2014). However, one should not forget score-based methods based on discrete scores similar to (4) proliferated in the early days for learning undirected graphs (e.g. Koller and Friedman, 2009, §20.7). More recently, extremely efficient algorithms have been developed for this problem using coordinate descent (Friedman et al., 2008) and Newton methods (Hsieh et al., 2014; Schmidt et al., 2009). As another example, deep neural networks are often learned using various descendants of stochastic gradient descent (SGD) (Bousquet and Bottou, 2008; Kingma and Ba, 2014; Bottou et al., 2016), although recent work has proposed other techniques such as ADMM (Taylor et al., 2016) and Gauss-Newton (Botev et al., 2017). One of the keys to the success of both of these models—and many other models in machine learning—was having a closed-form, tractable program for which existing techniques from the extensive optimization literature could be applied. In both cases, the application of principled optimization techniques led to significant breakthroughs. For undirected graphical models the major technical tool was convex optimization, and for deep networks the major technical tool was SGD.

Unfortunately, the general problem of DAG learning has not benefited in this way, and one of our main goals in the current work is to formulate score-based learning similarly as a closed-form, continuous program. Arguably, the challenges with existing approaches stem from the intractable form of the program (4). One of our main goals in the current work is to formulate score-based learning via a similar closed-form, continuous program. The key device in accomplishing this is a smooth characterization of acyclicity that will be introduced in the next section.

A new characterization of acyclicity

The values of $h$ quantify the “DAG-ness” of the graph;

$h$ and its derivatives are easy to compute.

We proceed in two steps: First, we consider the simpler case of binary adjacency matrices $B\in\{0,1\}^{d\times d}$ (Section 3.1). Note that since $\{0,1\}^{d\times d}$ is a discrete space, we cannot take gradients or do continuous optimization. For this we need the second step, in which we relax the function we originally define on binary matrices to real matrices (Section 3.2).

When does a matrix $B\in\{0,1\}^{d\times d}$ correspond to an acyclic graph? Recall the spectral radius $r(B)$ of a matrix $B$ is the largest absolute eigenvalue of $B$. One simple characterization of acyclicity is the following:

Suppose $B\in\{0,1\}^{d\times d}$ and $r(B)<1$. Then $B$ is a DAG if and only if

It essentially boils down to the fact that $\operatorname{tr}B^{k}$ counts the number of length-$k$ closed walks in a directed graph. Clearly an acyclic graph will have $\operatorname{tr}B^{k}=0$ for all $k=1,\dotsc,\infty$. In other words, $B$ has no cycles if and only if $f(B)=\sum_{k=1}^{\infty}\sum_{i=1}^{d}(B^{k})_{ii}=0$, then

Unfortunately, the condition that $r(B)<1$ is strong: although it is automatically satisfied when $B$ is a DAG, it is generally not true otherwise, and furthermore the projection is nontrivial. Alternatively, instead of the infinite series, one could consider the characterization based on finite series $\sum_{k=1}^{d}\operatorname{tr}B^{k}=0$, which does not require $r(B)<1$. However, this is impractical for numerical reasons: The entries of $B^{k}$ can easily exceed machine precision for even small values of $d$, which makes both function and gradient evaluations highly unstable. Therefore it remains to find a characterization that not only holds for all possible $B$, but also has numerical stability. Luckily, such function exists.

A binary matrix $B\in\{0,1\}^{d\times d}$ is a DAG if and only if

Similar to Proposition 1 by noting that $B$ has no cycles if and only if $(B^{k})_{ii}=0$ for all $k\geq 1$ and all $i$, which is true if and only if $\sum_{k=1}^{\infty}\sum_{i=1}^{d}(B^{k})_{ii}/k!=\operatorname{tr}e^{B}-d=0$. ∎

2 The general case: Weighted adjacency matrices

Unfortunately, the characterization (6) fails if we replace $B$ with an arbitrary weighted matrix $W$. However, we can replace $B$ with any nonnegative weighted matrix, and the same argument use to prove Proposition 2 shows that (6) will still characterize acyclicity. Thus, to extend this to matrices with both positive and negative values, we can simply use the Hadamard product $W\circ W$, which leads to our main result.

where $\circ$ is the Hadamard product and $e^{A}$ is the matrix exponential of $A$. Moreover, $h(W)$ has a simple gradient

and satisfies all of the desiderata (a)-(d).

The proof of (7) is similar to (6), and desiderata (c)-(d) follow from (8). To see why desiderata (b) holds, note that the proof of Proposition 1 shows that the power series $\operatorname{tr}(B+B^{2}+\cdots)$ simply counts the number of closed walks in $B$, and the matrix exponential simply re-weights these counts. Replacing $B$ with $W\circ W$ amounts to counting weighted closed walks, where the weight of each edge is $w_{ij}^{2}$. Thus, larger $h(W)>h(W^{\prime})$ means either (a) $W$ has more cycles than $W^{\prime}$ or (b) The cycles in $W$ are more heavily weighted than in $W^{\prime}$.

Moreover, notice that $h(W)\geq 0$ for all $W$ since each term in the series is nonnegative. This gives another interesting perspective of the space of DAGs as the set of global minima of $h(W)$. However, due to the nonconvexity, this is not equivalent to the first order stationary condition $\nabla h(W)=0$.

A key conclusion from Theorem 1 is that $h$ and its gradient only involve evaluating the matrix exponential, which is a well-studied function in numerical analysis, and whose $O(d^{3})$ algorithm (Al-Mohy and Higham, 2009) is readily available in many scientific computing libraries. Although the connection between trace of matrix power and number of cycles in the graph is well-known Harary and Manvel (1971), to the best of our knowledge, this characterization of acyclicity has not appeared in the DAG learning literature previously. We defer the discussion of other possible characterizations in the appendix. In the next section, we apply Theorem 1 to solve the program (3) to stationarity by treating it as an equality constrained program.

Optimization

Theorem 1 establishes a smooth, algebraic characterization of acyclicity that is also computable. As a consequence, the following equality-constrained program $(\mathsf{ECP})$ is equivalent to (3):

The main advantage of $(\mathsf{ECP})$ compared to both (3) and (4) is its amenability to classical techniques from the mathematical optimization literature. Nonetheless, since $\{W:h(W)=0\}$ is a nonconvex constraint, (9) is a nonconvex program, hence we still inherit the difficulties associated with nonconvex optimization. In particular, we will be content to find stationary points of (9); in Section 5.3 we compare our results to the global minimizer and show that the stationary points found by our method are close to global minima in practice.

In the follows, we outline the algorithm for solving (9). It consists of three steps: (i) converting the constrained problem into a sequence of unconstrained subproblems, (ii) optimizing the unconstrained subproblems, and (iii) thresholding. The full algorithm is outlined in Algorithm 1.

We will use the augmented Lagrangian method (e.g. Nemirovski, 1999) to solve $(\mathsf{ECP})$, which solves the original problem augmented by a quadratic penalty:

with a penalty parameter $\rho>0$. A nice property of the augmented Lagrangian method is that it approximates well the solution of a constrained problem by the solution of unconstrained problems without increasing the penalty parameter $\rho$ to infinity (Nemirovski, 1999). The algorithm is essentially a dual ascent method for (10). To begin with, the dual function with Lagrange multiplier $\alpha$ is given by

is the augmented Lagrangian. The goal is to find a local solution to the dual problem

Let $W_{\alpha}^{\star}$ be the local minimizer of the Lagrangian (11) at $\alpha$, i.e. $D(\alpha)=L^{\rho}(W_{\alpha}^{\star},\alpha)$. Since the dual objective $D(\alpha)$ is linear in $\alpha$, the derivative is simply given by $\nabla D(\alpha)=h(W_{\alpha}^{\star})$. Therefore one can perform dual gradient ascent to optimize (13):

where the choice of step size $\rho$ comes with the following convergence rate:

For $\rho$ large enough and the starting point $\alpha_{0}$ near the solution $\alpha^{\star}$, the update (14) converges to $\alpha^{\star}$ linearly.

In our experiments, typically fewer than 10 steps of the augmented Lagrangian scheme are required.

2 Solving the unconstrained subproblem

is the smooth part of the objective. Our goal is to solve the above problem to high accuracy so that $h(W)$ can be sufficiently suppressed.

In the special case of $\lambda=0$, the nonsmooth term vanishes and the problem simply becomes an unconstrained smooth minimization, for which a number of efficient numerical algorithms are available, for instance the L-BFGS (Byrd et al., 1995). To handle the nonconvexity, a slight modification (Nocedal and Wright, 2006, Procedure 18.2) needs to be applied.

When $\lambda>0$, the problem becomes composite minimization, which can also be efficiently solved by the proximal quasi-Newton (PQN) method (Zhong et al., 2014). At each step $k$, the key idea is to find the descent direction through a quadratic approximation of the smooth term:

where $\mathbf{\boldsymbol{g}}_{k}$ is the gradient of $f(\mathbf{\boldsymbol{w}})$ and $B_{k}$ is the L-BFGS approximation of the Hessian. Note that for each coordinate $j$, problem (17) has a closed form update $\mathbf{\boldsymbol{d}}\leftarrow\mathbf{\boldsymbol{d}}+z^{\star}e_{j}$ given by

Moreover, the low-rank structure of $B_{k}$ enables fast computation for coordinate update. As we describe in Appendix A, the precomputation time is only $O(m^{2}p+m^{3})$ where $m\ll p$ is the memory size of L-BFGS, and each coordinate update is $O(m)$. Furthermore, since we are using sparsity regularization, we can further speed up the algorithm by aggressively shrinking the active set of coordinates based on their subgradients (Zhong et al., 2014), and exclude the remaining dimensions from being updated. With the updates restricted to the active set $\mathcal{S}$, all dependencies of the complexity on $O(p)$ becomes $O(|\mathcal{S}|)$, which is substantially smaller. Hence the overall complexity of L-BFGS update is $O(m^{2}|\mathcal{S}|+m^{3}+m|\mathcal{S}|T)$, where $T$ is the number of inner iterations, typically $T=10$.

3 Thresholding

In regression problems, it is known that post-processing estimates of coefficients via hard thresholding provably reduces the number of false discoveries (Zhou, 2009; Wang et al., 2016). Motivated by these encouraging results, we threshold the edge weights as follows: After obtaining a stationary point $\widetilde{W}_{\mathsf{ECP}}$ of (10), given a fixed threshold $\omega>0$, set any weights smaller than $\omega$ in absolute value to zero. This strategy also has the important effect of “rounding” the numerical solution of the augmented Lagrangian (10), since due to numerical precisions the solution satisfies $h(\widetilde{W}_{\mathsf{ECP}})\leq\epsilon$ for some small tolerance $\epsilon$ near machine precision (e.g. $\epsilon=10^{-8}$), rather than $h(\widetilde{W}_{\mathsf{ECP}})=0$ strictly. However, since $h(\widetilde{W}_{\mathsf{ECP}})$ explicitly quantifies the “DAG-ness” of $\widetilde{W}_{\mathsf{ECP}}$ (see desiderata (b), Section 3), a small threshold $\omega$ suffices to rule out cycle-inducing edges.

Experiments

We compared our method against greedy equivalent search (GES) (Chickering, 2003; Ramsey et al., 2016), the PC algorithm (Spirtes et al., 2000), and LiNGAM (Shimizu et al., 2006). For GES, we used the fast greedy search (FGS) implementation from Ramsey et al. (2016). Since the accuracy of PC and LiNGAM was significantly lower than either FGS or NOTEARS, we only report the results against FGS here. This is consistent with previous work on score-based learning (Aragam and Zhou, 2015), which also indicates that FGS outperforms other techniques such as hill-climbing and MMHC (Tsamardinos et al., 2006). FGS was chosen since it is a state-of-the-art algorithm that scales to large problems.

2 Structure learning

3 Comparison to exact global minimizer

In order to assess the ability of our method to solve the original program given by (3), we used the GOBNILP program (Cussens, 2012; Cussens et al., 2017) to find the exact minimizer of (3). Since this involves enumerating all possible parent sets for each node, these experiments are limited to small DAGs. Nonetheless, these small-scale experiments yield valuable insight into how well NOTEARS performs in actually solving the original problem. In our experiments we generated random graphs with $d=10$, and then generated 10 simulated datasets containing $n=20$ samples (for high-dimensions) and $n=1000$ (for low-dimensions). We then compared the scores returned by our method to the exact global minimizer computed by GOBNILP along with the estimated parameters. The results are shown in Table 1. Surprisingly, although NOTEARS is only guaranteed to return a local minimizer, in many cases the obtained solution is very close to the global minimizer, as evidenced by deviations $\|\widehat{W}-W_{\mathsf{G}}\|$. Since the general structure learning problem is NP-hard, we suspect that although the models we have tested (i.e. ER and SF) appear amenable to fast solution, in the worst-case there are graphs which will still take exponential time to run or get stuck in a local minimum. Furthermore, the problem becomes more difficult as $d$ increases. Nonetheless, this is encouraging evidence that the nonconvexity of (9) is a minor issue in practice. We leave it to future work to investigate these problems further.

4 Real-data

We also compared FGS and NOTEARS on a real dataset provided by Sachs et al. (2005). This dataset consists of continuous measurements of expression levels of proteins and phospholipids in human immune system cells ($n=7466$ $d=11$, 20 edges). This dataset is a common benchmark in graphical models since it comes with a known consensus network, that is, a gold standard network based on experimental annotations that is widely accepted by the biological community. In our experiments, FGS estimated 17 total edges with an SHD of 22, compared to 16 for NOTEARS with an SHD of 22.

Discussion

We have proposed a new method for learning DAGs from data based on a continuous optimization program. This represents a significant departure from existing approaches that search over the discrete space of DAGs, resulting in a difficult optimization program. We also proposed two optimization schemes for solving the resulting program to stationarity, and illustrated its advantages over existing methods such as greedy equivalence search. Crucially, by performing global updates (e.g. all parameters at once) instead of local updates (e.g. one edge at a time) in each iteration, our method is able to avoid relying on assumptions about the local structure of the graph. To conclude, let us discuss some of the limitations of our method and possible directions for future work.

First, it is worth emphasizing once more that the equality constrained program (9) is a nonconvex program. Thus, although we overcome the difficulties of combinatorial optimization, our formulation still inherits the difficulties associated with nonconvex optimization. In particular, black-box solvers can at best find stationary points of (9). With the exception of exact methods, however, existing methods suffer from this drawback as well.GES (Chickering, 2003) is known to find the global minimizer in the limit $n\to\infty$ under certain assumptions, but this is not guaranteed for finite samples. The main advantage of NOTEARS then is smooth, global search, as opposed to combinatorial, local search; and furthermore the search is delegated to standard numerical solvers.

Second, the current work relies on the smoothness of the score function, in order to make use of gradient-based numerical solvers to guide the graph search. However it is also interesting to consider non-smooth, even discrete scores such as BDe (Heckerman et al., 1995). Off-the-shelf techniques such as Nesterov’s smoothing (Nesterov, 2005) could be useful, however more thorough investigation is left for future work.

Third, since the evaluation of the matrix exponential is $O(d^{3})$, the computational complexity of our method is cubic in the number of nodes, although the constant is small for sparse matrices. In fact, this is one of the key motivations for our use of second-order methods (as opposed to first-order), i.e. to reduce the number of matrix exponential computations. By using second-order methods, each iteration make significantly more progress than first-order methods. Furthermore, although in practice not many iterations ($t\sim 10$) are required, we have not established any worst-case iteration complexity results. In light of the results in Section 5.3, we expect there are exceptional cases where convergence is slow. Notwithstanding, NOTEARS already outperforms existing methods when the in-degree is large, which is known difficult spot for existing methods. We leave it to future work to study these cases in more depth.

Lastly, in our experiments, we chose a fixed, suboptimal value of $\omega>0$ for thresholding (Section 4.3). Clearly, it would be preferable to find a data-driven choice of $\omega$ that adapts to different noise-to-signal ratios and graph types. It is an intersting direction for future to study such choices.

The code is publicly available at https://github.com/xunzheng/notears.

References

Appendix A Details of Proximal Quasi-Newton

The detailed procedure of PQN is outlined in Algorithm 2.

Appendix B Sensitivity of threshold

We demonstrate the effect of threshold in Figure 4. For each setting, we computed the “ROC” curve for FDR and TPR with varying level of threshold, while ensuring the resulting graph is indeed a DAG. On the right, we also present the estimated edge weights of $\widetilde{W}_{\mathsf{ECP}}$ in decreasing order. One can first observe that in all cases most of the edge weights are equal or close to zero as expected. The remaining question is how to choose a threshold that separates out these (near zero) from signals (away from zero) so that best performance can be achieved. With enough samples, one can often notice a sudden change in the weight distribution as in Figure 4(a)(c). With insufficient samples, the breakpoint is less clear, and the optimal choice that balances between TPR and FDR is depends on the specific settings. Nonetheless, the predictive performance is less sensitive to threshold value as one can see from the slope of the decrease in the weights before getting close to zero. Indeed, in our experiments, we found a fixed threshold $\omega=0.3$ is a suboptimal yet reasonable choice across many different settings.

Appendix C Sensitivity of weight scale

We investigate the effect of weight scaling to the NOTEARS algorithm in Figure 5. In particular, we run experiments with $w_{ij}\in\alpha\cdot[0.5,2]\cup-\alpha\cdot[0.5,2]$ with $\alpha\in\{1.0,0.9,0.8,\dotsc,0.1\}$. On the left, we plot the smallest threshold $\omega$ required to obtain a DAG (see Section 4.3) for different scale $\alpha$. Overall, across different values of $\alpha$, the variation in the smallest $\omega$ required is minimal. We also hasten to point out that this also decreases the signal to noise ratio (SNR), which more directly affects the accuracy. Indeed, in the figure on the right, we can observe (as expected) some performance drop when using smaller value of $\alpha$.

Appendix D Experiments

We used simulated graphs from two well-known ensembles of random graphs:

Erdös-Rényi (ER). Random graphs whose edges are added independently with equal probability $p$. We simulated models with $d$, $2d$, and $4d$ edges (in expectation) each, denoted by ER-1, ER-2, and ER-4, respectively.

Scale-free networks (SF). Networks simulated according to the preferential attachment process described in Barabási and Albert (1999). We simulated scale-free networks with $4d$ edges and $\beta=1$, where $\beta$ is the exponent used in the preferential attachment process.

Gaussian noise ($\mathsf{Gauss}$). $z\sim\mathcal{N}(0,I_{d\times d})$.

Exponential noise ($\mathsf{Exp}$). $z_{j}\sim\operatorname{Exp}(1)$, $j=1,\ldots,d$.

Gumbel noise ($\mathsf{Gumbel}$). $z_{j}\sim\operatorname{Gumbel}(0,1)$, $j=1,\ldots,d$.

For each dataset, we ran FGS, PC, and LinGAM and NOTEARS to compare the performance in reconstructing the DAG $B$. We used the following implementations:

FGS and PC were implemented through the py-causal package, available at https://github.com/bd2kccd/py-causal. Both of these methods are written in highly optimized Java code.

LinGAM was implemented using the author’s Python code: https://sites.google.com/site/sshimizu06/lingam.

Since the accuracy of PC and LiNGAM was significantly lower than either FGS or NOTEARS, we only report the results against FGS. A few comments on FGS are in order: 1) FGS estimates a graph, so it does not output any parameter estimates; 2) Instead of returning a DAG, FGS returns a CPDAG (Chickering, 2003), which contains undirected edges; 3) FGS has a single tuning parameter that controls the strength of regularization. Thus, in our evaluations, we treated FGS favourably by treating undirected edges as true positives as long as the true graph had a directed edge in place of the undirected edge. For tuning parameters, we used the values suggested by the authors of the FGS code.

D.2 Metrics

We evaluated the learned graphs on four common graph metrics: 1) False discovery rate (FDR), 2) True positive rate (TPR), 3) False positive rate (FPR), and 4) Structural Hamming distance (SHD). Recall that SHD is the total number of edge additions, deletions, and reversals needed to convert the estimated DAG into the true DAG. Since we consider directed graphs, a distinction between True Positives (TP) and Reversed edges (R) is needed: the former is estimated with correct direction whereas the latter is not. Likewise, a False Positive (FP) is an edge that is not in the undirected skeleton of the true graph. In addition, Positive (P) is the set of estimated edges, True (T) is the set of true edges, False (F) is the set of non-edges in the ground truth graph. Finally, let (E) be the extra edges from the skeleton, (M) be the missing edges from the skeleton. The four metrics are then given by:

FDR $=(\mathit{R}+\mathit{FP})/\mathit{P}$

FPR $=(\mathit{R}+\mathit{FP})/\mathit{F}$

D.3 Further evaluations

Figure 7 and Figure 8 shows structure recovery results for $n=1000$ and $n=20$ for various random graphs and SEM noise types. Other than fixed $\omega$ as in the main paper, we also included the optimal choice of thresholding, marked as “best”. The trend is consistent with the main text: our method in general outperforms FGS, without tuning $\omega$ to the optimum for each setting.

Table 2 extends the global minimizer result for various random graph types. For each random graph and samples, we computed exact local scores as inputs to GOBNILP program, which finds the globally optimal structure for the given score. We can again observe that the difference between our estimate $\widehat{W}$ and global minimizer $W_{\mathsf{G}}$ is small across all cases.