A Core-Guided Approach to Learning Optimal Causal Graphs

1 Introduction

Discovery of causal relations from sample data is an important part of data analysis for various application fields, and thus an important area of AI research. Causal relation structures are generally represented as directed graphs, where nodes represent measurements and directed edges denote directed causal relations. It has been shown that many features of these structures can be determined even from passively observed data [Pearl, 2000; Spirtes et al., 2000]. The difficulty lies in determining presence of edges and their directionality, especially when latent confounders (unobserved common causes) and feedback (directed cycles) may be present.

For restricted settings without latent confounders and feedback cycles, various well-performing algorithms for causal discovery have been developed. When aiming for accuracy rather than very high scalability, the most accurate choices are exact score-based Bayesian network structure learning algorithms [Malone et al., 2015], which provide provably globally optimal graphs as solutions under well-defined (e.g. Bayesian marginal likelihood) scoring functions [Yuan and Malone, 2013; Bartlett and Cussens, 2017; Berg et al., 2014] by employing e.g. constraint optimization methods such as integer programming (IP) and maximum satisfiability (MaxSAT).

However, much less progress has been made for exact discovery algorithms when allowing for the presence of latent confounders or feedback. To handle these more general search spaces, the classic constraint-based causal discovery algorithms [Spirtes et al., 2000; Pearl, 2000] use

(in)dependence constraints obtained from statistical tests to narrow down the candidate graphs that may have produced the given data. The classic inexact algorithms, such as PC, CCD and FCI, scale up by pruning out independence tests based on earlier test results [Spirtes et al., 2000; Richardson, 1996]. Unfortunately, these methods are not very accurate in practice: early mistakes often jeopardize the learning results [Claassen and Heskes, 2012; Hyttinen et al., 2014].

Conflicting independence test results from sample data give rise to a combinatorial optimization problem which, when allowing for latent confounders and cycles, is over a drastically larger search space compared to more restricted settings such as Bayesian network structure learning. The first exact approach to this general setting recently developed by Hyttinen et al. [2014] is able find guaranteed-optimal causal structures in terms of weighted (in)dependence constraints, and exhibits higher accuracy [Hyttinen et al., 2014; Magliacane et al., 2016; Borboudakis and Tsamardinos, 2016]. The approach is based on applying MaxSAT optimization solvers-or, in this case equivalently, answer set programming [Hyttinen et al., 2014]. However, developing causal discovery algorithms with better running time performance without trading off accuracy for this general setting is a major challenge.

We take on this challenge by proposing several algorithmic techniques for speeding up the approach of Hyttinen et al. [2014], who apply a MaxSAT solver as a black box on their encoding of the causal discovery task. This leaves the actual search for a solution to an off-the-shelf MaxSAT solver that has no domain-specific knowledge on the problem domain in question. Instead, we show how to integrate domain-knowledge into a state-of-the-art MaxSAT solver for expediting the search for causal structures. Specifically, we (i) identify a number of generic cores describing small sets of conflicting constraints, (ii) develop an incremental core extraction technique, (iii) use dynamic partial encodings during search without influencing the optimality of found causal structures, and (iv) propose a problem-specific way based on linear programming relaxations for hardening constraints. The resulting solver, Dseptor, outperforms the approach of Hyttinen et al. [2014] in terms of running times on both simulated and real-world data. An implementation and an online appendix are available at the authors' homepages.

2 Exact Constraint-based Causal Discovery

We focus on the problem setting and optimization problem formulation of [Hyttinen et al., 2014]. We consider the class G of causal graphs G = (V,E) with set of nodes V , where the edge relation E is composed of directed causal edges and (symmetric) bi-directed edges (see Figure 1 for an example). Bi-directed edges X ↔ Y canonically represent latent confounders, e.g., structures X ← L → Y , where L /∈ V is an unmeasured common cause of two observed variables X and Y . A walk between X and Y is a sequence of consecutive edges in the graph (allowing for repeated edges and nodes). A node is a collider on a walk if both its adjacent edges on the walk point into the node. A walk in graph G is d-connecting w.r.t. a conditioning set S ⊆ V \ {X, Y } if every collider on the walk is in S and no other nodes on the walk are in S. Two nodes are d-connected given a conditioning set S if there is at least one d-connecting walk between them; otherwise they are d-separated. (This is equivalent to Pearl's standard definition [Studen´y, 1998].) In Figure 1, X and W are d-connected given Y through X ↔ Z → Y ← Z ← W. Nodes Y and Q are d-separated given W as all walks between them violate the d-connection criterion at node X.

Under the commonly used causal Markov and faithfulness assumptions [Spirtes et al., 2000], statistical dependence becomes equivalent to (a type of) reachability in the graph: two random variables are statistically dependent conditional on a set of variables S iff they are d-connected given S in the generating causal structure G. 1 That is, given enough samples from a model with structure in Figure 1, we would expect to find X statistically dependent on W given Y , and Y statistically independent of Q given W. Thus, throughout the rest of the paper we use X ⊥⊥ Y |S (X ̸⊥⊥ Y |S) to denote statistical independence (dependence) and d-separation (d-connection).

The aim of causal discovery is to recover as many properties of Gt as possible from the data, which we take to be an i.i.d. dataset sampled from a passive observational distribution that is Markov and faithful to an underlying "true" causal graph Gt. The (in)dependence constraints K are obtained by running statistical independence tests on the data. Since the tests produce both false positives and false negatives when run on finite sample data, exact constraint-based causal discovery can be viewed as the following abstract constrained optimization problem [Hyttinen et al., 2014].

INPUT: A set K of conditional (in)dependence constraints over V , and a non-negative weight (cost) w(k) for each k ∈ K.

TASK: Find a causal graph G∗ = (V,E∗ ) s.t.

In words, our goal is to find a single representative graph G∗ that minimizes the sum of the weights of the given conditional independence and dependence constraints not implied by G∗ . The constraints can be weighted according to their reliability, and conflicts among the constraints are well-resolved when the sum of the weights of the constraints not satisfied by the output graph are minimized. Given a representative, properties of its equivalence class can be examined in several ways, see e.g. [Hyttinen et al., 2013].

1 Spirtes [1995] shows this for the linear cyclic case. Pearl and Dechter [1996] give a proof for the discrete cyclic case, Neil [2000] points out the required additional assumptions. For the acyclic case, see e.g. [Spirtes et al., 2000].

Weights for (in)dependence constraints can be obtained in many ways. Claessen and Heskes [2012] calculate Bayesian probabilities by marginalizing over Bayesian network structures. Ways to turn frequentist p-values into reliability weights have recently been proposed [Magliacane et al., 2016; Triantafillou and Tsamardinos, 2015]. While the algorithms we develop are independent of the choice of weights, we obtain weights from local Bayesian model selection [Hyttinen et al., 2014; Margaritis and Bromberg, 2009] for our experiments.

Note that the optimization problem is computationally challenging. For obtaining good accuracy, a large number of (in)dependence constraints K are needed; we use all testable (in)dependence constraints ((n2 )2n-2 for n nodes). The dseparation condition for a solution satisfying a particular (in)dependence constraint is also quite intricate. On the other hand, this separation condition can be relatively naturally encoded declaratively as Boolean constraints. We give here an intuitive overview of the encoding of Hyttinen et al. [2014] in terms of MaxSAT to the extent necessary for the rest of this paper. In MaxSAT, problems are encoded using propositional clauses, using hard clauses (which have to be satisfied by all solutions) and weighted soft clauses (which, if not satisfied, incur a cost on a solution equal to their weights). In [Hyttinen et al., 2014] each (in)dependence constraint k ∈ K is encoded as a unit soft clause over a distinct Boolean variable representing k with weight w(k). Additional Boolean variables are used for representing the solutions searched over, i.e., the edge relation of causal graphs. The d-connecting walks are encoded as hard clauses, linking the edge relation with the (in)dependence constraints.

3 Implicit Hitting Set Approach

We build on the encoding of Hyttinen et al. [2014] represented as MaxSAT [Berg et al., 2015], integrating several domain-specific algoritmic improvements into the opensource MaxSAT solver LMHS [Saikko et al., 2016] to further improve the running times of the approach. LMHS implements the generic implicit hitting set approach [MorenoCenteno and Karp, 2013] for MaxSAT [Davies and Bacchus, 2013]. Analogously to MaxSAT, the causal discovery instances (D,K,w) considered here consist of a (hard) dseparability condition D, a set of (soft) (in)dependence constraints K, and a weight function w which characterizes the reliability of constraints in K.

In this context, the implicit hitting set approach of LMHS

(see Figure 2) works via an alternating sequence of calls to a Boolean satisfiability (SAT) solver and an integer programming (IP) solver. Each is used for a task it is well-suited for: the SAT solver for proving unsatisfiability, and the IP solver for efficiently handling a weighted objective function.

In detail, the SAT solver is used to accumulate a set C of (unsatisfiable) cores of the instance, i.e., subsets of constraints c ⊂ K which (together with D) are not satisfied by any causal graph. A core of an unsatisfiable instance in readily provided by modern SAT solvers without additional overhead once unsatisfiability is established. After each core extraction step, an IP solver is called on C to optimally solve the minimum-cost hitting set problem HSw(C) over the cores in C with weights given by w, i.e., to find a set H ⊆ ⋃ C minimizing ∑x∈H w(x) subject to H ∩ c ̸= ∅ for each c ∈ C, or in words, a minimum-cost set of constraints which intersects with each of the currently accumulated cores.

After a hitting set H is computed, the constraints in H are removed from the next SAT solver call. Concretely, this next call checks the satisfiability of a CNF formula consisting of the encoding of D and the encoding of the constraints of K not in H. Once the SAT solver finds a solution τ, i.e., finds a causal graph that satisfies K \ H and D, the solution given by the SAT solver is guaranteed to be an optimal solution, i.e., the optimal cost is ∑x∈H w(x). The intuition here, as shown for MaxSAT by Davies and Bacchus [2011], is that the optimal causal graphs of (D, K, w) all satisfy K \ H and D for some minimum-cost hitting set H over all unsatisfiable cores of the instance. In practice, such an H can be (and typically is) found without having to accumulate every unsatisfiable core. Further, the solution to each min-cost hitting set instance solved during search provides a lower bound on the cost of optimal causal graphs.

Example. Consider the following unit-weighted constraints:

The implicit hitting set approach starts by checking whether all of the constraints are simultaneously satisfiable, i.e., whether some causal graph satisfies K and D. No such graph exists, so the SAT solver will return an unsatisfiable core of the instance, for example, c1 = {X ̸⊥⊥ Z|W; Y ̸⊥⊥ Z|W; X ⊥⊥ Y |W; X ⊥⊥ Y |Z,W}. Thus, at least one of the four constraints has to be an error in the statistical tests (under the used assumptions). The search proceeds by solving the minimum-cost hitting set problem over this single core using an IP solver. One solution is

Table 1: Common domain-specific core types.

• (i) {X ⊥⊥ Z|S; X ⊥⊥ Y |S; X ̸⊥⊥ Z|Y, S}
• (ii) {X ̸⊥⊥ Z|S; Y ̸⊥⊥ Z|S; X ⊥⊥ Y |S; X ⊥⊥ Y |Z, S}
• (iii) {X ̸⊥⊥ Z|Y, S; Y ̸⊥⊥ Z|X, S; X ⊥⊥ Y |S; X ⊥⊥ Y |Z, S}
• (iv) {Y ̸⊥⊥ Z|S; X ̸⊥⊥ Z|S; Z ⊥⊥W|X, Y, S; X ⊥⊥ Y |Z, S; X ⊥⊥ Y |W, S}
• (v) {Y ̸⊥⊥ Z|S; X ̸⊥⊥ Z|S; Z ⊥⊥W|Y, S; X ⊥⊥ Y |Z, S; X ⊥⊥ Y |W, S}
• (vi) {X ̸⊥⊥ Y |Z, S; Y ̸⊥⊥ Z|X,W, S; W ̸⊥⊥ Y |Z, S; W ⊥⊥ X|Y, Z, S; X ⊥⊥ Z|W, S}
• (vii) {X ̸⊥⊥ Y |Z, S; Y ̸⊥⊥ Z|X,W, S; W ̸⊥⊥ Y |S; W ⊥⊥ X|Y, S; X ⊥⊥ Z|W, S}

X ⊥⊥ Y |Z,W, which gives a lower bound of 1 for the cost of optimal solutions. The search then iterates by checking whether a causal graph satisfies constraints other than X ⊥⊥ Y |Z,W. Again no such graph exists, and we obtain a core, say c2 = {X ⊥⊥ Y |W; X ̸⊥⊥ Y |W,Q; Y ⊥⊥ Q|W}. The IP solver can then provide the min-cost hitting set {X ⊥⊥ Y |W} of {c1, c2}. We then check for a graph in which all constraints other than X ⊥⊥ Y |W are satisfied. Now the SAT solver returns a solution, which can be interpreted as the graph in Figure 1. The cost of this graph is 1, since only one constraint, X ⊥⊥ Y |W, with unit weight is violated in the graph, and the search terminates. Optimality is proven without considering the core c3 = {X ̸⊥⊥ Z|Y,W; Y ̸⊥⊥ Z|X,W; X ⊥⊥ Y |W; X ⊥⊥ Y |Z,W}.

4 Domain-Specific Cores

As the first main contribution of this work, we identify several low-cardinality cores commonly encountered in causal discovery instances. Identifying general forms (or patterns) of such cores can be very beneficial for speeding up LMHS: we can find a number of cores by a simple pattern search through the constraints K (in negligible time), instead of making potentially very time-consuming SAT solver calls for each core. We can thereby use critical domain-specific knowledge in the LMHS search simply by adding the found domain specific cores directly to the core set maintained by LMHS.

During search, low-cardinality cores especially constrain the hitting set problems (and considerably improve the lower bound); however, their extraction can be time-consuming with a SAT solver. Thus, we ran plain LMHS on the encoding of Hyttinen et al. [2014] to search for examples of small cores. We then manually generalized them and proved their validity. Table 1 lists core types over at most four nodes, these core types are valid for any set S disjoint from the nodes mentioned. The induced cores were enough for solving all instances up to four nodes in our simulations, i.e., no further cores needed to be extracted using a SAT solver. The following theorem establishes that the core types (i)-(vii) in Table 1 are subset-minimal cores for causal discovery instances.

Theorem 1 Each of the sets of (in)dependence constraints

(i)-(vii) in Table 1 is a subset-minimal core: for each of the cores (i)-(vii), it holds that 1) there is no causal graph in G that satisfies all of the (in)dependence constraints in the core (assuming faithfulness), and 2) the constraints in every proper subset of the core are satisfied by a causal graph.

For example, the intuition behind core (i) is that if conditioning on node Y induces a dependence (d-connection) between X and Z (1st vs 3rd constraint), Y would have to be dependent on X (violating the 2nd constraint). Minimality of core (i) is verified by noting that any two constraints can be trivially simultaneously satisfied.

5 Incremental Core Extraction through Dynamic Partial Encoding

If we adhere exactly to the core loop shown in Figure 2, all constraints in K \ H are input at once to the SAT solver in a single SAT call for extracting a core. Here we propose instead finding a core by incrementally adding constraints from K \H one by one until the SAT solver reports unsatisfiability and yields a core; Alg. 1 outlines this approach as FINDCOREINCREMENTAL. This technique has several potential benefits (as also previously pointed out in [Ans´otegui et al., 2016] for the related so-called stratified approach in the context of pure SAT-based MaxSAT solvers). (i) We obtain a solution τ from each SAT solver call until the call on which the solver reports unsatisfiability. Each of the satisfiable calls thus gives an upper bound which can improve the currently known best upper bound U; improved upper bounds are important for the bounds-based constraint hardening described in Sect. 6 to have an impact. (ii) The core c obtained is minimal w.r.t. the order in which the soft constraints are added. This can result in smaller cores, which makes the core minimization procedures implemented within LMHS run faster. (iii) Many of the-compared to LMHS, additional-SAT solver calls on D,K′ are satisfiable and made on a relatively small set of constraints (K′ ); hence many of these calls can be expected to be relatively fast. Also, a number of SAT solver calls can be skipped, since the solution obtained during the preceding solver call often already satisfies the next constraint to be added.

Incremental core extraction allows for using domainspecific heuristics for the order in which constraints are added. For example, one can order by decreasing weights, following heuristics from inexact causal discovery algorithms [Triantafillou and Tsamardinos, 2015; Borboudakis and Tsamardinos, 2016] and MaxSAT [Ans´otegui et al., 2016], or by increasing conditioning set size, inspired by the FCI and PC algorithms [Spirtes et al., 2000]. However, these heuristics tend to consider a relatively high number of constraints before finding a core. We found that adding constraints one by one in a random order works well here. Intuitively, a random order diversifies the search for cores to different subsets of (in)dependence constraints, which can be very beneficial for proving improved lower bounds.

Dynamic Partial Encoding. We also employ an encoding technique which allows for switching off parts of the (in total large) d-separation encoding that are redundant w.r.t. the current set K′ of (in)dependence constraints considered in single

1: procedure FINDCOREINCREMENTAL(H,U)

• 2: K′ = ∅
• 3: while K′ ̸=K \H do
• 4: Add a new element k from K \ H to K′ .
• 5: If k is satisfied by previous solution τ : continue.
• 6:

(R, τ, c) = SATSOLVER(D, K′ )

• 7: If R = SAT: update upper bound U, continue.
• 8: If R = UNSAT: return core c.
• 9: Return SAT.

Algorithm 1: Incremental core extraction in Dseptor.

SAT solver call. While the full MaxSAT encoding of a causal discovery instance is input initially to LMHS, by definition each of the SAT solver calls take into account only a subset of the (in)dependence constraints: the constraints in the latest hitting set H are disregarded in the next SAT solver call. The incremental core extraction technique further decreases the number of enforced (in)dependence constraints per SAT solver call. However, the hard clauses imposing d-separation criterion related to a specific (in)dependence constraint remain active in the SAT solver-although the actual independence is not used or considered at all. The SAT solver may not always realize that such hard clauses are redundant. In hope of making the job of the SAT solver easier, we propose a way of dynamically switching off hard clauses unrelated to a given set of (in)dependence constraints.

Consider Figure 3, where the two triangles rooted at the (in)dependence constraints ki and kj, respectively, represent the set of hard clauses that are part of the d-separation encoding for ki and kj, i.e., the clauses in the cone of influence of the respective constraints. In our case, this structure is directly given by the encoding proposed in [Hyttinen et al., 2014]. Firstly, we add the logical implication ¬bi → ki per each (in)dependence constraint ki introducing a fresh blocking variable bi. Secondly, for clauses in the non-shared part of the cone marked with bi (similarly for bj), we augment each clause with the blocking variable bi: we replace each clause cl by logical implication ¬bi → cl. This means that bi = 0 implies that the clauses have to be satisfied, and bi = 1 relaxes (blocks) the clauses. For the shared part of the cones marked with bij, we again augment the clauses in that part with the blocking variable bij. We also add the logical implication (¬bi ∨ ¬bj) → ¬bij as a hard constraint. Thus bi = 0 (or bj = 0) enforces the clauses in the shared part as well. The idea here resembles the approach proposed for dealing with dependencies between blocking variables in other contexts [Lagniez and Biere, 2013]. Applying this strategy recursively on the d-separation encoding establishes that all clauses in the cone of influence of K′ are switched on iff all (in)dependence constraints in K′ are enforced to be satisfied during the incremental core extraction SAT calls.

6 Bounds-based Constraint Hardening

During search, we use the current bounds for boundsbased hardening of (in)dependence constraints. Our domainspecific approach is based on the following well-known technique of variable fixing for IP solving in operations research literature [Danzig et al., 1954; Crowder et al., 1983]. 2 By finding the min-cost hitting set over the current set of accumulated cores, conditioned on including a particular constraint k, we can obtain a conditional lower bound L¬k for the space of solutions in which k is not satisfied. If L¬k exceeds the current global upper bound U, k can be enforced as a hard constraint. Instead of solving the IP problem, a good bound L¬k can be obtained via the corresponding LP relaxation which is faster to solve. For hardening constraints early, this technique benefits crucially from the improved bounds obtained by incremental core extraction. Each hardened constraint simplifies the subsequent SAT and IP calls, and also speeds up incremental core extraction as only constraints that have not yet been hardened are considered there.

As a domain-specific extension, we also use this scheme to harden Boolean variables representing the causal graph solution, i.e., the existence of individual edges in the graph. Specifically, if the solution graph includes any type of an edge between nodes X and Y , then it leaves all constraints K′′ of the formX ⊥⊥ Y |S unsatisfied. If there are many sets S ⊆ V \ {X, Y } such that X ⊥⊥ Y |S, the sum of weights of such constraints can be large. Thus it makes sense to find the lower bound L′′ conditional on the existence of an edge between nodes X and Y-we can obtain such L′′ from the LP relaxation conditioned on including each constraint in K′′ . If L′′ is greater than the current global upper bound U, no edges can then be present between X and Y in an optimal graph and thus we can w.l.o.g. assign all edge variables between X and Y to be false. (Note that, in contrast, the inexact PC algorithm infers the absence of edges between X and Y upon finding a single (possibly unreliable) independence X ⊥⊥ Y |S [Spirtes et al., 2000].) Inferring hard absences of edges directly simplifies the subsequent SAT solver calls by tightening the solution space. The subsequently extracted cores are thereby often smaller, which makes the subsequent hitting set IPs also simpler.

7 Experiments

We implemented Dseptor by extending the MaxSAT solver LMHS with the domain-specific techniques proposed in Sects. 4-6; We use MiniSAT [E´en and S¨orensson, 2004] and IBM CPLEX as the internal SAT and IP solvers, respectively. We compare the performance of Dseptor to those of state-of-art MaxSAT solvers using the encoding of Hyttinen et al. [2014]. The experiments were run on 2.83-GHz Intel Xeon E5440 machines with 32-GB RAM and Debian GNU/Linux.

First we compare the performance of the solvers on synthetic data generated from 7-node cyclic (possibly) linear Gaussian models with correlated disturbances. The edges were sampled at random to obtain average degree 2, coefficients were from ±[0.2, 0.8]. We used BIC-based model selection to obtain the weights for the independence constraints over N = 1000 samples. Figure 4 gives a comparison of Dseptor with the MaxSAT solvers LMHS [Saikko et al., 2016], MaxHS [Davies and Bacchus, 2013], Maxino [Alviano et al., 2015], MSCG [Morgado et al., 2015], OpenWBO [Martins et al., 2014], QMaxSAT [Koshimura et al., 2012], and WPM3 [Ans´otegui et al., 2016], CPLEX run on a standard IP translation of MaxSAT [Davies and Bacchus, 2013; Ans´otegui and Gab`as, 2013], and the ASP solver Clingo (using its default branch-and-bound search) [Gebser et al., 2012] as used in [Hyttinen et al., 2014]. Dseptor outperforms all of the other solvers, solving 96% of the instances under the timeout of 600 seconds. Maxino comes in second (81%) ahead of plain LMHS and QMaxSAT. Clingo and many other solvers are not competitive here. Figure 4 also shows that all techniques proposed in Sects. 4-6 improve the performance: turning off each individual technique in Dseptor increases solving times.

2 Hardening using costs of residual formulas has also been used in pure SAT-based MaxSAT algorithms [Ans´otegui et al., 2013].

To examine the performances on more realistic sets of (in)dependence constraints, we looked at real-world datasets often used for benchmarking exact Bayesian network structure learning algorithms [Yuan and Malone, 2013; Bartlett and Cussens, 2017]. To comply with the current problem formulation, we considered suitable-sized (n) subsets of the variables in the datasets, the remaining variables becoming thus latent. We employed the BDEU score with equivalent sample size 10 to obtain constraint weights for this discrete data, and used a per-instance timeout of 7200 seconds. The results summarized in Figure 5, and detailed in Table 2 for Dseptor and its best competitors, LMHS (on which Dseptor is based), Maxino, and QMaxSAT, show that Dseptor clearly outperforms the competing approaches on these real-world benchmarks. Table 2 also gives the time taken by Dseptor to find an optimal graph (without proving its optimality yet) in parentheses. On a great majority of instances Dseptor finds an optimal solution very fast; the proof of optimality hence dominates the overall running times. This implies that Dseptor has very good anytime performance, being often able to give very good solutions in short time.