From Support Propagation to Belief Propagation in Constraint Programming (Extended Abstract)

1 Introduction

The distinctive driving force of Constraint Programming (CP) to solve combinatorial problems has been a privileged access to problem structure through the high-level models it uses. From that exposed structure in the form of so-called global constraints, powerful inference algorithms have shared information between constraints by propagating it through shared variables' domains, traditionally by removing unsupported values. The paper investigates a richer propagation medium made possible by recent work on model counting inside constraints.

Many forms of message passing algorithms have been investigated in AI. One of the earliest and best known, Belief Propagation (BP) also known as sum-product message passing, was proposed by Pearl [1982] to perform inference on graphical models. From a joint probability distribution it computes approximations of the marginal distributions onto individual variables (i.e. beliefs). Each message is a realvalued function over the domain of a variable that expresses a probability that the variable takes a given value in a model. BP is known to converge to the exact marginal distributions on tree topologies but may not converge in general.

In CP, constraint propagation can also be viewed as a message passing algorithm, announcing value deletions from domains through the constraint network and necessarily converging because some quantity, namely the size of the Cartesian product of the domains, decreases monotonically [Horsch and Havens, 2013; Werner, 2015]. From the perspective of message passing, the deleted values being propagated are messages taking the form of simpler Booleanvalued functions evaluating to false for these deleted values and to true otherwise, which is rather flat information since all non-deleted values are on an equal footing. One way to view a variable's filtered domain with respect to a constraint is as a set of variable-value pairs having non-zero frequency among its solution set [Dechter et al., 2010] - if instead we share the whole frequency distribution over individual variables we can discriminate between values from the perspective of each constraint and, for example, use it for branching in the search tree. Once we consider the frequency of a variable-value pair as its likelihood, or probability, of appearing in a solution to the given constraint, we come very close to belief propagation and other message passing algorithms for probabilistic inference. Beliefs about individual variablevalue assignments are exchanged between contraints and iteratively adjusted. It generalizes standard support propagation in CP and aims to converge to the true marginal distributions of the solutions over individual variables. Its advantage over standard belief propagation is that the higher-level models of CP featuring large-arity (global) constraints do not tend to create as many cycles, which are known to be problematic for convergence.

∗The full version of this paper was published as [Pesant, 2019].

Since in general we don't have an explicit description of the solution set of a constraint, that frequency distribution is not readily available and would have to be approximated, perhaps very coarsely. However recent work on solution counting is bringing it within reach to compute exact (or close) distributions for several families of constraints [Pesant, 2017]. Counting-based search [Pesant et al., 2012], a family of effective branching heuristics, drives search toward areas where there are many satisfying assignments to (some of) the constraints taken individually. These are expressed through the concept of solution density: given a constraint c(x1, . . . , xk), its number of solutions #c(x1, . . . , xk), respective finite domains D(xj) 1≤j≤k, a variable xi in the scope of c, and a value v ∈ D(xi),

defines the solution density of variable-value pair (xi, v) in constraint c, given as the fraction of solutions to c which include that assignment. Depending on the constraint and the combinatorial structure it encapsulates, such counting may be intractable. Nevertheless the recent work on countingbased search has provided efficient algorithms to count either exactly or approximately for several important families of constraints [Pesant et al., 2012; Brockbank et al., 2013; Pesant, 2015].

There is a well-established connection between probabilistic inference and weighted model counting [Chavira and Darwiche, 2008]. From the perspective of weighted counting, the concept of solution density has an underlying assumption that the likelihood of values in each domain follows a uniform distribution, thereby leading to all solutions having equal weight. We need to abandon this assumption and consider arbitrary distributions arising from the variable-to-constraint messages that reflect the belief acquired from the other constraints and which is iteratively adjusted.

The paper promotes a richer propagation medium in CP, made possible by extending the work on solution density to weighted counting inside constraints and close in spirit to message passing algorithms. It also contributes an instantiation of belief propagation that is less affected by cycles through the use of higher-order potentials corresponding to CP's global constraints, weighted-counting algorithms for some of the latter, and a publicly-available implementation for further research.

2 An Example

Consider the following example to illustrate our approach:

i. alldifferent(a, b, c) ii. a+ b + c +d = 7 iii. c ≤ d

three constraints with variables a, b, c, d ∈ {1, 2, 3, 4}. There are two solutions to that CP problem: (a = 2, b = 3, c = 1, d = 1) and (a = 3, b = 2, c = 1, d = 1). Even if we enforce domain consistency on each constraint, no filtering occurs. To solve it we would thus be left to branch on variables having identical domains.

Variable a takes value 2 in one of the two solutions, value 3 in one solution as well, and values 1 and 4 in no solution. If we look at the set of solutions as a multivariate discrete distribution, its projection onto a yields 〈0, 1, 1, 0〉 and under the assumption that either solution is equally likely the marginal probability distribution for a is then θa = 〈0, 1/2, 1/2, 0〉. Table 1 gives the marginal for each variable and Table 2, the marginals local to each constraint taken individually and over its own set of solutions. We see that from the point of view of the alldifferent constraint and for variable a (line

d

Table 3: Computed marginals after ten iterations

θia ) each value is equally likely since it appears in the same number of solutions to that constraint. Whereas for that same variable but from the point of view of the linear equality constraint (line θii a ) value 1 is ten times more likely than value 4. Note also that for variable d the two local marginals give conflicting beliefs.

Table 3 presents the computed marginals after ten iterations of belief propagation on this CP model using exact weighted counting on each constraint. Note how close these got to the true marginals in Table 1. To motivate our interest in taking advantage of the high-level modelling typically present in CP, made up of a relatively low number of high-arity constraints, consider Table 4 showing the impact of replacing the alldifferent constraint by its decomposition into three disequality constraints: a ̸= b, a ̸= c, b ̸= c. The computed marginals are much further from the true marginals.

3 Belief Propagation in a CP Solver

Our prototype 1 is built on top of MiniCP, a recent bare-bones open-source CP solver developed for academic purposes and written in Java [Michel et al., 2017]. Though few constraints and filtering algorithms are currently implemented, its small and clean architecture made it easier to implement the required architectural changes to the core of the solver without worrying about the potential impacts on a more complex system. We briefly describe the main changes.

In our setting a message corresponds to the frequency (marginal) distribution of a variable. In order to maintain the marginal distribution of a variable, each value in its domain is given a weight attribute. Variables are provided with simple methods to receive and send messages: incoming messages from constraints are combined by multiplying them and an outgoing message sent to a constraint excludes its own contribution to the distribution (this is termed the outside belief ).

1 Its current implementation is available at https://github.com/ PesantGilles/MiniCPBP

Constraints are also provided with methods to receive and send messages, and another method to perform the weighted counting.

Message passing is synchronized in two phases: in the first phase constraints receive messages from the variables and then update their beliefs by performing weighted counting based on these updated marginal distributions; in the second phase marginals for the variables are reset, constraints send back messages to the variables, and marginals are normalized. This goes on for a fixed number of iterations.

We solve a problem by backtrack search as usual. At each node of the search tree we perform belief propagation and then branch based on the computed marginals.

4 WeightedModel Counting

The main technical challenge of this approach is to design efficient weighted model counting algorithms for each family of constraints. Because each has a different combinatorial structure, the kind of algorithm required likely differs as well. We present some of them here.

4.1 Table

The table(x1, x2, . . . , xk, T) constraint restricts the sequence of variables 〈x1, x2, . . . , xk〉 to take on a combination of values appearing in T , the set of allowed tuples. The Compact Table algorithm [Demeulenaere et al., 2016] implementing the table constraint computes bit sets Txv ⊆ T identifying for each variable-value pair (x, v) the tuples in which it appears. It also maintains a bit set S of currently-supported tuples (those that could be a solution to the constraint given the current domains of the variables). Counting the solutions in which x = v then simply amounts to computing the cardinality of Txv ∩ S. For weighted counting we first compute the weight of each currently-supported tuple 〈v1, v2, . . . , vk〉 as the product of the outside belief of each pair (xi, vi) and sum the relevant weights for each Txv.

4.2 Sum and Regular

The sum(x1, x2, . . . , xk, y) constraint enforces y = ∑ki=1 xi. It can be used to state inequalities and weighted sums as well. As previously proposed to achieve domain consistency and to compute solution densities for knapsack constraints [Pesant and Quimper, 2008], conceptually we build a layered graph where each path from the first to the last layer represents a solution, with each component arc corresponding to an individual variable assignment. The only difference here is that now each arc carries a weight equal to the outside belief of the corresponding variable assignment and that instead of storing the number of incoming and outgoing paths at each node, we store weighted sums πin and πout of incoming and outgoing paths respectively where a path's weight is the product of the weights on its component arcs. We compute them through standard forward and backward passes and then use these weighted sums to compute the local beliefs for each variable-value pair.

The regular(x1, x2, . . . , xk,A) constraint restricts the sequence of variables 〈x1, x2, . . . , xk〉 to spell out a word belonging to the regular language recognized by automaton A. The same approach can be used for regular constraints since a similar layered graph has been proposed to achieve domain consistency and to count solutions [Pesant et al., 2012].

4.3 Among

The among(x1, x2, . . . , xk, V, c) constraint makes variable c correspond to the number of occurrences of values from set V in {x1, x2, . . . , xk}. By introducing binary indicator variables y1, y2, . . . , yk such that yi = 1 ⇔ xi ∈ V we rewrite an among constraint as sum(y1, y2, . . . , yk, c) and use its weighted counting algorithm.

4.4 Alldifferent

The alldifferent(x1, x2, . . . , xk) constraint restricts variables x1, x2, . . . , xk to take on distinct values. The one-to-one correspondence between solutions to an alldifferent constraint and maximum matchings in the associated bipartite graph is well known. Counting such matchings in turn corresponds to computing the permanent of the adjacency matrix A = (aij) of the bipartite graph,

where P denotes the set of all permutations of {1, 2, . . . , k}. Unfortunately this is a #P-complete problem (one of the socalled "hard to count - easy to decide" problems) [Valiant, 1979]. Nevertheless it is a well-studied problem, even for more general nonnegative matrices, for which several upper bounds have been proposed [Soules, 2003] and even polytime randomized approximation algorithms [Jerrum et al., 2001]. We will cast our weighted counting problem as that of computing the permanent of a nonnegative matrix built from the marginals.

We first construct matrix A by setting aij to the outside belief about variable xi taking value j. We interpret the permanent of A as a weighted counting of maximum matchings: if aij = 1 its participation in a matching is fully counted; the smaller aij is, the more it will be discounted; aij = 0 still means that no matching includes the corresponding assignment. For each variable xi and value j, the (unnormalized) local belief is equal to the weighted counting of matchings that include assignment xi = j, which we evaluate as the permanent of the sub-matrix of A obtained by removing row i and column j. We compute that permanent exactly for sub-matrices whose size does not exceed a given small threshold τ and otherwise use a quick-to-compute upper bound (bound U 3 from [Soules, 2003]) thus getting an approximate weighted counting.

mindom

1000 10000

5 Combinatorial Search

We investigate CP-based belief propagation's search guidance to find a solution to constraint satisfaction problems. We use the following two-way branching heuristic: max-strength assigns the variable-value pair exhibiting the largest positive difference between its marginal and the reciprocal of the domain size (i.e. its marginal strength) and disallows it upon backtracking. As a baseline we use minimum domain size variable ordering with random value ordering and report the average result over ten runs.

We report on three combinatorial problems from the XCSP3 repository 2 that can be modelled using the constraints from the previous section: Latin Square, Magic Square, and Primes. We stop an individual run after one hour of computation time.

5.1 Primes

The CP model for the Primes instances features 100 variables and many sum constraints. Figure 1 plots the number of instances solved against the number of failed search-tree nodes reached during search (left) and against the computation time (right). Each curve represents a particular configuration of the algorithm: mindom is our baseline branching heuristic (minimum domain size with random value ordering) with standard support propagation; all the others use the max-strength branching heuristic, with SBP/BP indicating that belief propagation is used with/without prior application of support propagation and the integer label (1, 2, 5, 10) representing the number of BP iterations used at each search-tree node.

We start with the number of fails. The SBP curves show that search guidance improves with the number of BP iterations, as expected. Remarkably for most instances the maxstrength branching heuristic guides directly to a solution, culminating at 10 iterations with 26 out of 32 instances being solved backtrack-free. Branching heuristic mindom does not guide nearly as well and only solves 4 instances backtrackfree. In general it exhibits several orders of magnitude more fails. Turning to computation time (right plot) we see that the gap between mindom and max-strength almost closes because of the computational cost of exact counting over the many sum constraints but the latter nevertheless solves a few more instances within the time limit.

2 http://www.xcsp.org/

instances solved

40

35

30

25

20

15

10

5

0

1E-2

1E+0

1E+2

1E+4

1E+6

fails

fails SBP 1

SBP 2

SBP 5

SBP 10

BP 2

BP 5

BP 10

mindom

1

10

100

1000

time (s)

Figure 2: Number of instances solved with respect to the number

of fails (left) and computation time (right) for the Magic Square in-

stances.

5.2 Magic Square

The CP model for these 9 × 9 partially-filled instances features sum and alldifferent constraints. Figure 2 shows that heuristic mindom requires on average tens of thousands of fails to solve any of them but manages to solve almost all of them given one order of magnitude more fails. SBP with 5 and 10 iterations requires a few orders of magnitude fewer fails than mindom. Because weighted counting on the alldifferent constraint is not exact until few unbound variables remain, not applying support propagation (the BP curves) deteriorates performance. Turning to computation time, even though mindom eventually catches up on the last few instances solved, for the most part it is about one order of magnitude slower than max-strength.

5.3 Latin Square

The CP model for these 30 × 30 partially-filled instances features alldifferent constarints. Figure 3 shows that heuristic mindom only manages to solve 17 out of 40 instances within the time limit on average. Again maxstrength's search guidance improves with the number of iterations and SBP 5 solves all but one instance. Because no constraint initially performs exact weighted counting here and therefore very little domain filtering will occur, not applying support propagation (the BP curves) is not advisable. Observing the foot of the SBP curves, where the number of fails is almost the same, gives an indication of the computational cost of each iteration of BP.

0

instances solved

40

35

30

25

20

15

10

5

0