As each generalized AXp for $\mathfrak{T}$ represents a prime implicant of the decision function $\tau(\mathbf{x})$ computed over literals $\tilde{o}_{ij}$, the above corollary claims that there are no overlapping rules in the result DS $\mathcal{R}$. This contrasts with other modern approaches to DS inference, where rule overlap is known to be a problem [29, 22]. Note that this approach still suffers from another common issue of DS models: namely, if DS $\mathcal{R}$ is computed for the training data $\mathcal{E}$, there may still be instances in $\mathbb{F}$ uncovered by $\mathcal{R}$.

**► Example 10.** Consider the running example BT model shown in Figure 1b. Its compiled DS representation computed by Algorithm 2 is shown in Figure 1a. Observe that there is no rule overlap in the DS computed. In fact, as the DS is computed by taking into account feature space $\mathbb{F}$, it computes the same classification function as the original BT model. ▼

**Feature Sorting.** Intuitively, how general and hence how applicable a rule is depends on how frequently the features used in it appear in the training data $\mathcal{E}$ labeled with the target class. Thus, a simple heuristic to apply when extracting a rule for prediction $c = \tau(\mathbf{v})$ is to sort the initial state of $\mathcal{Y} = \{\tilde{o}_{ij} \mid i \in \mathcal{F}, j \in [h_i]\}$ based on how frequently the corresponding literals $\tilde{o}_{ij}$ apply in examples $\mathcal{E}$ labeled with $c$. This feature sorting represented by line 3 in Algorithm 1 in practice (according to our experiments) results in significantly more general rules and so overall smaller DSes.

**Anytime Property.** Most widely used reasoning-based algorithms to infer DSes provide a solution only if the computation is completed; otherwise, no decision set is reported. In contrast to these, the proposed approach is an *anytime* algorithm, i.e. it can return a *valid* DS $\mathcal{R}$ even though the compilation process is interrupted before all the coverage set instances $C$ are covered. Furthermore, it can generate a more comprehensive DS $\mathcal{R}$, which covers more instances as it keeps going, i.e. after we have covered $C \subseteq \mathbb{F}$ we can continue running the algorithm for the (unseen) instances of $\mathbb{F}$.

## 4.3 Post-Hoc Model Reduction

The compiled DS $\mathcal{R}$ can be large (in terms of either the number of rules or the total number of literals) since each generalized AXp $\mathcal{Y} \in \mathcal{R}$ may need a significant number of literals to explain a prediction made by BT $\mathfrak{T}$, or/and many rules are required to explain all instances of $C$. Once the target DS is obtained, we can apply post-hoc heuristic methods for reducing its size and so making it more interpretable. The methods below are in a way inspired by the optimization problems studied in [18, 50]. Although these ideas are applicable to any DS inference method once the result model is devised, they do not look necessary for standard DS inference algorithms as they minimize the model while training. On the contrary, no minimization is applied in the rule enumeration process described above and so post-hoc model reduction plays a vital role in our approach to reduce the size of final DS models.

**Reducing the Number of Rules.** Given a set of rules $\mathcal{R}$, we can compute a minimum subset $\mathcal{R}^* \subseteq \mathcal{R}$ that is still equivalent to the BT $\mathfrak{T}$ wrt. the coverage set $C$ using discrete optimization, e.g. integer-linear programming (ILP). Concretely, the approach aims at selecting the smallest-size subset $\mathcal{R}^* \subseteq \mathcal{R}$ that covers all instances in $C$, where $\mathcal{R}$ is the compiled DS from $\mathfrak{T}$. Here, the size of $\mathcal{R}^*$ is measured as the total number of literals used. This can be done by solving the following *set cover problem* [28]. Namely, for each rule $\mathcal{Y}_j \in \mathcal{R}$, we introduce a Boolean variable $u_j$ such that $u_j = 1$ iff $\mathcal{Y}_j$ is included in $\mathcal{R}^*$.
38:10 From Formal Boosted Tree Explanations to Interpretable Rule Sets

Additionally, a Boolean variable $y_{ij}$ is used to indicate that $\mathcal{Y}_j$ covers $e_i \in \mathcal{C}$. As a result, the weighted set cover problem for minimizing the total number of literals used is as follows:

minimize

$$ \sum_{j=1}^{|\mathcal{R}|} (|\mathcal{Y}_j| + 1) \cdot u_j \qquad (4) $$

subject to

$$ \forall_{i \in [n]} \sum_{j=1}^{|\mathcal{R}|} y_{ij} \cdot u_j \geq 1 \qquad (5) $$

**Reducing the Number of Literals.** Additionally, one can minimize the total number of literals used in the rules of $\mathcal{R}$. Given a rule $\mathcal{Y} \in \mathcal{R}$, this can be done either lexicographically by maximizing rule accuracy followed by size minimization, or by optimizing both, or trading off misclassifications for rule size – in either case, a single MaxSAT call per rule to minimize can be made. The intuition is that if a rule $\mathcal{Y}$ misclassifies $k$ instances then its optimized version $\mathcal{Y}^* \subseteq \mathcal{Y}$ should not result in many more misclassifications on training data $\mathcal{E}$. Recall that a rule misclassifies an instance $\mathbf{v}_k \in \mathcal{C}$ if it matches $\mathbf{v}_k$ but assigns it to a wrong class.

Inspired by [18], we introduce a Boolean variable $p_k$, which is true iff rule $\mathcal{Y}$ covers $\mathbf{v}_k$ – this holds if $\mathcal{Y}$ does not use any literals incompatible with $\mathbf{v}_k$. If $\mathcal{Y}_{\mathbf{v}_k} = \{\tilde{o}_{ij} | i \in \mathcal{F}, j \in [h_i]\}$ are all the literals compatible with $\mathbf{v}_k$ then this can be modeled with constraints

$$ \forall_{k \in [|\mathcal{C}|]} p_k \leftrightarrow \bigwedge_{l \in \mathcal{Y} \setminus \mathcal{Y}_{\mathbf{v}_k}} \neg l \qquad (6) $$

Furthermore, let rule $\mathcal{Y}$ predict $c \in \mathcal{K}$ and let $\mathcal{C}_\oplus \subseteq \mathcal{C}$ contain all instances labeled with any other class. Thus, we can apply the objective below when minimizing rule $\mathcal{Y}$:

$$ \sum_{l \in \mathcal{Y}} l + \sum_{k \in [|\mathcal{C}_\oplus|]} W \cdot p_k \qquad (7) $$

If $W$ is large enough, say $|\mathcal{C}| + 1$, this lexicographically minimizes misclassifications and then literals. If $W$ is small, e.g. $1/\lambda \cdot |\mathcal{C}|$, this trades off $\lambda \cdot |\mathcal{C}|$ misclassifications for one literal.

5 Experimental Results

This section compares the proposed approach with the state-of-the-art DS learning algorithms on a variety of publicly available datasets in terms of accuracy, scalability, model and explanation size. The experiments are performed on an Intel Xeon 8260 CPU running Ubuntu 20.04.2 LTS, with the time limit of 3600s and the memory limit of 8GByte. Our experiments contain two parts, namely, exhaustive BT compilation and training-set BT compilation.

**Prototype implementation.** A prototype of the compilation-based approach to generating DSes was developed as a set of Python scripts using $\mathcal{C} = \mathcal{E}$, hereinafter referred to as *cpl*. The implementation of BT compilation exploits [17] and, therefore, makes use of the RC2 MaxSAT solver [19].<fnref n="4" /> The BTs to be compiled are computed by XGBoost [5]; the number of trees per class in a BT model is 50 and the maximum depth of each tree is 3. Post-hoc

<footnote n="4">⁴ Real weights in the objective function are not conventionally supported by MaxSAT solvers; the only other solver to support real weights besides RC2 is LMHS [47].</footnote>
literal reduction is done again with RC2 [19]. Let $cpl_l$ denote the implementation applying lexicographic optimization while $cpl_{l\lambda_1}$ trades off model accuracy for the number of literals used, with $\lambda_1 = 0.005$. Let $cpl_r$ denote the implementation with post-hoc rule reduction applied using the Gurobi ILP solver [14]. The configuration with both post-hoc lexicographic optimization and rule reduction is denoted $cpl_{lr}$. Finally, the proposed approach applying exhaustive compilation $C = \mathbb{F}$ is referred to as $cpl_f$.

**Competition.** Our approach is compared against: *twostg* a two-stage MaxSAT approach [18] for DSes perfectly accurate on the training data; *opt* another MaxSAT approach [50] for perfectly accurate DSes; *sp*$_{λ_1}$ a sparse alternative to *opt* by the same authors (with $λ_1 = 0.005$) optimizing like $cpl_{lλ_1}$; *imli*$_1$ and *imli*$_{16}$ using MaxSAT-based IMLI [12] to minimize the number of literals given a predefined number of rules (we use 1 or 16); *ids* a state-of-the-art approach [29] based on smooth local search;<fnref n="5" /> *ripper* a popular heuristic DS algorithm *RIPPER* [8]; and *CN2* (referred to as *cn2*) another heuristic algorithm [7, 6].⁶

**Datasets.** For the evaluation, 59 publicly available datasets from UCI Machine Learning Repository [9] and Penn Machine Learning Benchmarks [41] are considered. We apply 5-fold cross validation, resulting in 295 pairs of training and test (unseen) data. For the sake of a fair comparison, the datasets used are preprocessed so that each original feature $i \in \mathcal{F}$ is replaced with a number of non-intersecting feature intervals $x_i < d_{ij}$ defined by the XGBoost model (see Section 2). This guarantees that all competitors tackle the same problem instances.

## 5.1 Exhaustive BT Compilation

The first experiment compares exhaustive compilation, where $C = \mathbb{F}$ is the entire feature space. This is impractical except for 6 small benchmarks.

**Results.** Here we compare $cpl_f$ with the competition in terms of accuracy, the total number of literals used and explanation size. We present the results as cactus plots showing the number of datasets that e.g. reach a certain accuracy, or finish in a certain runtime, for each method. These experimental results are shown in Figures 2 and 3 as well as the average results across folds are described in Table 2 where only the results of the datasets completely solved by compared competitors are presented. Note that $cpl_f$ is nowhere near as scalable as the approaches described in the later experiments, but it is the most accurate approach to creating DSes we are aware of.

**Test accuracy.** An instance is considered misclassified if either there exists a rule of a wrong class that covers it, or it is not covered by any rule of the correct class. Thus, the test accuracy in this paper is calculated as $\frac{n-q}{n}$, where n is the total number of instances in the test data and q is the total number of misclassified instances. If an approach fails to train a model within the time limit, we assume its accuracy to be 0% for this dataset.

<footnote n="5">⁵ Since the original implementation performs poorly [22], here we consider the new implementation of IDS [10], which is claimed to be orders of magnitude faster than the original implementation.

⁶ Note that since RIPPER and IMLI compute a single class only given the training data, both of these competitors are augmented with a default rule predicting a class (1) different from the target class and (2) represented by the majority of training instances. Other algorithms, including our approach, incorporate a default rule that assigns a class based on the majority class in the training instances.</footnote>
<figure><img src="image_352.png" /><figcaption>■ **Figure 2** Accuracy of exhaustive compilation. The *standard* interpretation of cactus plots is assumed, i.e. a plot sorts the datapoints for each method by the *y*-axis value, and then shows them in increasing order independently of other methods. Thus, the order of datasets/folds differs for different methods. Also, the order of datasets for the same method differs in different subplots.</figcaption></figure>

<figure><img src="image_353.png" /><figcaption>■ **Figure 3** Succinctness of exhaustive compilation.</figcaption></figure>

As can be seen in Figure 2b and Table 2, the best accuracy is achieved by BTs and *cpl*<sub>*f*</sub>. In fact, these models share the same accuracy (this is also confirmed in Figure 2a), which should not come as a surprise given that *cpl*<sub>*f*</sub> replicates the behavior of the BT in the entire feature space F (see Proposition 8).

**Model Complexity.** In general, complexity of a DS model can be measured by the total number of literals used in this DS. The total number of literals used in DS models is compared in Figure 3a and Table 2. Though the accuracy of DSes trained by *cpl*<sub>*f*</sub> outperforms the other competitors, these models are significantly larger, which is no surprise given that *cpl*<sub>*f*</sub> computes many more rules with no post-hoc reduction applied.

**Explanation size.** Explanation size is defined as the number of literals required to explain an instance.<fnref n="7" /> This is arguably more important than the model size, since it defines “how hard” it is to understand an individual explanation. A small DS model tends to provide

<footnote n="7">⁷ See [51] for details.</footnote>
<table><caption>Table 2 Accuracy, number of literals used, and explanation size across folds.</caption><thead><tr><th rowspan="3">Approach</th><th colspan="6">Dataset</th></tr><tr><th colspan="6">Accuracy (%)</th></tr><tr><th>cardiotocography</th><th>hayes-roth</th><th>iris</th><th>new-thyroid</th><th>orbit</th><th>zoo</th></tr></thead><tbody><tr><td><strong>cpl<sub>f</sub></strong></td><td><strong>100.0</strong></td><td><strong>84.38</strong></td><td><strong>96.0</strong></td><td><strong>96.74</strong></td><td><strong>99.66</strong></td><td><strong>96.0</strong></td></tr><tr><td><strong>sp<sub>λ<sub>1</sub></sub></strong></td><td>100.0</td><td>73.44</td><td>94.0</td><td>91.63</td><td>99.43</td><td>89.05</td></tr><tr><td><strong>opt</strong></td><td>100.0</td><td>70.63</td><td>93.33</td><td>91.63</td><td>99.54</td><td>93.05</td></tr><tr><td><strong>twostg</strong></td><td>100.0</td><td>71.25</td><td>92.67</td><td>92.09</td><td>99.54</td><td>91.1</td></tr><tr><td><strong>cn<sub>2</sub></strong></td><td>100.0</td><td>62.5</td><td>92.67</td><td>93.02</td><td>99.54</td><td>89.1</td></tr><tr><td><strong>ripper</strong></td><td>45.3</td><td>66.25</td><td>57.33</td><td>80.93</td><td>94.11</td><td>60.33</td></tr><tr><td><strong>ids</strong></td><td>27.23</td><td>43.75</td><td>58.67</td><td>76.28</td><td>85.29</td><td>40.62</td></tr><tr><td><strong>imli<sub>16</sub></strong></td><td>27.23</td><td>38.75</td><td>25.34</td><td>69.77</td><td>70.55</td><td>43.33</td></tr><tr><td><strong>imli<sub>1</sub></strong></td><td>45.3</td><td>39.37</td><td>32.67</td><td>26.98</td><td>8.93</td><td>60.33</td></tr><tr><td colspan="7" style="text-align: center;"><strong>Number of literals used</strong></td></tr><tr><td><strong>cpl<sub>f</sub></strong></td><td>3120.0</td><td>76.0</td><td>214.0</td><td>3614.2</td><td>729.8</td><td>1422</td></tr><tr><td><strong>sp<sub>λ<sub>1</sub></sub></strong></td><td>21.0</td><td>33.5</td><td>9.0</td><td>15.4</td><td>10.0</td><td>23.2</td></tr><tr><td><strong>opt</strong></td><td>21.0</td><td>63.6</td><td>19.4</td><td>23.0</td><td>11.8</td><td>30.0</td></tr><tr><td><strong>twostg</strong></td><td>21.0</td><td>64.2</td><td>19.8</td><td>22.6</td><td>11.8</td><td>29.8</td></tr><tr><td><strong>cn<sub>2</sub></strong></td><td>21.0</td><td>116.2</td><td>27.2</td><td>36.6</td><td>13.2</td><td>40.8</td></tr><tr><td><strong>ripper</strong></td><td><strong>3.0</strong></td><td>12.8</td><td>5.0</td><td>8.2</td><td>4.0</td><td><strong>3.0</strong></td></tr><tr><td><strong>ids</strong></td><td>21.0</td><td>21.6</td><td>19.8</td><td>20.0</td><td>25.0</td><td>14.2</td></tr><tr><td><strong>imli<sub>16</sub></strong></td><td>5.0</td><td><strong>2.2</strong></td><td>7.4</td><td>7.4</td><td>6.4</td><td>5.0</td></tr><tr><td><strong>imli<sub>1</sub></strong></td><td><strong>3.0</strong></td><td><strong>2.2</strong></td><td><strong>3.0</strong></td><td><strong>4.2</strong></td><td><strong>3.0</strong></td><td><strong>3.0</strong></td></tr><tr><td colspan="7" style="text-align: center;"><strong>Explanation size</strong></td></tr><tr><td><strong>cpl<sub>f</sub></strong></td><td>7.26</td><td>3.76</td><td>3.02</td><td>4.9</td><td>3.18</td><td>5.4</td></tr><tr><td><strong>sp<sub>λ<sub>1</sub></sub></strong></td><td><strong>2.0</strong></td><td>6.31</td><td>2.45</td><td>4.13</td><td>2.86</td><td>3.64</td></tr><tr><td><strong>opt</strong></td><td><strong>2.0</strong></td><td>5.41</td><td>2.76</td><td>4.3</td><td>2.94</td><td>2.96</td></tr><tr><td><strong>twostg</strong></td><td><strong>2.0</strong></td><td>5.4</td><td>2.87</td><td>4.23</td><td>2.94</td><td>3.33</td></tr><tr><td><strong>cn<sub>2</sub></strong></td><td><strong>2.0</strong></td><td>6.94</td><td>3.02</td><td>4.47</td><td>3.02</td><td>4.05</td></tr><tr><td><strong>ripper</strong></td><td>2.73</td><td>10.15</td><td>4.3</td><td>4.3</td><td>3.15</td><td>2.59</td></tr><tr><td><strong>ids</strong></td><td>16.08</td><td>18.23</td><td>13.06</td><td>7.74</td><td>6.23</td><td>9.28</td></tr><tr><td><strong>imli<sub>16</sub></strong></td><td><strong>2.0</strong></td><td><strong>2.2</strong></td><td><strong>2.1</strong></td><td><strong>1.97</strong></td><td><strong>2.8</strong></td><td><strong>2.46</strong></td></tr><tr><td><strong>imli<sub>1</sub></strong></td><td>2.18</td><td><strong>2.2</strong></td><td>3.0</td><td>4.0</td><td>3.0</td><td><strong>2.2</strong></td></tr></tbody></table>

compact explanations but it is not always accurate. As can be seen in Figure 3b and Table 2 and similar to the total number of literals used in DSes, $cpl_f$ requires more literals to explain an instance than all competitors except *ids*.

A crucial observation to make here is that we test explanation size for each of the test instances available. Although test data are meant to extrapolate the overall unseen data, such approximation of the unseen feature space is not ideal. As a result, there may be numerous instances in $\mathbb{F}$ *uncovered* by all the approaches but $cpl_f$, in which case it will be the *only* approach providing a user not only with a prediction but also with a succinct explanation of the prediction made.

## 5.2 BT Compilation Targeting Training Data

Compilation to cover the training set $C = \mathcal{E}$ is much more efficient, and the main usage we expect of our algorithms.

**Scalability.** Figure 4a depicts scalability of all selected algorithms on the 295 considered datasets. Note that runtime of our approach includes BT training time. The best performance is demonstrated by the proposed implementation, i.e. $cpl$ and $cpl_*$, $* \in \{l, r, lr, l\lambda_1\}$, where all selected datasets are solved within the time limit. This is not surprising since the approach is an anytime algorithm that can always return a valid DS. As for other competitors, the heuristic method *ripper* and the MaxSAT approaches $imli_1$ as well as $imli_{16}$ also solve all considered datasets. Next is the heuristic algorithm *cn2*, where 235 datasets are solved
<figure><img src="image_354.png" /><figcaption>Figure 4 Summary of experimental results when the competitors aim at training a DS given training data $\mathcal{E}$ (i.e. $C = \mathcal{E}$).</figcaption></figure>

within the 3600s time limit. Followed by *ids*, which solves 166 considered datasets. The two-stage MaxSAT approach *twostg* successfully addresses 130 datasets, while the other MaxSAT algorithm for perfect decision sets *opt* and its sparse alternative $sp_{\lambda_1}$ solve 65 and 63 datasets respectively.

**Test Accuracy.** The accuracy among the selected approaches is shown in Figure 4b. The average accuracy among all selected datasets for BTs is 77.34%, beating all DS approaches. The highest accuracy among DSes is achieved by all the configurations of the proposed approach, i.e. *cpl* and *cpl*<sub>*</sub>, where the average accuracy ranges from 54.01% ($cpl_{l\lambda_1}$) to 57.49% ($cpl_{lr}$).<sup>8</sup> Unsurprisingly, the accuracy in $cpl_{l\lambda_1}$ is lower than the other configurations since $cpl_{l\lambda_1}$ trades off training accuracy on the number of literals in the computation process.

Next most accurate are the heuristic methods *cn2* (48.03%) followed by *ripper* (44.81%). The average accuracy of *imli*<sub>16</sub> and *imli*<sub>1</sub> is 35.47% and 29.7% respectively, while the average accuracy of *twostg* is 29.6% and *ids* is 26.78. Finally, the worst accuracy is demonstrated by *sp*<sub>λ₁</sub> and *opt* (18.84% and 18.27% on average respectively) as these tools fail to provide prediction information for many datasets within the time limit. We will omit further discussion of *sp* and *opt*<sub>λ₁</sub> since they solve so few datasets.

<footnote n="8"><sup>8</sup> Note that most datasets we used represent non-binary classification. Also, DSes are not to be compared with BTs. As Figure 4b shows (and as our work aims to demonstrate), our approach outperforms the state-of-the-art DS inference methods in terms of accuracy.</footnote>
<figure><img src="image_355.png" /><figcaption>(a) Acry: vs. *cn2*.</figcaption></figure>

<figure><img src="image_356.png" /><figcaption>(c) Acry: vs. *ripper*.</figcaption></figure>

<figure><img src="image_357.png" /><figcaption>(b) Xp size: vs. *cn2*.</figcaption></figure>

<figure><img src="image_358.png" /><figcaption>(d) Xp sz: vs. *ripper*.</figcaption></figure>

Figure 5 Comparison of *cpltr* vs. *cn2* and *ripper* in terms of accuracy and explanation size.

<figure><img src="image_359.png" /><figcaption>(a) Acry: vs. *twostg*.</figcaption></figure>

<figure><img src="image_360.png" /><figcaption>(c) Acry: vs. *imli*₁₆.</figcaption></figure>

<figure><img src="image_361.png" /><figcaption>(b) Xp sz: vs. *twostg*.</figcaption></figure>

<figure><img src="image_362.png" /><figcaption>(d) Xp sz: vs. *imli*₁₆.</figcaption></figure>

Figure 6 *cpltr* vs. *imli*₁₆ and *twostg* in terms of accuracy and explanation size.
**Model Complexity.** Figure 4c illustrates the comparison among selected approaches regarding the total number of literals used in each DS solution. The average number of literals are in order: *imli*₁ (2.77), *imli*₁₆ (8.26), *ids* (21.14), *ripper* (38.47), *cpl*₁λ₁ (118.47), *cpl*ᵢₜ (157.53), *cpl*ᵢ (213.27), *twostg* (265.98), *cpl*ᵢ (584.39) *cpl* (620.82), *cn*₂ (700.49). Clearly, rule reduction and literal reduction can significantly reduce the size of the model without significantly affecting accuracy. Note how our approaches while significantly larger than the least accurate competitors, are significantly smaller than the most accurate competitor *cn*₂.

**Explanation Size.** Figure 4d shows the explanation size for each competitor. The average explanation sizes are in order: *imli*₁ (2.61), *imli*₁₆ (3.00), *cpl*₁λ₁ (12.14), *ids* (15.28), *twostg* (17.5), *cpl*ᵢₜ (25.34), *cpl*ᵢ (26.18), *ripper* (29.08), *cn*₂ (81.93), *cpl*ᵢ (234.46), *cpl* (240.88). Figure 4d demonstrates that post-hoc literal reduction not only helps decrease the number of literals required to explain DS models, but also enables DSes to remain accurate, whereas rule reduction does not contribute to smaller explanations. With literal reduction applied our approaches are very competitive in terms of explanation size.

**Detailed Comparison.** While cactus plots allow us to compare many methods over a large suite of benchmarks, they do not allow direct comparison on individual benchmarks. We provide a detailed comparison of *cpl*ᵢₜ versus other decision set inference approaches in Figures 5 and 6, including *cn*₂, *ripper*, *twostg*, and *imli*₁₆.<fnref n="9" /> The scatter plots depicting explanation size are obtained for the datasets solvable by both competitors. Note that *cpl*ᵢₜ can generate more accurate DSes than the competitors. Also observe that the explanation size of DSes computed by *cpl*ᵢₜ is smaller than *cn*₂ and comparable with *twostg*. Although the explanation size of DSes in *cpl*ᵢₜ is larger than *ripper* and *imli*₁₆, the two approaches are less interpretable as they compute DSes representing only one class.

**Summary.** The experimental results were performed on various datasets, demonstrating that our approach computes DSes that outperform the state-of-the-art competitors in terms of accuracy and yield comparable explanation size to them.

# 6 Conclusions

This paper introduced a novel *anytime* approach to generating decision sets by means of on-demand extraction of generalized abductive explanations for boosted tree models. It can be used for exhaustive compilation of a BT model wrt. the entire feature space, or target a set of training instances. Augmented by a number of post-hoc model reduction techniques, the approach is shown to compute decision sets that are more accurate than decision sets computed by the state-of-the-art algorithms and comparable with them in terms of explanation size.

As the proposed approach targets generating a decision set by compiling a BT, a natural line of future work is to extend the proposed approach to compile BTs into the other interpretable models, i.e. decision trees and decision lists, making use of AXp extraction for BTs. Additionally, another future work is to apply AXp extraction to compile other accurate black box models, e.g. neural networks, into decision sets.

<footnote n="9">⁹ The average results across the folds are given in the appendix.</footnote>
References

1 ACM. Fathers of the deep learning revolution receive ACM A.M. Turing award. http://tiny.cc/9plzpz, 2018.

2 Jeremias Berg, Emir Demirovic, and Peter J. Stuckey. Core-boosted linear search for incomplete MaxSAT. In *CPAIOR*, pages 39–56, 2019.

3 Armin Biere, Marijn Heule, Hans van Maaren, and Toby Walsh, editors. *Handbook of Satisfiability*. IOS Press, 2021.

4 Leo Breiman, J. H. Friedman, R. A. Olshen, and C. J. Stone. *Classification and Regression Trees*. Wadsworth, 1984.

5 Tianqi Chen and Carlos Guestrin. XGBoost: A scalable tree boosting system. In *KDD*, pages 785–794, 2016.

6 Peter Clark and Robin Boswell. Rule induction with CN2: some recent improvements. In *EWSL*, pages 151–163, 1991.

7 Peter Clark and Tim Niblett. The CN2 induction algorithm. *Machine Learning*, 3:261–283, 1989.

8 William W. Cohen. Fast effective rule induction. In *ICML*, pages 115–123, 1995.

9 Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL: http://archive.ics.uci.edu/ml.

10 Jiri Filip and Tomas Kliegr. Pyids-python implementation of interpretable decision sets algorithm by lakkaraju et al, 2016. In *RuleML+RR*, 2019.

11 Nicholas Frosst and Geoffrey E. Hinton. Distilling a neural network into a soft decision tree. In *CEX@AI*\*IA, volume 2071 of *CEUR Workshop Proceedings*. CEUR-WS.org, 2017.

12 Bishwamittra Ghosh and Kuldeep S. Meel. IMLI: an incremental framework for MaxSAT-based learning of interpretable classification rules. In *AIES*, pages 203–210. ACM, 2019.

13 Jianping Gou, Baosheng Yu, Stephen J. Maybank, and Dacheng Tao. Knowledge distillation: A survey. *Int. J. Comput. Vis.*, 129(6):1789–1819, 2021.

14 Gurobi Optimization. Gurobi optimizer reference manual, 2022. URL: http://www.gurobi.com/.

15 Laurent Hyafil and Ronald L. Rivest. Constructing optimal binary decision trees is NP-complete. *Inf. Process. Lett.*, 5(1):15–17, 1976.

16 Alexey Ignatiev. Towards trustable explainable AI. In *IJCAI*, pages 5154–5158, 2020. doi:10.24963/ijcai.2020/726.

17 Alexey Ignatiev, Yacine Izza, Peter J. Stuckey, and João Marques-Silva. Using MaxSAT for efficient explanations of tree ensembles. In *AAAI*, pages 3776–3785, 2022.

18 Alexey Ignatiev, Edward Lam, Peter J. Stuckey, and Joao Marques-Silva. A scalable two stage approach to computing optimal decision sets. In *AAAI*, pages 3806–3814, 2021.

19 Alexey Ignatiev, Antonio Morgado, and Joao Marques-Silva. RC2: an efficient MaxSAT solver. *J. Satisf. Boolean Model. Comput.*, 11(1):53–64, 2019.

20 Alexey Ignatiev, Nina Narodytska, and Joao Marques-Silva. Abduction-based explanations for machine learning models. In *AAAI*, pages 1511–1519, 2019. doi:10.1609/aaai.v33i01.33011511.

21 Alexey Ignatiev, Nina Narodytska, and Joao Marques-Silva. On validating, repairing and refining heuristic ML explanations. *CoRR*, abs/1907.02509, 2019. arXiv:1907.02509.

22 Alexey Ignatiev, Filipe Pereira, Nina Narodytska, and João Marques-Silva. A SAT-based approach to learn explainable decision sets. In *IJCAR*, pages 627–645, 2018.

23 Yacine Izza, Alexey Ignatiev, and Joao Marques-Silva. On tackling explanation redundancy in decision trees. *J. Artif. Intell. Res.*, 75:261–321, 2022.

24 Yacine Izza, Alexey Ignatiev, Peter J. Stuckey, and Joao Marques-Silva. Delivering inflated explanations. *CoRR*, abs/2306.15272, 2023.

25 Yacine Izza and Joao Marques-Silva. On explaining random forests with SAT. In *IJCAI*, July 2021.
38:18 From Formal Boosted Tree Explanations to Interpretable Rule Sets

26 Anil P. Kamath, Narendra Karmarkar, K. G. Ramakrishnan, and Mauricio G. C. Resende. A continuous approach to inductive inference. *Math. Program.*, 57:215–238, 1992.

27 Kentaro Kanamori, Takuya Takagi, Ken Kobayashi, Yuichi Ike, Kento Uemura, and Hiroki Arimura. Ordered counterfactual explanation by mixed-integer linear optimization. In *AAAI*, pages 11564–11574. AAAI Press, 2021.

28 Richard M. Karp. Reducibility among combinatorial problems. In *Complexity of Computer Computations*, pages 85–103, 1972.

29 Himabindu Lakkaraju, Stephen H. Bach, and Jure Leskovec. Interpretable decision sets: A joint framework for description and prediction. In *KDD*, pages 1675–1684, 2016.

30 Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. *Nature*, 521(7553):436, 2015.

31 Zhendong Lei and Shaowei Cai. Solving (weighted) partial maxsat by dynamic local search for SAT. In *IJCAI*, pages 1346–1352, 2018.

32 Zhendong Lei and Shaowei Cai. Nudist: An efficient local search algorithm for (weighted) partial maxsat. *Comput. J.*, 63(9):1321–1337, 2020.

33 Zachary C. Lipton. The mythos of model interpretability. *Commun. ACM*, 61(10):36–43, 2018.
doi:10.1145/3233231.

34 Scott M. Lundberg and Su-In Lee. A unified approach to interpreting model predictions. In *NeurIPS*, pages 4765–4774, 2017. URL: https://proceedings.neurips.cc/paper/2017/hash/8a20a8621978632d76c43dfd28b67767-Abstract.html.

35 Chuan Luo, Shaowei Cai, Kaile Su, and Wenxuan Huang. CCEHC: an efficient local search algorithm for weighted partial maximum satisfiability. *Artif. Intell.*, 243:26–44, 2017.

36 Dmitry Malioutov and Kuldeep S. Meel. MLIC: A MaxSAT-based framework for learning interpretable classification rules. In *CP*, pages 312–327, 2018.

37 Joao Marques-Silva and Alexey Ignatiev. Delivering trustworthy AI through formal XAI. In *AAAI*, pages 12342–12350, 2022.

38 Hayden McTavish, Chudi Zhong, Reto Achermann, Ilias Karimalis, Jacques Chen, Cynthia Rudin, and Margo I. Seltzer. Fast sparse decision tree optimization via reference ensembles. In *AAAI*, pages 9604–9613, 2022.

39 Ryszard S Michalski. On the quasi-minimal solution of the general covering problem. In *International Symposium on Information Processing*, pages 125–128, 1969.

40 Nina Narodytska, Alexey Ignatiev, Filipe Pereira, and João Marques-Silva. Learning optimal decision trees with SAT. In *IJCAI*, pages 1362–1368, 2018.

41 Randal S. Olson, William G. La Cava, Patryk Orzechowski, Ryan J. Urbanowicz, and Jason H. Moore. PMLB: a large benchmark suite for machine learning evaluation and comparison. *BioData Min.*, 10(1):36:1–36:13, 2017.

42 Axel Parmentier and Thibaut Vidal. Optimal counterfactual explanations in tree ensembles. In *ICML*, volume 139 of *PMLR*, pages 8422–8431, 2021.

43 Marco Túlio Ribeiro, Sameer Singh, and Carlos Guestrin. "why should I trust you?": Explaining the predictions of any classifier. In *KDD*, pages 1135–1144, 2016. doi:10.1145/2939672.2939778.

44 Marco Túlio Ribeiro, Sameer Singh, and Carlos Guestrin. Anchors: High-precision model-agnostic explanations. In *AAAI*, pages 1527–1535, 2018. URL: https://www.aaai.org/ocs/index.php/AAAI/AAAI18/paper/view/16982.

45 Ronald L. Rivest. Learning decision lists. *Mach. Learn.*, 2(3):229–246, 1987.

46 Cynthia Rudin and Seyda Ertekin. Learning customized and optimized lists of rules with mathematical programming. *Mathematical Programming Computation*, 10:659–702, 2018.

47 Paul Saikko, Jeremias Berg, and Matti Järvisalo. LMHS: A SAT-IP hybrid MaxSAT solver. In *SAT*, pages 539–546, 2016.

48 Andy Shih, Arthur Choi, and Adnan Darwiche. A symbolic approach to explaining bayesian network classifiers. In *IJCAI*, pages 5103–5111, 2018. doi:10.24963/ijcai.2018/708.
49 Dylan Slack, Sophie Hilgard, Emily Jia, Sameer Singh, and Himabindu Lakkaraju. Fooling LIME and SHAP: adversarial attacks on post hoc explanation methods. In AIES, pages 180–186, 2020. doi:10.1145/3375627.3375830.

50 Jinqiang Yu, Alexey Ignatiev, Peter J. Stuckey, and Pierre Le Bodic. Computing optimal decision sets with SAT. In CP, pages 952–970, 2020.

51 Jinqiang Yu, Alexey Ignatiev, Peter J Stuckey, and Pierre Le Bodic. Learning optimal decision sets and lists with sat. JAIR, 72:1251–1279, 2021.

52 Jinqiang Yu, Alexey Ignatiev, Peter J. Stuckey, Nina Narodytska, and Joao Marques-Silva. Eliminating the impossible, whatever remains must be true: On extracting and applying background knowledge in the context of formal explanations. In AAAI, 2023.

## A Summaries of Results Across Folds

<figure><img src="image_363.png" /><figcaption>Figure 7 Experimental results of runtime and accuracy across folds.</figcaption></figure>

<figure><img src="image_364.png" /><figcaption>Figure 8 Experimental results of model complexity and explanation size across folds.</figcaption></figure>

Figures 7 and 8 illustrate the average experimental results across folds regarding scalability, accuracy, model complexity, and explanation size. Since 5-fold cross validation is used, these results for each dataset are obtained from the average of 5 pairs of training and test data. Here, observations similar to those described in Section 5 can be made, i.e. the best
38:20 From Formal Boosted Tree Explanations to Interpretable Rule Sets

scalability and accuracy among selected DS competitors are both demonstrated by *cpl* and
*cpl*<sub>*</sub>, * ∈ {l, r, lr, lλ<sub>1</sub>} , while *imli*<sub>1</sub> and *imli*<sub>16</sub> show the smallest model complexity and
explanation size.

<figure><img src="image_365.png" /><figcaption>Figure 9 *cpl<sub>lr</sub>* vs. CN2 and RIPPER across folds in terms of accuracy and explanation size.</figcaption></figure>
<figure><img src="image_366.png" /><figcaption>(a) Acry: vs. *twostg*.</figcaption></figure>

<figure><img src="image_367.png" /><figcaption>(c) Acry: vs. *imli*₁₆.</figcaption></figure>

<figure><img src="image_368.png" /><figcaption>(b) Xp sz: vs. *twostg*.</figcaption></figure>

<figure><img src="image_369.png" /><figcaption>(d) Xp sz: vs. *imli*₁₆.</figcaption></figure>

<caption>**Figure 10** *cpl*<sub>lr</sub> vs. *imli*<sub>16</sub> and *twostg* Across Folds in terms of accuracy and explanation size.</caption>

## B Detailed Comparisons Across Folds

In this appendix, we provide a detailed comparison of *cpl*<sub>lr</sub> versus other decision set inference approaches across folds.

Figure 9 and Figure 10 detail the comparisons of *cpl*<sub>lr</sub> with CN2, RIPPER, *imli*<sub>16</sub> and *twostg* in terms of average accuracy and explanation size across folds. As can be seen in Figure 9a, the accuracy of DSes generated by *cpl*<sub>lr</sub> is higher than the accuracy of CN2, where the average accuracy is 57.49% and 48.03%, respectively. Additionally, Figure 9b demonstrates that the explanation size of DSes produced by CN2 (81.93 on average) can be two orders of magnitude larger than the explanation size of *cpl*<sub>lr</sub> (25.88 on average).

Figure 9c illustrates that the average accuracy in RIPPER is 44.81%, which is 12.68% lower than the accuracy in *cpl*<sub>lr</sub>. Although Figure 9d depicts that RIPPER is comparable with *cpl*<sub>lr</sub> regarding explanation size (29.08 and 25.34 on average respectively), RIPPER is less interpretable as it computes DSes representing only one class.

As can be observed in Figure 10a, the accuracy of *twostg* (29.67% on average) is 27.82% lower than the accuracy in *cpl*<sub>lr</sub> while Figure 10b illustrates that the explanation size is comparable between the two approaches. Finally, Figure 10c demonstrates that the accuracy of *imli*<sub>16</sub> is 22.02% lower than the accuracy of *cpl*<sub>lr</sub> on average. However, as can be seen in Figure 10d, the explanation size of *imli*<sub>16</sub> is smaller than the explanation size of *cpl*<sub>lr</sub> but *imli*<sub>16</sub> generates DSes targeting only a single class, which significantly diminishes the interpretability of computed DSes.
<figure><img src="image_370.png" /></figure>
# Searching for Smallest Universal Graphs and Tournaments with SAT

Tianwei Zhang

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

Algorithms and Complexity Group, TU Wien, Austria

Stefan Szeider

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

Algorithms and Complexity Group, TU Wien, Austria

## Abstract

A graph is induced *k*-universal if it contains all graphs of order *k* as an induced subgraph. For over half a century, the question of determining smallest *k*-universal graphs has been studied. A related question asks for a smallest *k*-universal tournament containing all tournaments of order *k*.

This paper proposes and compares SAT-based methods for answering these questions exactly for small values of *k*. Our methods scale to values for which a generate-and-test approach isn't feasible; for instance, we show that an induced 7-universal graph has more than 16 vertices, whereas the number of all connected graphs on 16 vertices, modulo isomorphism, is a number with 23 decimal digits. Our methods include static and dynamic symmetry breaking and lazy encodings, employing external subgraph isomorphism testing.

**2012 ACM Subject Classification** Mathematics of computing → Extremal graph theory; Software and its engineering → Constraint and logic languages; Hardware → Theorem proving and SAT solving; Mathematics of computing → Graph enumeration

**Keywords and phrases** Constrained-based combinatorics, synthesis problems, symmetry breaking, SAT solving, subgraph isomorphism, tournament, directed graphs

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

**Supplementary Material** Software: https://doi.org/10.5281/zenodo.8147732

**Funding** The project leading to this publication has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 101034440, and was supported by the Vienna Science and Technology Fund (WWTF) within the project ICT19-065.

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

**Acknowledgements** This work was carried out in part while the second author visited the Simons Institute for the Theory of Computing, University of Berkeley, within the program *Extended Reunion: Satisfiability*. The authors thank Ciaran McCreesh for advising them to use and modify the Glasgow subgraph solver.

<figure><img src="image_374.png" /><figcaption>Figure 1 A smallest induced 4-universal graph and how all 11 graphs of order 4 embeds into it.</figcaption></figure>

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

© Tianwei Zhang and Stefan Szeider;
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. 39; pp. 39:1–39:20

Leibniz International Proceedings in Informatics

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

LIPICS Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
<figure><img src="image_377.png" /><figcaption>Figure 2 A smallest 4-universal tournament and how all 4 tournaments of order 4 embeds into it.</figcaption></figure>

1 Introduction

A graph $G$ is *induced k-universal* if it contains all *k*-vertex graphs as induced subgraphs; in other words, any *k*-vertex graph can be obtained from $G$ by deleting some of $G$'s vertices. This notion of a universal graph, which extends naturally to directed graphs and tournaments (orientations of complete graphs) was introduced by Rado in 1964 [25] and has since then attracted much attention in combinatorics. Figures 1 and 2 provide examples of a universal graph and a universal tournament, respectively. The fundamental question is to determine for each integer $k \ge 1$ the smallest $n$ such that an induced $k$-universal graph with $n$ vertices exists; this number is denoted $f(k) = n$. In 1965, Moon [21] obtained the bounds $2^{(k-1)/2} \le f(k) \le O(k \cdot 2^{k/2})$; Alstrup et al. [2] recently improved the upper bound to $f(k) \le 16 \cdot 2^{k/2}$. The tight asymptotic upper bound of $f(k) \le (1+o(1))2^{(k-1)/2}$ is due to Alon [1].

In contrast to these general bounds and asymptotic results, very little is known about exact $f(k)$ values, even for small $k$. Trimble [26] carried out a brute-force search, expanding an approach suggested by Preen [24], and could determine the values of $f(k)$ for $k \le 6$ and the interval $16 \le f(7) \le 18$. The brute-force search enumerates all candidate $n$-vertex graphs up to isomorphism and tests for each of them whether it contains all the $k$-vertex graphs as induced subgraphs. We refer to the $k$-vertex graphs that need to be checked as induced subgraphs as *pattern graphs*. Tools exist for the isomorph-free enumeration of all connected $n$-vertex graphs (we can assume that a candidate graph is connected), but the number of such graphs exceeds eleven million for $n = 10$, and one billion for $n = 11$ [23, A001349]. Thus, the feasibility of the brute-force approach quickly hits a rigid boundary. In particular, we cannot tighten the gap $16 \le f(7) \le 18$ in this way, as enumerating all graphs with 16 vertices is far out of reach, although we can easily enumerate all the 7-vertex pattern graphs; the exact number is 1044 [23, A000088].

This paper proposes new methods that allow us to break this boundary. The idea is to formulate the problem as a *synthesis problem* for the universal graph rather than the easier problem of testing a candidate graph for being universal. This way, we can avoid the bottleneck of enumerating all candidate graphs. We formulate the synthesis problem in propositional logic to harness the power of solvers for the propositional satisfiability problem (SAT) [8]. Suppose the sought-for universal graph $G$ has the numbers $1, \dots, n$ as its vertices. For each pair $1 \le i < j \le n$, we introduce a propositional variable $e_{i,j}$ (an *edge variable*) whose truth value determines whether there is an edge between vertices $i$ and $j$ in $G$.

We propose and compare various strategies and techniques of encoding the property of $G$ being induced $k$-universal. Each of the encodings must combine two fundamental properties: (i) to ensure that the found graph is indeed induced $k$-universal and (ii) to break symmetries to minimize the enormous search space. We have two fundamentally different approaches for properties (i) and (ii).

Our first approach for ensuring that the found graph is indeed induced *k*-universal uses a *direct SAT encoding*. For the relevant cases, the number of pattern graphs is reasonably small, and we can enumerate them up to isomorphism. For each pattern graph, we obtain
a formula which constrains the edge variables in such a way that the sought-for universal graph contains an induced subgraph that is isomorphic to the pattern graph. A conjunction over the formulas, one for each pattern graph, yields the SAT encoding.

Our second approach for ensuring universality utilizies the *Glasgow subgraph solver (GSS)* [19], a powerful constraint-based tool for testing subgraphs. However, during the SAT solver's search we don't have a candidate graph available for testing whether it accommodates all the required induced subgraphs. A partial truth assignment instead sets a subset of the edge variables to true or false, with some edge variables remaining undecided. Thus, a partial truth assignment gives rise to a *partially defined graph*. Our objective is to find out already for a partially defined graph whether it could be extended to a fully defined induced *k*-universal graph. Indeed, we could achieve this by accessing GSS's internal functionality.

Our first approach to breaking symmetries utilizes the *SAT modulo Symmetries (SMS)* framework [17]. Also SMS works on partially defined graphs as represented by a partial truth assignment to edge variables. Whenever the solver determines an edge, the partially defined graph represented by the current partial truth assignment is sent to an external propagator, which determines whether the partially defined graph can be extended to a fully defined graph that is canonical (i.e., a unique graph within its ismorphism class). If not, a clause is sent back to the solver.

Our second approach to symmetry breaking is based on the concept of *templates* which serve two purposes: symmetry breaking and search space partition. In an offline phase, we fix a small collection of highly symmetric pattern graphs and compute all possible ways these pattern graphs can interact as induced subgraphs in the universal graph, up to isomorphism. Each possible interaction gives rise to a template. By hardcoding a template, we get a formula that is satisfiable if and only if there exists an induced *k*-universal graph on *n* vertices in which the fixed collection of pattern graphs interact as prescribed in the template. *f(k) ≤ n* if an only if for at least one of the templates, the corresponding formula is satisfiable. In addition to the symmetry breaking aspect, the templates approach allows us to parallelize the search, running SAT calls for various templates independently.

We also evaluate our SAT-based approach to determining smallest *k*-universal tournaments. A tournament is a directed graph that can be obtained from a complete undirected graph by orienting its edges. The name tournament originates from such a directed graph's interpretation as the outcome of a round-robin tournament in which every player encounters every other player exactly once, and in which no draws occur. Tournaments are a central combinatorial object whose properties have been intensively studied. The notion of universality applies to tournaments in a natural way: we say that a tournament is *k*-universal if it contains each tournament on *k*-vertices as a subgraph. Let *t(k)* denote the smallest number of vertices a *k*-universal tournament can have. The study of *k*-universal tournaments goes back more than 50 years to Moon's book (entirely dedicated to tournaments [22]) where he determined the bounds $2^{(k-1)/2} \le t(k) \le O(k \cdot 2^{k/2})$. The improvements on the bounds parallel the ones for induced *k*-universal graphs: $f(k) \le 16 \cdot 2^{\lceil k/2 \rceil}$ by Alstrup et al. [2] and $f(k) \le (1+o(1))2^{(k-1)/2}$ by Alon [1].

The situation regarding exact values for $t(k)$ is even worse than for $f(k)$. While the number of tournaments up to isomorphism is for all $k \le 19$ [23, A000568], very little can be found on the value of $t(k)$ in the literature. The only obvious values are $t(1) = 1$ and $t(2) = 2$. In Moon's book [22], it is left as an exercise to the reader to determine the value of $t(3)$ and $t(4)$, indicating that these numbers might be easily deduced.

We implemented our SAT-based approaches to obtain new results on smallest induced *k*-universal graphs and smallest *k*-universal tournaments. In particular, we could improve the known lower bound for the order of an induced 7-universal graph by showing that no induced
39:4 Searching for Smallest Universal Graphs and Tournaments with SAT

<table><caption>**Table 1** Overview of results along the order *k* of the pattern graphs or tournaments. *G(k)/T(k)* give the number of pattern graphs/tournaments of order *k*, *f(k)/t(k)* give the order of a smallest (induced) *k*-universal graph or tournament, and *F(k)/T(k)* their number modulo isomorphism, respectively. New results are indicated in bold face.</caption><thead><tr><th rowspan="2">k</th><th colspan="3">Graphs</th><th colspan="3">Tournaments</th></tr><tr><th>G(k)</th><th>f(k)</th><th>F(k)</th><th>T(k)</th><th>t(k)</th><th>T(k)</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>2</td><td>2</td><td>3</td><td>2</td><td>1</td><td>2</td><td>1</td></tr><tr><td>3</td><td>4</td><td>5</td><td>5</td><td>2</td><td>4</td><td><b>3</b></td></tr><tr><td>4</td><td>11</td><td>8</td><td>438</td><td>4</td><td>5</td><td><b>1</b></td></tr><tr><td>5</td><td>34</td><td>10</td><td>22</td><td>12</td><td><b>8</b></td><td><b>1643</b></td></tr><tr><td>6</td><td>156</td><td>14</td><td>&gt;<b>36294</b></td><td>56</td><td><b>10</b></td><td><b>1088</b></td></tr><tr><td>7</td><td>1044</td><td>[17,18]</td><td>-</td><td>456</td><td>[13,15]</td><td>-</td></tr></tbody></table>

7-universal graph with 16 vertices exists, leaving the possibilities 17 and 18 for the smallest solution, and we could determine the exact order of smallest *k*-universal tournaments up to *k* = 6 and leaving the possibilities 13,14, and 15 for the smallest solution for *k* = 7. We also determined the precise number of optimal solutions for *k*-universal tournaments for *k* ≤ 6, up to isomorphism, which is of independent interest for combinatorial research. Our main results are summarized in Table 1, smallest *k*-universal tournaments that we identified are shown in Figures 6 and 7 in the appendix.

## 2 Preliminaries

We denote the set {1, ..., n} by [n] and the set {m, m + 1, ..., n} by [m, n] for m ≤ n.

### 2.1 Graphs

We consider simple undirected graphs *G*, denoting the vertex set and the edge set of *G* by *V*(*G*) and *E*(*G*), respectively. The order or a graph is the number |*V*(*G*)| of its vertices. An edge between vertices *u*, *v* ∈ *V*(*G*) is denoted *uv* or equivalently *vu*. A graph *G'* is a *subgraph* of a graph *G* if *V*(*G'* ) ⊆ *V*(*G*) and *E*(*G'* ) ⊆ *E*(*G*). A subgraph *G'* of *G* is *induced* if for any *u*, *v* ∈ *V*(*G'* ), *uv* ∈ *E*(*G*) implies *uv* ∈ *E*(*G'* ). Thus, an induced subgraph is determined by its set of vertices.

The neighborhood of a vertex *u* ∈ *V*(*G*) is denoted *N*<sub>*G*</sub>(*u*), i.e., *N*<sub>*G*</sub>(*u*) := { *x* ∈ *G* | *ux* ∈ *E*(*G*) }. Two vertices *v*, *v'* ∈ *G* are twins if they have the exact same neighbors excluding each other, i.e., *N*<sub>*G*</sub>(*u*) \ { *v* } = *N*<sub>*G*</sub>(*v*) \ { *u* }.

The complement graph $\bar{G}$ of a graph $G$ has $V(\bar{G}) = V(G)$ and $E(\bar{H}) := \{ uv \mid u \neq v \in V(H), uv \notin E(H) \}$.

We denote the complete graph of order *k* by $K_k$, the complete bipartite graph by $K_{l,r}$.

A monomorphism from $H$ to $G$ is an injective mapping $\phi: V(H) \to V(G)$ such that for all $u, v \in V(H)$ we have $\phi(u)\phi(v) \in E(G)$ if $uv \in E(H)$. Such a monomorphism $\phi$ is faithful if for all $u, v \in V(H)$ we have $uv \in E(H)$ if $\phi(u)\phi(v) \in E(G)$ [11]. We say $H$ embeds into $G$ if there is a faithful monomorphism from $H$ to $G$. A bijective faithful monomorphism is an isomorphism. We refer to $H$ as the pattern graph and $G$ as the target graph. For a monomorphims $\phi$ from $H$ to $G$ we put $\phi(V(H)) := \{ \phi(v) \mid v \in V(H) \}$.
Let $\mathcal{G}(k)$ denote the class of all graphs of order $k$ up to isomorphism. A graph $G$ is
induced $k$-universal if all $H \in \mathcal{G}(k)$ embed into $G$. $f(k)$ is the smallest order of an induced
$k$-universal graph, and $F(k)$ is the number of non-isomorphic $k$-universal graphs of order
$f(k)$.

## 2.2 Partially Defined Graphs

The notion of partially defined graphs was introduced in the context of SMS [17], for
capturing the combinatorial object represented by a partial truth assignment on the edge
variables. In this paper, we utilize this concept for two further purposes: to check with
an external subgraph solver whether the current branch of the search has the potential of
success (Section 4.5) and for the formulation of our template method for symmetry breaking
and search partition (Section 4.6).

A *partially defined graph* is a graph where some edges are undefined in the sense that their presence in the graph is open. Formally, a partially defined graph *G* is a graph whose edge set *E*(*G*) is split into disjoint sets *E*<sub>*d*</sub>(*G*) and *E*<sub>*u*</sub>(*G*), the sets of defined and undefined edges, respectively. *G* is *fully defined* if *E*<sub>*u*</sub>(*G*) = ∅. For a partially defined graph *G*, *X*(*G*) denotes the set of all fully defined graphs *G* can be extended to, i.e., *G*′ ∈ *X*(*G*) if and only if *V*(*G*)' = *V*(*G*), *E*<sub>*d*</sub>(*G*) ⊆ *E*(*G*)' ⊆ *E*(*G*).

We extend the notion of a faithful monomorphism and an isomorphism to partially defined
graphs. Let $H, G$ be partially defined graphs. An injective mapping $\phi : V(H) \to V(G)$ is a
faithful monomorphism from $H$ to $G$ if for all $u, v \in V(H)$ it holds that

1. $uv \in E_d(H)$ if and only if $\phi(u)\phi(v) \in E_d(G)$, and

2. $uv \in E_u(H)$ if and only if $\phi(u)\phi(v) \in E_u(G)$.

If there is a faithful monomorphism from $H$ to $G$ we say that $H$ embeds into $G$. We say that $H$ and $G$ are *isomorphic* if there is a *bijective* faithful monomorphism from $H$ to $G$.

## 2.3 Directed Graphs and Tournaments

We denote a directed graph (or digraph) $D$ by its vertex set $V(D)$ and arc set $A(D) \subseteq V(D) \times V(D)$. A digraph $D$ is an *oriented graph* if it has no directed digon $\{(u, v), (v, u)\}$. (Induced) subdigraphs, isomorphisms, and faithful monomorphisms are defined for digraphs analogously to graphs. A *tournament* is an oriented graph with a maximal number of arcs, i.e., adding any further arc creates a digon. A tournament is *transitive* if it contains no directed cycles. A tournament $D'$ is a subtournament of tournament $D$ if it is a subdigraph of $D$. If a tournament $H$ is isomorphic to a subtournament of $D$, then we say that $H$ embeds into $D$. Hence $H$ embeds into $D$ if and only if there exists a faithful monomorphism $\phi : V(H) \to V(D)$.

Let $\mathcal{T}(k)$ denote the class of all tournaments of order $k$ up to isomorphism. A tournament $D$ is $k$-universal if all $H \in \mathcal{T}(k)$ embed into $D$. Let $\mathcal{T}(k)$ denote the number of non-isomorphic tournaments of order $t(k)$ that are induced $k$-universal tournaments.

## 2.4 SAT

We consider propositional formals in conjunctive normal form (CNF). A CNF formula is
a conjunction of *clauses*, each clause is a disjunction of *literals*, a literal is a propositional
variable *x* or its negation **¬**(*x*). A propositional formula *F* is *satisfiable* if there exists a
mapping **τ** that assigns each variable a truth value ∈ {0, 1} such that each clause contains a
literal *x* with **τ」(*x*)** = 1 or a literal **¬**(*x*) with **τ」(*x*)** = 0. A truth assignment is *partial* for a CNF
formula if it is defined only for a subset of the formula's variables. SAT solvers are tools that decide whether a given CNF formula is satisfiable or not [8]. Today's most powerful complete SAT solvers follow the conflict driven clause learning (CDCL) paradigm. Modern CDCL SAT solvers like Cadical [4] can certify the unsatisfiability of a CNF formula by producing DRAT proofs (deletion, reverse asymmetric tautology) that can be independently checked and verified [12, 27].

### 3 A Lower Bound for Induced $k$ Universal Graphs

Trimble [26] showed that there is no graph of order less than $2k + 2$ that is universal for $\{K_k, \overline{K_k}, K_{3,3}, \overline{K_{3,3}}\}$, for all $k \ge 6$. We extend this argument to show the following more general proposition, which gives a strictly tighter lower bound for $k \ge 8$.

**Proposition 1.** $2k + 2\lfloor\frac{k}{2}\rfloor - 4 \le f(k)$ for all $k \ge 2$.

**Proof.** Suppose $G$ is an induced $k$-universal graph. Graphs $K_k$ and $\overline{K_k}$ are both elements of $G(k)$, and therefore there exist faithful monomorphisms $\phi_K : V(K_k) \to V(G)$ and $\phi_I : V(\overline{K_k}) \to V(G)$. Let $S_1 := \{\phi_K(v) \mid v \in V(K_k)\}$, $S_2 := \{\phi_I(v) \mid v \in V(\overline{K_k})\}$ and $S_3 := V(G)/(S_1 \cup S_2)$. Clearly, the subgraphs induced by $S_1$ and $S_2$ are a clique and an independent set, respectively, and therefore may overlap by no more than one vertex. This implies that $|S_1 \cup S_2| \ge 2k - 1$ and $|S_3| \le n - 2k + 1$. Define $H := K_{\lfloor\frac{k}{2}\rfloor, \lfloor\frac{k}{2}\rfloor+(k \bmod 2)}$. Let $G_1$ be an induced subgraph of $G$ that is isomorphic to $H$. Since there are no edges between the two cliques of $G_1$, it must be the case that at least one of the two cliques of $G_1$ does not intersect $S_1$. Denote this clique by $C$. Since $S_2$ is an independent set in $G$, $C$ can only intersect $S_2$ by at most one vertex. Since $C$ has at least $\lfloor\frac{k}{2}\rfloor$ vertices, it follows that $C$ intersects $S_3$ by at least $\lfloor\frac{k}{2}\rfloor - 1$ vertices. In other words, the subgraph induced by $S_3$ contains a clique $D$ of size $\lfloor\frac{k}{2}\rfloor - 1$.

Now consider the graph $\overline{H}$. Since $\overline{H}$ is an induced subgraph of $G$, it follows that $H = \overline{H}$ is an induced subgraph of $\overline{G}$. We can repeat the argument of the previous paragraph with the roles of $S_1$ and $S_2$ reversed to show that the subgraph induced by $S_3$ contains an independent set of size $\lfloor\frac{k}{2}\rfloor - 1$. Since this independent set can overlap with the clique $D$ in at most one vertex, it follows that the minimal size of $S_3$ is $2(\lfloor\frac{k}{2}\rfloor - 1) - 1$. In other words, for any $k$-universal graph of order $n$ where $S_1 \cap S_2 \ne \emptyset$, $n - 2k + 1 \ge 2(\lfloor\frac{k}{2}\rfloor - 1) - 1$, which implies that $f(k) \ge 2k + 2\lfloor\frac{k}{2}\rfloor - 4$. Hence, the proposition is shown. ◀

**Proposition 1** gives us $f(2) \ge 2$, $f(3) \ge 4$, $f(4) \ge 8$, $f(5) \ge 10$, $f(6) \ge 14$ and $f(7) \ge 16$. Note that for $k = 4, 5, 6$, the lower bound given by the proposition is also the exact value of $f(k)$. Even though it does not give a tighter lower bound for $f(7)$, it does improve on the existing knowledge for $k \ge 8$.

**Corollary 2.** Let $G$ be a $k$-universal graph. If there exist faithful monomorphisms $\phi_K : V(K_k) \to V(G)$ and $\phi_I : V(\overline{K_k}) \to V(G)$ with $\phi_K(V(K_k)) \cap \phi_I(V(\overline{K_k})) = \emptyset$, then the order of $G$ is at least $2k + 2\lfloor\frac{k}{2}\rfloor - 3$.

### 4 Encodings for Induced $k$ universal Graphs

Throughout this section, we fix an integer $n$ and consider a potential universal graph $G$ with $V(G) = [n]$ whose edges are represented by edge variables $e_{u,v}$, $1 \le u < v \le n$.
## 4.1 Encoding Universality

First, we define a CNF formula $M(H, n)$ that ensures that $(m_{u,v})_{u \in V(H), v \in V(G)}$ encodes an injective mapping from $V(H)$ to $V(G)$.

$$M(H, n) := \bigwedge_{u \in V(H)} \left( \bigvee_{x \in [n]} m_{u,x} \right) \wedge \bigwedge_{\substack{u \neq v \in V(H) \\ x \in [n]}} \left( \neg m_{u,x} \vee \neg m_{v,x} \right) \wedge \bigwedge_{\substack{u \in V(H) \\ x \neq y \in [n]}} \left( \neg m_{u,x} \vee \neg m_{u,y} \right)$$

Now, we specify a CNF formula $F(H, n)$ which is satisfiable if and only if $H$ embeds into $G$. $F(H, n)$ encodes the existence of a faithful monomorphism $\phi$ from $H$ to $G$ in terms of variables $m_{u,v}$, $u \in V(H)$, $v \in V(G)$, which are true if and only if $\phi(u) = v$.

$$F(H, n) := M(H, n) \wedge \bigwedge_{\substack{uv \in E(H) \\ 1 \le x < y \le n}} (\neg m_{u,x} \vee \neg m_{u,y} \vee e_{x,y}) \wedge \bigwedge_{\substack{u \neq v \in V(H) \\ \text{s.t. } uv \notin E(H) \\ 1 \le x < y \le n}} (\neg m_{u,x} \vee \neg m_{v,y} \vee \neg e_{x,y}).$$

The second and third conjunct ensure that the mapping preserves edges and non-edges, respectively. The formula

$$U(k, n) = \bigwedge_{H \in \mathcal{G}(k)} F(H, n)$$

is satisfiable if and only if there exists a $k$-universal graph on $n$ vertices. From a satisfying assignment we can read off a $k$-universal graph.

## 4.2 Symmetry-breaking Based on Twins

A basic symmetry-breaking can already be achieved by observing the symmetries of a pattern graph $H \in \mathcal{G}(k)$. Specifically, if $v, v' \in V(H)$ are twins and $\phi$ is a faithful monomorphism from $H$ to $G$, then $\phi' : V(H) \to V(G)$ defined as

$$\phi'(x) := \begin{cases} \phi(v') & \text{if } x = v, \\ \phi(v) & \text{if } x = v', \\ \phi(x) & \text{otherwise,} \end{cases}$$

is also a faithful monomorphism. This means that we can always require that the faithful monomorphism encoded by $(m_{u,v})_{u \in V(H), v \in V(G)}$ preserves the order of vertices (seen as integers) for twins. The following CNF formula encodes this additional requirement.

$$\text{Twin}(k) := \bigwedge_{\substack{H \in \mathcal{G}(k) \\ u < v \in V(H) \\ \text{s.t. } u, v \text{ are twins} \\ 1 \le x < y \le n}} \neg m_{v,x} \vee \neg m_{u,y}$$

## 4.3 Embedding $K_k$ and $\overline{K_k}$

From preliminary experiments we observed that fixing (or “hardcoding”) some edges/non-edges of the target graph can drastically speed up the solving process. Trimble [26] hardcodes $K_k$ and $\overline{K_k}$ in his experiments, and we agree that it is a good place to start. Since the two subgraphs can have at most one overlapping vertex, we have the following two formulas, each hardcoding one of the two possibilities. We assume $n \ge 2k$.

$$\text{KI}_{\text{disjoint}}(k) := \bigwedge_{1 \le x < y \le k} e_{x,y} \wedge \neg e_{x+k,y+k}, \quad \text{KI}_{\text{overlap}}(k) := \bigwedge_{1 \le x < y \le k} e_{x,y} \wedge \neg e_{x+k-1,y+k-1}.$$
39:8 Searching for Smallest Universal Graphs and Tournaments with SAT

In Section 4.6, we will extend this idea by hardcoding more pattern graphs.

To sum up, the following two encodings are shared by all methods.

$$U_d(k, n) := U(k, n) \wedge \text{Twin}(k) \wedge \text{KI}_{\text{disjoint}}(k),$$

$$U_o(k, n) := U(k, n) \wedge \text{Twin}(k) \wedge \text{KI}_{\text{overlap}}(k).$$

Note that $U(k, n)$ is satisfiable if and only if $U_d(k, n)$ or $U_o(k, n)$ is satisfiable.

## 4.4 SAT Modulo Symmetries (SMS)

If a SAT solver checks the satisfiability of $U_d(k, n)$ or $U_o(k, n)$, it does not break symmetries and considers isomorphic copies of graphs implicitly represented by the edge variables, making the approach impractical. Fortunately, the SAT modulo Symmetries (SMS) framework [17] offers a solution.

For two (fully defined) graphs $G_1, G_2$, let $G_1 \preceq G_2$ if and only the adjacency matrix of $G_1$ is lexicographically smaller or equal to the adjacency matrix of $G_2$ (we consider the matrix as the string obtained by concatenating the rows of the matrix).

SMS works as follows. Whenever the SAT solver decides on an edge variable $e_{u,v}$, the partially defined graph $G$ represented by the current truth assignment on the edge variables is sent to an external propagator. With a certain pre-determined probability the propagator checks a necessary condition for $\mathcal{X}(G)$ containing a $\preceq$-minimal graph. If the condition is not met, then this branch of the search can be terminated: a clause that excludes $G$ is learned and sent back to the SAT solver.

In our context, while the satisfiability of $U(k, n)$ and $\text{Twin}(k)$ are unaffected by our choice of the $\preceq$-minimal graph as the representative target graph for each isomorphic class, $\text{KI}_{\text{disjoint}}(k)$ and $\text{KI}_{\text{overlap}}(k)$ might cause an otherwise satisfying assignment to fail the minimality check.

The solution to this is to take advantage of the following functionality of SMS. When SMS checks the $\preceq$-minimality of a partially defined graph $G$, one can specify a partition $P$ of the vertex set such that $\preceq$-minimality is only checked among all partially defined graphs that can be obtained from $G$ by permuting vertices within an equivalence class of $P$. For $U_d(k, n)$, we give SMS the partition $\{[k], [k+1, 2k], [2k+1, n]\}$, and for $U_o(k, n)$, we give the partition $\{[k-1], [k, k], [k+1, 2k-1], [2k, n]\}$. The disadvantage of this adaption is that we are not able to trim down search branches as efficiently, but the advantage is that we can enjoy the speed-up brought by hardcoding a part of the target graph.

## 4.5 External Subgraph Isomorphism Testing

Next we lay out a lazy encoding strategy where we check during the SAT solver's run whether the partially defined graph represented by the current partial assignment to the edge variables can be extended to a $k$-universal graph. We utilize the *Glasgow subgraph solver (GSS)* [19] to achieve this. The GSS is a state-of-the-art solver for the subgraph isomorphism problem based on constraint programming. Given a pattern graph $H$ and a target graph $G$, the subgraph isomorphism problem asks whether there exists a monomorphism from $H$ to $G$. One of the key concepts underlying the GSS's functionality is to generate auxiliary graphs $H_1, \dots, H_m$ with $V(H) = V(H_i), 1 \le i \le m$, and $G_1, \dots, G_m$ with $V(G) = V(G_i), 1 \le i \le m$, such that a mapping from $V(G)$ to $V(H)$ is a (faithful) monomorphisms from $H$ to $G$ if and only if it is a monomorphisms from $H_i$ to $G_i$, for every $1 \le i \le m$. Now the solver searches for a mapping from $V(G)$ to $V(H)$ that is a monomorphism from $H_i$ to $G_i$, for every $1 \le i \le m$. For instance, the GSS can decide whether $H$ is isomorphic to an *induced* subgraph of $G$ by setting $H_1 = H, G_1 = G, H_2 = \bar{H}$, and $G_2 = \bar{G}$.
We utilize this functionality of the GSS to check whether a partially defined graph can be extended to an inducted k-universal graph.

For a partially defined graph $G$ we define two fully defined graphs $G_{-0}$ and $G_{-1}$ by setting
$V(G_{-0}) = V(G_{-1}) = V(G)$,

$$E(G_{-0}) = \{ uv \mid u \neq v \in V(G), uv \in E_d(G) \cup E_u(G) \}, \text{ and}$$

$$E(G_{-1}) = \{ uv \mid u \neq v \in V(G), uv \notin E_d(G) \}. $$

**Proposition 3.** Let $H$ be a graph and $G$ a partially defined graph $G$. $H$ embeds into some graph in $\mathcal{X}(G)$ if and only if there is a mapping from $V(H)$ to $V(G)$ that is a monomorphism both from $H$ to $G_{-0}$ and from $\bar{H}$ to $G_{-1}$,

Thus, we can use the GSS with a minimal change to its interface to check whether $H$ embeds into some graph in $\mathcal{X}(G)$. We integrate the subgraph solver into the external propagator of the SMS framework. Whenever the SAT solver decides on an edge variable $e_{u,v}$, the partially defined graph $G$ represented by the current truth assignment on the edge variables is sent to the propagator. The propagator checks whether all graphs from a set of randomly selected pattern graphs of a certain size are embeddable into $G$. If the check fails, this branch of the search can be terminated: a clause that excludes $G$ is learned and sent back to the SAT solver. We set up three parameters for this process to control the usage of the subgraph solver: *Sample size* specifies the number of randomly selected pattern graphs. *Threshold* specifies the minimal number of determined edges/non-edges that is required to present in $G$ for the embeddability check to be performed. Finally, *frequency* specifies the likelihood the embeddability check is performed and is indicated as an integer $f$. Everytime the SAT solver decides on an edge variable $e_{u,v}$, if the threshold is met, then there is a $1/f$ chance the embeddability check will be performed.

## 4.6 Templates

Next we propose a different approach, based on the concept of a template, which extends the simple method of hardcoding edges or non-edges for speeding up the solving process discussed in Section 4.3.

As it turned out, templates can be defined as partially defined graphs, but with a different purpose than the partially defined graphs used in SMS. Fix a small collection of graphs from $\mathcal{G}(k)$. We want each template to provide a distinct way in terms of how this small collection of graphs can be embedded at simultaneously in the universal graph. Ideally, these templates should allow us to hardcode as many edges/non-edges into the universal graph as possible for the benefit of speed-up. At the same time, there should not be too many of them, as we set a separate SAT instance of each template and need to show that all of them are unsatisfiable to establish a lower bound. In the end, we also want them to be easy to generate.

Let $\mathcal{H} \subseteq \mathcal{G}(k)$. An $(n, \mathcal{H})$-template is a partially defined graph $G$ such that there exists a family of faithful monomorphism $(\delta_H : V(H) \to V(G))_{H \in \mathcal{H}}$ with the additional condition that for all $x, y \in V(G)$, if there is no $H \in \mathcal{H}$ such that $x, y \in \delta_H(V(H))$, then $xy \in E_u(G)$. Let $\mathcal{H}^{(n)}$ denote the set of all $(n, \mathcal{H})$-templates modulo isomorphism.

**Proposition 4.** Let $k < n$ integers and $\mathcal{H} \subseteq \mathcal{G}(k)$. Then $f(k) \le n$ if and only if $\mathcal{X}(G)$ contains a k-universal graph for some $G \in \mathcal{H}^{(n)}$.

According to this proposition, we can limit the search for an induced k-universal graph to the graphs in $\mathcal{X}(G)$ for some $G \in \mathcal{H}^{(n)}$. It remains to generate $\mathcal{H}^{(n)}$ for suitable sets $\mathcal{H}$ for which
39:10 Searching for Smallest Universal Graphs and Tournaments with SAT

$\mathcal{H}^{(n)}$ remains reasonably small. We accomplish that by representing templates (i.e., partially defined graphs) as directed graphs and utilizing SMS for directed graphs as introduced by Kirchweger et al. [16], as follows.

Let us define a mapping $\mathcal{D}$ from the set of partially defined undirected graphs to the set of (fully defined) directed graphs: for each partially defined graph $G$, we put $V(\mathcal{D}(G)) := V(G)$ and $A(\mathcal{D}(G)) := \{uv \in E_d(G) | u > v\} \cup E_u(G)$. We also define the reverse mapping $\mathcal{D}'$ from directed graphs to partially defined undirected graphs: for each directed graph $D$, we put $V(\mathcal{D}'(D)) := V(D)$, $E_d(\mathcal{D}'(D)) := \{uv | uv \in A(D), vu \notin A(D)\} \cup \{uv | vu \in A(D), uv \notin A(D)\}$ and $E_u(\mathcal{D}'(D)) := \{uv | uv, vu \in A(D)\}$. For any partially defined graph $G$ we have that $\mathcal{D}'(\mathcal{D}(G))$ and $G$ are isomorphic.

We encode into CNF the condition for a directed graph $D$ to be $\mathcal{D}(G)$ for some template $G$ of $\mathcal{H}^{(n)}$, use SMS for digraphs to enumerate all such $D$ modulo isomorphism, and then apply $\mathcal{D}'$ to recover the corresponding $G$. We encode the condition for a directed graph $D$ to be $\mathcal{D}(G)$ for some template $G$ of $\mathcal{H}$ as follows.

$$
\begin{align*}
D(\mathcal{H}, n) &:= \bigwedge_{H \in \mathcal{H}} (M(H, n) \wedge \\
& \phantom{:=} \quad \bigwedge_{uv \in E(H)} (m_{u,x} \wedge m_{v,y}) \rightarrow (\neg a_{x,y} \wedge a_{y,x}) \wedge \\
& \phantom{:=} \quad \quad \bigwedge_{\substack{uv \notin E(H) \\ 1 \le x < y \le n}} (m_{u,x} \wedge m_{v,y}) \rightarrow (\neg a_{x,y} \wedge \neg a_{y,x})) \\
& \phantom{:=} \quad \wedge \bigwedge_{1 \le x < y \le n} \left( (a_{x,y} \wedge a_{y,x}) \vee \bigvee_{H \in \mathcal{H}} \left( \bigvee_{i \in V(H)} m_{i,x} \wedge \bigvee_{i \in V(H)} m_{i,y} \right) \right).
\end{align*}
$$

The first three conjuncts in $D(\mathcal{H}, n)$ enforce that $H$ is an induced subgraph of $G$, and the last conjunct ensures that for all pairs $u, v$ of $V(G)$ such that $u, v$ are not both in the range of $\delta_H$ for all $H \in \mathcal{H}$, where $\delta_H$ refers to the faithful monomorphism from $H$ to $G$, we have $uv \in U(G)$.

By trial and error, we determined that an ideal set of pattern graphs to generate templates
for the task of determining $f(7)$ is $\mathcal{H}_4 := \{K_7, \overline{K_7}, K_{3,4}, \overline{K_{3,4}}\}$. Running SMS on $D(\mathcal{H}_4, 16)$
and $D(\mathcal{H}_4, 17)$ gives us 350 and 2772 solutions, respectively. We recover the templates
from these solutions by applying $\mathcal{D}'$. Then, we filter out isomorphic templates using Nauty
with the method explained in the Nauty and Traces User's Guide [20], which converts
detecting isomorphisms between graphs with edge-colorings to that between graphs with
vertex-colorings. This leaves us with 40 and 359 templates, respectively. Hence we have
established the following result.

► **Proposition 5.** $|\mathcal{H}_4^{(16)}| = 40$ and $|\mathcal{H}_4^{(17)}| = 359$.

**5** **Encodings for *k* Universal Tournaments**

**5.1 Basic Encoding**

For a fixed integer $n$, we consider a potential $k$-universal tournament $G$ with $V(G) = [n]$ whose arcs are represented by arc variables $a_{u,v}$, $u \neq v \in [n]$. The truth value of $a_{u,v}$ indicates whether there is an arc from $u$ to $v$.

First, we specify a CNF formula $F'(H, n)$ which is satisfiable if and only if $H$ embeds
into $G$, i.e., if $H$ is isomorphic to a subtournament of $G$. The following CNF formula encodes
the existence of a faithful monomorphism $\phi$ from $H$ to $G$, in a similar way $F(H, n)$ encodes
the universality condition for universal graphs.

$$F'(H, n) := M(H, n) \wedge \bigwedge_{x \neq y \in [n]} ((\neg a_{x,y} \vee \neg a_{y,x}) \wedge (a_{x,y} \vee a_{y,x})) \wedge \bigwedge_{\substack{uv \in A(H) \\ 1 \leq x < y \leq n}} (\neg m_{u,x} \vee \neg m_{v,y} \vee a_{x,y}).$$

$M(H, n)$ ensure that $(m_{u,v})_{u \in V(H), v \in V(G)}$ encodes an injective mapping from $V(H)$ to $V(G)$. The second conjunct ensures that there is exactly one arc between any two distinct vertices. The third conjunct ensures that the mapping preserves arcs. The formula

$$U'(k, n) = \bigwedge_{H \in \mathcal{T}(k)} F'(H, n)$$

is satisfiable if and only if there exists a k-universal tournament on n vertices. From a satisfying assignment we can read off a k-universal tournament.

## 5.2 Approaches to Universal Tournaments

For the SMS approach, we use a recent extension of the framework to directed graphs [16].

For the SAT approach, there are two more features that we add to the basic encoding. The first feature is to hardcode an embedding of the transitive tournament onto $[k]$.

$$\text{TT}(k) = \bigwedge_{1 \leq i < j \leq k} a_{i,j}$$

The second feature is to utilize the perfect symmetry-breaking clauses from Lohn et al.'s work [18] who defined a series of CNF formulas $I(t)$, $1 \leq t \leq 8$, which they call isolators. These isolators restrict arc variables to achieve a perfect symmetry-breaking for tournaments in the sense that the satisfying assignments to $I(t)$ are in 1-1 correspondence with the tournaments of order t modulo isomorphism. Since the size of $I(t)$ grows quickly in t, we only utilize isolators for $t \leq 6$.

Given $l \leq r$, we write $I(l, r)$ to denote the isolator resulting from considering $[l, r]$ instead of $[r - l + 1]$ to be the underlying vertex set. To obtain the overall encoding TTI$(k, n)$, we start with TT$(k)$ and see how many vertices there are in $[k + 1, n]$. If $n - k \leq 6$, then we add $I(k + 1, n)$ to the encoding and finish. If $n - k > 6$, then we add $I(k + 1, k + 6)$ to the encoding and repeat this process with $[k + 7, n]$. For example, TTI$(6, 10) = \text{TT}(6) \wedge I(7, 10)$ and TTI$(7, 16) = \text{TT}(7) \wedge I(8, 13) \wedge I(14, 16)$.

## 6 Certificates

There is a long tradition of computer-assisted proofs in mathematics; a famous example is the proof of the Four-Color Theorem by Apple and Haken [3]. There are mathematical statements whose proof seems to require extensive brute-force search, which is impossible for a human to check [13]. For trusting mathematical statements established by such computational methods, certifying the computational result with a formal proof that can be independently checked is highly desirable. Proof checking should ideally be accomplished through a simple algorithm whose correctness can be trusted or even fully verified. Most of the SAT-based methods described in our paper can be independently verified. In particular, the runs of a CDCL SAT solver on our encodings can produce a DRAT proof, for which a verified proof checker DRAT-trim is available [12, 27]. This approach extends to SMS as follows [15]. During the SMS run, we collect all the generated symmetry-breaking clauses and add them to the clauses generated from the encoding. We can run a plain CDCL SAT solver on this extended set of clauses that generates a DRAT proof. The validity of the symmetry-breaking clauses can also be checked offline with a simple tool.
39:12 Searching for Smallest Universal Graphs and Tournaments with SAT

One could extend the our current verification tool chain in various ways. One can formally verify SAT encodings themselves [5, 6, 9] to provide additional trust in the obtained results. To certify the correctness of auxiliary computations provided by the Glasgow subgraph solver, one can utilize a Cutting Planes proof logging format [10]. What remains is the usage of Nauty [20] for generating pattern graphs and occasionally for filtering isomorphic copies, which does not support independent certification. The former usage, the generation of pattern graphs, can easily be replaced by SMS. For the latter, if an isomorphism between two objects has been found, it is easy to verify whether this is indeed a correct isomorphism. On the other hand, if Nauty fails to identify two objects as being isomorphic, this would only make our overall approach less efficient but would not alter the final results.

# 7 Experiments

## 7.1 Configurations

Based on the techniques discussed in Sections 4 and 5, we consider several configurations that we test in our experiments. We use a naming convention for configurations, where a name starting with UG indicates a configuration for universal graphs and a name starting with UT indicates a configuration for universal tournaments. All UG configurations contain the following basic encoding: the encoding for universality (Section 4.1), the symmetry breaking for twins (Section 4.2), and the hardcoded $K_k$, $\overline{K_k}$ (Section 4.3). We distinguish between the use of SMS and the use of SAT without SMS (plain SAT).

<table><tr><td>UG-SAT</td><td>basic encoding with plain SAT.</td></tr><tr><td>UG-SMS</td><td>basic encoding with SMS.</td></tr><tr><td>UG-SMS-GSS</td><td>basic encoding with SMS extended with the Glasgow Subgraph Solver.</td></tr><tr><td>UG-SAT-Temp</td><td>basic encoding together with templates from Section 4.6 with plain SAT.</td></tr><tr><td>UT-SMS</td><td>encoding for universality from Section 5.1 with SMS.</td></tr><tr><td>UT-SAT-TTI</td><td>encoding for universality with fixed transitive tournament and isolators from Section 5.2 and plain SAT.</td></tr></table>

Theoretically, it is also possible to run SMS extended with GSS without the basic encoding. However, the preliminary experiments we conduct show that this is impractical given the drastically increased running time. Hence, we do not use this method in our main experiments.

## 7.2 Setup

The setup for the experiments is as follows. We generate $G(k)$ and $T(k)$ with Nauty 2.8.6 [20]. For configurations based on SMS, we use a recent SMS version that invoke a recent version of the Cadical SAT solver through the new IPASIR-UP interface [7]. For configurations that do not use SMS we employ Cadical as provided by the PySAT package [14]. We made minor changes to GSS [19] to access its internal functionality. All non-standard tools can be found in the supplementary material.

The experiments are run on a Sun Grid Engine (SGE) with Ubuntu 18.04.6 LTS. The architecture of the nodes in the SGE are among the following: 2x Intel Xeon E5540 with 2.53 GHz Quad Core, 2x Intel Xeon E5649 with 2.53 GHz 6-core, 2x Intel Xeon E5-2630 v2 with 2.60GHz 6-core, 2x Intel Xeon E5-2640 v4 with 2.40GHz 10-core and 2x AMD EPYC 7402 with 2.80GHz 24-core.
## 7.3 Universal Graphs

### 7.3.1 Verifying $f(k)$ and $F(k)$ for $1 \le k \le 5$

We run the configurations UG-SMS and UG-SAT-Temp to verify values $f(k)$ and $F(k)$ for $1 \le k \le 5$ as obtained by Trimble [26]. To verify the values of $f(k)$, we run both configurations on $n = f(k) - 1$ and $n = f(k)$ to see whether the former gives a negative answer (unsat) and the latter a positive (sat). To verify the values of $F(k)$, we solve incrementally, i.e., every time a universal graph is discovered, a clause that prohibits the generation of this particular graph is added to the solver and the solver continues. For UG-SMS, we do not hardcode $K_k$ and $\overline{K_k}$ and start with $([n])$ as the initial partition. This way, the minimality check is enough to filter out non-isomorphic solutions. For UG-SAT-Temp, we detect and exclude isomorphic copies among the found universal graphs using Nauty. It turns out, that both configurations can verify each of the vaules within a couple of seconds.

### 7.3.2 Comparing the Efficiency of Computing $f(6)$

The case of $k = 6$ serves as a good benchmark for testing the efficiency of different methods, since it is neither too hard nor too easy. The value $f(6) = 14$ was determined by Trimble [26]. Hence, we test various configurations on $k = 6, n = 13$ and $k = 6, n = 14$. In the experiments described below, we run all instances 5 times and take the average running time.

For the case of $k = 6, n = 13$, we only test UG-SAT and UG-SMS since this case is computationally easy. The configurations use 0.61s and 1.39s, respectively.

For the case of $k = 6, n = 14$, we test UG-SAT, UG-SMS, and UG-SMS-GSS. The first two configurations use 13.9s and 521s, respectively. With UG-SMS-GSS, we try all possible combinations of the following parameter values: threshold ∈ {40, 50, 60, 70, 80}, frequency ∈ {5, 20, 50, 100, 500, 1000}, and sample size ∈ {5, 20, 60, 100, 156}. We record the total time used in solving and the time used by the embedability check.

<figure><img src="image_378.png" /><figcaption>■ **Figure 3** Data points for UG-SMS-GSS on $k = 6, n = 14$, grouped according to the value of *frequency*. We cast the data points in a two-dimensional space where the *x*-axis represents the percentage of time taken by the embedability check, and the *y*-axis represents the total time used. We draw a red horizontal line standing for the case of SMS, and a black horizontal line standing for the case of pure SAT.</figcaption></figure>
<table><caption>Table 2 Running times for UG-SMS-GSS on the unsatisfiable case with $k = 7$ and $n = 16$.</caption><thead><tr><th>case</th><th>sample size</th><th>frequency</th><th>threshold</th><th>total time/h</th><th>time used in embedability check/h</th></tr></thead><tbody><tr><td>1</td><td>1044</td><td>500</td><td>100</td><td>73.63</td><td>0.12</td></tr><tr><td>2</td><td>1044</td><td>1000</td><td>80</td><td>71.34</td><td>3.40</td></tr><tr><td>3</td><td>1044</td><td>500</td><td>70</td><td>82.36</td><td>11.63</td></tr><tr><td>4</td><td>1044</td><td>1000</td><td>110</td><td>76.34</td><td>0</td></tr></tbody></table>

To understand the parameters' influence on the running time, and how the three configurations compare with each other, we present our findings in Figure 3. We first observe that UG-SAT uses significantly less time than UG-SMS or UG-SMS-GSS. Second, the introduction of the embedability check slightly decreases the running time in some cases and drastically increases the running time in others. Third, the solving time positively correlates with the percentage of the time used by the embedability check. In general, UG-SAT is clearly the most efficient among the configurations and the benefit of UG-SMS-GSS compared to UG-SMS is not obvious.

In Figure 3, data points of different *frequency* are shown in different colors. We see that data points of the same color form loose clusters. A further check into the distribution in terms of sample size and threshold within each group does not yield any meaningful pattern. We also tried grouping the same data points according to *sample size* and *threshold*, but did not obvious any obvious patterns. This indicates that the *frequency* parameter might be the only useful one among the three in terms of tuning for efficiency, and giving it a reasonably large value (which corresponds to the infrequent invocation of the embedability check) is preferred.

We also tried to determine the value of $F(6)$ incrementally with UG-SAT. After a few days we terminated the program. At that state, it had generated 36294 non-isomorphic 6-universal graphs of order 14.

### 7.3.3 $f(7) > 16$

To show $f(7) > 16$, we use UG-SAT, UG-SMS, UG-SMS-GSS, and UG-SAT-Temp. All four configurations conclude that a 7-universal graph with 16 vertices does not exist. UG-SAT takes 39 hours and UG-SMS takes 68 hours. For UG-SMS-GSS, we narrow down the choices of the possible values for the parameters based on previous experiments, and test all combinations with *sample size* = 1044, *frequency* ∈ {500, 1000} and *threshold* {70, 80, 90, 100, 110}. Only 4 test cases out of 50 terminate within 97 hours. Table 2 provides some details. For UG-SAT-Temp, there are $|H_4^{(16)}| = 40$ templates to consider. The easiest template takes 0.4s, the hardest 38 minutes. The total time spent on all 40 templates is 3.6 hours, making the average time spent on each template 322.4s.

Since UG-SAT-Temp seems to be the most effective among the considered configurations, we also tried using it to determine whether there exists a 7-universal graph of order 17. Here we have $|H_4^{(17)}| = 359$ templates to consider. Unfortunately, this approach does not seem to be able to bring the problem down to the reach of modern SAT solvers. By the time of the writing, 25 out of the 359 cases still have not terminated. The 334 cases that have terminated spent 4237 CPU days, all returning with a negative answer. Figure 4 shows the time spent on each template in an ascending order. Even though we have no definitive result, the facts that (i) most cases could be shown unsatisfiable, and (ii) the SAT solver was not able to find a counter example on the remaining cases even after a long runtime, provide some evidence that a 7-universal graph of order 17 might not exist.
<table><caption>Table 3 Overview of running times in seconds for finding one or all solutions up to isomorphism.</caption><thead><tr><th rowspan="3">k</th><th rowspan="3">T(k)</th><th rowspan="3">t(k)</th><th rowspan="3">T(k)</th><th colspan="2">t(k) - 1</th><th colspan="4">t(k)</th></tr><tr><th colspan="2">UT-SAT-TTI</th><th colspan="2">UT-SAT-TTI</th><th colspan="2">UT-SMS</th></tr><tr><th>one</th><th>all</th><th>one</th><th>all</th></tr></thead><tbody><tr><td>5</td><td>12</td><td>8</td><td>1643</td><td>0.03</td><td>0.17</td><td>0.00</td><td>48.28</td><td>0.06</td><td>3.27</td></tr><tr><td>6</td><td>56</td><td>10</td><td>1088</td><td>6.41</td><td>30.68</td><td>11.83</td><td>N/a</td><td>105.36</td><td>47033.06</td></tr></tbody></table>

<figure><img src="image_379.png" /><figcaption>Figure 4 Time spent on each template ordered by their value. The blue bars represent the time spent by templates that have terminated, and the red bars show the time spent on templates that have not terminated by the time of writing.</figcaption></figure>

## 7.4 Universal Tournaments

### 7.4.1 $t(k)$ and $T(k)$ for $1 \le k \le 6$

We run the configurations UT-SMS and UT-SAT-TTI to calculate $t(k)$ and $T(k)$ for $1 \le k \le 6$. To calculate the value of $f(k)$, we start with the lower bound

$$t(k) \ge \max(\{k, t(k-1)\}, \min\{x \mid C_k^x \ge |T(k)|\})$$

where $C_k^x$ denotes the number of $k$-combinations from a given set of $x$ elements, i.e., $C_k^x = \frac{x!}{k!(x-k)!}$. We test TTI($k, n$) while increasing $n$ until we find a universal tournament. To calculate $T(k)$, we modify the two methods similarly as we did for verifying $F(k)$. Table 3 shows the time spent on the case $n = t(k) - 1$ and $n = t(k)$. The time consumed in all applicable cases for $0 \le k \le 4$ are under 0.005s, so we omit them in the table. The process of finding out all non-isomorphic 6-universal tournaments of order 10 with UT-SAT-TTI does not terminate within a couple of days.
### 7.4.2 $13 \le t(7) \le 15$

An obvious lower bound for $t(7)$ is 11 since $C_7^{10} = 120 < 456 = |T(7)|$. Neither UT-SMS nor UT-SAT-TTI gives a result on $k = 7, n = 11$ within half an hour. We generate separate SAT instances where we hardcode each of the four non-isomorphic tournaments of order 4 onto $[8, 11] \subseteq V(G)$ in addition to TTI(7, 11). The four instances terminate with a negative answer in 14.55s, 458.32s, 1399.09s, and 1212.89s, respectively. This way, we have shown that $t(7) > 11$. To show that $t(7) > 12$, we generate and test the following SAT instance for each $\vec{a} = (a_i)_{i \in [7]}, a_i \in [4]$,

$$ \text{TTI}(7, 12) \wedge \bigvee_{i \in [7]} \psi_{a_i}(i) $$

where

$$ \psi_1(i) := |\{j \in [8, 12] \mid e_{i,j} = 1\}| \le 1, \quad \psi_2(i) := |\{j \in [8, 12] \mid e_{i,j} = 1\}| = 2 $$

$$ \psi_3(i) := |\{j \in [8, 12] \mid e_{i,j} = 1\}| = 3, \text{ and } \psi_4(i) := |\{j \in [8, 12] \mid e_{i,j} = 1\}| \ge 4. $$

This gives us 16384 instances in total and all of them return with a negative answer. Thus, $t(7) > 12$. The sum of the time spent is around 2002 CPU days, making the average running time for each case 2.9 hours. More than half of the cases terminate within 40 minutes, and the hardest case takes almost 4 days.

The upper bound $t(7) \le 15$ is obtained with UT-SAT-TTI, which finds a 7-universal tournament of order 15 within 57.95s. As for the cases of TTI(7, 13) and TTI(7, 14), UT-SAT-TTI does not terminate within a few days, and it is not practical to generate separate SAT instances in the way we did for TTI(7, 12), given the large number of instances and the estimated time per case from our preliminary experimentation.

## 8 Conclusion

We have proposed several methods to compute the smallest order of induced $k$-universal graphs and tournaments as a synthesis problem that does not require the enumeration of all potential candidates. Our approaches make use of state-of-the-art CDCL SAT solvers. For the SAT encoding, we proposed a direct encoding and a lazy encoding that utilizes the Glasgow subgraph solver (GSS); for symmetry breaking, we proposed using the SAT modulo Symmetries (SMS) framework and a new templates approach.

We tested and compared these methods, which let us understand their strengths and weaknesses. The template approach was the fastest for unsatisfiable instances; SMS was particularly strong for enumerating all solutions; the GSS was able to speed up the SMS approach. We could improve the known lower bound and show that no induced 7-universal graph with 16 vertices exists, leaving possibilities 17 and 18 for the smallest solution.

We could determine the exact order of smallest $k$-universal tournaments up to $k = 6$ and the interval $[13, 15]$ for $k = 7$. We determined the precise number of optimal solutions, which is of independent interest for combinatorial research.

As future work, determining the exact order of a smallest induced 7-universal graph might be within reach by refining our methods and additional cubing. It would be interesting to see whether it is possible to combine SMS with the templates approach, as this might give a significant boost to possibly allow us to attack the $k$-universality problems for even larger values of $k$.
References

1 Noga Alon. Asymptotically optimal induced universal graphs. *Geom. Funct. Anal.*, 27(1):1–32, 2017. doi:10.1007/s00039-017-0396-9.

2 Stephen Alstrup, Haim Kaplan, Mikkel Thorup, and Uri Zwick. Adjacency labeling schemes and induced-universal graphs. In Rocco A. Servedio and Ronitt Rubinfeld, editors, *Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing, STOC 2015, Portland, OR, USA, June 14-17, 2015*, pages 625–634. ACM, 2015. doi:10.1145/2746539.2746545.

3 Kenneth Appel and Wolfgang Haken. The solution of the four-color-map problem. *Sci. Amer.*, 237(4):108–121, 152, 1977. doi:10.1038/scientificamerican1077-108.

4 Armin Biere, Katalin Fazekas, Mathias Fleury, and Maximillian Heisinger. CaDiCaL, Kissat, Paracooba, Plingeling and Treengeling entering the SAT Competition 2020. In Tomas Balyo, Nils Froleyks, Marijn Heule, Markus Iser, Matti Järvisalo, and Martin Suda, editors, *Proc. of SAT Competition 2020 – Solver and Benchmark Descriptions*, volume B-2020-1 of *Department of Computer Science Report Series B*, pages 51–53. University of Helsinki, 2020. URL: https://tuhat.helsinki.fi/ws/files/142452772/sc2020_proceedings.pdf.

5 Cayden R. Codel. Verifying SAT encodings in lean. Master's thesis, Computer Science Department School of Computer Science Carnegie Mellon University Pittsburgh, PA 15213, May 2022. URL: http://reports-archive.adm.cs.cmu.edu/anon/2022/CMU-CS-22-106.pdf.

6 Luís Cruz-Filipe, João Marques-Silva, and Peter Schneider-Kamp. Formally verifying the solution to the boolean Pythagorean triples problem. *Journal of Automated Reasoning*, 63(3):695–722, 2019. doi:10.1007/s10817-018-9490-4.

7 Katalin Fazekas, Aina Niemetz, Mathias Preiner, Markus Kirchweger, Stefan Szeider, and Armin Biere. IPASIR-UP: User propagators for CDCL. In Meena Mahajan and Friedrich Slivovsky, editors, *The 26th International Conference on Theory and Applications of Satisfiability Testing (SAT 2023), July 04-08, 2023, Alghero, Italy*, LIPIcs. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2023.

8 Johannes K. Fichte, Daniel Le Berre, Markus Hecher, and Stefan Szeider. The silent (r)evolution of SAT. *Communications of the ACM*, 66(6):64–72, June 2023. doi:10.1145/3560469.

9 Sofia Giljegård and Johan Wennerbeck. Puzzle solving with proof—writing a verified SAT encoding chain in HOL4. Master's thesis, Chalmers University of Technology and University of Gothenburg, 2021. URL: https://hdl.handle.net/20.500.12380/304104.

10 Stephan Gocht, Ciaran McCreesh, and Jakob Nordström. Subgraph isomorphism meets cutting planes: Solving with certified solutions. In Christian Bessiere, editor, *Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence, IJCAI 2020*, pages 1134–1140. ijcai.org, 2020. doi:10.24963/ijcai.2020/158.

11 Geňa Hahn and Claude Tardif. Graph homomorphisms: structure and symmetry. In Geňa Hahn and Gert Sabidussi, editors, *Graph Symmetry*, pages 107–166. Kluwer Academic Publishers, Dordrecht, 1997.

12 Marijn Heule, Warren A. Hunt Jr., Matt Kaufmann, and Nathan Wetzler. Efficient, verified checking of propositional proofs. In Mauricio Ayala-Rincón and César A. Muñoz, editors, *Interactive Theorem Proving - 8th International Conference, ITP 2017, Brasília, Brazil, September 26-29, 2017, Proceedings*, volume 10499 of *Lecture Notes in Computer Science*, pages 269–284. Springer Verlag, 2017. doi:10.1007/978-3-319-66107-0_18.

13 Marijn J. H. Heule and Oliver Kullmann. The science of brute force. *Communications of the ACM*, 60(8):70–79, 2017. doi:10.1145/3107239.

14 Alexey Ignatiev, Antonio Morgado, and Joao Marques-Silva. PySAT: A Python toolkit for prototyping with SAT oracles. In *SAT*, pages 428–437, 2018. doi:10.1007/978-3-319-94144-8_26.

15 Markus Kirchweger, Manfred Scheucher, and Stefan Szeider. A SAT attack on Rota's Basis Conjecture. In *Theory and Applications of Satisfiability Testing - SAT 2022 - 25th International Conference, Haifa, Israel, August 2-5, 2022, Proceedings*, 2022. doi:10.4230/LIPIcs.SAT.2022.4.
39:18 Searching for Smallest Universal Graphs and Tournaments with SAT

16 Markus Kirchweger, Manfred Scheucher, and Stefan Szeider. SAT-based generation of planar graphs. In Meena Mahajan and Friedrich Slivovsky, editors, *The 26th International Conference on Theory and Applications of Satisfiability Testing (SAT 2023), July 04-08, 2023*, Alghero, Italy, LIPIcs. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2023.

17 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 2021)*, Leibniz International Proceedings in Informatics (LIPIcs), pages 39:1-39:17. Dagstuhl, 2021. doi:10.4230/LIPICs.CP.2021.34.

18 Evan Lohn, Chris Lambert, and Marijn J. H. Heule. Compact symmetry breaking for tournaments. In Alberto Griggio and Neha Rungta, editors, *22nd Formal Methods in Computer-Aided Design, FMCAD 2022, Trento, Italy, October 17-21, 2022*, pages 179-188. IEEE, 2022. doi:10.34727/2022/isbn.978-3-85448-053-2_24.

19 Ciaran McCreesh, Patrick Prosser, and James Trimble. The Glasgow subgraph solver: Using constraint programming to tackle hard subgraph isomorphism problem variants. In Fabio Gadducci and Timo Kehrer, editors, *Graph Transformation - 13th International Conference, ICGT 2020, Held as Part of STAF 2020, Bergen, Norway, June 25-26, 2020, Proceedings*, volume 12150 of *Lecture Notes in Computer Science*, pages 316-324. Springer, 2020. doi:10.1007/978-3-030-51372-6_19.

20 Brendan D. McKay and Adolfo Piperno. Practical graph isomorphism, II. *J. Symbolic Comput.*, 60:94-112, 2014. doi:10.1016/j.jsc.2013.09.003.

21 Jhon W. Moon. On minimal *n*-universal graphs. *Proc. Glasgow Math. Assoc.*, 7:32-33 (1965), 1965. doi:10.1017/S2040618500035139.

22 John W. Moon. *Topics on tournaments*. Holt, Rinehart and Winston, New York-Montreal, Que.-London, 1968.

23 OEIS Foundation Inc. The On-Line Encyclopedia of Integer Sequences. Published electronically at http://oeis.org.

24 James Preen. What is the smallest graph that contains all non-isomorphic 4-node and 5-node connected graphs as induced subgraphs? online, October 2011. Mathematics Stack Exchange, https://math stackexchange com/q/75513, last edited 2020-04-25.

25 Richard Rado. Universal graphs and universal functions. *Acta Arithmetica*, 9(4):331-340, 1964. URL: http://eudml.org/doc/207488.

26 James Trimble. Induced universal graphs for families of small graphs. *CoRR*, abs/2109.00075, 2021. doi:10.48550/arXiv.2109.00075.

27 Nathan Wetzler, Marijn J. H. Heule, and Warren A. Hunt. DRAT-trim: Efficient checking and trimming using expressive clausal proofs. In *Theory and Applications of Satisfiability Testing - SAT 2014*, volume 8561 of *Lecture Notes in Computer Science*, pages 422-429. Springer Verlag, 2014. doi:10.1007/978-3-319-09284-3_31.
A Appendix

<figure><img src="image_380.png" /><figcaption>Figure 5 Smallest induced *k*-universal graphs and their incidence matrices.</figcaption></figure>

<figure><img src="image_381.png" /><figcaption>Figure 6 All smallest *k*-universal tournaments for *k* = 1, 2, 3, 4, and their incidence matrices.</figcaption></figure>
<figure><img src="image_382.png" /><figcaption>■ **Figure 7** The three displayed tournaments are examples of k-universal tournaments of order 8, 10, and 15 for $k = 5$, $k = 6$, and $k = 7$, respectively. We found many more with this property. We also show their incidence matrices. For $k = 5, 6$ we know that there are no k-universal tournaments of smaller order, for $k = 7$ we cannot rule out that there is one of order 14.</figcaption></figure>
FastMapSVM for Predicting CSP Satisfiability

Kexin Zheng<fnref n="1" />

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

University of Southern California, Los Angeles, CA, USA

Ang Li<fnref n="1" />

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

University of Southern California, Los Angeles, CA, USA

Han Zhang

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

University of Southern California, Los Angeles, CA, USA

T. K. Satish Kumar

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

University of Southern California, Los Angeles, CA, USA

## Abstract

Recognizing the satisfiability of Constraint Satisfaction Problems (CSPs) is NP-hard. Although several Machine Learning (ML) approaches have attempted this task by casting it as a binary classification problem, they have had only limited success for a variety of challenging reasons. First, the NP-hardness of the task does not make it amenable to straightforward approaches. Second, CSPs come in various forms and sizes while many ML algorithms impose the same form and size on their training and test instances. Third, the representation of a CSP instance is not unique since the variables and their domain values are unordered. In this paper, we propose FastMapSVM, a recently developed ML framework that leverages a distance function between pairs of objects. We define a novel distance function between two CSP instances using maxflow computations. This distance function is well defined for CSPs of different sizes. It is also invariant to the ordering on the variables and their domain values. Therefore, our framework has broader applicability compared to other approaches. We discuss various representational and combinatorial advantages of FastMapSVM. Through experiments, we also show that it outperforms other state-of-the-art ML approaches.

**2012 ACM Subject Classification** Computing methodologies → Machine learning

**Keywords and phrases** Constraint Satisfaction Problems, Machine Learning, FastMapSVM

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

**Funding** This work at the University of Southern California is supported by DARPA under grant number HR001120C0157 and by NSF under grant number 2112533. The views, opinions, and/or findings expressed are those of the author(s) and should not be interpreted as representing the official views or policies of the sponsoring organizations, agencies, or the U.S. Government.

# 1 Introduction

Constraints constitute a very natural and general means for formulating regularities in the real world. A fundamental combinatorial structure used for reasoning with constraints is that of the Constraint Satisfaction Problem (CSP). The CSP formally models a set of variables, their corresponding domains, and a collection of constraints between subsets of the variables. Each constraint restricts the set of allowed combinations of values of the participating variables. A solution of a given CSP instance is an assignment of values to all the variables from their respective domains such that all the constraints are satisfied. Technologies for efficiently solving CSPs bear immediate and important implications on how fast we can solve computational problems that arise in several other areas of research, including computer vision, spatial and temporal reasoning, model-based diagnosis, planning and scheduling, and language understanding.

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

<footnote n="1">¹ Corresponding Author</footnote>
Unfortunately, solving CSPs is NP-hard since they generalize the Satisfiability (SAT) problem. Although many technologies have been developed for solving CSPs in practice [4], they do not sufficiently harness the power of Machine Learning (ML) techniques. While there have been a lot of attempts to apply ML techniques to CSPs, none of these attempts have yielded spectacular results: They do not consistently produce high-quality outcomes. Examples of ML approaches used in the CSP domain include the application of Support Vector Machines (SVMs) [2], linear regression [22], decision tree learning [7, 6], clustering [15, 8], *k*-nearest neighbors [13], and others [9]. However, most of these approaches have had limited success for a variety of challenging reasons.

First, from a complexity theory perspective, the NP-hardness of the task does not make it amenable to straightforward ML approaches. For example, since a neural network (NN) is essentially a continuous differentiable form of a circuit, it is not straightforward to make NNs effective in the CSP domain. Second, CSPs come in various forms and sizes while many ML algorithms use an architectural framework that imposes the same form and size on their training and test instances. For example, an NN may have a fixed input layer that it uses for the training and test instances alike. Third, the representation of a CSP instance is not unique since the variables and their domain values are unordered. This poses a significant combinatorial challenge for ML algorithms since they have to learn the permutation invariance with respect to orderings on the variables and their domain values. For example, an NN may pose the overhead of having to be trained on all permutations of the same CSP instance to become effective.

In this paper, we consider the problem of predicting the satisfiability of CSP instances using ML. In ML terminology, this is essentially a binary classification problem defined on CSPs with the two possible classification labels “satisfiable” and “unsatisfiable”. This classification problem is a cornerstone task for addressing the combinatorics of CSPs. It also serves as a stepping stone for the task of solving CSPs. In fact, any ML framework expected to be viable for solving CSPs should likely first demonstrate its success on solving the aforementioned classification problem on CSPs.

We propose to solve the above classification problem on binary CSPs<fnref n="2" /> using a recently developed ML framework called FastMapSVM [19]. While most ML algorithms learn to identify characteristic features of individual objects in a class, FastMapSVM leverages a domain-specific distance function on pairs of objects. It does this by combining the strengths of FastMap [5] and SVMs. In its first stage, FastMapSVM invokes FastMap, an efficient linear-time algorithm that maps complex objects to points in a Euclidean space, while preserving pairwise distances between them. In its second stage, it invokes SVMs and kernel methods for learning to classify the points in this Euclidean space.

FastMapSVM has demonstrated success on classifying complex objects such as seismograms in Earthquake Science [19]. It offers several advantages over ML algorithms that reason about individual objects instead of pairs of objects. First, FastMapSVM enables domain experts to incorporate their domain knowledge using a distance function. This avoids relying on complex ML models to infer the underlying structure in the data entirely. Second, because the distance function encapsulates domain knowledge, FastMapSVM naturally facilitates interpretability and explainability. In fact, it even provides a perspicuous visualization of the objects and the classification boundaries between them. Third, FastMapSVM uses significantly smaller amounts of time and data for model training compared to other ML algorithms. Fourth, it extends the applicability of SVMs and kernel methods to domains with complex objects.

<footnote n="2">² Binary CSPs have at most two variables per constraint but are allowed to have non-Boolean variables. Binary CSPs are representationally as powerful as general CSPs with bounded arity of the constraints.</footnote>
<figure><img src="image_388.png" /><figcaption>Figure 1 The left panel shows the (0, 1)-matrix representation of a single constraint $C(X_1, X_2)$. The right panel shows the (0, 1)-matrix representation of an entire binary CSP instance.</figcaption></figure>

In applying FastMapSVM to the CSP domain, we define a novel distance function between two CSP instances. This distance function uses maxflow computations and is well defined for CSP instances of different sizes. It is also invariant to the ordering on the variables and their domain values in the CSP instances. Therefore, FastMapSVM has broader applicability compared to other ML approaches in the CSP domain. Moreover, since it uses the intelligence of SVMs, kernel methods, and maxflow computations, it is able to significantly outperform competing ML approaches. It is also able to outperform procedures that invest polynomial time in establishing local consistency—such as arc-consistency—to discover unsatisfiable CSP instances. This demonstrates that a trained FastMapSVM model acquires an intelligence beyond that of prominent polynomial-time procedures.<fnref n="3" /> We discuss various other representational and combinatorial advantages of FastMapSVM and, through experiments, we also demonstrate its superior performance.

## 2 Preliminaries and Definitions

A CSP instance is defined by a triplet $(\mathcal{X}, \mathcal{D}, \mathcal{C})$, where $\mathcal{X} = \{X_1, X_2 \dots X_N\}$ is a set of variables and $\mathcal{C} = \{C_1, C_2 \dots C_M\}$ is a set of constraints on subsets of them. Each variable $X_i$ is associated with a finite discrete-valued domain $D_i \in \mathcal{D}$, and each constraint $C_i$ is a pair $(S_i, R_i)$ defined on a subset of variables $S_i \subseteq \mathcal{X}$, called the *scope of* $C_i$. $|S_i|$ is referred to as the *arity* of the constraint. $R_i \subseteq D_{S_i}$ ($D_{S_i} = \times_{X_j \in S} D_j$) denotes all compatible tuples of $D_{S_i}$ allowed by the constraint. The absence of a constraint on a certain subset of the variables is equivalent to a constraint on the same subset of the variables that allows all combinations of values to them. A *solution* of a CSP instance is an assignment of values to all the variables from their respective domains such that all the constraints are satisfied. A binary CSP instance has at most two variables per constraint. Binary CSPs are representationally as powerful as general CSPs with bounded arity of the constraints [4].

A binary CSP is *arc-consistent* if and only if for all variables $X_i$ and $X_j$, and for every instantiation of $X_i$, there exists an instantiation of $X_j$ such that the direct constraint between them is satisfied. Similarly, a binary CSP is *path-consistent* if and only if for all variables

<footnote n="3">³ This is an important hallmark of an ML algorithm. In [20], a deep NN model is presented to recognize the satisfiability of CSP instances with Boolean variables and binary constraints. However, this class of CSP instances is equivalent to 2-SAT and can be solved in polynomial time, diminishing the advantages of an ML framework over polynomial-time reasoning.</footnote>
<figure><img src="image_389.png" /></figure>

<figure><img src="image_390.png" /><figcaption>■ **Figure 2** The figure, borrowed from [3], illustrates how coordinates are computed and recursion is carried out in FastMap.</figcaption></figure>

$X_i$, $X_j$ and $X_k$, and for every instantiation of $X_i$ and $X_j$ that satisfies the direct constraint between them, there exists an instantiation of $X_k$ such that the direct constraints between $X_i$ and $X_k$ and between $X_j$ and $X_k$ are also satisfied.

For a given binary CSP instance, we can build a matrix representation for it using a simple mechanism. First, we assume that the domain values of each variable are ordered in some way. (We can simply use the order in which the domain values of each of the variables are specified.) Under such an ordering, we can represent each binary constraint as a 2-dimensional matrix with all its entries set to either 1 or 0 based on whether the corresponding combination of values to the participating variables is allowed or not by that constraint. The left panel of Figure 1 shows the (0,1)-matrix representation of a binary constraint between two variables $X_1$ and $X_2$ with domain sizes of 5 each. The combination of values ($X_1 \leftarrow d_{12}$, $X_2 \leftarrow d_{21}$) is an allowed combination, and the corresponding entry in the matrix is therefore set to 1. However, the combination of values ($X_1 \leftarrow d_{14}$, $X_2 \leftarrow d_{22}$) is a disallowed combination, and the corresponding entry is therefore set to 0. In general, $d_{ip}$ denotes the p-th domain value of $X_i$ assuming an index ordering on the domain values of $X_i$.

The (0,1)-matrix representation of an entire binary CSP instance can be constructed simply by stacking up the matrix representations for the individual constraints into a bigger “block” matrix. The right panel of Figure 1 illustrates how a binary CSP instance on 3 variables $X_1$, $X_2$, and $X_3$ can be represented as a “mega-matrix” with 3 sets of rows and 3 sets of columns. Each block-entry inside this mega-matrix is the matrix representation of the direct constraint between the corresponding row and column variables. Therefore, the matrix representation of an entire binary CSP instance has $\sum_{i=1}^N |D_i|$ rows and $\sum_{i=1}^N |D_i|$ columns.

3 FastMap and FastMapSVM

In this section, we describe FastMap and FastMapSVM to set up the groundwork for our approach. Both these rely on a domain-specific distance function $D(\cdot, \cdot)$ on pairs of objects. $D(\cdot, \cdot)$ is required to be a non-negative symmetric function.
## 3.1 FastMap

FastMap [5] is a Data Mining algorithm that embeds complex objects—like audio signals, seismograms, DNA sequences, electrocardiograms, or magnetic-resonance images—into a K-dimensional Euclidean space, for a user-specified value of K and a user-supplied distance function $D(\cdot, \cdot)$ on pairs of objects. The Euclidean distance between any two objects in the embedding approximates the domain-specific distance between them. Therefore, similar objects, as quantified by $D(\cdot, \cdot)$, map to nearby points in Euclidean space while dissimilar objects map to distant points. Although FastMap preserves $O(N^2)$ pairwise distances between N objects, it generates the embedding in only $O(KN)$ time. Because of its efficiency, FastMap has already found numerous real-world applications, including in Data Mining [5], shortest-path computations [3], solving combinatorial optimization problems on graphs [12], and community detection and block modeling [11].

FastMap embeds a collection of complex objects in an artificially created Euclidean space that enables geometric interpretations, algebraic manipulations, and downstream ML algorithms. It gets as input a collection of complex objects $\mathcal{O}$, where $D(O_i, O_j)$ represents the domain-specific distance between objects $O_i, O_j \in \mathcal{O}$. It generates a Euclidean embedding that assigns a K-dimensional point $\mathbf{p}_i \in \mathbb{R}^K$ to each object $O_i$. A good Euclidean embedding is one in which the Euclidean distance $\|\mathbf{p}_j - \mathbf{p}_i\|_2$ between any two points $\mathbf{p}_i$ and $\mathbf{p}_j$ closely approximates $D(O_i, O_j)$.

In the first iteration, FastMap heuristically identifies the farthest pair of objects $O_a$ and $O_b$ in linear time. Once $O_a$ and $O_b$ are determined, every other object $O_i$ defines a triangle with sides of lengths $d_{ai} = D(O_a, O_i)$, $d_{ab} = D(O_a, O_b)$, and $d_{ib} = D(O_i, O_b)$, as illustrated in Figure 2a. The lengths of the sides of the triangle define its entire geometry, and the projection of $O_i$ onto the line $\overline{O_aO_b}$ is given by

$$x_i = (d_{ai}^2 + d_{ab}^2 - d_{ib}^2)/(2d_{ab}). \qquad (1)$$

FastMap sets the first coordinate of $\mathbf{p}_i$, the embedding of $O_i$, equal to $x_i$. In the subsequent K − 1 iterations, the same procedure is followed for computing the remaining K − 1 coordinates of each object; however, the distance function is adapted for each iteration. For example, for the first iteration, the coordinates of $O_a$ and $O_b$ are 0 and $d_{ab}$, respectively. Because these coordinates fully explain the true distance between these two objects, from the second iteration onward, the rest of $\mathbf{p}_a$ and $\mathbf{p}_b$'s coordinates should be identical. Intuitively, this means that the second iteration should mimic the first one on a hyperplane that is perpendicular to the line $\overline{O_aO_b}$. Figure 2b illustrates this. Although the hyperplane is never explicitly constructed, it conceptually implies that the distance function for the second iteration should be changed for all i and j in the following way:

$$D_{new}(O'_i, O'_j)^2 = D(O_i, O_j)^2 - (x_i - x_j)^2, \qquad (2)$$

in which $O'_i$ and $O'_j$ are the projections of $O_i$ and $O_j$, respectively, onto this hyperplane, and $D_{new}(\cdot, \cdot)$ is the new distance function. The distance function is recursively updated according to Equation 2 at the beginning of each of the K − 1 iterations that follow the first.

In each of the K iterations, FastMap heuristically finds the farthest pair of objects according to the distance function defined for that iteration. These objects are called pivots and are stored as reference objects. There are very few, that is, $\le 2K$, reference objects. Technically, finding the farthest pair of objects in any iteration takes $O(N^2)$ time. However, FastMap uses a linear-time “pivot changing” heuristic [5] to efficiently and effectively identify a pair of objects $O_a$ and $O_b$ that is very often the farthest pair. It does this by initially
choosing a random object $O_b$ and then choosing $O_a$ to be the farthest object away from $O_b$. It then reassigns $O_b$ to be the farthest object away from $O_a$, reassigns $O_a$ to be the farthest object away from $O_b$, and so on, until convergence or a maximum of *C* iterations, for a small constant $C \le 10$.

## 3.2 FastMapSVM

FastMapSVM [19] elegantly combines the strengths of FastMap and SVMs. In the first phase, it invokes FastMap for efficiently mapping complex objects to points in a Euclidean space, while preserving pairwise distances between them. In the second phase, it invokes SVMs and kernel methods for learning to classify the points in this Euclidean space. FastMapSVM has several advantages over other methods.

First, FastMapSVM leverages domain-specific knowledge via a distance function. There are many real-world domains in which feature selection for *individual* objects is hard. While a domain expert can occasionally identify and incorporate domain-specific features of the objects to be classified, doing so becomes increasingly hard with increasing complexity of the objects. Therefore, many existing ML algorithms for classification find it hard to leverage domain knowledge when used off the shelf. However, in many real-world domains with complex objects, a distance function on *pairs* of objects is well defined and easy to compute. In such domains, FastMapSVM is more easily applicable than other ML algorithms that focus on the features of individual objects. FastMapSVM also enables domain experts to incorporate their domain knowledge via a distance function instead of relying on complex ML models to infer the underlying structure in the data entirely. Examples of such real-world objects include audio signals, seismograms, DNA sequences, electrocardiograms, and magnetic-resonance images. While these objects are complex and may have many subtle features that are hard to recognize, there exists a well-defined distance function on *pairs* of objects that is easy to compute. For instance, individual DNA sequences have many complex and subtle features but the *edit distance*<fnref n="4" /> between two DNA sequences is well defined and easy to compute. Similarly, the Minkowski distance [1] is well defined for images and the cosine similarity [16] is well defined for text documents.

Second, FastMapSVM facilitates interpretability, explainability, and visualization. Many existing ML algorithms produce results that are hard to interpret or explain. For example, in NNs, a large number of interactions between neurons with nonlinearities makes a meaningful interpretation or explanation of the results very hard. In fact, the very complexity of the objects in the domain can hinder interpretability and explainability. FastMapSVM mitigates these challenges and thereby supports interpretability and explainability. While the objects themselves may be complex, FastMapSVM embeds them in a Euclidean space by considering only the distance function defined on *pairs* of objects. In effect, it simplifies the description of the objects by assigning Euclidean coordinates to them. Moreover, since the distance function is itself user-supplied and encapsulates domain knowledge, FastMapSVM naturally facilitates interpretability and explainability. In fact, it even provides a perspicuous visualization of the objects and the classification boundaries between them. This aids human interpretation of the data and results. It also enables a human-in-the-loop framework for refining the processes of learning and decision making. As a hallmark, FastMapSVM produces the visualization very efficiently since it invests only linear time in generating the Euclidean embedding.

<footnote n="4">⁴ The edit distance between two strings is the minimum number of insertions, deletions, or substitutions that are needed to transform one to the other.</footnote>
Third, FastMapSVM uses significantly smaller amounts of time and data for model training compared to other ML algorithms. While NNs and other ML algorithms store abstract representations of the training data in their model parameters, FastMapSVM stores explicit references to some of the original objects, referred to as pivots. While making predictions, objects in the test instances are compared directly to the pivots using the user-supplied distance function. Thereby, FastMapSVM obviates the need to learn a complex transformation of the input data and thus significantly reduces the amounts of time and data required for model training. Moreover, given $N$ training instances, that is, $N$ objects and their classification labels, FastMapSVM leverages $O(N^2)$ pieces of information via the distance function that is defined on every pair of objects. In contrast, ML algorithms that focus on individual objects leverage only $O(N)$ pieces of information.

Fourth, FastMapSVM extends the applicability of SVMs and kernel methods to complex objects. Generally speaking, SVMs are particularly good for classification tasks. When combined with kernel methods, they recognize and represent complex nonlinear classification boundaries very elegantly [17]. Moreover, soft-margin SVMs with kernel methods [14] can be used to recognize both outliers and inherent nonlinearities in the data. While the SVM machinery is very effective, it requires the objects in the classification task to be represented as points in a Euclidean space. Often, it is very difficult to represent complex objects as precise geometric points without introducing inaccuracy or losing domain-specific representational features. In such cases, deep NNs have gained more popularity compared to SVMs for the reason that it is unwieldy for SVMs to represent all the features of complex objects in Euclidean space. However, FastMapSVM revives the SVM approach by leveraging a distance function and creating a low-dimensional Euclidean embedding of the complex objects.

# 4 Distance Function on CSPs

In this section, we describe a distance function on binary CSPs. This distance function is based on maxflow computations and is illustrated in Figure 3. It is well defined for CSP instances $I_1$ and $I_2$ that may have different sizes. The maxflow computations are utilized in: (a) a single high-level “maximum matching of minimum cost” problem posed on the variables of $I_1$ and $I_2$, and (b) multiple low-level “maximum matching of minimum cost” problems posed on the domain values of pairs of variables, one from each of $I_1$ and $I_2$.

The high-level “maximum matching of minimum cost” problem is posed on a complete bipartite graph, in which the two partitions of the bipartite graph correspond to the variables of $I_1$ and $I_2$, respectively. If the number of variables in $I_1$ does not match the number of variables in $I_2$, dummy variables are added to the CSP instance with fewer variables. Figure 3 (top panel) illustrates this for $I_1$ and $I_2$ with variables $\{X_1, X_2, X_3, X_4\}$ and $\{X'_1, X'_2, X'_3\}$, respectively. A dummy variable $X'_4$ is added to $I_2$. The dummy variable has a single domain value that is designed to be consistent with all domain values of all other variables, since this does not change the CSP instance.

The distance between $I_1$ and $I_2$ is defined to be the cost of the “maximum matching of minimum cost” on the high-level bipartite graph. This bipartite graph has an edge between every $X_i$ in $I_1$ and every $X'_j$ in $I_2$. The cost annotating an edge between $X_i$ and $X'_j$ is itself set to be the cost of the “maximum matching of minimum cost” posed at the low level on the domain values of $X_i$ and $X'_j$. Figure 3 (bottom-left panel) shows the high-level bipartite graph and highlights an edge between $X_1$ and $X'_2$ for explanation of the low-level “maximum matching of minimum cost”.
40:8 FastMapSVM for Predicting CSP Satisfiability

<figure><img src="image_391.png" /><figcaption>■ **Figure 3** The top panel shows two CSP instances with variables {$X_1, X_2, X_3, X_4$} (left) and {$X'_1, X'_2, X'_3$} (middle), respectively. A dummy variable $X'_4$ with a singleton domain is added to the CSP instance with fewer variables (right). The bottom panel (left) shows how a “maximum matching of minimum cost” problem is posed on a complete bipartite graph with the variables of the two CSP instances in each partition. The cost annotating the edge between $X_i$ and $X'_j$ is itself derived from a “maximum matching of minimum cost” problem posed on the domain values of $X_i$ and $X'_j$. The bottom panel (right) shows this “maximum matching of minimum cost” problem for the variables $X_1$ and $X'_2$. It is posed on a complete bipartite graph with the domain values {$d_{11}, d_{12}, d_{13}$} and {$d'_{21}, d'_{22}, d'_{23}, d'_{24}, d'_{25}, d'_{26}$} in each partition. The cost annotating the edge between $d_{11}$ and $d'_{21}$ is the absolute value of the difference between the average compatibility of $d_{11}$ and the average compatibility of $d'_{21}$.</figcaption></figure>

The low-level “maximum matching of minimum cost” problem posed on the domain values of $X_i$ and $X'_j$ also uses a complete bipartite graph. The two partitions consist of the domain values of $X_i$ and $X'_j$, respectively. The cost annotating the edge between $d_{ip}$ and $d'_{jq}$ is the absolute value of the difference between the average compatibility of $d_{ip}$ and the average compatibility of $d'_{jq}$. Figure 3 (bottom-right panel) shows the low-level “maximum matching of minimum cost” problem posed on the domain values of $X_1$ and $X'_2$. The domains of $X_1$ and $X'_2$ are {$d_{11}, d_{12}, d_{13}$} and {$d'_{21}, d'_{22}, d'_{23}, d'_{24}, d'_{25}, d'_{26}$}, respectively. Consider the edge between $d_{11}$ and $d'_{21}$. The average compatibility of $d_{11}$ is the fraction of “1”s in the column “$d_{11}$” in the matrix representation of $\mathcal{I}_1$. This fraction is equal to 8/16. The average compatibility of $d'_{21}$ is the fraction of “1”s in the column “$d'_{21}$” in the matrix representation of $\mathcal{I}_2$ after adding the dummy variable $X'_4$. This fraction is equal to 4/11. Therefore, the cost annotating the edge between $d_{11}$ and $d'_{21}$ is equal to $|8/16 - 4/11|$.

The “maximum matching of minimum cost” problems in the high level and the low level are posed on bipartite graphs. Since the costs annotating the edges of the bipartite graphs in the high level and the low level are non-negative, the distance function is also non-negative.
Moreover, since the two partitions of any bipartite graph can be viewed interchangeably without affecting the “maximum matching of minimum cost”, the overall distance function is symmetric. It is also easy to observe that the distance between two identical CSP instances is always 0. These properties of the distance function satisfy all the requirements imposed on it by the FastMap component of FastMapSVM.

The high-level bipartite graph is invariant to the orderings on the elements within each partition. That is, it is invariant to the orderings on the variables of $\mathcal{I}_1$ and $\mathcal{I}_2$. Similarly, all low-level bipartite graphs are invariant to the orderings on the domain values of the participating variables. Therefore, the overall distance function is invariant to variable-orderings as well as domain value-orderings. This allows us to bypass data augmentation methods typically required for training other ML models.<fnref n="5" /> In the context of CSPs, a CSP instance is typically augmented by changing the ordering on its variables or the ordering on the domain values of individual variables. However, doing so generates an exponential number of CSP training instances within the same equivalence class. This drawback of traditional ML algorithms of having to learn equivalence classes is now intelligently addressed within the framework of FastMapSVM by utilizing a distance function that is invariant to both variable-orderings and domain value-orderings.

We note that the above distance function could have been defined in many other ways. For example, we could have introduced dummy domain values in the low-level “maximum matching of minimum cost” problems to equalize the domain sizes of the participating variables. We could have also chosen not to use dummy variables in the high-level “maximum matching of minimum cost” problem. In addition, we could have defined the costs annotating the edges of the bipartite graphs using many other characteristics of the CSPs. These variations of the distance function are not of fundamental importance to this paper. Instead, in this paper, we focus on the advantages of the FastMapSVM framework as a whole. The study of more refined distance functions is delegated to future work.

# 5 Experimental Results

In this section, we describe the comparative performance of FastMapSVM against other state-of-the-art ML approaches on predicting CSP satisfiability.

## 5.1 Experimental Setup

We evaluate FastMapSVM against three competing approaches. The first is a state-of-the-art deep graph convolutional neural network (DGCNN) [23]. The second is a state-of-the-art graph isomorphism network (GIN) [21]. Both these networks ingest a CSP instance in the form of a graph, as shown in Figure 4. In the graphical representation of a binary CSP instance, a vertex represents a domain value and is tagged with the name of the variable that it belongs to. Information in these tags is utilized by the DGCNN and the GIN. An edge between two vertices $v_1$ and $v_2$ with tags $X_i$ and $X_j$, respectively, represents the compatible combination ($X_i \leftarrow v_1, X_j \leftarrow v_2$) allowed by the direct constraint between $X_i$ and $X_j$.⁶ DGCNN and GIN do not require the CSP training and test instances to be of the same size.

The third is a polynomial-time algorithm based on establishing arc-consistency. This algorithm first establishes arc-consistency and then checks whether any variable's domain is annihilated. If so, it declares the CSP instance to be “unsatisfiable”. Otherwise, it declares

<footnote n="5">⁵ Data augmentation refers to transformations of data without changing their labels, known as label-preserving transformations. For example, to generate more training data serving object recognition tasks in computer vision applications, an image can be augmented by translating it or reflecting it</footnote>
40:10 FastMapSVM for Predicting CSP Satisfiability

<figure><img src="image_392.png" /><figcaption>Figure 4 The left panel shows the matrix representation of a binary CSP instance. The right panel shows its graphical representation. The vertices represent domain values and are clustered into four groups, corresponding to the four variables {$X_1, X_2, X_3, X_4$}.</figcaption></figure>

the CSP instance to be “satisfiable”. This algorithm is used in our evaluation to demonstrate that FastMapSVM’s capabilities go beyond that of a polynomial-time algorithm.⁷ Of course, a polynomial-time algorithm based on establishing path-consistency could also have been used. But we excluded this algorithm since arc-consistency already provides the required proof of concept and establishing path-consistency is prohibitively expensive.

In our experiments, we do not include results from the CSP-cNN framework of [20] for the following reasons. First, this framework has the drawback that it is applicable only to CSP training and test instances of the same size. (FastMapSVM does not have this drawback since the distance function does not require the CSP instances to be of the same size.) Second, the success of CSP-cNN critically depends on data augmentation methods: For each CSP training instance, a very large number of other training instances that permute the variables and their domain values also have to be generated. (This requirement is completely obviated by FastMapSVM since the distance function is invariant to such permutations.) Third, CSP-cNN has been shown to be successful only on Boolean binary CSP instances, that is, the polynomial-time solvable 2-SAT problems. (FastMapSVM does not have this limitation; below, we demonstrate its success on general binary CSP instances.)

We implemented FastMapSVM and arc-consistency in Python3 and ran them on a laptop with an Apple M2 chip with 16 GB memory. We ran DGCNN and GIN on a Linux system with an Intel(R) Xeon(R) Silver 4216 CPU at 2.10 GHz. The different platforms are inconsequential to the comparative performances of these algorithms with respect to effectiveness. For each dataset, we trained DGCNN and GIN for 100 epochs with a learning rate of 0.0001 and a minibatch size of 100 to obtain representative results.

## 5.2 Instance Generation

We generate the binary CSP instances for both training and testing using the Model A method in [18, 20]. We generate a CSP instance by first picking the number of variables *N* uniformly at random to be an integer within the range [1, 100]. Then, we pick the domain

<footnote>horizontally without changing its label [10].

⁶ The graphical representation of a binary CSP instance is obtained by parsing its matrix representation. Thus, we correctly represent the compatible tuples of domain values between every pair of variables, even if there does not exist a direct constraint between those variables.

⁷ This is done to avoid the pitfalls of [20], as mentioned before.</footnote>
size of each variable independently and uniformly at random to be an integer within the range [1, 10]. We use a probability parameter $P_1$ to independently determine the existence of a direct constraint between each pair of distinct variables. That is, for each pair of distinct variables $X_i$ and $X_j$, we introduce a direct constraint between them with probability $P_1$. Moreover, we use a probability parameter $P_2$ to determine the compatible tuples of a direct constraint. For a pair of variables $X_i$ and $X_j$ with a direct constraint between them, each tuple $(X_i \leftarrow d_{ip}, X_j \leftarrow d_{jq})$ is independently deemed to be compatible with probability $1 - P_2$. We set $P_1 = 1$ and $P_2 = 0.4$ to obtain representative results for all approaches.

Model A has a tendency to produce mostly unsatisfiable CSP instances with increasing $N$ [18, 20]. Therefore, we use a “hidden solution” method to generate satisfiable CSP instances whenever required. In this method, a set of hidden solutions of the CSP instance is chosen a priori.<fnref n="8" /> A hidden solution $(X_1 \leftarrow d_{1p_1}, X_2 \leftarrow d_{2p_2} \dots X_N \leftarrow d_{Np_N})$ is utilized as follows: While generating the direct constraints using Model A, a direct constraint between variables $X_i$ and $X_j$ reserves the tuple $(X_i \leftarrow d_{ip_i}, X_j \leftarrow d_{jp_j})$ as being compatible before the other tuples are set using the probability parameter $P_2$. Therefore, $(X_1 \leftarrow d_{1p_1}, X_2 \leftarrow d_{2p_2} \dots X_N \leftarrow d_{Np_N})$ satisfies all the direct constraints and, consequently, qualifies as a solution. Similarly, multiple hidden solutions can be utilized with the following modification in the generation procedure: A direct constraint between variables $X_i$ and $X_j$ reserves multiple tuples as being compatible. For generating satisfiable CSP instances, we pick the number of hidden solutions uniformly at random to be an integer within the range [1, 10]. We pick a hidden solution itself by assigning a domain value chosen independently and uniformly at random for each variable from its domain.

We generate three datasets: Dataset-1, Dataset-2, and Dataset-3. For each dataset, we generate 1000 training instances and 1000 test instances. Each training and test set has an equal number of satisfiable and unsatisfiable instances.

In Dataset-1, we generate the instances using Model A. Since Model A frequently generates unsatisfiable instances, we use a complete CSP solver to identify and collect such instances. We generate the satisfiable instances using the hidden solution method, as described above.

In Dataset-2, we generate the satisfiable instances as in Dataset-1. However, we design and generate the unsatisfiable instances to be more challenging. We do this by hiding two complementary pseudo-solutions $(X_1 \leftarrow d_{1p_1}, X_2 \leftarrow d_{2p_2} \dots X_N \leftarrow d_{Np_N})$ and $(X_1 \leftarrow d_{1q_1}, X_2 \leftarrow d_{2q_2} \dots X_N \leftarrow d_{Nq_N})$. We identify a pair of distinct variables $X_i$ and $X_j$ such that $d_{ip_i} \neq d_{iq_i}$ and $d_{jp_j} \neq d_{jq_j}$. All direct constraints between distinct variables $X_s$ and $X_t$ such that $\{X_s, X_t\} \neq \{X_i, X_j\}$ are generated as before by reserving the tuples $(X_s \leftarrow d_{sp_s}, X_t \leftarrow d_{tp_t})$ and $(X_s \leftarrow d_{sq_s}, X_t \leftarrow d_{tq_t})$ as being compatible. However, the direct constraint between $X_i$ and $X_j$ reserves the tuples $(X_i \leftarrow d_{ip_i}, X_j \leftarrow d_{jp_j})$ and $(X_i \leftarrow d_{iq_i}, X_j \leftarrow d_{jp_j})$ as being compatible and reserves the tuples $(X_i \leftarrow d_{ip_i}, X_j \leftarrow d_{jp_j})$ and $(X_i \leftarrow d_{iq_i}, X_j \leftarrow d_{jq_j})$ as being not compatible. We finally use a complete CSP solver to verify that the CSP instance is indeed unsatisfiable.<fnref n="9" />

In Dataset-3, we generate the satisfiable instances as in Dataset-1. However, we design and generate the unsatisfiable instances differently from in Dataset-2. We do this by first hiding two complementary pseudo-solutions $(X_1 \leftarrow d_{1p_1}, X_2 \leftarrow d_{2p_2} \dots X_N \leftarrow d_{Np_N})$ and $(X_1 \leftarrow d_{1q_1}, X_2 \leftarrow d_{2q_2} \dots X_N \leftarrow d_{Nq_N})$, as in Dataset-2. However, we gather all variables $X_{r_1}, X_{r_2}, \dots, X_{r_M}$ for which the two pseudo-solutions have different assignments

<footnote n="8">⁸ The CSP instance can have other solutions as well.</footnote>

<footnote n="9">⁹ This procedure frequently generates unsatisfiable instances, as required. However, satisfiable instances that are generated occasionally are filtered out by the CSP solver.</footnote>
of domain values, that is, $d_{rp_{r_m}} \neq d_{rq_{r_m}}$, for all $1 \le m \le \bar{M}$. For any two distinct variables $X_i$ and $X_j$ in $\{X_{r_1}, X_{r_2}, \dots, X_{r_{\bar{M}}}\}$, we reserve the tuples $(X_i \leftarrow d_{ip_i}, X_j \leftarrow d_{jq_j})$ and $(X_i \leftarrow d_{iq_i}, X_j \leftarrow d_{jp_j})$ as being not compatible. Finally, we pick two distinct variables $X_s$ and $X_t$ from $\{X_{r_1}, X_{r_2}, \dots, X_{r_{\bar{M}}}\}$ and overwrite the tuples $(X_s \leftarrow d_{sp_s}, X_t \leftarrow d_{tq_t})$ and $(X_s \leftarrow d_{sq_s}, X_t \leftarrow d_{tp_t})$ as being compatible and reserve the tuples $(X_s \leftarrow d_{sp_s}, X_t \leftarrow d_{tp_t})$ and $(X_s \leftarrow d_{sq_s}, X_t \leftarrow d_{tq_t})$ as being not compatible. As before, we use a complete CSP solver to verify that the CSP instance is indeed unsatisfiable.

## 5.3 Results

We show three sets of results pertaining to FastMapSVM. First, we show the 2-dimensional and the 3-dimensional embeddings that FastMapSVM produces to aid visualization. Second, we show the behavior of FastMapSVM with respect to the hyperparameter $K$, that is, the number of dimensions and with respect to the size of the training data. Third, we show the comparative performance of FastMapSVM against DGCNN, GIN, and arc-consistency.

FastMapSVM used the SVM classifier from the scikit-learn library. Its hyperparameter settings were determined by grid search. For Dataset-1 and Dataset-3, the hyperparameters were regularization parameter = 8, kernel = "rbf", and kernel coefficient = "scale". For Dataset-2, the hyperparameters were regularization parameter = 8, kernel = "poly", and kernel coefficient = "scale".

Figure 5 shows a perspicuous visualization of the CSP test instances for all three datasets. This visualization capability is unique to FastMapSVM. We note that while the accuracy, recall, precision, and the F1 score of FastMapSVM typically increase with increasing $K$, $K=2$ and $K=3$ are the only two values that support visualization. Still, in most cases, Figure 5 shows a clear separation between the satisfiable and unsatisfiable instances. Moreover, the separation is clearer in the 3-dimensional embeddings compared to their 2-dimensional counterparts.

Figure 6a shows the behavior of FastMapSVM with respect to the number of dimensions $K$ on Dataset-1. Its behavior on the other datasets is similar. The performance metrics, that is, the accuracy, recall, precision, and the F1 score, improve with increasing $K$. This is intuitively expected since the distances between the CSP instances can be embedded with lower distortion in higher dimensions. However, Figure 6a also shows that a point of diminishing returns is attained rather quickly at around $K=8$. This shows that $K=8, 9$, or $10$ is good enough for the CSP domain. Finally, Figure 6a also shows that the improvements in the performance metrics are significant between $K=2$ and $K=8$.

Figure 6b shows the behavior of FastMapSVM with respect to the size of the training data on Dataset-1. Its behavior on the other datasets is similar. The performance metrics improve with increasing size of the training data. Figure 6b also shows that the improvements in the performance metrics are significant between 128 and 256 training data instances. Further improvements are gradual between 256 and 1000 training data instances. This shows that FastMapSVM has the capability to achieve good performance from relatively small amounts of training data and training time.

Table 1 shows a comparison of all the competing methods on all three datasets with respect to all of the performance metrics. It uses $K=8$ for FastMapSVM. It also shows two versions of DGCNN and GIN: the "labeled" version and the "unlabeled" version.

We recollect that in the graphical representation of a binary CSP instance, a vertex represents a domain value and is tagged with the name of the variable that it belongs to. Information in these tags is available to be utilized by the DGCNN and the GIN. The labeled versions of DGCNN and GIN utilize this information while the unlabeled versions ignore this
<figure><img src="image_393.png" /><figcaption>(a) Dataset-1 CSP instances embedded in a 2-dimensional Euclidean space by FastMapSVM</figcaption></figure>

<figure><img src="image_394.png" /><figcaption>(c) Dataset-2 CSP instances embedded in a 2-dimensional Euclidean space by FastMapSVM</figcaption></figure>

<figure><img src="image_395.png" /><figcaption>(e) Dataset-3 CSP instances embedded in a 2-dimensional Euclidean space by FastMapSVM</figcaption></figure>

<figure><img src="image_396.png" /><figcaption>(b) Dataset-1 CSP instances embedded in a 3-dimensional Euclidean space by FastMapSVM</figcaption></figure>

<figure><img src="image_397.png" /><figcaption>(d) Dataset-2 CSP instances embedded in a 3-dimensional Euclidean space by FastMapSVM</figcaption></figure>

<figure><img src="image_398.png" /><figcaption>(f) Dataset-3 CSP instances embedded in a 3-dimensional Euclidean space by FastMapSVM</figcaption></figure>

■ **Figure 5** The figure shows the low-dimensional Euclidean embeddings produced by FastMapSVM for classifying CSP instances. Mostly, there is a clear separation of satisfiable instances (blue) and unsatisfiable instances (red).

information. Table 1 shows that the unlabeled versions perform better than their labeled counterparts. While this is a little surprising, it is likely that the unlabeled versions indirectly perform permutation reasoning on the tags (names of variables) much more efficiently.

Table 1 shows that FastMapSVM generally outperforms all other competing methods by a significant margin. Even on a particular dataset where it is not the top performer with respect to a particular performance metric, it is a close second. Overall, Dataset-1 seems to be the easiest for all methods and Dataset-2 seems to be the hardest for all methods.

In comparison to arc-consistency, FastMapSVM is significantly better on Dataset-2 and Dataset-3. On these datasets, arc-consistency declares all test instances as being “satisfiable”, leading to a perfect recall score but very poor precision, accuracy, and F1 scores. On the one
40:14 FastMapSVM for Predicting CSP Satisfiability

<figure><img src="image_399.png" /><figcaption>(a) Influence of the number of dimensions on the performance of FastMapSVM</figcaption></figure>

<figure><img src="image_400.png" /><figcaption>(b) Influence of the size of the training data on the performance of FastMapSVM</figcaption></figure>

■ **Figure 6** The figure shows the behavior of FastMapSVM with respect to the number of dimensions and with respect to the size of the training data. The performance metrics include the accuracy, recall, precision, and the F1 score.

hand, this shows that arc-consistency is ineffective in recognizing unsatisfiable CSP instances. On the other hand, it also shows that CSP instances generated as in Dataset-1 are insufficient to conclusively evaluate competing ML methods. In contrast, FastMapSVM performs well on all three datasets.

In comparison to DGCNN and GIN, FastMapSVM is significantly better on all three datasets. On the accuracy, recall, and F1 scores, FastMapSVM is better than DGCNN, which in turn is better than GIN. GIN generally has high precision scores but very poor recall scores. This shows that it is poor in identifying satisfiable instances but is mostly correct
<table><caption>**Table 1** The table shows all performance metrics for all competing methods on all datasets. “AC” represents arc-consistency.</caption><thead><tr><th>Dataset</th><th>Model</th><th>Accuracy</th><th>Recall</th><th>Precision</th><th>F1</th></tr></thead><tbody><tr><td rowspan="6">Dataset-1</td><td>FastMapSVM</td><td>96.7%</td><td>97.0%</td><td>96.4%</td><td>96.7%</td></tr><tr><td>AC</td><td><b>99.3%</b></td><td><b>100.0%</b></td><td>96.7%</td><td><b>98.3%</b></td></tr><tr><td>DGCNN (unlabeled)</td><td>94.2%</td><td>92.2%</td><td>96.0%</td><td>94.1%</td></tr><tr><td>DGCNN (labeled)</td><td>82.3%</td><td>82.0%</td><td>82.5%</td><td>82.2%</td></tr><tr><td>GIN (unlabeled)</td><td>84.1%</td><td>67.0%</td><td><b>99.4%</b></td><td>80.0%</td></tr><tr><td>GIN (labeled)</td><td>56.4%</td><td>52.6%</td><td>56.9%</td><td>54.7%</td></tr><tr><td rowspan="6">Dataset-2</td><td>FastMapSVM</td><td><b>82.9%</b></td><td>72.8%</td><td><b>91.2%</b></td><td><b>81.0%</b></td></tr><tr><td>AC</td><td>50.1%</td><td><b>100.0%</b></td><td>50.1%</td><td>66.8%</td></tr><tr><td>DGCNN (unlabeled)</td><td>73.4%</td><td>61.4%</td><td>80.8%</td><td>69.8%</td></tr><tr><td>DGCNN (labeled)</td><td>53.9%</td><td>51.2%</td><td>54.1%</td><td>52.6%</td></tr><tr><td>GIN (unlabeled)</td><td>71.9%</td><td>49.0%</td><td>90.4%</td><td>63.6%</td></tr><tr><td>GIN (labeled)</td><td>54.4%</td><td>51.6%</td><td>54.7%</td><td>53.1%</td></tr><tr><td rowspan="6">Dataset-3</td><td>FastMapSVM</td><td><b>95.4%</b></td><td>94.4%</td><td>96.3%</td><td><b>95.3%</b></td></tr><tr><td>AC</td><td>50.0%</td><td><b>100.0%</b></td><td>50.0%</td><td>66.7%</td></tr><tr><td>DGCNN (unlabeled)</td><td>90.3%</td><td>86.6%</td><td>93.5%</td><td>89.9%</td></tr><tr><td>DGCNN (labeled)</td><td>74.7%</td><td>71.8%</td><td>76.2%</td><td>73.9%</td></tr><tr><td>GIN (unlabeled)</td><td>78.4%</td><td>53.6%</td><td><b>98.5%</b></td><td>69.4%</td></tr><tr><td>GIN (labeled)</td><td>57.4%</td><td>53.8%</td><td>58.0%</td><td>55.8%</td></tr></tbody></table>

when it does so. FastMapSVM does not have this drawback. Moreover, on the accuracy and F1 scores, FastMapSVM outperforms the closest competitor (DGCNN) by larger margins with increasing hardness of the CSP instances, that is, in the order of Dataset-1, Dataset-3, and Dataset-2. Even on the metric of efficiency, FastMapSVM outperforms DGCNN and GIN.<fnref n="10" />

# 6 Conclusions and Future Work

In this paper, we introduced a novel ML framework, called FastMapSVM, for the task of predicting CSP satisfiability. FastMapSVM overcomes the hurdles faced by other ML approaches in the CSP domain. It leverages a distance function on CSPs that is defined via maxflow computations. FastMapSVM is applicable to CSP training and test instances of different sizes and is invariant to both variable-orderings and domain value-orderings. This allows it to bypass the onus of having to learn equivalence classes of CSP instances and, therefore, requires significantly smaller amounts of time and data for model training compared to other ML algorithms. FastMapSVM also uses the intelligence of SVMs, kernel methods, and maxflow computations, accounting for its superior empirical performance, even over state-of-the-art graph neural networks. Moreover, it facilitates a perspicuous visualization of the CSP instances, their distribution, and the classification boundaries between them. Overall, the FastMapSVM framework for CSPs has broader applicability and various representational and combinatorial advantages compared to other ML approaches.

There are many avenues for future work. These include the design of better distance functions on CSPs, the application of FastMapSVM to optimization variants of CSPs, and the general facilitation of integrating constraint reasoning and ML methods via FastMapSVM.

<footnote n="10">¹⁰ DGCNN and GIN ran on a different platform compared to FastMapSVM. However, the ballpark results are still conclusive.</footnote>
References

1 Abder-Rahman Ali, Micael S Couceiro, Aboul Ella Hassanien, and D Jude Hemanth. Fuzzy c-means based on Minkowski distance for liver CT image segmentation. *Intelligent Decision Technologies*, 2016.

2 Alejandro Arbelaez, Youssef Hamadi, and Michele Sebag. Continuous search in constraint programming. In *Proceedings of the 22nd IEEE International Conference on Tools with Artificial Intelligence*, 2010.

3 Liron Cohen, Tansel Uras, Shiva Jahangiri, Aliyah Arunasalam, Sven Koenig, and T. K. Satish Kumar. The FastMap algorithm for shortest path computations. In *Proceedings of the International Joint Conference on Artificial Intelligence*, 2018.

4 Rina Dechter. *Constraint processing*. Morgan Kaufmann, 2003.

5 Christos Faloutsos and King-Ip Lin. FastMap: A fast algorithm for indexing, data-mining and visualization of traditional and multimedia datasets. In *Proceedings of the ACM SIGMOD International Conference on Management of Data*, 1995.

6 Ian Philip Gent, Christopher Anthony Jefferson, Lars Kotthoff, Ian James Miguel, Neil Charles Armour Moore, Peter Nightingale, and Karen Petrie. Learning when to use lazy learning in constraint solving. In *Proceedings of the 19th European Conference on Artificial Intelligence*, 2010.

7 Alessio Guerri and Michela Milano. Learning techniques for automatic algorithm portfolio selection. In *Proceedings of the 16th European Conference on Artificial Intelligence*, 2004.

8 Serdar Kadioglu, Yuri Malitsky, Meinolf Sellmann, and Kevin Tierney. ISAC-instance-specific algorithm configuration. In *Proceedings of the 19th European Conference on Artificial Intelligence*, 2010.

9 Lars Kotthoff. Algorithm selection for combinatorial search problems: A survey. *Data Mining and Constraint Programming: Foundations of a Cross-Disciplinary Approach*, 2016.

10 Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. Imagenet classification with deep convolutional neural networks. *Communications of the ACM*, 2017.

11 Ang Li, Peter Stuckey, Sven Koenig, and T. K. Satish Kumar. A FastMap-based algorithm for block modeling. In *Proceedings of the International Conference on the Integration of Constraint Programming, Artificial Intelligence, and Operations Research*, 2022.

12 Jiaoyang Li, Ariel Felner, Sven Koenig, and T. K. Satish Kumar. Using FastMap to solve graph problems in a Euclidean space. In *Proceedings of the International Conference on Automated Planning and Scheduling*, 2019.

13 Eoin O'Mahony, Emmanuel Hebrard, Alan Holland, Conor Nugent, and Barry O'Sullivan. Using case-based reasoning in an algorithm portfolio for constraint solving. In *Proceedings of the Irish Conference on Artificial Intelligence and Cognitive Science*, 2008.

14 Arti Patle and Deepak Singh Chouhan. SVM kernel functions for classification. In *Proceedings of the International Conference on Advances in Technology and Engineering*, 2013.

15 Luca Pulina and Armando Tacchella. A multi-engine solver for quantified boolean formulas. In *Proceedings of the 13th International Conference on Principles and Practice of Constraint Programming*, 2007.

16 Faisal Rahutomo, Teruaki Kitasuka, and Masayoshi Aritsugi. Semantic cosine similarity. In *Proceedings of the International Student Conference on Advanced Science and Technology*, 2012.

17 V David Sánchez A. Advanced support vector machines and kernel methods. *Neurocomputing*, 2003.

18 Barbara M Smith and Martin E Dyer. Locating the phase transition in binary constraint satisfaction problems. *Artificial Intelligence*, 1996.

19 Malcolm White, Kushal Sharma, Ang Li, T. K. Satish Kumar, and Nori Nakata. FastMapSVM: Classifying complex objects using the FastMap algorithm and support-vector machines. *arXiv preprint arXiv:2204.05112*, 2022.
20 Hong Xu, Sven Koenig, and T. K. Satish Kumar. Towards effective deep learning for constraint satisfaction problems. In *Proceedings of the 24th International Conference on Principles and Practice of Constraint Programming*, 2018.

21 Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka. How powerful are graph neural networks? *arXiv preprint arXiv:1810.00826*, 2018.

22 Lin Xu, Frank Hutter, Holger H Hoos, and Kevin Leyton-Brown. SATzilla: portfolio-based algorithm selection for SAT. *Journal of Artificial Intelligence Research*, 2008.

23 Muhan Zhang, Zhicheng Cui, Marion Neumann, and Yixin Chen. An end-to-end deep learning architecture for graph classification. In *Proceedings of the AAAI Conference on Artificial Intelligence*, 2018.
<figure><img src="image_401.png" /></figure>
Improving Local Search for Pseudo Boolean
Optimization by Fragile Scoring Function and Deep
Optimization

Wenbo Zhou

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

School of Information Science and Technology, Northeast Normal University, Changchun, China
Key Laboratory of Applied Statistics of MOE, Northeast Normal University, Changchun, China
Key Laboratory of Symbolic Computation and Knowledge Engineering of MOE, Jilin University, Changchun, China

Yujiao Zhao

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

School of Information Science and Technology, Northeast Normal University, Changchun, China

Yiyuan Wang<fnref n="1" />

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

School of Information Science and Technology, Northeast Normal University, Changchun, China
Key Laboratory of Applied Statistics of MOE, Northeast Normal University, Changchun, China

Shaowei Cai

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

State Key Laboratory of Computer Science, Institute of Software, Chinese Academy of Sciences,
Beijing, China
School of Computer Science and Technology, University of Chinese Academy of Sciences, Beijing,
China

Shimao Wang

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

School of Information Science and Technology, Northeast Normal University, Changchun, China

Xinyu Wang

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

School of Information Science and Technology, Northeast Normal University, Changchun, China

Minghao Yin<fnref n="1" />

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

School of Information Science and Technology, Northeast Normal University, Changchun, China
Key Laboratory of Applied Statistics of MOE, Northeast Normal University, Changchun, China

Abstract

Pseudo-Boolean optimization (PBO) is usually used to model combinatorial optimization problems,
especially for some real-world applications. Despite its significant importance in both theory and
applications, there are few works on using local search to solve PBO. This paper develops a novel
local search framework for PBO, which has three main ideas. First, we design a two-level selection
strategy to evaluate all candidate variables. Second, we propose a novel deep optimization strategy
to disturb some search spaces. Third, a sampling flipping method is applied to help the algorithm
jump out of local optimum. Experimental results show that the proposed algorithms outperform
three state-of-the-art PBO algorithms on most instances.

2012 ACM Subject Classification Computing methodologies → Search methodologies

Keywords and phrases Local Search, Pseudo-Boolean Optimization, Deep Optimization

Digital Object Identifier 10.4230/LIPIcs.CP.2023.41

Supplementary Material Software (Source Code): https://github.com/yiyuanwang1988/DeepOpt-PBO, archived at swh:1:dir:a1d58cf93325bed3675f5872f42adf459a518a92

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

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

<footnote n="1">¹ corresponding author</footnote>
**Funding** This work is supported by the National Natural Science Foundation of China under Grant No.61806050, the Natural Science Research Foundation of Jilin Province of China under Grant Nos.YDZJ202201ZYTS412 and YDZJ202201ZYTS423, Ministry of Education under Grant No.2021BCF01002, CCF-Huawei Populus Grove Fund, and the Fundamental Research Funds for the Central Universities under Grant Nos.2412023YQ003, 2412022ZD015, 2412022ZD018, 2412022QD040 and 93K172022K10.

1 Introduction

The Boolean Satisfiability (SAT) problem is a prototypical NP-complete problem, whose
aim is to determine whether a given propositional formula is satisfiable or not. The SAT
problem plays a core role in many domains of computer science and artificial intelligence [14].
Many real-world problems can be encoded into SAT and its optimization version MaxSAT
and solved using their powerful solvers. However, due to the limited expressive power of
SAT and MaxSAT, their encodings often generate very large problem instances. For such
cases, pseudo-Boolean optimization (PBO) provides a more expressive and natural way to
express constraints than SAT and MaxSAT [23]. Besides, PBO is very close enough to SAT
to benefit from the recent advances in SAT solving [23].

The PBO consists of a set of pseudo-Boolean constraints and an objective function, whose
goal is to find a solution that minimizes the objective function and satisfies all pseudo-Boolean
constraints. Solvers for PBO can be divided into complete and incomplete solvers. Up to
now, many complete PBO solvers have been proposed, such as Sat4j [3], Open-WBO [21],
NaPS [24], RoundingSat [12, 10], RSCard [13], RS/lp [9] and PBO-IHS [26, 27].

Compared to complete PBO solvers, there are relatively fewer incomplete solvers on PBO. The reason may be that PBO is more complicated than SAT and how to find a suitable variable to flip is still difficult. As one of the most popular incomplete approaches, local search can find an approximate solution within a reasonable time [32, 33]. Lei et al. [18] proposed a novel local search algorithm called LS-PBO to handle PBO. Key features of LS-PBO include a converting method to obtain corresponding objective constraints, a weighting scheme to guide the search direction, and a well-designed scoring function to flip some candidate variables. Besides, for some special cases of PBO such as NK-Landscapes and MAX-kSAT, a new perturbation strategy called VIGbP was proposed [29]. According to the literature, the current best incomplete algorithm for PBO is LS-PBO.

In this work, to further improve the performance of local search algorithms on solving
PBO, we propose a local search framework for PBO based on three main ideas.

First, we present a two-level selection strategy to choose which variable to flip. In our proposed algorithm, we use the previous scoring function *score* [18] as a primary scoring function, which is defined as the decrease of the total penalty of related constraints and objective function. To address the issue about tie-breaking in the primary scoring function, we propose a fragile scoring function *hhscore* as the secondary scoring function. The proposed *hhscore* can greatly differentiate between two kinds of satisfied constraints by using the definition of satisfied threshold.

Second, we propose a novel deep optimization strategy (DeepOpt) to deeply probe some
regions using locked and unlocked operations during the local search. Recently, the DeepOpt
strategy was first proposed by Chen et al. [8] and has been successfully applied in solving
dominating set problems. In our proposed DeepOpt strategy, we preferentially select some
variables in unsatisfied constraints as unlocked operations. Moreover, we use some trigger
conditions to decide whether the algorithm calls DeepOpt or not.
Third, we design a sampling flipping method to modify a current candidate solution when the algorithm tramps into local optimum. In the proposed method, we adopt two sampling ways to collect some samples among all unsatisfied constraints. We further employ a probabilistic heuristic “best from multiple selections” (BMS) [6] to select a candidate variable. Besides, we also apply a new scoring function to simultaneously flip two variables.

By incorporating these three ideas and some other tricks, we develop two local search algorithms for PBO. Extensive experiments are carried out to evaluate our algorithms on the benchmarks used in the literature. Experimental results show that the proposed algorithms outperform three state-of-the-art PBO algorithms on almost all the benchmarks.

## 2 Preliminaries

Since a non-linear pseudo-Boolean (PB) constraint can be translated into an equivalent set of linear PB constraints [23], we only start with a review of the basics of linear PB constraints here. A linear PB constraint is defined over a finite set of Boolean variables. Boolean variable $x_i$ can take only two values false (0) and true (1). A literal $l_i$ over a Boolean variable $x_i$ is either $x_i$ or $\bar{x}_i = 1 - x_i$. A linear PB constraint is a 0-1 integer inequality.

$$ \sum_{i} a_i l_i \triangleright b \qquad (1) $$

where $a_i$ and $b$ are integer constants, $l_i$ are literals and $\triangleright \in \{=, >, \ge, <, \le\}$ is one of the classical relational operators. All PB constraints can be normalized into the following form.

$$ \sum_{i} a_i l_i \geq b \qquad (2) $$

where all the literals $l_i$ are distinct, and all the coefficients $a_i$ and the degree $b$ are non-negative integers.

A PB formula is a conjunction of PB constraints, denoted as $F = C_1 \land C_2 \land \dots \land C_m$, where $C_p$ ($p \in \mathbb{Z}$, $1 \le p \le m$) is a PB constraint. The PBO problem consists of a PB formula $F$ and an objective function $O: \sum_i c_i l_i$ where $c_i$ is a non-negative integer coefficient. Given a PB constraint $C_p: \sum_i a_i^p l_i^p \ge b^p$, the sum of its coefficients is defined as $sum(C_p) = \sum_i a_i^p$, its average coefficient is defined as $coeff(C_p) = sum(C_p)/|L(C_p)|$ where $L(C_p)$ is the set of literals in $C_p$, and its maximum coefficient $a_{max}^p$ is the maximum value of coefficients in $C_p$. The average coefficient of an objective function $O$ is defined as $coeff(O) = \sum_i c_i/|L(O)|$ where $L(O)$ is the set of literals in $O$.

Given a PB formula $F$, its complete assignment is a mapping that assigns 0 or 1 to each variable. Given a complete assignment of $F$, if a literal evaluates to true, we say it is a true literal and otherwise it is a false literal. A PB constraint $C_p$ is satisfied when the left and right terms of the constraint evaluate to integers which satisfy the relational operator. Otherwise, $C_p$ is unsatisfied. The sum of coefficients of true literals in $C_p$ is denoted as $SatL(C_p) = \sum_{i=1}^{p} \sum_{l_i \in C_p} a_i^p$. An assignment $\alpha$ of $F$ is feasible if and only if $\alpha$ satisfies all PB constraints in $F$. The value of the objective function of a feasible solution $\alpha$ is denoted as $obj(\alpha)$. The PBO problem aims to obtain a feasible solution for $F$ with the minimum objective value.

### 2.1 Review for Weighting and Scoring Function

Constraint weighting techniques have usually been used to guide and diversify the search process [31, 5, 15]. The weighting based scoring function score for PBO is recently proposed by Lei et al. [18]. Each PB constraint $C_p \in F$ and an objective function $O$ have the property
41:4 Improving Local Search for Pseudo Boolean Optimization

of weighting, denoted as $w(C_p)$ and $w(O)$, respectively. Before introducing the weighting, we
first present a basic concept. Given a PB formula $F$ and $m$ is the number of constraints, the
average constraint coefficient of $F$ is defined as $\text{avg\_coeff} = \sum_{p=1}^{m} w(C_p) \times \text{coeff}(C_p)/m$.
The property of weighting works as follows.

**Weighting Rule 1:** At first, $w(O) = 1$ and $w(C_p) = 1$ for an objective function $O$ and each PB constraint $C_p$.

**Weighting Rule 2:** For each unsatisfied constraint $C_p$, $w(C_p) = w(C_p) + 1$.

**Weighting Rule 3:** If the current $obj(\alpha)$ of the objective function is better than the current best-found value of objective function during the search process so far and $w(O) \times coeff(O) - avg\_coeff \le \zeta$, then $w(O) = w(O) + 1$. In our work, we use the same parameter value of $\zeta$ as [18], i.e., $\zeta = 100$.

Given an assignment $\alpha$ of $F$, if a PB constraint $C_p: \sum_i a_i^p l_i^p \ge b^p$ is unsatisfied, the
penalty of $C_p$ is defined as $w(C_p) \times (b^p - \sum_i a_i^p l_i^p)$. For the objective function $O: \sum_i c_i l_i$, the penalty of $O$ is defined as $w(O) \times \sum_i c_i l_i$. Based on the definition of penalty, we introduce two scoring functions including hard score $hscore$ and objective score $oscore$ as below. For a Boolean variable $x_i$, the respective $hscore(x_i)$ and $oscore(x_i)$ are the decrease of the total penalty of unsatisfied PB constraints and the objective function caused by flipping $x_i$. Combing the above scoring functions, we define the score of a variable $x_i$ as $score(x_i) = hscore(x_i) + oscore(x_i)$. Remark that after flipping some variables during the search process, the corresponding score values should be updated accordingly.

3 Two-Level Selection Strategy

In this section, we introduce a secondary scoring function to reinforce local search algorithms
for PBO and then propose a two-level selection strategy to decide a candidate variable.

**3.1 Fragile Scoring Function**

As a core guidance for the search process, scoring functions play an important role in local
search algorithms, which measure the benefits of a candidate variable. In local search
algorithms for PBO, such benefits *hscore* can be set as the distance between the sum of
coefficients of true literals and the corresponding degree, whereas *oscore* can be set as the
distance between the current and best-found objective values [18].

In our algorithm, we consider the sum of *hscore* and *oscore* (i.e., *score*) as the primary
scoring function, which is the same method as LS-PBO [18]. According to our preliminary
experiments, 9% candidate variables on average have the same largest *score* value during
the search. To address the issue about tie-breaking in the primary scoring function, previous
work uses the *age* information of variables as the secondary scoring function, where *age* is
defined as the number of steps since the last time it is flipped. But the experimental results
show that the use of only *age* cannot effectively guide the search process. Thus, to further
choose a variable among these variables with the same best *score* value, we design a novel
fragile scoring function denoted as *hhscore*.

Before introducing *hhscore*, we first introduce a necessary concept. For a PB constraint
$C_p$, $gap(C_p) = \min\{b^p + a_{max}^p, sum(C_p)\}$ is used to denote the *satisfied threshold* of $C_p$, which plays a key role in our proposed *hhscore*. Based on the satisfied threshold of PB constraints, we define a fragile satisfied PB constraint in the following.

► **Definition 1.** For a satisfied PB constraint $C_p$ (i.e., $\text{SatL}(C_p) \ge b^p$), the $C_p$ is a fragile satisfied PB constraint if and only if $\text{SatL}(C_p) < gap(C_p)$.
On the one hand, if a satisfied PB constraint $C_p$ is fragile, we think that flipping any variable in this PB constraint would probably make this constraint become unsatisfied. On the other hand, if a satisfied PB constraint is not fragile (i.e., $SatL(C_p) \ge gap(C_p)$), we think this PB constraint is solid, which means that flipping any true literal in $C_p$ would not make $C_p$ become unsatisfied with a high probability. Although we have also tried a similar property that measures the number of true literals in a satisfied PB constraint $C_p$, such as $gap(C_p) = b^p + 1$, we did not find it useful in our algorithm.

To maintain the information of “fragile” for each literal during the search process, we define inner as below.

**► Definition 2.** Suppose that a Boolean variable $x_i$ whose literal $l_i$ appears in some PB constraints (e.g., $C_p$) and the coefficient of $l_i$ in $C_p$ is $a_i^p$. The value of inner($x_i, C_p$) is calculated as follows.

(a) $C_p$ is an unsatisfied PB constraint, i.e., $SatL(C_p) < b^p$:

- If $l_i$ is a true literal in $C_p$, inner($x_i, C_p$) = 0;

- If $l_i$ is a false literal in $C_p$, inner($x_i, C_p$) = $max\{a_i^p - (b^p - SatL(C_p)), 0\}$.

(b) $C_p$ is a fragile satisfied PB constraint, i.e., $b^p \le SatL(C_p) < gap(C_p)$:

- If $l_i$ is a true literal in $C_p$, inner($x_i, C_p$) = $-min\{a_i^p, SatL(C_p) - b^p\}$;

- If $l_i$ is a false literal in $C_p$, inner($x_i, C_p$) = $min\{a_i^p, gap(C_p) - SatL(C_p)\}$.

(c) $C_p$ is satisfied but not a fragile PB constraint, i.e., $gap(C_p) \le SatL(C_p)$:

- If $l_i$ is a true literal in $C_p$, inner($x_i, C_p$) = $-max\{a_i^p - (SatL(C_p) - gap(C_p)), 0\}$;

- If $l_i$ is a false literal in $C_p$, inner($x_i, C_p$) = 0.

The value of inner($x_i, C_p$) is subject to three distinct factors including the state of $C_p$ (i.e., satisfied or unsatisfied), the value of $l_i$ (i.e., true or false), and the coefficient of $l_i$ (i.e., $a_i^p$). To make the readers easily understand the above Definition 2, we list all the situations corresponding to different inner values in Figure 1.

Considering all the PB constraints in which variable $x$'s literal appears, the proposed fragile scoring function hhscore of variable $x_i$ is defined as below.

$$hhscore(x_i) = \sum_{p=1}^{n_x} inner(x_i, C_p) \qquad (3)$$

where $n_x$ is the number of PB constraints including the literal of $x_i$.

## 3.2 Selection Rule

Combining the respective advantages of score and hhscore, we propose a two-level selection strategy as follows.

**Flipping Rule.** Flip a variable $x_i$ with the biggest score($x_i$) value, breaking ties by preferring the one with the biggest hhcore($x_i$) value, further ties are broken randomly.

The proposed secondary scoring function is inspired by subscore [7], but has two essential differences. First, our proposed scoring function can be considered as a general version of subscore because the value of satisfied threshold gap is equal to 2 for the SAT problem, which is the same function as subscore. Second, previous work uses a linear combination of subscore and another scoring function as the primary scoring function, whereas our work considers a two-level scoring function to guide the search process. In addition, experiments show that the trigger fraction of using hhscore at least once for all tested instances is about 99.6%.
<table><tr><td>(a) $SatL(C_p) < b^p$</td></tr><tr><td>$C_p$ unsatisfied</td></tr><tr><td>$l_i$ true</td></tr><tr><td>0</td></tr></table>

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

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

<table><tr><td>(a) $SatL(C_p) < b^p$</td></tr><tr><td>$C_p$ unsatisfied</td></tr><tr><td>$l_i$ false</td></tr><tr><td>$\max(a_i^p - (b^p - SatL(C_p)), 0)$</td></tr></table>

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

$$ (a1) \quad a_i^p - (b_p - SatL(C_p)) > 0 $$

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

$$ (a2) \quad a_i^p - (b_p - SatL(C_p)) \le 0 $$

<table><tr><td>(b) $b^p \le SatL(C_p) < gap(C_p)$</td></tr><tr><td>$C_p$ satisfied</td></tr><tr><td>$l_i$ true</td></tr><tr><td>$-\min(a_i^p, SatL(C_p) - b^p)$</td></tr></table>

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

$$ (b1) \quad a_i^p \le SatL(C_p) - b_p $$

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

$$ (b2) \quad a_i^p > SatL(C_p) - b_p $$

<table><tr><td>(b) $b^p \le SatL(C_p) < gap(C_p)$</td></tr><tr><td>$C_p$ satisfied</td></tr><tr><td>$l_i$ false</td></tr><tr><td>$\min(a_i^p, gap(C_p) - SatL(C_p))$</td></tr></table>

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

$$ (b3) \quad a_i^p \le gap(C_p) - SatL(C_p) $$

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

$$ (b4) \quad a_i^p > gap(C_p) - SatL(C_p) $$

<table><tr><td>(c) $gap(C_p) \le SatL(C_p)$</td></tr><tr><td>$C_p$ satisfied</td></tr><tr><td>$l_i$ true</td></tr><tr><td>$-\min(a_i^p - (SatL(C_p) - gap(C_p)), 0)$</td></tr></table>

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

$$ (c1) \quad a_i^p - (SatL(C_p) - gap(C_p)) > 0 $$

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

$$ (c2) \quad a_i^p - (SatL(C_p) - gap(C_p)) \le 0 $$

<table><tr><td>(c) $gap(C_p) \le SatL(C_p)$</td></tr><tr><td>$C_p$ satisfied</td></tr><tr><td>$l_i$ false</td></tr><tr><td>0</td></tr></table>

<figure><img src="image_421.png" /><figcaption>Figure 1 A graphical explanation of inner.</figcaption></figure>
■ **Algorithm 1** DeepOpt.

Input: PBO instance *F* and a perturbation initialization assignment *α*
Output: A perturbation solution *α* of *F*
1 *UnSatSet* := *unsat*(*α*) and *SatSet* := *F* \ *UnsatSet*;
2 *CandSet* := ∅;
3 **while** |*CandSet*| ≤ *γ* × *n* **do**
4 &nbsp;&nbsp;&nbsp;&nbsp;**if** *UnSatSet* ≠ ∅ **then**
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;select a random unsatisfied PB constraint *C*;
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*UnSatSet* := *UnSatSet* \ {*C*};
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** each variable *x* ∈ *C* **do**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** *x*'s literal is true &amp;&amp; with 50% probability **then** *α* := *α* with *x* flipped;
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*CandSet* := *CandSet* ∪ {*x*};
10 &nbsp;&nbsp;&nbsp;&nbsp;**else**
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;select a random satisfied PB constraint *C*;
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*SatSet* := *SatSet* \ {*C*};
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** each variable *x* ∈ *C* **do**
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** |*unsat*(*α*)| > *MaxHard* **then** **break**;
15 &nbsp;&nbsp;&nbsp;&nbsp;**for** step = 0; step < *L*<sub>opt</sub>; step++ **do**
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;select a variable *x* among |*CandSet*|/2 samples from *CandSet* based on **Flipping Rule**;
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*α* := *α* with *x* flipped;
18 &nbsp;&nbsp;&nbsp;&nbsp;**return** *α*;

4 Deep Optimization for PBO

Local search algorithms usually search the entire space and focus on exploring some promising
spaces using several heuristic strategies. Recently, a general perturbation mechanism called
deep optimization has been proposed by Chen et al. [8], which can deeply probe some regions
based on locked and unlocked operations and then converge to a new solution quickly. Note
that deep optimization is somewhat similar to some classic search frameworks such as large
neighborhood search [25]. According to this general framework, we propose a new deep
optimization approach DeepOpt for PBO.

The pseudo-code of our proposed DeepOpt is reported in Algorithm 1 and the corresponding trigger conditions of DeepOpt will be displayed in the next section. We use the unsat(α) function to denote the set of unsatisfied PB constraints under an assignment α. At first, the algorithm uses UnSatSet and SatSet to store unsatisfied and satisfied PB constraints under a perturbation initialization assignment, respectively (Line 1). Afterward, a candidate variable set CandSet is initialized to an empty set (Line 2), which stores all unlocked variables in the following search. Note that, in our work, we use CandSet to store all unlocked variables which means that these variables can be flipped in the following search phase, while the remaining variables can be seen as locked variables.

The DeepOpt usually consists of two phases including a selection phase (Lines 3–15) to
generate several candidate local search spaces and a search phase (Lines 16–18) to repair
these search spaces.

In the first phase, we flip several variables to achieve the purpose of early preparation.
There are two exit conditions in this phase. The first one is that |CandSet| is larger than
γ × n where n is the number of variables (Line 3), whereas the second one is that the number
of unsatisfied constraints under the current assignment is larger than MaxHard (Line 15)
where $\gamma$ and $\text{MaxHard}$ are two parameters in DeepOpt. At each iteration, if $\text{UnSatSet}$ is not empty, the algorithm preferentially chooses a random unsatisfied PB constraint $C$ from $\text{UnSatSet}$ (Lines 4–5) and puts all variables of $C$ into $\text{CandSet}$ (Line 9). For each variable $x \in C$, if $x$'s literal in $C$ is true, it occurs with 50% probability to flip $x$ (Line 8). If the algorithm selects a random satisfied constraint $C$, the algorithm only adds all variables of $C$ into $\text{CandSet}$ (Line 14). Note that according to our preliminary experiments, if the first phase only selects the candidate variables and does not flip these variables, then the performance of DeepOpt will become bad.

In the second phase, the algorithm applies the search process to perturb the assignment $\alpha$ until the limit of iterations $L_{opt}$ is reached (Line 16). The algorithm employs the BMS strategy [6] in the process of selecting a candidate variable, i.e., randomly choosing $|\text{CandSet}|/2$ variables to compose a candidate set. Afterward, the algorithm flips a variable $x$ among a candidate set based on the flipping rule (Lines 17–18). At last, the algorithm returns the final assignment $\alpha$ (Line 19).

## 5 DeepOpt-PBO Algorithm

In this section, we propose an effective local search algorithm DeepOpt-PBO for PBO, whose main framework is presented in Algorithm 2. The DeepOpt-PBO is mainly divided into the initialization and search phases.

In the initialization phase (Lines 1–3), the best solution $\alpha^*$ is set to $\emptyset$ and its objective value $obj^*$ is set to $+\infty$. To obtain an initial assignment $\alpha$, all variables are set to 0. The algorithm initializes the related weight information according to the weighting rule 1.

In the following, there is an outer loop (Lines 4–32) and an inner loop (Lines 7–31). During the search, whenever a better feasible assignment is obtained, $\alpha^*$ and $obj^*$ are updated accordingly (Line 12). After each inner loop, the algorithm uses the `RestartVar` function to restart a current assignment (Line 32), which will be presented in Section 5.2. Finally, the algorithm returns $\alpha^*$ and $obj^*$ when reaching a time limit.

In each inner loop (`step < L`), the algorithm searches for a local optimal assignment $\alpha$. The algorithm uses `GoodSet` to store variables whose score is larger than 0 (Line 15). If `GoodSet` is not empty, the algorithm flips a candidate variable based on the proposed flipping rule (Lines 16–18). Otherwise, it means that the algorithm tramps into local optimum. The algorithm uses the modified weighting rule 2 and weighting rule 3 to update the corresponding weight value (Lines 20–21). Specifically, we optimize previous weighting rule 2, resulting in a novel modified weighting rule 2, which will be mentioned in Section 5.2. Moreover, the algorithm will use a sampling technique `SampleFlip` to flip one or more variables to help the algorithm itself jump out of local optimum, which will be introduced in the next subsection (Line 22).

**Trigger Conditions of DeepOpt.** In the below part, we will introduce some trigger conditions of DeepOpt in our proposed algorithm. At first, three variables are defined as below.

1) `unhard` is used to denote the number of unsatisfied PB constraints. Before each inner loop, `unhard` is initialized to the number of unsatisfied constraints under the current assignment (i.e., $|\text{unsat}(\alpha)|$) (Line 6). In the inner loop, whenever $|\text{unsat}(\alpha)| < \text{unhard}$, `unhard` is updated accordingly (Lines 8–9). After calling the DeepOpt method, `unhard` will be updated by $|\text{unsat}(\alpha)|$ (Line 31).
Algorithm 2 DeepOpt-PBO.

Input : PBO instance *F* and cutoff time *cutoff*
Output : An assignment *α*\* of *F* and its objective value
1 &alpha;\* := /0/ and *obj*\* := +&infin;;
2 &alpha; := all variables are set to 0;
3 initialize the weight value of the objective function and each constraint to 1;
4 **while** elapsed time < cutoff **do**
5 &nbsp;&nbsp;&nbsp;&nbsp;step<sub>opt</sub> := 1 and *coff* := 1;
6 &nbsp;&nbsp;&nbsp;&nbsp;unhard := |*unsat*(*alpha*)|;
7 &nbsp;&nbsp;&nbsp;&nbsp;**for** step = 1; step < *L*; step++ **do**
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** |*unsat*(*alpha*)| < unhard **then**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;unhard := |*unsat*(*alpha*)|;
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;step<sub>opt</sub> := step<sub>opt</sub>/2 and step := step/2;
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** &alpha; is feasible &amp;&amp; *obj*(*alpha*) < *obj*\* **then**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&alpha;\* := &alpha; and *obj*\* := *obj*(*alpha*);
13 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;step<sub>opt</sub> := step := 1;
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** *coff* &ne; 1 **then** *coff* := *coff*/2;
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;GoodSet := {x | *score*(*x*) > 0};
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** GoodSet &ne; /0/ **then**
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;select a variable *x* in GoodSet based on Flipping Rule;
18 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&alpha; := &alpha; with *x* flipped;
19 &nbsp;&nbsp;&nbsp;&nbsp;**else**
20 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;update the weight of each constraint based on Modified Weighting Rule 2;
21 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;update the weight of the objective function based on Weighting Rule 3;
22 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;SampleFlip(*F*, &alpha;);
23 &nbsp;&nbsp;&nbsp;&nbsp;step<sub>opt</sub> := step<sub>opt</sub> + 1;
24 &nbsp;&nbsp;&nbsp;&nbsp;**if** step<sub>opt</sub>%(coff &times; MinStep) == 0 **then**
25 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** |*unsat*(*alpha*)| &le; MinHard &amp;&amp; with 50% probability **then**
26 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** *coff* &ne; &delta; **then** *coff* := *coff* &times; 2;
27 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;DeepOpt(*F*, &alpha;);
28 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** &alpha;\* &ne; /0/ **then**
29 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&alpha; := &alpha;\*;
30 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;DeepOpt(*F*, &alpha;);
31 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;unhard := |*unsat*(*alpha*)| and step<sub>opt</sub> := 1;
32 &nbsp;&nbsp;&nbsp;&nbsp;RestartVar(*F*, &alpha;);
33 &nbsp;&nbsp;&nbsp;&nbsp;**return** (*alpha*\*, *obj*\*);

2) step<sub>opt</sub> records the non-improvement steps after the initialization phase or the last DeepOpt operation. Before each inner loop, the algorithm sets step<sub>opt</sub> to 1 (Line 5). In each iteration of the inner loop, step<sub>opt</sub> is increased by 1 (Line 23). When the algorithm obtains a better assignment (Line 13) or the algorithm calls the DeepOpt method (Line 27 or 30), step<sub>opt</sub> will be set to 1 (Line 31). If |unsat(*α*)| < unhard, then step<sub>opt</sub> will be cut in half (Lines 8 and 10).

3) *coff* is used to control the frequency of using the DeepOpt method. Before each inner loop, *coff* is initialized to 1 (Line 5). When the algorithm finds a better assignment, the value of *coff* is divided by 2 to reduce the frequency of perturbations (Line 14). If |*unsat*(*α*)| is smaller than parameter *MinHard*, it means unsatisfied PB constraints are few enough to consider perturbations, avoiding tramping into a local optimum. The value of *coff* will be doubled with a 50% probability (Line 26). During this process, parameter *δ* controls the maximum value of *coff*.
41:10 Improving Local Search for Pseudo Boolean Optimization

The algorithm judges whether to call the DeepOpt method under each (coff × MinStep) iteration (Line 24). If |unsat(α)| ≤ MinHard, the algorithm calls the DeepOpt method with a 50% probability (Line 27). Otherwise, if α* is not empty, the algorithm will use α* as a perturbation initialization assignment and then employ the DeepOpt method (Lines 28–30).

■ **Algorithm 3** SampleFlip.

Input : PBO instance *F* and an assignment α
Output : A modified assignment α
1 **if** |*unsat*(*α*)| == 0 **then**
2 &nbsp;&nbsp;&nbsp;select a random variable *x* with *oscore*(*x*) > 0;
3 &nbsp;&nbsp;&nbsp;*α* := *α* with *x* flipped;
4 &nbsp;&nbsp;&nbsp;**return** *α*;
5 select a random unsatisfied constraint *C*;
6 *CSet* := *VSet* := ∅;
7 **if** |*unsat*(*α*)| ≥ *β* **then**
8 &nbsp;&nbsp;&nbsp;**for** *i* = 1 to *β* **do**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;select a random unsatisfied constraint *C*<sub>*i*</sub>;
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*CSet* := *CSet* ∪ {*C*<sub>*i*</sub>};
11 **if** |*L*(*C*)| == 1 **then**
12 &nbsp;&nbsp;&nbsp;**select** only variable *x* in *C* and *α* := *α* with *x* flipped;
13 **else if** |*L*(*C*)| == 2 **then**
14 &nbsp;&nbsp;&nbsp;/* Two variables *x*₁ and *x*₂ in *C* */
15 &nbsp;&nbsp;&nbsp;select a variable *x*<sub>*i*</sub> with the largest *score*<sub>*i*</sub>(*x*₁, *x*<sub>*i*</sub>) value, breaking ties randomly;
16 &nbsp;&nbsp;&nbsp;select a variable *x*<sub>*j*</sub> with the largest *score*<sub>*j*</sub>(*x*₂, *x*<sub>*j*</sub>) value, breaking ties randomly;
17 &nbsp;&nbsp;&nbsp;**if** *score*<sub>*i*</sub>(*x*₁, *x*<sub>*i*</sub>) > 0 || *score*<sub>*j*</sub>(*x*₂, *x*<sub>*j*</sub>) > 0 **then**
18 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** *score*<sub>*i*</sub>(*x*₁, *x*<sub>*i*</sub>) > *score*<sub>*i*</sub>(*x*₂, *x*<sub>*i*</sub>) **then**
19 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**α** := *α* with *x*₁ and *x*<sub>*i*</sub> flipped;
20 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else**
21 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**α** := *α* with *x*₂ and *x*<sub>*j*</sub> flipped;
22 **else**
23 &nbsp;&nbsp;&nbsp;**if** |*unsat*(*α*)| ≥ *β* **then**
24 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** each constraint *C*<sub>*p*</sub> ∈ *CSet* **do**
25 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**select** |*L*(*C*<sub>*p*</sub>)|/2 samples from *L*(*C*<sub>*p*</sub>) and put them into *VSet*;
26 &nbsp;&nbsp;&nbsp;**else**
27 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**select** |*L*(*C*)|/2 samples from *L*(*C*) and put them into *VSet*;
28 **select** a variable *x* from *VSet* based on **Flipping Rule**;
29 **α** := *α* with *x* flipped;
30 **return** *α*;

5.1 Sampling Flipping

The pseudo-code of *SampleFlip* is outlined in Algorithm 3. First, we introduce a new scoring function, which is used in our *SampleFlip* method. When flipping two variables *x*<sub>s</sub> and *x*<sub>z</sub> simultaneously, *score*<sub>t</sub>(*x*<sub>s</sub>, *x*<sub>z</sub>) is defined as the sum of the decrease of the total penalty of unsatisfied PB constraints and the objective function. Although it is easy to see the computation complexity of *score* is quite lower than *score*<sub>t</sub>, *score*<sub>t</sub> can find a better
flipping operation compared to *score*. The method of selecting more candidate elements simultaneously has also been used in some different NP-hard problems [28, 34]. Moreover, in the *SampleFlip* method, we only consider *score*<sub>*t*</sub> in a special case that the number of literals in the selected constraint *C* is 2 (i.e., |*L*(*C*)| = 2).

In the beginning, if there are no unsatisfied constraints, the algorithm selects and flips a random variable *x* with *oscore*(*x*) > 0 (Lines 1–4). There are two main important components, including clause sampling (Lines 5–10) and variable flipping (Lines 11–29) in the algorithm. In the phase of clause sampling, the algorithm adopts two kinds of sampling ways. The first sampling method is to select a random unsatisfied constraint *C* (Line 5). If |*unsat*(*α*)| > *β*, then the algorithm activates the second sampling method, i.e., adding *β* unsatisfied constraints into *CSet* (Lines 7–10).

In the phase of variable flipping, if *C* is a unit clause (i.e., |*L*(*C*)| = 1), the algorithm flips the only variable *x* in *C* (Lines 11–12). If *C* is a binary clause (i.e., *L*(*C*) = {*x*₁, *x*₂}), the algorithm tries to find two pairs of variables {*x*₁, *x*ᵢ} and {*x*₂, *x*ⱼ} with the largest *score*ᵢ value (Lines 14–15). If there exists a positive flipping operation among these two pairs, the algorithm flips a pair with the better *score*ᵢ value (Lines 16–19). Otherwise, the algorithm still flips only one variable among *x*₁ and *x*₂ based on the flipping rule (Lines 20–21). In the subsequent process, the algorithm will depend on the value of parameter *β* to collect some samples from *CSet* or *C* into *VSet* (Lines 23–27). At last, the algorithm selects and then flips the best variable *x* from *VSet* (Lines 28–29).

The intuitive explanations behind *SampleFlip* come from two aspects. First, single flipping mechanism for local search is easy to fall into a local optimum, while sampling strategy can explore the solution space more effectively in a look-ahead way. Second, the sampling strategy selects variables from various unsatisfied constraints, providing more search directions. In addition, we specially focus on the case of two literals, to keep running time low but a good performance.

## 5.2 Some Other Techniques

Based on the main part as shown above, we also introduce some additional heuristics to further improve the efficiency, including a weighting rule and a restart strategy.

**Modified Weighting Rule 2.** The weighting scheme plays an important role in the search process. By increasing the weight of unsatisfied constraints, we can accurately guide the search process toward more efficient directions. In each iteration of the inner loop, if *GoodSet* is empty, the visited times of each unsatisfied constraint will be increased by 1. For an unsatisfied constraint *C*<sub>*p*</sub>, when its visited times are larger than *b*<sup>*p*</sup>/*coeff*(*C*<sub>*p*</sub>), the value of *w*(*C*<sub>*p*</sub>) will be increased by 1, and the value of its visited times is reset to 0.

**Restart Strategy.** The second adopted technique is the restart strategy. Under a current assignment *α*, we can change the original value of each variable with a certain probability, i.e., from 1 (0) to 0 (1). Besides, the restart strategy will reset the weight value of each constraint according to the weighting rule 1. If the number of literals in an objective function *O* is larger than *β* (i.e., |*L*(*O*)| > *β*), then *w*(*O*) will be reset based on the weighting rule 1, too.

# 6 Experimental Evaluation

We first introduce the seven selected benchmarks, three state-of-the-art PBO competitors, and the adopted experimental setup. Then, we carry out extensive experiments to evaluate the performance of our proposed algorithm.
## 6.1 Experiment Preliminaries

Since the literature on local search algorithms for handling PBO is very sparse, we selected all used instances from [18, 10]. To be specific, we considered 3738 instances obtained from three application benchmarks and four standard benchmarks: (1) 24 instances from the minimum-width confidence band problem (MWCB) [2]. These MWCB instances were obtained based on the MIT-BIH arrhythmia database<fnref n="2" />. (2) 18 instances from the wireless sensor network optimization problem (WSNO) [16, 17]. We used the same encoding as previous work [18] to obtain an optimization version of WSNO. (3) 21 instances from the seating arrangements problem (SAP) [1]. These instances were originally proposed in the MaxSAT Evaluation 2017. (4) 1600 OPT-SMALL-INT instances from the most recent PB Competition in 2016 (PB2016)<fnref n="3" />. The PB2016 benchmark is often considered as the main target for comparing against some other PB solvers [10, 27]. (5) The 0-1 integer linear programming optimization benchmark (MIPLIB) contains 267 instances from the mixed integer programming library MIPLIB 2017⁴. (6) 1025 crafted combinatorial instances (CRAFT) are provided in the literature [30]. (7) The Knapsack benchmark (KNAP) consists of a total of 783 instances [22].

<table><caption>Table 1 Tuned parameters of our proposed algorithms.</caption><thead><tr><th>Parameter</th><th>Range</th><th>Final value</th></tr></thead><tbody><tr><td colspan="3"><strong>DeepOpt</strong></td></tr><tr><td>MinStep</td><td>{10<sup>3</sup>,10<sup>4</sup>,10<sup>5</sup>,10<sup>6</sup>}</td><td>10<sup>5</sup></td></tr><tr><td>δ</td><td>{64,128,256}</td><td>128</td></tr><tr><td>γ</td><td>{0.02,0.05,0.08,0.11}</td><td>0.05</td></tr><tr><td>MaxHard</td><td>{30,50,70,90}</td><td>50</td></tr><tr><td>MinHard</td><td>{5,10,15}</td><td>10</td></tr><tr><td>L<sub>opt</sub></td><td>{10,30,50,70}</td><td>50</td></tr><tr><td colspan="3"><strong>Some other parameters in our proposed algorithms</strong></td></tr><tr><td>β</td><td>{50,100,150,200}</td><td>100</td></tr><tr><td>L</td><td>{10<sup>2</sup>,10<sup>3</sup>,10<sup>4</sup>,10<sup>5</sup>}</td><td>10<sup>4</sup></td></tr></tbody></table>

According to the frequency of using the proposed restart strategy, we propose two versions of DeepOpt-PBO, resulting in DeepOpt-PBO-v1 and DeepOpt-PBO-v2. In detail, DeepOpt-PBO-v1 does not use the restart strategy (i.e., parameter *L* is set to INT_MAX). DeepOpt-PBO-v2 is run for half of the computation time given to the algorithm with the restart strategy, while the second half of the available computation time of DeepOpt-PBO-v2 is given to the algorithm without the restart strategy. We also attempted to use the restart strategy during the whole time, and the result is not satisfactory. According to our preliminary experiments by using the automatic configuration tool irace [20], Table 1 shows the parameter values. Specifically, since these benchmarks have different scales, we built a training set and randomly selected 10 instances from the corresponding tested benchmark. The tuning process is given a budget of 5000 runs for the training set with a time budget of 3600s per run.

The proposed algorithms are compared against three state-of-the-art PBO algorithms, including a local search algorithm LS-PBO [18] and two exact PBO solvers, i.e., RoundingSat [10] and PBO-IHS [27]. Although mixed integer programming solvers (e.g., SCIP [4]) and

<footnote n="2">² http://physionet.org/physiobank/database/mitdb/</footnote>

<footnote n="3">³ http://www.cril.univ-artois.fr/PB16/</footnote>

<footnote n="4">⁴ http://miplib.zib.de</footnote>