**Proof.** We will use mathematical induction on the length of the extension, $m$, on $S$ to prove the theorem. Let $A$ denote the VA clauses of $M_1$, and $A_0$ denote the VA clauses of $B_1$.

Base case is trivial as $B_0$ and $B_1$ are given as isomorphic when $m = 0$.

Induction step: Suppose the search path $S$ is extended $m$ times, where $m > 1$, so that $S_m$ is LNH-compliant and isomorphic to a subset $A_m \subseteq A$. Then by Lemma 8, $S_m$ can be extended by one VA clause with the cell term $t_{m+1}$, chosen by the said cell selection strategy, to $S_{m+1}$ which is LNH-compliant and isomorphic to $A_{m+1} \subseteq A$.

Note that a model finder may do propagations after a cell value assignment. That is, some cell terms can be assigned values inferred from existing VA clauses. Propagations can be viewed as part of the cell selection strategy and be handled the same way as regular cell value assignments.

We can therefore conclude by mathematical induction that $S$ can be extended to a complete search path when all cell terms in $\mathcal{F}$ are filled with values such that $S$ represents the model $M_0$, is LNH-compliant, and is isomorphic to $A_s \subseteq A$. Since $M_0$ and $M_1$ are of the same size, so $A_s$ and $A$ must necessarily be of the same size and hence must be equal. Therefore, $M_0$ is isomorphic to $M_1$.

Theorem 9 shows that isomorphic cubes always extend to isomorphic models. So, one of the isomorphic cubes may be discarded without losing any non-isomorphic model.

**Corollary 10.** On searching under the LNH with any cell selection strategy on a signature $\Sigma$ and an FOL formula $\mathcal{F}$ on $\Sigma$, if $M_1$ is a model in $\mathcal{F}$, then there is a complete search path $S$ under the said LNH that results in a model $M_0$ which is isomorphic to $M_1$.

Corollary 10 proves the completeness of the LNH in that every model in any search is isomorphic to some model found by searching under the LNH. An alternative proof of the corollary is given in [50].

# 4 Searching with Cubes

Cubes can be constructed to partition the search space into non-overlapping subtrees that can be processed in parallel. It is not necessary to search all the subtrees that originate from the collection of cubes that span the entire search space because isomorphic cubes in the same collection can be eliminated without losing non-isomorphic models. For example, suppose we want to search for models of order 3 or more on a function $f: D^2 \to D$ under the LNH with a cell selection strategy that selects $f(0,0)$ then $f(1,1)$ as the first 2 cell terms in the search process. There are at most 6 cubes of length 2 (listed below) under the said LNH and cell selection strategy, so together they must span the whole search space in the sense that every search path that starts with the cell terms $f(0,0)$ then $f(1,1)$ in the search tree must include one of the 6 cubes in it.

1. $\langle f(0,0) = 0; f(1,1) = 0 \rangle$.

2. $\langle f(0,0) = 0; f(1,1) = 1 \rangle$.

3. $\langle f(0,0) = 0; f(1,1) = 2 \rangle$.

4. $\langle f(0,0) = 1; f(1,1) = 0 \rangle$.

5. $\langle f(0,0) = 1; f(1,1) = 1 \rangle$.

6. $\langle f(0,0) = 1; f(1,1) = 2 \rangle$.

Since $\pi_{(0,1)}(\text{Cube } 1) = \{f(1,1) = 1, f(0,0) = 1\} = \pi_{id}(\text{Cube } 5)$, so Cubes 1 and 5 are isomorphic and one of them can thus be discarded without losing non-isomorphic models per Theorem 9. This example demonstrates the importance of keeping the LNH in the search –
<figure><img src="image_66.png" /></figure>

Note: $t_0$ denotes $f(0,0)$ and $t_1$ denotes $f(1,1)$. A dotted line with a cross is a branch pruned by the LNH, except for the branch ending on the VA clause $t_1 = 1$ (the shaded node), which is pruned by the isomorphic cubes removal algorithm.

<caption>Figure 2 Partial Search Tree Showing Cubes of Length 2.</caption>

it cuts the search space from potentially $n^2$ cubes down to 6. Theorem 9 allows us to further cut the number of cubes down to 5 (see Figure 2 for illustration). More isomorphic cubes can be removed with longer cubes (see Table 2).

The procedure of removing isomorphic cubes starts with generating a set of short cubes (typically of length 2 for a binary operation) that spans the entire search space. The model finder takes short cubes as inputs and runs with them as if they are generated by itself to generate longer cubes of predefined length $l$. Specifically, the model finder runs as usual, except that it emits the cubes of length $l$ when the depth of the search tree reaches $l$. After outputting the cube, the model finder backtracks as if it has reached the bottom of the search tree, and runs on a new branch as usual until all cubes of length $l$ are generated. Some models may be generated in this process due to propagation, and they are kept as part of the final outputs. Next, the cubes are compared for isomorphism and only one of any pair of isomorphic cubes is kept. This new set of non-isomorphic cubes of length $l$ will be used as inputs to the model finder in the next round of generation of longer cubes. The process is repeated until the desired length of cubes is reached.

For searching models defined by one operation of arity $k$, we use the sequence of lengths $l: k, 2^k, 3^k, 4^k, \dots$. This is to match the concentric cell selection strategy (see Example 4 for its definition) of the finite model finder such as Mace4. We will discuss the best cube length to use in Section 5.3.

Finally, non-isomorphic cubes of the target length can then be processed independently in parallel and their output models collected separately.

## 4.1 Invariants

To speed up the isomorphic cubes removal process, the same invariant-based algorithm described in [2] to remove isomorphic models can be applied to cubes. Invariants, such as number of distinct domain elements, are properties that must be identical for cubes to be isomorphic. For example, the cubes $A = \langle f(2,2) = 2; f(2,3) = 4 \rangle$ and $B = \langle f(3,3) = 2; f(1,2) = 2 \rangle$ cannot be isomorphic because $A$ contains an idempotent 2 but $B$ does not. Powerful and inexpensive invariants for binary operations include:
1. Number of y such that $x = (x * y) * x$.

2. Number of y such that $y = x * z$ for all $z \in D$.

3. Number of y such that $y = z * y$ for all $z \in D$.

4. Number of idempotents $x$ (i.e. $x * x = x$) for all $x \in D$.

5. Number of y such that $y * y = x$ for each $x \in D$.

First, invariant vectors (i.e. ordered lists of invariants) for cubes are calculated and used as hash keys to group cubes having the same invariant vectors into hash buckets. Then, cubes within the same bucket are tested for isomorphism. There is no need to test for isomorphism across buckets because isomorphic cubes must have the same invariant vectors. This saves tremendous amounts of testing time. Furthermore, buckets can be processed independently and in parallel to further speed up the process.

## 4.2 Work Stealing

In the basic form of this cube-based parallel algorithm, cubes are statically generated before the model enumeration process begins. It has the advantage of low runtime overheads as no synchronization among running finite model finders is needed. The preprocessing time for generating the cubes is also small for short to medium-length cubes. The disadvantage is that the workload may be uneven among the parallel processes. Some jobs may take a long time to finish when free workers sitting idle after finishing their jobs.

This problem can be solved with *work stealing* algorithms (also used in SAT [26]) in which a busy finite model searcher releases cubes that are not currently being worked on. For example, suppose a running model searcher is working on a cube $B_0 = \langle f(0,0) = 0 \rangle$, and its cell selection strategy picks the cell $f(1,1)$ to assign value next. Under the LNH, $f(1,1)$ may be assigned a value from $\{0,1,2\}$. If the model searcher is requested to spin out some work for other free workers, then it generates three cubes, $B_1 = \langle f(0,0) = 1; f(1,1) = 0 \rangle$, $B_2 = \langle f(0,0) = 1; f(1,1) = 1 \rangle$, and $B_3 = \langle f(0,0) = 1; f(1,1) = 2 \rangle$. It continues to work on the cube $B_0$ and releases $B_1$ and $B_2$ to other free workers.

# 5 Experimental Results

We integrate the cube-based algorithms into the finite model enumerator Mace4, which supports searching on FOL with the LNH and many cell selection strategies [35]. Parallelization is controlled outside Mace4. Only minor changes are made to Mace4 to

1. Accept cubes as inputs and continue searching for longer cubes or models from them.

2. Periodically check for signal for work stealing to spin off cubes for other workers.

The model searching logic in Mace4 remains intact. The concentric cell selection strategy (see Example 4 for its definition) is used in the experiments. A separate program removes isomorphic cubes by separating the cubes with equal invariants then check for isomorphisms (two cubes are isomorphic if one can be transformed to the other by a permutation).

We run the experiments on an Intel® Xeon® Silver 4110 CPU 2.0 GHz ×32 computer, with 64 GB of random access memory (RAM), using 30 parallel processes unless otherwise stated. All times reported are wall clock times.

We pick many disparate and challenging problems from the MarcieX database [3], which contains a collection of 158 most popular algebras. We also draw an example of semigroup subvariety from [1]. The definitions of the algebras used in the experiments in this section are listed in Table 1, in which all clauses are implicitly universally quantified.
8:10 Symmetries for Cube-And-Conquer in Finite Model Finding

<table><caption>Table 1 Definitions of Algebras Used in Experiments.</caption><thead><tr><th>Algebra</th><th>FOL Definition</th></tr></thead><tbody><tr><td>Semigroups</td><td>x * (y * z) = (x * y) * z.</td></tr><tr><td>Loops</td><td>x * y = x * z → y = z. y * x = z * x → y = z. x * 0 = x. 0 * x = x.</td></tr><tr><td>var{N<sub>2</sub><sup>1</sup> ∩ [x<sup>2</sup> = y<sup>2</sup>]}</td><td>x * (y * z) = (x * y) * z. (x * x) * x = x * x. x * y = y * x. x * x = y * y.</td></tr><tr><td>Tarski Algebras</td><td>(x * y) * y = (y * x) * x. x * (y * z) = y * (x * z). (x * y) * x = x.</td></tr><tr><td>Quasi-ordered Set</td><td>x &lt; y ∧ y &lt; z → x &lt; z. x &lt; x.</td></tr><tr><td>Involutive Lattices</td><td>(x * y) * z = x * (y * z). x * y = y * x. (x + y) + z = x + (y + z).<br/>x + y = y + x. (x * y) + x = x. (x + y) * x = x.<br/>-(x + y) = -x * -y. - - x = x.</td></tr></tbody></table>

In the tables showing experimental results in this section, the rows with cube length 0 show the results of running Mace4 in a single thread without the cube-based algorithms.

Table 2 shows the results of applying Theorem 9 to remove isomorphic cubes for the binary operation of the semigroups of order 7. Observe that the percentage reduction of the number of cubes increases as the cube length increases. The isomorphic cubes removal algorithm is therefore complementary to the LNH because the LNH removes a lot of short cubes but loses its effectiveness as the length of the cubes grows.

<table><caption>Table 2 #Cubes for Semigroups of Order 7.</caption><thead><tr><th rowspan="2">Cube Length</th><th colspan="2"># Cubes</th><th rowspan="2">Reduction (%)</th></tr><tr><th>w/o Removal of Isomorphic Cubes</th><th>w/ Isomorphic Cubes Removed</th></tr></thead><tbody><tr><td>2</td><td>6</td><td>5</td><td>16.7</td></tr><tr><td>4</td><td>34</td><td>28</td><td>17.6</td></tr><tr><td>9</td><td>1,568</td><td>888</td><td>43.4</td></tr><tr><td>16</td><td>56,206</td><td>12,036</td><td>78.2</td></tr><tr><td>25</td><td>1,028,171</td><td>59,056</td><td>94.3</td></tr></tbody></table>

We run Mace4 to enumerate semigroups defined by a single binary operation. The results show a speedup of over 100 times when cubes of length 25 are used, with over 96% of the isomorphic models suppressed (see Table 3). The results on semigroups are indicative of the algorithm's usefulness in general to the computational algebraists because algebraic structures related to semigroups are ubiquitous in algebra. Not only are there many well-known semigroup-related algebras, but also many semigroup varieties and subvarieties that are of high research interests [1].

Table 4 shows the results for loops (a quasigroup-related class of algebra) defined by a single non-associative binary operation. Here the reduction in the number of the output isomorphic models is not as pronounced. This is expected because the LNH works very well with the Latin square and removes a high percentage of the isomorphic models [48] before the isomorphic cubes removal takes place. For example, while only 0.16% of semigroups of order 7 generated by the LNH are non-isomorphic, 8.7% (106,228,849 out of 1,216,226,816) of the models generated for the loops of order 8 under the LNH are non-isomorphic. Nevertheless, the parallel algorithm provides 15 times improvement in speed for cube length of 16.
Table 5 shows the results of running the algorithms on the semigroup subvariety var{$N_2^1 \cap [x^2 = y^2]}$} (see p. 40 of [1] for its definition and discussions). With longer cubes, the algorithms speed up the process by 26 times with 30 threads. The results confirm that the proposed algorithms work remarkably well with semigroup-related algebras.

The Tarski algebras are unlike both the semigroups and the quasigroups in that its multiplication table is not associative and is not a Latin square [3]. It shows the cube-based algorithms perform better and better as the length of the cube increases (see Table 6).

The quasi-ordered set is defined by one binary relation. The isomorphic cubes algorithms work well on relations just as it works well on functions. As shown in Table 7, when cubes of length 36 are used, over 99% of the isomorphic models are suppressed, and the search process is sped up by over 200 times.

As an example to demonstrate the effectiveness of the algorithms on more complex algebras, consider the Involutive Lattice [3], which is defined by two associative binary operations and one unary operation. For Involutive Lattices of order 13, the search tree has a maximum depth of 351. Using cubes of length of 105, we obtain a speedup of 300 times, with almost 98% of the isomorphic cubes suppressed (see Table 8). The results show that the isomorphic cubes algorithms are highly effective for both simple and complex algebras.

<table><caption>Table 3 Running Cubes on Semigroups of Order 7.</caption><thead><tr><th rowspan="2">Cube Length</th><th rowspan="2">#Cubes</th><th rowspan="2">#Models (Millions)</th><th colspan="2">Time in min.</th></tr><tr><th>Gen.</th><th>Total Cubes</th></tr></thead><tbody><tr><td>0</td><td></td><td>1,021.1</td><td></td><td>235.2</td></tr><tr><td>2</td><td>5</td><td>717.7</td><td>0.0</td><td>12.5</td></tr><tr><td>4</td><td>28</td><td>611.1</td><td>0.1</td><td>9.4</td></tr><tr><td>9</td><td>888</td><td>360.2</td><td>0.1</td><td>5.2</td></tr><tr><td>16</td><td>12,036</td><td>158.2</td><td>0.2</td><td>2.8</td></tr><tr><td>25</td><td>59,056</td><td>39.5</td><td>0.9</td><td>1.7</td></tr></tbody></table>

<table><caption>Table 5 Running Cubes on var{$N_2^1 \cap [x^2 = y^2]}$} of Order 9.</caption><thead><tr><th rowspan="2">Cube Length</th><th rowspan="2">#Cubes</th><th rowspan="2">#Models (Millions)</th><th colspan="2">Time in min.</th></tr><tr><th>Gen.</th><th>Total Cubes</th></tr></thead><tbody><tr><td>0</td><td></td><td>313.0</td><td></td><td>72.0</td></tr><tr><td>2</td><td>1</td><td>156.5</td><td>0.0</td><td>2.9</td></tr><tr><td>4</td><td>1</td><td>156.5</td><td>0.1</td><td>2.8</td></tr><tr><td>9</td><td>2</td><td>156.5</td><td>0.1</td><td>2.8</td></tr><tr><td>16</td><td>5</td><td>120.9</td><td>0.1</td><td>2.3</td></tr><tr><td>25</td><td>16</td><td>55.5</td><td>0.2</td><td>1.3</td></tr><tr><td>36</td><td>70</td><td>13.0</td><td>0.3</td><td>0.8</td></tr><tr><td>49</td><td>331</td><td>1.5</td><td>1.0</td><td>1.1</td></tr></tbody></table>

<table><caption>Table 4 Running Cubes on Loops of Order 8.</caption><thead><tr><th rowspan="2">Cube Length</th><th rowspan="2">#Cubes</th><th rowspan="2">#Models (Millions)</th><th colspan="2">Time in min.</th></tr><tr><th>Gen.</th><th>Total Cubes</th></tr></thead><tbody><tr><td>0</td><td></td><td>1,216</td><td></td><td>564.0</td></tr><tr><td>2</td><td>1</td><td>1,216</td><td>0.0</td><td>47.4</td></tr><tr><td>4</td><td>2</td><td>1,216</td><td>0.1</td><td>47.3</td></tr><tr><td>9</td><td>18</td><td>1,216</td><td>0.1</td><td>46.2</td></tr><tr><td>16</td><td>3,583</td><td>1,214</td><td>0.1</td><td>45.3</td></tr></tbody></table>

<table><caption>Table 6 Running Cubes on Tarski Algebras of Order 13.</caption><thead><tr><th rowspan="2">Cube Length</th><th rowspan="2">#Cubes</th><th rowspan="2">#Models (Millions)</th><th colspan="2">Time in min.</th></tr><tr><th>Gen.</th><th>Total Cubes</th></tr></thead><tbody><tr><td>0</td><td></td><td>379.6</td><td></td><td>1,949.9</td></tr><tr><td>2</td><td>3</td><td>189.8</td><td>0.0</td><td>70.2</td></tr><tr><td>4</td><td>1</td><td>189.8</td><td>0.1</td><td>69.9</td></tr><tr><td>9</td><td>3</td><td>183.3</td><td>0.1</td><td>67.7</td></tr><tr><td>16</td><td>11</td><td>158.8</td><td>0.1</td><td>58.1</td></tr><tr><td>25</td><td>55</td><td>111.9</td><td>0.2</td><td>40.1</td></tr><tr><td>36</td><td>157</td><td>62.1</td><td>0.2</td><td>21.8</td></tr><tr><td>49</td><td>174</td><td>24.9</td><td>0.5</td><td>8.8</td></tr><tr><td>64</td><td>171</td><td>6.6</td><td>1.0</td><td>3.7</td></tr></tbody></table>
## 8:12 Symmetries for Cube-And-Conquer in Finite Model Finding

<table><caption>**Table 7** Running Cubes on Quasi-ordered Set of Order 8.</caption><thead><tr><th rowspan="2">Cube Length</th><th rowspan="2">#Cubes</th><th rowspan="2">#Models (Millions)</th><th colspan="2">Time in min.</th></tr><tr><th>Gen.</th><th>Total Cubes</th></tr></thead><tbody><tr><td>0</td><td></td><td>642.8</td><td></td><td>59.9</td></tr><tr><td>2</td><td>1</td><td>642.8</td><td>0.0</td><td>4.2</td></tr><tr><td>4</td><td>3</td><td>474.6</td><td>0.1</td><td>3.2</td></tr><tr><td>9</td><td>9</td><td>209.5</td><td>0.1</td><td>1.7</td></tr><tr><td>16</td><td>33</td><td>61.3</td><td>0.1</td><td>0.8</td></tr><tr><td>25</td><td>139</td><td>12.6</td><td>0.2</td><td>0.3</td></tr><tr><td>36</td><td>713</td><td>2.0</td><td>0.3</td><td>0.3</td></tr></tbody></table>

<table><caption>**Table 8** Running Cubes on Involutive Lattices of Order 13.</caption><thead><tr><th rowspan="2">Cube Length</th><th rowspan="2">#Cubes</th><th rowspan="2">#Models (Millions)</th><th colspan="2">Time in min.</th></tr><tr><th>Gen.</th><th>Total Cubes</th></tr></thead><tbody><tr><td>0</td><td></td><td>423.0</td><td></td><td>4,719.7</td></tr><tr><td>3</td><td>2</td><td>423.0</td><td>0.0</td><td>432.5</td></tr><tr><td>6</td><td>3</td><td>423.0</td><td>0.1</td><td>432.8</td></tr><tr><td>10</td><td>6</td><td>263.9</td><td>0.1</td><td>270.0</td></tr><tr><td>21</td><td>23</td><td>178.6</td><td>0.1</td><td>180.9</td></tr><tr><td>36</td><td>108</td><td>84.9</td><td>0.2</td><td>88.3</td></tr><tr><td>55</td><td>555</td><td>46.0</td><td>0.3</td><td>46.2</td></tr><tr><td>78</td><td>1,710</td><td>19.8</td><td>0.5</td><td>20.6</td></tr><tr><td>105</td><td>5,048</td><td>8.7</td><td>4.9</td><td>14.3</td></tr></tbody></table>

The reductions in time and number of models (on top of the LNH) are summarized in Figures 3 and 4. Note that the reduction in total time is over 90% even for short cubes. However, the biggest gain in both reduction in time and in isomorphic models is when longer cubes are used. Reduction in isomorphic models also helps tremendously in the post processing step to extract non-isomorphic models.

<figure><img src="image_67.png" /><figcaption>**Figure 3** Reduction in Number of Output Models.</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>**Figure 4** Reduction in Total Time with 30 Parallel Processes.</figcaption></figure>

### 5.1 Speedup of Finite Model Enumeration with Parallelization

As discussed, the cubes algorithms allow low-cost parallelization of the finite model enumeration process. Figure 5 and Table 9 show the performance of the parallel cubes algorithms with 1 to 16 parallel processes. Here, the reported times do not include isomorphic mode filtering; they are for Mace4 to generate models only. Note that when many processes compete for limited amount of RAM, swapping could slow down the processes substantially. This helps to explain why larger algebras, such as the Involutive Lattice of order 13, have their curves flattened out much faster than small algebras, such as the Semigroups of order 7. More processes also mean more work-stealing and higher overheads.
<table><caption>Table 9 Performance w/ Multiprocessing.</caption><thead><tr><th rowspan="2">Algebra</th><th rowspan="2">Order</th><th rowspan="2">Cube Length</th><th colspan="5">Time in seconds</th></tr><tr><th>1 Process</th><th>2 Processes</th><th>4 Processes</th><th>8 Processes</th><th>16 Processes</th></tr></thead><tbody><tr><td>Semigroups</td><td>7</td><td>25</td><td>6,626</td><td>3,397</td><td>1,757</td><td>940</td><td>425</td></tr><tr><td>Loops</td><td>7</td><td>16</td><td>202</td><td>108</td><td>50</td><td>35</td><td>21</td></tr><tr><td>Tarski algebras</td><td>13</td><td>64</td><td>1,766</td><td>973</td><td>552</td><td>273</td><td>250</td></tr><tr><td>var{N<sub>2</sub><sup>1</sup> ∩ [x<sup>2</sup> = y<sup>2</sup>]}</td><td>9</td><td>49</td><td>130</td><td>84</td><td>80</td><td>57</td><td>53</td></tr><tr><td>Quasi Ordered</td><td>8</td><td>36</td><td>123</td><td>77</td><td>51</td><td>37</td><td>25</td></tr><tr><td>Involutive Lattices</td><td>12</td><td>105</td><td>1,496</td><td>794</td><td>480</td><td>378</td><td>320</td></tr></tbody></table>

<figure><img src="image_69.png" /><figcaption>Figure 5 Performance w/ Multiprocessing.</figcaption></figure>

## 5.2 Isomorphic Cubes Removal Speeds up Isomorphic Models Filtering

As pointed out Section 1, reducing the number of Mace4 outputs also reduces the efforts needed to filter out isomorphic models. Table 10 shows, using involutive lattices as an example, the out-sized effect of the reduction of Mace4 outputs on the time to filter out the isomorphic models using the invariant-based isomorphic model filtering algorithm [2], with 30 parallel processes. With the reduction in number of Mace4 models, the isomorphic model filtering process is sped up by 2 orders of magnitude. The improvement in speed is observed to be better with models of higher orders. We would also point out that the isomorphic model filter generates the same non-isomorphic models with or without the cubes algorithms.

<table><caption>Table 10 Running Invariant-based Isomorphic Models Filter on Involutive Lattices.</caption><thead><tr><th rowspan="2">Order</th><th rowspan="2">#Non-iso Models</th><th colspan="2">w/o Cubes</th><th colspan="2">w/ Cubes</th></tr><tr><th>#Mace4 Output</th><th>Isomorphic Model Filter Time (s)</th><th>Cube Length</th><th>#Mace4 Output</th><th>Isomorphic Model Filter Time (s)</th></tr></thead><tbody><tr><td>9</td><td>122</td><td>72,470</td><td>29</td><td>78</td><td>3,670</td><td>1</td></tr><tr><td>10</td><td>389</td><td>575,463</td><td>496</td><td>105</td><td>13,789</td><td>4</td></tr><tr><td>11</td><td>906</td><td>4,771,035</td><td>28,424</td><td>105</td><td>97,680</td><td>135</td></tr><tr><td>12</td><td>3,047</td><td>43,851,030</td><td>N/A</td><td>105</td><td>971,416</td><td>2,802</td></tr></tbody></table>
## 5.3 Optimal Cube Length

In general, the search process using longer cubes finishes earlier with fewer isomorphic models. However, we observe that there are three limiting factors on the lengths of the cubes.

First, as the length of the cubes gets longer, more and more models are generated as a result of propagations. This reduces the impact of removing isomorphic cubes because they represent a progressively smaller proportion of the isomorphic models. It is observed that when more than *n* − 2 symbols out of the *n* domain elements are used in the cell terms, the number of (isomorphic) models will be substantial and extending the cube length does not bring enough reduction in isomorphic models to justify the increase in processing time.

Second, the isomorphic cubes removal time grows quite fast as the length of the cube grows. When the isomorphic cubes removal process takes more than a few minutes, further lengthening of the cubes will result in prohibitive overheads in the search process.

Lastly, when the final number of cubes is more than tens of thousands, the overheads in processing them becomes so high that the search becomes slower. This factor depends heavily on the number of processors available. More processors mean more parallel processes can be run without slowing down the whole search process.

One heuristic is to run cube generation until the number of cubes reaches some threshold or the runtime exceeds some threshold, then switch to model generation. The thresholds are system-dependent and can further be fine-tuned by experiments with algebras of interest.

## 6 Related Work

There is extensive research on *paralyzing SAT solving*, where the predominant approaches are search space partitioning and portfolios, c.f. [33]. We find inspiration in the *cube-and-conquer* approach proposed by Heule and colleagues [20–22], where the search space is partitioned by a lookahead solver into (many) cubes and then each subspace is solved by a CDCL SAT solver. In SAT, partitioning by a CDCL solver is nontrivial [32] and that is why the lookahead solver is useful for this task. Nevertheless, the use of the lookahead solver is not seen as an indispensable feature of the cube-and-conquer, as noted by Subercaseaux and Heule [46]. In our approach, we have a tight control over the decisions of the solver and we do not need a separate solver to perform the splitting. Additionally, we invest extra effort into search space splitting by identifying symmetries in the cubes.

The *adaptive prefix-assignment technique* [25] is a symmetry reduction algorithm used in SAT. The prefix is equivalent to a propositional cube, and the algorithm also tries to eliminate isomorphic cubes. In our case, we exploit symmetries specific to FOL-LNH and isomorphism at FOL level, which is absent in their algorithm (and in SAT in general).

*Parallel algorithms* can be characterized by how the search is done. There are two main search methods: embarrassingly parallel search (EPS) and work stealing search [6, 10, 25, 33, 40, 42]. In the former method, the task is decomposed into many sub-tasks that are queued up to be processed by free worker threads/processes. In the latter method, when a worker completes its task, it asks other workers for more work. The busy workers may split their tasks into smaller sub-tasks and pass some of them to the free workers. The main focus of this method is to keep all the CPUs running until all jobs are done, although for some cases, the work stealing scheme can affect efficiency [10]. The EPS method is a natural choice for the cube-based parallelization scheme because preprocessing can be performed to generate numerous non-isomorphic cubes by splitting the search space. However, a work-stealing procedure is essential in supplementing the EPS to balance uneven workloads [33].
Parallel algorithms can also help select the best strategy in solving a problem with the EPS method [39]. After a problem is decomposed into a large number of sub-tasks, a small number (e.g., 1%) of these sub-tasks are run in parallel using different strategies of the same solver or different solvers. The strategy that gives the best performance on the subset of sub-tasks will be used to run all sub-tasks. The same idea is used in the invariant-based isomorphic models removal algorithm [2]: it randomly generates a large number of invariants, then applies them to a small percentage of models to pick the best performing random invariants to apply to the whole set of models. This idea can be applied to the finite model finders that support multiple cell selection strategies to pick the best function order and cell selection strategy for any specific problem.

Finite model enumeration can be posed as a *constraint programming (CP)* task [27]. Some CP solvers, e.g., Minion [17] and Gecode [37], support parallelization [31]. In CP, the search space can be partitioned by adding constraints to rule in and/or out partitions. Each partition can be processed by a separate worker thread/process. Minion further implements a work stealing search scheme that also partitions the search space dynamically by splitting the existing constraint model after the search has started [15, 29]. However, to effectively add *symmetry-breaking* constraints such as lex-leaders to a CP solver often requires deep knowledge of the solver *and* the problem at hand (e.g., the semigroups in [15]) which may not be available when mathematicians first define and study a new algebraic structure.

Moreover, to use traditional CP solvers for finite model enumerations, mathematicians need to learn a new CP-specific language such as CHR [45] and Savile Row [38]. It is possible to use a translator to translate between languages, but that adds uncertainties to the fidelity and the optimality of the translated specifications. FOL remains one of the most popular languages among mathematicians due to its simple and intuitive syntax. Moreover, a popular automatic theorem prover, Prover9 [34], shares the same input language with Mace4. This adds more than just convenience to the process, as it also reduces the chances of discrepancies between Prover9 and Mace4 on the same problem.

A well-known issue with enumerating models defined with FOL are the isomorphic models included in the outputs. This is an inherent symmetry property of FOL [40]. There is extensive research on *symmetry-breaking* [4, 11–13, 28, 40, 43, 47]. Although complete symmetry-breaking is known to be computationally challenging [13, 47], many useful algorithms, such as the LNH and the XLNH [4, 5], have emerged in partial symmetry-breaking. The LNH can be considered a symmetry-breaking with interchangeable values in constraint satisfaction problems (CSP) [19]. The XLNH is more restrictive as it only works on unary operations. The LNH is implemented in many systems such as Falcon [50], SEM [51], FMSET [6], and Mace4. The isomorphic cubes algorithm, which removes more cubes as the cube length grows, complements the LNH.

Another important symmetry-breaking strategy is to steer the search engine away from the fruitless exploration of sub-search space by adding symmetry-breaking input clauses [13, 47]. The cube-based parallel algorithms are compatible with algorithms of this kind of strategy as long as they do not break the LNH.

Some finite model finders, such as SEMK [8] and SEMD [24], try to completely suppress isomorphic models in the search process. However, these isomorph-free algorithms are not easy to parallelize as global information is generated and consumed in many steps, requiring high-cost synchronizations between cooperating workers, especially when they run on different computers. The cube-based parallel algorithm, on the other hand, is an EPS method that requires no synchronizations between workers. The static removal of isomorphic cubes done in a preprocessing step is shown to be effective in suppressing isomorphic models even
before the actual search begins. The augmented work stealing algorithm is not high in
synchronization costs because it does not involve communications between running jobs. The
remaining isomorphic models from the cube-based algorithms can be efficiently removed by
the invariant-based isomorphic model filtering algorithm as a postprocessing step.

Another algorithm, DSYM [4], exploits local symmetries by finding symmetries (synonym to isomorphisms in their terminology) under invariant partial interpretations (which are invariant cubes) and without parallelism. It also works with the LNH and XLNH. DSYM is a predictive algorithm that works at the *parent* level and predicts which of its immediate children will be isomorphic cubes. It can be seen as a special case of the isomorphic cube algorithm because it removes isomorphic cubes having the same immediate parents, while the isomorphic cube algorithm removes all isomorphic cubes, irrespective of their parents. Nevertheless, for the cases that DSYM covers, it does so right before the cubes are generated, while the isomorphic cube algorithm only detects the symmetries right after the cubes are generated. A disadvantage of DSYM is that it is not clear how it can be effectively parallelized. Furthermore, DSYM only detects symmetries under the same subtree. The isomorphic cubes removal algorithm, on the other hand, detects both global and local symmetries the same way, and hence detects and removes more symmetries than DSYM. Moreover, DSYM uses only two invariants in testing isomorphism between cubes, while we use many invariants that are proven successful in the invariant-based isomorphic model removal algorithm in our isomorphic cubes removal process. Nevertheless, DSYM can be applied to the cube generation process as well as the final model generation process. That is, the isomorphic cube removal algorithm is compatible with DSYM, as with any other symmetry-breaking algorithm that works with the LNH.

7 Conclusions and Future Work

In this paper, we introduce an efficient parallel algorithm together with a novel symmetry-
removal mechanism for enumerating finite models. The approach is inspired by the cube-and-
conquer paradigm, successfully used in SAT solving, which partitions the search space into
cubes and then massively paralyzes. In contrast, our approach applies symmetries specific to
finite model finding.

In conclusion, this paper fulfills an important unmet need for an efficient algorithm for enumerating finite algebraic models in computational algebra by enhancing the existing finite model enumeration process with the parallel cubes algorithm and the isomorphic cubes removal algorithm that reduce both the runtime and the number of output isomorphic models. These new algorithms are so scalable that they can be used on a laptop as well as on a cluster of powerful computers, and they require minimal efforts to safely integrate into existing finite model finders. Very importantly, these algorithms can be used as a black-box without requiring the users to have any knowledge about the way they work.

Future research will focus on improving isomorphic cube removal, on best cell selection
strategy, and on predicting of optimal cube length.

References

1 João Araújo, João Pedro Araújo, Peter J. Cameron, Edmond W. H. Lee, and Jorge Raminhos. A survey on varieties generated by small semigroups and a companion website, 2019. arXiv: 1911.05817.

2 João Araújo, Choiwah Chow, and Mikoláš Janota. Boosting isomorphic model filtering with invariants. *Constraints*, 27(3):360–379, July 2022. doi:10.1007/s10601-022-09336-x.
3 João Araújo, David Matos, and João Ramires. MarcieDB: a model and theory database. https://marciedb.pythonanywhere.com, 2022.

4 Gilles Audemard, Belaid Benhamou, and Laurent Henocque. Predicting and detecting symmetries in FOL finite model search. *J. Autom. Reason.*, 36(3):177–212, 2006. doi:10.1007/s10817-006-9040-3.

5 Gilles Audemard and Laurent Henocque. The eXtended least number heuristic. In Rajeev Goré, Alexander Leitsch, and Tobias Nipkow, editors, *Automated Reasoning, First International Joint Conference, IJCAR*, volume 2083 of *Lecture Notes in Computer Science*, pages 427–442, Berlin, Heidelberg, 2001. Springer. doi:10.1007/3-540-45744-5_35.

6 Belaid Benhamou and Laurent Henocque. A hybrid method for finite model search in equational theories. *Fundam. Informaticae*, 39(1-2):21–38, 1999. doi:10.3233/FI-1999-391202.

7 R.D. Blumofe and C.E. Leiserson. Scheduling multithreaded computations by work stealing. In *Proceedings 35th Annual Symposium on Foundations of Computer Science*, pages 356–368, 1994. doi:10.1109/SFCS.1994.365680.

8 Thierry Boy de la Tour and Prakash Countcham. An isomorph-free SEM-like enumeration of models. *Electronic Notes in Theoretical Computer Science*, 125(2):91–113, 2005. Proceedings of the 5th International Workshop on Strategies in Automated Deduction (Strategies 2004). doi:10.1016/j.entcs.2005.01.003.

9 Stanley Burris and Hanamantagouda P. Sankappanavar. *A course in universal algebra*, volume 78 of *Graduate texts in mathematics*. Springer, New York, NY, 1981.

10 Geoffrey Chu, Christian Schulte, and Peter J. Stuckey. Confidence-based work stealing in parallel constraint programming. In Ian P. Gent, editor, *Principles and Practice of Constraint Programming - CP*, volume 5732, pages 226–241. Springer, 2009. doi:10.1007/978-3-642-04244-7_20.

11 Koen Claessen and Niklas Sörensson. New techniques that improve MACE-style finite model finding. In *Proceedings of the CADE-19 Workshop: Model Computation - Principles, Algorithms, Applications*, 2003.

12 Michael Codish, Alice Miller, Patrick Prosser, and Peter James Stuckey. Breaking symmetries in graph representation. In Francesca Rossi, editor, *IJCAI 2013, Proceedings of the 23rd International Joint Conference on Artificial Intelligence*, pages 510–516. IJCAI/AAAI, 2013. URL: http://www.aaai.org/OCS/index.php/IJCAI/IJCAI13/paper/view/6480.

13 James M. Crawford, Matthew L. Ginsberg, Eugene M. Luks, and Amitabha Roy. Symmetry-breaking predicates for search problems. In Luigia Carlucci Aiello, Jon Doyle, and Stuart C. Shapiro, editors, *Proceedings of the Fifth International Conference on Principles of Knowledge Representation and Reasoning (KR)*, pages 148–159. Morgan Kaufmann, 1996.

14 A. Distler and J. Mitchell. Smallsemi, a library of small semigroups in GAP, Version 0.6.12. https://gap-packages.github.io/smallsemi/, 2019. GAP package.

15 Andreas Distler, Christopher Jefferson, Tom Kelsey, and Lars Kotthoff. The semigroups of order 10. In Michela Milano, editor, *Principles and Practice of Constraint Programming - CP*, volume 7514, pages 883–899. Springer, 2012. doi:10.1007/978-3-642-33558-7_63.

16 The GAP Group. *GAP – Groups, Algorithms, and Programming, Version 4.11.1*, 2021. URL: https://www.gap-system.org.

17 Ian P. Gent, Christopher Jefferson, and Ian Miguel. Minion: A fast scalable constraint solver. In Gerhard Brewka, Silvia Coradeschi, Anna Perini, and Paolo Traverso, editors, *ECAI, 17th European Conference on Artificial Intelligence, Including Prestigious Applications of Intelligent Systems (PAIS), Proceedings*, volume 141 of *Frontiers in Artificial Intelligence and Applications*, pages 98–102, Amsterdam, Netherlands, 2006. IOS Press. URL: http://www.booksonline.iospress.nl/Content/View.aspx?pid=1654.

18 Ian P. Gent, Ian Miguel, Peter Nightingale, Ciaran McCreesh, Patrick Prosser, Neil C. A. Moore, and Chris Unsworth. A review of literature on parallel constraint solving. *Theory Pract. Log. Program.*, 18(5-6):725–758, 2018. doi:10.1017/S1471068418000340.
8:18 Symmetries for Cube-And-Conquer in Finite Model Finding

19 Pascal Van Hentenryck, Pierre Flener, Justin Pearson, and Magnus Ågren. Tractable symmetry breaking for CSPs with interchangeable values. In Georg Gottlob and Toby Walsh, editors, *IJCAI-03, Proceedings of the Eighteenth International Joint Conference on Artificial Intelligence*, pages 277–284. Morgan Kaufmann, 2003. URL: http://ijcai.org/Proceedings/03/Papers/041.pdf.

20 Marijn Heule, Oliver Kullmann, Siert Wieringa, and Armin Biere. Cube and conquer: Guiding CDCL SAT solvers by lookaheads. In Kerstin Eder, João Lourenço, and Onn Shehory, editors, *Hardware and Software: Verification and Testing - 7th International Haifa Verification Conference, HVC, Revised Selected Papers*, volume 7261, pages 50–65. Springer, 2011. doi:10.1007/978-3-642-34188-5_8.

21 Marijn J. H. Heule, Oliver Kullmann, and Armin Biere. Cube-and-conquer for satisfiability. In Youssef Hamadi and Lakhdar Sais, editors, *Handbook of Parallel Constraint Reasoning*, pages 31–59. Springer, 2018. doi:10.1007/978-3-319-63516-3_2.

22 Marijn J. H. Heule, Oliver Kullmann, and Victor W. Marek. Solving and verifying the Boolean Pythagorean triples problem via cube-and-conquer. In *Theory and Applications of Satisfiability Testing (SAT)*, 2016. doi:10.1007/978-3-319-40970-2_15.

23 Antti E. J. Hyvärinen, Matteo Marescotti, and Natasha Sharygina. Lookahead in partitioning SMT. In *Formal Methods in Computer Aided Design, FMCAD*, pages 271–279. IEEE, 2021. doi:10.34727/2021/isbn.978-3-85448-046-4_37.

24 Xiangxue Jia and Jian Zhang. A powerful technique to eliminate isomorphism in finite model search. In Ulrich Furbach and Natarajan Shankar, editors, *Automated Reasoning*, pages 318–331, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg.

25 Tommi Junttila, Matti Karppa, Petteri Kaski, and Jukka Kohonen. An adaptive prefix-assignment technique for symmetry reduction. *Journal of Symbolic Computation*, 99:21–49, 2020. doi:10.1016/j.jsc.2019.03.002.

26 Bernard Jurkowiak, Chu Min Li, and Gil Utard. A parallelization scheme based on work stealing for a class of SAT solvers. *J. Autom. Reason.*, 34(1):73–101, 2005. doi:10.1007/s10817-005-1970-7.

27 Majid Ali Khan. Efficient enumeration of higher order algebraic structures. *IEEE Access*, 8:41309–41324, 2020. doi:10.1109/ACCESS.2020.2976876.

28 Markus Kirchweger and Stefan Szeider. SAT modulo symmetries for graph generation. In Laurent D. Michel, editor, *27th International Conference on Principles and Practice of Constraint Programming, CP*, volume 210 of LIPIcs, pages 34:1–34:16. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021. doi:10.4230/LIPIcs.CP.2021.34.

29 Lars Kotthoff and Neil C. A. Moore. Distributed solving through model splitting. *ArXiv*, abs/1008.4328, 2010.

30 Kenneth Kunen. The structure of conjugacy closed loops. *Transactions of the American Mathematical Society*, 352(6):2889–2911, 2000.

31 Arnaud Malapert, Jean-Charles Régin, and Mohamed Rezgui. Embarrassingly parallel search in constraint programming. *J. Artif. Intell. Res.*, 57:421–464, 2016. doi:10.1613/jair.5247.

32 Ruben Martins, Vasco M. Manquinho, and Inês Lynce. Improving search space splitting for parallel SAT solving. In *22nd IEEE International Conference on Tools with Artificial Intelligence, ICTAI*. IEEE Computer Society, 2010. doi:10.1109/ICTAI.2010.56.

33 Ruben Martins, Vasco M. Manquinho, and Inês Lynce. An overview of parallel SAT solving. *Constraints An Int. J.*, 17(3):304–347, 2012. doi:10.1007/s10601-012-9121-3.

34 W. McCune. Prover9 and mace4. http://www.cs.unm.edu/~mccune/prover9/, 2005–2010.

35 William McCune. Mace4 reference manual and guide, August 2003. URL: https://www.cs.unm.edu/~mccune/prover9/mace4.pdf.

36 Gábor Nagy and Petr Vojtěchovský. LOOPS, computing with quasigroups and loops in GAP, Version 3.4.1. https://gap-packages.github.io/loops/, November 2018. Refereed GAP package.

37 Morten Nielsen. Parallel search in gecode. *Technical Report, Gecode*, 2006.
38 Peter Nightingale, Özgür Akgün, Ian P. Gent, Christopher Jefferson, Ian Miguel, and Patrick Spracklen. Automatically improving constraint models in savile row. *Artificial Intelligence*, 251:35-61, 2017. doi:10.1016/j.artint.2017.07.001.

39 Anthony Palmieri, Jean-Charles Régin, and Pierre Schaus. Parallel strategies selection. *CoRR*, abs/1604.06484, 2016. arXiv:1604.06484.

40 Giles Reger, Martin Riener, and Martin Suda. Symmetry avoidance in MACE-style finite model finding. In Andreas Herzig and Andrei Popescu, editors, *Frontiers of Combining Systems FroCoS*, volume 11715, pages 3-21, Switzerland AG, 2019. Springer. doi:10.1007/978-3-030-29007-8_1.

41 Jean-Charles Régin and Arnaud Malapert. Parallel constraint programming. In Youssef Hamadi and Lakhdar Sais, editors, *Handbook of Parallel Constraint Reasoning*, pages 337-379. Springer, 2018. doi:10.1007/978-3-319-63516-3_9.

42 Jean-Charles Régin, Mohamed Rezgui, and Arnaud Malapert. Embarrassingly parallel search. In Christian Schulte, editor, *Principles and Practice of Constraint Programming - 19th International Conference, CP 2013, Uppsala, Sweden, September 16-20, 2013. Proceedings*, volume 8124 of *Lecture Notes in Computer Science*, pages 596-610, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg. doi:10.1007/978-3-642-40627-0_45.

43 Francesca Rossi, Peter van Beek, and Toby Walsh, editors. *Handbook of Constraint Programming*, volume 2 of *Foundations of Artificial Intelligence*. Elsevier, 2006. URL: https://www.sciencedirect.com/science/bookseries/15746526/2.

44 Neil J. A. Sloane and The OEIS Foundation Inc. The on-line encyclopedia of integer sequences, 2020. URL: http://oeis.org/?language=english.

45 Jon Sneyers, Peter van Weert, Tom Schrijvers, and Leslie de Koninck. As time goes by: Constraint handling rules: A survey of chr research from 1998 to 2007. *Theory and Practice of Logic Programming*, 10(1):1-47, 2010. doi:10.1017/S1471068409990123.

46 Bernardo Subercaseaux and Marijn Heule. Toward optimal radio colorings of hypercubes via SAT-solving. In Ruzica Piskac and Andrei Voronkov, editors, *Proceedings of 24th International Conference on Logic for Programming, Artificial Intelligence and Reasoning*, volume 94 of *EPiC Series in Computing*, pages 386-404. EasyChair, 2023. doi:10.29007/qrmp.

47 Toby Walsh. Symmetry breaking constraints: Recent results. In Jörg Hoffmann and Bart Selman, editors, *Proceedings of the Twenty-Sixth AAAI Conference on Artificial Intelligence, July 22-26, 2012, Toronto, Ontario, Canada*. AAAI Press, 2012. URL: http://www.aaai.org/ocs/index.php/AAAI/AAAI12/paper/view/4974.

48 H. Zhang. Combinatorial designs by SAT solvers. *Handbook of Satisfiability*, pages 533-568, 2009. URL: https://cir.nii.ac.jp/crid/1571980076163512448.

49 Hantao Zhang and Jian Zhang. MACE4 and SEM: A comparison of finite model generators. In Maria Paola Bonacina and Mark E. Stickel, editors, *Automated Reasoning and Mathematics - Essays in Memory of William W. McCune*, volume 7788 of *Lecture Notes in Computer Science*, pages 101-130. Springer, 2013. doi:10.1007/978-3-642-36675-8_5.

50 Jian Zhang. Constructing finite algebras with FALCON. *Journal of Automated Reasoning*, 17:1-22, August 1996. doi:10.1007/BF00247667.

51 Jian Zhang and Hantao Zhang. SEM: a system for enumerating models. In *IJCAI*, pages 298-303, 1995. URL: http://ijcai.org/Proceedings/95-1/Papers/039.pdf.
<figure><img src="image_70.png" /></figure>
# Guiding Backtrack Search by Tracking Variables During Constraint Propagation

Gilles Audemard ☑️ ID

CRIL, Univ. Artois & CNRS, France

Christophe Lecoutre ☑

CRIL, Univ. Artois & CNRS, France

Charles Prud'homme ☑️ ID

TASC, IMT-Atlantique, LS2N-CNRS, France

## Abstract

It is well-known that variable ordering heuristics play a central role in solving efficiently Constraint Satisfaction Problem (CSP) instances. From the early 80's, and during more than two decades, the dynamic variable ordering heuristic selecting the variable with the smallest domain was clearly prevailing. Then, from the mid 2000's, some adaptive heuristics have been introduced: their principle is to collect some useful information during the search process in order to take better informed decisions. Among those adaptive heuristics, wdeg/dom (and its variants) remains particularly robust. In this paper, we introduce an original heuristic based on the midway processing of failing executions of constraint propagation: this heuristic called pick/dom tracks the variables that are directly involved in the process of constraint propagation, when ending with a conflict. The robustness of this new heuristic is demonstrated from a large experimentation conducted with the constraint solver ACE. Interestingly enough, one can observe some complementary between the early, midway and late forms of processing of conflicts.

**2012 ACM Subject Classification** Computing methodologies → Discrete space search

**Keywords and phrases** Variable Ordering Heuristics, Variable Weighting

**Digital Object Identifier** 10.4230/LIPIcs.CP.2023.9

**Funding** This work has benefited from the support of the National Research Agency under France 2030, MAIA Project ANR-22-EXES-0009.

## 1 Introduction

Backtrack search remains a classical approach for solving instances of the Constraint Satisfaction Problem (CSP), and the related Constraint Optimization Problem (COP). It is based on depth-first exploration, which is conducted by instantiating variables in sequence and backtracking when dead-ends occur. For efficiently exploring the search space, a filtering process is performed at each step of the search so as to reduce the domains of the variables; typically most of the constraints guarantee the property known as (generalized) arc consistency [5, 1, 26, 14]. The order in which variables are chosen during the depth-first traversal of the search space is decided by a *variable ordering heuristic H*. At each internal node of the search tree built by the backtrack search algorithm, the next variable *x* is selected by *H*, and a value is assigned to *x* according to a *value ordering heuristic*. Choosing the right heuristics for solving a given constraint network is a key issue since different heuristics can lead to drastically different search trees.

In modern constraint solvers, three main principles are considered for guiding search (i.e., performing depth-first exploration):

- First, one should start by assigning variables that belong to the most difficult part(s) of the problem instance. This principle is derived from the recognition that there is no point in traversing the easy part(s) of an instance and then backtracking repeatedly when

<figure><img src="image_71.png" /></figure>

© Gilles Audemard, Christophe Lecoutre, and Charles Prud'homme;
licensed under Creative Commons License CC-BY 4.0

29th International Conference on Principles and Practice of Constraint Programming (CP 2023).

Editor: Roland H. C. Yap; Article No. 9; pp. 9:1–9:17

<figure><img src="image_72.png" /></figure>

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
9:2 Guiding Backtrack Search by Tracking Variables During Constraint Propagation

it turns out that the first choices are incompatible with the remaining difficult part(s).
Here the underlying *fail-first* principle is [13]: “To succeed, try first where you are most
likely to fail”.

- Second, value selection should be based on the *succeed-first* or *promise* principle, which comes from the simple observation that to find a solution quickly, it is better to move at each step to the most promising subtree, primarily by selecting a value that is most likely to participate in a solution.

- Third, when starting to build the search tree, one should pay attention to the initial variable/value choices that are particularly important. Indeed, bad choices near the root of the search tree may turn out to be disastrous because they lead to exploration of very large fruitless subtrees. To make good initial choices, one strategy is to select the first branching decisions with special care, perhaps calling sophisticated and expensive procedures for this purpose. Another relevant strategy is to restart search several times, ideally learning some information each time in order to refine search guidance.

The current response (in solvers) to following these principles is as follows:

- Generic *adaptive* variable ordering heuristics that learn from conflicts during exploration are usually employed; a classical such heuristic being wdeg/dom [3], possibly combined with a mechanism, called last-conflict reasoning (lc) simulating a certain form of intelligent backjumps [18].

- Promising attempts to select values (or pairs variable-value) according to some elegant mechanisms such as Belief Propagation [23] have been introduced, but, unfortunately, we are not aware of any generic robust value ordering heuristic. One exception may be BIVS (Bound-Impact Value Selector) [8], but controlling its computation cost remains a difficult question. Consequently, it is rather frequent that the first (smallest) value be the default choice. Interestingly, for optimization, it is highly recommended to use in priority the value present in the last found solution, which is a technique known as solution(-based phase) saving [27, 7], clearly in concordance with the *promise* principle (as initially mentioned in [6, 10]).

- To address the issue of heavy-tailed runtime distributions [11], the search is restarted regularly, following a geometric progression (or the Luby sequence). Besides, by collecting nogoods [17] along the leftmost branch of the search tree at the end of each run (i.e., just before restarting), we have the guarantee of never exploring again the same parts of the search space (which is a nice feature when exhibiting all distinct solutions of a CSP instance).

This is the context of our contribution. More specifically, we focus our interest on the
first crucial component: variable ordering heuristics. Indeed, in this paper, we introduce
an original heuristic called *pick/dom* that learns from conflicts by identifying the variables
that are directly involved in the process of constraint propagation (when failing). We show
how to implement it in the variable-oriented propagation scheme. The robustness of this
new heuristic is demonstrated by conducting a large experimentation with the well-known
constraint solver ACE [15].

2 Preliminaries

A *Constraint Network* (CN) is composed of a finite set of *n* variables *X*, and a finite set of *e*
constraints *C*. Each variable *x* must be assigned a value from its current domain, denoted by
*dom*(*x*). Each constraint *c* represents a mathematical relation over an ordered set of variables,
called the *scope* of *c*, and denoted by **scp**(*c*). The *arity* of a constraint *c* is the size of its scope. The *degree* of a variable *x* is the number of constraints of *C* involving *x*. A *solution* to a CN *P* = (*X*, *C*) is the assignment of a value to each variable of *X* such that all constraints of *C* are satisfied. A constraint network is *satisfiable* iff it admits at least one solution. The *Constraint Satisfaction Problem (CSP)* is to determine whether a given constraint network is satisfiable, or not. A classical approach for solving this NP-complete problem is to perform a depth-first search with backtracking, while filtering domains after each taken decision. This procedure builds a binary search tree *T*: for each internal node *v* of *T*, a pair (*x*, *v*) is selected where *x* is a variable and *v* is a value in dom(*x*). Then, two cases are considered: the assignment *x* = *v* (positive decision) and the refutation *x* ≠ *v* (negative decision). The *future* variables of a constraint *c*, denoted by **fut**(*c*), are the variables in **scp**(*c*) at a given node of the search tree that have not been explicitly assigned by the search algorithm.

When an objective function (integer or real-valued function defined on a subset of variables of $\mathcal{X}$) is added to the constraint network, we obtain an instance of the *Constraint Optimization Problem (COP)*. Backtrack search for COP relies on an optimization strategy based on decreasingly updating the maximal bound (assuming minimization) whenever a solution is found; this is a kind of ramp-down strategy (related to Branch and Bound), whose principle is equivalent (still assuming a minimization problem) to adding a special objective constraint $\text{obj} < \infty$ to the constraint network (although it is initially trivially satisfied), and to update the limit of this constraint whenever a new solution is found.

## 3 Variable Ordering Heuristics

We provide in this section a quick overview of popular general-purpose search heuristics. The simple variable ordering heuristic dom [13], which selects variables in sequence of increasing size of domain, has long been considered as the most robust backtrack search heuristic. However, twenty years ago, *adaptive* heuristics were introduced: they take into account information collected along the part of the search space (tree) already explored.

In this paper, we shall mainly focus our attention to the popular *adaptive heuristics* based on constraint weighting (*wdeg*, *wdeg/dom*, *cacd*, *chs*), and failure rating (*frba/dom*). We will also refer to *impact*, *activity* and *counting-based heuristics*, which are defined as follows:

- **impact**, or **ibs** (Impact-Based Search), selects in priority the variable with the highest impact. The impact of a variable *x* gives a measure about the importance of *x* in reducing the search space [25]. The size of the search space of a CN *P* is the product of all current domain sizes:

$$ \text{size}(P) = \prod_{x \in \mathcal{X}} |\text{dom}(x)| $$

The impact *I* of a variable assignment *x* = *a* on *P* is computed as follows:

$$ I(x = a) = 1 - \frac{\text{size}(P')}{\text{size}(P)} $$

where $P' = \phi(P|_{x=a})$ denotes the CN obtained after assigning *x* to *a* and running the filtering process $\phi$ (e.g., enforcing arc consistency). Note that if $P'$ leads to a failure, then $I(x = a) = 1$. It is easy to see that this heuristic can be used for value selection as well.

- **activity**, or **abs** (Activity-Based Search), selects in priority the variable with the highest activity. The activity of a variable *x* is roughly measured by the number of times the domain of *x* is reduced during search [22]. This heuristic is motivated by the key role of propagation in constraint programming and relies on a decaying sum to forget the oldest statistics progressively. More formally, the activity $A(x)$ of a variable *x* is updated at
9:4 Guiding Backtrack Search by Tracking Variables During Constraint Propagation

each (new) node of the search tree (after a decision has been taken by the solver followed
by constraint propagation) regardless of the outcome (success or failure) by the following
two rules:

- $A(x) = A(x) * \gamma$, where $0 \le \gamma \le 1$ is an age decay parameter, if the domain of $x$ has not been affected (i.e., has not been reduced)

- $A(x) = A(x) + 1$ otherwise

The activities are initialized by making random probing in the search space.

- Counting-based search relies on computing the solution density of each variable-value assignment for a constraint in order to build an integrated variable-selection and value-selection heuristic [24]. Depending on the constraints, computing such information can carry a high computational cost although some mechanisms have been proposed to accelerate it [9].

Now, to introduce wdeg and wdeg/dom, we need to describe the way constraint propagation
is run each time a decision is taken by the backtrack search algorithm. Algorithm 1 describes
the constraint-oriented propagation scheme, which uses of a set of constraints for piloting
propagation. This simplifies the presentation here, whereas later, we will introduce the
new heuristic pick/dom in the context of the variable-oriented scheme. Initially, the set Q
contains the whole set of constraints of the constraint network. Then, each constraint c in
Q is picked in turn and a filtering process is applied from c: typically, this is for enforcing
arc-consistency (or a partial form) by calling filter(c) at Line 4. The call to this function
returns a subset of variables involved in c, denoted by X, whose domains have been modified
(i.e., such that at least one value has been removed from each of these domains). By means
of X, we can update Q so as to ensure constraint propagation is run until a fixed point is
reached. If ever the domain of one variable of X is empty, it simply means that a conflict
occurred (a dead-end has been identified) and so, a backtrack is required. This is triggered
by the returned Boolean value false, after having called the function incrementWeight
with the culprit constraint (responsible for the domain wipeout) passed as a parameter. In
the initial paper [3], the principle of constraint weighing is very simple: the weight of the
culprit constraint c is incremented by 1.

■ Algorithm 1 propagate((X, C): CN): Boolean.

1 $Q \leftarrow C$
2 while $Q \neq \emptyset$ do
3 &nbsp;&nbsp;&nbsp; pick and delete $c$ from $Q$
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$X \leftarrow \text{filter}(c)$ // $X$ are variables in $\text{scp}(c)$ with reduced domains
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if $\exists x \in X \mid \text{dom}(x) = \emptyset$ then
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{incrementWeight}(c)$
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\text{return false}$ // detected inconsistency
8 &nbsp;&nbsp;&nbsp;foreach $c' \in C \mid c' \neq c \land X \cap \text{scp}(c') \neq \emptyset$ do
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$Q \leftarrow Q \cup \{c'\}$
10 return true

To summarize, each constraint *c* admits a weight, initially set to 1, which is incremented whenever a domain wipeout occurs while filtering *c*. Importantly, it was observed experimentally that it was more effective to consider only the future variables involved in a culprit constraint. Technically, instead of associating a global weight *c.weight* with each constraint
c, one can introduce a local weight `c.weight[x]` to be associated with each variable `x` in `scp(c)`. Hence, when a conflict occurs, instead of incrementing the weight `c.weight` of the culprit constraint, one can decide to increment the local weight `c.weight[x]` of each future variable involved in `scp(c)`.

The heuristics `wdeg` and `wdeg/dom` are defined as follows:

- `wdeg` selects in priority the future variable with the highest "weighted degree". Each variable `x` is given a weighted degree, which is the sum of the weights over all constraints involving `x` and at least another future variable. For each future variable `x`, the score of `x` according to `wdeg` is:

$$ \sum_{c \in C: x \in scp(c) \land |\text{ fut}(c)| > 1} c.\text{weight}[x] $$

- `wdeg/dom` selects in priority the future variable with the highest ratio "weighted degree to current domain size". For each future variable `x`, the score of `x` according to `wdeg/dom` is:

$$ \frac{\sum_{c \in C: x \in scp(c) \land |\text{ fut}(c)| > 1} c.\text{weight}[x]}{|\text{dom}(x)|} $$

In classical forms of `wdeg` and `wdeg/dom`, counters are incremented by 1., which remains very simplistic and does not differentiate between constraints. This is why constraint weighting, in the so-called variant `cacd` [28], has been refined by exploiting as information both the "current arity" of the culprit constraint (i.e., the number of future variables) and the size of the current domains of the future variables. The increment values computed for the classical and `cacd` variants are precisely shown in Algorithm 2.

**Algorithm 2** incrementWeight(c: Constraint).

1. **foreach** $x \in \text{fut}(c)$ **do**
2. &nbsp;&nbsp;&nbsp;&nbsp;**if** variant `cacd` **then**
3. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;`c weights[x]` $\leftarrow$ `c weights[x]` + $\frac{1}{|\text{fut}(c)| \times (1 + |\text{dom}(x)|)}$
4. &nbsp;&nbsp;&nbsp;&nbsp;**else**
5. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;`c weights[x]` $\leftarrow$ `c weights[x]` + 1

Note that to break ties, which correspond to sets of variables that are considered as equivalent by the heuristic, one can use a second criterion. However, for adaptive heuristics, it is usual that the first encountered variable with the best score is selected.

Finally, we introduce two recent heuristics: the former, `chs` [12], exploits the history of search failures, while the latter, `frba/dom` [20], computes the proportion of failing assignments for each variable.

- `chs` (Conflict-History Search), selects in priority variables appearing in recent failures. All failures are registered with a timestamp. More precisely, `chs` maintains for each constraint `c`, a score `q(c)` and updates it at every domain wipeout with an exponential recency weighted average:

$$ q(c) = (1 - \alpha) \times q(c) + \alpha \times r(c) $$

where $\alpha = 0.4$ (decreasing as time goes by) and $r(c)$ is the reward given when a domain wipeout occurred. The reward is higher when the constraint frequently enters in conflict:

$$ r(c) = \frac{1}{\#\text{Conflicts} - \text{Conflict}(c) + 1} $$
9:6 Guiding Backtrack Search by Tracking Variables During Constraint Propagation

where #Conflicts is the total number of conflicts and Conflict(c) stores the last value of #Conflicts when c led to a failure. The conflict history score of a variable x which will be used in selecting the branching variable is given by:

$$
\frac{\sum_{c \in C: x \in \text{scp}(c) \land |\text{ fut}(c)| > 1} q(c) + \delta}{|\text{dom}(x)|}
$$

where δ is a positive real number close to 0 that avoids random selection at the beginning
of search.

- **frba/dom (Failure Rate Based Search)**, selects in priority the variables that most often lead to a conflict when assigning them. For this purpose, two counters are associated with each variable *x*: the former, **#Conflicts**(*x*), records the number of times a failure has been observed just by propagating an assignment involving *x*, and the latter, **#Assigns**(*x*), the number of times the variable *x* was assigned. The failure rate of a variable *x* is then:

$$
fr(x) = \frac{\#\text{Conflicts}(x)}{\#\text{Assigns}(x)}
$$

In addition, similarly to the factor used for chs, we compute:

$$
a(x) = \frac{1}{\#\text{Conflicts} - \text{Conflict}(x) + 1}
$$

where #Conflicts is the total number of conflicts and Conflict(x) stores the last value
of #Conflicts when x being assigned led to a failure.

The failure rate score of a variable *x* by *frba/dom* is then:

$$
\frac{fr(x) + a(x)}{|dom(x)|}
$$

4 Variable Tracking in Conflicting Propagation

In this section, we introduce the principle of Variable Tracking in Conflicting Propagation (VTCP). More specifically, we introduce a new variable ordering heuristic called **pick/dom**, whose principle is to track the variables that are used to trigger filtering operations during constraint propagation. However, it is important to note that this tracking is only used to update counters (called “pick degrees”) associated with variables when constraint propagation ends with a conflict. We show how to implement such variable tracking within the variable-oriented propagation scheme.

To record information about tracked variables, we just need to associate a counter pick[x]
with each variable x of the CN. Initially, this counter (“pick degree”) is set to 0. According
to the selected mode (see below) used to update these counters, recorded values may be in
real or integer forms. The heuristic pick/dom is simply defined as follows:

- **pick/dom** selects in priority the future variable with the highest ratio “pick degree to current domain size”. For each future variable *x*, the score of *x* according to **pickdom** is:

$$
\frac{\text{pick}[x]}{|\text{dom}(x)|}
$$

In the rest of this section, we show how pick degrees are computed.

In ACE, constraint propagation follows the variable-oriented scheme (as initially introduced in [21]): the set *Q* contains variables. The principle is that whenever a value is removed
from the domain of a variable $x$, this variable is added to $Q$. In a first step, by ignoring any statement related to $L$ and $\Delta_x$, we can rather easily recognize the variable-oriented propagation scheme in Algorithm 3: as long as there is a variable in $Q$, one of them, $x$ is picked, and we execute the filtering algorithms (propagators) attached to all constraints involving $x$, while updating $Q$ when necessary.

■ Algorithm 3 propagate((X, C): CN): Boolean.

1 L = <>
2 Q ← X
3 while Q ≠ ∅ do
4 pick and delete x from Q
5 Δ<sub>x</sub> ← 0
6 for c ∈ C | x ∈ scp(c) do
7 X ← filter(c, Δ<sub>x</sub>) // Δ<sub>x</sub> is updated during call
8 if ∃y ∈ X | dom(y) = ∅ then
9 append (x, Δ<sub>x</sub>) to L
10 incrementPick(L)
11 return false // detected inconsistency
12 Q ← Q ∪ X
13 if Δ<sub>x</sub> > 0 then
14 append (x, Δ<sub>x</sub>) to L
15 return true

■ Algorithm 4 incrementPick(L: Sequence).

1 foreach *i* ranging from 1 to |L| do
2   (x_i, Δ_{x_i}) ← L[i]
3   switch VARIANT do
4     case 0 do
5       increment ← 1
6     case 1 do
7       increment ← Δ_{x_i}
8     case 2 do
9       increment ← 100 × (Δ_{x_i} / (Σ_{j=1}^{|L|} Δ_{x_j}))
10    case 3 do
11     increment ← (n−depth/n) × 100 × (Δ_{x_i} / (Σ_{j=1}^{|L|} Δ_{x_j}))
12   pick[x_i] ← pick[x_i] + increment

Now, let us consider first the structure *L*, which is a list, initially empty, keeping track of
any variable *x* that plays a role (i.e., triggers some effective filtering) during propagation,
together with an information indicating the degree Δ*x* of this role. Notice that, since *L* is
a list, the same variable may occur several times. Concerning the local variable Δ*x*, it is
initialized to 0 in Line 5. When the loop starting at Line 6 ends, Δ*x* indicates how many
9:8 Guiding Backtrack Search by Tracking Variables During Constraint Propagation

values have been deleted from the domain of *x* in the different calls to Function *filter* at Line 7. In practice, it is possible to handle $\Delta_x$ in a non-intrusive way by introducing a global variable whose value is incremented whenever a value is deleted (whatever the domain is). If at Line 13, the value of $\Delta_x$ is 0, it means that no filtering/reduction was performed at all since the time *x* was picked. This is then a useless “pick”, which is the reason why we do not update the structure *L* at Line 14. Importantly, the list *L* is only exploited if Algorithm 3 returns *false* (because a conflict is detected). Before returning *false*, the last picked variable is added to *L* (because we have the guarantee of some filtering) and the function *incrementPick* is called in order to update some picked degrees.

The way picked degrees are updated is shown in Algorithm 4. Four modes (denoted by values ranging from 0 to 3) are possible. In mode 0, the picked degree of any occurrence of a variable present in *L* is incremented by 1. In mode 1, the increment is given by $\Delta_x$, the impact of *x* after having been picked. In mode 2, each time constraint propagation is run, 100 points are shared according to the relative impacts of the variables present in *L*. In mode 3, a coefficient is applied to 100, depending on the current depth of the solver. As a first extreme case, the current depth is 0, which means that we are at the root node, and so, 100 points are spread. As a second extreme case, the current depth is *n*, meaning that we are a leaf, and 0 point is shared. The rationale is that we give more importance to nodes near the top of the search tree.

5 Experimental Results

In our experiments, we have compared general-purpose variable ordering heuristics based on constraint weighting, failure rating and variable tracking during conflicting propagation, with the constraint solver ACE [15]. More specifically, we have compared the four variants of *pick/dom* with *wdeg-cacd* [28], *wdeg/dom* [3] and *chs* [12], as well as the recently introduced *frba/dom*. From now on, for simplicity, *pick/dom*, *wdeg-cacd* and *frba/dom* will be referred as *pick* (while appending the mode in subscript text), *cacd* and *frba*, respectively. Note that *ibs* and *abs* are not retained in our experiments because they are usually outperformed when used in ACE. Concerning the value ordering heuristic, it systematically chooses the smallest value in domains.

We have considered two benchmarks, denoted by *xcsp-csp* and *xcsp-cop*, which are respectively composed of all CSP and COP instances selected for the main tracks of the XCSP³ competitions (In 2019, instances were randomly selected from existing series, and there were no competitions held in 2020 and 2021) organized in 2017, 2018 [16] and 2022 [2] (most of them generated by the Python library PyCSP³ [19]). They correspond to two full sets of 942 and 1,034 instances in format XCSP³ [4], for exactly 77 and 50 problems, respectively. A time limit of 1,200 seconds was given per instance.

**Ranking.** Results will be partly analyzed from the scoring function used for the 2022 XCSP³ competition. For self-containedness, we recall it now. The number of points won by a solver *S* is decided as follows:

- for CSP, this is the number of times *S* is able to solve an instance, i.e., to decide the satisfiability of an instance (either exhibiting a solution, or indicating that the instance is unsatisfiable)

- for COP, this is, roughly speaking, the number of times *S* gives the best known result, compared to its competitors. More specifically, for each instance *I*:
- if *I* is unsatisfiable, 1 point is won by *S* if *S* indicates that the instance *I* is unsatisfiable, 0 otherwise,

- if *S* provides a solution whose bound is less good than another one (found by another competing solver), 0 point is won by *S*,

- if *S* provides an optimal solution, while indicating that it is indeed the optimality, 1 point is won by *S*,

- if *S* provides (a solution with) the best found bound among all competitors, this being possibly shared by some other solver(s), while indicating no information about optimality: 1 point is won by *S* if no other solver proved that this bound was optimal, 0.5 otherwise.

## 5.1 Global Overview of Results

We start this experimental section with a global overview of the obtained results. The scores of tested heuristics on *xcsp-csp* and *xcsp-cop* are given in Table 1 and Table 2. First, let us make some comments on CSP. Here, the default version of ACE is the best one (cacd). The heuristics are relatively close, and the differences mainly come from the number of (solved) SAT instances. On our benchmark, the heuristic *frba* appears to be the worst one. Concerning the *pick* variants, they are quite close, even if *pick*<sub>3</sub> is the one that is most often the fastest heuristic.

<table><caption>■ Table 1 Ranking on *xcsp-csp*. For each heuristic, we report the number of SAT/UNSAT instances, the total number of solved instances and the number of times a heuristic is the fastest for solving an instance. We also report the result for the virtual best solver/heuristic (VBS).</caption><thead><tr><th>Heuristic</th><th>Solved</th><th>SAT</th><th>UNSAT</th><th>Fastest</th></tr></thead><tbody><tr><td>VBS</td><td>676</td><td>481</td><td>195</td><td>676</td></tr><tr><td><b>cacd</b></td><td><b>646</b></td><td><b>457</b></td><td><b>189</b></td><td>411</td></tr><tr><td><b>chs</b></td><td>636</td><td>449</td><td>187</td><td>390</td></tr><tr><td><b>pick<sub>3</sub></b></td><td>632</td><td>441</td><td><b>189</b></td><td><b>442</b></td></tr><tr><td><b>pick<sub>0</sub></b></td><td>631</td><td>442</td><td>188</td><td>439</td></tr><tr><td><b>pick<sub>1</sub></b></td><td>630</td><td>441</td><td><b>189</b></td><td>427</td></tr><tr><td><b>wdeg/dom</b></td><td>626</td><td>442</td><td>183</td><td>377</td></tr><tr><td><b>pick<sub>2</sub></b></td><td>625</td><td>437</td><td>188</td><td>434</td></tr><tr><td><b>frba</b></td><td>618</td><td>431</td><td>187</td><td>391</td></tr></tbody></table>

For COP, all *pick* variants are the most efficient heuristics: indeed, the gap with the other heuristics is significant. The variant *pick*<sub>3</sub> is the best one in terms of the number of proved optima and the number of found best bounds. The heuristics based on constraint weighting got quite similar results, and even if the heuristic *frba* appears to be also outperformed, it is important to report that this heuristic is very good in proving optima.

In order to provide a deeper analysis of the results of our experiments, we propose, in the next section, to focus only on the following three heuristics:

- *cacd* as a robust representative of the heuristics based on constraint weighting. The results of *cacd*, *wdeg/dom* and *chs* are quite similar on COP but *cacd* appears to be the most efficient heuristic on CSP. In addition, this is the default search strategy of ACE.

- *frba* as a recent proposed heuristic based on failure rating.
9:10 Guiding Backtrack Search by Tracking Variables During Constraint Propagation

■ Table 2 Ranking on *xcsp-cop*. For each heuristic, we report the number of proved optima, the number of times the best bound has been found and the score as computed for the 2022 XCSP³ competition. We also report the results for the virtual best solver/heuristic (VBS).

<table><thead><tr><th>Heuristic</th><th>Score</th><th>Optimum</th><th>Best Bound</th></tr></thead><tbody><tr><td>VBS</td><td>959.00</td><td>372</td><td>956</td></tr><tr><td><b>pick<sub>3</sub></b></td><td><b>624.50</b></td><td><b>347</b></td><td><b>627</b></td></tr><tr><td><b>pick<sub>2</sub></b></td><td>618.00</td><td>343</td><td>621</td></tr><tr><td><b>pick<sub>1</sub></b></td><td>617.50</td><td>336</td><td>621</td></tr><tr><td><b>pick<sub>0</sub></b></td><td>612.50</td><td>341</td><td>617</td></tr><tr><td><b>wdeg/dom</b></td><td>569.50</td><td>312</td><td>583</td></tr><tr><td><b>cacd</b></td><td>557.50</td><td>315</td><td>570</td></tr><tr><td><b>frba</b></td><td>557.00</td><td>341</td><td>560</td></tr><tr><td><b>chs</b></td><td>554.50</td><td>324</td><td>563</td></tr></tbody></table>

■ *pick*<sub>3</sub> as the best variant of variable tracking in conflict propagation. Actually, it is the best heuristic on COP, able to find, most often, best bounds and also to prove them. On CSP, it is also the most efficient variant of VTCP and the fastest one.

5.2 Comparing Best Heuristics

In this section, we compare the three selected heuristics, namely, *cacd*, *frba* and *pick*<sub>3</sub>. We start this comparison on the *xcsp-csp* benchmark with Figure 1, which shows some classical scatter plots (permitting to compare two algorithms with rather good precision).

<figure><img src="image_73.png" /><figcaption>■ Figure 1 Scatter plots for CSP instances. Each dot represents an instance: its value on the *x*-axis (resp. *y*-axis) represents the time needed by the heuristic labelling the *x*-axis (resp. *y*-axis) to solve it. Blue (resp. Orange) dots corresponds to SAT (resp. UNSAT) instances. The dots below the diagonal represent then the instances where the *y*-axis heuristic is faster than the *x*-axis one.</figcaption></figure>

Here, it is clear that *frba* is less efficient than *cacd* and *pick*<sub>3</sub>. Even if *cacd* is the best heuristic on CSP (see Table 1), notably on the hardest instances, an instance-by-instance comparison between *cacd* and *pick*<sub>3</sub> looks less obvious, as the dots are uniformly distributed over both parts separated by the diagonal.

As an alternative for visualizing these results, a Venn diagram is depicted in Figure 4a. In such a diagram, each circle represents the instances solved by a heuristic. An overlapping region represents a set of instances solved in an *equivalent manner* by two heuristics, or
three in the case of the central region. For CSP, two heuristics are equivalent if they require the same amount of time (with a tolerance of one second) to find the same result (SAT or UNSAT). A region with no overlap emphasizes the instances that are better solved by a single heuristic. Here, the ranking is clear: *pick*₃ is the winner, followed by *cacd* and finally *frba*. To summarize the results for CSP: although *cacd* is the most robust heuristic when considering the number of solved instance (with a timeout set to 1,200 seconds), *pick*₃ is usually the fastest heuristic.

Next, we consider the *xcsp-cop* benchmark. A first analysis can be made from the scatter plots in Figure 2. Each plot is based only on the instances whose optimality has been proved by at least one of the two compared heuristics. On the one hand, one can see in Figure 2a and Figure 2b that *frba* and *pick*₃ are better at proving optimality than *cacd*. On the other hand, Figure 2c does not show any dominance between *frba* and *pick*₃.

<figure><img src="image_74.png" /><figcaption>Figure 2 Scatter plots for COP instances. Each dot represents an instance: its value on the *x*-axis (resp. *y*-axis) represents the time needed by the heuristic labelling the *x*-axis (resp. *y*-axis) to prove its optimum. The dots below the diagonal represent then the instances where the *y*-axis heuristic is faster than the *x*-axis one.</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>Figure 3 Plots for COP Instances. Above each figure, the number of instances (partially) solved by at least one of the two heuristics as well as the total number of instances is indicated. When a heuristic is better than the other, this contributes to its area, representing the difference in resolution time. The *vol.* value indicates the volume of this surface, the *perc.* value indicates the ratio of the volume to the sum of the two volumes. The difference between two heuristics takes into account the best bounds found and possibly the proof of optimality.</figcaption></figure>

Figure 3 allows a complementary analysis. Each graphics shows a pairwise comparison with two areas. When a heuristic is better than the other, the difference in solving time is
computed. An heuristic *a* is better than a heuristic *b* if it proves optimality whereas *b* does not, or if it finds a better bound than *b*, or, finally, if they find the same bound and it is faster than *b*. In the two first cases, the resolution time of *b* is set to the time limit, namely 1,200. Each computed time difference contributes to the area of the best heuristic. Hence, a larger area means that the attached heuristic offers better performances than the other one. Each area is annotated with *vol*. and *perc*. which respectively denote the volume of the area and the ratio, in percentage, of the volume to the sum of the two volumes. Finally, strictly equivalent instances are indicated between the 2 vertical bars. In Figure 3a, one can observe that *cacd* is slightly more robust than *frba*. Interestingly, in Figure 3b and Figure 3c, the competitiveness of *pick*<sub>3</sub> is clearly visible: the area of *pick*<sub>3</sub> is almost twice as large as that of *cacd* or *frba*.

<figure><img src="image_76.png" /><figcaption>(a) Venn Diagram on CSP Instances.</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>(b) Venn Diagram on COP Instances.</figcaption></figure>

■ **Figure 4** Venn diagrams for *cacd*, *frba* and *pick*<sub>3</sub> on CSP and COP instances. Each circle represents the instances solved by a heuristic. An overlapping region represents a set of instances solved in an equivalent manner by two heuristics, or three in the case of the central region. Two heuristics are considered as being equivalent if they found the same result in the same amount of time. For CSP, this is the time for proving (un)satisfiability. For COP, this is for getting the best bound or proving optimality. A region not overlapped emphasizes the instances that are better solved by a single heuristic.

These results are confirmed by the Venn diagram in Figure 4b. The circles, and especially
the one for pick₃, stand out from the centre. This Venn diagram also suggests that each
heuristic would be more suitable for certain problems. This is confirmed in Table 3. First,
frba performs particularly well on five problems, namely Cutstock, DC, ItemsetMining
OpenStacks and TravelingSalesman. Indeed, it is able to close more instances and provides
better bounds than cacd and pick₃. In turn, cacd behaves better on six problems: Cyc-
licBandwith, Ramsey, StillLife and Taillard, and to a lesser extent PseudoBoolean and
QueenAttacking. As for pick₃, the set of problems where it dominates is clearly larger:
Auction, BinPacking, ChessBoardColoration, CoinsGrid, EchelonStock2, FAPP, GraphColor-
ing, MultiAgentPathFinding, NurseRostering, OPD, QuadraticAssignment, RCPSP, RLFAP,
Spot, SteelMillSlab, SumColoring TAL, TravelingTournament, Triangular and Warehouse.
This is just under half of the problems. For some problems the gap is quite large, namely
BinPacking, FAPP or RCPSP.

Finally, we show in Table 4 some details for some chosen (representative and singular)
instances of some problems. This may be helpful for testing and reproducing the results we have obtained.

<table><caption>**Table 3** Details per problem (COP). For each problem, the number of instances is displayed, and for each heuristic, the couple 'number of proved optima : number of best bounds found' is provided. Best results are highlighted (when at least one heuristic is outperformed by the other(s)).</caption><thead><tr><th>Problem</th><th>cacd</th><th>frba</th><th>pick<sub>3</sub></th><th>Problem</th><th>cacd</th><th>frba</th><th>pick<sub>3</sub></th></tr></thead><tbody><tr><td>AirCraft (13)</td><td>4:10</td><td>5:6</td><td>4:10</td><td>NursingWork (12)</td><td>1:2</td><td>1:4</td><td>1:4</td></tr><tr><td>Auction (16)</td><td>2:6</td><td>2:4</td><td>2:14</td><td>OPD (17)</td><td>8:9</td><td>9:14</td><td>9:15</td></tr><tr><td>BACP (24)</td><td>4:24</td><td>4:24</td><td>4:24</td><td>OpenStacks (16)</td><td>12:16</td><td>14:16</td><td>12:16</td></tr><tr><td>BinPacking (51)</td><td>1:17</td><td>0:25</td><td>2:33</td><td>PeacableA (14)</td><td>4:9</td><td>4:9</td><td>4:8</td></tr><tr><td>BusScheduling (10)</td><td>1:5</td><td>1:5</td><td>1:8</td><td>PizzaVoucher (10)</td><td>4:7</td><td>6:8</td><td>5:10</td></tr><tr><td>CVRP (10)</td><td>0:2</td><td>0:2</td><td>0:8</td><td>PseudoB (30)</td><td>13:23</td><td>13:18</td><td>13:22</td></tr><tr><td>ChessBoard (17)</td><td>3:13</td><td>2:13</td><td>3:14</td><td>Quadratic (36)</td><td>6:26</td><td>6:14</td><td>6:27</td></tr><tr><td>ClockTriplet (10)</td><td>2:9</td><td>2:6</td><td>2:9</td><td>QueenAtt (17)</td><td>7:11</td><td>5:3</td><td>8:6</td></tr><tr><td>CoinsGrid (10)</td><td>2:7</td><td>2:4</td><td>3:9</td><td>Rack (4)</td><td>2:4</td><td>3:4</td><td>3:4</td></tr><tr><td>CrosswordDes (13)</td><td>3:5</td><td>3:3</td><td>3:5</td><td>Ramsey (17)</td><td>5:17</td><td>4:16</td><td>5:14</td></tr><tr><td>CutStock (17)</td><td>6:9</td><td>8:9</td><td>7:9</td><td>RCPSP (43)</td><td>31:34</td><td>34:35</td><td>35:43</td></tr><tr><td>CyclicBand (12)</td><td>4:9</td><td>2:3</td><td>2:6</td><td>RLFAP (50)</td><td>10:22</td><td>12:32</td><td>12:48</td></tr><tr><td>DC (26)</td><td>7:14</td><td>10:25</td><td>12:18</td><td>Spot (10)</td><td>2:7</td><td>2:6</td><td>2:8</td></tr><tr><td>EchelonStock2 (10)</td><td>0:9</td><td>0:7</td><td>0:10</td><td>SteelMill (17)</td><td>2:2</td><td>2:2</td><td>2:7</td></tr><tr><td>FAPP (18)</td><td>2:2</td><td>2:7</td><td>3:13</td><td>StillLife (47)</td><td>15:46</td><td>14:21</td><td>13:29</td></tr><tr><td>FastFood (17)</td><td>17:17</td><td>17:17</td><td>17:17</td><td>SumColoring (14)</td><td>4:8</td><td>4:7</td><td>4:10</td></tr><tr><td>Filters (8)</td><td>8:8</td><td>8:8</td><td>8:8</td><td>Taillard (51)</td><td>5:17</td><td>5:6</td><td>5:13</td></tr><tr><td>GolombRuler (47)</td><td>19:31</td><td>18:27</td><td>18:35</td><td>TAL (20)</td><td>9:16</td><td>10:15</td><td>10:18</td></tr><tr><td>GraphCol (28)</td><td>14:24</td><td>15:23</td><td>17:24</td><td>TemplateDes (15)</td><td>12:13</td><td>12:13</td><td>11:12</td></tr><tr><td>ItemsetMining (15)</td><td>3:10</td><td>7:12</td><td>6:10</td><td>TravelingTour (14)</td><td>2:7</td><td>2:2</td><td>2:10</td></tr><tr><td>Knapsack (31)</td><td>22:31</td><td>26:31</td><td>26:31</td><td>TravelingSale (29)</td><td>3:16</td><td>3:17</td><td>3:15</td></tr><tr><td>LowAuto (31)</td><td>12:21</td><td>12:23</td><td>14:19</td><td>Triangular (10)</td><td>0:3</td><td>1:5</td><td>1:9</td></tr><tr><td>Mario (10)</td><td>10:10</td><td>10:10</td><td>10:10</td><td>VRP (17)</td><td>0:9</td><td>0:1</td><td>0:9</td></tr><tr><td>MultiAgentP (20)</td><td>4:6</td><td>8:9</td><td>9:17</td><td>WarOrPeace (10)</td><td>6:10</td><td>6:10</td><td>6:10</td></tr><tr><td>NurseRost (41)</td><td>2:16</td><td>3:3</td><td>3:29</td><td>WareHouse (9)</td><td>0:3</td><td>0:0</td><td>0:9</td></tr></tbody></table>

# 6 Discussion

The three different ways of exploiting conflicts for guiding search, as experimented in the last section, show somewhat complementary behaviors. This can be explained by the fact that information is extracted at different moments: at the very beginning of the process conducting to a conflict (i.e., at the time of the decision), during constraint propagation, or at the time the last propagator (filtering algorithm) is solicited. One can then refer to such approaches as *early* (E), *midway* (M) and *late* (L) operational treatment of conflicts. This is illustrated in Figure 5 where a new decision $x = a$ is taken, when solving a CN $P$, in the continuity of two previously taken decisions $v = a$ and $w \neq b$. In our scenario, running constraint propagation $\phi$ on (the current state of) $P$ after having assigned the value $a$ to $x$, i.e., $\phi(P|x=a)$, reveals a new conflicting (dead-end) situation (denoted by $\perp$). The early processing of this new conflict consists in considering the variable $x$ involved in the decision as the main culprit. This is the principle behind the heuristic `frba/dom`. The *midway* processing of this conflict consists in considering all variables having played a role (i.e., having been picked) during propagation as having contributed to the failure. This is the principle underlying the heuristic `pick/dom`. The late processing of this conflict consists in considering the last constraint (here, $c_w^2$) provoking a domain-wipeout (i.e., removing the last value of a domain) as the object of interest. This is the principle of constraint weighting, as in `wdeg/dom`.

One related heuristic, which is based on some form of *midway* strategy, is `abs`. However, VTCP and `abs` collect information according to different strategies. `abs` updates activity
<figure><img src="image_78.png" /><figcaption>Figure 5 Illustration of pivotal moments for collecting information about conflicts: this corresponds to early (E), midway (M) and late (L) processing of conflicts.</figcaption></figure>

counters at each node, whereas VTCP only considers conflicting nodes. **abs** systematically updates the activity counters of each variable, whereas VTCP only increments the counters of variables at the origin of calls to effective propagators.

Importantly, we do believe that the experimental results we have obtained are significant, for several reasons. First, they can be reproduced in an open-source constraint solver (with a repository available in GitHub). Second, the number of models and instances used for our experiments is very large, involving more than 120 problems of various nature. Third, ACE is a competitive constraint solver and (although not officially engaged) showed good performances in (notably the COP main track) of the 2022 XCSP³ competition.

Finally, it is true that the importance of variable tracking in conflicting propagation looks more limited when solving CSP instances. Actually, solving a (satisfiable) CSP instance involves one single phase: finding a solution, whereas solving a COP instance involves two subsequent phases: moving down towards an optimal solution, and proving optimality. We believe that **pick/dom** is rather efficient for the first phase of COP solving (for the second phase, a learning mechanism like in Picat [29] or OR-Tools becomes central).

## 7 Conclusion

In this paper, we have introduced a new way of exploiting conflicts during backtrack search so as to build a well-informed variable ordering heuristic. In contrast to existing heuristics relying on the early and late treatments of failing nodes, this new heuristic, **pick/dom**, consists in a midway processing of conflicts by tracking variables during constraint propagation. The robustness of **pick/dom** has been demonstrated from a vast experimentation campaign involving more than 120 problems (and around 2,000 instances). Interestingly, the three different forms of exploiting conflicts, based on different moments when to collect information, entail somewhat complementary behaviours of the solver. This opens some perspectives for building a still more robust solver by combining these conflict-based heuristics in a clever way. Identifying features or properties of problem instances (e.g., tree width, backbone size, presence of strong communities, structure of variable arrays, etc.) which are favorable for a certain form of conflict processing is an issue which would deserve to be addressed.
**Table 4** Details per instance. For each selected instance and each heuristic, the best bound, the time (in seconds) to find it (TB) and the time (in seconds) to prove the optimum (TO, if proven) are provided. Best results are highlighted. DC-A, DC-B, DC-C and DC-D stands for DC-midori-xor4-d1-t0-r05-v128, DC-midori-xor4-d1-t0-r08-v128, DC-skinny-xor0-d0-t0-r09-v64-z0, and DC-midori-xor4-d1-t0-r09-v128, respectively. Some other names have been shortened for more visibility.

<table><thead><tr><th rowspan="2">Instance</th><th colspan="3">cacd</th><th colspan="3">frba</th><th colspan="3">pick<sub>3</sub></th></tr><tr><th>Bound</th><th>TB</th><th>TO</th><th>Bound</th><th>TB</th><th>TO</th><th>Bound</th><th>TB</th><th>TO</th></tr></thead><tbody><tr><td>DC-A</td><td>5</td><td>733</td><td>739</td><td>5</td><td>9</td><td>23</td><td>5</td><td>21</td><td>28</td></tr><tr><td>DC-B</td><td>41</td><td>59</td><td>-</td><td>8</td><td>174</td><td>-</td><td>8</td><td>316</td><td>484</td></tr><tr><td>DC-C</td><td>41</td><td>5</td><td>-</td><td>41</td><td>9</td><td>-</td><td>41</td><td>3</td><td>1,087</td></tr><tr><td>DC-D</td><td>53</td><td>753</td><td>-</td><td>9</td><td>699</td><td>-</td><td>43</td><td>179</td><td>-</td></tr><tr><td>Fapp-m2s-21-0500</td><td>180,499K</td><td>1,195</td><td>-</td><td>63,832K</td><td>870</td><td>881</td><td>177,187K</td><td>1,008</td><td>-</td></tr><tr><td>Fapp-m2s-02-0250</td><td>221,762K</td><td>1,191</td><td>-</td><td>221,520K</td><td>1,158</td><td>-</td><td>221,591K</td><td>1,197</td><td>-</td></tr><tr><td>Fapp-m2s-test03-0400</td><td>453,213K</td><td>1,199</td><td>-</td><td>449,289K</td><td>1,192</td><td>-</td><td>450,098K</td><td>1,200</td><td>-</td></tr><tr><td>QueenAttacking-09</td><td>0</td><td>94</td><td>98</td><td>1</td><td>119</td><td>-</td><td>0</td><td>439</td><td>442</td></tr><tr><td>QueenAttacking-11</td><td>2</td><td>245</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1</td><td>708</td><td>736</td></tr><tr><td>QueenAttacking-13</td><td>11</td><td>158</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>StillLife-11-14</td><td>81</td><td>264</td><td>-</td><td>79</td><td>652</td><td>-</td><td>81</td><td>625</td><td>-</td></tr><tr><td>StillLife-wastage-12</td><td>76</td><td>12</td><td>872</td><td>76</td><td>23</td><td>-</td><td>76</td><td>549</td><td>-</td></tr><tr><td>StillLife-wastage-37</td><td>681</td><td>1,189</td><td>-</td><td>619</td><td>1,182</td><td>-</td><td>585</td><td>398</td><td>-</td></tr><tr><td>Auction-cnt-d100</td><td>829K</td><td>200</td><td>-</td><td>825K</td><td>1,099</td><td>-</td><td>849K</td><td>1,197</td><td>-</td></tr><tr><td>Auction-sum-d100</td><td>840K</td><td>69</td><td>-</td><td>824K</td><td>77</td><td>-</td><td>854K</td><td>102</td><td>-</td></tr><tr><td>Auction-sum-d500</td><td>3,368K</td><td>186</td><td>-</td><td>3,264K</td><td>224</td><td>-</td><td>3,344K</td><td>89</td><td>-</td></tr><tr><td>CVRP-A-n32-k5</td><td>1,095</td><td>15</td><td>-</td><td>1,006</td><td>257</td><td>-</td><td>835</td><td>492</td><td>-</td></tr><tr><td>CVRP-A-n36-k5</td><td>1,050</td><td>77</td><td>-</td><td>1,039</td><td>1,112</td><td>-</td><td>892</td><td>74</td><td>-</td></tr><tr><td>CVRP-A-n34-k5</td><td>-</td><td>-</td><td>-</td><td>915</td><td>57</td><td>-</td><td>813</td><td>385</td><td>-</td></tr><tr><td>NurseRostering-01</td><td>607</td><td>1</td><td>-</td><td>607</td><td>1</td><td>11</td><td>607</td><td>1</td><td>10</td></tr><tr><td>NurseRostering-02</td><td>1,024</td><td>5</td><td>-</td><td>928</td><td>62</td><td>-</td><td>833</td><td>383</td><td>-</td></tr><tr><td>NurseRostering-19</td><td>68,621</td><td>1,200</td><td>-</td><td>-</td><td>-</td><td>-</td><td>60,805</td><td>1,197</td><td>-</td></tr><tr><td>Rlfap-graph-05-opt</td><td>2,882</td><td>883</td><td>-</td><td>221</td><td>2</td><td>10</td><td>221</td><td>62</td><td>70</td></tr><tr><td>Rlfap-graph-06-opt</td><td>46,647</td><td>1,198</td><td>-</td><td>34,004</td><td>1,198</td><td>-</td><td>8,882</td><td>769</td><td>-</td></tr><tr><td>Rlfap-scen-06-opt</td><td>11,211</td><td>1,166</td><td>-</td><td>10,102</td><td>491</td><td>-</td><td>3,389</td><td>124</td><td>-</td></tr><tr><td>Triangular-10</td><td>20</td><td>89</td><td>-</td><td>20</td><td>4</td><td>138</td><td>20</td><td>0</td><td>142</td></tr><tr><td>Triangular-22</td><td>50</td><td>46</td><td>-</td><td>50</td><td>161</td><td>-</td><td>52</td><td>167</td><td>-</td></tr><tr><td>Triangular-38</td><td>95</td><td>10</td><td>-</td><td>97</td><td>224</td><td>-</td><td>97</td><td>337</td><td>-</td></tr><tr><td>SteelMillSlab-m2s-3-0</td><td>68</td><td>1,059</td><td>-</td><td>-</td><td>-</td><td>-</td><td>43</td><td>995</td><td>-</td></tr><tr><td>SteelMillSlab-m2s-3-2</td><td>306</td><td>794</td><td>-</td><td>-</td><td>-</td><td>-</td><td>171</td><td>1,177</td><td>-</td></tr><tr><td>SteelMillSlab-m2s-4-0</td><td>161</td><td>1,109</td><td>-</td><td>-</td><td>-</td><td>-</td><td>94</td><td>1,186</td><td>-</td></tr></tbody></table>
References

1 K.R. Apt. *Principles of Constraint Programming*. Cambridge University Press, 2003.

2 G. Audemard, C. Lecoutre, and E. Lonca. Proceedings of the 2022 XCSP3 competition. *CoRR*, abs/2209.00917, 2022. doi:10.48550/arXiv.2209.00917.

3 F. Boussemart, F. Hemery, C. Lecoutre, and L. Sais. Boosting systematic search by weighting constraints. In *Proceedings of ECAI'04*, pages 146–150, 2004.

4 F. Boussemart, C. Lecoutre, G. Audemard, and C. Piette. XCSP3: an integrated format for benchmarking combinatorial constrained problems. CoRR, abs/1611.03398, 2016. URL: http://arxiv.org/abs/1611.03398, arXiv:1611.03398.

5 R. Dechter. *Constraint processing*. Morgan Kaufmann, 2003.

6 R. Dechter and J. Pearl. Network-based heuristics for constraint satisfaction problems. *Artificial Intelligence*, 34(1):1–38, 1988.

7 E. Demirovic, G. Chu, and P. Stuckey. Solution-based phase saving for CP: A value-selection heuristic to simulate local search behavior in complete solvers. In *Proceedings of CP’18*, pages 99–108, 2018.

8 J.-G. Fages and C. Prud’homme. Making the first solution good! In *Proceedings of ICTAI’17*, pages 1073–1077, 2017.

9 S. Gagnon and G. Pesant. Accelerating counting-based search. In *Proceedings of CPAIOR’18*, pages 245–253, 2018.

10 P.A. Geelen. Dual viewpoint heuristics for binary constraint satisfaction problems. In *Proceedings of ECAI’92*, pages 31–35, 1992.

11 C. Gomes, B. Selman, N. Crato, and H. Kautz. Heavy-tailed phenomena in satisfiability and constraint satisfaction problems. *Journal of Automated Reasoning*, 24:67–100, 2000.

12 D. Habet and C. Terrioux. Conflict history based heuristic for constraint satisfaction problem solving. *Journal of Heuristics*, 27(6):951–990, 2021.

13 R.M. Haralick and G.L. Elliott. Increasing tree search efficiency for constraint satisfaction problems. *Artificial Intelligence*, 14:263–313, 1980.

14 C. Lecoutre. *Constraint networks: techniques and algorithms*. ISTE/Wiley, 2009.

15 C. Lecoutre. ACE, a generic constraint solver. CoRR, abs/2302.05405, 2023. doi:10.48550/arXiv.2302.05405.

16 C. Lecoutre and O. Roussel. Proceedings of the 2018 XCSP3 competition. CoRR, abs/1901.01830, 2019. arXiv:1901.01830.

17 C. Lecoutre, L. Sais, S. Tabary, and V. Vidal. Recording and minimizing nogoods from restarts. *Journal on Satisfiability, Boolean Modeling and Computation (JSAT)*, 1:147–167, 2007.

18 C. Lecoutre, L. Sais, S. Tabary, and V. Vidal. Reasonning from last conflict(s) in constraint programming. *Artificial Intelligence*, 173(18):1592–1614, 2009.

19 C. Lecoutre and N. Szczepanski. PyCSP3: modeling combinatorial constrained problems in Python. CoRR, abs/2009.00326, 2020. arXiv:2009.00326.

20 H. Li, M. Yin, and Z. Li. Failure based variable ordering heuristics for solving CSPs. In *Proceedings of CP’21*, 2021.

21 J.J. McGregor. Relational consistency algorithms and their application in finding subgraph and graph isomorphisms. *Information Sciences*, 19:229–250, 1979.

22 L. Michel and P. Van Hentenryck. Activity-based search for black-box constraint programming solvers. In *Proceedings of CPAIOR’12*, pages 228–243, 2012.

23 G. Pesant. From support propagation to belief propagation in constraint programming. *Journal of Artificial Intelligence Research*, 66:123–150, 2019.

24 G. Pesant, C.-G. Quimper, and A. Zanarini. Counting-based search: Branching heuristics for constraint satisfaction problems. *Journal of Artificial Intelligence Research*, 43:173–210, 2012.

25 P. Refalo. Impact-based search strategies for constraint programming. In *Proceedings of CP’04*, pages 557–571, 2004.

26 F. Rossi, P. van Beek, and T. Walsh, editors. *Handbook of Constraint Programming*. Elsevier, 2006.
27 J. Vion and S. Piechowiak. Une simple heuristique pour rapprocher DFS et LNS pour les COP. In *Proceedings of JFPC'17*, pages 39–45, 2017.

28 H. Wattez, C. Lecoutre, A. Paparrizou, and S. Tabary. Refining constraint weighting. In *Proceedings of ICTAI'19*, pages 71–77, 2019.

29 N.-F. Zhou. An XCSP3 Solver in Picat. In *XCSP3 Competition 2022 Proceedings*, 2022.
<figure><img src="image_79.png" /></figure>
# Incremental Constrained Clustering by Minimal Weighted Modification

Aymeric Beauchamp ✉️ ID

University of Orléans, INSA Centre Val de Loire, LIFO EA 4022, France

Thi-Bich-Hanh Dao ✉️ ID

University of Orléans, INSA Centre Val de Loire, LIFO EA 4022, France

Samir Loudni ✉️ ID

TASC (LS2N-CNRS), IMT Atlantique, France

GREYC, University of Caen Normandy, France

Christel Vrain ✉️ ID

University of Orléans, INSA Centre Val de Loire, LIFO EA 4022, France

## Abstract

Clustering is a well-known task in Data Mining that aims at grouping data instances according to their similarity. It is an exploratory and unsupervised task whose results depend on many parameters, often requiring the expert to iterate several times before satisfaction. Constrained clustering has been introduced for better modeling the expectations of the expert. Nevertheless constrained clustering is not yet sufficient since it usually requires the constraints to be given before the clustering process. In this paper we address a more general problem that aims at modeling the exploratory clustering process, through a sequence of clustering modifications where expert constraints are added on the fly. We present an incremental constrained clustering framework integrating active query strategies and a Constraint Programming model to fit the expert expectations while preserving the stability of the partition, so that the expert can understand the process and apprehend its impact. Our model supports instance and group-level constraints, which can be relaxed. Experiments on reference datasets and a case study related to the analysis of satellite image time series show the relevance of our framework.

**2012 ACM Subject Classification** Theory of computation → Constraint and logic programming; Computing methodologies → Semi-supervised learning settings

**Keywords and phrases** Incremental constrained clustering, Constrained optimization problem, User feedback

Digital Object Identifier 10.4230/LIPIcs.CP.2023.10

Related Version Extended Version: hal-04158825

**Supplementary Material** Software (Source Code): https://github.com/aymericb213/IAC
archived at swh:1:dir:8f611496af5ad016912ca3d55379b50b36541f89

**Funding** This work was supported by the French national research project HERELLES under grant agreement ANR-20-CE23-0022.

**Acknowledgements** The authors want to thank the anonymous reviewers for their comments and suggestions which helped to improve this paper.

## 1 Introduction

Clustering is a popular task in Data Mining in which data instances (i.e. points of a dataset) are grouped into distinct clusters according to their similarity. Over time, many strategies have been explored to compute data partitions, each of them having their own strengths and biases. Constrained clustering [16] aims to find relevant clusters by stating some desired properties in the form of constraints, thus alleviating the aforementioned biases. The most

<figure><img src="image_80.png" /></figure>

© Aymeric Beauchamp, Thi-Bich-Hanh Dao, Samir Loudni, and Christel Vrain;

licensed under Creative Commons License CC-BY 4.0

29th International Conference on Principles and Practice of Constraint Programming (CP 2023).

Editor: Roland H. C. Yap; Article No. 10; pp. 10:1–10:22

<figure><img src="image_81.png" /></figure>

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
## 10:2 Incremental Constrained Clustering by Minimal Weighted Modification

commonly used constraints for clustering state that two points must be clustered together (*must-link*) or apart (*cannot-link*) [32]. Other constraints exist as for instance limiting the size or the diameter of the clusters [1]. They are usually given by a subject matter expert (SME) that possesses domain knowledge on the data; they model his/her knowledge and expectations about the result by expert constraints.

In practice, it may be difficult for an expert to express constraints solely on the basis of the data. It is usually easier to him/her to provide feedback on the current partition to refine it. This leads to a *human-in-the-loop* clustering process, where new constraints given or validated by the expert are incrementally integrated, modifying the result until user satisfaction. However, this raises several non trivial questions such as how to elicit the constraints, how to integrate them, or how to further exploit them. By incorporating human feedback and domain knowledge, the resulting clusters are more likely to align with the expert expectations while making them more intuitive and interpretable. This incremental setting mimics the natural step-by-step progression of an exploratory task such as clustering, where the user needs to iterate several times before satisfaction. In such a process, the result at each step should (1) be computed fast enough, (2) exploit the expert constraints efficiently (3) be similar to the result of the previous step, in order not to disturb the expert.

A naive way to integrate new expert constraints is to restart a constrained clustering algorithm. However, this presents at least two weaknesses: it starts from scratch without ever considering intermediate results and the new constraints can lead to a partition very different from the one previously shown to the expert. To cope with this problem, we propose a first generic framework that allows for truly interactive and iterative constrained clustering that fulfills the conditions mentioned above. Our main contributions are:

- an incremental constrained clustering framework designed for human interaction, combining active constraint selection and clustering modification;

- a new constraint programming model for minimal clustering modification, which ensures the stability of the partition.

The paper is organized as follows. We review in Section 2 related work on incremental constrained clustering and minimal clustering modification and present our method in Section 3. Experiments on reference and satellite image time series datasets are presented in Section 4 and perspectives are discussed in Section 5.

## 2 Related Work

The first works on using constraints in clustering are extensions of classic algorithms for handling *must-link/cannot-link* constraints [33, 10, 35]. They either search for a solution satisfying all the constraints [33] or a compromise between constraint satisfaction and clustering quality [6]. Thereafter methods allowing to integrate more general constraints, using declarative frameworks such as SAT [12], ILP [29] or CP [9] have been proposed. When new constraints are given by the user, all these methods require to restart from scratch, without any guarantee that the new partition will be similar to the previous one. Therefore they are not suited for an incremental setting.

There is a growing body of works related to incremental constrained clustering. In [7], the idea of gathering feedback from an existing result rather than expecting the user to provide insightful constraints by themselves is demonstrated. Later, the authors of [11] studied the problem of adding or removing a constraint from a constraint set satisfied by a partition. They described conditions under which the problem is easy to solve and an algorithm working under these conditions. In [26], a cluster refinement framework uses subclustering to find
representatives to present to the user before learning an embedding using the feedback. More recently, [21] proposed an incremental variant of a collaborative constrained clustering system that integrates constraint satisfaction in its objective function. Among declarative approaches, an ILP model for minimal clustering modification (MCM) is proposed in [20]. It constrains cluster diameters with the objective of removing undesirable properties from a partition. Another ILP model [28] computes a membership score for each point to each cluster and optimizes a criterion based on this score. It can modify the assignment of points to clusters, even if they are not involved in constraints. However those models are restrictive since they cannot handle conflicting constraints, which make the problem unsolvable.

Our approach seeks to preserve the general cluster structure of an existing partition while satisfying new constraints. To the best of our knowledge, this is the first approach explicitly tackling both quality and stability. It is based on CP, therefore it can integrate different kinds of constraints that can be either soft or hard, with some control over constraint relaxation and the ability to handle conflicting constraints. We use subclustering in a similar way to [26], albeit for a different purpose as it allows generalizing the changes decided by our CP model as well as finding cluster representatives used in our objective function.

3 Incremental and Active Clustering Framework

<figure><img src="image_82.png" /><figcaption>Figure 1 Schematic view of the incremental clustering cycle.</figcaption></figure>

In this section, we describe our proposed incremental and active clustering (IAC) framework. Figure 1 gives a general overview of IAC. The incremental constrained clustering loop starts from an initial partition computed by any clustering algorithm. This partition is then shown to the user to collect his/her general feedback (satisfied or not). If he/she is not satisfied, he/she can modify this partition in two ways: by manually providing a set of constraints and/or by inferring it through an active constraint selection method. He/she can also set the proportion of constraints to satisfy as well as the scope of modifications. The clustering modification step updates the current partition according to these new constraints, optionally generalizing modifications to unconstrained data instances. The output is a new partition satisfying the constraints while preserving its stability, i.e. similar to the previous one. This process is repeated until the user is satisfied by the resulting partition. It is noteworthy that our framework is generic as any active constraint selection method and
## 10:4 Incremental Constrained Clustering by Minimal Weighted Modification

any modification algorithm could be used as long as the constraints generated match the constraints handled in modification. We formulate the problem of clustering modification in a declarative way and present a CP model. This has the benefit of being able to integrate several types of constraints for the user feedback.

### 3.1 Minimal Weighted Clustering Modification

We consider the minimal weighted clustering modification (MWCM) problem: given a partition $\mathcal{P}$ of $N$ data instances (numbered from 1 to $N$) into a number $K$ of clusters, and a set of user constraints $\mathcal{C}$, the objective is to find a new partition $\mathcal{P}'$ such that the constraints are satisfied while minimizing some function $f$ modeling the difference between $\mathcal{P}$ and $\mathcal{P}'$. For solving this problem, Algorithm 1 shows the different steps that are detailed below.

**Algorithm 1** Minimal Weighted Clustering Modification.

**Input:** Dataset $\mathcal{X}$, partition $\mathcal{P}$, constraints $\mathcal{C}$, anchor generation rate $\alpha$, super-instance rate $\beta$, constraint satisfaction rate $\delta$

**Output:** modified partition $\mathcal{P}'$

1: anchors $\leftarrow$ COMPUTEREPRESENTATIVES($\mathcal{X}, \mathcal{P}, \alpha$) ▷ See section 3.1.1

2: $X \leftarrow$ COMPUTECOPINSTANCES($\mathcal{X}, \mathcal{P}, \mathcal{C}, \beta$) ▷ Instances used in COP (Section 3.1.2)

3: $\mathcal{D} \leftarrow$ DISTANCEMATRIX($X, anchors$) ▷ See section 3.1.1

4: $p \leftarrow$ GETCONSTRAINEDPARTITION($X, \mathcal{P}$) ▷ Cluster membership of the constrained instances

5: mods $\leftarrow$ SOLVEMODEL($\mathcal{D}, p, \mathcal{C}, \delta$) ▷ Solves the COP in Section 3.2

6: **return** APPLYMODIFICATIONS(mods, $\mathcal{P}$) ▷ Updates $\mathcal{P}$ and generalizes modifications

#### 3.1.1 Objective function and anchors

A straightforward candidate for $f$ is to count the number of instances that have changed their cluster membership between $\mathcal{P}$ and $\mathcal{P}'$ [20]:

$$ \arg \min \sum_{i=1}^{N} \mathbb{I}(\mathcal{P}[i] \neq \mathcal{P}'[i]) \qquad (1) $$

where $\mathcal{P}[i]$ denotes the number $c \in [1, K]^1$ of the cluster containing instance $i \in [1, N]$ and $\mathbb{I}$ the indicator function that returns 1 if the expression given as argument is true, and 0 otherwise. The main drawback of Equation (1) is that it does not take into account the structure of the clusters. For example, it is reasonable to consider that putting two instances in a nearby cluster is more akin to the idea of minimal modification than putting one instance into a faraway cluster. Therefore we propose an alternate objective function that integrates a distance-based weighting of the modifications:

$$ \arg \min \sum_{i=1}^{N} \mathbb{I}(\mathcal{P}[i] \neq \mathcal{P}'[i]) \mathcal{D}[i, \mathcal{P}'[i]] \qquad (2) $$

where $\mathcal{D}$ is a distance matrix of dimensions $N \times K$ such that $\mathcal{D}[i, c]$ is the distance of instance $i$ to cluster $c$. This objective function measures the changes between $\mathcal{P}$ and $\mathcal{P}'$ and keeps the cumulative distances resulting from these changes small. This objective function is integrated into the model for MWCM that will be presented in Section 3.2.

<footnote n="1">¹ We use the notation $[1, K]$ for the set $\{1, ..., K\}$.</footnote>
A simple way to compute $\mathcal{D}[i, c]$ is to set $\mathcal{D}[i, c] = d(i, \mu_c)$, where $\mu_c$ is the representative (medoid or centroid) of $c$ and $d$ is a distance measure, typically the Euclidean distance. This has however a known limitation: the modifications made by the model implicitly treat all clusters as spherical. It can be counterproductive if the user seeks more complex shaped clusters. To overcome this, we use anchors [17] such that each cluster will be represented by a subset of its instances, in order to better represent its structure. Anchors are computed by dividing each cluster $c$ of the input partition into smaller sub-clusters using Single-Link hierarchical clustering. For each of these sub-clusters, an anchor is defined as the instance minimizing the sum of its distances with the other instances of the sub-cluster. Using anchors, $\mathcal{D}[i, c]$ represents the distance of instance $i$ to its closest anchor belonging to cluster $c$. A parameter $\alpha$ defining the proportion of instances per cluster (in percentage) that will become anchors is used. For example, a rate of 0% means that we only use medoids, while at a rate of 100%, all instances are anchors (cf. Figure 2).

<figure><img src="image_83.png" /><figcaption>Figure 2 Anchor positions for different values of $\alpha$ computed from the partition generated by Kmeans on lsun dataset. The anchors are represented bigger than normal instances.</figcaption></figure>

### 3.1.2 Generalizing constraints with super-instances

In a real use case, we assume that the expert will only react on a small number of instances per iteration. As a result, the clustering modification could become unnoticeable when the dataset size is large compared to the number of constraints, which is a fairly common case. Hence, exploiting expert feedback to *generalize* the modifications to relevant unconstrained instances is an important issue w.r.t. the challenge of asking a reasonable number of queries. This issue has been highlighted in [34], and is especially relevant in our incremental setting. Furthermore, this generalization must be controlled to ensure the expert can grasp the scope of potential modifications. We assume that, in most cases, the expert will want to modify a zone around the selected instances and not only the instances themselves. Bearing this intuition in mind, making use of nearest neighbors or a proximity radius seem adequate, but these methods lack predictability: if an instance is within the radius of two constrained points who were reassigned to different clusters, determining how to resolve the generalization is not obvious.

We propose to use super-instances, i.e. virtual instances grouping several real data points, to generalize the modifications. Note that *generalization* does not mean that new constraints are generated, rather that the super-instances are passed directly to the model instead of the data instances (see Appendix C). Thus any modification in the cluster membership of a super-instance amounts to changing the membership of every real data instance that compose it. The generalization scope is controlled as follows: *the less a cluster is divided, the stronger the impact of a modification*. Thus the scope depends on the number of super-instances, determined by a rate $\beta$ proportional to the cluster size. As such, setting $\beta$ to e.g. 10%
## 10:6 Incremental Constrained Clustering by Minimal Weighted Modification

means that each cluster will be split into a number of super-instances equal to 10% of its size, with 1 being equivalent to not generalizing at all. Super-instances are determined by sub-clustering each cluster of the current partition into small groups, each group representing a super-instance. We empirically found that complete-link hierarchical clustering is adequate, despite its memory usage which makes it unsuitable on large datasets. Other alternatives include the density-based OPTICS as in [26], or the Furthest Point First (FPF) algorithm [15]. However, user constraints defined on instances need to be transferred to super-instances. This may raise potential conflicts. To avoid this pitfall, we ensure that every super-instance contains no more than one constrained data point. If this is not the case, we split the super-instance using the constrained instances as centers of the new split super-instances. An illustrative case is given in Appendix A.

## 3.2 Constraint Optimization Problem Formulation

Taking advantage of declarative approaches, we formulate the problem of finding a similar partition satisfying user constraints as a Constraint Optimization Problem (COP). In the following, we use the term instance to denote a data instance or a super-instance if it is used.

**Variables and Objective Function.** Only instances that are subject to the constraints will be concerned by the COP. Function `GETCONSTRAINEDPARTITION` in Algorithm 1 produces the subset $X$ containing the constrained instances. For each instance in $i \in X$, we define a variable $G_i$ with the domain $[1, K]$, where $G_i = c$ means instance $i$ is assigned to cluster $c$ in the new partition $\mathcal{P}'$. Using Eq. (2), the objective function is:

$$ \arg\min \sum_{i \in X} \mathbb{I}(G_i \neq \mathcal{P}[i]) \mathcal{D}[i, G_i] \qquad (3) $$

**User constraints.** Several instance-level and group-level constraints can be expressed in our model, as below. Must-link (ML)/cannot-link (CL) constraints on two instances $i, j$ stating that the instances must/cannot be in the same cluster, can be expressed by $G_i = G_j$ for ML and $G_i \neq G_j$ for CL. We also compute the transitive closure on ML/CL constraints [25], which derives supplementary constraints according to three rules: (i) if $ML(a, b)$ and $ML(b, c)$, then $ML(a, c)$; (ii) similarly $ML(a, b)$ and $CL(b, c)$ imply $CL(a, c)$; (iii) in a binary clustering case ($K = 2$), $CL(a, b)$ and $CL(b, c)$ imply $ML(a, c)$.

Triplet constraint $(a, p, n)$ [23] states that a reference instance $a$ is more similar to instance $p$ than to instance $n$. Instance $p$ is therefore called positive instance and $n$ negative. This constraint can be expressed using an implication constraint as follows:

$$ G_a = G_n \implies G_a = G_p \qquad (4) $$

Span-limited constraints [27] restricting the span of a set of instances $S \subseteq X$. A specific span-limited constraint states that the instances of $S$ must be assigned only to clusters from a given subset $C \subseteq [1, K]$. It can be expressed using the count global constraint<fnref n="2" />:

$$ \text{count}(c, [G_i | i \in S], =, 0) \quad \forall c \notin C \qquad (5) $$

A generic span-limited constraint specifies that the instances of $S$ must be assigned to at most a number $\gamma$ of clusters. It can be expressed using the constraint `atmost_nvalue`<fnref n="3" /> [5]:

$$ \text{atmost\_nvalue}(\gamma, [G_i | i \in S]) \qquad (6) $$

<footnote n="2">² https://sofdem.github.io/gccat/gccat/Ccount.html</footnote>

<footnote n="3">³ https://sofdem.github.io/gccat/gccat/Catmost_nvalue.html</footnote>
More thematic expert feedback can also be integrated as implication constraints. In our model, they are of the form $P \Rightarrow Q$, where P and Q are conjunctions of simpler constraints such as ML/CL.

**Cluster creation.** Our model allows $P'$ to have more clusters than $P$, by defining the domain of $G_i$ as $[1, K']$, with $K' > K$. This enables assigning an instance to a new cluster if some constraints are unsatisfied. For example, if $K = 2$ and the expert states three constraints $G_i \neq G_j, G_j \neq G_l$ and $G_l \neq G_i$, then it is necessary to create a new cluster. Given its conditions of apparition, it will typically be very small. In order to prevent the model from creating clusters because it would be optimal to do so w.r.t. the objective function, we set the distance $D[i, k']$ to a value greater than all other distances for $K < k' \le K'$.

**Relaxing constraints.** During the incremental process, the expert could make mistakes when trying to improve the partition, which would result in adding conflicting constraints, thus leading to an over-constrained CP model. We reify the user constraints to gain control over constraint satisfaction. Each constraint $c \in C$ is associated with a Boolean variable $S_c$ such that $S_c = 1$ iff c is satisfied. The satisfaction rate $\delta$ sets a lower/upper bound or the exact value of the number of constraints the model must satisfy:

$$ \sum_{c \in C} S_c \overset{\geq}{\leq} \delta \cdot |C| \qquad (7) $$

Constraint relaxation can both solve problems with conflicting constraints and ignore - or warn the user about - constraints that would modify instances far from their new cluster. Relaxation however increases runtime due to the additional Boolean variables. As is, our framework automatically detects ML/CL conflicts and reduces the satisfaction rate accordingly.

**Managing the constraint store.** Algorithm 1 solves this COP with the user constraints collected at each iteration. The incremental setting raises the issue of managing the constraints between iterations. It is possible to store every constraint received since the beginning of the process to ensure that all expert feedbacks are respected. In the experiments, we choose instead to treat the constraint set given at each iteration independently. In this way, if the expert adds a constraint in conflict with another one given previously, we consider that the user is simply rescinding some of his/her feedback. The expert can also mark some constraints as mandatory so that they are kept satisfied throughout the process. Appropriately managing the constraint store is a potential future research lead.

## 3.3 Active Constraint Selection

Obtaining constraints manually can be costly. This motivates active constraint selection methods [2, 25, 37], which select the most informative constraints to query. To evaluate the interest of exploiting an active constraint selection approach within our framework, we use NPU [37], a neighborhood-based sampling strategy. Neighborhoods $N$ are groups of instances whose cluster assignment is certain, they represent the underlying clusters. NPU iteratively builds the neighborhoods by selecting the most informative instance $x^*$ and querying its relation with respect to existing neighborhoods. The informativeness of an instance $x$ is defined by the ratio $H(N|x)/E[q(x)]$, where $H(N|x)$ is the entropy measure of the uncertainty to assign $x$ to a neighborhood in $N$, and $E[q(x)]$ the expected number of queries needed to discover its neighborhood.
10:8 Incremental Constrained Clustering by Minimal Weighted Modification

Exploiting this informativeness, we adapt the NPU framework to the incremental setting. The neighborhoods N, which are initially empty, are constructed and kept throughout the iterations. For each informative instance x*, queries are put on the membership relation between x* and each neighborhood N ∈ N. Once the user answers favorably, a ML constraint is created with the queried neighborhood, otherwise a CL is created. If no ML is achieved, a new neighborhood is created for x* (see the detailed algorithm in Appendix B). In relation to constraint management between iterations, it must be pointed out that the preservation of the neighborhoods between iterations prevents selecting constraints conflicting with those selected in earlier iterations. Indeed, an instance stored in a neighborhood is never picked again by NPU, and is only used to generate constraints with an instance that has not been presented to the user before. The only factor that could cause previous constraints to be involuntarily relaxed is a high generalization scope. We found no such occurrences in our experiments with the values we tested for β (see Section 4.2.1). Fig. 3 illustrates the framework walkthrough on a toy example, with noticeable separation between non-spherical clusters that KMEANS is unable to recover.

<figure><img src="image_84.png" /><figcaption>Figure 3 Illustration of IAC with (α = 20%, β = 30%) over 10 iterations on lsun dataset, starting from a KMEANS partition (Fig. 3a). Subsequent figures show the evolution of the partition after t iterations of IAC with NPU. Adjusted Rand Index with ground truth is reported on each figure.</figcaption></figure>

4 Experiments

In this section, the experiments aim to answer the following research questions (RQ):

1. What effect do IAC parameters (α and β) have on clustering results ? (Section 4.2.1)

2. How does our CP model scale with the number and type of constraints ? (Section 4.2.2)

3. How does the constraint relaxation of IAC compare with other methods that use soft constraints ? (Section 4.2.3)

4. How effective is IAC in an active constraint selection context ? (Section 4.2.4)

5. What is the performance of our framework in terms of clustering quality, partition similarity and runtime when compared to state of the art methods ? (Section 4.2.4)

6. How effective is our framework on a real use case with human feedback ? (Section 4.3)

4.1 Experimental Methodology

**Evaluation measures.** For all experiments we use datasets for which a *ground-truth* labeling is known. The produced partition is then compared to the known partition using an external measure. A high value of the measure indicates a good partition, meaning that the clustering algorithm has successfully identified the already known structure. We consider three measures: Adjusted Rand Index (ARI) [19], Adjusted Mutual Information (AMI) [31] and Folkes-Mallows Index (FMI) [13]. ARI measure is defined by:

$$
ARI(\mathcal{P}, \mathcal{P}') = \frac{2(ab - cd)}{(a+d)(d+b) + (a+c)(c+b)} \quad (8)
$$
where a (resp. b) is the number of instances pairs clustered together (resp. apart) in $\mathcal{P}$ and $\mathcal{P}'$, and c (resp. d) is the number of pairs clustered together (resp. apart) in $\mathcal{P}$ and apart (resp. together) in $\mathcal{P}'$. AMI is a variation of Normalized Mutual Information corrected for chance:

$$AMI(\mathcal{P}, \mathcal{P}') = \frac{MI(\mathcal{P}, \mathcal{P}') - E(MI(\mathcal{P}, \mathcal{P}'))}{\max(H(\mathcal{P}), H(\mathcal{P}')) - E(MI(\mathcal{P}, \mathcal{P}'))} \quad (9)$$

where MI is the measure of mutual information, H the entropy of a partition, and E the expected mutual information if the partitions are random. Finally FMI is the geometric mean of precision and recall, using one partition as a reference for the other:

$$FMI(\mathcal{P}, \mathcal{P}') = \sqrt{\frac{TP}{TP + FP} \cdot \frac{TP}{TP + FN}} \quad (10)$$

We also measure the runtime of the modification step (reclustering or MWCM) since it is crucial to take it into account in the incremental setting. We set a timeout after 1 hour of modification time.

**Experimental protocol.** We have implemented our CP model in Python 3.11 using the CPMpy library [18], interfacing with CP-SAT solver from or-tools<fnref n="4" />. We chose this library because it allows to easily use landmark ML libraries such as scikit-learn together with CP. Code for reproducing the experiments is available at the repository given at the summary of this paper. The implementation of NPU and all clustering algorithms we considered for comparison (COPKMeans, PCKMeans and MPCKMeans) are from the active-semi-supervised-clustering<fnref n="5" /> library. All experiments were run on a computer with two 48-core Intel Xeon processors at 4 GHz and 64 GB of RAM running Ubuntu 20.04.

For each dataset, we first generate an initial partition with KMEANS [24] with K set to the true number of clusters. Queries correspond to pairwise constraints. We emulate user feedback using the ground truth labeling of the data as an oracle i.e. a must-link constraint is added if the selected pair of instances belong to the same class, and a cannot-link otherwise. For RQs 1, 4 and 5, we perform 10 iterations of selection-modification loop. At each iteration, we use the current partition to select a batch of 10 queries with NPU, get feedback from the oracle, and apply either a constrained clustering algorithm to the full dataset or our CP model for cluster modification. In order to smooth out the random effects occurring in the partition initialisation and in constraint selection, we repeat each experiment 90 times.

To evaluate the overall performance over the 11 successive partitions (including the initial partition) obtained for each run, we compute for each metric the area under the budget curve (AUBC) [38]) for different fixed budgets of queries to ask the user. Given the budget curve, the AUBC is calculated by the trapezoid method, and the higher value reflects better performance of the evaluated method under varying budgets. For each metric, we compute two types of AUBC: $AUBC_{quality}$ when comparing the successive partitions to the ground truth partition, and $AUBC_{similarity}$ when comparing two consecutive intermediate partitions. Since $AUBC_{similarity}$ values are defined over the interval [0, 0.9], we perform a min-max normalization so that all metrics are defined over the [0, 1] range. To statistically compare the performance of different algorithms and/or configurations of the same algorithm for different parameter settings on multiple data sets, we resort to Bayesian pairwise comparison [4] using

<footnote n="4">⁴ https://developers.google.com/optimization</footnote>

<footnote n="5">⁵ https://github.com/datamole-ai/active-semi-supervised-clustering</footnote>
10:10 Incremental Constrained Clustering by Minimal Weighted Modification

<table><caption>**Table 1** Dataset Characteristics, with *N* the number of instances, *A* the number of features and *K* the number of clusters or classes.</caption><thead><tr><th colspan="4">UCI</th><th colspan="4">FCPS</th></tr><tr><th>Name</th><th>(N, A, K)</th><th>Name</th><th>(N, A, K)</th><th>Name</th><th>(N, A, K)</th><th>Name</th><th>(N, A, K)</th></tr></thead><tbody><tr><td>iris</td><td>(150, 4, 3)</td><td>ionosphere</td><td>(351, 34, 2)</td><td>lsun</td><td>(400, 2, 3)</td><td>chainlink</td><td>(1000, 3, 2)</td></tr><tr><td>wine</td><td>(178, 13, 3)</td><td>yeast</td><td>(1484, 8, 10)</td><td>target</td><td>(770, 2, 6)</td><td>wingnut</td><td>(1016, 2, 2)</td></tr><tr><td>sonar</td><td>(208, 60, 2)</td><td>statlog</td><td>(2310, 19, 7)</td><td>atom</td><td>(800, 3, 2)</td><td>engytime</td><td>(4096, 2, 2)</td></tr><tr><td>glass</td><td>(214, 9, 6)</td><td>Letters</td><td>(20000, 16, 26)</td><td></td><td></td><td></td><td></td></tr><tr><td>ecoli</td><td>(336, 7, 8)</td><td>MNIST</td><td>(70000, 784, 10)</td><td></td><td></td><td></td><td></td></tr></tbody></table>

baycomp<fnref n="6" /> library. The principle is to use Bayes' rule to update a prior statistical distribution representing the null hypothesis (both compared algorithms have the same performance), with a likelihood function modeling the experimental observations. We then get a posterior distribution, reflecting how the prior belief has changed, taking the observations into account. Using the Markov chain Monte Carlo method, the posterior is sampled 50,000 times to estimate the probability of one algorithm being better than the other as well as the probability of being in the region of practical equivalence (or rope). In practice, querying the posterior distribution allows to simulate repeating the whole experimental process and to quantify the likelihood of our results. We choose to fix to 1% the difference of performance between the methods as the rope.

## 4.2 UCI and FCPS Datasets

In this section, we report experimental results on ten real-world datasets from the UCI repository<fnref n="7" /> and on six synthetic datasets from the FCPS [30] suite designed to address specific challenges to the clustering algorithms such as lack of linear separability, classes defined by data density rather than data spacing, no cluster structure at all, etc. A summary of the basic characteristics is given in Table 1. We used the versions of datasets available under the library clustering-benchmarks<fnref n="8" /> [14].

### 4.2.1 Parameter Settings of IAC

To answer RQ1, we evaluate the effects of different parameter settings of the clustering modification step of our IAC framework: the anchor generation rate $\alpha \in \{0\%, 5\%, 20\%\}$ and the super-instance generation rate $\beta \in \{10\%, 30\%, 50\%, 100\%\}$. This makes a total of 12 configurations of parameter combinations to evaluate. Recall that $\alpha = 0\%$ means that we only compute the cluster medoids, while $\beta = 100\%$ means no generalization by super-instances is performed. For each configuration and each metric, we perform Bayesian pairwise comparison according to AUBC<sub>quality</sub> and AUBC<sub>similarity</sub> values for each of the three metrics ARI, AMI and FMI over all the datasets and count the number of wins. More precisely, given a pairwise comparison between two configurations $conf_1$ and $conf_2$, using a Bayesian hierarchical model [8], we get three probabilities: the probability that $conf_1$ has higher scores than $conf_2$, the probability that differences are within the region of practical equivalence (rope), or that $conf_2$ has higher scores. If ($p_{conf_1} > p_{conf_2} + p_{rope}$), then we count this comparison as winning for $conf_1$.

<footnote n="6">⁶ https://baycomp.readthedocs.io/en/latest/index.html</footnote>

<footnote n="7">⁷ https://archive.ics.uci.edu/ml/index.php</footnote>

<footnote n="8">⁸ https://clustering-benchmarks.gagolewski.com/weave/suite-v1.html</footnote>
<figure><img src="image_85.png" /><figcaption>Figure 4 Pairwise plot comparing the number of wins for each configuration (α, β) using a Bayesian hierarchical model, according to AUBC<sub>quality</sub> (horizontal) or AUBC<sub>similarity</sub> (vertical). Best values are close to the top left, with configurations on the Pareto front in orange.</figcaption></figure>

The table of wins is available in Appendix D. Configurations using anchors (i.e. α ≠ 0) ensure the highest AUBC<sub>quality</sub> values in almost all configurations. Additionally, the best values are obtained with α = 20%. For a fixed value of α, AUBC<sub>quality</sub> values increase significantly with the decrease of β, 10% being the best value. This result suggests that a large scope of generalization does not strongly impact the modification of the clustering and thus its quality. AUBC<sub>similarity</sub> values have the opposite behavior: β = 10% is the worst setting for similarity. Note that the impact of α seems negligible compared to β. It slightly improves similarity when β is low, while its contribution seems negligible in the absence of generalization. Generalizing modifications understandably degrades similarity, albeit not dramatically. A Pareto front made of three configurations on every metric (see Fig. 4) emerges from the results: one best in quality (20%, 10%), another best in similarity (20%, 100%), and a compromise setup (20%, 30%).

Figure 5 plots the posterior distribution of pairwise comparison of configurations in a simplex. The distribution is shown as a triangle with regions corresponding to different samples of the distribution, e.g. Figure 5a compares the posterior distribution obtained with α = 0% (medoids only) and β = 10% with the one obtained without generalization. This indicates that there is a 88.9% probability that using medoids combined with SI for this value of β is better than using medoids without SI.

<figure><img src="image_86.png" /><figcaption>Figure 5 Simplex view of Bayesian comparison of two configurations w.r.t AUBC<sub>quality</sub> of ARI. Each sample is plotted according to probabilities p<sub>conf1</sub> (left), p<sub>rope</sub> (top) and p<sub>conf2</sub> (right).</figcaption></figure>

### 4.2.2 Impact of Number and Types of User Constraints on Runtime

In this section, we answer **RQ2** by studying two points: the computational efficiency of our CP model and the expressiveness of our approach (see Sect. 3.2). For these experiments, we consider 3 sizes of constraint set (10, 100, 1000). For each test case, we randomly generate sets of four types of constraints (pairwise, triplet, span-limited specific and generic). We compute an initial KMeans partition and run our CP model only once for each set of constraints and
10:12 Incremental Constrained Clustering by Minimal Weighted Modification

<figure><img src="image_87.png" /><figcaption>Figure 6 Evolution of running time of our CP model for the two largest datasets when varying the number and type of constraints, with 95% confidence interval. CPU times are in log-scale.</figcaption></figure>

we report the average CPU times over 90 runs. For pairwise constraints, we disable the computation of transitive closure to keep the number of constraints unchanged. Span-limited constraints are created by randomly choosing 10 instances and finding the ground truth set of clusters - or number of clusters, in the generic case - to which the group belongs.

**Runtime analysis.** Figure 6 shows the results we obtained for *letters* and *MNIST* datasets with $\alpha = 0\%$ and $\beta = 100\%$. Our CP model can process 10 ML/CL constraints in less than 0.05 seconds, while for triplet and specific span-limited constraints the runtime reaches 0.035 and 0.15 seconds respectively. This seems very reasonable in an incremental context. For the *yeast* dataset, with 100 pairwise constraints, it takes 0.35 seconds, whereas for triplet and span-limited constraints the runtime increases up to 1.36 seconds. With 1000 constraints, the runtime is more than 44s for specific span-limited, 67s for ML/CL, and over 180s for triplet; generic span-limited constraints took up too much memory to finish. However, in practice, the number of constraints expected from the user is in the tens rather than the hundreds or thousands. Triplet and span-limited constraints show a substantial increase in runtime compared to ML/CL constraints.

**Mixed constraint types.** One of the advantages of our approach is its ability to easily combine different types of constraints without the need to create a specialized algorithm. To demonstrate this ability, we compared two composite settings:

- **mlcl+triplet**: generate pairwise and triplet constraints in equal proportions, e.g. 50 ML/CL and 50 triplet constraints for the 100 constraints case.

- **all**: similar to **mlcl+triplet**, except a pairwise (resp. triplet) constraint is replaced by a specific (resp. generic) span-limited constraint.

As far as we are aware, such problems cannot be solved by any existing techniques. Problems with ML/CL and triplet constraints take much less time to solve than those involving the three types of constraints, the latter being much more time-consuming. Surprisingly, with 1000 constraints, the ML/CL+Triplet combination takes less time compared to the case where only triplet constraints are involved. All these results underline the relevance of carefully selecting a small number of constraints to guarantee a good compromise between efficiency and quality of the final clustering.
<table><caption>**Table 2** Comparative study for clustering with pairwise constraints relaxation for the mk2 dataset. Metrics are ARI with ground truth (Quality), ARI with unconstrained KMeans (Similarity), runtime and number of constraints relaxed in the solution ($n_r$).</caption><thead><tr><th rowspan="2"></th><th colspan="4">Test case with conflicts</th><th colspan="4">Test case with δ = 94%</th></tr><tr><th>Quality</th><th>Similarity</th><th>Time</th><th>n<sub>r</sub></th><th>Quality</th><th>Similarity</th><th>Time</th><th>n<sub>r</sub></th></tr></thead><tbody><tr><td>IAC+Anchors</td><td>0.576</td><td><strong>0.075</strong></td><td>5.262</td><td>49.7</td><td>0.309</td><td><strong>0.177</strong></td><td>3.051</td><td>60.1</td></tr><tr><td>IAC+Anchors+SI</td><td><strong>0.760</strong></td><td>0.024</td><td>4.868</td><td>49.83</td><td>0.393</td><td>0.108</td><td>2.972</td><td>60.1</td></tr><tr><td>PCK-Means</td><td>0.081</td><td>0.051</td><td><strong>3.639</strong></td><td>149.3</td><td>0.375</td><td>0.045</td><td><strong>2.834</strong></td><td>66.5</td></tr><tr><td>MPCK-Means</td><td>0.078</td><td>0.017</td><td>29.27</td><td>159.9</td><td><strong>0.406</strong></td><td>0.019</td><td>26.98</td><td>61.2</td></tr></tbody></table>

### 4.2.3 Relaxing Constraints

We address the research question **RQ3** by comparing IAC with existing constrained clustering methods for constraint relaxation. We selected two methods: *Pairwise Constrained K-Means* (PCK-Means) [3] allows the violation of ML and CL constraints, and *Metric PCK-Means* (MPCK-Means) [6] combines PCK-Means with distance-metric learning [36]. As the alternatives only support pairwise constraints, we compare performance for problems that involve only ML/CL constraints. To that end, we consider two settings: in the first one, we generate 950 constraints at random based on the ground truth and add 50 conflicting constraints so that some constraints must be relaxed to solve the problem; in the second, we generate 1000 constraints at random without any explicit conflict constraint and set the satisfaction rate $\delta$ of IAC to the mean satisfaction rate of PCK-Means and MPCK-Means. We run two variants of IAC: IAC+Anchors ($\alpha = 20\%$, $\beta = 100\%$) and IAC+Anchors+SI ($\alpha = 20\%$, $\beta = 30\%$). Table 2 shows the performance of the different approaches for the mk2 dataset, measured by ARI, runtime and number of constraints relaxed. Each value in the table is the average of 90 runs with different sets of constraints. As previously, for each run, we compute an initial KMeans partition and run our CP model only once for each set of constraints. Results show that in the presence of conflicts, our method violates fewer constraints than the other methods. Furthermore, the runtimes of IAC are comparable to those of PCK-Means. In contrast, MPCK-Means is significantly more expensive. In terms of clustering accuracy (measured by ARI), our approach clearly outperforms the compared to the alternatives. When imposing a satisfaction rate to the model, PCK-Means and MPCK-Means achieve comparable or better quality than IAC, but similarity stays low. However, IAC achieves again better performance in number of constraints relaxed.

### 4.2.4 Comparing IAC with alternatives in the incremental setting

In this section, we address research questions **RQ4** and **RQ5**, and conduct experiments to evaluate the interest and performance of our IAC framework for active constraint section context. NPU can be used in combination with any semi-supervised clustering algorithm, we use the same ones as in the previous section, including COPK-Means algorithm [33]. This leads to several combinations, and for each combination we perform 10 iterations of selection-modification loop. For each iteration, we use the current partition to select a batch of 10 queries with NPU and at random, get feedback from the user, and perform either a reclustering or MCM. In this experiment, queries correspond to pairwise constraints. In light of the results of Section 4.2.1, we choose the compromise configuration between quality and similarity ($\alpha = 20\%$, $\beta = 30\%$) for this experiment.
10:14 Incremental Constrained Clustering by Minimal Weighted Modification

<figure><img src="image_88.png" /><figcaption>Figure 7 ARI scores of the partition at each iteration of selection-modification, compared to ground truth (left) or to the previous partition (right), for the glass dataset. All ARI scores are the mean and 95% confidence interval over 90 runs. Higher is better.</figcaption></figure>

**Clustering quality comparison.** Figure 7 shows the evolution of the ARI scores over 10 iterations of incremental and active clustering modification for the **glass** dataset. IAC+NPU produces increasingly better clusterings as more iterations are given (cf. Fig. 7a), while keeping a high similarity throughout the iterations (cf. Fig. 7b). None of the competitors produces a clustering with a high ARI. Interestingly, both PCK-Means and COPK-Means with NPU are able to find good clusterings while MPCK-Means is not, even after a relatively large number of iterations. However, the behaviour of the competitors are more chaotic in terms of similarity. We observe similar results for the other datasets (due to lack of space, all other results are available via our link in the summary). These results also show that methods with a random selection of constraints produce typically worse results. We validate these observations by Bayesian comparison w.r.t. AUBC<sub>quality</sub> and AUBC<sub>similarity</sub> values for each metric. Fig. 8 show that IAC has a high probability to perform better than the alternatives on ARI; we have similar results for AMI and FMI. Interestingly, using NPU leads to better similarity (see Fig. 9). This suggests that the use of an active constraint selection strategy brings another advantage besides improving the quality of the clustering. However, in an online context, runtime is particularly important as it requires user interaction and selecting the next query can be very costly, superseding the time taken for modification. For the biggest datasets (Letters and MNIST), only methods using random selection finish before timeout. We can conclude from these results that our model for minimal clustering modification is effectively a better way for active incremental constrained clustering compared to the naive approach to incrementality.

<figure><img src="image_89.png" /><figcaption>Figure 8 Bayesian comparisons with IAC using NPU w.r.t. AUBC<sub>quality</sub> values for ARI.</figcaption></figure>
<figure><img src="image_90.png" /><figcaption>**Figure 9** Bayesian comparison of IAC with or without NPU w.r.t AUBC<sub>similarity</sub> for all metrics.</figcaption></figure>

**Runtime comparison.** In Fig. 10, the runtime of modification (reclustering or MWCM) is shown for each dataset. The runtime of our model is comparable to those of the competitors, although it seems to have better scaling on the largest datasets. Empirically, IAC is an order of magnitude faster on yeast, statlog and letters datasets that have a large number of clusters. It is also noteworthy that MPCK-Means is much slower than other methods due to metric learning, yet this increase does not translate into better quality or similarity than IAC. The limiting factor for the modification step of IAC is the computation of anchors and super-instances, whose scaling is worse than solving the COP in itself.

<figure><img src="image_91.png" /><figcaption>**Figure 10** Evolution of the runtime of tested methods on our benchmark datasets, in seconds (log scale). Methods using NPU are hatched.</figcaption></figure>

## 4.3 Tree Cut Data

**Introducing the data.** Our case study for research question **RQ6** concerns the analysis of satellite image time series (SITS) composed of 11 images of dimensions 724 × 337 of a zone of the Vosges mountains in eastern France, taken irregularly on the span of 3 years from 2016 to 2018. Each pixel is associated to a series of NDVI (*Normalized Difference Vegetation Index*) values denoting the level of vegetation at each timestamp. At our disposition are labels that separate the SITS into three classes: vegetation, artificial structures, and tree cut zones. This last class has several properties: it was precisely labeled by domain experts, whereas the two other classes are more approximately defined. It is also a very small class as shown in Fig. 11a, containing only 639 instances (less than 0.3% of the data), which makes it hard to detect with unsupervised learning. Image-wise, 10 zones have been identified as places where trees have been cut within the time interval. Lastly, the evolution of these
## 10:16 Incremental Constrained Clustering by Minimal Weighted Modification

zones (a sharp decrease of NDVI value followed by a slow return to normal) is similar to that of field harvesting or grassland mowing, which complicates the problem further. In these conditions, expert intervention is paramount.

**Problem definition.** A set of 179 ML/CL constraints has been collected in [21] from domain experts, focused on the two largest tree cut zones (204 and 147 instances, i.e. more than half of tree cuts) as shown in Fig.11b. We define the problem as recovering these areas with binary constrained clustering. Following [22], we clustered the dataset with K-Means and $K = 15$, only retaining the cluster covering the areas the most as the “positive” cluster of our problem. In Fig. 11c, this cluster is displayed in colors, while the “negative” cluster is composed of all pixels not colored. We then used this binary partition as input of IAC, and selected the unsatisfied user constraints in the partition to improve it. In our experiments, 79 constraints were unsatisfied. IAC was set to iterate until all constraints are satisfied. We set $\beta$ to 100% as the large dataset size means that super-instance computation takes hours to complete, which is not compatible with a real life setting.

<figure><img src="image_92.png" /><figcaption>Figure 11 Some views of the use case; pictures (c) and (d) show the “positive” cluster, before and after modification. Highlighted therein are the true positives (green), false negatives (yellow), and false positives (purple). Best viewed in colors.</figcaption></figure>

**Results.** The modified clear cut cluster is displayed in Fig. 11d. The recovering of tree cut can be observed as the green zone of true positives enlarges in this figure : the left area progressed from 37 to 92 true positives, covering almost half the area. The right area gained 10 true positives. The modification was made within 22 seconds.

## 5 Conclusion

We have developed IAC, a framework for clustering modification that can be used in an incremental setting where an expert iteratively adds constraints, either manually or using an active method. A CP model for minimal weighted clustering modification ensures that the general cluster structure is preserved to maintain some continuity between iterations, as
shown in experiments on reference datasets and on a real use case. It can also efficiently exploit an active query strategy to converge faster, and handle contradictory constraints that the user may give as input. The runtime of IAC is dependent of the constraint selection step, which requires further experiments with more active methods and/or to develop a new one suited for the incremental setting. It would also be interesting to explore the use of multiple CP models for modification, such as [20] for minimal modification with cluster-level constraints. Lastly, there remain open questions about the potential reuse of relaxed constraints at a later iteration: What constraints to choose ? When to propose them to the user ? The conception of a strategy answering these interrogations is worth considering.

References

1 Arindam Banerjee and Joydeep Ghosh. Scalable clustering algorithms with balancing constraints. *Data Min. Knowl. Discov.*, 13(3):365–395, 2006.

2 Sugato Basu, Arindam Banerjee, and Raymond J. Mooney. Active semi-supervision for pairwise constrained clustering. In *Proceedings of the Fourth SIAM International Conference on Data Mining*, pages 333–344. SIAM, 2004. doi:10.1137/1.9781611972740.31.

3 Sugato Basu, Arindam Banerjee, and Raymond J. Mooney. Active Semi-Supervision for Pairwise Constrained Clustering. In *ICDM*, pages 333–344, 2004.

4 Alessio Benavoli, Giorgio Corani, Janez Demšar, and Marco Zaffalon. Time for a change: A tutorial for comparing multiple classifiers through Bayesian analysis, 2017.

5 Christian Bessière, Emmanuel Hébrard, Brahim Hnich, Zeynep Kiziltan, and Toby Walsh. Filtering Algorithms for the NValue Constraint. *Constraints*, 11(4):271–293, 2006.

6 Mikhail Bilenko, Sugato Basu, and Raymond J. Mooney. Integrating constraints and metric learning in semi-supervised clustering. In *Proceedings of the 21st International Conference on Machine Learning*, pages 11–18, 2004.

7 David Cohn, Rich Caruana, and Andrew Mccallum. Semi-Supervised Clustering with User Feedback. Technical report, Cornell University, 2001. doi:10.1201/9781584889977.ch2.

8 Giorgio Corani, Alessio Benavoli, Janez Demšar, Francesca Mangili, and Marco Zaffalon. Statistical comparison of classifiers through Bayesian hierarchical modelling. *Machine Learning*, 106(11):1817–1837, November 2017. doi:10.1007/s10994-017-5641-9.

9 Thi-Bich-Hanh Dao, Khanh-Chuong Duong, and Christel Vrain. Constrained clustering by constraint programming. *Artificial Intelligence*, 244:70–94, 2017.

10 Ian Davidson and S. S. Ravi. Agglomerative Hierarchical Clustering with Constraints: Theoretical and Empirical Results. In *Knowledge Discovery in Databases: PKDD 2005*, Lecture Notes in Computer Science, pages 59–70, 2005.

11 Ian Davidson, S. S. Ravi, and Martin Ester. Efficient incremental constrained clustering. In *Proceedings of the 13th ACM SIGKDD international conference on Knowledge discovery and data mining*, pages 240–249, 2007.

12 Ian Davidson, S. S. Ravi, and Leonid Shamis. A SAT-based Framework for Efficient Constrained Clustering. In *Proceedings of the 2010 SIAM International Conference on Data Mining*, pages 94–105, 2010.

13 E. B. Fowlkes and C. L. Mallows. A Method for Comparing Two Hierarchical Clusterings. *Journal of the American Statistical Association*, 78(383):553–569, 1983.

14 Marek Gagolewski. A Framework for Benchmarking Clustering Algorithms, 2022.

15 Teofilo F. Gonzalez. Clustering to minimize the maximum intercluster distance. *Theoretical Computer Science*, 38:293–306, 1985.

16 Germán González-Almagro, Daniel Peralta, Eli De Poorter, José-Ramón Cano, and Salvador García. Semi-Supervised Constrained Clustering: An In-Depth Overview, Ranked Taxonomy and Future Research Directions, 2023.
10:18 Incremental Constrained Clustering by Minimal Weighted Modification

17 Mathieu Guilbert, Christel Vrain, Thi-Bich-Hanh Dao, and Marcilio C. P. de Souto. Anchored Constrained Clustering Ensemble. In *International Joint Conference on Neural Networks*, IJCNN, 2022.

18 Tias Guns. Increasing modeling language convenience with a universal n-dimensional array, cppy as python-embedded example. In *Modref*, volume 19, 2019.

19 Lawrence Hubert and Phipps Arabie. Comparing partitions. *Journal of Classification*, 2(1):193–218, 1985.

20 Chia-Tung Kuo, S. S. Ravi, Thi-Bich-Hanh Dao, Christel Vrain, and Ian Davidson. A framework for minimal clustering modification via constraint programming. In *Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence*, AAAI’17, pages 1389–1395, 2017.

21 Baptiste Lafabregue, Pierre Gançarski, Jonathan Weber, and Germain Forestier. Incremental constrained clustering with application to remote sensing images time series. In *2022 IEEE International Conference on Data Mining Workshops (ICDMW)*, 2022.

22 Thomas Lampert, Baptiste Lafabregue, Thi-Bich-Hanh Dao, Nicolas Serrette, Christel Vrain, and Pierre Gancarski. Constrained Distance-Based Clustering for Satellite Image Time-Series. *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, 12(11):4606–4621, 2019.

23 Eric Yi Liu, Zhaojun Zhang, and Wei Wang. Clustering with relative constraints. In *Proceedings of the 17th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*, KDD ’11, pages 947–955, New York, NY, USA, August 2011. Association for Computing Machinery. doi:10.1145/2020408.2020564.

24 James MacQueen. Some Methods For Classification And Analysis Of Multivariate Observations. In *Proceedings of the fifth Berkeley symposium on mathematical statistics and probability*, volume 1, pages 281–297, 1967.

25 Pavan Kumar Mallapragada, Rong Jin, and Anil K. Jain. Active query selection for semi-supervised clustering. In *19th International Conference on Pattern Recognition*, pages 1–4, 2008.

26 Logan Adam Mitchell. INCREMENT - Interactive Cluster Refinement. PhD thesis, Brigham Young University, 2016.

27 Nguyen-Viet-Dung Nghiem, Christel Vrain, and Thi-Bich-Hanh Dao. Knowledge integration in deep clustering. In *Machine Learning and Knowledge Discovery in Databases - European Conference, ECML PKDD 2022, Proceedings, Part I*, volume 13713 of *Lecture Notes in Computer Science*, pages 174–190. Springer, 2022. doi:10.1007/978-3-031-26387-3_11.

28 Nguyen-Viet-Dung Nghiem, Christel Vrain, Thi-Bich-Hanh Dao, and Ian Davidson. Constrained Clustering via Post-processing. In *Discovery Science*, Lecture Notes in Computer Science, pages 53–67, 2020.

29 Abdelkader Ouali, Samir Loudni, Yahia Lebbah, Patrice Boizumault, Albrecht Zimmermann, and Lakhdar Loukil. Efficiently finding conceptual clustering models with integer linear programming. In *Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence*, IJCAI’16, pages 647–654, 2016.

30 Alfred Ultsch and Jörn Lötsch. The Fundamental Clustering and Projection Suite (FCPS): A Dataset Collection to Test the Performance of Clustering and Data Projection Algorithms. *Data*, 5(1):13, 2020.

31 Nguyen Xuan Vinh, Julien Epps, and James Bailey. Information theoretic measures for clusterings comparison: Is a correction for chance necessary? In *Proceedings of the 26th Annual International Conference on Machine Learning*, 2009.

32 Kiri Wagstaff and Claire Cardie. Clustering with instance-level constraints. In *Proceedings of the Seventeenth International Conference on Machine Learning*, ICML ’00, pages 1103–1110, San Francisco, CA, USA, 2000. Morgan Kaufmann Publishers Inc.

33 Kiri Wagstaff, Claire Cardie, Seth Rogers, and Stefan Schrödl. Constrained k-means clustering with background knowledge. In *Proceedings of the Eighteenth International Conference on Machine Learning (ICML 2001)*, pages 577–584, 2001.
34 Kiri L. Wagstaff. Value, Cost, and Sharing: Open Issues in Constrained Clustering. In *Knowledge Discovery in Inductive Databases*, pages 1–10, 2007.

35 Xiang Wang and Ian Davidson. Flexible constrained spectral clustering. In *Proceedings of the 16th ACM SIGKDD international conference on Knowledge discovery and data mining*, KDD '10, pages 563–572, 2010.

36 Eric Xing, Michael Jordan, Stuart J Russell, and Andrew Ng. Distance metric learning with application to clustering with side-information. In *Advances in Neural Information Processing Systems*, volume 15. MIT Press, 2002.

37 Sicheng Xiong, Javad Azimi, and Xiaoli Z. Fern. Active Learning of Constraints for Semi-Supervised Clustering. *IEEE Trans. on Knowledge and Data Engineering*, 26(1):43–54, 2014.

38 Xueying Zhan, Huan Liu, Qing Li, and Antoni B. Chan. A Comparative Survey: Benchmarking for Pool-based Active Learning. In *Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21*, volume 5, pages 4679–4686, August 2021.
A Super-instances generation

This is an example case of generating super-instances and splitting to prevent emergent
conflicts. In Fig. 12a, the clusters of the partition in Fig. 2 are divided into super-instances.
However, some super-instances contain multiple constrained instances, e.g. the green one in
the bottom left, which could lead to conflicts. In Fig. 12b, these super-instances have been
splitted.

<figure><img src="image_93.png" /><figcaption>(a) Result of generating super-instances through *complete-link* hierarchical clustering on every cluster.</figcaption></figure>

<figure><img src="image_94.png" /><figcaption>(b) Final super-instances after splitting. Greyed instances are unconstrained and not used in the CSP.</figcaption></figure>

**Figure 12** Exemple preprocessing for super-instance generation on lsun dataset. Instances sharing a color are represented by the same super-instance. ML constraints are in red, CL constraints in dashed blue.
# B Incremental NPU

This is the variant of NPU we use in the selection step of IAC. In Algorithm 2, the main modifications are the removal of the reclustering step of the original NPU, and the output of the set of constraints obtained from the queries for the modification step. Algorithm 3 is unchanged and is shown to give a complete view of the algorithm.

**Algorithm 2 Incremental NPU.**

**Input** : Dataset $\mathcal{D}$, partition $\mathcal{P}$, oracle
**Output** : constraint set $\mathcal{C}$
1: $\mathcal{C} \leftarrow \emptyset$; $l \leftarrow 1$; $\mathcal{N} \leftarrow N_1$ | $N_1 = \{random(\mathcal{D})\}$
2: $x^* \leftarrow MostInformative(\mathcal{D}, \mathcal{P}, \mathcal{N})$
3: **for** each $N_i \in \mathcal{N}$ in decreasing order of $P(x^* \in N_i)$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;Query $x^*$ against any $x_i \in N_i$ to the oracle
5: &nbsp;&nbsp;&nbsp;&nbsp;**if** ($x^*, x_i, ML$) **then**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mathcal{C} \leftarrow (x^*, x_i, ML)$
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$N_i = N_i \cup x^*$$
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;break
9: &nbsp;&nbsp;&nbsp;&nbsp;**else**
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\mathcal{C} \leftarrow (x^*, x_i, CL)$$
11: &nbsp;&nbsp;&nbsp;&nbsp;**if** no ML is returned **then**
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$l \leftarrow l + +$$; $N_l = x^*$; $\mathcal{N} \leftarrow \mathcal{N} \cup N_l$
13: &nbsp;&nbsp;&nbsp;&nbsp;**return** $\mathcal{C}$

**Algorithm 3 MostInformative.**

**Input** : Dataset $\mathcal{D}$, partition $\mathcal{P}$, set of neighborhoods $\mathcal{N}$
**Output** : most informative data point $x^*$
1: Learn a random forest classifier using $\mathcal{P}$ as labels
2: Compute the similarity matrix $M$ s.t. $M[i, j]$ is the number of leaves where $i$ and $j$ are together normalized by the number of trees of the RF
3: **for** each $x \in \mathcal{U} = \mathcal{D} \setminus \mathcal{N}$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;**for** $i = 1$ to $l$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\frac{1}{|N_i|} \sum_{x_j \in N_i} M(x, x_j)$$
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$p(x \in N_i) = \frac{\sum_{p=1}^l \frac{1}{|N_p|} \sum_{x_j \in N_p} M(x, x_j)}{\sum_{p=1}^l \frac{1}{|N_p|} \sum_{x_j \in N_p} M(x, x_j)}$$
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$H(\mathcal{N}|x) = - \sum_{i=1}^l p(x \in N_i) \log_2 p(x \in N_i)$$
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$E(x) = \sum_{i=1}^l i * p(x \in N_i)$$
9: &nbsp;&nbsp;&nbsp;&nbsp;**return** $\arg\max_{x \in \mathcal{U}} \frac{H(\mathcal{N}|x)}{E(x)}$
## C Modifications and generalization

For each modified instance (or super-instance), we store its initial cluster membership and its new cluster membership. This allows the framework to keep a history of modifications and to easily retrieve the partition at any given iteration. Considering generalization, we also keep track of the composition of each constrained super-instance. When the COP produces a solution, Algorithm 4 transmits the modifications from the super-instance to the real data.

Algorithm 4 APPLYMODIFICATIONS.

Input : dataset $\mathcal{X}$, super-instances $S$, modifications $\mathcal{M}$, partition $\mathcal{P}$

Output : modified partition $\mathcal{P}'$

1: $\mathcal{P}' \leftarrow \mathcal{P}$

2: for each $sp \in S$ do

3: points $\leftarrow \{x \in \mathcal{X} | x \in sp\}$

4: for each $p \in points$ do

5: Update the membership of $p$ in $\mathcal{P}'$ with the corresponding value in $\mathcal{M}$
return $\mathcal{P}'$

## D Bayesian pairwise comparison of IAC configurations

Table 3 Number of wins for each configuration $(\alpha, \beta)$ using a Bayesian hierarchical model. Values in parentheses indicate the number of cases where the probability that a configuration has a higher score is greater than 95%. Values of $\alpha$ and $\beta$ are in percentage (%).

<table><thead><tr><th rowspan="2">(α, β)</th><th colspan="3">AUBC<sub>quality</sub></th><th colspan="3">AUBC<sub>similarity</sub></th></tr><tr><th>ARI</th><th>AMI</th><th>FMI</th><th>ARI</th><th>AMI</th><th>FMI</th></tr></thead><tbody><tr><td>(0, 10)</td><td>9 (1)</td><td>10 (1)</td><td>9 (1)</td><td>0</td><td>0</td><td>0</td></tr><tr><td>(0, 30)</td><td>2 (0)</td><td>1 (0)</td><td>1 (0)</td><td>3 (1)</td><td>2 (0)</td><td>1 (1)</td></tr><tr><td>(0, 50)</td><td>0</td><td>0</td><td>0</td><td>6 (3)</td><td>6 (2)</td><td>4 (2)</td></tr><tr><td>(0, 100)</td><td>0</td><td>4 (1)</td><td>3</td><td>9 (9)</td><td>9 (9)</td><td>9 (9)</td></tr><tr><td>(5, 10)</td><td>10 (3)</td><td>10 (3)</td><td>10 (3)</td><td>1 (0)</td><td>1 (0)</td><td>1 (0)</td></tr><tr><td>(5, 30)</td><td>7 (3)</td><td>5 (3)</td><td>7 (3)</td><td>4 (2)</td><td>4 (1)</td><td>3 (1)</td></tr><tr><td>(5, 50)</td><td>1 (0)</td><td>2 (0)</td><td>1 (0)</td><td>6 (4)</td><td>6 (3)</td><td>4 (2)</td></tr><tr><td>(5, 100)</td><td>5 (2)</td><td>4 (0)</td><td>5 (1)</td><td>9 (9)</td><td>9 (9)</td><td>9 (9)</td></tr><tr><td>(20, 10)</td><td>10 (7)</td><td>10 (6)</td><td>10 (6)</td><td>1 (1)</td><td>2 (1)</td><td>1 (0)</td></tr><tr><td>(20, 30)</td><td>8 (3)</td><td>6 (3)</td><td>7 (3)</td><td>4 (2)</td><td>4 (2)</td><td>3 (1)</td></tr><tr><td>(20, 50)</td><td>4 (1)</td><td>3 (1)</td><td>4 (1)</td><td>6 (5)</td><td>6 (4)</td><td>6 (4)</td></tr><tr><td>(20, 100)</td><td>5 (2)</td><td>7 (3)</td><td>5 (1)</td><td>9 (9)</td><td>9 (9)</td><td>9 (9)</td></tr></tbody></table>
# Simplifying Step-Wise Explanation Sequences

Ignace Bleukx ⬒ ID

<figure><img src="image_95.png" /></figure>

<figure><img src="image_96.png" /></figure>

DTAI, KU Leuven, Belgium

Jo Devriendt ⬒ ID

<figure><img src="image_97.png" /></figure>

<figure><img src="image_98.png" /></figure>

DTAI, KU Leuven, Belgium

Emilio Gamba ⬒ ID

<figure><img src="image_99.png" /></figure>

<figure><img src="image_100.png" /></figure>

Data Analytics Lab, VUB, Brussels, Belgium

DTAI, KU Leuven, Belgium

Bart Bogaerts ⬒ ID

<figure><img src="image_101.png" /></figure>

<figure><img src="image_102.png" /></figure>

Artificial Intelligence Lab, VUB, Brussels, Belgium

Tias Guns ⬒ ID

<figure><img src="image_103.png" /></figure>

<figure><img src="image_104.png" /></figure>

DTAI, KU Leuven, Belgium

Data Analytics Lab, VUB, Brussels, Belgium

## Abstract

Explaining constraint programs is useful for debugging an unsatisfiable program, to understand why a given solution is optimal, or to understand how to find a unique solution. A recently proposed framework for explaining constraint programs works well to explain the unique solution to a problem step by step. It can also be used to step-wise explain why a model is unsatisfiable, but this may create redundant steps and introduce superfluous information into the explanation sequence. This paper proposes methods to simplify a (step-wise) explanation sequence, to generate simple steps that together form a short, interpretable sequence. We propose an algorithm to greedily construct an initial sequence and two filtering algorithms that eliminate redundant steps and unnecessarily complex parts of explanation sequences. Experiments on diverse benchmark instances show that our techniques can significantly simplify step-wise explanation sequences.

**2012 ACM Subject Classification** Theory of computation → Constraint and logic programming

**Keywords and phrases** explanation, deduction, constraint programming, propagation

**Digital Object Identifier** 10.4230/LIPIcs.CP.2023.11

**Supplementary Material** Software (Source code): https://github.com/ML-KULeuven/SimplifySeq
archived at swh:1:dir:a0de0a65a6c2c3a69d84a1de4047a360d8854195

**Funding** This research was partly funded by the Flemish Government (AI Research Program), the Research Foundation - Flanders (FWO) project G070521N; the European Research Council (ERC) under the EU Horizon 2020 research and innovation programme (Grant No 101002802, CHAT-Opt) and the European Union's Horizon 2020 research and innovation programme under grant agreement No 101070149, project Tuples.

## 1 Introduction

As AI agents become more expressive and powerful, a growing need arises for *explainable AI*: methods that can explain the decisions of an automated agent. Such explanations can be needed for legal reasons,<fnref n="1" /> but are also essential to provide trust to users. Considerable research has gone into developing explanation techniques for black-box machine learning methods [17, 27], including techniques based on formal models [33].

<figure><img src="image_105.png" /></figure>

<footnote n="1">¹ E.g., the GDPR – https://gdpr.eu</footnote>
11:2 Simplifying Step-Wise Explanation Sequences

But also explaining formal models themselves, such as constraint programs and satisfiability problems [37], is of high importance. Although individual constraints typically have a clear meaning, their interaction can be highly non-trivial, which led to the development of explanation methods for constraint programs [21]. A prominent branch of *explainable constraint solving* is occupied with explaining why a set of constraints is unsatisfiable. Most of these methods [30, 32, 26, 28, 31, 25, 29] extract a *minimal unsatisfiable subset* (MUS): a (minimal) subset of the constraints that renders the problem unsatisfiable. Such a MUS provides a user with a (potentially large) subset of constraints that yield inconsistency. Recently, there has been work on guiding users with respect to *what* can be done to restore feasibility [22, 39], but there is a lack of tools to better explain *why* a problem is inconsistent.

In a sense, traces or proof-logs of a solver (as common in SAT [24] and recently also finding its way into richer formalisms [2, 18, 3, 8, 4]) provide an explanation of why a model is unsatisfiable. Still, they would quickly overwhelm a user, as it involves constraint reformulations, auxiliary variables, and branching decisions. Other solver-generated explanations can be extracted from highly effective *Lazy Clause Generation* solvers [34] combining constraint propagation and SAT solving. In this setting, every propagation is *explained* by adding a clause to the SAT solver. However, these kinds of explanations are by nature restricted to explaining the propagation of a single constraint. Instead, we start from *step-wise explanations* [7, 15] where each explanation step in a sequence refines a partial assignment, using a minimal set of *constraints* as well as *facts* derived in previous steps.

These step-wise approaches were developed in the context of explaining the unique solution of satisfiable logic puzzles. Each step explains why a certain variable in the unique solution took that specific value. An example is explaining how to solve a Sudoku-puzzle [13, 20] where each step derives the value of a cell in the solution. In this context, the number of explanation steps is at most the number of variables in the problem.

But this explanation framework can be applied more broadly to constraint satisfaction problems (CSPs), in particular also to those that are unsatisfiable. A step-wise explanation of an over-constraint model allows to *debug* it, by listing steps similar to a debugger for programming languages. Each step shows a (preferably small) subset of constraints causing the removal of allowed values in variables domains, up to where a conflict is derived.

In contrast to explanation sequences for satisfiable problems, in the unsatisfiable case, only a subset of the variables and derived values contribute to the conflict and should therefore be explained. Therefore, directly applying the step-wise explanations framework to this new setting results in overly complex explanations.

Furthermore, finding *the shortest* sequence of explanations is a hard problem. A recent paper [6] touches upon the complexity of finding the shortest sequence of arc-consistency propagations steps. In that setting, the goal is to explain the full result of the arc-consistency algorithm. In the general case, where other/stronger propagation algorithms are used to find propagation steps, the problem may become even harder.

In this paper, we do not consider finding *the best* explanation sequence. Instead, we investigate how to find good – interpreted here as short and with small steps – step-wise explanations in the context of explaining why a CSP is unsatisfiable. The proposed techniques also apply to explaining the objective value of an optimization problem, explaining the solution(s) of constraint problems, and can be of use in interactive configuration problems [23].

For this, we contribute the following:

1. For the first time, we consider the quality of explanation *sequences* as a whole;

2. We formalize the properties that good explanation sequences, and explanation steps from which they are built, should adhere to;
3. We propose a new normal form for explanation sequences;

4. In sections 5 and 6 we propose new algorithms to simplify sequences in this normal form;

5. We show our methods significantly simplify explanation steps *and* sequences compared to current approaches.

## 2 Preliminaries

We now formalize constraint satisfaction problems and step-wise explanation sequences as used throughout this paper.

### 2.1 Constraint Satisfaction Problems (CSPs)

► **Definition 1 (CSP).** A CSP is a triple $(\mathcal{X}, \mathcal{D}, \mathcal{C})$ [37] with

■ $\mathcal{X}$ a set of variables;

■ $\mathcal{D}$ a set of domains $D_x$ of allowed values for each variable $x$ of $\mathcal{X}$, i.e., $\mathcal{D} = \{D_x \mid x \in \mathcal{X}\}$;

■ $\mathcal{C}$ a set of constraints, each over a subset of the variables.

A *full assignment* to a CSP $(\mathcal{X}, \mathcal{D}, \mathcal{C})$ is a mapping such that each variable takes a value from its domain. A *constraint* is a function mapping full assignments to true or false, typically described by a formula (e.g., $x + y \ge 1$). A constraint is *satisfied* by a full assignment if it maps the constraint to true. A *solution* to a CSP is a full assignment satisfying all constraints in $\mathcal{C}$. A CSP is *unsatisfiable* if it has no solution. A set of constraints $\mathcal{S}$ is a *logical consequence* of another set $\mathcal{S}'$ if all solutions to $\mathcal{S}'$ are solutions to $\mathcal{S}$ – written as $\mathcal{S}' \models \mathcal{S}$. Constraints can be arranged in a *constraint graph* where nodes represent variables that are connected by an edge if they co-occur in a constraint.

For the remainder of this paper, we assume the set of variables and their domains are known. A *positive literal* is an equality $x = v$ and a *negative literal* is an inequality $x \neq v$ for variable $x$ and value $v$. Negative literals represent the constraint that a variable cannot be assigned to some value from its domain. We will often employ sets of negative literals, where we use $x \in \mathcal{R}$ as a shorthand for $\{x \neq v \mid v \in D_x \setminus \mathcal{R}\}$. For example, with $D_x = \{0, 1, 2, 3\}$, $x \in \{1, 2\}$ means the negative literals $\{x \neq 0, x \neq 3\}$. In other words, $x \in \mathcal{R}$ denotes a set of negative literals enforcing that $x$ can *only* take values remaining in $\mathcal{R}$. Note that a positive literal $x = v$ is equivalent to the set of negative literals $x \in \{v\}$. For the remainder of this paper, literals are assumed to be **negative** unless explicitly mentioned otherwise. With $\perp$ we denote the *trivial inconsistency*, i.e., the singleton set containing the literal *false*.

For simplicity, this paper only uses *integer domains* for variables. An integer domain is represented as either a finite range with a lower and upper bound, or an enumerated set. E.g., the range $[0..3]$ and the set $\{0, 1, 2, 3\}$ are identical. All ideas and algorithms in this paper apply to non-integer finite domains as well.

Given a set of constraints partitioned in *soft* constraints and a set of *hard* constraints, a *minimal unsatisfiable subset* (MUS) is a subset of the soft constraints that is unsatisfiable in conjunction with the hard constraints, and for which all strict subsets are satisfiable in conjunction with these hard constraints. MUS-calculation techniques are a well-studied research field and several algorithms for this exist [25, 32].

### 2.2 Step-wise explanations

The *step-wise explanation framework* was introduced in the context of first-order logic and Boolean satisfiability [10, 9, 7, 15]. We here reinterpret it from a finite-domain CSP perspective.
11:4 Simplifying Step-Wise Explanation Sequences

Given a CSP (*X*, *D*, *C*), an explanation step is a triple (*E*, *S*, *N*), where the input *E* and the output *N* are disjoint sets of literals, and *S* ⊆ *C* is the constraint subset, such that *E* ∪ *S* |= *N*. Informally, an explanation step consists of a subset of constraints which, together with some input literals, imply “new” output literals.

Given a CSP (*X*, *D*, *C*), a target set of literals *T* and a given set of literals *G*, an explanation sequence of length *n* from *G* to *T* is a sequence $\langle(\mathcal{E}_i, \mathcal{S}_i, \mathcal{N}_i)\rangle_{1\le i\le n}$ of explanation steps where $\mathcal{E}_i \subseteq \mathcal{G} \cup \bigcup_{1\le j<i} \mathcal{N}_j$ and $\mathcal{T} \subseteq \mathcal{G} \cup \bigcup_{1\le j\le n} \mathcal{N}_j$. Informally, an explanation sequence is a sequence of explanation steps where each step derives some new literals and can do so using $\mathcal{G}$ as well as previously derived outputs. Eventually, the sequence should derive all literals in $\mathcal{T}$. As a whole, an explanation sequence explains how a target set is entailed by the union of a given set and the constraints of a CSP.

For satisfiable CSPs with a unique solution, we can let $\mathcal{T}$ be a unique variables assignment.
In general, we can set $\mathcal{T}$ to the intersection of all solutions. An explanation sequence where
$\mathcal{G} = \emptyset$ and $\mathcal{T}$ is inconsistent explains a CSP's unsatisfiability.

We say an explanation step *derives* or *explains* a literal *l* if *l* ∈ *N*, and a sequence derives or explains a literal if one of its steps does.

The maximal output for an input $\mathcal{E}$ and constraint subset $\mathcal{S}$ is the set of all literals implied by $\mathcal{E} \cup \mathcal{S}$, minus $\mathcal{E}$. The maximal output can be calculated using a FULLPROPAGATE algorithm which finds the set of literals that are true in all solutions to the constraints. Such a FULLPROPAGATE function is implemented in multiple systems such as the Answer Set Programming system Clasp [16], the IDP system [11], and more recently in the pseudo-Boolean solver Exact [12, 14].

**► Definition 2 (Maximal sequence).** An explanation step $(\mathcal{E}_i, \mathcal{S}_i, \mathcal{N}_i)$ in a sequence is maximal if (1) $\mathcal{E}_i$ is the union of all previously derived literals and the given set, i.e., $\mathcal{E}_i = \mathcal{G} \cup \bigcup_{j<i} \mathcal{N}_j$, and (2) $\mathcal{N}_i$ is the maximal output of $\mathcal{E}_i$ and $\mathcal{S}_i$. An explanation sequence where all steps are maximal is a maximal sequence.

Given a sequence of constraint sets $\langle \mathcal{S}_i \rangle_{1\le i\le n}$ and a given set $\mathcal{G}$, the maximal step-
wise explanation sequence is the unique sequence $\langle (\mathcal{E}_i, \mathcal{S}_i, \mathcal{N}_i) \rangle_{1\le i\le n}$ where all explanation
steps are maximal. Clearly, maximal sequences contain much more input/output literals
than a user would care about. Moreover, calculating the maximal output of a step using
FULLPROPAGATE is an expensive operation. In the general case, any sound propagation
algorithm can be used to calculate the output of a step, but using a maximal one will provide
us with a useful normal form.

**► Example 3.** Consider the following unsatisfiable CSP and explanation sequences which we will use as a running example:

$\mathcal{X} = \{x, y, z, v, w\}$

$\mathcal{D} = \{D_x = D_y = [1..3], D_z = D_v = D_w = [0..3]\}$

$\mathcal{C} = \{x + y + z \ge 7, x + y + w \le 4, x \le v, v + z \le 3\}$

3 Greedy initial sequence construction

Naively, generating explanation sequences involves constructing sequences of explanation
steps that neatly match each other's input and output. In existing work [7, 15], explanation
sequences are constructed using an iterative loop that greedily searches for the best next
explanation step to add. Each individual step in the explanation sequence is optimal with
respect to an assumed cost function or heuristic. This cost function for example takes into
account the number of constraints and the number of input literals used for each step.
<table><caption>Table 1 Two explanation sequences for the same unsatisfiability for $\mathcal{G} = \emptyset$ and $\mathcal{T} = \perp$.</caption><thead><tr><th colspan="4">(a)</th></tr><tr><th><i>i</i></th><th><i>ε</i></th><th><i>S</i></th><th><i>N</i></th></tr></thead><tbody><tr><td>1</td><td>∅</td><td>x + y + z ≥ 7</td><td>z ≠ 0</td></tr><tr><td>2</td><td>∅</td><td>x + y + z ≥ 7<br>x + y + w ≤ 4</td><td>z ∈ [3..3]<br>w ≠ 3</td></tr><tr><td>3</td><td>z ∈ [3..3]</td><td>x ≤ v<br>v + z ≤ 3</td><td>⊥</td></tr></tbody></table>

<table><thead><tr><th colspan="4">(b)</th></tr><tr><th><i>i</i></th><th><i>ε</i></th><th><i>S</i></th><th><i>N</i></th></tr></thead><tbody><tr><td>1</td><td>∅</td><td>x + y + z ≥ 7<br>x + y + w ≤ 4</td><td>z ∈ [3..3]</td></tr><tr><td>2</td><td>∅</td><td>x ≤ v</td><td>v ∈ [1..3]</td></tr><tr><td>3</td><td>v ∈ [1..3]<br>z ∈ [3..3]</td><td>v + z ≤ 3</td><td>⊥</td></tr></tbody></table>

To find a set of constraints and literals that explains a new literal $n$, given a sequence of $i-1$ previous steps, we can extract an unsatisfiable subset from $C \cup G \cup \bigcup_{j<i} N_j \cup \{-n\}$. This is precisely what is done in the algorithms presented in [7], by enumerating subsets $S$ of increasing size and finding a small MUS in this way for each literal to explain.

Such a construction method works well for explaining unique solutions, where the target contains all consequences of the CSP. However, for arbitrary target set, a literal might be derivable with a simple step, but deriving that literal doesn't bring us any closer to explaining the target. Nevertheless, literals that do not appear in the target may well be useful to build other intermediate steps with. Balancing which literals (and their associated explanation steps) to include in the sequence, or to exclude, will be a crucial theme in the next sections.

We propose a greedy sequence construction algorithm inspired by [7]. Given a CSP $(\mathcal{X}, \mathcal{D}, \mathcal{C})$, a given set $\mathcal{G}$ and a target set $\mathcal{T}$, Algorithm 1 CONSTRUCT-GREEDY describes our construction algorithm. It keeps track of an input set $\mathcal{E}$ that contains the union of $\mathcal{G}$ and all literals explained so far. Using this set of input literals, it iteratively tests if constraint subsets $\mathcal{S} \subseteq \mathcal{C}$ of growing size $|\mathcal{S}| = i$ can produce an explanation step. For this, the algorithm calculates the maximal output $\mathcal{N}$ for $\mathcal{E} \cup \mathcal{S}$ for which $\mathcal{E} \cup \mathcal{S} \models \mathcal{N}$ and $\mathcal{N} \cap \mathcal{E} = \emptyset$ using a FULLPROPAGATE routine. If this output set $\mathcal{N}$ is not empty, i.e., a new literal could be derived, it adds the step $(\mathcal{E}, \mathcal{S}, \mathcal{N})$ to the sequence, extends $\mathcal{E}$ with $\mathcal{N}$ and resets the constraint subset size $i$. This process is repeated until the target has been explained ($\mathcal{T} \subseteq \mathcal{N}$) and the sequence is finished.

Notice the algorithm is guaranteed to terminate as $\mathcal{S}$ grows in every iteration. Eventually, $\mathcal{S}$ will be equal to $C$ if no smaller subset was able to derive a new literal. Of course, we here assume $\mathcal{G} \cup \mathcal{C} \models \mathcal{T}$.

**Example 4.** The table below shows a maximal explanation sequence with given set $\mathcal{G} = \emptyset$ and target set $\mathcal{T} = \perp$ for the following unsatisfiable CSP:

- $\mathcal{X} = \{p, q, r, s\}$

- $\mathcal{D} = \{D_p = D_q = D_r = D_s = [0..3]\}$

- $\mathcal{C} = \{p + q \le 1, q + 2r \le 4, 3s + p + r \le 1, \text{alldiff}(p, q, r, s)\}$

## 3.1 Properties

**Lemma 5.** For each step $(\mathcal{E}, \mathcal{S}, \mathcal{N})$ in a sequence generated by CONSTRUCT-GREEDY, and for each $\mathcal{S}' \subsetneq \mathcal{S}$, the maximal output for $\mathcal{S}'$ and $\mathcal{E}$ is empty.

Otherwise CONSTRUCT-GREEDY would have picked $\mathcal{S}'$ to form a step with $\mathcal{E}$. A direct consequence of this is that the sequences constructed by CONSTRUCT-GREEDY are *atomic*:

**Property 1 (Atomic).** An explanation step $(\mathcal{E}, \mathcal{S}, \mathcal{N})$ is atomic if no explanation sequence from $\mathcal{E}$ to $\mathcal{N}$ exists in which each step uses a strict subset of $\mathcal{S}$. An explanation sequence is atomic if all its steps are.
11:6 Simplifying Step-Wise Explanation Sequences

■ Algorithm 1 CONSTRUCT-GREEDY.

Input: CSP (*X*, *D*, *C*), given *G*, target *T*
1 *Seq* ← empty sequence
2 *E* ← *G*, *i* ← 0
3 while true do
4   **forall** *S* ⊆ *C* where |*S*| = *i* and *CONNECTED*(*S*) **do**
5     **N** ← FULLPROPAGATE(*E* ∪ *S*) \ *E*
6     **if** **N** ≠ ∅ **then**
7       **extend** *Seq* with (*E*, *S*, **N**)
8       *E* ← *E* ∪ **N**, *i* ← 0
9       **if** *T* ⊆ *E* **then**
10        **return** *Seq*
11      **break**
12   **i** ← *i* + 1

■ Table 2 A maximal sequence that could have been constructed by CONSTRUCT-GREEDY.

<table><thead><tr><th><i>i</i></th><th><i>E</i></th><th><i>S</i></th><th><i>N</i></th></tr></thead><tbody><tr><td>1</td><td>∅</td><td><i>p</i> + <i>q</i> ≤ 1</td><td><i>p</i> ∈ [0..1]<br><i>q</i> ∈ [0..1]</td></tr><tr><td>2</td><td><i>p</i> ∈ [0..1]<br><i>q</i> ∈ [0..1]</td><td><i>q</i> + 2<i>r</i> ≤ 4</td><td><i>r</i> ≠ 3</td></tr><tr><td>3</td><td><i>p</i> ∈ [0..1]<br><i>q</i> ∈ [0..1]<br><i>r</i> ≠ 3</td><td>3<i>s</i> + <i>p</i> + <i>r</i> ≤ 1</td><td><i>r</i> ≠ 2<br><i>s</i> ∈ [0..0]</td></tr><tr><td>4</td><td><i>p</i> ∈ [0..1]<br><i>q</i> ∈ [0..1]<br><i>r</i> ∈ [0..1]<br><i>s</i> ∈ [0..0]</td><td>alldiff(<i>p</i>, <i>q</i>, <i>r</i>, <i>s</i>)</td><td>⊥</td></tr></tbody></table>

The third step of the first sequence of Example 3 is not atomic: it can be split into two steps over $x \le v$ and $v + z \le 3$ each, as is done in the second sequence of the same example.

We believe this property is desirable for easy-to-understand explanation sequences as it requires each explanation step to be as small as possible. Note that explanation sequences generated by the literature [7, 15] also exhibit atomicity.

3.2 Efficiency optimizations

The two main bottlenecks in this algorithm are *enumerating large subsets of constraints* and *calculating the maximal output* $\mathcal{N}$ for a given $\mathcal{E}$ and $\mathcal{S}$. To improve runtime, we employ two efficiency optimizations. First, we cache calls to the FULLPROPAGATE routine. When a maximal output call for some $\mathcal{E}$ and $\mathcal{S}$ is repeated, the output is just taken from the cache. This is useful as the domains of many variables are unchanged between calls and the propagation only considers the literals from $\mathcal{E}$ with regard to the variables occurring in $\mathcal{S}$.

A second optimization inspects the constraint graph of each enumerated subset $\mathcal{S}$ using the CONNECTED function call on line 4 in the algorithm. This function checks if the constraint graph is connected. If not, the graph can be split up into two or more components $\mathcal{S}_c$ and the maximal output $\mathcal{N}$ for $\mathcal{S}$ and $\mathcal{E}$ is exactly the union of the maximal outputs $\mathcal{N}_i$ for each component $\mathcal{S}_c$ and input $\mathcal{E}$. By Lemma 5, all smaller $\mathcal{S}_c$ did not imply any new literal, so $\mathcal{S}$ cannot imply any either. Naturally, this means CONSTRUCT-GREEDY can skip subset $\mathcal{S}$.
## 3.3 Comparison to literature

The idea of CONSTRUCT-GREEDY is similar to the first literature approach [7]: iterate over increasingly larger subsets $S \subseteq C$ and check whether an explanation step $(\mathcal{E}, S, \mathcal{N})$ exists.

In contrast to [7], CONSTRUCT-GREEDY always greedily picks the first explanation step that can derive at least one unexplained literal to add to the explanation sequence, instead of generating multiple candidates and picking the optimal one under some quality function. The reasons for this difference are twofold: a simpler algorithm and less computational effort per step. Compared to logic puzzles, which are aimed to be solved by humans and contain at most hundreds of literals, for arbitrary CSPs the literals can easily grow into the millions (e.g., for large integer domains), so efficiency is more of a concern here. This is also the reason we do not employ the second literature approach [15], whose algorithm computes the optimal next explanation step according to a linear cost function directly, without explicit enumeration, by exploiting the hitting set duality between MUS's and correction subsets.

Still, by iterating over small constraint subsets first, CONSTRUCT-GREEDY adds steps that employ a small set of constraints, which is preferred to easily understand the explanation [7].

Note that CONSTRUCT-GREEDY constructs maximal sequences. This contrasts with the literature approaches which minimize the input literals for each step during construction. The motivation for this difference is that we first want to minimize the number of steps and the number of constraints for each step. The more literals are derived by a step, the more literals future steps can use as input. Thereby potentially needing fewer constraints to derive new literals. Moreover, by deriving more literals in a single step, the quicker we might reach the target set.

One can also minimize the number of input literals for each step in the sequence during post-processing of the sequence, thereby simulating the sequences produced by literature approaches in terms of input literals [7, 15]. In Section 7, this is precisely what is denoted as Filter simple.

# 4 Simple post-processing

In this section, we consider two straightforward ways to simplify explanation sequences requiring little to no computational effort.

## 4.1 Looseness

For an explanation sequence to be interpretable, we deem the length of the sequence to be an important metric. Each step in an explanation sequence requires effort by a user to process and understand. Naturally, this means no *loose* steps should be part of the sequence:

> **► Property 2 (Loose).** Given an explanation sequence $\langle (\mathcal{E}_i, S_i, N_i) \rangle_{1 \le i \le n} \text{ of } \mathcal{T} \text{ from } \mathcal{G},$ the step $(\mathcal{E}_j, S_j, N_j)$ is called *loose* if $\langle (\mathcal{E}_i, S_i, N_i) \rangle_{i \ne j}$ still forms an explanation sequence of $\mathcal{T}$ from $\mathcal{G}$. An explanation sequence is *loose* if one of its steps is *loose*.

A step is *loose* if it can be left out of the sequence. This is the case for instance when none of a step's output literals is used as input for later steps or is part of the target set.

The first step in the sequence 1a of Example 3 is *loose*, as it can just be removed without impacting the sequence. In sequence 1b of that same example, no *loose* steps are present as at least one output of each step is used later on.

Existing approaches [7, 15] focused mostly on explaining a unique solution of a CSP, where every step derives at least one yet unexplained literal of that unique solution. This means every step is guaranteed to derive a literal in the target set. However, [7] also describes