{'fraction_non_alphanumeric': 0.04374470830141515, 'fraction_numerical': 0.01261944119662944, 'mean_word_length': 3.9717378232110643, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Turbocharging Treewidth-Bounded Bayesian Network Structure Learning

We present a new approach for learning the structure of a treewidth-bounded Bayesian Network (BN). The key to our approach is applying an exact method (based on MaxSAT) locally, to improve the score of a heuristically computed BN. This approach allows us to scale the power of exact methods---so far only applicable to BNs with several dozens of nodes---to large BNs with several thousands of nodes. Our experiments show that our approach outperforms a state-of-the-art heuristic method.


Introduction
Bayesian network structure learning is the notoriously difficult problem of discovering a Bayesian network (BN) that optimally represents a given set of training data [4]. Since the exact inference on a BN is exponential in the BN's treewidth [14], one is particularly interested in learning BNs of bounded treewidth. However, learning a BN of bounded treewidth that optimally fits the data (i.e., with the largest possible score) is, in turn, an NP-hard task [13]. This predicament caused the research on treewidth-bounded BN structure learning to split into two branches: 1. Heuristic Learning (see, e.g., [5,17,23,24]), which is scalable to large BNs with thousands of nodes (but with a score that can be far from optimal), and 2. Exact Learning (see, e.g., [2,13,19]), which learns optimal BNs (but is scalable only to a few dozen nodes).
In this paper, we combine heuristic and exact learning and take the best of both worlds. The basic idea for our approach is to first compute a BN with a heuristic method (the global solver), and then to apply an exact method (the local solver) to parts of the heuristic solution. The parts are chosen small enough that they allow an optimal solution reasonably quickly with the exact method. Although the basic idea sounds compelling and reasonably simple, its realization requires several conceptual contributions and new results.
For the global solver, any heuristic algorithm for treewidth-bounded BN learning, such as the robust k-MAX algorithm by Scanagatta et al. [24], can be used. The task for the local solver is significantly more complex than treewidth-bounded BN structure learning, as several additional constraints need to be incorporated.
It is not sufficient that the BN computed by the local solver is acyclic. We need fortified acyclicity constraints that prevent cycles that run through the other parts of the BN, which have not been changed by the local solver. Similarly, it is not sufficient that the local BN is of bounded treewidth. We need fortified treewidth constraints that prevent the local BN from introducing links between a diverse set of nodes that, together with the other parts of the BN, which have not been changed by the local solver, increase the treewidth.
Given these additional requirements, we propose a new local solver BN-SLIM (for SAT-based Local Improvement Method), which satisfies the fortified constraints. We formulate a fortified version of the treewidth-bounded BN structure learning problem. In Theorem 1, we show that we can express the fortified constraints with the addition of certain virtual arcs and virtual edges. The virtual arcs represent directed paths that run outside the local instance, with these virtual arcs we can ensure fortified acyclicity. The virtual edges represent essential parts of a global tree decomposition using which we can ensure bounded treewidth.
The new formulation of the local problem is well-suited to be expressed as a MaxSAT (Maximum Satisfiability) problem and hence allows us to harvest the power of state-of-the-art MaxSAT solvers. A distinctive feature of the encoding is that, in contrast to the virtual edges, the virtual arcs are conditional and depend on the solution computed by the local solver.

Results
We implement BN-SLIM and evaluate it empirically on a large set of benchmark data sets, consisting between 64 and 10,000 nodes and for treewidth bounds 2, 5, and 8. As mentioned above, we use k-MAX [24] as the global solver. k-MAX improves over the k-greedy algorithm [23], which was the first algorithm for treewidth-bounded structure learning that scaled to thousands of nodes. k-MAX is an anytime algorithm that provides better and better solutions over time. However, as time progresses, improvements become less and less frequent. In our tests, we noticed that around 95% of all improvements occurred in the first 30 minutes when running for an hour. This is a good point to hand it over to BN-SLIM, which acts as a turbocharger. Our results show that BN-SLIM is capable of improving the global solution provided by k-MAX in nearly all cases; the improvement is often dramatic.
In our experimental setup, we split a total timeout of 60 minutes into an initial phase of 30 minutes, where the global solver computes the heuristic solution. In the remaining 30 minutes, we perform many local improvements with the local solver BN-SLIM. For comparison, we also record as the baseline, the solution found by k-MAX if left running for the remaining 30 minutes. We compare the final solution obtained by BN-SLIM with the baseline solution computed by k-MAX. According to the ∆BIC metric, which was used by Scanagatta et al. [24] for comparing treewidth-bounded BN structure learning algorithms, the results are "extremely positive" in favor of BN-SLIM in a vast majority of the experiments.

Related work
The first SAT-encoding for finding the treewidth of a graph was proposed by Samer and Veith [21]. Fichte et al. [6] proposed the first SAT-based local improvement method for treewidth, using the Samer-Veith encoding as the local solver. Recently, SAT encodings have been proposed for other graph and hypergraph width measures [7,8,15,25].
Several exact approaches to treewidth-bounded BN structure learning have been proposed. Korhonen and Parviainen [13] proposed a dynamic-programming approach, and Parviainen et al. [19] proposed a Mixed-Integer Programming approach. Berg et al. [2] proposed a MaxSAT approach by extending the basic Samer-Veith encoding for treewidth. Our approach for BN-SLIM uses a similar general strategy, but we encode acyclicity differently. Moreover, BN-SLIM deals with the fortified constraints in terms of virtual edges and virtual arcs.
Since the exact methods are limited to small domains, Nie et al. [17,18] suggested heuristic approaches that scale up to some hundreds of variables. The k-greedy algorithm proposed by Scanagatta et al. [23] at NIPS'16 provided a breakthrough, consistently yielding better DAGs than its competitors and scaling up to several thousands of variables. As mentioned above, k-MAX [24] is a more recent improvement over k-greedy.

Structure learning
We consider the problem of learning the structure (i.e., the DAG) of a BN from complete data set of N instances D 1 , . . . , D N over a set of n categorical random variables X 1 , . . . , X n . The goal is to find a DAG D = (V, E) where V is the set of nodes (one for each random variable) and E is the set of arcs (directed edges). The value of a score function determines how well a DAG D fits the data; D, together with local parameters, forms the BN.
We assume that the score is decomposable, i.e., being constituted by the sum of the scores of the individual variables. Hence we can assume that the score is given in terms of a local score function f that assigns each node v ∈ V and each subset P ∈ P v a real number f P (v). The score of a DAG D = (D, E) is then This setting accommodates several popular scores like BDeu, BIC, and AIC [1,11,26]. Notice that if a parent set of v has a score smaller than or equal to the score of the empty parent set, then we can safely disregard this parent set.

Treewidth
Treewidth is a graph invariant that provides a good indication of how costly probabilistic inference on a BN is. Treewidth is defined on undirected graphs and applies to BNs via the moralized graph where T is a tree and χ is a function that assigns each tree node t a set χ(t) of vertices of G such that the following conditions hold: T1 For every edge e of G there is a tree node t such that e ⊆ χ(t).
T2 For every vertex v of G, the set of tree nodes t with v ∈ χ(t) induces a non-empty subtree of T .
The sets χ(t) are called bags of the decomposition T , and χ(t) is the bag associated with the tree node t. The width of a tree decomposition (T, χ) is the size of a largest bag minus 1. The treewidth of G, denoted by tw(G), is the minimum width over all tree decompositions of G.
The treewidth-bounded BN structure learning problem takes as input a set V of nodes, a local score function f on V , and an integer W , and it asks to compute a DAG D = (V, E) of treewidth ≤ W , such that f (D) maximal.

Local improvement
Consider an instance (V, f, W ) of the treewidth-bounded BN structure learning problem, and assume we have computed a solution D = (V, E) heuristically, together with a tree decomposition T = (T, χ) of width ≤ W of the moralized graph M (D).
We select a subtree S ⊆ T such that the number of vertices in V S := t∈V (S) χ(t) is at most some budget B. The budget is a parameter that we specify beforehand, such that the subinstance induced by V S is small enough to be solved optimally by an exact method, which we call the local solver. The local solver computes for each v ∈ V S a new parent set, optimizing the score of the resulting DAG The local solver ensures that the following conditions are met: We assume that the local solver certifies C2 by computing a tree decomposition The two conditions stated above are not sufficient to ensure that D new is acyclic and that treewidth of M (D new ) remains bounded by W . Acyclicity can be violated by cycles formed by the combination of the new incoming arcs of vertices in S together with old arcs that are kept from D. The treewidth can increase by a number that is linear in |V S |.
Hence, we need additional side conditions, which we will formulate using the following additional concepts.
Let us call a vertex v ∈ V S a boundary vertex if there exists a node t ∈ V (T ) \ V (S) such that v ∈ χ(t), i.e., it occurs in some bag outside S. We call the other vertices in V S internal vertices, and the vertices We can now formulate the side conditions.
C3 S new is a tree decomposition of the extended moral graph M ext .
is acyclic. We note that condition C4 implies that in D new , all parents of an internal vertex are internal. . . , T r be the connected components of T \ V (S), i.e., the T i 's are the subtrees of T that we get when deleting the subtree S. Let V i = t∈V (Ti) χ(t), 1 ≤ i ≤ r, and observe that each external vertex x belongs to exactly one of the sets V 1 , . . . , V r . Let B i = V S ∩ V i , 1 ≤ i ≤ r, be the set of boundary vertices in V i . We observe that all the vertices in B i are mutually adjacent boundary vertices and occur together in a bag χ(s i ) of s i ∈ V (S) and in a bag χ(t i ), for t i ∈ V (T i ), as we can take s i and t i to be the two neighboring tree nodes of T with s i ∈ V (S) and t i ∈ V (T i ). We also observe that each B i forms a clique in the extended moral graph M ext .
Recall that by assumption, the local solver provides a tree decomposition Additionally, by condition C3, S new is also a tree decomposition of M ext , and hence, by a basic property of tree decompositions (see, e.g., [3, Lem. 3.1]), there must exist a bag Hence we can define T new as the tree we get by connecting the disjoint trees S new , T 1 , . . . , T r with the edges To prove the claim, we show that T new satisfies the properties T1 and T2.
If neither u nor v is external, then (u, v) ∈ E new S , and since S new is a tree decomposition of D new S , {u, v} ⊆ χ new (s) for some s ∈ V (S new ). If v is external but u isn't, then the arc (u, v) was already present in D, as the parents of external vertices didn't change. Hence, since T is a tree decomposition of M (D), ∈ V S , then u, v are two adjacent boundary vertices, hence {u, v} is a virtual edge which, by C3, means {u, v} ⊆ χ new (s) for some s ∈ V (S new ). We conclude that T1 holds.
Condition T2. Let v ∈ V . If v is external, then there is exactly one i ∈ {1, . . . , r}, such that v ∈ V i = t∈V (Ti) χ(t). Since we do not change the tree decomposition of T i , condition T2 carries over from T to T new . Similarly, if v is internal, then v does not appear in any bag χ new (t) for t ∈ V (T new ) \ V (S new ), hence condition T2 carries over from S new to T new . It remains to consider the case where v is a boundary vertex. The nodes t ∈ V (S new ) with v ∈ χ(t) are connected, because S new satisfies T2, and for . Hence all the nodes t ∈ V (T new ) with v ∈ χ(t) are connected, and T2 also holds for boundary vertices.
To conclude the proof of the claim, it remains to observe the width of T new cannot exceed the widths of T or S new , hence the width of T new is at most W .
We observe that the vertices on the path from x j to x j on C all belong to some V i = t∈V (Ti) χ(t). Hence v j and v j are adjacent boundary vertices, and E → virt contains all the arcs (v j , v j ), 1 ≤ j ≤ p. However, the cycle C contains also the paths from v j to v j+1 (mod p) , for 1 ≤ j ≤ p, which only run through vertices in V S . These paths, together with the virtual arcs (v j , v j ) form a cycle C which lies in (V S , E new S ∪ E → virt ). This contradicts C5 which requires that this digraph be acyclic. Hence the claim holds. We observe that by taking D new = D we have a solution that satisfies all the required conditions and maintains the score.

Implementing the local improvement
In this section, we first discuss how the set S representing the subinstance is constructed. Then we provide a detailed explanation of the MaxSAT encoding that is responsible for solving the subinstance.

Constructing the subinstance
For this section, we follow the same notation as used in the previous section. To construct the subinstance, we initialize the subtree S with a node r picked at random from V (T ). We then expand S by performing a bread-first search from V (S) and adding a new node to S as long as the size of V S does not exceed the budget. Next, we compute E virt for the chosen S. Finally, we prune the parent sets of each vertex so as to only retain those parent sets which satisfy conditions C3 and C4. This can be done by first checking, for each parent set, if the required node t is present V (T ) \ V (S), and if it does, we record the set of virtual arcs that are imposed by this parent set as long as none of the virtual arcs are self-loops. For each v ∈ S and P ∈ P v , we denote by A → virt (v, P ) the set of imposed virtual arcs when v has the parent set P in D new .
We denote by P v , the collection of parent sets of node v that remain after this pruning process. Notice that, under this pruning, all remaining parent sets P ∈ P v satisfy C4. Also note that, since E → virt is conditional on the chosen parent sets, it cannot be precomputed.
Further, since we intend to solve the subinstance using a MaxSAT encoding, we need to ensure that the score of each parent set is non-negative. Recall that P v only contains those non-empty parent sets whose score is at least that of the empty parent set. Thus, we may assume that the empty parent set has the lowest score among all the parents of a certain vertex. Consequently, we can adjust the score function by setting for v ∈ S and P ∈ P v , which implies that f P (v) ≥ 0 for all v ∈ S and P ∈ P v .

MaxSAT encoding
We now describe the weighted partial MaxSAT instance that encodes conditions C1-C5. We build on top of the SAT encoding proposed by Samer and Veith [21]. The only difference in our case is that there are no explicit edges and hence we do not require the corresponding clauses. Instead, the edges of the moralized graph are dependent on and decided by other variables that govern the DAG structure. For convenience, let n denote the size of the subinstance, i.e., n := |S|. A part of the encoding is based on the elimination ordering of a tree decomposition (refer to e.g., [21,Sec. 2]).
The main variables used in our encoding are • variables par P v represent for each node v ∈ S the chosen parent set P , • n(n − 1)/2 variables acyc u,v represent the topological ordering of D new S , • n(n − 1)/2 variables ord u,v represent the elimination ordering of the tree decomposition, • n 2 variables arc u,v represent the arcs in the moralized graph M ext , along with the fill-in edges.
Since acyc u,v and ord u,v represent linear orderings, we enforce transitivity of these variables by means of the clauses To prevent self-loops in the moralized graph, we add the clauses For each node v ∈ S, and parent set P ∈ P v , the variable par P v is true if and only if P is the parent set of v. Since each node must have exactly one parent set, we introduce the cardinality constraint P ∈Pv Next, for each node v, parent set P , and u ∈ P , if P is the parent set of v then u must precede v in the topological ordering. Hence we add the clause par P v → acyc u,v for v ∈ S, P ∈ P v , and u ∈ P.
Similarly, for each node v, parent set P , and u ∈ P , if P is the parent set of v then we must add an arc in the moralized graph respecting the elimination ordering between u and v, as follows for v ∈ S, P ∈ P v , and u ∈ P.
Next, we encode the moralization by adding an arc between every pair of parents of a node, using the following clauses for v ∈ S, P ∈ P v , and u, w ∈ P.
Now, we encode the fill-in edges, with the following clauses Lastly, to bound the treewidth, we add a cardinality constraint on the number of outgoing arcs for each node as follows To complete the basic encoding, for every node v ∈ S, and every parent set P ∈ P v we add a soft clause weighted by the score of the parent set as follows To speed up the solving, we encode that for every pair of nodes, at most one of the arcs between them can exist. We add the following redundant clauses Now, we describe the additional clauses required to satisfy the fortified constraints, and thus conditions C3 and C5. For every virtual edge {u, v} ∈ E virt , we introduce a forced arc depending on the elimination ordering using the following pair of clauses This takes care of the fortified treewidth constraints, satisfying C3 and ensuring that the edge {u, v} ⊆ χ(s) for some s ∈ V (S). Finally, we add the clauses that encode the forced arcs E → virt . For each v ∈ S, P ∈ P v , and (u, v) ∈ A → virt (v, P ), we add the clause which forces the virtual arc (u, v) if P is the parent set of v in D new , thereby handling the fortified acyclicity constraints and ensuring that C5 is satisfied. This concludes the definition of the MaxSAT instance, to which we will refer as Φ D,f (S). We refer to the weight of a satisfying assignment τ of In particular, let K 0 be the weight of the truth assignment which corresponds to the parent sets of S as defined by the input DAG D. We summarize these observations in the following theorem.

Experiments
In this section, we describe the experiments conducted to analyze the performance of the local improvement algorithm. The current state-of-the-art heuristic algorithm for solving the treewidth-bounded BN structure learning problem is the k-MAX algorithm by Scanagatta, et al. [24], therefore, we compare our BN-SLIM approach with k-MAX. Accordingly, we found it fair to follow closely the experimental setup (including data sets, timeouts, comparison metrics) used by the authors of k-MAX to compare it with previous approaches.

Experimental Setup
We run all our experiments on a 4-core Intel Xeon E5540 2.53 GHz CPU, with each process having access to 8GB RAM. We use UWrMaxSat as the MaxSAT-solver primarily due to its anytime nature (available at the 2019 MaxSAT Evaluation webpage 1 ). We tried other solvers but found that UWrMaxSat works best for our use case. We use the BNGenerator package [12] in conjunction with the BBNConvertor tool [9] to generate and reformat random Bayesian Networks. We also make use of the implementation of the k-MAX algorithm available as a part of the BLIP package [22]. We implement the local improvement algorithm in Python 3.6.9, using the NetworkX 2.4 graph library [10]. The source code is available at https://www.ac.tuwien.ac.at/files/resources/software/bnslim.zip   andes  223  r0  2000  r5  4000  r10  10000  diabetes  413  r1  2000  r6  4000  r11  10000  pigs  441  r2  2000  r7  4000  r12  10000  link  724  r3  2000  r8  4000  r13  10000  munin  1041  r4  2000  r9  4000  r14  10000   Table 2: Synthetic data sets, n denotes the number of variables

Data sets
We consider 18 real data sets which have been previously used in [24], introduced in [16,27], and are publicly available 2 . These data sets are split into three subsets yielding 54 inputs. We also use 20 synthetic data sets, the first five of which are commonly used in literature as benchmarks 3 . The remaining 15 data sets are generated randomly using the BNGenerator tool, containing more variables than all the other data sets. Lastly, 5000 samples are drawn from each of these synthetic data sets. Thus we use a total of 74 data sets. Overall, the collection of data sets provides a broad variety of the different parameters and the nature of the data itself. We precompute the parent set score functions for all the data sets using independence selection (available in the BLIP package), and this cache is used as the input to both BN-SLIM and k-MAX.

Evaluation
We first conduct a preliminary analysis to find out the best values for the different parameters of BN-SLIM, which leads us to use a budget of 10 and a timeout of 2s per MaxSAT call. We run both BN-SLIM and k-MAX on the data sets for a total of one hour under the same hardware configuration. Further, since BN-SLIM needs an initial heuristic solution, we enlist k-MAX for this purpose. We denote by BN-SLIM(k-MAX), the algorithm which applies BN-SLIM on an initial solution provided by k-MAX. Out of the total one hour allotted to BN-SLIM(k-MAX), 30 minutes are used to generate the initial heuristic solution, and the remaining 30 minutes are used by BN-SLIM(k-MAX) to improve this initial solution.
For evaluating the performance of our algorithm, we use the same metric used by Scanagatta, et al., i.e., ∆BIC, which is the difference between the BIC scores of two solutions. Given graph D, the BIC score approximates the logarithm of the marginal likelihood of D. Thus, given two graphs D 1 and D 2 , the difference in their BIC scores approximates the ratio of their respective marginal likelihoods which is the Bayes Factor (BF). A positive ∆BIC score signifies positive evidence towards D 1 and a negative ∆BIC score signifies positive evidence towards D 2 . The values of ∆BIC can be mapped to a scale [20,Sec. 4.3] such that |∆BIC| > 10 counts as extremely positive/negative evidence, 10 > |∆BIC| > 6 counts as strongly positive/negative evidence, 6 > |∆BIC| > 2 counts as positive/negative evidence, and |∆BIC| ≤ 2 provides insignificant evidence and is labeled as neutral. The ∆BIC values of BN-SLIM(k-MAX) against k-MAX over all the data sets and three different treewidths are shown in Fig. 1.

Results
It can be observed from Fig. 1, that BN-SLIM(k-MAX) secures extremely positive evidence for significantly many datasets, across all tested treewidths, with the smaller treewidths being more favorable. The experiments demonstrate the effectiveness of the BN-SLIM approach and the combined power as a heuristic method of BN-SLIM(k-MAX) i.e., BN-SLIM on top of k-MAX.

Conclusion
With BN-SLIM, we have presented a powerful method for improving the outcome of treewidth-bounded BN structure learning heuristics. We have demonstrated the robustness and performance by applying BN-SLIM to the solution provided by the state-of-the-art heuristic k-MAX. Due to the complementary nature of the two approaches k-MAX (sampling) and BN-SLIM (reasoning), both methods work remarkably well together and exhibit the best of both worlds. Taking a 60-minute relay race, k-MAX starts running, but after 30 minutes, it runs out of steam, and improvements become rare. It then passes the baton to BN-SLIM, which now uses the remaining time for significant improvements and brings a highly optimized solution over the finish line. The experiments demonstrate that this collaborating team outperforms a single 60-minute k-MAX run.
The highly encouraging experimental outcome suggests several avenues for future work, which include the development of more sophisticated subinstance selection schemes, the inclusion of variable fidelity sampling (crude for global solver, fine-grained for local solver), as well as more complex collaboration protocols between local and global solver in a distributed setting.