Fulya Trösser

Simon De Givry

George Katsirelos

Structured Set Variable Domains in Bayesian Network Structure Learning

Keywords: 2012 ACM Subject Classification Computing methodologies → Learning in probabilistic graphical models; Theory of computation → Discrete optimization Keywords and phrases Combinatorial Optimization, Bayesian Networks, Decision Trees Digital Object Identifier 10.4230/LIPIcs.CP.2022.40

Constraint programming is a state of the art technique for learning the structure of Bayesian Networks from data (Bayesian Network Structure Learning -BNSL). However, scalability both for CP and other combinatorial optimization techniques for this problem is limited by the fact that the basic decision variables are set variables with domain sizes that may grow super polynomially with the number of random variables. Usual techniques for handling set variables in CP are not useful, as they lead to poor bounds. In this paper, we propose using decision trees as a data structure for storing sets of sets to represent set variable domains. We show that relatively simple operations are sufficient to implement all propagation and bounding algorithms, and that the use of these data structures improves scalability of a state of the art CP-based solver for BNSL.

Introduction

Bayesian Networks (BNs) are directed probabilistic graphical models, which can describe a normalized joint probability distribution over a potentially large set of random variables, by exploiting conditional independence to decompose the function. Learning the structure of BNs from data (the Bayesian Network Structure Learning problem, BNSL) is a challenging combinatorial optimization problem. There exist constraint-based approaches to learn BNs, which use local conditional independence tests, and score-based approaches, which use a decomposable score function to score each potential structure and aim to find the structure that minimizes this score. The former are known to be efficient, but have trouble with noisy data. The latter yield a known to be NP-hard problem [START_REF] Maxwell | Learning bayesian networks is NP-Complete[END_REF], which additionally has proved very challenging in practice.

There exist complete methods for score-based BNSL based on dynamic programming [START_REF] Silander | A simple approach for finding the globally optimal bayesian network structure[END_REF],

heuristic search [START_REF] Yuan | Learning optimal bayesian networks: A shortest path perspective[END_REF][START_REF] Fan | An improved lower bound for bayesian network structure learning[END_REF], maximum satisfiability [2], branch-and-cut [START_REF] Bartlett | Integer linear programming for the bayesian network 299 structure learning problem[END_REF] and constraint programming [22,[START_REF] Trösser | Improved acyclicity reasoning for bayesian network structure learning with constraint programming[END_REF]. Branch-and-cut and constraint programming have proven to be the most successful of these methods. However, scaling them up remains challenging. One challenge

40:2

Structured Set Variable Domains in Bayesian Network Structure Learning has to do with the decomposition of the scoring functions: these assign a score to each potential set of parents of each vertex and the score of a specific structure is the sum of the scores of each parent set. This means that the objective function must have a term for each potential parent set, a potentially exponential number of terms. There are various methods by which this number is made manageable, but it is still among the greatest obstacles to scalability. Moreover, the best solvers, ILP-based GOBNILP [START_REF] Bartlett | Integer linear programming for the bayesian network 299 structure learning problem[END_REF], and CP-based ELSA [START_REF] Trösser | Improved acyclicity reasoning for bayesian network structure learning with constraint programming[END_REF] also explicitly have this set of parent sets in other parts of the model as well, in the case of ELSA as domains of variables.

Here, we propose exploiting the fact that these domains are structured, i.e., that each value is a set. Specifically, we show that we can represent potential parent sets as paths on decision trees and that using these decision trees we can answer queries more efficiently than by traversing a list of domain values. This feature has not been exploited in BNSL in the past and allows us to solve large instances more efficiently.

Background

Bayesian Networks

A Bayesian Network is a directed graphical model B = G, P where G = V, E is a directed acyclic graph (DAG) called the structure of B and P are its parameters. A BN describes a normalized joint probability distribution. Each vertex of the graph corresponds to a random variable and presence of an edge between two vertices denotes direct conditional dependence. Each vertex v i is also associated with a Conditional Probability Distribution

P (v i | parents(v i ))
. The CPDs are the parameters of B.

Learning a BN from a set of multivariate discrete data using the score based method uses a decomposable scoring function (such as BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF][START_REF] Lam | Using new data to refine a bayesian network[END_REF] or BDeu [START_REF] Buntine | Theory refinement on bayesian networks[END_REF][START_REF] Heckerman | Learning bayesian networks: The combination of knowledge and statistical data[END_REF]) which assigns, based on the data, a score to each potential parent set of each vertex. The BNSL problem is the problem of finding the structure G which minimizes this scoring function.

The number of candidate parent sets can in principle be exponentially large, but it is typically kept in check. For one, the BIC scoring function [START_REF] Schwarz | Estimating the dimension of a model[END_REF][START_REF] Lam | Using new data to refine a bayesian network[END_REF] guarantees that the number of candidate parent sets grows only logarithmically with the size of the data set. Second, there exist dedicated pruning rules [START_REF] Polpo De Campos | Properties of bayesian dirichlet scores to learn bayesian network structures[END_REF][START_REF] Cassio P De Campos | Entropy-based pruning for learning bayesian networks using BIC[END_REF] which reduce the set further. As a last resort, an upper bound can be placed on the cardinality of parent sets. This is necessary especially in larger instances, where it is necessary to limit cardinality to as low as 3 in some cases.

CP-based BNSL

x v,S = 1 ∀v ∈ P (1b) v∈C,S∈P S -C (v) x v,S ≥ 1 ∀C ⊆ P (1c) x v,S ∈ {0, 1} ∀v ∈ P, S ∈ P S(v) (1d)
This is an exponentially large ILP, but on the flip side, the constraints (1c), called cluster constraints are facets of the polytope [START_REF] Cussens | Bayesian network structure learning with integer programming: Polytopes, facets and complexity[END_REF]. Hence, following GOBNILP, ELSA starts with none of the cluster constraints in the linear relaxation and then adds only those that can improve the dual bound. This is an NP-hard problem. GOBNILP solves this NP-hard problem to find violated cluster constraints, while ELSA uses a polynomial time algorithm which can identify a strict subset of all improving cluster constraints. The central element of the algorithm used in ELSA to find cluster constraints uses algorithm 1 on the domains restricted only to values which have reduced cost 0 in the current dual solution of the linear relaxation.

Both in finding improving cluster constraints and in enforcing GAC on the acyclicity constraint, the main bottleneck is line 8 of algorithm 1, which tests whether there exists in D(v) a value which is a subset of a given set. As domain sizes grow drastically faster than the number of random variables, it is crucial to optimize this step. In practical terms, even given the mitigations mentioned earlier, the average domain size can be in the thousands for larger instances.

Related work

A typical approach to dealing with large domain sizes in constraint programming is to enclose the set of domain values with an underestimation and an overestimation and reason with those instead. Sometimes, this can even be achieved without any loss in strength of inference. This is the case, for example, when representing only the bounds of a variables that are only used in linear inequalities. In the case where the values of a domain are sets, the variable is called a set variable. Its domain can be represented with the subset bound scheme [START_REF] Gervet | Conjunto: Constraint logic programming with finite set domains[END_REF], which underestimates by a set indicating all elements which appear in all remaining domain values and overestimates by a set indicating all elements which appear in any remaining domain value. The length-lex scheme uses lexicographic and cardinality information to get a tighter under-and over-estimation [START_REF] Gervet | Length-lex ordering for set CSPs[END_REF]. However, detecting infeasibility of the acyclicity constraint is crucial for the performance of CPBayes and even more for ELSA. Hence, over-estimating the actual domain in our case would lead to poor performance.

C
Hawkins et al. [START_REF] Hawkins | Solving set constraint satisfaction problems using ROBDDs[END_REF] followed an approach which is closer to our own, by using ROBDDs (reduced ordered binary decision diagrams) to represent domains. ROBDDs are diagrams like decision trees, but they require the same variable ordering in each branch and isomorphic subgraphs are merged, so that the underlying graph is a DAG rather than a tree. They can be significantly more compact than decision trees. However, Hawkins et al. used them in a setting where all constraints can be expressed as operations on ROBDDs. They do not deal with costs of the domain values, and in particular with reduced cost filtering.

Decision Trees as domain store

The set of sets that are in a domain can be seen as the set of solutions of a propositional formula, in which we have a propositional variable for each element of the universe. Therefore, knowledge compilation languages such as ROBDDs can be used to represent a domain.

There exist several queries and operations performed on the domains in ELSA, but not all are critical to optimize, as they are not performed often enough to dominate the runtime.

In particular, we want to address the test in line 8, which asks whether the domain contains a set which is a subset of another given set. Therefore, the main queries that need to be supported efficiently by a domain store for our purposes are:

1. Does there exist a domain value S such that S ⊆ T for some T ?

2. Does there exist a domain value S with reduced cost 0 such that S ⊆ T for some T ?

And the main operations, which also have to support backtracking, are:

1. Pruning a single value S

Updating the reduced cost of a value

The main issue that disqualifies ROBDDs and other reduced representations for us is that operation 1, reduced cost filtering, may remove arbitrary values, shattering the shared suffixes that an ROBDD exploits, which means that pruning may result in increasing the size of the representation and is not even guaranteed to be in linear time. Instead, we use decision trees here, in particular binary decision trees with implied literals, inspired by a similar technique in BDDs [START_REF] Lai | Reduced ordered binary decision diagram with implied literals: A new knowledge compilation approach[END_REF]. The main use of decision trees is in machine learning for classification, but their use as a data structure for representing sets of sets (or, equivalently, a knowledge compilation language) is straightforward.

We give below definitions for the specific case of binary decision trees and binary classification, as that is all we need.

Definition 1 (Binary decision tree).

Let A be a set of features {a 1 , . . . , a n } with Boolean along that path, it holds that I(l(n)) = l(e). We say that T classifies I as l(n) and that I and the path from the root to n are consistent with each other, or simply that I and n are consistent with each other.

To see how we can use binary decision trees as a data structure for a set of sets, observe that we can set the features to be the variables of the indicator function of the sets in the domain and the classes as in-set and not-in-set.

This allows us to further optimize the representation. Since we only care about the in-set class, from now on we assume that all nodes and arcs that do not appear on a path from the root to a leaf n with l(n) = in-set are removed from the decision tree.

Additionally, we can eliminate some nodes by adding implied literals in each node of the tree.

Definition 2 (Binary decision tree with implied literals). A binary decision tree with implied literals is a decision tree in which each node n (internal or leaf) is additionally labeled with a

set of literals lit(n, C i ) ⊆ {a = v | a ∈ A, v ∈ {true, f alse}} for i ∈ {1, 2}. An instantiation
I is consistent with a path to a leaf t with l(t) = C i if it is consistent with all the arcs it follows and all implied literal labels lit(n, C i ) for each node n on the path from the root to t.

In our case, we abbreviate lit(n, in-set) to lit(n), as we ignore the class not-in-set.

Decision trees with implied literals allow us to collapse chains, i.e., paths along which every node has outdegree 1, into a single node. Hence, they are not more compact than those without implied literals by more than a linear factor, but they have almost no overhead and, in preliminary experiments, we found them to provide some performance improvement.

In machine learning, the objective is not only to construct models that perform well on the training set, but that also generalize. Hence, it is not only acceptable, but also desirable to misclassify some samples in training sets, if that means a smaller and hence more general decision tree. In our setting, however, where we use decision trees to model a Boolean function, we accept no error. So no two sets that belong to different classes, i.e., one in in-set and one in not-in-set, are allowed to both be consistent with the same leaf node.

We place an additional constraint on the decision trees we construct, which is that each leaf node must be consistent with exactly one positive instantiation. This ensures that there exists a bijection between leaves of the tree and values in the domain. This is not as significant a constraint as it might seem at first. A leaf node n that is consistent only with positive instantiations but more than one of them is expanded into a full binary tree of depth k, where k is the number of variables (features) which have not appeared on the path from the root to n. However, for the queries that we care about, this means only that the corresponding algorithm will have to traverse an additional k nodes before answering, and, crucially, will only arrive at this point when it is guaranteed that it will give a positive answer. Even that overhead can be eliminated with some care. Indeed, while traversing the decision tree, we can determine that we have reached such a node n if the number of possible instantiations that are consistent with n is equal to the number of leaves reachable from n. The former is 2 n-lvl , where lvl is the distance from the root to n. The latter can be computed on construction and updated as values are removed. If these are equal, we know that the subtree contains all possible subsets and we can answer our query without more search. We give more detail later.

Constructing decision trees.

Constructing a minimum decision tree is NP-hard with respect to several metrics [START_REF] Laurent | Constructing optimal binary decision trees is NP-complete[END_REF]. We use the information gain heuristic [START_REF] Quinlan | Induction of decision trees[END_REF] to choose which variable to branch on in each node.

C P 2 0 2 2 40:6

Structured Set Variable Domains in Bayesian Network Structure Learning

It is a natural side effect of computing the information gain that we learn how many of the sets that are consistent with a node n contain the literals a = true and a = f alse for all a ∈ A. If either of these is 0, then its negation is added to the implied literal label for n and a is not considered as a candidate for branching. We also experimented with optimizing the in-memory layout for better cache behavior. Compared to the van-Emde Boas layout [START_REF] Van | Preserving order in a forest in less than logarithmic time[END_REF], a depth-first, false-child first layout performed better.

Maintaining a decision tree during search

It is fairly straightforward to update a decision tree for a pruning. In order to prune a value, we remove the unique leaf node that corresponds to it. Once we remove a leaf, its parent may no longer be able to reach any more leaves, hence we propagate this removal upwards.

We associate each removed node with the decision level in which it was removed, so on backtracking we add them back to restore the tree to its correct state.

This guarantees that the tree representation of a domain only shrinks down a branch of the branch and bound tree. Hence, the tree can remain static and we only mask nodes that do not lead to any leaves that correspond to unpruned values, which is simple to implement.

Reduced costs

ELSA solves the linear relaxation ( 1) from scratch at every node, and then strengthens it by discovering new violated cluster inequalities using the acyclicity checker (algorithm 1). Both these algorithms require an efficient implementation of the subset query on the subset of values which have reduced cost 0. In contrast to the domain itself, however, this set is reset to the empty set at the beginning of every node and grows monotonically until it admits an acyclic solution. Here again, the fact that there exists a bijection between values and leaves of the tree allows us to represent the set of 0-cost values as a subset of the full decision tree.

Every time the reduced cost of a value reaches 0, the unique leaf it corresponds to, as well as all its parents, are added to the set of visible nodes for these queries. This is implemented as an additional mask on top of that which hides pruned values.

Subset queries

To answer the query "does the domain contain a value S such that S ⊆ T ?", we perform a depth first traversal of the tree. At each node n, we check l(n). If l(n) / ∈ T , we only allow DFS to follow the outgoing arc labeled with false. If l(n) ∈ T , we allow DFS to follow both outgoing arcs. If the label lit(n) contains a literal p / ∈ T , we backtrack. If we reach a leaf, we stop and report success. If we exhaust the search without reaching a leaf, we report failure.

When this procedure reaches a node which is the root of a complete subtree of depth k, with no additional implied literal labels, it is guaranteed to terminate after visiting exactly k nodes and report success. Indeed, since this is a complete subtree, one of the outgoing arcs is always available to the depth first search, and it will reach a leaf after k more steps.

This procedure can be used to answer subset queries either on the entire domain, masking away only pruned values, or on those values which have reduced cost 0, masking away both pruned values and those whose reduced cost is greater than 0.

Experimental Results

We implemented decision trees as the domain representation on top of ELSA. The default implementation of a subset query in ELSA iterates over all domain values and returns if F. Trösser, S. de Givry, and G. Katsirelos 40:7

it finds one that is a subset of T . We replaced this by the depth-first traversal described in section 4 and denote this solver 1 ELSA IG . We compare against the previous version of ELSA 2 , GOBNILP 3 , and CPBayes 4 .

The datasets come from the UCI Machine Learning Repository 5 , the Bayesian Network Repository 6 , and the Bayesian Network Learning and Inference Package 7 . We have 51 medium datasets with |V | < 64, and 18 large datasets with 64 ≤ |V | < 128.

Local scores were computed from the datasets using B. Malone's code 8 . BDeu and BIC scores were used for medium datasets (less than 64 variables) and only BIC score for large datasets (above 64 variables). The maximum number of parents was limited to 5 for large datasets (except for accidents.test with maximum of 8), a high value that allows learning even complex structures [START_REF] Scanagatta | Learning bayesian networks with thousands of variables[END_REF]. For example, jester.test has 100 random variables, a sample size of 4, 116 and 770, 950 parent set values. For medium datasets, no restriction was applied except for some BDeu scores, where we limit sets to 6 or 8 to complete the computation of the local scores within 24 hours of CPU-time [START_REF] Lee | An experimental analysis of anytime algorithms for bayesian network structure learning[END_REF].

For the experiments, we modified the C++ source of CPBayes v1.1 just to allow us to run it with datasets having more than 64 variables. All computations were performed on a single core of Intel(R) Xeon(R) Gold 6248R CPU @ 3.00GHz and 1 TB of RAM with a 1-hour CPU time limit for the 51 medium datasets, as well as 3 of the large datasets: In Table 1, we show the time needed to find the optimal solution and prove optimality for all these solvers. We see that, while the use of decision trees has little effect, either positive or negative, for the smaller instances, it makes a great difference in the larger instances.

In particular, ELSA IG is the only solver that can prove optimality for the baudio.test 1 Comparison of GOBNILP, CPBayes, ELSA, and ELSA IG in terms of total running (and search) time in seconds. Time limit for the datasets above the line is 1 hour, and for the rest it is 10 hours. Datasets are sorted by increasing total domain size for each time limit category. For CPBayes as well as all variants of ELSA we report in parentheses time spent in search, after preprocessing finishes. † indicates a timeout.

and
where neither ELSA nor ELSA IG can match it. It seems, however, that ELSA IG is overall 285 the best performer. issues with ROBDDs or even DNNFs to achieve even more significant speedups.

P 2 0 2 2 40: 4 Structured

 24 Set Variable Domains in Bayesian Network Structure Learning

F

  domains and C 1 , C 2 be two classes. A binary decision tree T over the features A is a directed rooted binary tree. Each internal node n of T is labeled with a feature l(n) ∈ A and each arc e (of the at most two outgoing arcs) from n is labeled with l(e) ∈ {true, f alse} and are mutually exclusive. Each leaf node t is labeled with l(t) ∈ {C 1 , C 2 }. Given an instantiation I of the features, there is a unique path from the root to a leaf t so that for each arc e = (n, c)

  kdd.ts, kdd.test, and kdd.valid. For the remaining 15 large datasets, we had a 10-hour CPU time limit. For the preprocessing phase, we used two different settings depending on problem size n = |V |: l min = 20, l max = 26, r min = 50, r max = 500 if n ≤ 64, else l min = 20, l max = 20, r min = 15, r max = 30, where l min , l max are partition lower bound sizes and r min , r max are the number of restarts for the local search.

287

  We have shown that, in the BNSL problem, we can exploit the structure of domains to get 288 a significant speedup in learning the structure of BNs of larger datasets. Specifically, we 289 have shown that by treating domains as sets of sets instead of sets of values, and using 290 decision trees to represent these sets, we can answer subset queries significantly faster. This 291 is unlike the typical approach to handling large domains in CP, which uses over-and under-292 approximations. Although the current implementation shows some significant improvements, 293 answering subset queries is still the most time consuming operation performed by the solver. 294 Moreover, the fact remains that decision trees as a knowledge compilation language are fairly 295 weak in terms of conciseness. It remains an open question whether we can overcome the 296

F. Trösser, S. de Givry, and G. Katsirelos 40:3 Algorithm 1 Acyclicity checker

	1 acycChecker (P, D)
	2 order ← {}
	3 changes ← true
	4 while changes do
	5	changes ← f alse
	6	foreach v ∈ P \ order do
	8 8	if ∃S ∈ D(v) s.t. S ⊆ order then
	10 10	order ← order + v
	11	changes ← true
	12 return order
	of the ILP (1), which was proposed by Bartlett and Cussens [1] for the GOBNILP solver.
		min	σ v (S)x v,S	(1a)
		v∈P,S⊆V \{v}
		s.t.
		sets of X. These are unsurprisingly called parent set variables. There exists an acyclicity
		constraint over these which requires that their instantiation yields an acyclic graph. ELSA
		enforces GAC on this constraint. The central part of the GAC algorithm is algorithm 1,
		acycChecker. acycChecker determines in time O(n 2 d) whether the current set of domains
		admits an acyclic solution, based on the property that in any acyclic graph, for any subset
		of vertices C, at least one of the vertices v ∈ C has a parent set that does not intersect C.
		In addition, ELSA computes lower bounds by approximately solving the linear relaxation

ELSA

[START_REF] Trösser | Improved acyclicity reasoning for bayesian network structure learning with constraint programming[END_REF]

is a CP-based solver for the BNSL, based on the CPBayes solver

[22]

. The constraint model used in ELSA has several features that we do not discuss here. Instead, we focus on the part that is relevant to our contribution. For each random variable X, there exists a corresponding CSP variable P X whose domain is the set of candidate parent S∈P S

(v)

0 2 2 40:8 Structured Set Variable Domains in Bayesian Network Structure Learning

  jester.valid datasets. For the only instances where ELSA is significantly worse than CPBayes, bnetflix.ts, bnetflix.test, and bnetflix.valid, ELSA IG either closes the gap back down (bnetflix.valid) or is faster yet than CPBayes (bnetflix.ts and bnetflix.test). However, ELSA IG regresses with respect to ELSA in the accidents dataset and in plants.test. Part of the reason for this is that the benefit of the decision trees in terms of the reduction of the cost in answering the subset queries is comparatively reduced, hence the other overheads of decision trees dominate. For example, in bnetflix.ts, where ELSA IG significantly outperforms ELSA, ELSA looks at an average of 3315 values to answer each subset test, while ELSA IG visits just 90 nodes of the decision tree. On the other hand, in accidents.test, ELSA looks at an average of 80 values to answer each subset test, while ELSA IG visits 20 nodes of the decision tree. This difference is not enough to overcome other overheads.With respect to GOBNILP, ELSA IG mostly outperforms it, but there are some instances

		n	sum |D| GOBNILP	CPBayes	ELSA	ELSA IG
	carpo100_BIC	60	423	0.5	76.7 (27.5)	52.6 (0.1)	52.5 (0.0)
	insurance1000_BIC 27	506	0.6	31.6 (0.0)	32.8 (0.0)	37.2 (0.0)
	spectf_BIC	45	610	1.4	4.2 (3.5)	0.8 (0.0)	1.0 (0.1)
	sponge_BIC	45	618	1.6	5.1 (3.3)	1.8 (0.0)	2.1 (0.0)
	insurance1000_BDe 27	792	0.6	34.8 (0.0)	34.3 (0.0)	39.2 (0.0)
	alarm1000_BIC	37	1002	1.3	191.1 (159.1)	34.4 (1.0)	37.9 (1.9)
	flag_BDe	29	1324	4.0	16.6 (15.6)	1.0 (0.2)	1.3 (0.2)
	autos_BIC	26	2391	11.9	18.4 (0.0)	19.2 (0.0)	19.9 (0.1)
	soybean_BIC	36	5926	48.9	51.9 (1.7)	50.8 (3)	49.6 (0.0)
	wdbc_BIC	31	14613	86.3	459.4 (398.0)	56.0 (2.4)	61.7 (1.7)
	autos_BDe	26	25238	1005.2	239.5 (0.1)	145.8 (0.8)	177.1 (0.3)
	kdd.ts	64	43584	508.8	†	1452.3 (274.6)	1355.2 (141.3)
	steel_BIC	28	93026	† 1265.6 (1196.1)	124.2 (71.8)	100.6 (45.7)
	kdd.test	64	152873	3178.0	†	1594.3 (224.4)	1519.6 (48.9)
	mushroom_BDe	23	438185	†	167.0 (4.9)	182.6 (58.9)	150.1 (16.7)
	bnetflix.ts	100 446406	†	1086.9 (876.3)	2103.1 (1900.9)	557.9 (358.4)
	plants.test	111 520148	†	† 28049.6 (26312.9)	35961.7 (33712.7)
	jester.ts	100 531961	†	†	21550.5 (21003.7)	7951.4 (7301.6)
	accidents.ts	100 568160	1932.2	†	2302.2 (930.0)	†
	plants.valid	111 684141	†	† 17801.6 (14080.2)	19819.2 (14547.9)
	jester.test	100 770950	†	†	30186.8 (29455.0)	9644.5 (8742.8)
	baudio.test	100 1016403	†	†	† 31077.1 (29028.1)
	bnetflix.test	100 1103968	† 5794.5 (5486.2)	10333.1 (10096.5)	1448.8 (1137.7)
	bnetflix.valid	111 1325818	† 998.1 (451.0)	10871.7 (10527.7)	1476.5 (1041.5)
	accidents.test	100 1425966	14453.1	†	3641.7 (680.7)	8434.1 (4723.0)
	jester.valid	100 1463335	†	†	† 31949.5 (30624.2)
	accidents.valid	100 1617862	27730.5	†	†	†
	1 Available at https://gkatsi.github.io/elsa-cp22.tar.gz
	2 Available at https://gkatsi.github.io/elsa-ijcai21.tar.gz
	3 Version 1.6.3 with CPLEX 12.7.1
	4 Retrieved from http://cs.uwaterloo.ca/~vanbeek/Publications/CPBayes.zip
	5 http://archive.ics.uci.edu/ml

6

http://www.bnlearn.com/bnrepository 7 https://ipg.idsia.ch/software.php?id=132 8 http://urlearning.org C P 2

Table

Funding This work has been partly funded by the "Agence nationale de la Recherche" (ANR-16-CE40-0028 Demograph project and ANR-19-PIA3-0004 ANTI-DIL chair of Thomas Schiex).

Acknowledgements Thanks to the GenoToul (Toulouse, France) Bioinformatics platform for computational support. George Katsirelos; licensed under

: https://gkatsi.github.io/elsa-cp22.tar.gz


Abdelkader Beldjilali

Pierre Montalbano

David Allouche

George Katsirelos

Simon De Givry

Parallel Hybrid Best-First Search

Keywords: 2012 ACM Subject Classification Computing methodologies → Parallel algorithms Combinatorial Optimization, Parallel Branch-and-Bound, CFN Digital Object Identifier 10.4230/LIPIcs.CP.2022.36

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Cost Function Networks (CFNs), also known as Weighted Constraint Satisfaction Problems (WCSPs) [START_REF] Meseguer | Soft constraints processing[END_REF] is a mathematical framework which has been derived from Constraint Satisfaction Problems by replacing constraints with cost functions. In a CFN, we are given a set of variables with an associated finite domain and a set of local cost functions. Each cost function involves some variables and associates a non-negative integer cost to each of the possible combinations of values they may take. The usual WCSP problem considered is to assign all variables in a way that minimizes the sum of all costs. This minimization problem is NP-hard, and exact methods usually rely on Branch and Bound (B&B) algorithms exploring a binary search tree with soft local consistency maintained at each node in order 36:2

Parallel Hybrid Best-First Search to improve the problem lower bound (represented by c ∅ ) and prune domain values with a forbidden cost (represented by a maximum cost k) [START_REF] Cooper | Soft arc consistency revisited[END_REF].

Contraint Programming (CP) exact approaches usually rely on Depth-First Search (DFS) methods while Integer Linear Programming (ILP) approaches explore the tree in a best-first manner by exploiting strong bounds. We are interested in hybrid methods combining depthfirst and best-first with possibly weaker bounds but faster to compute. This is the case of the Hybrid Best-First Search (HBFS) method [2]. HBFS is a B&B algorithm for solving WCSPs.

Dealing with parallel computers or grids to speed-up solving time of exact methods has been explored in many different ways. For grids, with slow network interconnection, MapReduce is a general approach exploiting problem decomposition into independent subproblems solved in parallel (map on the grid processors) and then sequentially reduced at the end of the resolution. In CP, this decomposition approach is called Embarrassingly Parallel Search (EPS) [START_REF] Malapert | Embarrassingly parallel search in constraint programming[END_REF]. MapReduce has been applied also in the context of non serial dynamic programming in Graphical Models [START_REF] Otten | And/or branch-and-bound on a computational grid[END_REF] and CFNs [START_REF] Allouche | Towards parallel non serial dynamic programming for solving hard weighted csp[END_REF]. Message-passing approaches, on the other hand, take advantage of the low-latency communication of supercomputers, consisting of a large number of multiprocessor servers interconnected at high speed and low latency.

This allows for finer granularity in B&B parallelization. According to a recent survey [START_REF] Gent | A review of literature on parallel constraint solving[END_REF],

parallelizing the search based on message-passing and parallel B&B in CP are difficult problems and still poorly explored. In CP, for example, COMET [START_REF] Michel | Parallelizing constraint programs transparently[END_REF] uses work-stealing where workers which have run out of work take unexpanded nodes from other workers, leaving them less work to do and keeping all workers busy. In ILP, a recent review on parallel B&B was proposed in [START_REF] Ralphs | Parallel solvers for mixed integer linear optimization[END_REF]. We selected the Master-Worker protocol as the basis for our approach.

Other approaches rely on portfolios.

In this work we describe a parallel version of HBFS. We give an empirical evaluation on combinatorial optimization academic problems from Operations Research and real-life Graphical Model problems occuring in genetics and biology. Our experimental study analyses solving time and speed-ups of the parallel version compared to the original sequential HBFS.

We also compare our approach with a parallel ILP solver (IBM Ilog cplex). Moreover, we performed experiments on a high-performance cluster to study the scalability of our algorithm and compare with EPS.

Hybrid Best-First Search

The sequential version of HBFS [2] is a B&B method for CFNs that combines Best-First Search (BFS) and Depth-First Search (DFS). Like BFS, HBFS provides an anytime global lower bound on the optimum, while also providing anytime upper bounds, like DFS. Hence, it provides feedback on the progress of search and solution quality in the form of an optimality gap. Besides, it exhibits highly dynamic behavior that allows it to perform on par with methods like Limited Discrepancy Search [START_REF] Harvey | Limited discrepancy search[END_REF] and frequent restarting [START_REF] Gomes | Boosting combinatorial search through randomization[END_REF][START_REF] Demirovic | Solution-based phase saving for CP: A value-selection heuristic to simulate local search behavior in complete solvers[END_REF] When the problem is solved, the Master kills all the workers and returns the optimum value.

According to a round robin schema, the Master sends open nodes to every idle worker in a balanced way, ensuring a natural load balancing between the workers as soon as the number of open nodes in the frontier is larger than the number of workers. Moreover, an initial backtrack limit of Z i = 1 associated to each Worker i favors the production of open nodes at the beginning of the search. Each Z i is bounded by N as in sequential HBFS so that no worker takes too long.

The pseudo-code of the Master (resp. Worker) is given in Algorithm 2 (resp. Alg. 3).

In the implementation, we avoid to send the same solution twice to a Worker. Moreover, workers send their solution only if it improves compared to the last solution sent by the Master. This strategy allows to shorten messages in the Master-Worker protocol.

Improving the ramp-up phase

We observed that at the beginning of the search the first active worker may take a long time to build its list of open nodes when it reaches the initial backtrack limit (equal to one of open nodes when HBFS stops backtracking. During this period, called the ramp-up phase (where some workers have not been assigned at least one task), no parallelism is exploited.

Experimental Results

We implemented in C++ our parallel HBFS in the CFN solver toulbar2.2 We used the boost MPI library for the Master-Worker communication protocol. We kept default parameters of toulbar2 except no dichotomic branching in order to explore a binary search tree with DFS (option -d:). The variable ordering heuristic is dom/wdeg [START_REF] Boussemart | Boosting systematic search by weighting constraints[END_REF] combined with last conflict [START_REF] Lecoutre | Reasoning from last conflict(s) in constraint programming[END_REF]. The value ordering heuristic exploits the last solution found if any [START_REF] Demirovic | Solution-based phase saving for CP: A value-selection heuristic to simulate local search behavior in complete solvers[END_REF] or else EDAC existential value [START_REF] De Givry | Existential arc consistency: Getting closer to full arc consistency in weighted csps[END_REF]. EDAC is also used as soft local consistency during search.

Instances were preprocessed by VAC [START_REF] Cooper | Soft arc consistency revisited[END_REF] and the resulting CFNs saved to files before the experiments to reduce the setup sequential time of paralllel HBFS. We compared both the sequential and parallel version of HBFS and also with the integer programming solver cplex (version 20.1 with non-premature stop parameters EPAGAP=EPGAP=EPINT=0). We set the number of threads used by cplex to the desired number of cores.

Experiments were performed either on medium-scale computers (24-core Intel Xeon E5-2687W v4 at 3 GHz and 256 GB) with 1-hour timeout or on a large-scale cluster with more than 10, 000 cores (36-core per node of Intel Skylake 6140 at 2. We tested the methods on four benchmarks selected from [START_REF] Hurley | Multi-Language Evaluation of Exact Solvers in Graphical Model Discrete Optimization[END_REF] with a total of 134 instances: two academic benchmarks taken in Operations Research, uncapacitated warehouse location problem (Warehouses) with 15 instances [START_REF] Kratica | Solving the simple plant location problem by genetic alg[END_REF] and DIMACS maximum clique problem with 62 instances (MaxClique) 3 and two real-life Graphical Model benchmarks, linkage analysis problem occuring in genetics (Linkage) with 22 instances coming from UAI Evaluation 2008 4and computational protein design problem in biology (CPD) with 35 instances [START_REF] Allouche | Computational protein design as an optimization problem[END_REF]. We applied the tuple encoding to convert Linkage and CPD to integer linear programs [START_REF] Hurley | Multi-Language Evaluation of Exact Solvers in Graphical Model Discrete Optimization[END_REF]. For a comparison on MaxClique with another parallel branch and bound implementation, see [START_REF] Mccreesh | The shape of the search tree for the maximum clique problem and the implications for parallel branch and bound[END_REF]. Figure 2 Comparison on a medium-scale computer between sequential versus parallel HBFS with or without burst mode. The x-axis represents normalized time (with 0.2 corresponding to 720 seconds). The y-axis corresponds to normalized lower and upper bounds on 134 instances (with 1 corresponding to the optimum or best known cost, see the text description).

Comparison of parallel HBFS with its sequential version

We compared the anytime behavior of sequential (HBFS-1) and parallel HBFS (with 10 or 20 cores) with or without burst mode (see Sec. 3.1) on a medium-scale computer. We summarize the evolution of lower (clb in Alg. 1 and 2) and upper bounds (cub) for each method over all instances in Fig. 2. Specifically, for each instance we normalize all costs as follows: the initial lower bound c ∅ produced by EDAC is 0; the best but potentially suboptimal solution found by any method is 1; the worst solution is 2. This normalization is invariant to translation and scaling. Additionally, we simply normalize time from 0 to 1, corresponding to 1 hour. A point x, y on the lower bound line for method M in Fig. 2 means that after normalized runtime x, method M has proved on average over all instances a normalized lower bound of y and similarly for the upper bound. First, we observed that all parallel versions significantly outperformed the sequential HBFS lower bound curve. Concerning upper bound curves, the burst mode gave a clear advantage to parallel HBFS especially at the beginning of the search. In the sequel of the paper, we always report results of parallel HBFS with burst mode. As shown in the figure, increasing the number of cores from 10 to 20 slightly improved the bounds.

In Table 1 we report the number of instances solved by sequential and parallel HBFS for each benchmark. Parallel HBFS solved 1 more instance than the single core version in Linkage and 1 (resp. 3) in MaxClique using 10 (resp. 20) cores. We made local comparisons of solving times (shown in parentheses) by averaging on the subset of instances solved by the three methods (HBFS-1, HBFS-10, HBFS-20). It allows us to display overall speed-up of parallel approaches by giving the ratio of total sequential over parallel time. Parallel HBFS obtained near linear speed-up on MaxClique. Recall that 1 core is used by the master and the rest by the workers in the Master-Worker approach preventing us from full linear speed-up.

On CPD and Linkage the speed-up was halved. For Warehouses, only 50% of reduction in overall time was observed. This can be explained partly by the limited number of search nodes (Table 4 in Supplementary Material).We also observed that the evaluation of right branches made by the first active worker starting from the root node took most of the time. This is due to the fact that a first solution has been found by the worker resulting in more propagation on the right branches especially near the root. This pathological phenomenon did not appear on the other benchmarks.

78.26

Table 1 Number of solved instances within 1 hour (except for sequential HBFS-1 run on the cluster with a larger timeout of 10 hours) and average time in seconds in parentheses. To compute the mean we only consider for a given method (toulbar2 HBFS or cplex) the instances solved with any number of cores on the same computer (server with 3 GHz cores or cluster with 2.3 GHz cores).

Comparison of parallel HBFS with integer programming

In Table 1 we also report the number of instances solved and their average solving time (as explained above) by cplex using multithreading. It clearly dominates HBFS on Linkage (Supp. Fig. 5).For Warehouses, the differences are less important still in favor of cplex.

For MaxClique, although the global picture shows that it solved six more instances than HBFS with 20 cores, both methods performed well on different subsets of instances (e.g., HBFS-20 solved two instances -brock400_4 and sanr400_0.7 -unsolved by cplex-20 whereas cplex-20 solved eight instances unsolved by HBFS-20). For CPD, the CFN approach largely dominates the integer programming approach for all the instances. Concerning anytime curves shown in Fig. 3 (see also Supp. Fig. 4 and5), the CFN approach is also significantly superior to cplex on average in producing good upper bounds faster, HBFS-20 being the best method. Concerning overall speed-up, cplex had difficulties to benefit from parallelism on CPD, Linkage, and Warehouses where it usually develops a small amount of search nodes (less than 7, 059 nodes except on Linkage/pedigree19 and pedigree40), resulting in poor speed-up except in a few cases. The speed-up is better on MaxClique but seems to stagnate when going from 10 to 20 cores (it was even slower on four instances).

Comparison of parallel HBFS with EPS on a cluster

The EPS approach is a two-phase procedure. First, the problem to be solved is decomposed into a list of l independent subproblems. Next, all the subproblems are solved in parallel (with at most p workers running at the same time) based on a particular scheduling strategy with no communication between the workers. For optimization problems, we need to provide a good initial upper bound. Otherwise the search tree can be much larger than needed. In the first phase, we used the original HBFS method to collect l subproblems. As soon as Each open node ν defines an independent subproblem with partial assignment ν.δ. In order to collect open nodes more rapidly we fix the (maximum) backtrack limit Z = N = 1. Ideally l should be 30 × p with p the number of available cores [START_REF] Malapert | Embarrassingly parallel search in constraint programming[END_REF]. In the second phase, we schedule on the cluster the subproblems that are solved by the original HBFS method using a simple scheduling heuristic based on increasing |ν.δ|.

In Table 2 we report for nine difficult instances their optimum value, the upper bound found at the end of EPS Phase-1, the actual number of generated subproblems, the average solving time of all subproblems, the maximum solving time, the number of failed subproblems (timeout of 1 hour) and the overall solving time of EPS Phase-2 using 180 cores on the cluster.

We compare with HBFS using the same number of cores. Our EPS strategy failed on 4/9 instances. In parentheses, we indicate the maximum depth ν.depth of failed subproblems.

Clearly, finding the right number l of not-too-difficult subproblems corresponding to partial assigments greater than a given depth is a challenging task. In our experiments, we tried with different values for l ∈ [50, 6000], selecting the largest threshold value with a Phase-1 duration being less than 1 second for Linkage (l = 6000), 6 seconds for MaxClique (l = 6000) and 44 seconds for CPD (l = 1000). On the opposite, we did not tune any specific parameter for our parallel HBFS method.

In Table 1 we also report the overall speed-up of HBFS-180 compared to HBFS-1 on the cluster. HBFS-180 got a two-order-of-magnitude speed-up on Linkage.

Conclusion

Although the speed-up offered by the parallel version of HBFS was very instance dependent, we observed significant gain on several instances, outperforming in some cases state-of-theart solvers like cplex. Even if the scalability of our approach must be subject of deeper investigation, due to the minimal size of the information shared between the Master and the Workers, our approach is very likely compliant with a larger number of cores.

A more challenging task which remains as future work is to exploit the structure of CFNs by parallelizing Backtrack with Tree Decomposition (BTD-HBFS) [2]. Shared memory protocols may be more suitable for this task to make learnt nogoods available to all Workers.

On the practical side, our parallel HBFS could ran in conjunction with a parallel large neighborhood search strategy [START_REF] Ouali | Iterative Decomposition Guided Variable Neighborhood Search for Graphical Model Energy Minimization[END_REF] offering even better anytime lower and upper bounds. Comparison on a medium-scale computer between toulbar2 using parallel HBFS (with burst mode) and cplex using multiple threads. The x-axis represents normalized time (with 0.2 corresponding to 720 seconds). The y-axis corresponds to normalized lower and upper bounds on 134 instances (with 1 corresponding to the optimum or best known cost, see the text description

Function 3 4 NodesRecompute

 34 HBFS(clb,cub): integer ; /* Returns the optimum value */ 1 open := {ν(δ = ∅, lb = clb)} ; /* Initializes the open list with a root node */ 2 while (open ̸ = ∅ and clb < cub) do ν :=pop(open) ; /* Chooses a node with minimum lower bound and maximum depth */ Restores state ν.δ, leading to assignment Aν , maintaining soft local consistency ; := NodesRecompute + ν.depth ; 5 cub :=DFS(Aν ,cub,Z) ; /* Increase Nodes and put all right open branches in open */ 6 clb := max(clb, lb(open)) ; if (NodesRecompute > 0) then 7 if (NodesRecompute/Nodes > β and Z ≤ N ) then Z := 2 × Z; 8 else if (NodesRecompute/Nodes < α and Z ≥ 2) then Z := Z/2; return cub; Algorithm 1 Hybrid Best-First Search. Initial call: HBFS(c ∅ ,k) with Z = 1. DFS increases a counter Nodes at each branching decision. It can backtrack (taking right branches) up to a limit of Z backtracks. When this limit is reached, all the unexplored right branches are placed in open. HBFS controls the balance between best-first search (partially exploring more open nodes) and depth-first search (complete exploration from a given starting node). Best-first search requires recomputing the state ν.δ of a node (line 3) which can be costly in practice. HBFS uses a simple rule to limit this recomputation effort (measured by NodesRecompute at line 4). It tries to keep the ratio NodesRecompute Nodes in the interval [α, β] by increasing (by a power of two) the backtrack limit Z if the ratio value is above β or decreasing Z if it is below alpha (lines 7-8). Initially, Z is set to 1. In order to avoid exponential DFS behavior, HBFS limits the maximum value taken by Z to N . We kept the same value α = 5%, β = 10%, N = 2 14 in our experiments as in the original paper [2]. C P 2 0 2 2 Parallel HBFS The parallel version of HBFS is based on the Master-Worker parallel paradigm [21] where the Master is in charge of the open node frontier and dispatches the current best (with minimum lower bound) open node plus the current best solution found so far to the next available Worker. The Worker performs a bounded DFS starting from the received node and returns to the Master the resulting list of open nodes (see Fig. 1, with a DFS limit here of 4 backtracks). The Worker also returns the best solution found during its restricted search if any. Only the Master has a global view of the whole search and reports optimality gaps ( cub-clb cub ) until the proof of optimality is reached: when the current best lower bound in the frontier of open nodes, including active worker starting nodes, is equal or greater than the cost of the best solution found so far or the frontier is empty and there are no active workers.

Function 11 A

 11 HBFS-Master(clb, cub, S): integer ; /* S queue of workers, return the optimum */ open := {ν(δ = ∅, lb = clb)} ; /* Initializes the open list with a root node */ I := S ; /* Queue of idle workers */ A := ∅ ; /* Maps active workers to open nodes currently being processed */ while ((open ̸ = ∅ or A ̸ = ∅) and clb < cub) do while (open ̸ = ∅ and I ̸ = ∅) do ν :=pop(open) ; /* Chooses a node with minimum lower bound and maximum depth */ i :=popFront(I) ; /* Unqueue the first idle worker */ A := A ∪ {(i, ν)} ; Send ν and best solution cub to Worker i ; 9 Receive a list of open nodes V and solution cub ′ by worker j ; /* Wait for message */ push(open, V) ; /* Adds worker open nodes to the Master open list */ cub := min(cub, cub ′ ) ; /* Checks if a better solution as been found */ 10 pushBack(I, j) ; /* Pushes Worker j at the end of the idle worker queue I */ := A \ {(j, A[j])} ; /* Removes Worker j from active workers */ clb := max(clb, min(lb(open), min{lb(ν) for (i, ν) ∈ A})) ; /* Global lower bound */ return cub; Algorithm 2 Parallel HBFS-Master. Initial call for p workers: HBFS-Master(c ∅ ,k,(1, . . . , p)).

  We modified our communication protocol to send a message to the master as soon as an open-node has been collected or a new solution has been found by a worker inside its DFS subroutine (line 12). Such messages are received by the Master (line 9) which does not change the Worker state to idle (lines 10 and 11) until it receives a closing-node message by the Worker (sent at line 15). By doing so, it allows the Master to distribute open nodes to idle workers earlier before the first active worker has finished its initial DFS. We call this modified Master-Worker protocol the burst mode. However, the Worker can potentially send O(nd) more messages and it disallows data compression of the open list messages. 1

  HBFS has more than l open nodes in its frontier it stops and returns the current upper bound (cub) and the list of open nodes (without those having a lower bound lb(ν) ≥ cub).

Figure 3

 3 Figure 3Comparison on a medium-scale computer between toulbar2 using parallel HBFS (with burst mode) and cplex using multiple threads. The x-axis represents normalized time (with 0.2 corresponding to 720 seconds). The y-axis corresponds to normalized lower and upper bounds on 134 instances (with 1 corresponding to the optimum or best known cost, see the text description).

  is the depth of the node, ν.depth. HBFS always chooses the next open node to explore with minimum lower bound ν.lb (best-first principle) and, in case of ties, maximum depth ν.depth (depth-first principle) in the frontier. The minimum of all open node lower bounds, denoted lb(open), is a valid global lower bound (kept in clb at line 6) for the problem. HBFS also maintains the current upper bound (cub) as the cost of the best solution found so far by DFS (line 5). The search ends when the open list is empty or contains nodes with a lower bound greater than or equal to cub (line 2).

	Figure 1 A tree that is partially explored by DFS with a backtrack limit Z = 4. Nodes with a
	bold border are leaves, nodes with no border are placed in the open list after the backtrack bound is
	exceeded. Nodes are numbered in the order they are visited.
	DFS returns a new list of open nodes to be inserted in the BFS frontier.
	The pseudo-code of HBFS is given in Algorithm 1. The main procedure is in charge of the
	BFS frontier of open nodes. Here a node ν corresponds to a sequence of decisions ν.δ. The
	root node has an empty decision sequence (line 1). When a node is explored by DFS (line 5),
	an unassigned variable is chosen and a branching decision to either assign the variable to
	a chosen value (left branch, positive decision) or remove the value from the domain (right
	branch, negative decision) is taken. The number of decisions taken to reach a given node

in terms of quickly finding good solutions. As in BFS, HBFS maintains a frontier of open search nodes. It expands each open node using DFS with a limit on its number of backtracks. Each bounded ν

  Parallel HBFS-Worker. Initial call for Worker i: HBFS-Worker(k,i) with Zi = 1.

	/* or closing-node mes. in burst
	mode */
	Algorithm 3

). It can be explained by the fact that if it found a new solution then this improved upper bound will possibly imply more work in subsequent propagation made later when assessing the lower bound of each open node. This has the effect to slow-down the construction of the list Procedure HBFS-Worker(cub,rank) ; /* rank: Worker ID */ while (true) do openi := ∅ ; /* local open list of Worker i */ Receive an open node ν and solution cub ′ by Master ; /* Wait for message */ cub := min(cub, cub ′ ) ; /* Updates cub and best solution if any */ Restores state ν.δ, leading to assignment Aν , maintaining soft local consistency ; NodesRecompute := NodesRecompute + ν.depth ; 12 cub :=DFS(Aν ,cub,Zi) ; /* Increase Nodes ; put all right open branches in open i */ if (NodesRecompute > 0) then 13 if (NodesRecompute/Nodes > β and Zi ≤ N ) then Zi := 2 × Zi; 14 else if (NodesRecompute/Nodes < α and Zi ≥ 2) then Zi := Zi/2; 15 Send openi and best solution cub to the Master ;

Table 2

 2 ). EPS and HBFS-180 results on hard instances (with n variables and maximum domain size d). A '-' indicates that some (see #failed) subproblems could not be solved in less than 3, 600sec.

	instance	n	d	opt.	cub	l	av. time max. t. #fail(depth) EPS-180 HBFS-180
	linkage/pedigree19	259	5	4625	5684	5114	20.57	-	1 (4)	-	69.1
	linkage/pedigree40	274	6	7300	8838	5641	101.99	-	49 (21)	-	1680
	linkage/pedigree51	295	5	6406	6802	5798	0.61	497.38	0	499	5.7
	cpd/1BRS	38	178 4007610 4007679	956	2.94	38.90	0	44	37.5
	cpd/1CDL	38	170 3590514 3590825 1001	6.66	79.04	0	79	18.3
	cpd/1GVP	52	170 5196719 5196841	979	14.59	170.66	0	171	17.0
	maxcl./brock400_1	400	2	373	379	6010	63.95	-	12 ( 149 )	-	1812
	maxcl./brock400_2	400	2	371	379	5975	65.27	-	18 ( 149 )	-	880
	maxcl./san400_0.5_1 400	2	387	392	6073	5.07	414.96	0	3652	1220
											C P 2 0 2 2

In non-burst mode, all right branches share a common prefix in their ν.δ and only the deepest δ information need to be sent to the Master.

https://toulbar2.github.io/toulbar2 version 1.2.

We removed the largest instances keller6 and p_hat1500-1,2,3 from the original 66 DIMACS instances.

Linkage instances were further preprocessed by variable elimination limited to at most 8 neighbors[START_REF] Favier | Pairwise decomposition for combinatorial optimization in graphical models[END_REF].

.fr/degivry/Beldjilali22Supp.pdf Funding This work has been partially funded by the French "Agence Nationale de la Recherche", through grant ANR-19-P3IA-0004. It was performed using HPC resources from CALMIP (Grant 2022-P21010).

Carbon footprint

The experiments in this paper took approximately 17, 000 hours and emitted 68kg of CO2, with an estimate of 4g/h per core.

Supplementary Material https://miat.inrae