Bayesian Network Structure Learning with Integer Programming: Polytopes, Facets and Complexity (Extended Abstract)

1 Introduction

The study of probabilistic graphical models is a central topic in modern artificial intelligence research. Bayesian networks [Koller and Friedman, 2009] form a central class of probabilistic graphical models. A central problem related to Bayesian networks (BNs) is that of learning them from data. An essential part of this learning problem is aimed at learning the structure of a Bayesian network-represented as a directed acyclic graph-that accurately represents the (hypothetical) joint probability distribution underlying the data.

There are two principle approaches to Bayesian network learning: constraint-based and score-based. In the constraint-based approach [Spirtes et al., 1993; Colombo et al., 2012] the goal is to learn a network which is consistent with conditional independence relations which have been inferred from the data. The score-based approach to Bayesian network structure learning (BNSL)-focused on here-treats the BNSL problem as a combinatorial optimization problem of finding a BN structure that optimises a score function for given data.

∗This paper is an extended abstract of an article in the Journal of Artificial Intelligence Research [Cussens et al., 2017]. The authors gratefully acknowledge financial support from: UK Medical Research Council Grant G1002312 (JC, MB); Senior Postdoctoral Fellowship SF/14/008 from KU Leuven (JC); UK NC3RS Grant NC/K001264/1 (JC); Academy of Finland under grants 251170 COIN Centre of Excellence in Computational Inference Research, 276412, and 284591 (MJ); Research Funds of the University of Helsinki (MJ); and Icelandic Research Fund grant 152679-051 (JHK). Part of the work was done while JHK was at University of Helsinki and at Reykjavik University.

Learning an optimal BN structure is a computationally challenging problem: even the restriction of the BNSL problem where only BDe scores [Heckerman et al., 1995] are allowed is known to be NP-hard [Chickering, 1996]. Due to NP-hardness, much work on BNSL has focused on developing approximate, local search style algorithms [Tsamardinos et al., 2006] that in general cannot guarantee that optimal structures in terms of the objective function are found. Recently, despite its complexity, several advances in exact approaches to BNSL have surfaced [Koivisto and Sood, 2004; Silander and Myllym¨aki, 2006; Cussens, 2011; de Campos and Ji, 2011; Yuan and Malone, 2013; van Beek and Hoffmann, 2015], ranging from problem-specific dynamic programming branch-and-bound algorithms to approaches based on A∗ -style state-space search, constraint programming, and integer linear programming (IP), which can, with certain restrictions, learn provably-optimal BN structures with tens to hundreds of nodes.

As shown in a recent study [Malone et al., 2014], perhaps the most successful exact approach to BNSL is provided by the GOBNILP system [Cussens, 2011]. GOBNILP implements a branch-and-cut approach to BNSL, using state-of-the-art IP solving techniques together with specialised BNSL cutting planes. The focus of this work is on providing further understanding of the IP approach to BNSL from the theoretical perspective.

Viewed as a constrained optimisation problem, a central source of intractability of BNSL is the acyclicity constraint imposed on BN structures. In the IP approach to BNSL-as implemented by GOBNILP-the acyclicity constraint is handled in the branch-and-cut framework via deriving specialised cutting planes called cluster constraints which were originally introduced by Jaakkola et al. [2010]. These cutting planes are found by solving a sequence of so-called subIPs arising from solutions to linear relaxations of the underlying IP formulation of BNSL without the acyclicity constraint. Finding these cutting planes is an example of a separation problem for a linear relaxation solution, so called since the cutting plane will separate that solution from the set of feasible solutions to the original (unrelaxed) problem. Understanding fundamental aspects of these cutting planes and the sub-IPs used to find them is important not only from a purely theoretical perspective, but also since it holds out the promise of further improving the efficiency of state-of-the-art approaches to solving BNSL exactly. This is the focus of and underlying motivation for this work.

In this extended abstract we summarise the main results presented in [Cussens et al., 2017]. First, we study the computational complexity of the separation problem solved via sub-IPs. We establish that the sub-IPs are themselves NPhard to solve. From the practical perspective, this both gives a theoretical justification for applying an exact IP solver to solve the sub-IPs within GOBNILP, and motivates further work on improving the efficiency of the sub-IP solving via either improved exact techniques and/or further approximate algorithms.

Second, we study the facets of the convex hull of acyclic digraphs. As a key theoretical result, we show that cluster constraints are in fact facet-defining inequalities of this convex hull. From the more practical perspective, achieved via exhaustive computation, we provide a complete enumeration of facets for low-dimensional convex hulls. Mapping to practice, explicit knowledge of such facets has the potential to further speed up state-of-the-art BNSL solving by integrating (some of) these facets explicitly into the search.

2 BN Structure Learning

We focus on the Integer Programming (IP) approach to exact score-based Bayesian network structure learning. The essence of this technique is to encode the BNSL problem as an equivalent IP formulation. IP solving techniques can then be used to find an optimal solution to the problem and hence an optimal BN.

We restrict attention to decomposable score functions, where the score is defined locally by the parent set choices for each random variable, i, in the vertex set V . Specifically, for i ∈ V and J ⊆ V \ {i}, let i ← J denote the the pair (i, J), called a family. In our framework, we assume that the score function gives a local score ci←J for each family i ← J. A global score c(B) for each candidate structure (V,B) is then defined as

and the task is to find an acyclic digraph (V,B) maximising c(B) over all acyclic digraphs over V .

We allow for the possibility that not all J ⊆ V \ {i} are allowed as candidate parent sets for i, writing P(i) to denote the set of parent sets which are allowed. For example, we may only wish to consider parent sets of small cardinality either to reflect prior knowledge or to make solving easier. In fact, restricting parent set cardinality does not make BNSL easier in this sense: any BNSL problem can be converted into one where no parent set contains more than two parents, at the MAXIMISE

SUBJECT TO

Figure 1: The IP formulation of the BNSL problem.

expense of creating a polynomial number of auxiliary vertices [Cussens et al., 2017].

3 The GOBNILP System

GOBNILP [Cussens, 2011] is a program that has been developed to find optimal Bayesian networks by solving an IP formulation of the BNSL problem using the SCIP IP system [Achterberg, 2007]. The search is implemented as a branch-and-cut approach, the essentials of which are outlined next.

The IP formulation [Bartlett and Cussens, 2017] with which we work is shown in Figure 1. The binary IP variables used, xi←J, correspond to node i having the set of nodes J as parents in the network. If set to 1, it means that J is the parent set of i in the network; if set to 0, J is not the parent set of i.

The objective function follows from the global score. For any particular valid assignment of values to variables, the objective function sums the local scores for the chosen families (i.e. those where xi←J = 1) and ignores the others. This corresponds directly to the definition of the global score for that network, as given in Equation 1.

Equation 3 enforces the constraint that each node has exactly one parent set, while inequalities (4) prevent any solution in which there would exist a cycle in the resulting graph.

Inequalities (4) are just one way that cycles can be prevented [Cussens, 2010; Peharz and Pernkopf, 2012; Cussens et al., 2013]. These particular constraints are known as cluster constraints [Jaakkola et al., 2010] as each corresponds to the constraint that any cluster (set) of nodes must have at least one node with no parents in that cluster. An illustration of how acyclic graphs satisfy all cluster constraints but cyclic graphs do not is given in Figure 2.

GOBNILP begins with a relaxed version of the IP problem in which (i) the integrality constraints on all variables are removed, and (ii) all cluster constraints are removed. This relaxed problem is first solved, giving an optimal relaxed solution. In almost all cases, this relaxed solution is not a valid solution to the original IP problem due to violating cluster constraints

One or more cluster constraints that are violated by the relaxed solution are identified and added to the relaxed problem, and this modified relaxed problem is solved again.

The process of solving the relaxed problem and then adding violated cluster constraints repeats until the relaxed solution no longer breaks any cluster constraints. If the solution thus obtained assigns integer values (i.e. 0 or 1) to all variables, then an optimal solution to the problem has been found. Alternatively (and more commonly), if the relaxed solution contains one or more variables with fractional values then the solver branches on one of the variables with a fractional value, forming one subproblem in which this variable is fixed to 0 and one in which it is fixed to 1. The algorithm continues in this manner (cutting and then branching) for each of these problem instances.

4 Complexity of the GOBNILP Sub-IPs

As GOBNILP adds cluster constraints as cutting planes, it is necessary to detect which of these constraints are violated by the solution to the relaxed problem, or alternatively to determine that the solution respects all cluster constraints. In GOBNILP one of the techniques used to do this is through a sub-IP - a second IP that yields a solution that corresponds to a violated cluster constraint if one exists or is infeasible if the current solution obeys all cluster constraints.

Suppose that the solution to the current relaxed problem is x∗ . For each variable x∗i←J > 0, the sub-IP will contain an associated binary variable yi←J. If the solution to the sub-IP sets one of these variables to 1, it indicates that the corresponding variable in the original IP is included in the cluster constraint found. Additionally, for each i ∈ V , a binary variable yi is created. Variables of this type which are set to 1 in the solution indicate that the associated node is a member of the cluster which defines the cut.

As shown in [Cussens, 2011; Bartlett and Cussens, 2017], the sub-IP in Figure 3 either detects a cluster constraint to add as a cut if this is possible, or is infeasible if no such cut exists.

GOBNILP spends a considerable proportion of its run time solving this problem. It is therefore an important issue to understand the complexity of this problem. In particular, is there an algorithm with polynomial time complexity that can either detect a cluster constraint to add as a cutting plane or determine with certainty that none exists? A polynomial-time algorithm could potentially be used to significantly speed up

SUBJECT TO

Figure 3: The GOBNILP Sub-IP.

GOBNILP. While certain polynomial-time solvable cases can be identified - specifically, when x∗ is integer-valued and so represents a (possibly cyclic) directed graph - the complexity of solving the sub-IPs in the general case has been an open question.

We show that the task solved by the sub-IPs, which we call the weak separation problem for BNSL, is in fact NP-hard. Formally, we define the weak separation problem for BNSL as follows: given a solution x∗ to a relaxed problem, find a separating cluster C ⊆ V , |C| > 1, for which

or establish that no such C exists.

Theorem 1. The weak separation problem for BNSL is NPhard, even when restricted to instances where each parent set has size at most 2, that is, J ∈ P(i) for all i ∈ V only if |J| ≤ 2.

Theorem 1 is proven by a reduction from vertex cover. Informally, the main idea of the proof is to represent each edge {u, v} ∈ E in the original vertex cover instance G = (V, E) by assigning weight x∗s←{u,v} = 1/ |E| in the new instance, where s is a node not appearing in V , as illustrated in Figure 4. Clearly, U ⊆ V is a vertex cover in G if and only if the total weight of terms x∗s←{u,v} such that U intersects the parent set is 1. In other words, the summation constraint J ∩ C ̸= ∅ in inequality (11) is used to model the fact that a vertex cover is a hitting set for edges. Finally, we duplicate this construction for multiple new nodes s, and use terms x∗v←{s} to control the size of the vertex cover; for full details, see [Cussens et al., 2017].

5 BNSL Facets

Let G be the set of (family-variable encoded) acyclic digraphs for some vertex set V . PF(V, PV ) denotes the convex hull of G: the set of all convex combinations of elements of G, or equivalently the smallest convex set containing G. We call PF(V,PV ) the family-variable polytope 1 . It is a polytope with G as its vertices.

There is a set of linear inequalities which provide a minimal representation of the family-variable polytope. Such inequalities are known as facet-defining. To get an intuition of what a facet is, consider the cube with vertices

This is a 3-d polytope with 6 vertices and 6 facets defined by the 6 facet-defining inequalities

As we discuss in Section 6, when |V | ≤ 4 it is possible to compute the set of all facet-defining inequalities of the family-variable polytope. For bigger vertex sets this is infeasible, however it remains important to identify facet-defining inequalities since they make the best cuts. In [Cussens et al., 2017], we show that cluster constraints are facet-defining. In fact, the stronger result that all k-cluster constraints, as defined in [Cussens, 2011], are facet-defining is given.

Theorem 2. For any C ⊆ V and any positive integer κ < |C|, the valid inequality

is facet-defining for the family variable polytope PF(V, PV ).

We also show in [Cussens et al., 2017] that any facetdefining inequality for a vertex set V can be 'lifted' to generate a facet-defining inequality for a bigger vertex set V ′ :

Theorem 3. Let PV (i) := 2V \{i} for all i ∈ V . Let

1 In the full length article [Cussens et al., 2017] we also consider two other related polytopes: the digraph and the cluster polytopes.

be a facet-defining inequality of the family variable polytope PF(V, PV ) which is not a lower bound on a variable. Let V ′ be a node set such that V ⊆ V ′ . Then

is facet-defining for PF(V ′ , PV ′ ) and is not a lower bound on a variable.

Theorems 2 and 3 concern family-variable polytopes when all parent sets are allowed for each vertex. Since in practical BNSL it is often necessary (or desirable) to limit the choice of parent sets in some way, it is important to be able to determine when a facet-defining inequality remains facet-defining when some parent sets are disallowed. The following theorem of [Cussens et al., 2017] helps us do this.

Theorem 4. Let πx ≤ π0 define a facet for the familyvariable polytope PF(V,P). Suppose that πi←J = πi←J ′ for some i ∈ V , J, J ′ ∈ P(i) with J ⊊ J ′ , J ̸= ∅. Let ˘π be π with the component πi←J ′ removed. Let ˘ P be identical to P except that J ′ is removed from P(i). Then ˘πx ≤ π0 defines a facet for the polytope PF(V, ˘ P).

Given a facet-defining inequality of an all-parent-setsallowed polytope PF(V, PV ) and a parent set cardinality limit κ, Theorem 4 establishes that if the coefficients for all family variables xi←J ′ with |J ′ | > κ are not strictly larger than the coefficient for some family variable xi←J with J ⊊ J ′ so that |J| ≤ κ, then the inequality also defines a facet for the polytope with family variables restricted by κ. In [Cussens et al., 2017] this is confirmed computationally for the case where |V | = 4 and κ = 2. It follows that a normal (k = 1) cluster constraint is facet-defining for any limit κ on the size of parent sets.

6 Enumeration of Low-Dimensional Facets

Complementing the theoretical results, we provide a complete enumeration of facets for low-dimensional familyvariable polytopes. Mapping to practice, explicit knowledge of such facets has the potential to further speed up state-ofthe-art BNSL solving by integrating (some of) these facets explicitly into the search.

We provide a complete enumeration (see [Cussens et al., 2017, Appendix A] for details) of the facet-defining inequalities over 2-4 nodes and confirm the enumeration is consistent with the theoretical results. The enumeration was done using the cdd polyhedral computation software [Fukuda, 2016].

• For node set of size 2, i.e., V = {a, b}, there are 3 acyclic digraphs and three facet-defining inequalities: the two lower bounds xa←b ≥ 0 and xa←b ≥ o and the 1-cluster constraint xa←b + xb←a ≤ 1.
• For |V | = 3, there are 25 acyclic digraphs and and the convex hull of these 25 acyclic digraphs (i.e. the familyvariable polytope) has 17 facet-defining inequalities
• For |V | = 4, there are 543 acyclic digraphs, and 135 facet-defining inequalities of the family variable polytope.