Amina Kemmar

Samir Loudni

Yahia Lebbah

Patrice Boizumault

Thierry Charnois

PREFIX-PROJECTION Global Constraint for Sequential Pattern Mining

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Mining useful patterns in sequential data is a challenging task. Sequential pattern mining is among the most important and popular data mining task with many real applications such as the analysis of web click-streams, medical or biological data and textual data. For effectiveness and efficiency considerations, many authors have promoted the use of constraints to focus on the most promising patterns according to the interests given by the final user. In line with [START_REF] Pei | Mining sequential patterns with constraints in large databases[END_REF], many efficient ad hoc methods have been developed but they suffer from a lack of genericity to handle and to push simultaneously sophisticated combination of various types of constraints. Indeed, new constraints have to be hand-coded and their combinations often require new implementations.

Recently, several proposals have investigated relationships between sequential pattern mining and constraint programming (CP) to revisit data mining tasks in a declarative and generic way [START_REF] Coquery | A SAT-based approach for discovering frequent, closed and maximal patterns in a sequence[END_REF][START_REF] Métivier | A constraint programming approach for mining sequential patterns in a sequence database[END_REF][START_REF] Kemmar | Mining relevant sequence patterns with cp-based framework[END_REF][START_REF] Negrevergne | Constraint-based sequence mining using constraint programming[END_REF]. The great advantage of these approaches is their flexibility. The user can model a problem and express his queries by specifying what constraints need to be satisfied. But, all these proposals are not effective enough because of their CP encoding. Consequently, the design of new efficient declarative models for mining useful patterns in sequential data is clearly an important challenge for CP.

To address this challenge, we investigate in this paper the other side of the cross fertilization between data-mining and constraint programming, namely how the CP framework can benefit from the power of candidate pruning mechanisms used in sequential pattern mining. First, we introduce the global constraint PREFIX-PROJECTION for sequential pattern mining. PREFIX-PROJECTION uses a concise encoding and its filtering relies on the principle of projected databases [START_REF] Pei | PrefixSpan: Mining sequential patterns by prefix-projected growth[END_REF]. The key idea is to divide the initial database into smaller ones projected on the frequent subsequences obtained so far, then, mine locally frequent patterns in each projected database by growing a frequent prefix. This global constraint utilizes the principle of prefix-projected database to keep only locally frequent items alongside projected databases in order to remove infrequent ones from the domains of variables. Second, we show how the concise encoding allows for a straightforward implementation of the frequency constraint (PREFIX-PROJECTION constraint) and constraints on patterns such as size, item membership and regular expressions and the simultaneous combination of them. Finally, experiments show that our approach clearly outperforms CP approaches and competes well with ad hoc methods on large datasets for mining frequent sequential patterns or patterns under various constraints. It is worth noting that the experiments show that our approach achieves scalability while it is a major issue of CP approaches.

The paper is organized as follows. Section 2 recalls preliminaries. Section 3 provides a critical review of ad hoc methods and CP approaches for sequential pattern mining. Section 4 presents the global constraint PREFIX-PROJECTION. Section 5 reports experiments we performed. Finally, we conclude and draw some perspectives.

Preliminaries

This section presents background knowledge about sequential pattern mining and constraint satisfaction problems.

Sequential Patterns

Let I be a finite set of items. The language of sequences corresponds to L I = I n where n ∈ N + . Definition 1 (sequence, sequence database). A sequence s over L I is an ordered list s 1 s 2 . . . s n , where s i , 1 ≤ i ≤ n, is an item. n is called the length of the sequence s. A sequence database SDB is a set of tuples (sid, s), where sid is a sequence identifier and s a sequence.

Definition 2 (subsequence, relation).

A sequence α = α 1 . . . α m is a subsequence of s = s 1 . . . s n , denoted by (α s), if m ≤ n and there exist integers 1 ≤ j 1 ≤ . . . ≤ j m ≤ n, such that α i = s ji for all 1 ≤ i ≤ m. We also say that α is contained in s or s is a super-sequence of α. For example, the sequence BABC is a super-sequence of AC : AC BABC . A tuple (sid, s) contains a sequence α, if α s.

The cover of a sequence p in SDB is the set of all tuples in SDB in which p is contained. The support of a sequence p in SDB is the number of tuples in SDB which contain p. Definition 3 (coverage, support). Let SDB be a sequence database and p a sequence. cover SDB (p)={(sid, s) ∈ SDB | p s} and sup SDB (p) = #cover SDB (p). Definition 4 (sequential pattern). Given a minimum support threshold minsup, every sequence p such that sup SDB (p) ≥ minsup is called a sequential pattern [START_REF] Agrawal | Mining sequential patterns[END_REF]. p is said to be frequent in SDB.

SPM under Constraints

In this section, we define the problem of mining sequential patterns in a sequence database satisfying user-defined constraints Then, we review the most usual constraints for the sequential mining problem [START_REF] Pei | Mining sequential patterns with constraints in large databases[END_REF]. Problem statement. Given a constraint C(p) on pattern p and a sequence database SDB, the problem of constraint-based pattern mining is to find the complete set of patterns satisfying C(p). In the following, we present different types of constraints that we explicit in the context of sequence mining. All these constraints will be handled by our concise encoding (see Sections 4.2 and 4.5).

-The minimum size constraint size(p, min ) states that the number of items of p must be greater than or equal to min .

-The item constraint item(p, t) states that an item t must belong (or not) to a pattern p.

-The regular expression constraint [START_REF] Garofalakis | Mining sequential patterns with regular expression constraints[END_REF] reg(p, exp) states that a pattern p must be accepted by the deterministic finite automata associated to the regular expression exp.

Projected Databases

We now present the necessary definitions related to the concept of projected databases [START_REF] Pei | PrefixSpan: Mining sequential patterns by prefix-projected growth[END_REF].

Definition 6 (prefix, projection, suffix). Let β = β 1 . . . β n and α = α 1 . . . α m be two sequences, where m ≤ n. -Sequence α is called the prefix of β iff ∀i ∈ [1..m], α i = β i .
-Sequence β = β 1 . . . β n is called the projection of some sequence s w.r.t. α, iff (1) β s, (2) α is a prefix of β and (3) there exists no proper super-sequence β of β such that β s and β also has α as prefix.

-Sequence γ = β m+1 . . . β n is called the suffix of s w.r.t. α. With the standard concatenation operator "concat", we have β = concat(α, γ).

Definition 7 (projected database). Let SDB be a sequence database, the α-projected database, denoted by SDB| α , is the collection of suffixes of sequences in SDB w.r.t. prefix α.

[14] have proposed an efficient algorithm, called PrefixSpan, for mining sequential patterns based on the concept of projected databases. It proceeds by dividing the initial database into smaller ones projected on the frequent subsequences obtained so far; only their corresponding suffixes are kept. Then, sequential patterns are mined in each projected database by exploring only locally frequent patterns.

Example 2. Let us consider the sequence database of Table 1 with minsup = 2. PrefixSpan starts by scanning SDB 1 to find all the frequent items, each of them is used as a prefix to get projected databases. For SDB 1 , we get 3 disjoint subsets w.r.t. the prefixes A , B , and C . For instance, SDB 1 | A consists of 3 suffix sequences: {(1, BCBC ), (2, BC ), (3, B )}. Consider the projected database SDB 1 | <A> , its locally frequent items are B and C. Thus, SDB 1 | <A> can be recursively partitioned into 2 subsets w.r.t. the two prefixes AB and AC . The AB -and AC -projected databases can be constructed and recursively mined similarly. The processing of a αprojected database terminates when no frequent subsequence can be generated.

Proposition 1 establishes the support count of a sequence γ in SDB| α [START_REF] Pei | PrefixSpan: Mining sequential patterns by prefix-projected growth[END_REF]: Proposition 1 (Support count). For any sequence γ in SDB with prefix α and suffix

β s.t. γ = concat(α, β), sup SDB (γ) = sup SDB|α (β).
This proposition ensures that only the sequences in SDB grown from α need to be considered for the support count of a sequence γ. Furthermore, only those suffixes with prefix α should be counted.

CSP and Global Constraints

A Constraint Satisfaction Problem (CSP) consists of a set X of n variables, a domain D mapping each variable X i ∈ X to a finite set of values D(X i ), and a set of constraints C. An assignment σ is a mapping from variables in X to values in their domains: ∀X i ∈ X, σ(X i ) ∈ D(X i ). A constraint c ∈ C is a subset of the cartesian product of the domains of the variables that are in c. The goal is to find an assignment such that all constraints are satisfied. Domain consistency (DC). Constraint solvers typically use backtracking search to explore the space of partial assignments. At each assignment, filtering algorithms prune the search space by enforcing local consistency properties like domain consistency. A constraint c on X is domain consistent, if and only if, for every X i ∈ X and for every d i ∈ D(X i ), there is an assignment σ satisfying c such that σ(X i ) = d i . Such an assignment is called a support. Global constraints provide shorthands to often-used combinatorial substructures. We present two global constraints. Let X = X 1 , X 2 , ..., X n be a sequence of n variables. Let V be a set of values, l and u be two integers s.t. 0 ≤ l ≤ u ≤ n, the constraint Among(X, V, l, u) states that each value a ∈ V should occur at least l times and at most u times in X [START_REF] Beldiceanu | Introducing global constraints in CHIP[END_REF]. Given a deterministic finite automaton A, the constraint Regular(X, A) ensures that the sequence X is accepted by A [START_REF] Pesant | A regular language membership constraint for finite sequences of variables[END_REF].

This section provides a critical review of ad hoc methods and CP approaches for SPM.

3.1 Ad hoc Methods for SPM GSP [START_REF] Srikant | Mining sequential patterns: Generalizations and performance improvements[END_REF] was the first algorithm proposed to extract sequential patterns. It uses a generate-and test approach. Later, two major classes of methods have been proposed: -Depth-first search based on a vertical database format e.g. cSpade incorporating contraints (max-gap, max-span, length) [START_REF] Zaki | Sequence mining in categorical domains: Incorporating constraints[END_REF], SPADE [START_REF] Zaki | SPADE: An efficient algorithm for mining frequent sequences[END_REF] or SPAM [START_REF] Ayres | Sequential pattern mining using a bitmap representation[END_REF].

-Projected pattern growth such as PrefixSpan [START_REF] Pei | PrefixSpan: Mining sequential patterns by prefix-projected growth[END_REF] and its extensions, e.g. CloSpan for mining closed sequential patterns [START_REF] Yan | CloSpan: Mining closed sequential patterns in large databases[END_REF] or Gap-BIDE [START_REF] Li | Efficient mining of gap-constrained subsequences and its various applications[END_REF] tackling the gap constraint.

In [START_REF] Garofalakis | Mining sequential patterns with regular expression constraints[END_REF], the authors proposed SPIRIT based on GSP for SPM with regular expressions. Later, [START_REF] Trasarti | Sequence mining automata: A new technique for mining frequent sequences under regular expressions[END_REF] introduces Sequence Mining Automata (SMA), a new approach based on a specialized kind of Petri Net. Two variants of SMA were proposed: SMA-1P (SMA one pass) and SMA-FC (SMA Full Check). SMA-1P processes by means of the SMA all sequences one by one, and enters all resulting valid patterns in a hash table for support counting, while SMA-FC allows frequency based pruning during the scan of the database. Finally, [START_REF] Pei | Mining sequential patterns with constraints in large databases[END_REF] provides a survey for other constraints such as regular expressions, length and aggregates. But, all these proposals, though efficient, are ad hoc methods suffering from a lack of genericity. Adding new constraints often requires to develop new implementations.

CP Methods for SPM

Following the work of [START_REF] Guns | Itemset mining: A constraint programming perspective[END_REF] for itemset mining, several methods have been proposed to mine sequential patterns using CP. Proposals. [START_REF] Coquery | A SAT-based approach for discovering frequent, closed and maximal patterns in a sequence[END_REF] have proposed a first SAT-based model for discovering a special class of patterns with wildcards1 in a single sequence under different types of constraints (e.g. frequency, maximality, closedness). [START_REF] Métivier | A constraint programming approach for mining sequential patterns in a sequence database[END_REF] have proposed a CSP model for SPM. Each sequence is encoded by an automaton capturing all subsequences that can occur in it. [START_REF] Kemmar | Mining relevant sequence patterns with cp-based framework[END_REF] have proposed a CSP model for SPM with wildcards. They show how some constraints dealing with local patterns (e.g. frequency, size, gap, regular expressions) and constraints defining more complex patterns such as relevant subgroups [START_REF] Novak | Supervised descriptive rule discovery: A unifying survey of contrast set, emerging pattern and subgroup mining[END_REF] and top-k patterns can be modeled using a CSP. [START_REF] Negrevergne | Constraint-based sequence mining using constraint programming[END_REF] have proposed two CP encodings for the SPM. The first one uses a global constraint to encode the subsequence relation (denoted global-p.f), while the second one encodes explicitly this relation using additional variables and constraints (denoted decomposed-p.f).

All these proposals use reified constraints to encode the database. A reified constraint associates a boolean variable to a constraint reflecting whether the constraint is satisfied (value 1) or not (value 0). For each sequence s of SDB, a reified constraint, stating whether (or not) the unknown pattern p is a subsequence of s, is imposed: (S s = 1) ⇔ (p s). A great consequence is that the encoding of the frequency measure is straightforward: f req(p) = s∈SDB S s . But such an encoding has a major drawback since it requires (m = #SDB) reified constraints to encode the whole database. This constitutes a strong limitation of the size of the databases that could be managed.

Most of these proposals encode the subsequence relation (p s) using variables P os s,j (s ∈ SDB and 1 ≤ j ≤ ) to determine a position where p occurs in s. Such an encoding requires a large number of additional variables (m× ) and makes the labeling computationally expensive. In order to address this drawback, [START_REF] Negrevergne | Constraint-based sequence mining using constraint programming[END_REF] have proposed a global constraint exists-embedding to encode the subsequence relation, and used projected frequency within an ad hoc specific branching strategy to keep only frequent items before branching over the variables of the pattern. But, this encoding still relies on reified constraints and requires to impose m exists-embedding global constraints.

So, we propose in the next section the PREFIX-PROJECTION global constraint that fully exploits the principle of projected databases to encode both the subsequence relation and the frequency constraint. PREFIX-PROJECTION does not require any reified constraints nor any extra variables to encode the subsequence relation. As a consequence, usual SPM constraints (see Section 2.2) can be encoded in a straightforward way using directly the (global) constraints of the CP solver.

PREFIX-PROJECTION Global Constraint

This section presents the PREFIX-PROJECTION global constraint for the SPM problem.

A Concise Encoding

Let P be the unknown pattern of size we are looking for. The symbol 2 stands for an empty item and denotes the end of a sequence. The unknown pattern P is encoded with a sequence of variables P 1 , P 2 , . . . , P s.t. ∀i ∈ [1 . . . ], D(P i ) = I ∪ {2}. There are two basic rules on the domains:

1. To avoid the empty sequence, the first item of P must be non empty, so (2 ∈ D 1 ). 2. To allow patterns with less than items, we impose that ∀i ∈ [1..( -1)], (P i = 2) → (P i+1 = 2).

Definition and Consistency Checking

The global constraint PREFIX-PROJECTION ensures both subsequence relation and minimum frequency constraint. Proof: This is a direct consequence of proposition 1. We have straightforwardly sup SDB (σ) = sup SDB|σ ( ) = #SDB| σ . Thus, suffixes of SDB| σ are supports of σ in the constraint PREFIX-PROJECTION (P, SDB, minsup), provided that #SDB| σ ≥ minsup. 2

The following proposition characterizes values in the domain of unassigned (i.e. future) variable P i+1 that are consistent with the current assignment of variables P 1 , ..., P i . Proposition 3. Let σ2 = d 1 , . . . , d i be a current assignment of variables P 1 , . . . , P i , P i+1 be a future variable. A value d ∈ D(P i+1 ) appears in a solution for PREFIX-PROJECTION (P, SDB, minsup) if and only if d is a frequent item in SDB| σ :

#{(sid, γ)|(sid, γ) ∈ SDB| σ ∧ d γ} ≥ minsup
Proof: Suppose that value d ∈ D(P i+1 ) occurs in SDB| σ more than minsup. From proposition 1, we have sup SDB (concat(σ, d )) = sup SDB|σ ( d ). Hence, the assignment σ ∪ d satisfies the constraint, so d ∈ D(P i+1 ) participates in a solution. 2 Anti-monotonicity of the frequency measure. If a pattern p is not frequent, then any pattern p satisfying p p is not frequent. From proposition 3 and according to the anti-monotonicity property, we can derive the following pruning rule: Proposition 4. Let σ = d 1 , . . . , d i be a current assignment of variables P 1 , . . . , P i . All values d ∈ D(P i+1 ) that are locally not frequent in SDB| σ can be pruned from the domain of variable P i+1 . Moreover, these values d can also be pruned from the domains of variables P j with j ∈ [i + 2, . . . , ].

Proof: Let σ = d 1 , . . . , d i be a current assignment of variables P 1 , . . . , P i . Let

d ∈ D(P i+1 ) s.t. σ = concat(σ, d ). Suppose that d is not frequent in SDB| σ .
According to proposition 1, sup SDB|σ ( d ) = sup SDB (σ ) < minsup, thus σ is not frequent. So, d can be pruned from the domain of P i+1 . Suppose that the assignment σ has been extended to concat(σ, α), where α corresponds to the assignment of variables P j (with j > i). If d ∈ D(P i+1 ) is not frequent, it is straightforward that sup SDB|σ (concat(α, d )) ≤ sup SDB|σ ( d ) < minsup. Thus, if d is not frequent in SDB| σ , it will be also not frequent in SDB| concat(σ,α) . So, d can be pruned from the domains of P j with j ∈ [i + 2, . . . , ]. 2 Example 3. Consider the sequence database of Table 1 with minsup = 2. Let P = P 1 , P 2 , P 3 with D(P 1 ) = I and D(P 2 ) = D(P 3 ) = I ∪ {2}. Suppose that σ(P 1 ) = A, PREFIX-PROJECTION(P, SDB, minsup) will remove values A and D from D(P 2 ) and D(P 3 ), since the only locally frequent items in SDB 1 | <A> are B and C. Proposition 4 guarantees that any value (i.e. item) d ∈ D(P i+1 ) present but not frequent in SDB| σ does not need to be considered when extending σ, thus avoiding searching over it. Clearly, our global constraint encodes the anti-monotonicity of the frequency measure in a simple and elegant way, while CP methods for SPM have difficulties to handle this property. In [START_REF] Negrevergne | Constraint-based sequence mining using constraint programming[END_REF], this is achieved by using very specific propagators and branching strategies, making the integration quite complex (see [START_REF] Negrevergne | Constraint-based sequence mining using constraint programming[END_REF]).

Building the projected databases.

The key issue of our approach lies in the construction of the projected databases. When projecting a prefix, instead of storing the whole suffix as a projected subsequence, one can represent each suffix by a pair (sid, start) where sid is the sequence identifier and start is the starting position of the projected suffix in the sequence sid. For instance, let

}. This is the principle of pseudo-projection, adopted in PrefixSpan, exploited during the filtering step of our PREFIX-PROJECTION global constraint. Algorithm 1 details this principle. It takes as input a set of projected sequences P rojSDB and a prefix α. The algorithm processes all the pairs (sid, start) of P rojSDB one by one (line 2), and searches for the lowest location of α in the sequence s corresponding to the sid of that sequence in SDB (lines 6-8).

In the worst case, PROJECTSDB processes all the items of all sequences. So, the time complexity is O( × m), with m = #SDB and is the length of the longest sequence in SDB. The worst case space complexity of pseudo-projection is O(m), since we need to store for each sequence only a pair (sid, start), while for the standard projection the space complexity is O(m× ). Clearly, the pseudo-projection takes much less space than the standard projection.

Filtering

Ensuring DC on PREFIX-PROJECTION(P, SDB, minsup) is equivalent to finding a sequential pattern of length ( -1) and then checking whether this pattern remains a frequent pattern when extended to any item d in D(P ). Thus, finding such an assignment (i.e. support) is as much as difficult than the original problem of sequential pattern mining. [START_REF] Yang | Computational aspects of mining maximal frequent patterns[END_REF] has proved that the problem of counting the number of maximal3 frequent patterns in a database of sequences is #P-complete, thereby proving the NP-hardness of the problem of mining maximal frequent sequences. The difficulty is due to the exponential number of candidates that should be parsed to find the frequent patterns. Thus, finding, for every variable P i ∈ P and for every d i ∈ D(P i ), an assignment σ satisfying PREFIX-PROJECTION(P, SDB, minsup) s.t. σ(P i ) = d i is of exponential nature.

So, the filtering of the PREFIX-PROJECTION constraint maintains a consistency lower than DC. This consistency is based on specific properties of the projected databases Algorithm 2: FILTER-PREFIX-PROJECTION(SDB, σ, i, P , minsup) Data: SDB: initial database; σ: current prefix σ(P1), . . . , σ(Pi) ; minsup: the minimum support threshold; PSDB: internal data structure of PREFIX-PROJECTION for storing pseudo-projected databases begin 1

if (see Proposition 3), and anti-monotonicity of the frequency constraint (see Proposition 4), and resembles forward-checking regarding Proposition 3. PREFIX-PROJECTION is considered as a global constraint, since all variables share the same internal data structures that awake and drive the filtering.

(i ≥ 2 ∧ σ(Pi) = 2)
Algorithm 2 describes the pseudo-code of the filtering algorithm of the PREFIX-PROJECTION constraint. It is an incremental filtering algorithm that should be run when some i first variables are assigned according to the following lexicographic ordering P 1 , P 2 , . . . , P of variables of P . It exploits internal data-structures enabling to enhance the filtering algorithm. More precisely, it uses an incremental data structure, denoted PSDB, that stores the intermediate pseudo-projections of SDB, where PSDB i (i ∈ [0, . . . , ]) corresponds to the σ-projected database of the current partial assignment σ = σ(P 1 ), . . . , σ(P i ) (also called prefix) of variables P 1 , . . . , P i , and PSDB 0 = {(sid, 1)|(sid, s) ∈ SDB} is the initial pseudo-projected database of SDB (case where σ = ). It also uses a hash table indexing the items I into integers (1 . . . #I) for an efficient support counting over items (see function getFreqItems).

Algorithm 2 takes as input the current partial assignment σ = σ(P 1 ), . . . , σ(P i ) of variables P 1 , . . . , P i , the length i of σ (i.e. position of the last assigned variable in P ) and the minimum support threshold minsup. It starts by checking if the last assigned variable P i is instantiated to 2 (line 1). In this case, the end of sequence is reached (since value 2 can only appear at the end) and the sequence σ(P 1 ), . . . , σ(P i ) constitutes a frequent pattern in SDB; hence the algorithm sets the remaining ( -i) unassigned variables to 2 and returns true (lines 2-4). Otherwise, the algorithm computes incrementally PSDB i from PSDB i-1 by calling function PROJECTSDB (see Algorithm 1). Then, it checks in line 6 whether the current assignment σ is a legal prefix for the constraint (see Proposition 2). This is done by computing the size of PSDB i . If this size is less than minsup, we stop growing σ and we return false. Otherwise, the algorithm computes the set of locally frequent items F I in PSDB i by calling function getFreqItems (line 8).

Function getFreqItems processes all the entries of the pseudo-projected database one by one, counts the number of first occurrences of items a (i.e. SupCount[a]) in each entry (sid, start), and keeps only the frequent ones (lines [START_REF] Novak | Supervised descriptive rule discovery: A unifying survey of contrast set, emerging pattern and subgroup mining[END_REF][START_REF] Pei | PrefixSpan: Mining sequential patterns by prefix-projected growth[END_REF][START_REF] Pei | Mining sequential patterns with constraints in large databases[END_REF][START_REF] Pesant | A regular language membership constraint for finite sequences of variables[END_REF][START_REF] Srikant | Mining sequential patterns: Generalizations and performance improvements[END_REF][START_REF] Trasarti | Sequence mining automata: A new technique for mining frequent sequences under regular expressions[END_REF][START_REF] Yan | CloSpan: Mining closed sequential patterns in large databases[END_REF][START_REF] Yang | Computational aspects of mining maximal frequent patterns[END_REF][START_REF] Zaki | Sequence mining in categorical domains: Incorporating constraints[END_REF]. This is done by using ExistsItem data structure. After the whole pseudo-projected database has been processed, the frequent items are returned (line 22), and Algorithm 2 updates the current domains of variables P j with j ≥ (i + 1) by pruning inconsistent values, thus avoiding searching over not frequent items (lines 9-11).

(m × + m × ( + d) + × d) = O(m × + m × d + × d).
The space complexity of the filtering algorithm lies in the storage of the PSDB internal data structure. In the worst case, we have to store pseudo-projected databases. Since each pseudo-projected database requires O(m), the worst case space complexity is O(m × ). 2

Encoding of SPM Constraints

Usual SPM constraints (see Section 2.2) can be reformulated in a straightforward way. Let P be the unknown pattern.

-Minimum size constraint: size(P, min ) ≡ i= min i=1

(P i = ) -Item constraint: let V be a subset of items, l and u two integers s.t. 0 ≤ l ≤ u ≤ . item(P, V ) ≡ t∈V Among(P, {t}, l, u) enforces that items of V should occur at least l times and at most u times in P . To forbid items of V to occur in P , l and u must be set to 0.

-Regular expression constraint: let A reg be the deterministic finite automaton encoding the regular expression exp. reg(P, exp) ≡ Regular(P, A reg ).

Experimental Evaluation

This section reports experiments on several real-life datasets from [START_REF] Fournier-Viger | SPMF: A Java Open-Source Pattern Mining Library[END_REF][START_REF] Béchet | Sequential pattern mining to discover relations between genes and rare diseases[END_REF][START_REF] Trasarti | Sequence mining automata: A new technique for mining frequent sequences under regular expressions[END_REF] of large size having varied characteristics and representing different application domains (see Table 2). Our objective is (1) to compare our approach to existing CP methods as well as to state-of-the-art methods for SPM in terms of scalability which is a major issue of existing CP methods, (2) to show the flexibility of our approach allowing to handle different constraints simultaneously.

Experimental protocol. The implementation of our approach was carried out in the Gecode solver 4 . All experiments were conducted on a machine with a processor Intel X5670 and 24 GB of memory. A time limit of 1 hour has been used. For each dataset, we varied the minsup threshold until the methods are not able to complete the extraction of all patterns within the time limit. was set to the length of the longest sequence of SDB. The implementation and the datasets used in our experiments are available online 5 . We compare our approach (indicated by PP) with:

1. two CP encodings [START_REF] Negrevergne | Constraint-based sequence mining using constraint programming[END_REF], the most efficient CP methods for SPM: global-p.f and decomposed-p.f; 2. state-of-the-art methods for SPM : PrefixSpan and cSpade; 3. SMA [START_REF] Trasarti | Sequence mining automata: A new technique for mining frequent sequences under regular expressions[END_REF] for SPM under regular expressions. We used the author's cSpade implementation 6 for SPM, the publicly available implementations of PrefixSpan by Y. Tabei7 and the SMA implementation 8 for SPM under regular expressions. The implementation 9 of the two CP encodings was carried out in the Gecode solver. All methods have been executed on the same machine. (a) Comparing with CP Methods for SPM. First we compare PP with the two CP encodings global-p.f and decomposed-p.f (see Section 3.2). CPU times (in logscale for BIBLE, Kosarak and PubMed) of the three methods are shown on Fig. 1. First, decomposed-p.f is the least performer method. On all the datasets, it fails to complete the extraction within the time limit for all values of minsup we considered. Second, PP largely dominates global-p.f on all the datasets: PP is more than an order of magnitude faster than global-p.f. The gains in terms of CPU times are greatly amplified for low values of minsup. On BIBLE (resp. PubMed), the speed-up is 84.4 (resp. 33.5) for minsup equal to 1%. Another important observation that can be made is that, on most of the datasets (except BIBLE and Kosarak), global-p.f is not able to mine for patterns at very low frequency within the time limit. For example on FIFA, PP is able to complete the extraction for values of minsup up to 6% in 1, 457 seconds, while global-p.f fails to complete the extraction for minsup less than 10%.

To complement the results given by Fig. 1, Table 3 reports for different datasets and different values of minsup, the number of calls to the propagate routine of Gecode (column 5), and the number of nodes of the search tree (column 6). First, PP explores less nodes than global-p.f. But, the difference is not huge (gains of 45% and 33% on FIFA and BIBLE respectively). Second, our approach is very effective in terms of number of propagations. For PP, the number of propagations remains small (in thousands for small values of minsup) compared to global-p.f (in millions). This is due to the huge number of reified constraints used in global-p.f to encode the subsequence relation. On the contrary, our PREFIX-PROJECTION global constraint does not require any reified constraints nor any extra variables to encode the subsequence relation. (b) Comparing with ad hoc Methods for SPM. Our second experiment compares PP with state-of-the-art methods for SPM. Fig. 2 shows the CPU times of the three methods. First, cSpade obtains the best performance on all datasets (except on Protein). However, PP exhibits a similar behavior as cSpade, but it is less faster (not counting the highest values of minsup). The behavior of cSpade on Protein is due to the vertical representation format that is not appropriated in the case of databases having large sequences and small number of distinct items, thus degrading the performance of the mining process. Second, PP which also uses the concept of projected databases, clearly outperforms PrefixSpan on all datasets. This is due to our filtering algorithm combined together with incremental data structures to manage the projected databases.

On FIFA, PrefixSpan is not able to complete the extraction for minsup less than 12%, while our approach remains feasible until 6% within the time limit. On Protein, PrefixSpan fails to complete the extraction for all values of minsup we considered. These results clearly demonstrate that our approach competes well with state-of-the-art methods for SPM on large datasets and achieves scalability while it is a major issue of existing CP approaches.

(c) SPM under size and item constraints. Our third experiment aims at assessing the interest of pushing simultaneously different types of constraints. We impose on the PubMed dataset usual constraints such as the minimum frequency and the minimum size constraints and other useful constraints expressing some linguistic knowledge such as the item constraint. The goal is to retain sequential patterns which convey linguistic regularities (e.g., gene -rare disease relationships) [START_REF] Béchet | Sequential pattern mining to discover relations between genes and rare diseases[END_REF]. The size constraint allows to remove patterns that are too small w.r.t. the number of items (number of words) to be relevant patterns. We tested this constraint with min set to 3. The item constraint imposes that the extracted patterns must contain the item GENE and the item DISEASE. As no ad hoc method exists for this combination of constraints, we only compare PP with global-p.f. Fig. 3 shows the CPU times and the number of sequential patterns extracted with and without constraints. First, pushing simultaneously the two constraints enables to reduce significantly the number of patterns. Moreover, the CPU times for PP decrease slightly whereas for global-p.f (with and without constraints), they are almost the same. This is probably due to the weak communication between the m exists-embedding reified global constraints and the two constraints. This reduces significantly the quality of the whole filtering. Second (see Table 4), when considering the two constraints, PP clearly dominates global-p.f (speed-up value up to 51.5).

Moreover, the number of propagations performed by PP remains very small as com- pared to global-p.f. Fig. 3c compares the two methods under the minimum size constraint for different values of min , with minsup fixed to 1%. Once again, PP is always the most performer method (speed-up value up to 53.1). These results also confirm what we observed previously, namely the weak communication between reified global constraints and constraints imposed on patterns (i.e., size and item constraints).

(d) SPM under regular constraints. Our last experiment compares PP-REG against two variants of SMA: SMA-1P (SMA one pass) and SMA-FC (SMA Full Check). Two datasets are considered from [START_REF] Trasarti | Sequence mining automata: A new technique for mining frequent sequences under regular expressions[END_REF]: one synthetic dataset (data-200k), and one real-life dataset (Protein). On the synthetic dataset, our approach is very effective. For RE14, our method is more than an order of magnitude faster than SMA. On Protein, the gap between the 3 methods shrinks, but our method remains effective. For the particular case of RE2, the Regular constraint can be substituted by restricting the domain of the first and third variables to {S, T } and {R, K} respectively (denoted as PP-SRE), thus improving performances.

Conclusion

We have proposed the global constraint PREFIX-PROJECTION for sequential pattern mining. PREFIX-PROJECTION uses a concise encoding and provides an efficient filtering based on specific properties of the projected databases, and anti-monotonicity of the frequency constraint. When this global constraint is integrated into a CP solver, it enables to handle several constraints simultaneously. Some of them like size, item membership and regular expression are considered in this paper. Another point of strength, is that, contrary to existing CP approaches for SPM, our global constraint does not require any reified constraints nor any extra variables to encode the subsequence relation. Finally, although PREFIX-PROJECTION is well suited for constraints on sequences, it would require to be adapted to handle constraints on subsequence relations like gap. Experiments performed on several real-life datasets show that our approach clearly outperforms existing CP approaches and competes well with ad hoc methods on large datasets and achieves scalability while it is a major issue of CP approaches. As future work, we intend to handle constraints on set of sequential patterns such as closedness, relevant subgroup and skypattern constraints.

Proposition 5 .

 5 In the worst case, filtering with PREFIX-PROJECTION global constraint can be achieved in O(m × + m × d + × d). The worst case space complexity of PREFIX-PROJECTION is O(m × ). Proof: Let be the length of the longest sequence in SDB, m = #SDB, and d = #I. Computing the pseudo-projected database PSDB i can be done in O(m × ): for each sequence (sid, s) of SDB, checking if σ occurs in s is O( ) and there are m sequences. The total complexity of function GETFREQITEMS is O(m × ( + d)). Lines (9-11) can be achieved in O( × d). So, the whole complexity is O

Fig. 1 :

 1 Fig. 1: Comparing PP with global-p.f for SPM on real-life datasets: CPU times.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Comparing PREFIX-PROJECTION with state-of-the-art algorithms for SPM.

  For data-200k, we used two RE: RE10 ≡ A * B(B|C)D * EF * (G|H)I * and RE14 ≡ A * (Q|BS * (B|C))D * E(I|S) * (F |H)G * R. For Protein, we used RE2 ≡ (S|T ) . (R|K) (where . represents any symbol). Fig. 4 reports CPU-times comparison.

Fig. 4 :

 4 Fig. 4: Comparing PREFIX-PROJECTION with SMA for SPM under RE constraint.

Table 1 :

 1 SDB 1 : a sequence database example.Example 1. Table1represents a sequence database of four sequences where the set of items is I = {A, B, C, D}. Let the sequence p = AC . We have cover SDB1 (p) = {(1, s 1 ), (2, s 2 )}. If we consider minsup = 2, p = AC is a sequential pattern because sup SDB1 (p) ≥ 2.

	sid Sequence
	1	ABCBC
	2	BABC
	3	AB
	4	BCD

Definition 5 (sequential pattern mining (SPM)). Given a sequence database SDB and a minimum support threshold minsup. The problem of sequential pattern mining is to find all patterns p such that sup SDB (p) ≥ minsup.

  Algorithm 1: PROJECTSDB(SDB, P rojSDB, α)

		Data: SDB: initial database; P rojSDB: projected sequences; α: prefix
		begin
	1	SDB|α ← ∅ ;
	2	for each pair (sid, start) ∈ P rojSDB do
	3	s ← SDB[sid] ;
	4	posα ← 1; poss ← start ;
	5	while (posα ≤ #α ∧ poss ≤ #s) do
	6	if (α[posα] = s[poss]) then
	7	posα ← posα + 1 ;
	8	poss ← poss + 1 ;
	9	if (posα = #α + 1) then
	10	SDB|α ← SDB|α ∪ {(sid, poss)}
	11	return SDB|α ;

us consider the sequence database of Table 1. As shown in example 2, SDB| A consists of 3 suffix sequences: {(1, BCBC ), (2, BC ), (3, B )}. By using the pseudoprojection, SDB| A can be represented by the following three pairs: {(1, 2), (2, 3),

Table 2 :

 2 Dataset Characteristics.

	dataset	#SDB #I avg (#s) maxs∈SDB (#s)	type of data
	Leviathen 5834 9025 33.81	100	book
	Kosarak	69999 21144 7.97	796	web click stream
	FIFA	20450 2990 34.74	100	web click stream
	BIBLE	36369 13905 21.64	100	bible
	Protein	103120 24	482	600	protein sequences
	data-200K 200000 20	50	86	synthetic dataset
	PubMed	17527 19931	29	198	bio-medical text

Table 3 :

 3 PP vs. global-p.f.

				CPU times (s)	#PROPAGATIONS	#NODES
				PP global-p.f	PP global-p.f	PP global-p.f
		20	938	8.16	129.54 1884	11649290 1025	1873
		18	1743	13.39	222.68 3502	19736442 1922	3486
	FIFA	16 14	3578 7313	24.39 44.08	396.11 7181 704 14691	35942314 3923 65522076 8042	7151 14616
		12	16323	86.46	1271.84 32820	126187396 18108	32604
		10	40642	185.88	2761.47 81767	266635050 45452	81181
		10	174	1.98	105.01	363	4189140 235	348
		8	274	2.47	153.61	575	5637671 362	548
	BIBLE	6 4	508 1185	3.45 5.7	270.49 1065 552.62 2482	8592858 669 15379396 1575	1016 2371
		2	5311	15.05	1470.45 11104	39797508 7048	10605
		1	23340	41.4	3494.27 49057	98676120 31283	46557
		5	2312	8.26	253.16 4736	15521327 2833	4619
		4	3625	11.17	340.24 7413	20643992 4428	7242
	PubMed	3 2	6336 13998	16.51 28.91	536.96 12988 955.54 28680	29940327 7757 50353208 17145	12643 27910
		1	53818	77.01	2581.51 110133	124197857 65587	107051
		99.99	127	165.31	219.69	264	26731250 172	221
		99.988	216	262.12	411.83	451	44575117 293	390
	Protein	99.986 99.984	384 631	467.96 753.3	909.47 1443.92 1322 805	80859312 514 132238827 845	679 1119
		99.982	964	1078.73	2615 2014	201616651 1284	1749
		99.98	2143	2315.65	-4485	-2890	-

Table 4 :

 4 PP vs. global-p.f under minimum size and item constraints.

	Dataset minsup (%) #PATTERNS	CPU times (s) PP global-p.f	#PROPAGATIONS PP global-p.f	#NODES PP global-p.f
		5	279	6.76	252.36 7878	12234292 2285	4619
		4	445	8.81	339.09 12091	16475953 3618	7242
	PubMed	3 2	799 1837	12.35 20.41	535.32 20268 953.32 43088	24380096 6271 42055022 13888	12643 27910
		1	7187	49.98	2574.42 157899	107978568 52508	107051
	data-200k (RE10)		data-200k (RE14)		Protein (RE2)

A wildcard is a special symbol that matches any item of I including itself.

We indifferently denote σ by d1, . . . , di or by σ(P1), . . . , σ(Pi) .

A sequential pattern p is maximal if there is no sequential pattern q such that p q.

http://www.gecode.org

https://sites.google.com/site/prefixprojection4cp/

http://www.cs.rpi.edu/~zaki/www-new/pmwiki.php/Software/

https://code.google.com/p/prefixspan/

http://www-kdd.isti.cnr.it/SMA/

https://dtai.cs.kuleuven.be/CP4IM/cpsm/

Acknowledgments. The authors would like to thank the anonymous referees for their valuable comments. This work is partly supported by the ANR (French Research National Agency) funded projects Hybride ANR-11-BS002-002.


Abdelkader Ouali

David Allouche

Simon De Givry

Samir Loudni

Yahia Lebbah

Francisco Eckhardt

Lakhdar Loukil

Iterative Decomposition Guided Variable Neighborhood Search for Graphical Model Energy Minimization

Graphical models factorize a global probability distribution/energy function as the product/sum of local functions. A major inference task, known as MAP in Markov Random Fields and MPE in Bayesian Networks, is to find a global assignment of all the variables with maximum a posteriori probability/minimum energy. A usual distinction on MAP solving methods is complete/incomplete, i.e. the ability to prove optimality or not. Most complete methods rely on tree search, while incomplete methods rely on local search. Among them, we study Variable Neighborhood Search (VNS) for graphical models. In this paper, we propose an iterative approach above VNS which uses (partial) tree search inside its local neighborhood exploration. The resulting hybrid method offers a good compromise between completeness and anytime behavior than existing tree search methods while still being competitive for proving optimality. We further propose a parallel version of our method improving its anytime behavior on difficult instances coming from a large graphical model benchmark. Last we experiment on the challenging minimum energy problem found in Computational Protein Design, showing the practical benefit of our parallel version.

INTRODUCTION

Probabilistic graphical models (Koller and Friedman, 2009) are formed by variables linked to each other by stochastic relationships. They enable to model complex systems with heterogeneous data and to capture uncertainty. Graphical models have been applied in a wide range of areas such as image analysis, speech recognition, bioinformatics, and ecology.

We focus on models with discrete variables like Markov Random Field and Bayesian Network. Our goal is to find a global assignment of all the variables with maximum a posteriori probability. This optimization task defines an NP-complete problem (Shimony, 1994). Solving methods can be categorized in two groups: exact and local search methods. Exact methods rely on tree search, variable elimination, linear programming, or a combination of them [START_REF] Marinescu | [END_REF][START_REF] Otten | [END_REF]Allouche et al., 2015). Graph-cut and message-passing algorithms like loopy belief propagation and variational approaches (Kolmogorov, 2006;[START_REF] Wainwright | [END_REF]Sontag et al., 2008;2012;Wang and Koller, 2013) are exact only in some particular cases (e.g., Potts model or tree structure). Local search methods are stochastic algorithms like Gibbs sampling, Guided Local Search (Park, 2002;Hutter et al., 2005), and Stochastic Greedy Search [START_REF] Mengshoel | [END_REF]. Some of them have theoretical asymptotic proof of convergence, i.e., the optimal solution is guaranteed to be found if infinite time is available. In practice, they may exhibit a better anytime behavior than exact methods on large and difficult problems (Marinescu et al., 2003;Hutter et al., 2005;[START_REF] Mengshoel | [END_REF], i.e., they produce better solutions in less time.

A few attempts have been done to combine exact and local search methods. A simple way is to run sequentially a local search algorithm then tree search, where solutions found by local search will be used as initial upper bounds for branch and bound exact methods. Another approach is to design a local search framework where the neighborhood exploration is performed by tree search in a systematic or non-systematic way. VNS/LDS+CP [START_REF] Loudni | [END_REF] combines a metaheuristic, Variable Neighborhood Search (VNS) [START_REF] Mladenović | [END_REF] with Limited Discrepancy Search (LDS) (Harvey and Ginsberg, 1995), a partial tree search method (Section 2). We propose in Section 3 an iterative variant of VNS/LDS+CP called UDGVNS, adapted to graphical models and able to prove optimality when the neighborhood size is maximal and with unbounded tree search. In Section 4, we describe a coarse-grained parallel version called UPDGVNS with asynchronous cooperative execution of UDGVNS processes with centralized information exchange as in (Davidovic and Crainic, 2012;[START_REF] Ouali | [END_REF]. We report experimental results in Sections 5-6.

PRELIMINARIES

2.1 GRAPHICAL MODEL Definition 1 A probabilistic graphical model (or Gibbs distribution) (Koller and Friedman, 2009) is a triplet (X , D, F) with X = {X 1 , . . . , X n }, a set of n random variables, D = {D 1 , . . . , D n }, a set of finite domains of values of maximum size d = max n i=1 |D i |, and F, a set of potential functions. Each variable X i takes values in D i . An assignment of X is a set x = (x 1 , . . . , x n ), with x i ∈ D i . The set of all possible assignments of X is denoted ∆ = n i=1 D i . Let A = {D 1 , . . . , D n } with D i ⊆ D i represent a restricted set of ∆ called a partial assignment. If S is a subset of V = {1, . . . , n}, X S , x S and ∆ S are respectively the subset of random variables {X i , i ∈ S}, the assignment (x i , i ∈ S) obtained from x, and the set of all possible assignments of X S . Given a set S of part of V , the set F = {f S } S∈S of maps from ∆ S to R + , called potential functions, is said to factorize a joint probability distribution P iff:

P(x) = 1 Z f S ∈F fS(xS) (1)
where Z = x∈∆ f S ∈F f S (x S ) is the normalizing constant, also called partition function.

Among the various tasks, the Most Probable Explanation (MPE) problem is to find the most likely assignment x ∈ ∆ to all the variables in X maximizing P(x). By taking the opposite of the logarithm of P(x), i.e., log P(x) = f S ∈Flog f S (x S )log Z, we obtain an additive model with ϕ(x S ) =log f S (x S ) called energy functions. Finding a solution of minimum energy is equivalent to MPE. In the rest of the paper, we consider energy minimization. When ϕ(x S ) maps to N + ∪ {∞}, the corresponding deterministic graphical model is called a Cost Function Network (CFN) (Meseguer et al., 2006). Finding a solution of minimum cost is the same as doing energy minimization on the equivalent probabilistic model [START_REF] Hurley | [END_REF].

TREE DECOMPOSITION

Definition 2 A tree decomposition of a connected graphical model is a pair (C T , T ) where T = (I, A) is a tree with nodes set I and edges set A and

C T = {C i | i ∈ I} is a family of subsets of X , called clusters, such that: (i) ∪ i∈I C i = X , (ii) ∀ f S ∈ F, ∃ C i ∈ C T s.t. S ⊆ C i , (iii) ∀ i, j, k ∈ I, if j is on the path from i to k in T , then C i ∩ C k ⊆ C j .
Definition 3 The intersection of two clusters C i and C j is called a separator, and noted sep(C i , C j ).

Definition 4 A graph of clusters for a tree decomposition (C T , T ) is an undirected graph G = (C T , E) that has a vertex for each cluster C i ∈ C T , and there is an edge

(C i , C j ) ∈ E when sep(C i , C j ) = ∅.
As finding an optimal tree decomposition is NP-hard, we use fast approximate algorithms like min-fill heuristic.

DFBB & LIMITED DISCREPANCY SEARCH

Depth-First Branch and Bound (DFBB) methods explore a search tree in a systematic way by recursively choosing the next unassigned variable to assign and by choosing a value in its domain for the assignment (the branch part) until a better solution is found or it can be proved that the subtree rooted at the current search node has no better solutions and it can be pruned (the bound part). DFBB depends on its variable and value ordering heuristics for branching in order to find good solutions rapidly and to reduce the size of the search tree to be explored.

Limited Discrepancy Search (LDS) (Harvey and Ginsberg, 1995) is a heuristic method that explores the search tree in a non-systematic way by making a limited number of wrong decisions w.r.t. its value ordering heuristic. We assume a binary search tree where at each search node either the selected variable is assigned to its chosen preferred value (left branch) or the value is removed from the domain (right branch). Each value removal corresponds to a wrong decision made by the search, it is called a discrepancy. The number of discrepancies is limited by a parameter denoted . See Algorithm 1 where lb(A) gives a lower bound on the minimum energy min x∈ D i ∈A Dilog P(x) of the partial assignment A.

In order to produce good quality solutions as time passes, a simple strategy is to iterate LDS with an increasing number of discrepancies going from min to max . Each iteration does at most discrepancies along the path from the root search node to a terminal node. Let h be the maximum height of the explored search tree, the number of terminals with discrepancies is bounded by ( h ). Thus, LDS has a time complexity in O(h +1 ) (for < h/2 [START_REF] Larrosa | [END_REF]) and a linear space complexity (thanks to the depth-first search principle as in DFBB). Notice that h ≤ n + ( right branches followed by n left branches). Without any discrepancy limit, we have at most n(d -1) value removals (right branches)

Algorithm 1: Pseudo-code of mono LDS Procedure LDS( , A, ub : In/Out, x : In/Out, opt : In/Out) if (∃Di ∈ A, |Di| > 1) then Choose an unassigned variable Xi ∈ X , |Di| > 1 ; Choose a value xi ∈ Di ; A ← (A \ {Di}) ∪ {{xi}} ; // left branch if (lb(A ) < ub) then LDS( , A , ub, x, opt) ; if ( > 0) then A ← (A \ {Di}) ∪ {Di \ {xi}} ; // right b if (lb(A ) < ub) then LDS( -1, A , ub, x, opt) else opt ← false ; else ub ← lb(A), x ← A ; // new sol. found
to reach a terminal node, where d is the maximum domain size of all variables. In order to be able to prove optimality, we set max = n(d -1), i.e., the last iteration is complete. In practice, optimality may be proved with a much smaller1 , and we just have to check if the discrepancy limit was reached during the search (falsified condition at line 1 of Alg. 1) to detect an incomplete search (flag opt set to false at line 2).

VARIABLE NEIGHBORHOOD SEARCH

Variable Neighborhood Search (VNS) [START_REF] Mladenović | [END_REF]) is a metaheuristic that uses a finite set of pre-selected neighborhood structures N k , k = 1, 2, ..., k max to escape from local minima by systematically changing the neighborhood structure if the current one does not improve the current incumbent solution. VNS repeatedly performs three major steps. In the first one, called shaking, a solution x is randomly generated in the neighborhoods of x denoted N k (x). In the second one, a local search method is applied from x to obtain a local optimum x . In the third one, called neighborhood change, if x is better than x, x is replaced with x and k is set to 1; otherwise, k is increased by one.

VNS METHODS IN GRAPHICAL MODELS

The use of VNS scheme for solving deterministic graphical models started with VNS/LDS+CP [START_REF] Loudni | [END_REF]. This approach is related to LNS (Shaw, 1998), but it adjusts dynamically the neighborhood size when search seems to stagnate as in VNS. VNS/LDS+CP algorithm. Let N k be the neighborhood structure of size k; N k denotes the set of all subsets of k variables among X . Let x be the current solution. First, x is partially destroyed by un-assigning a subset of k variables and an exploration of its (large) neighborhood is performed until the solution is repaired with a new one. The core of VNS/LDS+CP is its reconstruction phase. It runs LDS with a fixed discrepancy to explore the neighborhood of the solution. One advantage of this choice is its exploration speed that improves the quality profile and allows a more balanced exploration of the search tree. Decomposition Guided VNS. More recently, Fontaine et al. [START_REF] Fontaine | [END_REF] investigated the incorporation of tree decomposition within VNS/LDS+CP in order to efficiently guide the exploration of large neighborhoods. They proposed DGVNS, a first local search approach that exploits the graph of clusters provided by a tree decomposition to build relevant neighborhood structures. For both DGVNS and VNS/LDS+CP, the reconstruction phase is performed using LDS with a fixed discrepancy , but neighborhood structures are managed in a different way. Instead of the neighborhood structures N k used by VNS/LDS+CP, DGVNS uses neighborhood structures N k,c , where k is the neighborhood size and C c is the cluster where the variables will be selected from. If (k > |C c |), we complete the set of candidate variables to be unassigned by adding clusters adjacent to C c . The neighborhood change in DGVNS is performed in the same way as in VNS. However, DGVNS considers successively all the clusters. This ensures a better diversification by covering a large number of different regions.

AN ITERATIVE DGVNS METHOD

We present UDGVNS (see Algorithm 2), unifying two complete and incomplete search methods. As done by iterative LDS, it restores the completeness of DGVNS by applying successive calls with an increasing discrepancy limit, controlling if tree search was partial using opt flag and the fact that the current neighborhood kept some variables assigned in A (test at line 9). In UDGVNS, optimality can be proven in two cases: (i) when the current neighborhood corresponds to the whole problem and the discrepancy value is greater than or equal to the maximum number of right branches or (ii) by examining the bounds at the root node (ub = lb(D) at lines 5 and 8). In this case, the search space is implicitly explored by the algorithm, therefore optimality is proven.

The initial solution of UDGVNS is obtained at line 4 by a modified version of LDS, denoted LDS r , that stops immediately after a first solution is found2 . Thus LDS r (∞, D, ub, x, opt) will return a solution with nonzero probability (i.e., satisfying the constraints).

UDGVNS has two main parameters for tuning its compromise between optimality proof and anytime behavior: the discrepancy limit ( ) and the local search neighbor-Algorithm 2: Unified DGVNS algorithm.

Function UDGVNS( min, max, + , kmin, kmax, + k , ub :

In/Out, x : In/Out) : boolean let (CT , T ) be a tree decomposition of (X , D, F) ; opt ← true ;

LDS r (∞, D, ub, x, opt) ; // initial solution

if (ub = lb(D)) then opt ← true; c ← 1 ; // current cluster index r ← 0 ;
// number of iterations ← min ;

// initial discrepancy limit while (¬opt ∧ ≤ max) do i ← 0 ; // nb. of failed neighbor.

k ← kmin ; // init. neighbo. size while (¬opt ∧ k ≤ kmax) do A ← getNeighborhood(x, Cc, k) ; ub ← ub, opt ← true ; LDS r ( , A, ub , x , opt) ; // nei. search if (ub = lb(D)) then opt ← true; else if (A = D) then opt ← false; if (ub < ub) then x ← x , ub ← ub ; // new best sol i ← 0, k ← kmin ; r ← 0, ← min ; else i ← i + 1 ; if (k < kmax) then k ← min(kmax, kmin + k i) ; else k ← ∞; c ← 1 + c mod |CT | ; // get next clus. r ← r + 1 ; if ( < max) then ← min( max, min + r) ; else ← ∞;
return opt ; hood size (k). Recall that inside DGVNS, the maximum tree height for LDS r is h ≤ k + ≤ k(d -1), so LDS r has time complexity in O((k min( , d)) +1 ). The complexity increases exponentially with and polynomially with k. As soon as a better solution is found by the current neighborhood search (line 7), we stop the search in order to reinitialize the two parameters to their minimum value, as it is faster to explore small neighborhoods (lines 11-12). For each parameter, and k, we have tested three updating rules: increase by one at each iteration (+ /k = +), multiply by two at each iteration (a

+ /k b = mult2(a, b) = a × 2 b )
, and apply a Luby strategy [START_REF] Luby | [END_REF]

(a + /k b = Luby(a, b) = a × luby(1 + b)) 3 .
Operator + k controls the compromise between intensification and diversification. The goal of the Luby strategy is to exponentially increase the number of small neighborhoods explored compared to the number of larger ones. Whereas classical VNS algorithms will stuck on large problems, trying to diversify the search by explor-Operator + controls the compromise between incomplete and complete search. Using a fast growing strategy emphasis completeness whereas a slow growth should favor anytime behavior. We noticed that it is worthwhile to cover all the variables by the union of the explored neighborhoods in order to not miss some important variables. We tested a forth strategy for k which consists in a slow increment (by +1) at the beginning

until k = max i∈I (|C i |) + |C T | -1 then it jumps directly to k = k max .
This ensures that k grows slowly until the largest cluster has been totally explored by at least one neighborhood search. Then, when k = k max = |X |, UDGVNS does a restart, looking for an improved starting solution, using LDS r applied on the whole problem. If it fails to find a better solution, a larger discrepancy can be selected and UDGVNS continues its intensification process starting with a small neighborhood size (line 6).

A PARALLEL VERSION OF UDGVNS

This section describes how UDGVNS can been parallelized. We denote the resulting algorithm UPDGVNS for Unified Parallel DGVNS. It relies on asynchronous cooperative execution of several UDGVNS (i.e. workers) processes with centralized information exchange. In the proposed parallelization one global best solution is shared among the worker processes. Initially, the master initiates the search by launching npr worker processes in parallel with the same initial solution (line 4 done in the master). Each worker process obtains from the master a copy of the current best solution x, the index c of the cluster to be processed and performs destroy and repair operations on its local copy. As soon as a new solution is found by a worker, it is sent to the master and a new cluster is requested (line 13 done in the master) to restart a new exploration starting from the best available overall solution. Like UDGVNS, each worker process controls locally how the discrepancy limit and the neighborhood size k evolve during successive explorations (except that ub comes from the master at line 10). The whole process stops as soon as a worker terminates (line 14) with an optimality proof or it reaches max and k max .

UAI EVALUATION RESULTS

Benchmarks description. We performed experiments on probabilistic and deterministic graphical models coming from a large multi-paradigm benchmark repository [START_REF] Hurley | [END_REF] 4 . Among the 3016 available instances, we selected all the instances that were used in previous UAI competitions, in image analysis, or in Cost Function Network optimization. It includes 319 instances from UAI 2008 Evaluation and Probabilistic Inference Challenge (PIC 2011) 5 , 1461 instances from the Computer Vision and Pattern Recognition (CVPR) OpenGM2 benchmark 6 [START_REF] Kappes | [END_REF], and 281 instances from the Cost Function Library 7 . In order to have a fair comparison between solvers, we preprocessed all the instances by polynomial-time problem reformulations and simplifications that remove variables (using bounded and functional variable elimination (Favier et al., 2011) , 2015). IBM ILOG CPLEX 12.7.0.0 (using parameters EPAGAP, EPGAP, and EPINT set to zero to avoid premature stop) was reported as being very competitive on some image analysis [START_REF] Kappes | [END_REF] and Markov Random Field problems [START_REF] Hurley | [END_REF]. CPLEX explores its search tree using best-first search. It applies several heuristics methods to find good solutions before and during the search. We also compared with message-passing algorithms: LIBDAI 12 (Mooij, 2010), winner of the UAI 2010 Evaluation at 20sec. and 1hour time limits, MPLP2 13 (Sontag et al., 2008;2012), and TRW-S 14 (Kolmogorov, 2006) Anytime upper bound profiles. In order to summarize the evolution of upper bounds as time passes, we took a subset of 114 difficult instances, unsolved in 1 hour by our DFBB, LDS, and VNS methods (whereas CPLEX could solve 17 instances). Specifically, for each instance I we normalize all energies as follows: the best, potentially suboptimal solution found by any algorithm is 1, the worst solution is 2. This normalization is invariant to translation and scaling. Figure 1b shows the upper bound behavior for different VNS strategies compared to DFBB and LDS (see Supp. material for a zoom). The ranking of best methods is the opposite of the cactus plot order, except for (k add1/jump, mult2) which comes in second position. The = 3 strategy got the best upper bounds in average, but still very close to the other VNS strategies, except may-be (k Luby, mult2) and (k mult2, mult2). We conclude that our new iterative UDGVNS method (especially (k add1/jump, mult2)) offers a good compromise between anytime behavior and optimality proof. These results also show that variable neighborhood search is by far superior to classical systematic (DFBB) or nonsystematic (LDS) tree search methods, improving by more than 20% (resp. 10%) the quality of the solu-tions. In the following figures, we assume UDGVNS with (k add1/jump, mult2) strategy.

Figure 1c compare UDGVNS with state-of-the-art methods. Message-passing algorithms like TRW-S and LIB-DAI gave the worst results. They could not find any solution for 20 (resp. 19) instances (mostly in UAI/Linkage and CFN/SPOT5 categories, both containing hard constraints). The same problem occurred for MPLP2 on 5 instances (SPOT5), but it obtained much better results on the remaining instances. CPLEX ran out of memory on two instances without producing any solution due to the heavy local polytope encoding (CFN/Warehouse/capb,capmq5). All other methods got better results in average and produced at least one solution per instance. According to its initial phase setting, DAOOPT provides different anytime behaviors, very close to the best solutions in 1 hour. UDGVNS performed the best, slightly better than INCOP+TOULBAR2, improving by 1.7% (resp. 2.3%) on average after 1 hour (20 min). Parallelization. Finally, in order to evaluate the impact of core numbers, we consider the anytime upper bound behavior of the parallel release: UPDGVNS using (k add1/jump, mult2) with min = 3, taken from the best strategies enlightened by UDGVNS. We made a comparison with CPLEX using 10 and 30 cores (results in Supp. material), improving its anytime behavior, but still being far from the other competitors (30 cores gave solutions 10% higher than UPDGVNS after 1 hour). We could not compare with the parallel version of DAOOPT as it is based on a different cluster engine (condor) and it does not parallelize its initial phase.

Figure 1d shows that UPDGVNS (with 10 or 30 cores) provides slightly better upper bounds than UDGVNS (using 1 core) in less than 20 min. The results seems to be, in average, poorly sensitive to the cores number, due to the fact than the 10-core curve is extremely close from the 30-core one. The 10 and 30-core UPDGVNS curves converge quickly in less than 2min. INCOP+TOULBAR2 quickly drops out around 1 min and never reaches the same quality level. UDGVNS converges slower but still going down after 20 min. DAOOPT (1200sec setting) gave results 10% in average worse than UPDGVNS with either 10 and 30 processors.

COMPUTATIONAL PROTEIN DESIGN

Benchmarks description. One main challenge in Computational Protein Design (CPD) lies in the exploration of the amino-acid sequence space, while considering, to some extent, side chain flexibility. The exorbitant size of the search space urges for the development of efficient exact deterministic search methods enabling identification of low-energy sequence-conformation models, corresponding to the Global Minimum Energy Conformation (GMEC). The GMEC corresponds to a maximum probability mass due to the Boltzmann relation between molecular energy and probability, which is, for that matter, equivalent in the Markov Random Field modeling with Maximum A Posteriori probability (MAP-MRF) estimation.

CPD has been recently modeled as CFN, and various intensive benchmarking had demonstrated that in most cases CFN techniques are more efficient than other discrete optimization methods (Allouche et al., 2014), especially for full design of instances lesser than 100 Amino acid (Simoncini et al., 2015) . In this paper we would like to evaluate VNS methods capability to solve difficult instances. Accordingly, we tried to generate new larger ones supposed to be more difficult to solve than those generated in (Simoncini et al., 2015). For this aims, in a first stage, we selected proteins structure in the PDB databases with size between 100 and 300 amino-acid.

On the basis of a database query described in the supplementary material, we obtained a first set including 438 PDB references, which has been sorted according to the average volume per variable calculated using the radius of gyration of each instance (see Supp. material). At last, we extracted for benchmarking only the 20th first elements and used as benchmarking set. Each protein structure was fully redesigned according to (Simoncini et al., 2015) protocol. On the basis of energy matrix generated with a modified release of PYROSETTA.4 script (Simoncini et al., 2015) in order to use the last released Rosetta force-field (aka Beta November 2016) [START_REF] Alford | [END_REF]. The instances (see Supp. material) contain from 130 up to n = 282 variables with maximum domain size from 383 to 438, and between 1706 and 6208 cost/energy functions, and besides the tree width ranges from 21 to 68 and from 0.16 to 0.34 for a normalized tree width.

In this section, we evaluate the effectiveness of UDGVNS and UPDGVNS on CPD instances in terms of optimality proof (cf. part (a)) and solution quality vs. CPU time (cf. part (b)). All computations were performed on a cluster of 96-core AMD Opteron 6174 at 2.2 GHz and 256 GB of RAM. The parallelization has been done in MPI.

Experimental Protocol. We compared UDGVNS and UPDGVNS with VNS/LDS+CP and FIXBB (CPD dedicated algorithm provided by Rosetta package). We also compared with TOULBAR2. We tested TOULBAR2 using Virtual Arc Consistency [START_REF] Cooper | [END_REF] as preprocessing combined with hybrid best-first search with tree decomposition (BTD-HBFS) using min-fill heuristic. Dead end elimination is turned off, according to (Simoncini et al., 2015) (the resulting command line is: toulbar2-B=1-A-dee:-O=-3). The TOULBAR2 experiments use a 24-hour CPU time limit.

Concerning the VNS methods, in order to evaluate variability due to the random selection of neighborhoods, a set of 10 runs per instance with different seeds has been performed with a time limit of 1-hour per run.

For the parallel strategy, the number of processes npr is set to 96 (i.e. maximum number of available processes). For UDGVNS and UPDGVNS, following the results observed in Section 5, we considered the following two settings for operators + k and + : (k add1, = 3) which yields the best anytime performances and (k add1/jump, mult2) which offers the best compromise between both anytime performances and optimality proof. k min , k max and min have been respectively set to 4, n (the total number of variables) and 3 (they correspond to the same parameter settings as those described in [START_REF] Fontaine | [END_REF]).

(a) Optimality results. Our first set of experiments aims at evaluating the efficiency of UDGVNS and UPDGVNS in terms of optimality proof comparing with TOUL-BAR2.

As indicated before, we used the setting (k add1/jump, mult2). For all methods (see Supp. material), the CPU-time is the time to find and prove an optimal solution. TOULBAR2 clearly outperforms the two VNS methods: UDGVNS and UPDGVNS were able to prove optimality only on 3 instances among the 10 ones closed by TOULBAR2. However, the VNS methods are faster (except on 5dbl), despite the fact that BTD-HBFS benefits from the lower bounds reported by HBFS in individual clusters to improve the anytime behavior and from the global pruning lower bounds of BTD. This greatly improves the overall performance of TOULBAR2 compared to VNS methods. (b) CPU time and solution quality results. Our second set of experiments aims at evaluating VNS capability with respect to finding the optimal solution or a solution of better quality on instances for which optimal solutions are unavailable. For this aims, we selected (k add1, = 3) as setting for operators + k and + .

Table 1 shows a comparative evaluation of VNS methods with FIXBB and TOULBAR2. For each instance and each VNS method, we report the number of successful runs to reach the optimum (or the best found solution for unsolved instances) within a 3600-second time limit, the average CPU times (in seconds) over the 10 runs (for unsuccessful runs, the CPU time is set to the time limit) ± the standard deviation. The Energy gap ∆E between the best VNS solution and the two external references FIXBB and TOULBAR2 are also given. For TOULBAR2, reported CPU times correspond to times to find an optimal solution (for solved instances) or a best one (for unsolved ones) within the 24-hour time limit.

VNS methods vs. FIXBB. Rosetta Modeling suite is one of the most popular software package used in CPD field. It provides a Monte Carlo based Simulated Annealing algorithm called FIXBB. In this work, the best solutions exhibited by 1000 FIXBB cycles performed on each CPD instance have been used as base-line to compare solution quality of the solution provided by VNS methods when TOULBAR2 BTD-HBFS fails to solve the instance.

The FIXBB CPU times are two orders of magnitude higher than the 1h time limit imposed for VNS evaluation. They are not reported as they exceed 100-hours in sequential mode, even if FIXBB cycles are independent and thus be easy to parallelize. As we can see in Table 1 (see column ∆E of ( 4)), the solution quality of FIXBB is in all cases inferior to the best solution found by the VNS methods. The Energetic gap ∆E between FIXBB solution and the best VNS overall solution ranges between +0.16 and +5.20 Rosetta Energy Unit (R.E.U). As shown in (Simoncini et al., 2015) such a level of energy difference can strongly impact the designed protein solution (i.e. the corresponding sequences of the two methods can be far in terms of hamming distance).

VNS methods vs. TOULBAR2. The comparison between best solutions found by VNS methods and TOULBAR2 shows that, excepted VNS/LDS+CP, both UDGVNS and UPDGVNS provide in all case the same or even better solution than TOULBAR2 (see column ∆E of (5) in Table 1). The Energetic gap ∆E in the worst case reaches 17.86 R.U.E. Concerning the number of successful runs reported over the 10 runs, VNS/LDS+CP seems less robust respectively than UDGVNS and the parallel release UPDGVNS. This last one provides in all cases the best solution over all. Table 1 also compares the VNS methods in terms of speedups. We observe that speedup values are fluctuating from one instance to another, very likely due to the tree decomposition resulting from the 3D shape of each instance. For VNS/LDS+CP and UDGVNS, it range between 0.62 and 2.86 over the 14 instances solved by both methods. As expected, when tree decomposition and parallelization are used, not only the speed of resolution increases but the reliability too (speedup values between 4.1 and 18.5). Moreover, the comparison between UDGVNS and UPDGVNS shows significant accelerations (between 2.06 and 9.22), thus confirming the practical interest of parallelization in addition to the exploitation of problem decomposition. We can observe (in Supp. material) that both UDGVNS and UPDGVNS find optimal solutions more quickly than TOULBAR2 within the one hour time limit (with the exception of UDGVNS

CONCLUSION

In this paper we proposed a unified view of VNS methods including various LDS and neighborhood evolution strategies. Experiments performed on difficult instances, coming from a large graphical model benchmark, showed that our hybrid method has a much better anytime behavior than existing tree search methods and still being competitive for proving optimality. We further proposed a parallel version of our method improving its anytime behavior. It remains as future work to automatically choose the best parameter settings per instance.

Supplementary Material

July 11, 2017

1 Cactus plot of UDGVNS versus cplex, daoopt, and incop+toulbar2 We report daoopt time from (OD17), obtained on a cluster of dual 2.67 GHz Intel Xeon X5650 6-core CPUs and 24 GB of RAM.

CPD instances filtering

In this paper we would like to evaluate VNS methods capability to solve difficult instances. Accordingly, we tried to generate new larger ones supposed to be more difficult to solve than those generated in (SAdG + 15). For this aims, in a first stage, we selected protein structures in the PDB databases (http://www.rcsb.org/pdb/home/home.do) with sizes range between 100 and 300 amino-acid. The resulting query has been filtered with the following criteria: resolution has to be lower than 2.5 A, membrane proteins, proteins complex, as well as proteins with disulfur bright have been removed, in addition, with proteins including non natural amino acid. We also discarded proteins with missing residues, out of the N and C terminal part of the sequence in order to select a protein subset without any hole. Proteins with identity sequences higher than 90% are unselected. The corresponding remaining set includes 438 PDB references. Therefore, we sorted the 438 putative instances according to the average volume per variable and we extracted for benchmarking only the 20th first elements and used as benchmarking set. Each protein structure was fully redesigned according to (SAdG + 15) protocol. On the basis of energy matrix generated with a modified release of pyrosetta.4 script (SAdG + 15) in order to use the last released Rosetta force-field (aka Beta November 2016) (A + 17). The instances characteristic contain from 130 up to n = 282 variables width maximum domain size from 383 to 438, and between 1706 and 6208 cost functions. The tree width ranges from 21 to 68, and from 0.16 to 0.34 for a normalized tree width.

In order to select well decomposable instances, after preliminary reorientation of each protein according to his inertial moment, we filtered the resulting PDB set by 3D geometrical criteria. By construction, due to the cutoff distance used for pairwise energy calculation and his related constraints in the model, Globular protein (.i.e spherical one) will correspond to CFN very nearby as click. In practice, the spherical shape which can be detected by a symmetric repartition of the Radius of gyration (Rg) component Rg(x), Rg(y), Rg(z). Indeed, Rg and his subcomponent are defined as the root mean square distance from each atom of the protein to their centroid or in the orthogonal plan to the related axis (x,y,z). Thus, for selecting non spherical protein , we first calculated the Rg i respectively around x , y and z axis as and we filtered the proteins characterized by a ratio:

min(Rg(x)/Rg + Rg(y)/Rg + Rg(z)/Rg) (1)
The min(Rg(i)/Rg) is one figure to detect putative well structured instances.

Another structural criterion that can be used for structured instances filtering is a decreasing sort according to the approximated means space volume occupied associated to each variable. This volume is also related to the Gyration radius by the following equation:

V = 4 3 * π * Rg 3 |X| (2)
This criterion does not give at first sight a direct information about fold shape, but interestingly, both structural criteria produce the same 15th first instances set, with only few re-ranking, very likely because protein compactness got maximum bound and indeed a minimum average volume per amino-acid when the protein fold is close to sphere or vice versa for linear protein fold. After the 15th instances, the two critters remains less correlated. As a matter of fact, in the 5 next instances, the covering between the two sorted list is 2 over 5. On Account of computer power for benchmarking, we can not bench the full protein set. Table 3: Tables (3a) and (3b) report, for each instance, the CPU times spent by VNS methods (within the 3600-seconds time limit) to obtain the best solution computed by toulbar2 (5). A '-' indicates that the corresponding solver was not able to compute a solution of equal/better quality than toulbar2.

  (d) Comparing the anytime behavior of parallel UPDGVNS.

Figure 1 :

 1 Figure 1: Comparative evaluation on probabilistic and deterministic graphical models

Figure 1 :Figure 2 :

 12 Figure 1: Number of instances solved by each method as time passes (UDGVNS = UDGVNS(k add1/jump, mult2)).

  (1) : VNS(k add1, = 3) (2): UDGVNS(k add1, = 3) (3): UPDGVNS(npr, k add1, = 3). Table (3c) reports the CPU-times required by UDGVNS, UPDGVNS and toulbar2 to prove the optimum within the 24-hours time limit. A '-' indicates that the corresponding solver failed to prove optimality. (2 * ): UDGVNS(k add1/jump, mult2)(3 * ): UPDGVNS(npr, k add1/jump, mult2).

  ), values (using dead-end elimination (de Givry et al., 2013)), and fixed-value potentials, after an ini-Experimental protocol. LDS r employs a randomized (for breaking ties) dynamic variable ordering heuristic 9 . Its value ordering heuristic chooses the EAC value as the preferred value and lower bounds are deduced by enforcing EDAC, as explained in(Larrosa et al., 2005). In the following, we setk min = 4, k max = n = |X |, l min = 1, and l max = n(d -1). Here, DFBB corresponds to UDGVNS(∞, ∞, +, ∞, ∞, +), LDS to UDGVNS(1, ∞, mult2, |X |, |X |, +), and DGVNS to UDGVNS(3, 3, +, k min , k max , +).We compared with state-of-the-art exact solvers. DAOOPT 10 won PIC 2011. It has a time-bounded initial phase of lower bound computation based on Message Passing Linear Programming algorithm(Sontag et al., 2008; 2012;[START_REF] Ihler | [END_REF] and mini-bucket elimination[START_REF] Dechter | [END_REF] with iterative min-fill heuristic. It also finds good initial upper bounds using stochastic local search GLS +(Hutter et al., 2005). We tested three parameter settings as suggested in(Otten et al., 2012), controlling the time spent to compute initial lower and upper bounds. In the 3600sec setting, SLS is run 20 times with 10 seconds per run. The best solution found is used as an initial upper bound for an AND/OR exhaustive tree search. INCOP+TOULBAR2 11[START_REF] Hurley | [END_REF] won the UAI 2010 Evaluation at 20 min. time limit.

	from 4 (CVPR-GeomSurf-3-gm13) to 48,566 (CVPR-
	ColorSeg-8-crops-small) with mean value n ≈ 403.4 (instead of 1, 316.5 before preprocessing). For DGVNS
	methods, we built tree decompositions using min-fill
	heuristic. Because the number of clusters m = |C T | can be very large (m ≈ 256.7), we merged any pair of connected clusters (C i , C j ) when the separator size
	is too large compared to the individual cluster sizes
	(|sep(C i , C j )| > 0.7 min(|C i |, |C j |)), resulting in m ≈ 19.9 and mean tree width max i∈I (|C i |) ≈ 76.9 (instead of 49.6 without merging). In order to experiment sequen-
	tial and parallel methods on the most difficult instances,

4 genoweb.toulouse.inra.fr/ ˜degivry/evalgm 5 graphmod.ics.uci.edu/uai08/Evaluation/Report/ Benchmarks, www.cs.huji.ac.il/project/PASCAL 6 hci.iwr.uni-heidelberg.de/opengm2 7 costfunction.org/benchmark tial lower bound computation by Equivalence Preserving Transformations (Cooper et al., 2010) (enforcing Virtual Arc Consistency (VAC) as a message-passing algorithm). The resulting instances are smaller while preserving the same optimum. We used TOULBAR2 with options -A -z=2 for this preprocessing step. We kept 1669 non-trivial instances (with more than one variable) for the experiments. The number of variables n ranges we selected a subset of instances unsolved in 1 hour by all our DFBB, LDS, and VNS algorithms. We took at most twenty instances per problem category (avoiding over-representation issues by some categories), resulting in a selection of 114 difficult instances 8 . TOULBAR2 takes a starting solution from the best result of three runs of the IDWalk (Neveu et al., 2004) local search algorithm (100,000 local moves per run). It is followed by an exhaustive hybrid best-first search (Allouche et al.

  . Note that LIBDAI and TRW-S are applied on the original instances rather than the preprocessed ones as we found they produced better results without applying VAC first. All solvers read problems in uai tabular format, except CPLEX which uses the local polytope formulation (called support encoding in[START_REF] Hurley | [END_REF]). All computations were per-

	formed on a cluster of 48-core AMD Opteron 6176 at
	2.3 GHz and 384 GB of RAM with a 1-hour CPU time
	limit 15 .
	Optimality results. The efficiency of DFBB, LDS, and
	VNS methods to prove optimality is shown in the cac-
	tus plot of Figure 1a. DFBB was slightly more effi-
	cient than LDS and solved 1442 (resp. 1433) instances
	among 1669 in 1-hour time limit. They are followed
	by three UDGVNS strategies with (k mult2, mult2)
	(1430 solved), (k Luby, mult2) (1425 solved) and
	(k add1/jump, mult2) (1421 solved), remaining
	very close in terms of performance. Next, a set
	of three less-and-less efficient UDGVNS strategies rise:
	(k add1, mult2) (1384 solved), (k add1, Luby)
	(1361 solved) and (k add1, add1) (1333 solved),
	showing the importance of faster discrepancy increase

to speed up optimality proofs. The worst strategy here was using a fixed discrepancy level ( = 3 as originally proposed in

[START_REF] Fontaine | [END_REF]

) which solved 1128 instances. For comparison, DAOOPT (3600sec setting) solved 1409 instances, CPLEX solved 1423, and IN-COP+TOULBAR2 1440 instances (results in Supp. material).

Table 1 :

 1 Comparative evaluation on CPD instances. In bold instances solved by TOULBAR2. ncl is the number of clusters. ∆E is in Rosetta Energy unit. It has been obtained by the difference of solutions costs divided by the cost shift used during modeling, in this case 10 8 . on 5dbl and 5e10). For UDGVNS, speedup values range from 0.81 to 18.3 with a mean of 4.23 over all the solved instances. For UPDGVNS, the ratio in terms of speedup is greatly amplified (between 5.25 and 56.74 with a mean of 18.48 over all the solved instances).

	Instance	ncl	Succ.	(1)	Time (s)	Succ.	(2)	Time (s)	Succ.	(3) Time (s)	(4) ∆E	(5) Time (s)	∆E	(1/2)	Speed-up (1/3)	(2/3)
	5dbl	87 10/10	2,761± 67	10/10		963± 53	10/10	149± 18	0.27	783	0	2.86 18.53	6.46
	5jdd 3r8q	168 157	0/10 0/10		TO TO	10/10 10/10	3,248 ± 162 3,478± 82	10/10 10/10	646± 72 397± 47	4.08 4.19	20,662 12,762	0.04 0	--	--	5.02 8.76
	4bxp	108 10/10	1,213± 32	10/10	1,352± 40	10/10	216± 24	0.26	2,966	0	0.89	5.61	6.26
	1f00	177	0/10		TO	10/10	2,575± 29	10/10	542± 30	4.38	9,749	0	-		4.75
	2x8x	131 10/10 3,312± 111	4/10 2,801± 1,071	10/10	521± 58	4.16	69,213	3.61	1.18	6.35	5.37
	1xaw	66	0/10		TO	10/10		581± 38	10/10	260± 57	2.73	2,804	0	-	-	2.23
	5e10	74 10/10	1,667± 86	10/10	1,412± 21	10/10	132± 24	0.26	1,171	0	1.18 12.62	10.7
	1dvo	82 10/10		940± 22	10/10		940± 22	10/10	197± 25	2.90	34,142	0.18	1	4.77	4.77
	1ytq	67 10/10 2,235± 211	10/10	1,304± 40	10/10	280± 26	1.67	17,063	0.31	1.71	7.98	4.65
	2af5	140	0/10		TO	10/10	2,659± 75	10/10 894± 247	4.37	86,029	0.60	-	-	2.97
	1ng2	86 10/10	1,065± 69	10/10		547± 18	10/10	260± 26	1.14	38,730	5.93	1.94	4.09	2.10
	3sz7	79	0/10		TO	10/10	2,952± 276	10/10	320± 50	3.11	82,625	0.54	-	-	9.22
	2gee	110 10/10	1,647± 8	10/10	1,276± 19	10/10	286± 38	1.68	5,021	0	1.29	5.75	4.46
	5e0z	73 10/10		621± 15	10/10		995± 18	10/10	105± 8	0.16	999	0	0.62	5.91	9.47
	1yz7	87 10/10	2,148± 9	10/10		945± 24	10/10	457± 85	2.91	83,816	3.20	2.27	4.70	2.06
	3lf9	72 10/10	1,636± 31	10/10	894± 41.17	10/10	216± 21	2.41	2,667	0	1.82	7.57	4.13
	3e3v	91	0/10		TO	10/10	2,327± 621	10/10	263± 38	2.57	81,574	0.15	-	-	8.84
	1is1	107 10/10	3,178± 51	10/10	2,274± 421	10/10	329± 48	3.45	63,832	0.42	1.39	9.65	6.91
	5eqz	89 10/10	1,849± 49	10/10		697± 6	10/10	225± 16	2.20	12,768	0	2.65	8.21	3.097
	4uos	118 10/10 2,304± 988	10/10	2,109± 62	10/10	469± 79	5.20	58,589 17.86	1.09	4.91	4.5
	TO: TimeOut (1): VNS/LDS+CP(k add1, = 3) (2): UDGVNS(k add1, = 3) (3): UPDGVNS(npr, k add1, = 3) (4): FIXBB	(5): TOULBAR2
						For unsolved in-
	stances, we measured the CPU times spent by the three
	VNS methods (within the 3600-second time limit) to
	reach the best solution produced by TOULBAR2. Re-
	sults (in Supplementary material) show a clear order-
	ing in terms of CPU times across different solvers, from
	TOULBAR2, VNS/LDS+CP, UDGVNS, and UPDGVNS.
	The speedup values are significantly improved, in par-
	ticular with UPDGVNS (between 32.3 and 354.6 with a
	mean of 168.7 over all the unsolved instances). These
	results confirm the superiority of VNS methods in terms
	of anytime performance as compared to TOULBAR2.

3

  Anytime upper bound of UDGVNS and UPDGVNS versus cplex

		1.6		UPDGVNS (30 cores)
				UPDGVNS (10 cores)
		1.5		UDGVNS (1 core) cplex (30 cores)
	Normalized upper bounds	1.2 1.3 1.4		cplex (10 cores) cplex (1 core)
		1.1
		1
		0	20min	40min	1hour
			Wall-clock real time
	Figure 3: Anytime upper bound with 1 , 10 and 30 processors respectively for cplex and UPDGVNS
	(UPDGVNS = UPDGVNS(k add1/jump, mult2)).
	Linkage (optimum / worst solution) (4625/21439) (5258/166553) (6651/104904) (6406/629929) pedigree19 pedigree31 pedigree44 pedigree51 cplex (1 core) 790 59.3 6.35 36.23 cplex (10 cores) 191 9.00 2.48 9.43 cplex (30 cores) 75 7.17 2.69 5.34 daoopt (1 core) 375,110 16,238 95,830 101,788 daoopt (20 cores) 27,281 1,055 6,739 6,406 daoopt (100 cores) 7,492 201 1,799 1,578 udgvns (1 core) -(4949) -(5258) -(6722) -(6406) updgvns (10 cores) -(4762) 3,341 -(6651) -(6406) updgvns (30 cores) -(4626) 1,775 -(6651) -(6406)

4 Solving time and anytime behavior of U(P)DGVNS versus cplex and daoopt on UAI-Linkage category

Table 1 :

 1 Best CPU time (in seconds) for sequential versions and best wall-clock time for multiplecore ones to find and prove optimality on Pedigree instances. A "-" means no proof of optimality in 1 hour (except daoopt with no time limit) (in parenthesis, unnormalized upper bound founds after 1 hour).

Table 2 :

 2 Characteristics of PDB instances: pdbid is the code reference in PDB database, |X | is the number of variable, d is the maximum domain size, e is number of cost functions, tw is the Min-fill tree width and tw/|X | a normalized tree width by |X |. The two last columns correspond to structural criteria respectively defined in (1) and (2). Grey color corresponds to the re-ranking when one min(Rg x /Rg) or V is used for PDB list sorting.6 Comparing best solutions of toulbar2 vs vns methods

	pdbid |X | 5dbl 130 384 1,706 21 d e tw tw/|X | min(Rgx/Rg) V ( 3 /var) 0.16 0.150 1,212.49 5jdd 263 406 5,220 41 0.16 0.239 655.58 3r8q 271 418 5,518 43 0.16 0.341 472.88 4bxp 170 439 2,636 33 0.19 0.316 457.81 1f00 282 430 6,208 51 0.18 0.269 439.28 2x8x 235 407 4,745 44 0.19 0.354 404.42 1xaw 107 412 1,623 28 0.26 0.308 378.04 5e10 133 400 2,286 34 0.26 0.294 344.68 1dvo 152 389 2,587 51 0.34 0.420 343.38 1ytq 181 415 3,449 54 0.30 0.392 332.69 2af5 292 410 5,693 68 0.23 0.427 330.23 1ng2 176 397 3,135 60 0.34 0.473 309.16 3sz7 151 450 2805 49 0.32 0.403 304.87 2gee 188 397 3,715 38 0.20 0.367 293.25 5e0z 136 420 2,367 36 0.26 0.362 279.00 1yz7 176 418 3,538 49 0.28 0.414 276.35 3e3v 154 436 2,976 37 0.26 0.367 251.97 3lf9 120 416 2,133 31 0.24 0.323 251.51 1is1 185 431 3,740 48 0.26 0.459 245.58 5eqz 138 434 2,567 33 0.24 0.338 241.93 4uos 188 383 4,161 44 0.23 0.347 234.11

It was less than or equal to 128 for all instances completely solved by LDS and VNS methods in Section 5.

It stops the recursive LDS procedure at line

of Alg. 1 and set the optimality flag opt to false.

Recall luby(i) = {1, 1, 2, 1, 1, 2,

4, 1, 1, 2, 1, 1, 2, 4, 8, . . .}.ing larger neighborhoods, VNS using Luby will spend more time on small neighborhoods in order to locally improve the current solution, favoring intensification. The mult2 strategy reduces the number of neighborhood explorations at a given discrepancy limit, in order to try larger discrepancy limits more rapidly. If the problem is solvable by a complete search within the time limit, it will also speed-up the optimality proof.

UAI DBN, Grid, Linkage, ObjectDetection, CVPR ChineseChars,

Using parameter -pe parallel smp min(2x, 30) on a SUN Grid Engine for a method exploiting x core(s) to ensure half-load of the nodes on the cluster.

Acknowledgements This work has been partially funded by the french "Agence nationale de la Recherche", reference ANR-16-C40-0028. We are grateful to the genotoul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul) and the HPC HAYTHAM of university of Oran 1 for providing computing and storage resources. We thank Mathieu Fontaine for his contribution to the code of DGVNS.