11:8 Simplifying Step-Wise Explanation Sequences

nested explanations. These clarify a single explanation step ($\mathcal{E}, S, \mathcal{N}$) by creating a new explanation sequence from $S$ with $\mathcal{E}$ and the negation of $\mathcal{N}$ as the given set, and $\perp$ as target set. The construction of nested explanations in [7] may lead to loose steps.

Loose steps are easy to detect in a sequence and can simply be removed from the explanation sequence until none remain.

4.2 Pertinence

After constructing a sequence and filtering out loose steps, the output of a step may still contain irrelevant literals as not *all* output literals have to be used later in the sequence. This requires a user to “waste” effort in understanding why useless literals are implied by the step. Clearly, this is undesired, motivating the following definition.

**► Property 3 (Pertinent).** An explanation sequence is pertinent if (i) $\mathcal{G}$ and the $N_i$ are pairwise disjoint, and (ii) all output literals are part of the target or the input of some following step.

The first condition states that each literal should only be derived once (and obviously, given literals should not be rederived). The second condition states that everything we derived should be used. To transform an explanation sequence into a pertinent one, we can simply loop over the sequence and remove any outputs not satisfying any of the conditions in the definition. Removing loose steps from a pertinent sequence corresponds to removing steps with an empty output set. In Example 3, the first step of sequence 1a is not pertinent, as the literal $w \neq 3$ in its output is not part of any input or of the target. For satisfiable CSPs, the algorithms described in the literature [7] are guaranteed to produce a pertinent explanation sequence, but not so for unsatisfiable CSPs.

5 Relaxation-based filtering

Informally, a sequence is pertinent if it derives only those literals that are (indirectly) needed for the target. However, a step in a sequence may also contain input literals that are not needed to derive its output. Naturally, this is an undesired property of explanation steps as it requires the user to process more literals than needed.

**► Property 4 (Sparse).** An explanation step $(\mathcal{E}, S, \mathcal{N})$ is sparse if no step $(\mathcal{E}', S, \mathcal{N})$ with $\mathcal{E}' \subsetneq \mathcal{E}$ exists. An explanation sequence is sparse if all its steps are.

The sequence of Example 4 contains non-sparse steps, e.g., the second step does not need the literals $p \in [0..1]$ as input to derive $z \neq 3$. The sequences in Example 7 are both sparse.

Previous approaches [7, 15] construct explanation steps by minimizing the input $\mathcal{E}$, making the resulting steps and sequences sparse. In contrast, CONSTRUCT-GREEDY constructs a maximal sequence which typically has non-sparse steps due to non-useful output literals of previous steps. We address this by a *relaxation-based filtering* algorithm, which calculates a so-called *relaxation* for each step.

**► Definition 6 (Relaxation).** A relaxation of a step $(\mathcal{E}, S, \mathcal{N})$ is a step $(\mathcal{E}', S, \mathcal{N})$ with $\mathcal{E}' \subseteq \mathcal{E}$.

An explanation step is sparse if and only if it allows no strict relaxations. Given a set of implied literals $\mathcal{N}$, let $\neg \mathcal{N}$ be the constraint denoting that at least one of the literals in $\mathcal{N}$ is false, i.e., $\neg \mathcal{N} = \bigvee_{x \neq v \in \mathcal{N}} x = v$. Then, $(\mathcal{E}, S, \mathcal{N})$ is sparse if for every subset $\mathcal{E}' \subsetneq \mathcal{E}$ it holds that $\mathcal{E}' \cup S \cup \{\neg \mathcal{N}\}$ is satisfiable, as in that case, $\mathcal{E}' \cup S \not= \mathcal{N}$ and $\mathcal{E}'$ cannot form a relaxation for $(\mathcal{E}, S, \mathcal{N})$.
<table><caption>Algorithm 2 RELAXATION-BASED FILTERING.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Input:</td><td>Explanation sequence Seq, given set G, target set T</td></tr><tr><td>1</td><td>Req &larr; T \ G</td></tr><tr><td>2</td><td>for (E<sub>i</sub>, S<sub>i</sub>, N<sub>i</sub>) &isin; reverseSeq) do</td></tr><tr><td>3</td><td>N<sub>i</sub> &larr; Req &cap; N<sub>i</sub></td></tr><tr><td>4</td><td>if N<sub>i</sub> is empty then</td></tr><tr><td>5</td><td>drop (E<sub>i</sub>, S<sub>i</sub>, N<sub>i</sub>) from Seq</td></tr><tr><td>6</td><td>else</td></tr><tr><td>7</td><td>E<sub>i</sub> &larr; MUS(soft: E<sub>i</sub>, hard: S<sub>i</sub> &cup; {&not;N<sub>i</sub>})</td></tr><tr><td>8</td><td>N<sub>i</sub> &larr; FULLPROPAGATE(E<sub>i</sub> &cup; S<sub>i</sub>)</td></tr><tr><td>9</td><td>Req &larr; ((Req \ N<sub>i</sub>) &cup; E<sub>i</sub>) \ G</td></tr></tbody></table>

To calculate sparse steps, we need to find a relaxation for each step. For this, we can calculate a MUS with $\mathcal{E}$ as soft constraints and $S \cup \{\neg \mathcal{N}\}$ as hard constraints. This MUS is a subset of $\mathcal{E}$ that yields a sparse step.

A straightforward algorithm to make all steps sparse would be to iterate over the steps and calculate such a MUS. However, the disadvantage of this is that sequences may no longer be pertinent: if the input of some step is reduced, then the output of a previous step may no longer be needed. And conversely, naively making a sequence pertinent may make it no longer sparse, as a smaller output for a step may reduce the number of input literals needed.

Instead, we propose a method that makes a sequence sparse by relaxing its steps from back to front. Additionally, this will make some steps loose, allowing the method to remove these steps. Hence, we call it *RELAXATION-BASED FILTERING*.

Algorithm 2 shows the pseudocode. It loops from back to front over the sequence and keeps a set of required literals *Req* that still need to be explained by some step earlier in the sequence. *Req* is initialized as the target set *T*. For each step in the sequence, the algorithm first keeps only those output literals that are also in *Req*, since keeping more would make the sequence non-pertinent (they are not required for the rest of the sequence). Next, it relaxes the step by calculating a MUS.

In the general case, many such MUSes may exist, and choosing one MUS over the other can have a profound impact on the sequence. This is clearly visible in Example 7 where the difference between both sequences arises from choosing a different MUS to relax the last step.

In this setting, we want a MUS allowing us to keep the set of required literals *Req* as small as possible. The reason for this is twofold: we want to keep future relaxations as small as possible and increase the number of detected loose steps. To compute a MUS satisfying this preference, we can split the calculation into two parts: first, minimize the extra input literals needed, and next minimize the subset of *Req*. Other techniques for this exist and include the OCUS-algorithm from [15].

The above process yields a sparse step deriving all required literals the step originally derived. However, the sparse input and constraints may imply more output literals than the originally *new* literals. To detect this, we make $\mathcal{N}$ maximal in line 8 of the algorithm. Any output this step can derive in this way can be removed from *Req*. But any input the step uses should be added to *Req*. This is what happens in the last line.

*RELAXATION-BASED FILTERING* guarantees sparsity by design since the input of each step is generated by a MUS call. Transforming the sequence into a pertinent one can be done trivially by removing literals from each $\mathcal{N}$ if they are not used later in the sequence or if they are derived as well at some point earlier in the sequence, as argued in Section 4.2.
## 11:10 Simplifying Step-Wise Explanation Sequences

Ensuring pertinence of this sequence in this way will not break sparsity in case the original sequence was maximal (which is guaranteed by our construction). To see this note that by definition in a maximal sequence, each literal is derived as soon as it can possibly be derived. Hence, any literal in $\mathcal{N}$ after Line 3 of the algorithm will not be derived at any earlier point in the sequence. Since these literals are required, they are part of some later input, and hence will not be removed from the output. For this reason, the step will stay sparse after making the sequence pertinent.

If the input sequence was maximal, then after RELAXATION-BASED FILTERING and making the sequence pertinent, there will also be no loose steps. To see this, note if RELAXATION-BASED FILTERING keeps a step, it is because some of its newly derived literals were required. If the input sequence is maximal, all literals are derived as early in the sequence as possible, meaning that this literal will not be derived earlier (and hence the pertinence making of the sequence will never remove it).

**► Example 7.** Below are two filtered versions of the maximal sequence in Example 4

<table><caption>**Table 3** Two sparse explanation sequences that could be produced using RELAXATION-BASED FILTERING by filtering the maximal sequence of example 4.</caption><thead><tr><th>i</th><th>ε</th><th>S</th><th>N</th><th>i</th><th>ε</th><th>S</th><th>N</th></tr></thead><tbody><tr><td>1</td><td>∅</td><td>p + q ≤ 1</td><td>p ∈ [0..1]<br>q ∈ [0..1]</td><td rowspan="2">1</td><td>∅</td><td>3s + p + r ≤ 1</td><td>p ∈ [0..1]<br>r ∈ [0..1]<br>s ∈ [0..1]</td></tr><tr><td>2</td><td>∅</td><td>3s + p + r ≤ 1</td><td>r ≠ 2, r ≠ 3</td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td>p ∈ [0..1]<br>q ∈ [0..1]<br>r ∈ [0..1]</td><td>alldiff(p, q, r, s)</td><td>⊥</td><td>2</td><td>p ∈ [0..1]<br>r ∈ [0..1]<br>s ∈ [0..1]</td><td>alldiff(p, q, r, s)</td><td>⊥</td></tr></tbody></table>

# 6 Deletion-based filtering

The techniques proposed so far are essential to make explanation sequences more easily understandable by users. The experiments in Section 7 will show that RELAXATION-BASED FILTERING significantly reduces the number of steps for many sequences. However, more steps could be removed if we allow to change the input and output of other, later steps. In the left sequence of Example 7, step 1 can be omitted from the sequence, as shown on the right in that same example. However, depending on the exact MUSes (and hence, relaxations) calculated by RELAXATION-BASED FILTERING, it may miss this, and the step remains redundant.

**► Property 5 (Redundant).** Given an explanation sequence $\langle (\mathcal{E}_i, S_i, \mathcal{N}_i) \rangle_{1 \le i \le n}$ of $\mathcal{T}$ from $\mathcal{G}$, the step $(\mathcal{E}_j, S_j, \mathcal{N}_j)$ is redundant if the maximal explanation sequence arising from $\mathcal{G}$ and $\langle S_i \rangle_{i \ne j}$ is an explanation sequence of $\mathcal{T}$. An explanation sequence is redundant if one of its steps is redundant.

Intuitively, a redundant step is a step for which all useful literals it derives could also be derived by later steps in the sequence. The first and second steps in Example 4 are redundant, as there exists a maximal sequence with the constraints of steps three and four: it is exactly the right sequence of Example 7.

The literature approaches [7, 15] may very well lead to redundant sequences, as once a step is added to the sequence under construction, no check is later made whether the step could be subsumed by one added later. Notice this also holds for the application of explaining the unique solution of a logic puzzle, whereas loose steps could not occur in that setting.
■ Algorithm 3 DELETION-BASED FILTERING.

**Input:** Explanation sequence *Seq* with maximal inputs and outputs, given set *G*,
target set *T*

1 **for** (*E*<sub>*i*</sub>, *S*<sub>*i*</sub>, *N*<sub>*i*</sub>) ∈ *reverse*(*Seq*) **do**
2 | let *Sub* be a copy of *Seq* from *i* + 1 to end
3 | **if** *TRYDELETION*(*E*<sub>*i*</sub>, *Sub*) **then**
4 | | shrink *Seq* to size *i* − 1
5 | | append *Sub* to *Seq*
6 **Function** *TRYDELETION*(*E*<sub>*i*</sub>, *Sub*):
7 | *N*<sub>*i*</sub> ← ∅
8 | **for** (*E*<sub>*j*</sub>, *S*<sub>*j*</sub>, *N*<sub>*j*</sub>) ∈ *Sub* **do**
9 | | *E*<sub>*j*</sub> ← *G* ∪ *E*<sub>*j*−1</sub> ∪ *N*<sub>*j*−1</sub>
10 | | *N*<sub>*j*</sub> ← FULLPROPAGATE(*E*<sub>*j*</sub> ∪ *S*<sub>*j*</sub>)
11 | | **if** *T* ⊆ *N*<sub>*j*</sub> **then**
12 | | | **return** True
13 | | **return** False

Algorithm DELETION-BASED FILTERING filters redundant steps. For this, it assumes a maximal sequence as input – if not, the sequence can be made maximal by a linear iteration over all steps from front to back. Next, DELETION-BASED FILTERING iterates over the sequence from back to front. In every iteration, it leaves out the current step, calculates the resulting maximal explanation sequence, and checks if it still explains the target. Calculating the resulting maximal sequence requires a call to a full propagation routine for $\mathcal{E} \cup S$ for each later step to derive the new outputs.

In the worst case where no step can be removed from a sequence of length *n*, this results in *n*(*n* − 1) / 2 calls to a full propagation routine. Although the pseudo-Boolean solver Exact allows for stateful computation of these maximal outputs, this remains a computationally expensive task and several optimizations are needed to make the algorithm work in practice, as we now explain.

**6.1 Necessary condition on the existence of a sequence**

For a maximal sequence with given set $\mathcal{G}$ and target set $\mathcal{T}$, DELETION-BASED FILTERING has to check for a step $(\mathcal{E}_i, S_i, \mathcal{N}_i)$ whether the maximal sequence determined by $\langle S_j \rangle_{i<j \le n}$ and $\mathcal{E}_i$ as given set, still derives $\mathcal{T}$. A necessary condition is that $\mathcal{E}_i \cup \bigcup_{i<j \le n} S_j \models \mathcal{T}$ – if not, $\mathcal{T}$ cannot be explained with just the constraints from the remaining steps in the sequence. This can be checked via a simple solve call with the constraints $\mathcal{E}_i \cup \bigcup_{i<j \le n} S_j \cup \neg \mathcal{T}$. In case this returns a solution, step *i* cannot be dropped from the sequence.

6.2 Partial propagation

Full propagation – calculating the maximal output of a step – is an expensive operation.
However, weaker and more efficient propagation algorithms are researched [5, 38] and are im-
plemented in several state-of-the-art constraint solvers. When DELETION-BASED FILTERING
has to check whether a target set $\mathcal{T}$ can be derived by the maximal sequence arising from
$\langle S_i \rangle_{j \le i \le n}$ and given set $\mathcal{E}_i$, it can first compute a “partial maximal” sequence using a com-
putationally cheaper propagation algorithm. If the target set is explained according to this
weaker sequence, it also will be explained by the maximal sequence, so the step can be
dropped and DELETION-BASED FILTERING can continue to the next iteration.
## 6.3 Caching

As explained in Section 3.2, we can cache the result of calls to a full propagation algorithm. In addition, we can cache the result of computing a maximal sequence (containing multiple maximal outputs), and we can cache the result of checking the necessary condition for a maximal sequence as described in Section 6.1. These caching strategies become particularly useful given the monotonicity of entailment, as formalized in the following lemma.

**Lemma 8.** For any set of constraints $\mathcal{S}$ and any two sets of literals $\mathcal{E} \subseteq \mathcal{E}',$ if $\mathcal{E} \cup \mathcal{S} \models \mathcal{N}$ then $\mathcal{E}' \cup \mathcal{S} \models \mathcal{N}$.

This entails that if a maximal sequence determined by $\langle S_j \rangle_{i<j \le n}$ with given set $\mathcal{E}_i$ entails a target set $\mathcal{T}$, so will the maximal sequence using the same constraints with a superset of $\mathcal{E}_i$. Conversely, if that maximal sequence does *not* entail $\mathcal{T}$, neither will the maximal sequence with any subset of $\mathcal{E}_i$.

We propose to generate explanation sequences using *CONSTRUCT-GREEDY* and filter these sequences using *DELETION-BASED FILTERING* so no *redundant* steps remain. Lastly, these non-redundant sequences need to undergo *RELAXATION-BASED FILTERING* so all steps are made sparse and can be made pertinent. In Appendix A, we construct a formal argument to prove why this pipeline ensures *atomicity*, *sparsity*, *pertinence* and *non-redundancy* for general explanation sequences.

Notice that filtering explanation sequences using *DELETION-BASED FILTERING* can potentially filter more steps compared to *RELAXATION-BASED FILTERING* but may increase the complexity of individual steps as more literals are used/explained at a time. Nevertheless, the constraint sets of each step are not altered by any of the algorithms we proposed. These are entirely determined by the construction algorithm.

# 7 Experimental results

We evaluate and review the algorithms and ideas presented in previous sections using three benchmark sets consisting of unsatisfiable CSP instances.

**Sudoku.** 50 9x9 sudoku instances in which an empty cell is given a wrong value, such that this variable assignment is non-trivial, i.e., it does not directly falsify any constraint. The original puzzles were generated by the QQwing [35] tool using the “Intermediate” difficulty setting. The sudoku constraints are modelled with AllDifferent global constraints.

**Jobshop.** 50 jobshop instances generated by the approach in [40] with 5 machines and 5 tasks per job. The time horizon is equal to 50 units. Unsatisfiability arises from restricting the makespan to a better-than-optimal value. All instances are modelled using Cumulative constraints [1] as is usual in scheduling problems.

**Debug.** Using the diverse set of CSP models from Håkan Kjellerstrand,<fnref n="2" /> we construct a set of 202 unsatisfiable CSPs by introducing artificial errors in the models to mimic a modelling error by a user – an unsatisfiability bug. We follow a similar approach to [29]. The errors are introduced such that each constraint separately remains satisfiable.

<footnote n="2">² http://www.hakank.org/cpmpy/</footnote>
These three benchmark families correspond to three use cases where step-wise explanations
can be of help. **Sudoku** simulates the counterfactual explanation provided during an
interactive collaboration (e.g., [23]) by a human agent and an automated system. Here the
user wants to know why the system derived that a certain value was no longer possible for a
certain value. **Jobshop** corresponds to explaining to a user why an inferred objective value
is optimal for some optimization problem. **Debug** simulates the situation where a user is
modelling a CSP but discovers the model has no solutions and the user needs to "debug" it.

For every instance, we compute three explanation sequences using different random seeds.
As metrics, we consider the number of steps in sequences, as well as the complexity of
individual steps. This first metric is different from what literature approaches consider [7, 15]
as in their use-case of explaining a unique solution, they focus only on the complexity of
individual steps.

RQ1. How many redundant steps do explanation sequences of unsatisfiable CSPs contain?

RQ2. How well can the proposed filtering techniques simplify explanation sequences?

RQ3. How much time do the proposed algorithms take to run?

RQ4. How does simplifying the set of input constraints affect the explanation sequence?

## 7.1 Experimental setup

All experiments were run on a single core of an Intel(R) Xeon(R) Silver 4214 CPU with 128GB of RAM. FULLPROPAGATE is handled by the stateful pruneDomains functionality of the Exact solver [12, 14] v1.0.0.<fnref n="3" /> Partial propagation is calculated with OR-Tools [36] v9.6 presolve routine. All algorithms are implemented using a custom branch of CPMpy [19] v0.9.12, embedded in Python v3.10.9. To calculate MUSes we employ the mus tool provided by CPMpy . Cardinal-minimal MUSes are computed using the hitting-set approach of [25]. All code and benchmarks are available on GitHub: https://github.com/ML-KULeuven/SimplifySeq

## 7.2 Approaches under investigation

Throughout this section, we investigate all algorithms and combinations of them presented in
this paper. All sequences are generated using CONSTRUCT-GREEDY and are indicated using
**Original** in each of the plots. By making each step sparse, we can mimic the explanation
sequences generated by approaches described in the literature [7, 15]. This can simply be
done by finding a minimal unsatisfiable subset with hard constraints $S \cup \{\neg N\}$ and soft
constraints $\mathcal{E}$. Afterwards filtering *loose* steps and making the sequence *pertinent* using
the approach described in Section 4 is denoted as **Filter simple** and is considered current
state-of-the-art. We consider two versions of RELAXATION-BASED FILTERING differing in
the type of MUS extracted: a MUS and a smallest MUS (SMUS), shown as **MUS Relax**
resp. **SMUS Relax**. Lastly, we investigate sequences exposed to the entire pipeline using
DELETION-BASED FILTERING and both versions of RELAXATION-BASED FILTERING. These
are shown as **Del + MUS Relax** and **Del + SMUS Relax**.

## 7.3 Redundant steps

The first research question is straightforward: how many redundant steps are there in naively
generated sequences, and how many can we filter using our proposed techniques? From
Figure 1 it is clear greedy explanation sequences contain a lot of redundant steps. For all

<footnote n="3">³ https://gitlab.com/JoD/exact git commit 34c4fa46</footnote>
11:14 Simplifying Step-Wise Explanation Sequences

<figure><img src="image_106.png" /><figcaption>Figure 1 Comparison of different filter algorithms in terms of **number of steps** in the sequence.</figcaption></figure>

benchmarks, the deletion-based method leads the way in terms of the number of steps filtered from the sequence. The price to be paid is a significant increase in computational effort for some benchmarks and more complex individual steps as outlined in the next sections. Simple filtering of loose steps as described in Section 4 is able to reduce the number of steps, but falls short compared to fully-fledged relaxation-based methods. This is most prominently visible in the Jobshop benchmark set. Interestingly, both variants of RELAXATION-BASED FILTERING (MUS Relax and SMUS Relax) show similar performance in terms of steps left in the sequence after filtering.

7.4 Complexity of individual steps

<figure><img src="image_107.png" /><figcaption>(a) Comparison in terms of **average number of input literals** for each step in an explanation sequence.</figcaption></figure>

<figure><img src="image_108.png" /><figcaption>(b) Comparison in terms of **average number of output literals** for each step in an explanation sequence.</figcaption></figure>

<caption>Figure 2 Complexity of individual steps in terms of literals, algorithms producing maximal sequences are omitted from these plots.</caption>

Next, we investigate how different filtering methods compare in terms of input and output literals for each step. The simple filtering method is able to keep the number of input literals relatively low. As expected, it again falls short compared to relaxation-based methods as **Filter simple** does not modify the set of input literals based on the output literals that are actually useful. While the deletion-based methods are able to remove more steps, it increases the number of input literals for the remaining steps compared to relaxation-based methods. This is similar to the number of output literals, shown in Figure 2b. Deletion-based methods require steps to derive more output while relaxation-based methods work better for this metric. We argue this larger number of output literals is less of a concern compared to the increase of input literals as steps can be split up so they derive fewer literals at a time.
The left-hand side of Table 4 shows that the number of constraints is low for each individual step, ranging from 1 to a few constraints propagated at a time. Notice how during filtering, some steps occupying little constraints are deleted, while the largest step for each sequence is always present after filtering still.

## 7.5 Explaining a MUS

Table 4 Effect of explaining a MUS of the unsatisfiable problem instead of all constraints in the problem. Metrics are shown before filtering → after filtering using DELETION-BASED FILTERING.

<table><thead><tr><th rowspan="2"></th><th colspan="3">Explaining all constraints</th><th colspan="3">Explaining a MUS</th></tr><tr><th>#steps</th><th>Max |S|</th><th>Avg |S|</th><th>#steps</th><th>Max |S|</th><th>Avg |S|</th></tr></thead><tbody><tr><td>Sudoku</td><td>79.62 → 26.33</td><td>1.00 → 1.00</td><td>1.00 → 1.00</td><td>46.10 → 27.87</td><td>2.04 → 2.04</td><td>1.07 → 1.08</td></tr><tr><td>Jobshop</td><td>118.03 → 35.69</td><td>1.02 → 1.02</td><td>1.00 → 1.00</td><td>96.55 → 38.10</td><td>1.35 → 1.35</td><td>1.01 → 1.01</td></tr><tr><td>Debug</td><td>18.70 → 5.91</td><td>1.32 → 1.32</td><td>1.10 → 1.17</td><td>6.86 → 5.51</td><td>1.48 → 1.48</td><td>1.20 → 1.26</td></tr></tbody></table>

In the literature, the most common way to explain why a model is unsatisfiable is to first extract a MUS. What happens to step-wise explanations when we first extract a MUS and build the explanation sequence from that MUS?

The result of this experiment is summarized in Table 4. Here we observe that when extracting a MUS, the explanation sequence is much shorter before any filtering is applied. Interestingly, explaining the full CSP may actually lead to shorter sequences after filtering, as observed for **Sudoku** and **Jobshop**. Equally remarkable, limiting the explanation sequence to a MUS increases the number of constraints used in individual steps. This is most clear from **Sudoku** where the maximal number of constraints in a step increases from 1 constraint to 2.04 constraints on average. Intuitively, limiting the constraints with which to build an explanation sequence to a smaller unsatisfiable subset reduces the ease with which explanations can be built. E.g., the **Sudoku** problem contains many redundant constraints that are implied by some subset of other constraints. Nevertheless, these redundant constraints can be very useful to elegantly explain something to a user.

## 7.6 Runtime and optimizations

Figure 3 displays the runtime for each of the filtering algorithms for all benchmark sets. For **Jobshop**, the deletion-based method is computationally more expensive compared to any of the relaxation-based variants. This is not the case for **Sudoku** and **Debug**. As **Cumulative** constraints are harder to fully propagate compared to **AllDifferent** constraints, the quadratic number of full-propagation calls is much more prominent for the **Jobshop** benchmark. For all benchmarks, calculating a smallest-minimal MUS requires more computational effort while not impacting the sequence enough to justify this extra cost.

The optimizations implemented for DELETION-BASED FILTERING as described at the end of Section 6 are most effective with the **Debug** benchmark, where only 7% of the theoretic number of full-propagation calls are executed. For **Sudoku** and **Jobshop**, this is 19% resp. 34%. Overall, the proposed optimizations are necessary to reach the observed performance.

# 8 Conclusion and Future work

We investigated the problem of step-wise explaining unsatisfiable problems, with applications in debugging unsatisfiability, explaining optimality, and explanations in interactive configuration. To ensure a *simple* step-wise explanation, both the number of steps and the
11:16 Simplifying Step-Wise Explanation Sequences

<figure><img src="image_109.png" /><figcaption>Figure 3 ECDF plot relating solved instances to runtime for different filtering algorithms.</figcaption></figure>

number of constraints and literals in each explanation step must be kept low. We proposed the formal properties of *atomicity*, *pertinence*, *sparsity*, and (ir)redundancy, to which an explanation sequence should strive to adhere. For this, we proposed a workflow combining greedy construction, deletion-based filtering and relaxation-based filtering that guarantees all properties are satisfied.

We presented extensive experimental results comparing current approaches to our proposed techniques. This showed filtering of explanation sequences is especially essential in the challenging setting of explaining unsatisfiability and is key in creating easy-to-understand explanation sequences. Our methods are able to considerably shorten explanation sequences with deletion-based filtering leading to sequences without any redundant steps. Compared to relaxation-based techniques which may produce longer sequences, deletion-based methods may lead to more complex steps after filtering.

For future work, observe that none of the filtering algorithms considered in this paper alter the set of constraints for explanation steps, nor the order of explanation steps. These are entirely determined by the construction algorithm. Changing either of these aspects might further simplify an explanation sequence. We leave these areas for future work, where both new sequence construction methods and more elaborate post-processing are viable options to explore.

Finally, the small number of constraints in most explanation steps shows that computing explanation sequences seem a worthwhile approach to help a user understand unsatisfiable constraints. Still, the expressivity of explanation steps deserves further study with our results opening the door for further investigation into what domain experts perceive as “simple” explanation sequences. User studies or learning techniques can investigate how it helps them to understand the interplay of complex constraints.

References

1 Abderrahmane Aggoun and Nicolas Beldiceanu. Extending CHIP in order to solve complex scheduling and placement problems. In Jean-Paul Delahaye, Philippe Devienne, Philippe Mathieu, and Pascal Yim, editors, *JFPL'92, 1eres Journées Francophones de Programmation Logique, 25-27 Mai 1992, Lille, France*, page 51, 1992.

2 Mario Alviano, Carmine Dodaro, Johannes Klaus Fichte, Markus Hecher, Tobias Philipp, and Jakob Rath. Inconsistency proofs for ASP: the ASP - DRUPE format. *Theory Pract. Log. Program.*, 19(5-6):891–907, 2019. doi:10.1017/S1471068419000255.

3 Haniel Barbosa, Andrew Reynolds, Gereon Kremer, Hanna Lachnitt, Aina Niemetz, Andres Nötzli, Alex Ozdemir, Mathias Preiner, Arjun Viswanathan, Scott Viteri, Yoni Zohar, Cesare Tinelli, and Clark W. Barrett. Flexible proof production in an industrial-strength SMT solver. In Jasmin Blanchette, Laura Kovács, and Dirk Pattinson, editors, *Automated Reasoning - 11th International Joint Conference, IJCAR 2022, Haifa, Israel, August 8-10, 2022, Proceedings,*
volume 13385 of *Lecture Notes in Computer Science*, pages 15–35. Springer, 2022. doi: 10.1007/978-3-031-10769-6_3.

4 Jeremias Berg, Bart Bogaerts, Jakob Nordström, Andy Oertel, and Dieter. Vandesande. Certified core-guided MaxSAT solving. In *Proceedings of the 29th International Conference on Automated Deduction (CADE)*, 2023. Accepted for publication.

5 Christian Bessiere. Constraint propagation. In Francesca Rossi, Peter van Beek, and Toby Walsh, editors, *Handbook of Constraint Programming*, volume 2 of *Foundations of Artificial Intelligence*, pages 29–83. Elsevier, 2006. doi:10.1016/S1574-6526(06)80007-6.

6 Christian Bessiere, Clément Carbonnel, Martin C. Cooper, and Emmanuel Hebrard. Complexity of Minimum-Size Arc-Inconsistency Explanations. In Christine Solnon, editor, *28th International Conference on Principles and Practice of Constraint Programming (CP 2022)*, volume 235 of *Leibniz International Proceedings in Informatics (LIPIcs)*, pages 9:1–9:14, Dagstuhl, Germany, 2022. Schloss Dagstuhl – Leibniz-Zentrum für Informatik. doi:10.4230/LIPIcs.CP.2022.9.

7 Bart Bogaerts, Emilio Gamba, and Tias Guns. A framework for step-wise explaining how to solve constraint satisfaction problems. *Artif. Intell.*, 300:103550, 2021. doi:10.1016/j.artint.2021.103550.

8 Bart Bogaerts, Stephan Gocht, Ciaran McCreesh, and Jakob Nordström. Certified symmetry and dominance breaking for combinatorial optimisation. In *Thirty-Sixth AAAI Conference on Artificial Intelligence, AAAI 2022, Thirty-Fourth Conference on Innovative Applications of Artificial Intelligence, IAAI 2022, The Twelveth Symposium on Educational Advances in Artificial Intelligence, EAAI 2022 Virtual Event, February 22 - March 1, 2022*, pages 3698–3707. AAAI Press, 2022. URL: https://ojs.aaai.org/index.php/AAAI/article/view/20283.

9 Jens Claes, Bart Bogaerts, Rocío Sáez Canoy, Emilio Gamba, and Tias Guns. Zebratutor: Explaining how to solve logic grid puzzles. In Katrien Beuls, Bart Bogaerts, Gianluca Bontempi, Pierre Geurts, Nick Harley, Bertrand Lebichot, Tom Lenaerts, Gilles Louppe, and Paul Van Ecke, editors, *Proceedings of the 31st Benelux Conference on Artificial Intelligence (BNAIC 2019) and the 28th Belgian Dutch Conference on Machine Learning (Benelearn 2019)*, Brussels, Belgium, November 6-8, 2019, volume 2491 of *CEUR Workshop Proceedings*. CEUR-WS.org, 2019. URL: http://ceur-ws.org/Vol-2491/demo96.pdf.

10 Jens Claes, Bart Bogaerts, Emilio Gamba, Rocío Sáez Canoy, and Tias Guns. Human-oriented solving and explaining of logic grid puzzles, November 2019. BNAIC 2019 ; Conference date: 07-11-2019 Through 08-11-2019.

11 Broes De Cat, Bart Bogaerts, Maurice Bruynooghe, Gerda Janssens, and Marc Denecker. Declarative logic programming. In Michael Kifer and Yanhong Annie Liu, editors, *Declarative Logic Programming: Theory, Systems, and Applications*, chapter Predicate Logic As a Modeling Language: The IDP System, pages 279–323. Association for Computing Machinery and Morgan & Claypool, New York, NY, USA, 2018. doi:10.1145/3191315.3191321.

12 Jo Devriendt. Exact solver, 2023. URL: https://gitlab.com/JoD/exact.

13 William Dumez, Simon Vandevelde, and Joost Vennekens. Step-wise explanations of sudokus using IDP. In *BNAIC/BeNeLearn*, November 2022.

14 Jan Elffers and Jakob Nordström. Divide and conquer: Towards faster pseudo-Boolean solving. In Jérôme Lang, editor, *Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, IJCAI 2018, July 13-19, 2018, Stockholm, Sweden*, pages 1291–1299. ijcai.org, 2018. doi:10.24963/ijcai.2018/180.

15 Emilio Gamba, Bart Bogaerts, and Tias Guns. Efficiently explaining CSPs with unsatisfiable subset optimization. In Zhi-Hua Zhou, editor, *Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21*, pages 1381–1388. International Joint Conferences on Artificial Intelligence Organization, August 2021. Main Track. doi:10.24963/ijcai.2021/191.

16 Martin Gebser, Benjamin Kaufmann, and Torsten Schaub. Conflict-driven answer set solving: From theory to practice. *Artif. Intell.*, 187:52–89, 2012. doi:10.1016/j.artint.2012.04.001.
11:18 Simplifying Step-Wise Explanation Sequences

17 Leilani H. Gilpin, David Bau, Ben Z. Yuan, Ayesha Bajwa, Michael A. Specter, and Lalana Kagal. Explaining explanations: An overview of interpretability of machine learning. In Francesco Bonchi, Foster J. Provost, Tina Eliassi-Rad, Wei Wang, Ciro Cattuto, and Rayid Ghani, editors, 5th IEEE International Conference on Data Science and Advanced Analytics, DSAA 2018, Turin, Italy, October 1-3, 2018, pages 80-89. IEEE, 2018. doi:10.1109/DSAA.2018.00018.

18 Stephan Gocht, Ciaran McCreesh, and Jakob Nordström. An auditable constraint programming solver. In Christine Solnon, editor, 28th International Conference on Principles and Practice of Constraint Programming, CP 2022, July 31 to August 8, 2022, Haifa, Israel, volume 235 of LIPIcs, pages 25:1-25:18. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2022. doi:10.4230/LIPIcs.CP.2022.25.

19 Tias Guns. Increasing modeling language convenience with a universal n-dimensional array, cppy as python-embedded example. In Proceedings of the 18th workshop on Constraint Modelling and Reformulation at CP (Modref 2019), volume 19, 2019.

20 Tias Guns, Milan Pesa, Maxime Mulamba, Ignace Bleukx, Emilio Gamba, and Senne Berden. Sudoku assistant - an AI-powered app to help solve pen-and-paper sudokus, 2022.

21 Sharmi Dev Gupta, Begum Genc, and Barry O'Sullivan. Explanation in constraint satisfaction: A survey. In Zhi-Hua Zhou, editor, Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI 2021, Virtual Event / Montreal, Canada, 19-27 August 2021, pages 4400-4407. ijcai.org, 2021. doi:10.24963/ijcai.2021/601.

22 Sharmi Dev Gupta, Begum Genc, and Barry O'Sullivan. Finding counterfactual explanations through constraint relaxations. CoRR, abs/2204.03429, 2022. doi:10.48550/arXiv.2204.03429.

23 Pieter Van Hertum, Ingmar Dasseville, Gerda Janssens, and Marc Denecker. The KB paradigm and its application to interactive configuration. Theory Pract. Log. Program., 17(1):91-117, 2017. doi:10.1017/S1471068416000156.

24 Marijn J. H. Heule. Proofs of unsatisfiability. In Armin Biere, Marijn Heule, Hans van Maaren, and Toby Walsh, editors, Handbook of Satisfiability - Second Edition, volume 336 of Frontiers in Artificial Intelligence and Applications, pages 635-668. IOS Press, 2021. doi:10.3233/FAIA200998.

25 Alexey Ignatiev, Alessandro Previti, Mark H. Liffiton, and João Marques-Silva. Smallest MUS extraction with minimal hitting set dualization. In Gilles Pesant, editor, Principles and Practice of Constraint Programming - 21st International Conference, CP 2015, Cork, Ireland, August 31 - September 4, 2015, Proceedings, volume 9255 of Lecture Notes in Computer Science, pages 173-182. Springer, 2015. doi:10.1007/978-3-319-23219-5_13.

26 Ulrich Junker. Quickxplain: Conflict detection for arbitrary constraint propagation algorithms. In IJCAI'01 Workshop on Modelling and Solving problems with constraints, volume 4. Citeseer, 2001.

27 Pat Langley, Ben Meadows, Mohan Sridharan, and Dongkyu Choi. Explainable agency for intelligent autonomous systems. In Satinder Singh and Shaul Markovitch, editors, Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, February 4-9, 2017, San Francisco, California, USA, pages 4762-4764. AAAI Press, 2017. URL: http://aaai.org/oocs/index.php/IAAI/IAAI17/paper/view/15046.

28 Niklas Lauffer and Ufuk Topcu. Human-understandable explanations of infeasibility for resource-constrained scheduling problems. In ICAPS 2019 Workshop XAIP, 2019.

29 Kevin Leo and Guido Tack. Debugging unsatisfiable constraint models. In Domenico Salvagnin and Michele Lombardi, editors, Integration of AI and OR Techniques in Constraint Programming - 14th International Conference, CPAIOR 2017, Padua, Italy, June 5-8, 2017, Proceedings, volume 10335 of Lecture Notes in Computer Science, pages 77-93. Springer, 2017. doi:10.1007/978-3-319-59776-8_7.
30 Mark H. Liffiton, Alessandro Previti, Ammar Malik, and João Marques-Silva. Fast, flexible MUS enumeration. *Constraints An Int. J.*, 21(2):223-250, 2016. doi:10.1007/s10601-015-9183-0.

31 Mark H. Liffiton and Karem A. Sakallah. Algorithms for computing minimal unsatisfiable subsets of constraints. *J. Autom. Reason.*, 40(1):1-33, 2008. doi:10.1007/s10817-007-9084-z.

32 João Marques-Silva. Minimal unsatisfiability: Models, algorithms and applications (invited paper). In *40th IEEE International Symposium on Multiple-Valued Logic, ISMVL 2010, Barcelona, Spain, 26-28 May 2010*, pages 9-14. IEEE Computer Society, 2010. doi:10.1109/ISMVL.2010.11.

33 João Marques-Silva. Logic-based explainability in machine learning. In Leopoldo E. Bertossi and Guohui Xiao, editors, *Reasoning Web. Causality, Explanations and Declarative Knowledge - 18th International Summer School 2022, Berlin, Germany, September 27-30, 2022, Tutorial Lectures*, volume 13759 of *Lecture Notes in Computer Science*, pages 24-104. Springer, 2022. doi:10.1007/978-3-031-31414-8_2.

34 Olga Ohrimenko, Peter J. Stuckey, and Michael Codish. Propagation via lazy clause generation. *Constraints An Int. J.*, 14(3):357-391, 2009. doi:10.1007/s10601-008-9064-x.

35 Stephen Ostermiller. QQwing. URL: https://qqwing.com/.

36 Laurent Perron and Vincent Furnon. Or-tools, November 2022. URL: https://developers.google.com/optimization/.

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

38 Christian Schulte and Peter J. Stuckey. Efficient constraint propagation engines. *ACM Trans. Program. Lang. Syst.*, 31(1):2:1-2:43, 2008. doi:10.1145/1452044.1452046.

39 Ilankaikone Senthooran, Matthias Klapperstück, Gleb Belov, Tobias Czauderna, Kevin Leo, Mark Wallace, Michael Wybrow, and Maria Garcia de la Banda. Human-centred feasibility restoration. In Laurent D. Michel, editor, *27th International Conference on Principles and Practice of Constraint Programming, CP 2021, Montpellier, France (Virtual Conference), October 25-29, 2021*, volume 210 of *LIPIcs*, pages 49:1-49:18. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021. doi:10.4230/LIPIcs.CP.2021.49.

40 Eric Taillard. Benchmarks for basic scheduling problems. *European journal of operational research*, 64(2):278-285, 1993.

## A Satisfying all desirable properties

To transform a maximal sequence generated by CONSTRUCT-GREEDY to one that is *sparse*, *pertinent*, *atomic* and contains no *redundant* (and hence no *loose*) steps, we propose to generate sequences with CONSTRUCT-GREEDY, then filter those with DELETION-BASED FILTERING, finally followed by RELAXATION-BASED FILTERING.

The application of RELAXATION-BASED FILTERING in the last step guarantees sparsity. Modifying the output of steps to guarantee pertinence will not break sparsity as the output of DELETION-BASED FILTERING is maximal (see the argument of Section 5).

After the application of DELETION-BASED FILTERING, all redundant steps have been removed. In other words, for any *i*, there exists no strict subsequence $\langle S_i \rangle_{1 \le i \le n, i \ne j}$ for which its maximal sequence with the same given set $\mathcal{G}$ explains the same target set $\mathcal{T}$. Hence, RELAXATION-BASED FILTERING will not actually filter any more steps, as also shown in the experiments in Section 7. Moreover, RELAXATION-BASED FILTERING does not change any of the constraint subsets $S_i$, so the maximal sequence after RELAXATION-BASED FILTERING is equal to the one after DELETION-BASED FILTERING. Therefore, RELAXATION-BASED FILTERING preserves the guarantee of DELETION-BASED FILTERING that no redundant steps exist.

To argue atomicity, we will consider three sequences
11:20 Simplifying Step-Wise Explanation Sequences

■ $\langle (\mathcal{E}_i, S_i, \mathcal{N}_i) \rangle_{1 \le i \le n}$ denotes the sequence obtained after construction,

■ $\langle (\mathcal{E}'_i, S_i, \mathcal{N}'_i) \rangle_i$ not deleted the sequence obtained after DELETION-BASED FILTERING, and

■ $\langle (\mathcal{E}''_i, S_i, \mathcal{N}''_i) \rangle_i$ not deleted the sequence obtained after RELAXATION-BASED FILTERING.

Moreover, we will write $\mathcal{I}_i$ (resp. $\mathcal{I}'_i, \mathcal{I}''_i$) for $\mathcal{G} \cup \bigcup_{j<i} \mathcal{N}_j$ (resp. $\mathcal{G} \cup \bigcup_{j<i} \mathcal{N}'_j, \mathcal{G} \cup \bigcup_{j<i} \mathcal{N}''_j$). Intuitively, $\mathcal{I}_i$ represents the set of all facts available prior to step $i$. By definition of an explanation sequence, we have that $\mathcal{E}_i \subseteq \mathcal{I}_i$ (and similar for the other two sequences).

To argue that $(\mathcal{E}''_i, S_i, \mathcal{N}''_i)$ is still atomic, we first prove four claims.

1. $\mathcal{I}_i \supseteq \mathcal{I}'_i \supseteq \mathcal{I}''_i$;

2. $\mathcal{N}'_i \subseteq \mathcal{I}_{i+1}$ and $\mathcal{N}''_i \subseteq \mathcal{I}_{i+1}$

3. For each literal $l$ in $\mathcal{N}_i \setminus \mathcal{I}_i$ (i.e., every newly derived literal in step $i$), and every subset $S \subseteq S_i$, $\mathcal{I}_i \cup S \not\models l$;

4. If step $i$ is not deleted, then $\mathcal{N}'_i \cap (\mathcal{N}_i \setminus \mathcal{I}_i) \neq \emptyset$ and $\mathcal{N}''_i \cap (\mathcal{N}_i \setminus \mathcal{I}_i) \neq \emptyset$, i.e., at least one literal that was newly derived at step $i$ is still derived at that step.

Claim 1 follows from the fact that, by construction, the first two sequences are maximal (in their inputs as well as outputs), and hence derive everything that is derivable at each point. Since the constraint sets in the three sequences are the same, nothing extra can then be derived in the trimmed sequences. Claim 2 again follows directly from the fact that the initial sequence is maximal: what is derived at step $i$ can be at most $\mathcal{I}_i \cup \mathcal{N}_i$, which equals $\mathcal{I}_{i+1}$. Claim 3 holds by construction: our greedy construction algorithm only generates steps for which there is no smaller set of constraints that can derive something. Claim 4 holds since, using Claim 2, if that intersection would be empty, $\mathcal{N}'_i$ would consist only of literals of $\mathcal{I}_i$, in which case DELETION-BASED FILTERING would clearly delete this step.

Now, take any $i$. We continue to show that $(\mathcal{E}''_i, S_i, \mathcal{N}''_i)$ is in fact atomic. Using Claim 4, $\mathcal{N}''_i \cap (\mathcal{N}_i \setminus \mathcal{I}_i) \neq \emptyset$. From Claim 2 it follows that $\mathcal{N}''_i$ is the union of two sets: $O := \mathcal{N}''_i \cap \mathcal{I}_i$ and $D := \mathcal{N}''_i \cap (\mathcal{N}_i \setminus \mathcal{I}_i)$, where $O$ is the set of literals that were originally “old” at step $i$ (derived before $i$) and $D$ are those that were actually derived at step $i$. By Claim 4, $D$ is non-empty. Now assume towards contradiction that our step $(\mathcal{E}''_i, S_i, \mathcal{N}''_i)$ is not atomic. In that case, there is an explanation sequence that derives $\mathcal{N}''_i$ from $\mathcal{E}''_i$ in which each step uses a strict subset of $S_i$. From that sequence, consider the first step that derives an element from $D$. By definition, the input of that step can consist at most of $\mathcal{E}''_i \cup O$. Since $\mathcal{E}''_i \cup O \subseteq \mathcal{E}_i$, we find an explanation step that violates Claim 3, and we can conclude our proof by contradiction.
# Towards More Efficient Local Search for Pseudo-Boolean Optimization

Yi Chu

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

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

Institute of Software, Chinese Academy of Sciences, Beijing, China

Shaowei Cai<fnref n="1" />

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

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

State Key Laboratory of Computer Science, Institute of Software, Chinese Academy of Sciences, Beijing, China

Chuan Luo

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

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

School of Software, Beihang University, Beijing, China

Zhendong Lei

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

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

State Key Laboratory of Computer Science, Institute of Software, Chinese Academy of Sciences, Beijing, China

Cong Peng

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

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

Finovation in CCBFT, Beijing, China

## Abstract

Pseudo-Boolean (PB) constraints are highly expressive, and many combinatorial optimization problems can be modeled using pseudo-Boolean optimization (PBO). It is recognized that stochastic local search (SLS) is a powerful paradigm for solving combinatorial optimization problems, but the development of SLS for solving PBO is still in its infancy. In this paper, we develop an effective SLS algorithm for solving PBO, dubbed *NuPBO*, which introduces a novel scoring function for PB constraints and a new weighting scheme. We conduct experiments on a broad range of six public benchmarks, including three real-world benchmarks, a benchmark from PB competition, an integer linear programming optimization benchmark, and a crafted combinatorial benchmark, to compare *NuPBO* against five state-of-the-art competitors, including a recently-proposed SLS PBO solver *LS-PBO*, two complete PB solvers *PBO-IHS* and *RoundingSat*, and two mixed integer programming (MIP) solvers *Gurobi* and *SCIP*. *NuPBO* has been exhibited to perform best on these three real-world benchmarks. On the other three benchmarks, *NuPBO* shows competitive performance compared to state-of-the-art competitors, and it significantly outperforms *LS-PBO*, indicating that *NuPBO* greatly advances the state of the art in SLS for solving PBO.

**2012 ACM Subject Classification** Theory of computation → Randomized local search

**Keywords and phrases** Pseudo-Boolean Optimization, Stochastic Local Search, Scoring Function, Weighting Scheme

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

**Supplementary Material** Software (Source Code): https://github.com/filyouzicha/NuPBO

**Funding** This work is supported by the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDA0320000 and XDA0320300, the National Natural Science Foundation of China under Grant 62302492, and Grant 62202025, the Research Program of CCBFT, Grant No. KT2100040, and CCF-Huawei Populus Grove Fund under Grant CCF-HuaweiSY20231.

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

<footnote n="1">¹ Corresponding author.</footnote>
# 1 Introduction

Recent progress in the theory and application of the Boolean satisfiability (SAT) and maximum satisfiability (MaxSAT) problems has led to the development of high-performance complete solvers [11, 23, 1, 32, 2, 4, 31] and stochastic local search (SLS) solvers [8, 28, 7, 6, 27, 25, 5, 45]. SAT and MaxSAT solvers can address challenging problems in a wide variety of fields, and they are usually designed to deal with formulas encoded in conjunctive normal form (CNF).

For problems involving cardinality constraints, CNF solvers usually become ineffective, since expressing such constraints in CNF would dramatically increase the size of the formula and introduce many auxiliary variables and clauses [30]. Linear pseudo-Boolean (PB) constraints provide a more natural and direct way to express cardinality constraints than CNF. Meanwhile, linear PB constraints stay close to CNF and can benefit from advancements in SAT solving [33]. In practice, PB constraints occur in many areas, including VLSI design, economics, computer vision, and manufacturing [43, 44, 33]. The pseudo-Boolean optimization (PBO) problem is to find a satisfying assignment to a set of PB constraints that minimizes a given objective function.

## 1.1 Related Work

Existing pseudo-Boolean solvers are primarily based on complete methods. A number of PB solvers are based on resolution: they express the PB constraints in CNF and then call conflict-driven clause learning (CDCL) solvers, such as *MINISAT*+[13], *Open-WBO*[29], and *NaPS* [34]; alternatively, they deal with the PB constraints but derive new information only in the form of clauses [17]. CDCL is somewhat limited in its reasoning in that it is based on a resolution-proof system, for which exponential lower bounds are known for simple combinatorial principles [20, 15]. Another method requires going beyond resolution and using cutting planes, which can be found in recent PB solvers such as *Sat4j* [24], *RoundingSat* [14, 12] and *RoundingSat-Card* [15]. The success of the implicit hitting set (IHS) method in MaxSAT motivates another work, i.e., the implementation of the *PBO-IHS* solver for solving PBO [36, 37]. In addition, since PB constraints can be considered as 0-1 linear constraints, mixed integer programming (MIP) solvers can be directly applied to solving PBO. Representative and high-performance MIP solvers include *SCIP* [16] - one of the fastest non-commercial solvers, and *Gurobi* [19] - one of the most powerful commercial solvers.

Stochastic local search (SLS) is recognized to be one of the most powerful techniques for solving computationally hard problems in many areas of computer science, operations research, and engineering, and it has shown great success in solving SAT and MaxSAT [22]. In the book [21], a model for local search to solve constraint problems is presented. Somewhat surprisingly, there are only a few research works on using SLS for solving PBO [3, 39, 26].

Since the introduction of dynamic local search (DLS) methods [35, 9, 40, 41], weighting schemes play critical roles in the development of high-performance SLS algorithms. Modern SLS solvers for MaxSAT employ a scoring function defined as the weighted cost of unsatisfied clauses and incorporate a clause weighting scheme to adjust the weights during the search. In particular, most SLS solvers for MaxSAT focus on improving the scoring function through carefully designed weighting schemes. In recent years, the introduction of new weighting schemes has led to breakthroughs in the SLS algorithms for the (weighted) partial MaxSAT ((W)PMS) problem. The newly proposed weighting scheme, Weighting-PMS, in the SATLike [25] algorithm significantly improves the performance of SLS for (W)PMS. Currently, the state-of-the-art SLS algorithms for (W)PMS all employ the Weighting-PMS technique [5, 45].
It is intuitive that the scoring function commonly used in SLS algorithms for MaxSAT (i.e., measuring the total weight of all unsatisfied clauses) does not work well for the PBO problem, because it does not take into account the unsatisfied degree of PB constraints. Indeed, this issue has already attracted attention in the literature, and a scoring function that considers the unsatisfied degree of PB constraints is proposed in [42] and can effectively handle linear PB-constrained problems.

A recent SLS solver called *LS-PBO* [26] has been proposed for PBO and currently represents the state of the art in SLS for PBO. For *LS-PBO*, its scoring function measures the sum of the product between the degree of violation of all unsatisfied constraints and the weights of the constraints. However, such scoring function does not consider the balance between the degree of violation among different constraints. In addition, weighting schemes have not been applied to PBO until the introduction of *LS-PBO*. The weighting scheme in *LS-PBO* resembles the existing one named Weighting-PMS, which aims to increase the weights of unsatisfied constraints and the objective function when the algorithm falls into a local optimum, and to set an upper bound on the weight of the objective function. However, in the context of PBO solving, the research of designing weighting schemes is still in its infancy, which urgently calls for more powerful weighting schemes for PBO.

## 1.2 Contributions

In this work, we focus on improving the performance of SLS for solving PBO. In particular, we propose two main ideas. The first idea is a novel scoring function that considers the violation degree of unsatisfied constraints and utilizes a smooth function to balance the violation degree of different constraints. For each constraint, its smooth function is instantiated as the average of the coefficients of all variables appearing in the corresponding constraint. Since our scoring function is equipped with a weighting scheme, our second idea is a novel weighting scheme for PBO. Rather than setting an upper bound on the weight of the objective function, we adopt a weighting scheme with a stricter condition for updating the weight of the objective function.

On the basis of these two ideas, we develop a new SLS algorithm, named *NuPBO*. We conduct experiments on 6 benchmarks, which include 3 benchmarks encoded from real-world applications, and 3 standard benchmarks. On these 6 benchmarks, *NuPBO* is compared to 5 solvers, including *LS-PBO* [26], *PBO-IHS* [37], *RoundingSat* [12], *Gurobi* [19], and *SCIP* [16]. On the 3 application benchmarks, *NuPBO* achieves improvement over *LS-PBO*, and significantly outperforms other competitors. On the other 3 benchmarks, *NuPBO* exhibits competitive performance compared to its competitors, including the commercial solver *Gurobi*.
This represents a significant advance in the research of SLS solvers for PBO. In addition, we evaluate the effectiveness of the underlying ideas on all benchmarks.

## 2 Preliminaries

Given a set of $n$ Boolean variables $V = \{x_1, x_2, \dots, x_n\}$, a literal is either a variable $x_i$ or its negation $\neg x_i$. A clause is a disjunction of literals, i.e., $c_i = l_{i1} \lor l_{i2} \lor \dots \lor l_{ik}$, where $k$ denotes the length of clause $c_i$. A CNF formula $F$ is a conjunction of clauses. An assignment is a mapping that assigns a Boolean value (*True* (i.e., 1) or *False* (i.e., 0)) to each variable. Given an assignment $\alpha$, a clause $c$ is satisfied if at least one literal in $c$ is *True*; otherwise, $c$ is unsatisfied. Given a CNF formula $F$, the Boolean satisfiability (SAT) problem is to decide whether an assignment exists such that all clauses are satisfied, and the maximum satisfiability (MaxSAT) problem is to find an assignment that maximizes the number of satisfied clauses.
12:4 Towards More Efficient Local Search for Pseudo-Boolean Optimization

A linear pseudo-Boolean constraint (LPB constraint, PB constraint for short) has the following form:

$$ \sum_{j=1}^{n} a_{j} l_{j} \triangleright b, \quad a_{j}, b \in \mathbb{Z} \qquad (1) $$

where b is called the degree of the constraint, $l_j$ is a literal, $a_j$ is the coefficient of $l_j$, n is the length of the constraint, $\triangleright$ is one of the classical relational operators (>=, >, $\ge$, < or $\le$), and $\mathbb{Z}$ is the integer set.

For each Boolean variable, $x_i = 1 - \neg x_i$. It is important to note that for an equality constraint, there need to be two normalized constraints to represent it. Therefore, all PB constraints can be normalized into the following form:

$$ \sum_{j=1}^{n} a_{j} l_{j} \geq b, \quad a_{j}, b \in \mathbb{N}_{0}^{+} \qquad (2) $$

where $\mathbb{N}_0^+$ is the non-negative integer set [33].

In the following sections, we assume that the PB constraints are of the normalized form. A PB formula $F$ is a conjunction of PB constraints. An assignment is a mapping that assigns a Boolean value to each variable. Given an assignment $\alpha$, a PB constraint $c$ is satisfied if the corresponding inequality holds under $\alpha$; otherwise, $c$ is unsatisfied. If an assignment $\alpha$ satisfies all constraints in $F$, then we say $\alpha$ is a feasible solution (or solution for short).

A pseudo-Boolean optimization (PBO) instance consists of a PB formula $F$ and a linear Boolean objective function $\sum_{j=1}^{n} e_{j} l_{j} + d$, $e_j \in \mathbb{N}^+$, $d \in \mathbb{Z}$, and the task is to find an assignment that satisfies all PB constraints in $F$ and minimizes the objective function. Given an assignment $\alpha$, we use $obj(\alpha)$ to denote the value of the objective function. Given a solution $\alpha$, the cost of the solution $\alpha$ is equal to $obj(\alpha)$. We say a solution $\alpha_1$ is better than another solution $\alpha_2$, if $obj(\alpha_1) < obj(\alpha_2)$.

The average coefficient of a PB constraint $c$ is denoted as $avg_{coe}(c) = (\sum_{j=1}^{n} a_j)/n$. The average coefficient of an objective function $o$ is denoted as $avg_{coe}(o) = (\sum_{j=1}^{n} e_j)/n$. Because PB constraints must be satisfied in a PBO problem, the PB constraints are referred to as hard constraints.

Given an assignment $\alpha$, we define the value of violation of a hard constraint $c$ as

$$ \mathrm{viol}(c) = \max \left( 0, b - \sum_{j=1}^{n} a_{j} l_{j} \right) $$

In other words, if the hard constraint $c$ is satisfied under $\alpha$, then $viol(c) = 0$; otherwise (i.e., $c$ is unsatisfied), $viol(c)$ is the integer distance of $c$ from being satisfied. In a PBO instance, a solution is an assignment, under which all hard constraints are satisfied.

The SLS algorithms that employ constraint weighting schemes usually maintain weight for each constraint. We use $w(c)$ to represent the weight of each hard constraint $c$, and $w(o)$ to represent the weight of the objective function $o$.

3 Main Ideas

In general, the search directions of SLS algorithms are guided by the scoring function. It is recognized that the effectiveness of the scoring function could be enhanced through working with a weighting scheme. In this section, we first propose a new scoring function, and then design a new weighting scheme to work with it.
<table><caption>**Table 1** The violation and objective value under all assignments of the PBO instance $I_1$ in Example 3.1.</caption><thead><tr><th rowspan="2"></th><th colspan="7">Assignment (x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub>)</th></tr><tr><th>(0, 0, 0)</th><th>(0, 0, 1)</th><th>(0, 1, 0)</th><th>(0, 1, 1)</th><th>(1, 0, 0)</th><th>(1, 0, 1)</th><th>(1, 1, 0)</th><th>(1, 1, 1)</th></tr></thead><tbody><tr><th>viol(obj)</th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>viol(c<sub>1</sub>)</th><td>0</td><td>0</td><td>0</td><td>0</td><td>4</td><td>3</td><td>3</td><td>2</td></tr><tr><th>viol(c<sub>2</sub>)</th><td>1</td><td>0</td><td>0</td><td>0</td><td>4</td><td>2</td><td>3</td><td>1</td></tr><tr><th>viol(c<sub>3</sub>)</th><td>29</td><td>13</td><td>14</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><th>obj(α)</th><td>1</td><td>1</td><td>0</td><td>0</td><td>2</td><td>2</td><td>1</td><td>1</td></tr></tbody></table>

## 3.1 A New Scoring Function

Given a PBO instance, which consists of $n$ PB constraints (hard constraints) and one objective function, we assume that the current assignment is $\alpha$.

**Scoring Function in LS-PBO.** Before introducing our new scoring function, we first describe the existing scoring function proposed in LS-PBO [26], which is presented as follows.

- For a hard constraint $c$, the penalty function is defined as $\text{penalty}(c) = w(c) \times \text{viol}(c)$; then the hard score of a variable $x$ is defined as the decrement of the sum of the penalty function values of all hard constraints caused by flipping $x$, which is denoted by $\text{hscore}(x)$.

- For the objective function $o$, the value of the objective function is $\text{obj}(\alpha)$, and the penalty function is defined as $\text{penalty}(o) = w(o) \times \text{obj}(\alpha)$; then the soft score of a variable $x$ is defined as the decrement of the penalty function value of the objective function caused by flipping $x$, which is denoted by $\text{sscore}(x)$.

- The score of a variable $x$ is defined as $\text{score}(x) = \text{hscore}(x) + \text{sscore}(x)$.

**An Intuitive View.** The above scoring function has a drawback, which is due to its underlying penalty function. The penalty function above considers the weights and *viol* (resp. *obj*) values of hard (resp. soft) constraints. In this way, it measures the importance of a variable in a constraint $c$ by the coefficient of the corresponding variable in $c$. Nevertheless, it should be noted that, as the value of $\text{score}(x)$ is determined by all hard constraints and the objective function involving $x$, the above penalty function may overemphasize the importance of $x$ in constraints with relatively large coefficients, resulting in an unreasonable value of its score. In the following, we illustrate our intuition through a simple PBO instance.

**► Example 3.1.** Let us consider a PBO instance $I_1$, which consists of three hard constraints and an objective function: $c_1: 4 \neg x_1 + x_2 + x_3 \ge 4$, $c_2: 3 \neg x_1 + x_2 + 2x_3 \ge 4$, $c_3: 29x_1 + 15x_2 + 16x_3 \ge 29$, minimize: $x_1 + \neg x_2$. The values of *viol* of hard constraints and the value of *obj* of the objective function under all assignments are presented in Table 1. Solutions for $I_1$ are those resulting in the zero value of *viol* for all hard constraints. From Table 1, the optimal solution is $\alpha^* = (0, 1, 1)$, and $\text{obj}(\alpha^*)$ is 0.

Given instance $I_1$, consider a scoring function without any weighting scheme, or equivalently, the weight of each constraint is 1, i.e., $w(c_1) = 1$, $w(c_2) = 1$, $w(c_3) = 1$, $w(o) = 1$; the initial assignment $\alpha = (0, 0, 0)$. In accordance with the definition of the scoring function in LS-PBO, $\text{score}(x_1) = 21$, $\text{score}(x_2) = 17$, and $\text{score}(x_3) = 17$. Actually, in order to optimize the assignment, SLS algorithms tend to select the variable to be flipped as the one with the largest score, so in this situation, $x_1$ is picked. After flipping $x_1$, the assignment becomes $\alpha = (1, 0, 0)$, and the score value of each variable becomes $\text{score}(x_1) = -21$, $\text{score}(x_2) = 3$, $\text{score}(x_3) = 3$. Then, no matter whether $x_2$ or $x_3$ is flipped, the Hamming distance between
the current assignment and the optimal solution is the same as that between the initial assignment and the optimal solution, which is two. The search is not progressing in the direction towards the optimal solution.

In practice, PBO instances encoded from real-world problems are much more complex than the given illustrative example. If SLS algorithms conduct the search in incorrect directions, it would be difficult to identify a promising search space that is more likely to contain the optimal solution or those close to optimality.

As presented in Table 1, when we focus on the value of *viol(c<sub>3</sub>)*, for those cases where the value of *viol(c<sub>3</sub>)* is not 0, its value is much larger than the *viol* value of other hard constraints. Considering that each hard constraint has a *penalty* value directly proportional to its *viol* value, utilizing scoring function aims to guide the search towards the area with a lower sum of *penalty* values. Consequently, through making use of such scoring function, the algorithms would prefer the falsified literal with the largest coefficient to be *True* (in instance *I*<sub>1</sub>, under the assumption that the current assignment is (0, 0, 0), this falsified literal is *x*<sub>1</sub> in the hard constraint *c*<sub>3</sub>).

**Our New Scoring Function.** In our opinion, a good scoring function for PBO should balance the *viol* values of different constraints. To this end, we propose to smooth the *penalty* values of constraints. For simplicity, we denote the smoothing function of a hard constraint *c* as *smooth(c)*, and the smoothing function of the objective function *o* as *smooth(o)*. Based on the idea of balancing the *viol* value, we propose the following, new scoring function:

- For a hard constraint *c*, the penalty function is defined as $\text{penalty}(c) = \frac{w(c) \times \text{viol}(c)}{\text{smooth}(c)}$; then the hard score of a variable *x* is defined as the decrement of the sum of the penalty function values of all hard constraints caused by flipping *x*, which is denoted by $\text{hscore}(x)$.

- For the objective function *o*, the value of the objective function is $\text{obj}(\alpha)$, and the penalty function is defined as $\text{penalty}(o) = \frac{w(o) \times \text{obj}(\alpha)}{\text{smooth}(o)}$; then the soft score of a variable *x* is defined as the decrement of the penalty function value of the objective function caused by flipping *x*, which is denoted by $\text{sscore}(x)$.

- The score of a variable *x* is defined as $\text{score}(x) = \text{hscore}(x) + \text{sscore}(x)$.

In order to instantiate the above scoring function, we propose to use a method for smoothing by using the average of the constraint coefficients, i.e., $\text{smooth}(c) = \text{round}(\text{avg}_{\text{coe}}(c))$, $\text{smooth}(o) = \text{round}(\text{avg}_{\text{coe}}(o))$ (round is a rounding function). Consider the PBO instance $I_1$ in Example 3.1, which has $\text{smooth}(c_1) = 2$, $\text{smooth}(c_2) = 2$, $\text{smooth}(c_3) = 20$, and $\text{smooth}(o) = 1$. Assume that each hard constraint and the objective function have a weight of 1 and the current assignment $\alpha = (0,0,0)$. Based on the new scoring function, $\text{score}(x_1) = -3.05$, $\text{score}(x_2) = 2.25$, and $\text{score}(x_3) = 1.3$. Hence, SLS algorithms would select variable $x_2$ to be flipped. Flipping $x_2$ would change the current assignment $\alpha$ to (0,1,0), and the score value of each variable would become $\text{score}(x_1) = -3.3$, $\text{score}(x_2) = -2.25$, $\text{score}(x_3) = 0.7$. Afterward, SLS algorithms would select variable $x_3$ to be flipped. If $x_3$ is flipped, assignment $\alpha$ becomes (0, 1, 1), which is the optimal solution of instance $I_1$.

According to this illustrative example, it can be observed that, by using the average of constraint coefficients to smooth the *penalty* value, the issue of the large difference among coefficients of a variable in various constraints can be alleviated, resulting in a more effective scoring function.
## 3.2 A New Weighting Scheme

Combinatorial optimization problems with both hard and soft constraints require effective weighting schemes that balance the weights of hard and soft constraints. A potential problem was pointed out in a previous study [10]: the excessive weight given to soft constraints may make it difficult to satisfy all the hard constraints, thereby hindering the algorithm's capability of finding solutions. Moreover, an existing study [38] demonstrates that designing a weighting scheme for problems with hard constraints is challenging as it requires weighting unsatisfied constraints while maintaining the distinction between hard and soft constraints.

To alleviate the above problem, the weighting scheme proposed in LS-PBO sets an upper bound $\zeta$ (an integer parameter) to the maximum value of the objective function weight. We use `unsat_hard_set` to denote the set of unsatisfied hard constraints. For a PBO instance $I$, the average of the product of the $avg_{coe}(c)$ and $w(c)$ of all constraints $c$ is denoted as $wavg_{coe}(I)$, that is, $wavg_{coe}(I) = (\sum_{i=1}^{m}(avg_{coe}(c_i) \times w(c_i))) / m$. Assuming that the current assignment is $\alpha$, the best solution that has been found is $\alpha^*$, and its corresponding objective function value is $obj(\alpha^*)$.

The weighting scheme adopted in LS-PBO is described as follows:

- Initialization phase: at the start of the local search process, the weight of each hard constraint $c$ is initialized as 1, i.e., $w(c) := 1$; the weight of the objective function $o$ is also initialized as 1, i.e., $w(o) := 1$.

- Update phase: when the search is trapped in a local optimum (i.e., there is no variable whose score value is greater than 0), for each $c$ in `unsat_hard_set`, $w(c) := w(c) + 1$; if $obj(\alpha) \ge obj(\alpha^*)$ and $w(o) \times avg_{coe}(o) - wavg_{coe}(I) \le \zeta$, $w(o) := w(o) + 1$, where $\zeta$ is a parameter introduced by LS-PBO.

In fact, the setting of $\zeta$ greatly affects the performance of LS-PBO, and if the average of the coefficients of the objective function is much greater than the average of the coefficients of the hard constraints, the weight of the objective function basically would not be updated. In addition, varying the timing of weighting constraints could be a promising strategy to improve the performance of the weighting scheme.

We propose to deal with these problems by modifying the condition of updating weights. Specifically, we propose a stricter condition for increasing objective function weight. Our proposed weighting scheme is as follows:

- Initialization phase: at the start of the local search process, the weight of each hard constraint $c$ is initialized as 1, i.e., $w(c) := 1$; **the weight of the objective function $o$ is initialized as 0**, i.e., $w(o) := 0$.

- Update phase: when the search is trapped in a local optimum (i.e., there is no variable whose score value is greater than 0), for each $c$ in `unsat_hard_set`, $w(c) := w(c) + 1$; **if `unsat_hard_set` is empty**, $w(o) := w(o) + 1$.

In the beginning, the weight of the objective function is initialized as 0, so that the algorithm would first focus on finding solutions. If the search is trapped in a local optimum, the weight of the objective function is increased only when the current assignment $\alpha$ is a solution (all hard constraints are satisfied under $\alpha$). Accordingly, if the algorithm frequently visits solutions, then the objective function would have a greater chance to increase its weight. Otherwise, there would be limited opportunities to increase the objective function weight.
**Algorithm 1** The NuPBO Algorithm.

**Input:** A PBO instance *I*, cutoff time.

**Output:** The best solution ($\alpha^*$) found and its objective function value $obj^*$, or "No solution found".

1 $\alpha^* := \emptyset$; $obj^* := +\infty$;

2 **while** no terminating criteria are met **do**

3 &nbsp;&nbsp;&nbsp;&nbsp; $\alpha :=$ an initial assignment;

4 &nbsp;&nbsp;&nbsp;&nbsp; **for** each *c* in hard constraints **do**

5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor w(c) \rfloor := 1$;

6 &nbsp;&nbsp;&nbsp;&nbsp; for the objective function *o*, $w(o) := 0$;

7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $L = 10000000$;

8 &nbsp;&nbsp;&nbsp;&nbsp; **for** step=0; step<L; step++ **do**

9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $\alpha$ is feasible and $obj^* > obj(\alpha)$ **then**
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor \alpha^* := \alpha; obj^* := obj(\alpha); L = step + 10000000 \rfloor$;

11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $D := \{x | score(x) > 0\} \neq \emptyset$ **then**

12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor v := \text{a variable in } D \text{ with the highest } score \rfloor$;

13 &nbsp;&nbsp;&nbsp;&nbsp; **else**
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; update constraints weights by the new weighting scheme described in Section 3.2;

15 &nbsp;&nbsp;&nbsp;&nbsp; **if** $\exists$ unsatisfied hard constraints **then**

16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor c := \text{a random unsatisfied hard constraint} \rfloor$;

17 &nbsp;&nbsp;&nbsp;&nbsp; $\lfloor v := \text{the variable whose literal is false with highest } score \text{ in } c \rfloor$;

18 &nbsp;&nbsp;&nbsp;&nbsp; **else**
19 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\lfloor v := \text{a randomly chosen variable with } sscore > 0 \rfloor$;

20 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\alpha := \alpha$ with *v* flipped;

21 **if** $\alpha^* \neq \emptyset$ **then** return $\alpha^*$ and $obj^*$;

22 **else** return No solution found;

## 4 The NuPBO Algorithm

In this section, we develop a new SLS algorithm named NuPBO, which is based on the main ideas proposed in Section 3. NuPBO adopts the Dynamic Local Search (DLS) framework as does LS-PBO. The pseudo-code of NuPBO is outlined in Algorithm 1. We use $\alpha^*$ and $obj^*$ to denote the best-found solution and the corresponding objective function value (i.e., the cost of the best-found solution), while $\alpha$ denotes the current assignment which is maintained during the search.

In the beginning, $\alpha^*$ is initialized as an empty set, and $obj^*$ is initialized as $+\infty$ (Line 1). NuPBO then iteratively calls the local search process until reaching a terminating criterion (e.g., reaching a preset cutoff time, or achieving a feasible assignment $\alpha$ whose corresponding $obj$ value is equal to 0) (Lines 2-20).

In the local search process, an initial assignment is generated by assigning each Boolean variable to a default value 0 (as in LS-PBO) (Line 3). NuPBO then initializes the weights of all hard constraints as 1, and sets the weight of objective function to 0 according to our proposed weighting scheme. After initialization, NuPBO conducts the search process (Lines 8-20). During the search process, whenever NuPBO finds a solution whose $obj$ is lower than $obj^*$, then $\alpha^*$ and $obj^*$ are updated accordingly.

In each search step, NuPBO selects a variable and flips it based on two situations: (I) If the set *D* of decreasing variables (i.e., the variable *x* with $score(x) > 0$) is not empty, a variable with the highest score is selected from *D*, breaking ties by preferring the variable
that has been flipped least recently. (II) When *D* is empty, which indicates that the search is trapped in a local optimum, then *NuPBO* updates the weights of constraints according to our new weighting scheme. Then, if there exist unsatisfied hard constraints, an unsatisfied hard constraint *c* is randomly picked. As we know, for a MaxSAT instance, if a clause *cm* is unsatisfied, then it would become satisfied after flipping any variables in *cm*. For a PBO instance, if a hard constraint *cp* is unsatisfied (i.e., *violcp*) > 0), only flipping variables whose literals are *False* under the current assignment can reduce *violcp* (Assuming under the current assignment *α*, *x*₁ is 1, then the literal *x*₁ is *True*, while the literal *¬x*₁ is *False*). Therefore, *NuPBO* picks the variable whose literal is *False* with the highest *score* in *c*. Otherwise (i.e., all the hard constraints are satisfied), *NuPBO* randomly chooses a variable whose *sscore* is greater than 0.

Finally, when any terminating criterion is met, *NuPBO* stops and reports the best solution *α*\* and *obj*\* if a solution is found; otherwise, it reports "No solution found".

# 5 Experimental Evaluations

In this section, we introduce experimental preliminaries and then conduct extensive experiments on 6 PBO benchmarks. First, we compare *NuPBO* with 5 state-of-the-art PBO solvers. Second, we conduct experiments to show that combining *NuPBO* with complete solvers can lead to better portfolios. Third, we report experimental results to demonstrate the effectiveness of our main ideas. Finally, we examine the stability of the SLS solvers by running each SLS solver 10 times with seeds ranging from 1 to 10.

## 5.1 Experimental Preliminaries

**Benchmarks.** We evaluate *NuPBO* on 6 benchmarks, which are described as follows:

- PB16: the OPT-SMALLINT-LIN benchmark from the latest 2016 pseudo-Boolean competition. As a mainstream benchmark for evaluating the performance of PBO solvers, it consists of 1600 instances of various categories.<fnref n="2" />

- MIPLIB: 0-1 integer linear programming optimization problems. This benchmark contains 291 instances of various categories, provided in the literature [12].<fnref n="3" />

- CRAFT: crafted combinatorial benchmarks whose coefficients are small integers. This benchmark contains 955 instances of various categories, provided in the literature [12].<fnref n="4" />

- MWCB: the Minimum-Width Confidence Band Problem. This benchmark contains 24 instances.

- SAP: the Seating Arrangements Problem. This benchmark contains 21 instances.

- WSNO: the Wireless Sensor Network Optimization Problem. This benchmark consists of 18 instances.

For the benchmarks of MWCB, SAP, and WSNO, the descriptions, the downloading websites, and the methods of converting the real-world applications into PBO instances and the encoded PBO instances are presented in the literature [26].<fnref n="5" />

<footnote n="2">² http://www.cril.univ-artois.fr/PB16/bench/PB16-used.tar</footnote>

<footnote n="3">³ https://zenodo.org/record/3870965</footnote>

<footnote n="4">⁴ https://zenodo.org/record/4036016</footnote>

<footnote n="5">⁵ https://lcs.ios.ac.cn/%7ecaisw/Resource/LS-PBO/</footnote>
12:10 Towards More Efficient Local Search for Pseudo-Boolean Optimization

**State-of-the-Art Competitors.** We compare NuPBO with 5 state-of-the-art solvers, including one SLS solver (i.e., LS-PBO) and 4 complete solvers. The 4 complete solvers include 2 PB solvers (i.e., PBO-IHS and RoundingSat) and 2 MIP solvers (i.e., Gurobi and SCIP):

- LS-PBO [26]: the state-of-the-art SLS algorithm for solving PBO. Adopt the parameter setting recommended by its authors. It outperforms Gurobi and RoundingSat on many real-world application benchmarks.⁵

- PBO-IHS [37]: a recent IHS PBO solver building upon RoundingSat [14].<fnref n="6" />

- RoundingSat [12]: a recent PBO solver combining core-guided search with cutting planes reasoning.<fnref n="7" />

- Gurobi [19]: one of the most powerful commercial MIP solvers (Version 9.1.2). The default configuration is used, along with a single thread.<fnref n="8" />

- SCIP [16]: one of the fastest non-commercial solvers for MIP (Version 7.0.3, using SoPlex 5.0.2 as its internal LP solver).<fnref n="9" />

**Experimental Setup.** LS-PBO and NuPBO are implemented in C++, and compiled with g++ (version 8.5.0) using the option "-O3". Installation procedures for other solvers follow their detailed guidelines. All the experiments are carried out on a workstation under the operating system CentOS, with the AMD EPYC7702 2.0GHz CPU.

In these experiments, we adopt two cutoff times of 300 CPU seconds (300s) and 3600 CPU seconds (1h). Each solver performs one run within a given cutoff time on each instance, and we record the cost of the best solution found by solver $S_j$ on instance $I_k$, denoted as $sol_{S_jI_k}$. The cost of the best solution found among all solvers in the same table within the same cutoff time on instance $I_k$ is denoted as $best_{I_k}$. For each solver $S$ solving a benchmark $B_i$ within a cutoff time, we use 3 *metrics* to evaluate the performance of $S$.

- \#win.: the number of instances where the corresponding $best_{I_k}$ can be obtained by solver $S$ on $B_i$ (i.e., the number of winning instances).

- avg\_score: in our experiments, the competition score of solver $S_j$ on instance $I_k$ is represented by $score_{S_jI_k} = \frac{best_{I_k}+1}{sol_{S_jI_k}+1}$, which measures the gap between $sol_{S_jI_k}$ and $best_{I_k}$. If solver $S_j$ could not report a solution on instance $I_k$, then $score_{S_jI_k} = 0$. We use avg\_score to denote the average competition score of a solver on a benchmark. The competition score of each solver on each instance is the metric to measure the performance of solvers in the incomplete track of recent MaxSAT Evaluations (2017-2023).

- \#feas.: the number of instances where solver $S$ obtains solutions on $B_i$.

In our experiments, avg\_score is calculated by ignoring the instances that are proven to have no solution by the complete solvers. Based on the preliminary experiments, we conclude that at least 123 instances in the PB16 benchmark and at least 17 instances in the MIPLIB benchmark do not have solutions. All instances in the CRAFT, MWCB, SAP, and WSNO benchmarks have solutions.<fnref n="10" />

The number of instances in each benchmark is indicated by \#inst.. For each of the above three metrics, if a solver obtains a larger metric value on a benchmark, then the solver exhibits better performance on the benchmark. The results highlighted in **bold** indicate the best performance for the corresponding metric.

<footnote n="6">⁶ https://bitbucket.org/coreo-group/pbo-ihs-solver/</footnote>

<footnote n="7">⁷ https://doi.org/10.5281/zenodo.4043124</footnote>

<footnote n="8">⁸ https://www.gurobi.com/products/gurobi-optimizer/</footnote>

<footnote n="9">⁹ https://www.scipopt.org/index.php#download</footnote>

<footnote n="10">¹⁰ Note that the definition of the competition score (metric) in this subsection has no relationship to the definition of the score in the scoring function in subsection 3.1.</footnote>
<table><caption>Table 2 Experimental results of NuPBO and all the competitors on all the benchmarks (top: cutoff 300s, bottom: cutoff 1h) (Benc, i.e., Benchmark. avg<sub>s</sub>, i.e., avg<sub>score</sub>).</caption><thead><tr><th rowspan="3">Benc</th><th rowspan="3">#inst.</th><th colspan="2">SLS solvers</th><th colspan="2">PB solvers</th><th colspan="2">MIP solvers</th></tr><tr><th>NuPBO</th><th>LS-PBO</th><th>PBO-IHS</th><th>RoundingSat</th><th>Gurobi</th><th>SCIP</th></tr><tr><th>#win.<br>avg.<br>#feas.</th><th>#win.<br>avg.<br>#feas.</th><th>#win.<br>avg.<br>#feas.</th><th>#win.<br>avg.<br>#feas.</th><th>#win.<br>avg.<br>#feas.</th><th>#win.<br>avg.<br>#feas.</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">cutoff 300s</td></tr><tr><td>PB16</td><td>1600</td><td>1049<br>1351</td><td>0.8800</td><td>700<br>1183</td><td>0.7516</td><td>930<br>1401</td><td>0.8712</td><td>964<br>1392</td><td>0.8916</td><td><b>1158</b><br>1365</td><td><b>0.9011</b></td><td>887<br>1244</td><td>0.7918</td></tr><tr><td>MIPLIB</td><td>291</td><td>130<br>242</td><td><b>0.8480</b></td><td>75<br>222</td><td>0.7375</td><td>84<br>230</td><td>0.7350</td><td>89<br>245</td><td>0.7855</td><td><b>166</b><br>235</td><td>0.8013</td><td>114<br>221</td><td>0.7360</td></tr><tr><td>CRAFT</td><td>955</td><td><b>894</b><br>943</td><td>0.9868</td><td>816<br>930</td><td>0.9682</td><td>809<br>930</td><td>0.9433</td><td>825<br>936</td><td>0.9639</td><td>893<br>955</td><td><b>0.9961</b></td><td>764<br>921</td><td>0.9583</td></tr><tr><td>MWCB</td><td>24</td><td><b>24</b><br>24</td><td><b>1.0000</b></td><td>0<br>24</td><td>0.9448</td><td>0<br>19</td><td>0.4496</td><td>0<br>24</td><td>0.6247</td><td>0<br>12</td><td>0.4004</td><td>0<br>3</td><td>0.0761</td></tr><tr><td>SAP</td><td>21</td><td><b>21</b><br>21</td><td><b>1.0000</b></td><td>0<br>21</td><td>0.9750</td><td>0<br>0</td><td>0.0000</td><td>0<br>0</td><td>0.0000</td><td>0<br>1</td><td>0.0395</td><td>0<br>0</td><td>0.0000</td></tr><tr><td>WSNO</td><td>18</td><td><b>18</b><br>18</td><td><b>1.0000</b></td><td>11<br>18</td><td>0.9026</td><td>2<br>5</td><td>0.1738</td><td>4<br>18</td><td>0.6624</td><td>4<br>5</td><td>0.2431</td><td>0<br>6</td><td>0.1631</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">cutoff 1h</td></tr><tr><td>PB16</td><td>1600</td><td>1064<br>1360</td><td>0.8897</td><td>814<br>1278</td><td>0.8164</td><td>990<br>1412</td><td>0.8875</td><td>1012<br>1401</td><td>0.9100</td><td><b>1229</b><br>1396</td><td><b>0.9354</b></td><td>1012<br>1304</td><td>0.8565</td></tr><tr><td>MIPLIB</td><td>291</td><td>127<br>243</td><td>0.8519</td><td>83<br>231</td><td>0.7673</td><td>100<br>239</td><td>0.7903</td><td>87<br>247</td><td>0.8099</td><td><b>199</b><br>253</td><td><b>0.9023</b></td><td>142<br>238</td><td>0.8104</td></tr><tr><td>CRAFT</td><td>955</td><td>891<br>955</td><td><b>0.9992</b></td><td>844<br>932</td><td>0.9714</td><td>886<br>947</td><td>0.9841</td><td>873<br>955</td><td>0.9902</td><td><b>930</b><br>955</td><td>0.9987</td><td>824<br>930</td><td>0.9704</td></tr><tr><td>MWCB</td><td>24</td><td><b>23</b><br>24</td><td><b>0.9998</b></td><td>1<br>24</td><td>0.9690</td><td>0<br>24</td><td>0.5620</td><td>0<br>24</td><td>0.7116</td><td>0<br>24</td><td>0.7437</td><td>0<br>17</td><td>0.5058</td></tr><tr><td>SAP</td><td>21</td><td><b>21</b><br>21</td><td><b>1.0000</b></td><td>0<br>21</td><td>0.9785</td><td>0<br>0</td><td>0.0000</td><td>0<br>0</td><td>0.0000</td><td>0<br>1</td><td>0.0451</td><td>0<br>0</td><td>0.0000</td></tr><tr><td>WSNO</td><td>18</td><td><b>18</b><br>18</td><td><b>1.0000</b></td><td>15<br>18</td><td>0.9985</td><td>5<br>14</td><td>0.5989</td><td>11<br>18</td><td>0.8660</td><td>4<br>13</td><td>0.4904</td><td>0<br>4</td><td>0.0842</td></tr></tbody></table>

## 5.2 Comparisons with State-of-the-Art Solvers

The comparative results of NuPBO and all the competitors on all the benchmarks are shown in Table 2. We first analyze the results with a cutoff time of 300s.

- In terms of the number of winning instances, NuPBO gives the best performance on 4 benchmarks, including CRAFT and the 3 real-world application benchmarks, and ranked second on the PB16 and MIPLIB benchmarks (with Gurobi being the best).

- In terms of avg<sub>score</sub>, NuPBO outperforms all the competitors on 4 benchmarks, including MIPLIB and the 3 real-world application benchmarks. On the PB16 benchmark, its avg<sub>score</sub> ranks third after Gurobi and RoundingSat. The avg<sub>score</sub> value of NuPBO ranks second on the CRAFT benchmark behind Gurobi.

- In terms #feas., NuPBO and LS-PBO find solutions for all instances in the 3 real-world application benchmarks, while PBO-IHS, RoundingSat, and Gurobi are respectively the best for finding solutions on PB16, MIPLIB, and CRAFT. Although the value of #feas. of NuPBO ranks second on MIPLIB and CRAFT benchmarks, and ranks fourth on PB16 benchmark, NuPBO performs considerably better than LS-PBO, an SLS solver for PBO.
12:12 Towards More Efficient Local Search for Pseudo-Boolean Optimization

<table><caption>Table 3 Experimental results of VBS<sub>all</sub>, VBS<sub>exclude_lspbo</sub>, and VBS<sub>exclude_nupbo</sub> on all the benchmarks (top: cutoff 300s, bottom: cutoff 1h).</caption><thead><tr><th rowspan="2">Benchmark</th><th rowspan="2">#inst.</th><th colspan="2">VBS<sub>all</sub></th><th colspan="2">VBS<sub>exclude_lspbo</sub></th><th colspan="2">VBS<sub>exclude_nupbo</sub></th></tr><tr><th>#win.</th><th>avg<sub>score</sub></th><th>#win.</th><th>avg<sub>score</sub></th><th>#win.</th><th>avg<sub>score</sub></th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">cutoff 300s</td></tr><tr><td>PB16</td><td>1600</td><td>1434</td><td>0.9709</td><td>1434</td><td>1430</td><td>0.9703</td><td>1434</td></tr><tr><td>MIPLIB</td><td>291</td><td>254</td><td>0.9270</td><td>254</td><td>254</td><td>0.9270</td><td>254</td></tr><tr><td>CRAFT</td><td>955</td><td>955</td><td>1.0000</td><td>955</td><td>955</td><td>1.0000</td><td>955</td></tr><tr><td>MWCB</td><td>24</td><td>24</td><td>1.0000</td><td>24</td><td>24</td><td>1.0000</td><td>24</td></tr><tr><td>SAP</td><td>21</td><td>21</td><td>1.0000</td><td>21</td><td>21</td><td>1.0000</td><td>21</td></tr><tr><td>WSNO</td><td>18</td><td>18</td><td>1.0000</td><td>18</td><td>18</td><td>1.0000</td><td>18</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">cutoff 1h</td></tr><tr><td>PB16</td><td>1600</td><td>1440</td><td>0.9749</td><td>1440</td><td>1438</td><td>0.9747</td><td>1440</td></tr><tr><td>MIPLIB</td><td>291</td><td>262</td><td>0.9562</td><td>262</td><td>262</td><td>0.9562</td><td>262</td></tr><tr><td>CRAFT</td><td>955</td><td>955</td><td>1.0000</td><td>955</td><td>955</td><td>1.0000</td><td>955</td></tr><tr><td>MWCB</td><td>24</td><td>24</td><td>1.0000</td><td>24</td><td>23</td><td>0.9998</td><td>24</td></tr><tr><td>SAP</td><td>21</td><td>21</td><td>1.0000</td><td>21</td><td>21</td><td>1.0000</td><td>21</td></tr><tr><td>WSNO</td><td>18</td><td>18</td><td>1.0000</td><td>18</td><td>18</td><td>1.0000</td><td>18</td></tr></tbody></table>

With the cutoff time of 1h, NuPBO outperforms LS-PBO on the 3 real-world application benchmarks. On the other 3 benchmarks, NuPBO shows competitive performance compared to Gurobi, and significantly outperforms the state-of-the-art SLS solver LS-PBO in terms of all metrics of #win., avg<sub>score</sub>, and #feas..

5.3 Complementarity between SLS Solvers and Complete Solvers

In this subsection, we conduct experiments to investigate the complementarity between SLS solvers and complete solvers when solving PBO.

To investigate the complementarity between SLS solvers and complete solvers, we construct three perfect portfolio selectors: given a set of base solvers Θ, for each instance, the solution of the perfect portfolio selector constructed on Θ is the best among the entire collection of solutions reported by all solvers in Θ. These three perfect portfolio selectors are built based on Θ₁={LS-PBO, NuPBO, PBO-IHS, RoundingSat, Gurobi, SCIP}, Θ₂={NuPBO, PBO-IHS, RoundingSat, Gurobi, SCIP}, and Θ₃={LS-PBO, PBO-IHS, RoundingSat, Gurobi, SCIP} dubbed VBS<sub>all</sub>, VBS<sub>exclude_lspbo</sub> and VBS<sub>exclude_nupbo</sub>, respectively. Then we conduct experiments to evaluate the performance of these three perfect portfolio selectors on all benchmarks. The related results are presented in Table 3.

The comparison between VBS<sub>all</sub> and VBS<sub>exclude_lspbo</sub> reveals the number of instances where only the LS-PBO solver can achieve the optimal solution among all solvers. Similarly, comparing VBS<sub>all</sub> and VBS<sub>exclude_nupbo</sub> shows the number of instances where only the NuPBO solver can obtain the optimal solution among all solvers. As shown in Table 3, taking the PB16 benchmark with a cutoff time of 300 seconds as an example, out of the 1600 instances, there are 4 instances where only LS-PBO can achieve the optimal solution among all solvers, and there are 78 instances where only NuPBO can obtain the optimal solution among all solvers. Additionally, in 1 instance, only NuPBO was able to find a feasible solution. The results in Table 3 demonstrate that, compared to the state-of-the-art SLS solver LS-PBO, NuPBO is able to enhance the complementarity between SLS solvers and complete solvers, which indicates that a portfolio selector, which combines NuPBO and complete solvers, could advance the state of the art in PBO solving.
<table><caption>**Table 4** Experimental results of NuPBO, NuPBO-alt(s), and NuPBO-alt(w) on all the benchmarks (top: cutoff 300s, bottom: cutoff 1h).</caption><thead><tr><th rowspan="2">Benchmark</th><th rowspan="2">#inst.</th><th colspan="3">NuPBO</th><th colspan="3">NuPBO-alt(s)</th><th colspan="3">NuPBO-alt(w)</th></tr><tr><th>#win.</th><th>avg<sub>score</sub></th><th>#feas.</th><th>#win.</th><th>avg<sub>score</sub></th><th>#feas.</th><th>#win.</th><th>avg<sub>score</sub></th><th>#feas.</th></tr></thead><tbody><tr><td colspan="11" style="text-align:center; font-weight:bold;">cutoff 300s</td></tr><tr><td>PB16</td><td>1600</td><td><b>1141</b></td><td><b>0.9026</b></td><td><b>1351</b></td><td>1024</td><td>0.8288</td><td>1253</td><td>1034</td><td>0.8822</td><td>1323</td></tr><tr><td>MIPLIB</td><td>291</td><td><b>182</b></td><td><b>0.8687</b></td><td><b>242</b></td><td>147</td><td>0.8347</td><td><b>242</b></td><td>112</td><td>0.8312</td><td>239</td></tr><tr><td>CRAFT</td><td>955</td><td><b>941</b></td><td><b>0.9874</b></td><td><b>943</b></td><td><b>941</b></td><td>0.9864</td><td>942</td><td>848</td><td>0.9648</td><td>925</td></tr><tr><td>MWCB</td><td>24</td><td><b>24</b></td><td><b>1.0000</b></td><td>24</td><td>0</td><td>0.9466</td><td>24</td><td>0</td><td>0.9700</td><td>24</td></tr><tr><td>SAP</td><td>21</td><td><b>13</b></td><td><b>0.9990</b></td><td>21</td><td>9</td><td>0.9974</td><td>21</td><td>0</td><td>0.9782</td><td>21</td></tr><tr><td>WSNO</td><td>18</td><td>17</td><td>0.9995</td><td>18</td><td><b>18</b></td><td><b>1.0000</b></td><td>18</td><td>15</td><td>0.9704</td><td>18</td></tr><tr><td colspan="11" style="text-align:center; font-weight:bold;">cutoff 1h</td></tr><tr><td>PB16</td><td>1600</td><td><b>1170</b></td><td><b>0.9076</b></td><td><b>1360</b></td><td>1113</td><td>0.8842</td><td>1329</td><td>1090</td><td>0.8911</td><td>1336</td></tr><tr><td>MIPLIB</td><td>291</td><td><b>188</b></td><td><b>0.8754</b></td><td>243</td><td>149</td><td>0.8520</td><td>244</td><td>121</td><td>0.8620</td><td><b>245</b></td></tr><tr><td>CRAFT</td><td>955</td><td><b>953</b></td><td><b>&gt;0.9999</b></td><td><b>955</b></td><td><b>953</b></td><td><b>&gt;0.9999</b></td><td><b>955</b></td><td>864</td><td>0.9683</td><td>926</td></tr><tr><td>MWCB</td><td>24</td><td><b>24</b></td><td><b>1.0000</b></td><td>24</td><td>0</td><td>0.9459</td><td>24</td><td>0</td><td>0.9655</td><td>24</td></tr><tr><td>SAP</td><td>21</td><td>11</td><td>0.9986</td><td>21</td><td><b>13</b></td><td><b>0.9988</b></td><td>21</td><td>0</td><td>0.9824</td><td>21</td></tr><tr><td>WSNO</td><td>18</td><td>18</td><td>1.0000</td><td>18</td><td>18</td><td>1.0000</td><td>18</td><td>18</td><td>1.0000</td><td>18</td></tr></tbody></table>

## 5.4 Analysis on the Underlying Ideas

In order to demonstrate the effectiveness of our two main ideas in our NuPBO solver, we conduct comparative experiments on 3 solvers. We develop two alternative versions of NuPBO, by replacing its scoring function and weighting scheme with the ones proposed in LS-PBO, dubbed NuPBO-alt(s) and NuPBO-alt(w), respectively. The weighting scheme proposed in LS-PBO introduces a parameter $\zeta$ for NuPBO-alt(w), which is set to 100 as recommended by LS-PBO's authors [26].

The comparative results of the cutoff time of 300 seconds and 1 hour are shown in Table 4. From Table 4, NuPBO outperforms its alternative versions on the majority of instances. We first discuss the effectiveness of the new scoring function.

**Regarding the New Scoring Function.** With the cutoff time of 300s, in terms of the metrics of #win. and avg<sub>score</sub>, NuPBO exhibits the best performance on 5 benchmarks. In terms of the metric of #feas., NuPBO performs better than NuPBO-alt(s) on the PB16 benchmark and CRAFT benchmark. On the remaining 4 benchmarks, the number of feasible solutions obtained by NuPBO is equal to that obtained by NuPBO-alt(s).

With the cutoff time of 1h, NuPBO exhibits significantly better performance than NuPBO-alt(s) on 3 benchmarks including PB16, MIPLIB, and MWCB. On the remaining 3 benchmarks, namely CRAFT, SAP and WSNO, the performance of NuPBO is comparable to that of NuPBO-alt(s).

To examine the intuition in Section 3.1, we conduct an experiment to analyze the relationship between the instance feature and the performance difference between NuPBO and NuPBO-alt(s). Due to the difficulty of counting the coefficients of all variables in different constraints within an instance, we use the Gini coefficient [18] of the degree of hard constraints as the instance feature, denoted by $Gini_d$. For a PBO instance $I_1$, if the degree values of all hard constraints in $I_1$ are arranged in ascending order, $Gini_d$ can be calculated as follows: $Gini_d = \frac{2}{n^2d} \sum_{i=1}^{n} i(d_i - \bar{d})$, where n is the number of hard constraints, i is the
<table><caption>Table 5 Experimental results of NuPBO, LS-PBO, NuPBO-alt(s), and NuPBO-alt(w) with seeds ranging from 1 to 10 on all the benchmarks (left: cutoff 300s, right: cutoff 1h).</caption><thead><tr><th rowspan="2">Benchmark</th><th rowspan="2">#inst.</th><th colspan="3">cutoff 300s</th><th colspan="3">cutoff 1h</th></tr><tr><th>avg<sub>avgsol</sub></th><th>avg<sub>stdev</sub></th><th>avg<sub>stdev</sub>/avg<sub>avgsol</sub></th><th>avg<sub>avgsol</sub></th><th>avg<sub>stdev</sub></th><th>avg<sub>stdev</sub>/avg<sub>avgsol</sub></th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">NuPBO</td></tr><tr><td>PB16</td><td>1600</td><td>34567.32</td><td>259.79</td><td>0.75%</td><td>34721.64</td><td>123.30</td><td>0.36%</td></tr><tr><td>CRAFT</td><td>955</td><td>3035410.02</td><td>0.10</td><td>&lt;0.01%</td><td>3000465.70</td><td>0.09</td><td>&lt;0.01%</td></tr><tr><td>MIPLIB</td><td>291</td><td>59271045.04</td><td>240944.55</td><td>0.41%</td><td>58116554.74</td><td>299542.17</td><td>0.52%</td></tr><tr><td>MWCB</td><td>24</td><td>197890.62</td><td>1504.80</td><td>0.76%</td><td>193513.34</td><td>888.91</td><td>0.46%</td></tr><tr><td>SAP</td><td>21</td><td>1039.04</td><td>3.74</td><td>0.36%</td><td>1033.64</td><td>3.00</td><td>0.29%</td></tr><tr><td>WSNO</td><td>18</td><td>1301.21</td><td>225.55</td><td>17.33%</td><td>1158.61</td><td>0.00</td><td>0.00%</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">LS-PBO</td></tr><tr><td>PB16</td><td>1600</td><td>30844.62</td><td>143.86</td><td>0.47%</td><td>34305.81</td><td>213.45</td><td>0.62%</td></tr><tr><td>CRAFT</td><td>955</td><td>3074484.99</td><td>1.69</td><td>&lt;0.01%</td><td>3071190.03</td><td>1.11</td><td>&lt;0.01%</td></tr><tr><td>MIPLIB</td><td>291</td><td>53327323.12</td><td>1262751.16</td><td>2.37%</td><td>50874231.31</td><td>924287.05</td><td>1.82%</td></tr><tr><td>MWCB</td><td>24</td><td>209821.48</td><td>1582.87</td><td>0.75%</td><td>201482.90</td><td>1525.52</td><td>0.76%</td></tr><tr><td>SAP</td><td>21</td><td>1066.74</td><td>4.61</td><td>0.43%</td><td>1059.17</td><td>3.43</td><td>0.32%</td></tr><tr><td>WSNO</td><td>18</td><td>1448.88</td><td>299.06</td><td>20.64%</td><td>1174.76</td><td>44.64</td><td>3.80%</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">NuPBO-alt(s)</td></tr><tr><td>PB16</td><td>1600</td><td>36842.12</td><td>112.50</td><td>0.31%</td><td>35779.84</td><td>129.39</td><td>0.36%</td></tr><tr><td>CRAFT</td><td>955</td><td>3038627.52</td><td>0.10</td><td>&lt;0.01%</td><td>3000465.70</td><td>0.10</td><td>&lt;0.01%</td></tr><tr><td>MIPLIB</td><td>291</td><td>65612369.74</td><td>911077.41</td><td>1.39%</td><td>63282671.88</td><td>90491.84</td><td>0.14%</td></tr><tr><td>MWCB</td><td>24</td><td>210377.43</td><td>1688.62</td><td>0.80%</td><td>205368.10</td><td>1267.23</td><td>0.62%</td></tr><tr><td>SAP</td><td>21</td><td>1039.50</td><td>3.49</td><td>0.34%</td><td>1034.15</td><td>2.68</td><td>0.26%</td></tr><tr><td>WSNO</td><td>18</td><td>1295.50</td><td>192.98</td><td>14.90%</td><td>1158.65</td><td>0.12</td><td>0.01%</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">NuPBO-alt(w)</td></tr><tr><td>PB16</td><td>1600</td><td>34535.79</td><td>312.20</td><td>0.90%</td><td>34594.61</td><td>214.89</td><td>0.62%</td></tr><tr><td>CRAFT</td><td>955</td><td>3094400.67</td><td>1.45</td><td>&lt;0.01%</td><td>3094398.61</td><td>1.15</td><td>&lt;0.01%</td></tr><tr><td>MIPLIB</td><td>291</td><td>59063484.68</td><td>259177.44</td><td>0.44%</td><td>57971941.39</td><td>350572.11</td><td>0.60%</td></tr><tr><td>MWCB</td><td>24</td><td>205188.29</td><td>1467.44</td><td>0.72%</td><td>201127.28</td><td>929.57</td><td>0.46%</td></tr><tr><td>SAP</td><td>21</td><td>1061.47</td><td>4.08</td><td>0.38%</td><td>1054.01</td><td>3.41</td><td>0.32%</td></tr><tr><td>WSNO</td><td>18</td><td>1293.97</td><td>143.69</td><td>11.10%</td><td>1159.25</td><td>2.02</td><td>0.17%</td></tr></tbody></table>

## 5.5 Stability of Local Search Solvers

In order to examine the stability of all four SLS solvers adopted in our experiments, each of the four SLS solvers runs 10 times with seeds ranging from 1 to 10 on all instances from all 6 benchmarks.

For a given solver *S* and an instance *I*: *sol*<sub>*SIJ*</sub> denotes the cost of the best solution found by solver *S* with seed *J* on instance *I*, *avgsol* denotes the average cost of best solutions obtained by solver *S* over all 10 runs on instance *I*, while *stdev* denotes the standard deviation of the cost of the best solutions obtained by solver *S* over all 10 runs on instance *I*. On a benchmark *B* consisting of multiple instances: *avg<sub>avgsol</sub>* represents the average value of *avgsol* obtained by solver *S* over all instances where solutions are obtained, while *avg<sub>stdev</sub>* stands for the average value of *stdev* obtained by solver *S* over all instances where solutions are found. The calculation of *avg<sub>avgsol</sub>* is based on instances where solutions are found,
12:16 Towards More Efficient Local Search for Pseudo-Boolean Optimization

different solvers may find solutions on different subsets of instances for a given benchmark and cutoff time. In addition, for a given benchmark, it is possible that a solver finds solutions on more instances within a cutoff time of 1h than adopting a cutoff time of 300s. Moreover, according to the above definition of $avg_{avgsol}$, we note that the value of $avg_{avgsol}$ cannot be used to compare the performance of different solvers.

The experimental results presented in Table 5 demonstrate that, with the cutoff time of 300s, all four SLS solvers exhibit stable performance on 5 out of 6 benchmarks, while on the WSNO benchmark, the performance is less stable compared to the other benchmarks. In addition, the values of $avg_{stdev}/avg_{avgsol}$ for NuPBO are less than 1% on all 5 benchmarks, which clearly indicates that NuPBO can achieve stable performance. With the cutoff time of 1h, all four SLS solvers perform stably on the 6 benchmarks.

6 Conclusions and Future Work

This paper is devoted to improving the performance of SLS solvers for solving the PBO problem via a new scoring function and a new weighting scheme. First, we introduced our new scoring function. Furthermore, we proposed a new weighting scheme that effectively determines when to increase the weight of the objective function. Based on these two main ideas, we developed a new SLS solver named NuPBO. Extensive experimental results demonstrate that NuPBO significantly outperforms LS-PBO on all testing benchmarks. NuPBO outperforms all its competitors on 3 real-world application benchmarks and shows competitive performance compared to state-of-the-art competitors on solving PB16, MIPLIB, and CRAFT benchmarks. In addition, NuPBO enhances the complementarity between SLS solvers and complete solvers on all testing benchmarks.

For future work, we would like to develop more efficient heuristic strategies and explore the effect of instance features on the performances of different categories of PBO solvers.

References

1 Carlos Ansótegui, Maria Luisa Bonet, and Jordi Levy. SAT-based MaxSAT algorithms. *Artificial Intelligence*, 196:77–105, 2013.

2 Carlos Ansótegui and Joel Gabàs. WPM3: An (in)complete algorithm for weighted partial MaxSAT. *Artificial Intelligence*, 250:37–57, 2017.

3 VL Beresnev, EN Goncharov, and AA Mel'nikov. Local search with a generalized neighborhood in the optimization problem for pseudo-boolean functions. *Journal of Applied and Industrial Mathematics*, 6:22–30, 2012.

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

5 Shaowei Cai and Zhendong Lei. Old techniques in new ways: Clause weighting, unit propagation and hybridization for maximum satisfiability. *Artificial Intelligence*, 287:103354, 2020.

6 Shaowei Cai, Chuan Luo, Jinkun Lin, and Kaile Su. New local search methods for partial MaxSAT. *Artificial Intelligence*, 240:1–18, 2016.

7 Shaowei Cai, Chuan Luo, and Kaile Su. Scoring functions based on second level score for k-sat with long clauses. *J. Artif. Intell. Res.*, 51:413–441, 2014. doi:10.1613/jair.4480.

8 Shaowei Cai, Chuan Luo, John Thornton, and Kaile Su. Tailoring local search for partial maxsat. In *Proceedings of AAAI 2014*, pages 2623–2629, 2014.

9 Byungki Cha and Kazuo Iwama. Performance test of local search algorithms using new types of random CNF formulas. In *Proceedings of IJCAI 1995*, pages 304–311, 1995.

10 Byungki Cha, Kazuo Iwama, Yahiko Kambayashi, and Shuichi Miyazaki. Local search algorithms for partial MAXSAT. In *Proceedings of AAAI 1997*, pages 263–268, 1997.
11 Jessica Davies and Fahemi Bacchus. Solving MAXSAT by solving a sequence of simpler SAT instances. In *Proceedings of CP 2011*, pages 225–239, 2011.

12 Jo Devriendt, Stephan Gocht, Emir Demirovic, Jakob Nordström, and Peter J. Stuckey. Cutting to the core of pseudo-boolean optimization: Combining core-guided search with cutting planes reasoning. In *Proceedings of AAAI 2021*, pages 3750–3758, 2021.

13 Niklas Eén and Niklas Sörensson. Translating pseudo-boolean constraints into SAT. *Journal on Satisfiability, Boolean Modeling and Computation*, 2(1-4):1–26, 2006.

14 Jan Elffers and Jakob Nordström. Divide and conquer: Towards faster pseudo-boolean solving. In Jérôme Lang, editor, *Proceedings of IJCAI 2018*, pages 1291–1299, 2018.

15 Jan Elffers and Jakob Nordström. A cardinal improvement to pseudo-Boolean solving. In *Proceedings of AAAI 2020*, pages 1495–1503, 2020.

16 Gerald Gamrath, Daniel Anderson, Ksenia Bestuzheva, Wei-Kun Chen, Leon Eifler, Maxime Gasse, Patrick Gemander, Ambros Gleixner, Leona Gottwald, Katrin Halbig, Gregor Hendel, Christopher Hojny, Thorsten Koch, Pierre Le Bodic, Stephen J. Maher, Frederic Matter, Matthias Miltenberger, Erik Mühmer, Benjamin Müller, Marc E. Pfetsch, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Christine Tawfik, Stefan Vigerske, Fabian Wegscheider, Dieter Weninger, and Jakob Witzig. The SCIP Optimization Suite 7.0. Technical report, Optimization Online, March 2020. URL: http://www.optimization-online.org/DB_HTML/2020/03/7705.html.

17 Martin Gebser, Benjamin Kaufmann, and Torsten Schaub. Conflict-driven answer set solving: From theory to practice. *Artificial Intelligence*, 187:52–89, 2012.

18 Corrado Gini. Concentration and dependency ratios. *Rivista di politica economica*, pages 769–792, 1997.

19 Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2021. URL: https://www.gurobi.com.

20 Armin Haken. The intractability of resolution. *Theoretical Computer Science*, 39:297–308, 1985.

21 Pascal Van Hentenryck and Laurent Michel. *Constraint-based local search*. The MIT press, 2009.

22 Holger H. Hoos, Laetitia Jourdan, Marie-Eléonore Kessaci, Thomas Stützle, and Nadarajen Veerapen. Special issue on "stochastic local search: Recent developments and trends". *International Transactions in Operational Research*, 27(1):697–698, 2020.

23 Miyuki Koshimura, Tong Zhang, Hiroshi Fujita, and Ryuzo Hasegawa. QMaxSAT: A partial Max-SAT solver. *Journal on Satisfiability, Boolean Modeling and Computation*, 8(1-2):95–100, 2012.

24 Daniel Le Berre and Anne Parrain. The sat4j library, release 2.2. *Journal on Satisfiability, Boolean Modeling and Computation*, 7(2-3):59–64, 2010.

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

26 Zhendong Lei, Shaowei Cai, Chuan Luo, and Holger H. Hoos. Efficient local search for pseudo boolean optimization. In *Proceedings of SAT 2021*, pages 332–348, 2021.

27 Chuan Luo, Shaowei Cai, Kaile Su, and Wenxuan Huang. CCEHC: An efficient local search algorithm for weighted partial maximum satisfiability. *Artificial Intelligence*, 243:26–44, 2017.

28 Chuan Luo, Shaowei Cai, Wei Wu, Zhong Jie, and Kaile Su. CCLS: An efficient local search algorithm for weighted maximum satisfiability. *IEEE Transactions on Computers*, 64(7):1830–1843, 2015.

29 Ruben Martins, Vasco M. Manquinho, and Inês Lynce. Open-WBO: A modular maxsat solver. In *Proceedings of SAT 2014*, pages 438–445, 2014.

30 Rafiq Muhammad and Peter J. Stuckey. A stochastic non-CNF SAT solver. In *Proceedings of PRICAI 2006*, pages 120–129, 2006.

31 Alexander Nadel. Anytime weighted MaxSAT with improved polarity selection and bit-vector optimization. In *Proceedings of FMCAD 2019*, pages 193–202, 2019.
32 Nina Narodytska and Fahiem Bacchus. Maximum satisfiability using core-guided MaxSAT resolution. In *Proceedings of AAAI 2014*, pages 2717–2723, 2014.

33 Olivier Roussel and Vasco Manquinho. Pseudo-boolean and cardinality constraints. In *Handbook of satisfiability*, pages 1087–1129. IOS Press, 2021.

34 Masahiko Sakai and Hidetomo Nabeshima. Construction of an ROBDD for a PB-constraint in band form and related techniques for pb-solvers. *IEICE Transactions on Information & Systems*, 98-D(6):1121–1127, 2015.

35 Bart Selman and Henry Kautz. Domain-independent extensions to GSAT: solving large structured satisfiability problems. In *Proceedings of IJCAI 1993*, pages 290–295, 1993.

36 Pavel Smirnov, Jeremias Berg, and Matti Järvisalo. Pseudo-boolean optimization by implicit hitting sets. In *Proceedings of CP 2021*, pages 51:1–51:20, 2021.

37 Pavel Smirnov, Jeremias Berg, and Matti Järvisalo. Improvements to the implicit hitting set approach to pseudo-boolean optimization. In *Proceedings of SAT 2022*, pages 13:1–13:18, 2022.

38 John Thornton and Abdul Sattar. Dynamic constraint weighting for over-constrained problems. In *Proceedings of PRICAI 1998*, pages 377–388, 1998.

39 Renato Tinós, Michal W Przewozniczek, and Darrell Whitley. Iterated local search with perturbation based on variables interaction for pseudo-boolean optimization. In *Proceedings of GECCO 2022*, pages 296–304, 2022.

40 Chris Voudouris and Edward Tsang. Partial constraint satisfaction problems and guided local search. *Proc., Practical Application of Constraint Technology (PACT'96)*, London, pages 337–356, 1996.

41 Christos Voudouris, Edward PK Tsang, and Abdullah Alsheddy. Guided local search. In *Handbook of metaheuristics*, pages 321–361. Springer, 2010.

42 Joachim P. Walser. Solving linear pseudo-boolean constraint problems with local search. In *Proceedings of AAAI 1997*, pages 269–274, 1997.

43 Robert Wille, Hongyan Zhang, and Rolf Drechsler. ATPG for reversible circuits using simulation, boolean satisfiability, and pseudo boolean optimization. In *Proceedings of ISVLSI 2011*, pages 120–125, 2011.

44 Yuhang Zhang, Richard Hartley, John Mashford, and Stewart Burn. Superpixels via pseudo-boolean optimization. In *Proceedings of ICCV 2011*, pages 1387–1394, 2011.

45 Jiongzhi Zheng, Kun He, Jianrong Zhou, Yan Jin, Chu-Min Li, and Felip Manya. BandMaxSAT: A local search maxsat solver with multi-armed bandit. In *Proceedings of IJCAI 2022*, pages 1901–1907, 2022.
# Boosting Decision Diagram-Based Branch-And-Bound by Pre-Solving with Aggregate Dynamic Programming

Vianney Coppé

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

UCLouvain, Louvain-la-Neuve, Belgium

Xavier Gillard

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

UCLouvain, Louvain-la-Neuve, Belgium

Pierre Schaus

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

UCLouvain, Louvain-la-Neuve, Belgium

## Abstract

Discrete optimization problems expressible as dynamic programs can be solved by branch-and-bound with decision diagrams. This approach dynamically compiles bounded-width decision diagrams to derive both lower and upper bounds on unexplored parts of the search space, until they are all enumerated or discarded. Assuming a minimization problem, relaxed decision diagrams provide lower bounds through state merging while restricted decision diagrams obtain upper bounds by excluding states to limit their size. As the selection of states to merge or delete is done locally, it is very myopic to the global problem structure. In this paper, we propose a novel way to proceed that is based on pre-solving a so-called aggregate version of the problem with a limited number of states. The compiled decision diagram of this aggregate problem is tractable and can fit in memory. It can then be exploited by the original branch-and-bound to generate additional pruning and guide the compilation of restricted decision diagrams toward good solutions. The results of the numerical study we conducted on three combinatorial optimization problems show a clear improvement in the performance of DD-based solvers when blended with the proposed techniques. These results also suggest an approach where the aggregate dynamic programming model could be used in replacement of the relaxed decision diagrams altogether.

**2012 ACM Subject Classification** Mathematics of computing → Combinatorial optimization

**Keywords and phrases** Discrete Optimization, Decision Diagrams, Aggregate Dynamic Programming

Digital Object Identifier 10.4230/LIPIcs.CP.2023.13

## 1 Introduction

On top of their use for Boolean encodings [27], formal verification [25], model checking [15], computer-aided design [29] and much more, *decision diagrams* (DDs) have recently emerged as a tool for discrete optimization. They provide a compact way to encode a set of solutions to a problem. Still, for large problems, DDs representing the whole solution space – called *exact* DDs – can quickly become intractable to compute. Two variants of DDs can be used instead: *restricted* [10] and *relaxed* [1, 8] DDs that respectively encode a subset and superset of the set of solutions. When compiled based on a *dynamic programming* (DP) model, these approximate DDs allow to compute bounds on the objective function for any subproblem while controlling the size of the DD compiled. Restricted DDs aim to find good admissible solutions by iteratively extending a bounded set of promising candidates while dropping others, in a beam search fashion. On the other hand, relaxed DDs rely on a problem-dependent state merging scheme to maintain an acceptable DD size while preserving all solutions of the problem. In [9], Bergman et al. presented a *branch-and-bound* algorithm solely based on these two ingredients, thus introducing a new general-purpose discrete optimization framework and solver.

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

© Vianney Coppé, Xavier Gillard, and Pierre Schaus;
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. 13; pp. 13:1–13:17

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

LIPICS

Leibniz International Proceedings in Informatics

Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
13:2 Boosting DD-Based Branch-And-Bound with Aggregate Dynamic Programming

In addition to exploiting the compactness of DP models, the main novelty of this
approach is its unique way of deriving lower and upper bounds. In the last few years, some
algorithmic improvements have been suggested to further strengthen these bounds. Assuming
a minimization problem, Gillard et al. [19] showed how user-defined lower bound formulas
can be integrated to prune DDs during their compilation and thus concentrate the search
on promising parts of the search space. They also proposed a way to compute tighter lower
bounds for all nodes contained in a relaxed DD through *local bounds*. Rudich et al. [30]
introduced a *peeling* operator that splits a relaxed DD in two: one part containing all paths
traversing a selected exact node and the other containing all remaining paths. It allows both
to warm-start the compilation of subsequent relaxed DDs and to strengthen the bounds of
the nodes inside the relaxed DD on which the peeling has been performed. More recently, [16]
generalized the ideas of [19] and introduced the use of a cache storing new thresholds that
further enhance the pruning power of the solver. Other factors impacting the quality of the
bounds provided by relaxed DDs have been studied, including variable orderings [7, 11, 26]
and alternative compilation schemes [24]. Yet, all these approaches rely on a problem-specific
state merging operator at the heart of the relaxation, which does not yield tight relaxations
for all problems, as our computational experiments show.

After covering the necessary background about DD-based optimization, this paper presents an alternate relaxation scheme for deriving good bounds by incorporating ideas from *aggregate dynamic programming* [2, 3] to the DD-based discrete optimization framework. The underlying idea of the approach is to deduce information about an original problem instance by creating and solving an aggregate – relaxed – version of it. This is achieved by *aggregating* the states of the DP model as to obtain a much smaller DP state space. If this aggregation is adequately specified, one can compute a lower bound for any original subproblem by finding the optimal solution of its aggregate version. Furthermore, this optimal aggregate solution can be *disaggregated* and transposed in the original problem to find good heuristic solutions. In practice, the aggregation-based lower bounds are used as additional pruning within the compilation of relaxed and restricted DDs. Moreover, aggregate solutions are translated into node selection heuristics to steer the compilation of restricted DDs toward resembling solutions to the original problem, which are thus expected to be good.

Throughout the paper, the framework is illustrated on three different combinatorial
problems: the Talent Scheduling Problem, the Pigment Sequencing Problem and the Aircraft
Landing Problem. They are then used for the experimental evaluation of the framework,
the results of which show that the aggregation-based bound brings additional pruning and
enables solving more instances. Furthermore, the aggregation-based node selection heuristic
improves the quality of the solutions found early in the search and thus contributes to
speeding up the overall resolution. Finally, we show that a DD-based solver using only the
aggregation-based bound as relaxation performs almost equally well, which is a promising
direction for problems for which defining a merging operator is difficult or inefficient.

Although this paper is – to the best of our knowledge – the first to combine aggregate
dynamic programming with the DD-based branch-and-bound paradigm proposed by Bergman
et al, there has already been some hybridization work to combine discrete optimization with
DDs and other methods. For instance, in [12], Cappart et al. propose to use reinforcement
learning to guess the variable ordering that should be used to derive the best possible bounds
from the compiled approximate DDs. Other attempts combined DDs with Lagrangian
relaxation [13, 23] or MIP [5, 22, 31, 32]. On a slightly different note, a method has
been proposed where restricted DDs are used to generate good neighborhoods in a *large
neighborhood search* framework [20].
# 2 Preliminaries

## 2.1 Discrete Optimization

A discrete optimization problem $\mathcal{P}$ involves finding the best possible solution $x^*$ from a finite set of feasible solutions $Sol(\mathcal{P}) = D \cap C$. This set is determined by the domain $D = D_0 \times \dots \times D_{n-1}$ from which the variables $x = (x_0, \dots, x_{n-1})$ each take on a value, i.e. $x_j \in D_j$, and by a set of constraints $C$ imposed on the value assignments. The quality of the solutions is evaluated according to an objective function $f(x)$ that must be minimized. Formally, the problem is defined as $\min \{f(x) | x \in D \cap C\}$ and any optimal solution $x^*$ must satisfy $x^* \in Sol(\mathcal{P})$ and $\forall x \in Sol(\mathcal{P}) : f(x^*) \le f(x)$. We describe below three optimization problems that will be utilized in the paper as illustrations for the aggregation-based framework.

**Talent Scheduling Problem.** The Talent Scheduling Problem (TalentSched) is a film shoot scheduling problem that considers a set $N = \{0, \dots, n-1\}$ of scenes and a set $A = \{0, \dots, m-1\}$ of actors. Each scene $i \in N$ involves a required set $R_i \subseteq A$ of actors for a duration $D_i \in \mathbb{N}$. Moreover, each actor $k \in M$ has a pay rate $C_k$ and is paid without interruption from their first to their last scheduled scene. The objective of TalentSched is to find a permutation of the scenes that minimizes the total cost of the film shoot.

**Pigment Sequencing Problem.** The Pigment Sequencing Problem (PSP) is a single-machine production planning problem that aims to minimize the stocking and changeover costs while satisfying a set of orders. There are different item types $I = \{0, \dots, n-1\}$ with a given stocking cost $S_i$ to pay for each time period between the production and the deadline of an order. For each pair $i, j \in I$ of item types, a changeover cost $C_{ij}$ is incurred whenever the machine switches the production from item type $i$ to $j$. Finally, the demand matrix $Q$ contains all the orders: $Q_p^i \in \{0, 1\}$ indicates whether there is an order for item type $i \in I$ at time period $p$ with $0 \le p < H$ and $H$ the time horizon.

**Aircraft Landing Problem.** The Aircraft Landing Problem (ALP) requires to schedule the landing of a set of aircrafts $N = \{0, \dots, n-1\}$ on a set of runways $R = \{0, \dots, r-1\}$. The aircrafts have a target $T_i$ and latest $L_i$ landing time. Moreover, the set of aircrafts is partitioned in disjoint sets $A_0, \dots, A_{c-1}$ corresponding to different aircraft classes in $C = \{0, \dots, c-1\}$. For each pair of aircraft classes $a, b \in C$, a minimum separation time $S_{a,b}$ between the landings is given. The goal is to find the schedule that minimizes the total waiting time of the aircrafts – the delay between their target time and scheduled landing time – while respecting their latest landing time.

## 2.2 Dynamic Programming

Dynamic programming (DP) is a divide-and-conquer strategy introduced by Bellman [4] for solving discrete optimization problems with an inherent recursive structure. It works by recursively decomposing the problem in smaller and overlapping subproblems. The cornerstone of the approach is the caching of intermediate results that allows each distinct subproblem to be solved only once. A DP model of a discrete optimization problem $\mathcal{P}$ can be defined as a labeled transition system consisting of:

- the control variables $x_j \in D_j$ with $j \in \{0, \dots, n-1\}$.

- a set of state-spaces $S = \{S_0, \dots, S_n\}$ among which one distinguishes the initial state $r$, the terminal state $t$ and the infeasible state $\hat{0}$.
13:4 Boosting DD-Based Branch-And-Bound with Aggregate Dynamic Programming

a set t of transition functions s.t. $t_j : S_j \times D_j \rightarrow S_{j+1}$ for $j = 0, \dots, n-1$ taking the system from one state $s^j$ to the next state $s^{j+1}$ based on the value d assigned to variable $x_j$, or to $\perp$ if assigning $x_j = d$ is infeasible. These functions should never allow one to recover from infeasibility, i.e. $t_j(\hat{0}, d) = \hat{0}$ for any $d \in D_j$.

a set h of transition value functions s.t. $h_j : S_j \times D_j \rightarrow \mathbb{R}$ representing the immediate reward of assigning some value $d \in D_j$ to the variable $x_j$ for $j = 0, \dots, n-1$.

a root value $v_r$.

On that basis, the objective function $f(x)$ of $\mathcal{P}$ is formulated as follows:

$$ \text{minimize} \quad f(x) = v_r + \sum_{j=0}^{n-1} h_j(s^j, x_j) $$

subject to $s^{j+1} = t_j(s^j, x_j)$, for all $j = 0, \dots, n-1$, with $x_j \in D_j$

$$ s^j \in S_j, j = 0, \dots, n \text{ and } x \in C. \tag{1} $$

**TalentSched.** A DP model for TalentSched was introduced in [17] that we slightly adapt here to make it suitable for the relaxation discussed in Section 2.3.1. States of this model are pairs $(M, P)$ where $M$ and $P$ are disjoint sets of scenes that respectively must or might still be scheduled. The only case where $P$ is non-empty happens when a state is relaxed.

Control variables: $x_j \in N$ with $0 \le j < n$ decides which scene is shot in j-th position.

State spaces: $S = \{(M, P) | M, P \subseteq N, M \cap P = \emptyset\}$. The root state is $r = (N, \emptyset)$ and the terminal states are of the form $(\emptyset, P)$.

Transition functions:

$$ t_j(s^j, x_j) = \begin{cases} (s^j.M \setminus \{x_j\}, s^j.P \setminus \{x_j\}) & \text{if } x_j \in s^j.M, \\ (s^j.M \setminus \{x_j\}, s^j.P \setminus \{x_j\}), & \text{if } x_j \in s^j.P \text{ and } |s^j.M| < n-j, \\ \hat{0}, & \text{otherwise.} \end{cases} $$

A scene from $P$ can only be selected if there are more spots left than scenes in $M$.

Transition value functions: let $a(Q) = \cup_{i \in Q} R_i$ be the required set of actors for a set of scenes $Q$. Given a state $s = (M, P)$, the set of actors that are guaranteed to be on-location is computed as $o(s) = a(s.M) \cap a(N \setminus (s.M \cup s.P))$ because they are required both for a scene that must still be scheduled and for another that is guaranteed to be scheduled. In the transition value functions, we add all the actors from $R_{x_j}$ to this set and sum the individual costs: $h_j(s^j, x_j) = D_{x_j} \sum_{k \in o(s^j) \cup R_{x_j}} C_k$.

Root value: $v_r = 0$.

**PSP.** The PSP was already tackled with a DD-based approach in [16, 20]. We hereby recall the DP model from [16] that allows the machine to be idle at some time periods. In this model, the decisions are made backwards – this allows to define transition functions that only lead to feasible production schedules. If variable $x_j$ decides the type of item to produce at period $j$, the reverse variable ordering $x_{H-1}, \dots, x_0$ is thus used. To simplify the transition functions, let us denote by $P_r^i$ the time period at which the r-th item of type i must be delivered, i.e. $P_r^i = \min\{0 \le q < H \mid \sum_{p=0}^q Q_p^i \ge r\}$ for all $i \in N, 0 \le r \le \sum_{0 \le p < H} Q_p^i$. Moreover, we define a dummy item type $\perp$ used for idle periods and $N' = N \cup \{\perp\}$.

States are pairs $(i, R)$ with $i$ the item type that the machine is currently set to produce and $R$ a vector that gives the remaining number $R_i$ of demands to satisfy for each type $i$.

Control variables: $x_j \in N'$ with $0 \le j < H$ decides the item type to produce at period $j$.
■ State space: $S = \{s \mid s.i \in N', \forall i \in N, 0 \le s.R_i \le \sum_{0 \le p < H} Q_p^i\}$. The root state is given by $r = \langle \perp, (\sum_{0 \le p < H} Q_p^0, \dots, \sum_{0 \le p < H} Q_p^{n-1}) \rangle$ and the terminal states are of the form $\langle i, (0, \dots, 0) \rangle$ with $i \in N'$.

■ Transition functions:

$$t_j(s^j, x_j) = \begin{cases} \langle t_j^i(s^j, x_j), t_j^R(s^j, x_j) \rangle, & \text{if } x_j \neq \perp \text{ and } s^j.R_{x_j} > 0 \text{ and } j \le P_{s^j.R_{x_j}}^{x_j}, \\ \langle t_j^i(s^j, x_j), t_j^R(s^j, x_j) \rangle, & \text{if } x_j = \perp \text{ and } \sum_{i \in N} s^j.R_i < j + 1, \\ \hat{0}, & \text{otherwise.} \end{cases}$$

where

$$\begin{aligned} t_j^i(s^j, x_j) &= \begin{cases} x_j, & \text{if } x_j \neq \perp \\ s^j.i, & \text{otherwise.} \end{cases} \\ t_j^R(s^j, x_j) &= \begin{cases} (s^j.R_0, \dots, s^j.R_{x_j} - 1, \dots, s^j.R_{n-1}), & \text{if } x_j \neq \perp \\ s^j.R, & \text{otherwise.} \end{cases} \end{aligned}$$

In the transition function, the first condition ensures that there remains at least one item to produce for the chosen type and that the current time period $j$ is earlier than its deadline. The second condition ensures that idle periods cannot be scheduled when the remaining quantity to produce is equal to the number of periods left.

■ Transition value functions: the changeover and stocking costs are computed as:

$$h_j(s^j, x_j) = \begin{cases} C_{x_j s^j i}, & \text{if } x_j \neq \perp \text{ and } s^j.i \neq \perp \\ 0, & \text{otherwise.} \end{cases} + \begin{cases} S_{x_j} \cdot (j - P_{s^j.R_{x_j}}^{x_j}), & \text{if } x_j \neq \perp \\ 0, & \text{otherwise.} \end{cases}$$

■ Root value: $v_r = 0$.

**ALP.** We reproduce here the DP model presented in [28] where states are pairs $(Q, ROP)$, with $Q$ a vector that gives the remaining number of aircrafts of each class to schedule and $ROP$ a runway occupation profile: a vector containing pairs $(l, c)$ that respectively give the time and aircraft class of the latest landing scheduled on each runway. Similarly to the PSP modeling, we denote by $\perp$ either a dummy aircraft class or a dummy runway.

■ Control variables: we use pairs of variables $(x_j, y_j) \in (C \times R) \cup \{(\perp, \perp)\}$ with $0 \le j < n$ that represent the decision to place an aircraft of class $x_j$ on runway $y_j$, or to schedule nothing at all in case of $(\perp, \perp)$.

■ State spaces:

$S = \{(Q, ROP) \mid \forall i \in C: Q_i \ge 0, \forall k \in R: ROP_k.l \ge 0, ROP_k.c \in C \cup \{\perp\}\}$. The root state is $r = ((|A_0|, \dots, |A_{c-1}|), ((0, \perp), \dots, (0, \perp)))$ and the terminal states are of the form $((0, \dots, 0), ROP)$.

■ Transition functions: if $A_i^k$ gives the aircraft from class $i$ that must be scheduled when there are $k$ aircrafts left from this class, we can define the function computing the earliest landing time given a state $s$, a class $x$ and a runway $y$:

$$E(s, x, y) = \begin{cases} T_{A_x^{s^j.Q_x}}, & \text{if } s.ROP_y.l = 0 \text{ and } s.ROP_y.c = \perp, \\ \max(s.ROP_y.l + \min_{i \in C} S_{i,x}, T_{A_x^{s^j.Q_x}}, & \text{if } s.ROP_y.l > 0 \text{ and } s.ROP_y.c = \perp, \\ \max(s.ROP_y.l + S_{s.ROP_y.c,x}, T_{A_x^{s^j.Q_x}}, & \text{otherwise.} \end{cases}$$

This allows us to define the transition functions as:

$$t_j(s^j, x_j, y_j) = \begin{cases} (t_j^Q(s^j, x_j, y_j), t_j^{ROP}(s^j, x_j, y_j)), & \text{if } x_j \neq \perp \text{ and } s^j.Q_{x_j} > 0 \\ (t_j^Q(s^j, x_j, y_j), t_j^{ROP}(s^j, x_j, y_j)), & \text{if } x_j = \perp \text{ and } \sum_{i \in C} s^j.Q_i = 0, \\ \hat{0}, & \text{otherwise.} \end{cases}$$
13:6 Boosting DD-Based Branch-And-Bound with Aggregate Dynamic Programming

where

$$
\begin{align*}
t_j^Q(s^j, x_j, y_j) &= \begin{cases} \langle s^j.Q_0, \dots, s^j.Q_{x_j} - 1, \dots, s^j.Q_{c-1} \rangle, & \text{if } x_j \neq \perp \\ s^j.Q, & \text{otherwise.} \end{cases} \\
t_j^{ROP}(s^j, x_j, y_j) &= \begin{cases} \langle s^j.ROP_0, \dots, (E(s^j, x_j, y_j), x_j), \dots, s^j.ROP_{r-1} \rangle, & \text{if } x_j \neq \perp \\ s^j.ROP, & \text{otherwise.} \end{cases}
\end{align*}
$$

The first condition of the transition function ensures that there remains at least one aircraft of the chosen class and that its earliest landing time is not greater its latest landing time. The second condition only allows us to schedule dummy aircrafts when there are no aircrafts left to schedule.

Transition value functions: the waiting time of the aircraft is computed as:

$$
h_j(s^j, x_j, y_j) = \begin{cases} E(s^j, x_j, y_j) - T_{A_{x_j}^{s^j.Q_{x_j}}}, & \text{if } x_j \neq \perp \\ 0, & \text{otherwise.} \end{cases}
$$

Root value: $v_r = 0$.

Because the runways are identical and independent, there are many symmetries in this model.
This can be mitigated by sorting the ROP of every state by increasing latest landing time,
breaking ties according to the previous aircraft class scheduled.

2.3 Decision Diagrams

When used to manipulate the DP model of a discrete optimization problem $\mathcal{P}$, DDs are
graphical encodings that represent a set of solutions of the problem. More precisely, a
DD $\mathcal{B} = (U, A, \sigma, l, v)$ is a layered directed acyclic graph composed of a set of nodes $U$
interconnected by a set of arcs $A$. Starting from a single node $u_r$ corresponding to a DP state
given by the function $\sigma(u_r)$, the process of iteratively extending a set of partial solutions is
called the *compilation* of a DD and is described by Algorithm 1. Note that the highlighted
portions concern the ingredients introduced in Section 3 and can be ignored for now. The
algorithm begins by initializing a layer $L_i$ that only contains the *root node* $u_r$, assuming its
state $\sigma(u_r)$ belongs to the *i*-th stage of the DP model. The subsequent layers of the DD are
then constructed sequentially by applying each valid transition of the DP model to every
node of the last completed layer at lines 8–16. Each layer thus corresponds to a stage of
the DP model and contains a single node for each state reached in order to preserve the
compactness of the model. The arcs $a \in A$ materialize the transitions that exist between
the states of consecutive stages. In particular, the arc $a = (u \xrightarrow{d} u')$ connecting nodes
$u \in L_j, u' \in L_{j+1}$ represents the transition between $\sigma(u)$ and $\sigma(u')$. The decision associated
with this transition is stored by the *label* $l(a) = d \in D_j$ and the transition value is given by
the arc *value* $v(a)$.

The algorithm completes when the last layer $L_n$ is generated, constituted by a single
node $t$ called the *terminal node*. The DD thus constructed contains a set of $u_r \rightsquigarrow t$
paths that can be combined with any previously discovered $r \rightsquigarrow u_r$ path, connecting
the *root* of the problem to $u_r$. Any $r \rightsquigarrow t$ path $p = (a_0, \dots, a_{n-1})$ represents a solution
given by $x(p) = (l(a_0), \dots, l(a_{n-1}))$. The objective value of such solution can also be
retrieved from the sequence of arcs by accumulating their values, and adding the root
value: $v(p) = v_r + \sum_{j=0}^{n-1} v(a_j)$. The set of solutions contained in the DD is denoted as
$\text{Sol}(\mathcal{B}) = \{x(p) \mid \exists p : r \rightsquigarrow t, p \in \mathcal{B}\}$. A DD rooted at a node $u_r$ is *exact* if it perfectly
represents the set of solutions of the corresponding subproblem $\mathcal{P}|_{u_r}$, i.e. $\text{Sol}(\mathcal{B}) = \text{Sol}(\mathcal{P}|_{u_r})$
and $v(p) = f(x(p)), \forall p \in \mathcal{B}$. The best value among the $u_1 \rightsquigarrow u_2$ paths in $\mathcal{B}$ is denoted
$v^*(u_1 \rightsquigarrow u_2 | \mathcal{B})$, and in particular $v^*(u | \mathcal{B}) = v^*(r \rightsquigarrow u | \mathcal{B})$.
**Algorithm 1** Compilation of DD $\mathcal{B}$ rooted at node $u_r$ with maximum width $W$.
1: $i \leftarrow$ index of the layer containing $u_r$
2: $L_i \leftarrow \{u_r\}$
3: $\tilde{P} \leftarrow \Delta(\tilde{p})$ with $\tilde{p}$ the optimal solution for $\pi(\sigma(u_r))$ // retrieve disaggregate solution
4: **for** $j = i$ **to** $n-1$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;**if** $|L_j| > W$ **then**
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;restrict or relax the layer to get $W$ nodes with Algorithm 2
7: &nbsp;&nbsp;&nbsp;&nbsp;$L_{j+1} \leftarrow \emptyset$
8: &nbsp;&nbsp;&nbsp;&nbsp;**for all** $u \in L_j$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$v_{rlb}(\sigma(u)) \leftarrow \max\{v_{rlb}(\sigma(u)), v_{agg}(\pi(\sigma(u)))\}$ // inject aggregation-based bound
10: &nbsp;&nbsp;&nbsp;&nbsp;**if** $v^*(u | \mathcal{B}) + v_{rlb}(\sigma(u)) \ge \bar{v}$ **then** // rough lower bound pruning w.r.t. incumbent
11: &nbsp;&nbsp;&nbsp;&nbsp;**continue**
12: &nbsp;&nbsp;&nbsp;&nbsp;**for all** $d \in D_j$ **do**
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;create node $u'$ with state $\sigma(u') = t_j(\sigma(u), d)$ or retrieve it from $L_{j+1}$
14: &nbsp;&nbsp;&nbsp;&nbsp;create arc $a = (u \xrightarrow{d} u')$ with $v(a) = h_j(\sigma(u), d)$ and $l(a) = d$
15: &nbsp;&nbsp;&nbsp;&nbsp;score(a) $\leftarrow 1$ if $l(a) \in \tilde{P}_j$, 0 otherwise
16: &nbsp;&nbsp;&nbsp;&nbsp;add $u'$ to $L_{j+1}$ and add $a$ to $A$

**► Example 1.** Let us define a TalentSched instance with 4 scenes with durations $D = \langle 3, 5, 2, 4 \rangle$ and 4 actors with pay rates $C = \langle 10, 20, 30, 40 \rangle$. The actor requirements for each scene are given by $R = \langle \{0, 3\}, \{0, 1, 3\}, \{0, 2, 3\}, \{0, 1, 2, 3\} \rangle$. Figure 1 shows the exact DD compiled for this instance with the DP model recalled in Section 2.2. Note that for each state $s = (M, P)$ corresponding to a node in the DD, we only show the set $M$ since $P$ is always empty in exact nodes. An optimal solution of the problem is $\langle 0, 2, 3, 1 \rangle$, which gives an objective value of 106.

As the reader might have guessed, the compilation of an exact DD for a combinatorial optimization problem suffers from the curse of dimensionality as much as the corresponding DP model. This is why DD-based discrete optimization rarely relies on exact DDs but rather on *restricted* and *relaxed* DDs. These two variants follow two distinct compilation schemes that allow to maintain the number of nodes of each layer – called the *width* – under a given parameter $W$. In Algorithm 1, this logic is performed at line 5 where the width of the current layer is compared with $W$. If needed, the layer is then either restricted or relaxed at line 6 by calling Algorithm 2.

### 2.3.1 Approximate Decision Diagrams

As stated by Algorithm 2, restricted DDs simply remove surplus nodes from the layer until it is reduced to $W$ nodes. A heuristic is used to evaluate the nodes and drop the least promising ones. Restricted DDs thus generate a subset of the solutions of the corresponding problem, i.e. $Sol(\bar{\mathcal{B}}) \subseteq Sol(\mathcal{P})$ and $v(p) = f(x(p)), \forall p \in \bar{\mathcal{B}}$ for a restricted DD $\bar{\mathcal{B}}$. They thus provide upper bounds on the objective value.

As opposed to restricted DDs, a relaxed DD $\mathcal{B}$ yields lower bounds by representing a superset of the solutions of the corresponding problem: $Sol(\bar{\mathcal{B}}) \supseteq Sol(\mathcal{P})$ and $v(p) \le f(x(p)), \forall p \in \bar{\mathcal{B}}$. This is achieved through a problem-specific *state merging* operator $\oplus(\sigma(\mathcal{M}))$ that defines an approximate representation that includes all states $\sigma(\mathcal{M}) = \{\sigma(u) | u \in \mathcal{M}\}$ corresponding to the merged nodes $\mathcal{M}$ and preserves all their outgoing transitions, although
13:8 Boosting DD-Based Branch-And-Bound with Aggregate Dynamic Programming

<figure><img src="image_126.png" /><figcaption>Figure 1 The exact DD for the TalentSched instance given in Example 1. Nodes are annotated with their state and the best prefix value. Arcs are labeled with the associated decision in bold and transition value. The arcs constituting one of the optimal solutions are highlighted in bold.</figcaption></figure>

**Algorithm 2** Restriction or relaxation of layer $L_j$ with maximum width $W$.

1: **while** $|L_j| > W$ **do**

2: $\mathcal{M} \leftarrow$ select nodes from $L_j$ according to their score

3: $L_j \leftarrow L_j \setminus \mathcal{M}$

4: create node $\mu$ with state $\sigma(\mu) = \oplus(\sigma(\mathcal{M}))$ and add it to $L_j$ // for relaxation only

5: **for all** $u \in \mathcal{M}$ **and** arc $a = (u' \xrightarrow{d} u)$ incident to $u$ **do**

6: replace $a$ by $a' = (u' \xrightarrow{d} \mu)$ and set $v(a') = \Gamma_{\mathcal{M}}(v(a), u)$

it may also introduce infeasible transitions. In Algorithm 2, a *meta-node* is created for this merged state at line 4 and the arcs pointing to the deleted nodes are redirected to this merged node at line 6. The operator $\Gamma_{\mathcal{M}}$ permits to adjust the value of these arcs if needed. In all three formulations given below, this operator is the identity function.

**TalentSched.** The merging operator is defined by $\oplus(\mathcal{M}) = (\oplus_M(\mathcal{M}), \oplus_P(\mathcal{M}))$ where $\oplus_M(\mathcal{M}) = \bigcap_{s \in \mathcal{M}} s.M$ and $\oplus_P(\mathcal{M}) = (\bigcup_{s \in \mathcal{M}} s.M \cup s.P) \setminus (\bigcap_{s \in \mathcal{M}} s.M)$. The definition of $\oplus_P(\mathcal{M})$ ensures that the resulting set of scenes that might be scheduled contains any scene that must or might be scheduled in any of the states, except those that still must be scheduled for all states.

**PSP.** A valid merging operator is $\oplus(\mathcal{M}) = (\perp, \langle \min_{s \in \mathcal{M}} s.R_0, \dots, \min_{s \in \mathcal{M}} s.R_{n-1} \rangle)$. The configuration of the machine is always reset to the dummy item type $\perp$ as there is little chance that merged states agree on it. For each item type, the remaining number of demands is computed by taking the minimum value among all merged states, meaning that any demand satisfied by at least one state is considered satisfied in the merged state.
<figure><img src="image_127.png" /><figcaption>Figure 2 Respectively on the left and the right, a restricted and relaxed DD for the TalentSched instance given in Example 1, compiled with *W* set to 3 and 4. Merged nodes are circled twice.</figcaption></figure>

**ALP.** The merging operator is again defined separately for each component of the states: $\oplus(\mathcal{M}) = (\oplus_Q(\mathcal{M}), \oplus_{ROP}(\mathcal{M}))$. First, the minimum remaining quantity of aircrafts for each class is stored in the merged state: $\oplus_{\mathcal{M}}(\mathcal{M}) = \langle \min_{s \in \mathcal{M}} s.Q_0, \dots, \min_{s \in \mathcal{M}} s.Q_{c-1} \rangle$. For the ROP, the minimum latest landing time on each runway is kept and the last aircraft classes scheduled are reset to $\perp$: $\oplus_{ROP}(\mathcal{M}) = \langle (\min_{s \in \mathcal{M}} s.ROP_0.l, \perp), \dots, (\min_{s \in \mathcal{M}} s.ROP_{r-1}.l, \perp) \rangle$.

**► Example 2.** Figure 2 shows approximate DDs for the TalentSched instance introduced in Example 1. Despite having a maximum width of 3, the best solution contained in the restricted DD is the optimal solution previously found. With a maximum width of 4, the relaxed DD provides a global lower bound of 86. The path corresponding to this lower bound is given by the assignment $\langle 0, 2, 3, 0 \rangle$, which is infeasible because scene 0 is scheduled twice.

### 2.3.2 Branch-and-Bound

In [9], a branch-and-bound algorithm based only on restricted and relaxed DDs was introduced. It maintains a queue of open nodes that represent the set of subproblems that remain to process. For each of them, a restricted DD is compiled in an attempt to improve the incumbent solution. Then, a relaxed DD is constructed in order to both decompose the given subproblem into even smaller ones and to compute a lower bound for each of them. These nodes are then added to the branch-and-bound queue for further exploration, unless the lower bound permits their direct elimination. Ultimately, the optimality of the best solution discovered during the search is confirmed once the queue has been emptied.
### 2.3.3 Rough Lower Bound

The *rough lower bound* (RLB) [19] is an additional optional modeling component that can be specified to speed up the resolution of any optimization problem. For any node *u*, the RLB gives a lower bound on the best value one can obtain when solving the corresponding subproblem $\sigma(u)$, i.e. $\underline{v}_{rlb}(\sigma(u)) \le v^*(u \rightsquigarrow t | \mathcal{B})$ with $\mathcal{B}$ the exact DD for the problem. It is used at line 10 of Algorithm 1 to filter nodes *a priori* by comparing this lower bound with the incumbent value $\bar{v}$. Since the RLB is computed for each node of the approximate DDs compiled throughout the branch-and-bound, it needs to be computationally cheap.

The RLB has the potential both to focus the compilation of restricted DDs on promising parts of the search space and to strengthen the bounds obtained through relaxed DDs. Furthermore, the branch-and-bound algorithm uses the RLB to make pruning decisions, if it happens to be tighter than the bound obtained with relaxed DDs.

**Example problems.** In our computational experiments, we use the lower bound given by Theorem 1 in [17] for TalentSched and the same RLB as in [20] for the PSP. We do not detail them in this article for the sake of conciseness.

## 3 Aggregate Dynamic Programming for Decision Diagrams

As stated in the introduction, optimizations techniques based on DP and DDs can prove highly effective [6, 13, 14, 18, 19]. In some cases, however, the state space of the DP models is simply too large and the bounds derived from restricted and relaxed DDs are of little to no use. This can be imputed either to the node selection heuristic or to the relaxation scheme. The *MinLP* heuristic traditionally used favors keeping nodes with the best prefix values. This locally-optimal selection policy may result in the elimination of all nodes that lead to the optimal solution, or even to any feasible solution, particularly in cases of highly constrained problems. In the latter case, the compilation of a restricted DD is a pure waste of time: no feasible solution is found at the end of the compilation, and not even a bound on the objective value can be exploited to reduce the optimality gap. The same phenomenon is detrimental to the usefulness of compiled relaxed DDs whose bounds might be of low quality when the node selection heuristic is oblivious to the global structure of the problem. Indeed, the merging operator yields a loose representation when applied to an arbitrary set of nodes for most problems. In the absence of a perfect heuristic, this situation will occur under certain conditions. It inspired our pursuit of a more globally-focused approach that could enhance the usefulness of the compiled DDs. This section presents a framework for integrating *aggregate dynamic programming* ideas with DD-based optimization that aims to address some of these shortcomings. Instead of relaxing the original problem by reasoning on merged states, it proposes to use problem instance and state aggregation operators that yield a simpler and relaxed version of the problem, which can be solved exactly. Solutions of the aggregated problem can provide bounds that capture the global problem structure, as well as guidance for the compilation of restricted DDs. This section details the role and meaning of the components of the framework one by one.

### 3.1 Preprocessing: Problem Instance Aggregation

The goal of this preprocessing step is to create an aggregate and simpler problem instance by reducing one or more dimensions of the problem. The *instance aggregation operator* $\Pi$ must be defined such that the aggregate problem instance $\mathcal{P}' = \Pi(\mathcal{P})$ is a relaxation of the original
problem instance $\mathcal{P}$. In practice, assuming the problem reasons over a set of elements, a clustering algorithm can be used to create clusters of such elements. Then, the aggregate problem instance can be obtained by considering aggregate elements that encompass all elements in a given cluster and by adapting the instance data accordingly. Formally, if a set $E$ of elements is clustered into $K$ clusters, we define two mapping functions: $\Phi: E \to \{0, \dots, K-1\}$ that gives the cluster for each original element and $\Phi^{-1}: \{0, \dots, K-1\} \to 2^E$ that gives the set of original elements for a given cluster.

**TalentSched.** In [17], it is proved that there always exists an optimal solution to the problem in which scenes with the same set of actors are scheduled together. This gives us the opportunity to aggregate the problem by creating $K$ clusters of scenes that require a similar set of actors, which is plausible to occur in real film shoots. Scenes belonging to the same clusters can then be aggregated by taking the intersection of their actor requirements and adding up their durations. Formally, we write $\Pi(\mathcal{P} = (N, A, R, D, C)) = (\Pi_N(N), A, \Pi_R(R), \Pi_D(D), C)$ with $\Pi_N(N) = \{0, \dots, K-1\}$. The aggregate actor requirements are computed as $\Pi_R(R) = R'$ with $R'_i = \cap_{j \in \Phi^{-1}(i)} R_j$ for all $i \in \Pi_N(N)$ and the aggregate durations as $\Pi_D(D) = D'$ with $D'_i = \sum_{j \in \Phi^{-1}(i)} D_j$ for all $i \in \Pi_N(N)$.

**PSP.** The number of item types considered in a PSP instance dramatically impacts the size of the state space – for instance, the case with only one item type can be solved greedily. Therefore, and because it is not unlikely that the machine will produce several sets of similar items, we propose to cluster item types that have similar stocking and changeover costs. The instance aggregation operator is thus $\Pi(\mathcal{P} = (I, S, C, H, Q)) = (\Pi_I(I), \Pi_S(S), \Pi_C(C), H, \Pi_Q(Q))$, where the aggregate set of item types is given by $\Pi_I(I) = \{0, \dots, K-1\}$. Their stocking costs are computed as $\Pi_S(S) = S'$ with $S'_k = \min_{i \in \Phi^{-1}(k)} S_i$ for all $k \in \Pi_I(I)$ and the changeover costs as $\Pi_C(C) = C'$ with $C'_{kl} = \min_{i \in \Phi^{-1}(k), j \in \Phi^{-1}(l)} C_{ij}$ for all $k, l \in \Pi_I(I)$. The aggregate demand matrix is defined as $\Pi_Q(Q) = Q'$ with $Q'_p^k = \sum_{i \in \Phi^{-1}(k)} Q_p^i$. However, as the demand matrix is only supposed to contain unit demands, one must redistribute surplus demands in $Q'$ to the left.

**ALP.** Similarly to the item types of the PSP, the aircraft classes can be aggregated to reduce the complexity of the problem. We thus propose to cluster them based on their minimum separation time with other classes and define the instance aggregation operator as $\Pi(\mathcal{P} = (N, R, C, A, S, T, L)) = (N, R, \Pi_C(C), \Pi_A(A), \Pi_S(S), T, \Pi_L(L))$. The set of aircrafts, their target landing time and the number of runways is kept. The aggregate set of classes is given by $\Pi_C(C) = \{0, \dots, K-1\}$ and their corresponding set of aircrafts is computed as $\Pi_A(A) = A'$ with $A'_i = \cup_{j \in \Phi^{-1}(i)} A_j$ for all $i \in \Pi_C(C)$. The smallest separation times between aggregate classes are kept, as formalized by $\Pi_S(S) = S'$ with $S'_{kl} = \min_{i \in \Phi^{-1}(k), j \in \Phi^{-1}(l)} S_{ij}$ for all $k, l \in \Pi_C(C)$. Finally, the aggregation operator adapts the latest landing times of all the aircrafts so that any aircraft with a given target landing time has a greater latest landing time than all other aircrafts of the same class with a smaller target landing time: $\Pi_L(L) = L'$ with $L'_i = \max\{L_j \mid \Phi(i) = \Phi(j), T_i \le T_j\}$ for all $i \in A$. This property is assumed to hold for the original problem instance, and must be preserved so that aircrafts from the same class can be scheduled sequentially in the DP model.

**► Example 3.** Let us apply the problem instance aggregation to our running example by creating $K = 2$ aggregate scenes. Assuming the following clustering is found: $\Phi(0) = 0, \Phi(1) = 1, \Phi(2) = 0, \Phi(3) = 1$ or equivalently $\Phi^{-1}(0) = \{0, 2\}, \Phi^{-1}(1) = \{1, 3\}$. We thus compute the aggregate scene durations as: $D' = \langle D_0 + D_2, D_1 + D_3 \rangle = \langle 5, 9 \rangle$ and the aggregate actor requirements as: $R' = \langle \{0, 3\} \cap \{0, 2, 3\}, \{0, 1, 3\} \cap \{0, 1, 2, 3\} \rangle = \langle \{0, 3\}, \{0, 1, 3\} \rangle$.
## 3.2 State Aggregation and Lower Bound

A second mapping function accompanies the problem instance aggregation operator: the state aggregation operator $\pi: S \to S'$ that projects each state of the state space $S$ of the original problem in the aggregate state space $S'$. The role of this operator is to translate each original state to its aggregate version by adapting the state information to fit the aggregate problem data. Let us denote by $\mathcal{B}$ and $\mathcal{B}'$ the exact DD for problem $\mathcal{P}$ and $\Pi(\mathcal{P})$, respectively. If the aggregation operators $\Pi$ and $\pi$ are defined such that $v^*(u \rightsquigarrow t | \mathcal{B}) \ge v^*(u' \rightsquigarrow t' | \mathcal{B}')$ for all $u \in \mathcal{B}, u' \in \mathcal{B}'$ with $\pi(\sigma(u)) = \sigma(u')$ and $\pi(\sigma(t)) = \sigma(t')$, then $v^*(u' \rightsquigarrow t' | \mathcal{B}')$ can be used as a lower bound in the original problem, which we will denote by $v_{agg}(\pi(\sigma(u)))$.

Assuming the aggregate problem can be pre-solved exactly and the solution of each subproblem is stored, this aggregation-based lower bound can be retrieved very quickly. One way to exploit it is to incorporate it in the RLB as shown at line 9 of Algorithm 1 so that it is used as often as possible. Another possibility would be to use the aggregate state space to replace the state merging scheme in relaxed DDs. Once a layer with greater width than $W$ is reached, all the states contained in the nodes of the layer could be mapped to the aggregate state space to pursue the compilation in a lower dimensional space.

**TalentSched.** The state compression operator for TalentSched is somewhat complex because we can only map to states where complete aggregate scenes have yet to be scheduled. As a result, if a state $s$ contains scenes in $s.P$ that can optionally be scheduled, we map it to a dummy aggregated state. The same logic is applied when $s.M$ only contains a subset of the scenes that compose an aggregate scene.

$$ \pi(s) = \begin{cases} (\emptyset, \emptyset), & \text{if } s.P \neq \emptyset, \\ (\emptyset, \emptyset), & \text{if } \exists i \in \Pi_N(N) : (\Phi^{-1}(i) \cap s.M) \neq \emptyset \land \Phi^{-1}(i) \not\subseteq s.M, \\ (M', \emptyset), & \text{otherwise, with } M' = \{i \in \Pi_N(N) \mid \Phi^{-1}(i) \subseteq s.M\}. \end{cases} $$

**PSP.** If we extend the definition of $\Phi$ such that $\Phi(\perp) = \perp$, the state aggregation operator can be defined as $\pi(s) = (\Phi(s.i), R)$ with $R_i = \sum_{j \in \Phi^{-1}(i)} s.R_j$ for all $i \in \Pi_I(I)$. The item type is projected to its corresponding aggregate type, and the remaining number of items to produce for each type is separately accumulated within each cluster.

**ALP.** Again, assuming $\Phi(\perp) = \perp$, the state aggregation operator is defined by $\pi(s) = (Q', ROP')$ with the remaining quantities of aircrafts aggregated as $Q'_i = \sum_{j \in \Phi^{-1}(i)} s.Q_j$ for all $i \in \Pi_C(C)$. For the ROP, one only needs to adapt the class of the last aircraft scheduled on each runway $ROP'_i = (s.ROP_0.l, \Phi(s.ROP_0.c))$ for all $i \in R$.

If lower bounds for original states are obtained only by pre-solving the aggregate problem, it is unlikely that the solution of an aggregate subproblem mapped with the state aggregation operator will be available, since the aggregate separation times between aircraft classes lead to very different landing times. However, a lower bound for an aggregate state $s^1 = (Q^1, ROP^1)$ can be provided by the solution of any state $s^2 = (Q^2, ROP^2)$ such that $Q^1 = Q^2$ and $ROP^1.c = ROP^2.c$ and $ROP^1.l \ge ROP^2.l$ for all $i \in R$.

**► Example 4.** Let us compute the aggregation-based lower bound for the root state of the running example $r = (\{0, 1, 2, 3\}, \emptyset)$ given its aggregate version $\pi(r) = (\{0, 1\}, \emptyset)$ and the clustering performed in Example 3. The aggregate version is trivial to solve since the objective function is symmetrical and there are only two scenes to schedule. We thus have $v_{agg}(r) = D'_0 \times (C_0 + C_3) + D'_1 \times (C_0 + C_1 + C_3) = 5 \times (1+4) + 9 \times (1+2+4) = 88$, which is a slightly better lower bound than the one obtained with the relaxed DD of Example 2.
## 3.3 Solution Disaggregation and Node Selection Heuristic

In order to exploit the solution of the aggregate version of a subproblem to find good heuristic solutions for the original subproblem, we need to specify the correspondence between decisions in the aggregate problem with decisions in the original problem. We therefore define a last modeling component, called the *decision disaggregation operator* $\delta(d) : D'_k \to 2^{D_i} \times \dots \times 2^{D_j}$ that maps the instantiation of a variable $x'_k$ in the aggregate problem to a vector of possible corresponding assignments for variables $x_i, \dots, x_j$ in the original problem.

Finally, we define the *path disaggregation operator* that transforms a sequence of decisions in the aggregate problem to a sequence of sets of possible decisions in the original problem: $\Delta(p = (a_k, \dots, a_{n'-1})) = \delta(l(a_k)) \dots \delta(l(a_{n'-1}))$ where $n'$ is the supposed number of aggregate variables and $\cdot$ denotes the concatenation of two vectors. Using this operator, we can compute a score for each decision made during the compilation of restricted DDs. At line 3 of Algorithm 1, we first retrieve the optimal value assignment of the aggregate subproblem and apply the path disaggregation operator on it. Then, a binary score is attributed to each arc at line 15, depending on its compatibility with the disaggregated solution. At line 2 of Algorithm 2, the maximum score obtained along any path up to each node can then be used to order nodes from most to least promising, favoring nodes with incoming paths that are highly compatible with the disaggregated solution. By doing so, the width of restricted DDs is controlled in the same way as before, enabling the preference of solutions even when no feasible solution with the maximum possible score is available.

**TalentSched.** Each aggregate scene corresponds to a set of original scenes, we thus need to map each aggregate decision to a sequence of original decisions: $\delta(i) = V$ where $V_j = \Phi^{-1}(i)$ for all $0 \le j < |\Phi^{-1}(i)|$. It corresponds to any of the scenes from the cluster $i$, duplicated $|\Phi^{-1}(i)|$ times so that they are all scheduled one after another, preferably.

**PSP.** The operator is much simpler to define for the PSP, since each decision concerns the production of one unit of a chosen aggregate item type. It can thus be interpreted as the decision of producing one unit of any item type in the corresponding cluster: $\delta(i) = \langle \Phi^{-1}(i) \rangle$.

**ALP.** The only difference with the PSP is that decisions also contain the runway on which the aircraft is scheduled to land, which remains the same: $\delta(a, r) = \langle \{ (a', r) | a' \in \Phi^{-1}(a) \} \rangle$.

► **Example 5.** As computed in Example 4, the schedule $\langle 0, 1 \rangle$ is optimal for the aggregate problem. By disaggregating this solution, we get $\langle \{0, 2\}, \{0, 2\}, \{1, 3\}, \{1, 3\} \rangle$. We can notice that the optimal schedule $\langle 0, 2, 3, 1 \rangle$ found in Example 1 is compatible with the disaggregated solution and would thus be favored by the aggregation-based node selection heuristic.

# 4 Computational Experiments

The impact of the aggregation-based bounds and heuristics was evaluated experimentally by extending the generic DD-based solver DDO [21] and injecting the modeling of the three discrete optimization problems presented throughout the paper. The version of DDO used includes the improvements introduced in [16, 19]. For each problem, random instances were generated with the following main parameters:

- TalentSched: number of scenes $n \in \{20, 22, 24, 26, 28\}$, number of actors $m \in \{10, 15\}$ and average fraction of actors required for each scene $\rho \in \{0.3, 0.4\}$.
13:14 Boosting DD-Based Branch-And-Bound with Aggregate Dynamic Programming

<figure><img src="image_128.png" /><figcaption>Figure 3 Number of instances solved over time for each configuration and problem studied.</figcaption></figure>

PSP: number of item types $n = 10$, horizon $H \in \{100, 150, 200\}$ and fraction of time periods with a demand $\rho \in \{0.9, 0.95, 1\}$.

ALP: number of aircrafts $n \in \{25, 50, 75, 100\}$, number of runways $r \in \{1, 2, 3, 4\}$, number of aircraft classes $c = 4$ and mean inter-arrival time $40/r$ for generating the target landing times according to a Poisson arrival process.

Furthermore, the instance generation tries to emulate an increasing number of groups of actor requirements, item types and aircraft classes that lend themselves more or less to aggregation. Each instance was presolved in its aggregate state space after aggregating its data according to $k$-means clustering for PSP and ALP and a custom hierarchical clustering for TalentSched that tries to maximize the remaining costs induced by the actor requirements. TalentSched instances can be presolved exactly with 20 aggregate scenes and PSP instances similarly with 5 aggregate item types. On the other hand, not all ALP instances reduced to 2 aggregate aircraft have a reasonable number of states so we employ a relaxed DD with maximum width 10000 for the presolving part instead. Note that the present approach does not compete with the state-of-the-art for TalentSched as it lacks much of the custom symmetry-breaking logic introduced in [17] and similarly for ALP regarding the dominance-breaking constraints presented in [28]. Six different configurations were created by combining the default DD-based solved DDO on one hand and a version using only restricted DDs and no relaxed DDs, denoted rDDO, on the other hand, with the aggregation-based bounds (AggB) and heuristics (AggH). Ten minutes were allotted for each configuration to solve each instance.

Figure 3 presents the cumulative number of instances solved with respect to the solving time. For TalentSched, it appears that any configuration of rDDO performs better than any of DDO. This suggests that the bounds provided by the relaxed DDs are looser than the RLB while being more expensive to compute. It confirms our intuition that the state merging scheme yields bounds with a limited impact for some problems, probably because the state information gets very dilute when many states are merged together. In this case, the RLB computation is also quite involved – see [17]. Still, adding the AggB and the AggH to either configurations improves the results by a small margin, although not that significant. This can be contrasted with the results obtained for the two other problems, which show a clear improvement when the AggB and the AggH are added to either configurations. Furthermore, in cases where rDDO alone yields the worst results, incorporating AggB leads to results that are similar to or better than those achieved by DDO. Combining it with the AggH performs better than DDO in both cases and almost equally well than DDO+AggB+AggH.

The impact of the AggB and the AggH can also be measured in terms of end gap $\frac{UB-LB}{UB}$. Figure 4 compares the end gap obtained for each instance by DDO and DDO+AggB+AggH. It shows that except for a few instances, DDO+AggB+AggH is always closer to terminating the search than DDO, especially for PSP. To validate the relevance of the AggH, we also
<figure><img src="image_129.png" /><figcaption>Figure 4 Comparison of the end gap obtained for each instance by DDO and DDO+AggB+AggH.</figcaption></figure>

<figure><img src="image_130.png" /><figcaption>Figure 5 Comparison of the value of the first solution found by DDO and DDO+AggB+AggH, and of the iteration at which the solution is found for ALP.</figcaption></figure>

compare the value of the first solution found by DDO and DDO+AggB+AggH on Figure 5(a). For TalentSched and PSP, the quality of the first solution is always better when using the AggH. However, there is no clear trend for the ALP. Unlike TalentSched and PSP, for which a solution is always found at the first iteration, the landing time windows of ALP make it difficult to find a feasible solution. This explains both the end gaps close to one in Figure 4 and the $\infty$ values in Figure 5(a), which represent the absence of a feasible solution. We thus compare on Figure 5(b) the iteration at which the first solution is found. We observe that DDO+AggB+AggH finds a feasible solution much earlier than DDO in most cases. This showcases well the benefits of a node selection heuristic with a more global awareness.

# 5 Conclusion

This paper explained how ideas from aggregate dynamic programming can be incorporated in DD-based optimization solvers. We proposed to derive lower bounds and node selection heuristics from a pre-solved aggregate version of the original problem at hand, and explained how these can be seamlessly added to the DD-based optimization framework. Computational experiments on three different problems showed that they provide lower bounds that further strengthen the current approach, and that could even be used as a replacement for relaxed DDs in some cases. Furthermore, the aggregation-based node selection heuristics were shown very valuable as they manage to steer the compilation of relaxed DDs toward better solutions earlier in the search. When applying this idea to a highly constrained problem, the heuristics proved to quickly lead to feasible solutions that were hard to find otherwise. These results suggest that aggregation-based bounds and heuristics capture global problem structures well, as opposed to the greedy MinLP heuristic traditionally used to compile approximate DDs.
13:16 Boosting DD-Based Branch-And-Bound with Aggregate Dynamic Programming

References

1 Henrik Reif Andersen, Tarik Hadzic, John N Hooker, and Peter Tiedemann. A constraint store based on multivalued decision diagrams. In *International Conference on Principles and Practice of Constraint Programming*, pages 118–132. Springer, 2007.

2 Sven Axsäter. State aggregation in dynamic programming—an application to scheduling of independent jobs on parallel processors. *Operations Research Letters*, 2(4):171–176, 1983.

3 James C Bean, John R Birge, and Robert L Smith. Aggregation in dynamic programming. *Operations Research*, 35(2):215–220, 1987.

4 Richard Bellman. The theory of dynamic programming. *Bulletin of the American Mathematical Society*, 60(6):503–515, November 1954. URL: https://projecteuclid.org:443/euclid.bams/1183519147.

5 David Bergman and Andre A. Cire. On finding the optimal bdd relaxation. In Domenico Salvagnin and Michele Lombardi, editors, *Integration of AI and OR Techniques in Constraint Programming*, volume 10335 of LNCS, pages 41–50. Springer, 2017.

6 David Bergman, Andre A Cire, Ashish Sabharwal, Horst Samulowitz, Vijay Saraswat, and Willem-Jan van Hoeve. Parallel combinatorial optimization with decision diagrams. In *Integration of AI and OR Techniques in Constraint Programming: 11th International Conference, CPAIOR 2014, Cork, Ireland, May 19-23, 2014. Proceedings 11*, pages 351–367. Springer, 2014.

7 David Bergman, Andre A Cire, Willem-Jan van Hoeve, and John N Hooker. Variable ordering for the application of bdds to the maximum independent set problem. In *International conference on integration of artificial intelligence (AI) and operations research (OR) techniques in constraint programming*, pages 34–49. Springer, 2012.

8 David Bergman, Andre A Cire, Willem-Jan van Hoeve, and John N Hooker. Optimization bounds from binary decision diagrams. *INFORMS Journal on Computing*, 26(2):253–268, 2014.

9 David Bergman, Andre A Cire, Willem-Jan van Hoeve, and John N Hooker. Discrete optimization with decision diagrams. *INFORMS Journal on Computing*, 28(1):47–66, 2016.

10 David Bergman, Andre A Cire, Willem-Jan van Hoeve, and Tallys Yunes. Bdd-based heuristics for binary optimization. *Journal of Heuristics*, 20(2):211–234, 2014.

11 Quentin Cappart, David Bergman, Louis-Martin Rousseau, Isabeau Prémont-Schwarz, and Augustin Parjadis. Improving variable orderings of approximate decision diagrams using reinforcement learning. *INFORMS Journal on Computing*, 34(5):2552–2570, 2022.

12 Quentin Cappart, Emmanuel Goutierre, David Bergman, and Louis-Martin Rousseau. Improving optimization bounds using machine learning: Decision diagrams meet deep reinforcement learning. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 33, pages 1443–1451, 2019.

13 Margarita P Castro, Andre A Cire, and J Christopher Beck. An mdd-based lagrangian approach to the multicommodity pickup-and-delivery tsp. *INFORMS Journal on Computing*, 32(2):263–278, 2020.

14 Margarita P Castro, Chiara Piacentini, Andre Augusto Cire, and J Christopher Beck. Solving delete free planning with relaxed decision diagram based heuristics. *Journal of Artificial Intelligence Research*, 67:607–651, 2020.

15 Edmund M Clarke, Orna Grumberg, and David E Long. Model checking and abstraction. *ACM transactions on Programming Languages and Systems (TOPLAS)*, 16(5):1512–1542, 1994.

16 Vianney Coppé, Xavier Gillard, and Pierre Schaus. Decision diagram-based branch-and-bound with caching for dominance and suboptimality detection, 2023. arXiv:2211.13118.

17 Maria Garcia de la Banda, Peter J Stuckey, and Geoffrey Chu. Solving talent scheduling with dynamic programming. *INFORMS Journal on Computing*, 23(1):120–137, 2011.
18 Xavier Gillard. *Discrete optimization with decision diagrams: design of a generic solver, improved bounding techniques, and discovery of good feasible solutions with large neighborhood search*. PhD thesis, UCL-Université Catholique de Louvain, 2022.

19 Xavier Gillard, Vianney Coppé, Pierre Schaus, and André Augusto Cire. Improving the filtering of branch-and-bound mdd solver. In *International Conference on Integration of Constraint Programming, Artificial Intelligence, and Operations Research*, pages 231–247. Springer, 2021.

20 Xavier Gillard and Pierre Schaus. Large neighborhood search with decision diagrams. In *International Joint Conference on Artificial Intelligence*, 2022.

21 Xavier Gillard, Pierre Schaus, and Vianney Coppé. Ddo, a generic and efficient framework for mdd-based optimization. In *Proceedings of the Twenty-Ninth International Conference on International Joint Conferences on Artificial Intelligence*, pages 5243–5245, 2021.

22 Jaime E Gonzalez, Andre A Cire, Andrea Lodi, and Louis-Martin Rousseau. Integrated integer programming and decision diagram search tree with an application to the maximum independent set problem. *Constraints*, pages 1–24, 2020.

23 John N. Hooker. Improved job sequencing bounds from decision diagrams. In Thomas Schiex and Simon de Givry, editors, *Principles and Practice of Constraint Programming*, volume 11802 of LNCS, pages 268–283. Springer, 2019.

24 Matthias Horn, Johannes Maschler, Günther R Raidl, and Elina Rönnberg. A\*-based construction of decision diagrams for a prize-collecting scheduling problem. *Computers & Operations Research*, 126:105125, 2021.

25 Alan J. Hu. *Techniques for efficient formal verification using binary decision diagrams*. PhD thesis, Stanford University, Department of Computer Science, 1995.

26 Anthony Karahalios and Willem-Jan van Hoeve. Variable ordering for decision diagrams: A portfolio approach. *Constraints*, 27(1):116–133, 2022.

27 C.-Y. Lee. Representation of switching circuits by binary-decision programs. *The Bell System Technical Journal*, 38(4):985–999, 1959.

28 Alexander Lieder, Dirk Briskorn, and Raik Stolletz. A dynamic programming approach for the aircraft landing problem with aircraft classes. *European Journal of Operational Research*, 243(1):61–69, 2015.

29 Shin-ichi Minato. *Binary decision diagrams and applications for VLSI CAD*, volume 342. Springer Science & Business Media, 1995.

30 Isaac Rudich, Quentin Cappart, and Louis-Martin Rousseau. Peel-And-Bound: Generating Stronger Relaxed Bounds with Multivalued Decision Diagrams. In Christine Solnon, editor, *28th International Conference on Principles and Practice of Constraint Programming (CP 2022)*, volume 235 of *Leibniz International Proceedings in Informatics (LIPIcs)*, pages 35:1–35:20. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2022.

31 Christian Tjandraatmadja. *Decision Diagram Relaxations for Integer Programming*. PhD thesis, Carnegie Mellon University Tepper School of Business, 2018.

32 Christian Tjandraatmadja and Willem-Jan van Hoeve. Target cuts from relaxed decision diagrams. *INFORMS Journal on Computing*, 31(2):285–301, 2019. doi:10.1287/ijoc.2018.0830.
<figure><img src="image_131.png" /></figure>
# Fast Matrix Multiplication Without Tears:
A Constraint Programming Approach

Arnaud Deza<fnref n="1" /> ☐

Department of Mechanical and Industrial Engineering, University of Toronto, Canada

Chang Liu<fnref n="1" /> ☐

Department of Mechanical and Industrial Engineering, University of Toronto, Canada

Pashootan Vaezipoor ☐

Department of Computer Science, University of Toronto, Canada

Elias B. Khalil ☐

Department of Mechanical and Industrial Engineering, University of Toronto, Canada

## Abstract

It is known that the multiplication of an $N \times M$ matrix with an $M \times P$ matrix can be performed using fewer multiplications than what the naive $NMP$ approach suggests. The most famous instance of this is Strassen's algorithm for multiplying $2 \times 2$ matrices in 7 instead of 8 multiplications. This gives rise to the constraint satisfaction problem of fast matrix multiplication, where a set of $R < NMP$ multiplication terms must be chosen and combined such that they satisfy correctness constraints on the output matrix. Despite its highly combinatorial nature, this problem has not been exhaustively examined from that perspective, as evidenced for example by the recent deep reinforcement learning approach of AlphaTensor. In this work, we propose a simple yet novel Constraint Programming approach to find algorithms for fast matrix multiplication or provide proof of infeasibility otherwise. We propose a set of symmetry-breaking constraints and valid inequalities that are particularly helpful in proving infeasibility. On the feasible side, we find that exploiting solver performance variability in conjunction with a sparsity-based problem decomposition enables finding solutions for larger (feasible) instances of fast matrix multiplication. Our experimental results using CP Optimizer demonstrate that we can find fast matrix multiplication algorithms for matrices up to $3 \times 3$ with $R = 23$ in a short amount of time.

**2012 ACM Subject Classification** Mathematics of computing

**Keywords and phrases** fast matrix multiplication, computer-assisted proofs, constraint programming, constraint satisfaction problem

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

## Supplementary Material

Software (Source Code): https://github.com/khalil-research/Matrix-Mult-CP/tree/main

# 1 Introduction

Matrix multiplication is a fundamental operation in linear algebra with applications in virtually every computational domain. As a result, extensive research has been dedicated to the development of faster matrix multiplication algorithms.

The elementary way of multiplying two $N \times N$ matrices requires $N^3$ multiplications. For example, multiplying two $2 \times 2$ matrices naively requires a total of $2^3 = 8$ multiplications. In 1969, Strassen [16] constructed an algorithm that finds the product of two $2 \times 2$ matrices in only 7 multiplications. This discovery has had significant implications as it opened up

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

<footnote n="1">¹ These authors contributed equally.</footnote>
14:2 Fast Matrix Multiplication Without Tears: A Constraint Programming Approach

the door for potentially faster algorithms for large-scale matrix or tensor computations. Strassen's algorithm has later been proved to be both canonical [4] (no smaller rank exists) and essentially unique [6] (all other solutions of the same rank are equivalent up to symmetry).

Currently, the best-known algorithm for multiplying $3 \times 3$ matrices requires $R = 23$ multiplications, compared to the naive elementary method that requires 27 multiplications. A known theoretical lower bound of $R = 19$ exists [2], however, it remains unclear whether $19 \le R \le 22$ is truly attainable. This is a testament to the difficulty of the *fast matrix multiplication* (FMM) problem, which has been intractable for existing methods even for tiny matrices.

In the literature, the general approach to finding FMM algorithms starts by representing matrix multiplication as a tensor operation using the multiplication tensor $T_N$ followed by finding exact or approximate low-rank decompositions that represent $T_N$. The factor matrices that are used in the low-rank decomposition encode FMM algorithms. A rank-7 decomposition (i.e., a multiplication algorithm that uses 7 multiplication operations) of a $2 \times 2$ matrix multiplication using Strassen's algorithm is shown in Figure 1. Existing methods for finding such factor matrices have several limitations. The most successful and common methods include local search [13] techniques for low-rank approximation, which cannot guarantee optimality. A more recent successful approach [7] searches for low-rank decomposition using *reinforcement learning* (RL) and was successful in finding faster algorithms for $N = 4$. However, this method is not exhaustive and hence cannot prove the infeasibility of a given rank.

In this work, we propose a novel approach to finding FMM algorithms by formulating the tensor decomposition problem, for the first time, as a *constraint satisfaction problem* (CSP) that is solved using *Constraint Programming* (CP). We believe that this is a very natural formulation of this highly combinatorial problem. CP is advantageous for FMM in that it is a flexible framework that can bring to bear a wide range of search and logical inference techniques that have been developed over the last few decades. It provides the ability to prove infeasibility when it is not possible to multiply two matrices using a given number of multiplications.

Besides a base CP formulation for FMM, we propose a set of symmetry-breaking constraints and valid inequalities that are useful for infeasibility proofs. On the feasible side, we show that “performance variability” w.r.t. solver random seeds can be exploited in conjunction with a sparsity-based decomposition of FMM for faster solving. Our experimental results, while limited to matrices of size up to $3 \times 3$, demonstrate the effectiveness of the aforementioned constraints and techniques. The CP approach to FMM is uniquely positioned to close open questions such as whether it is possible to multiply two $3 \times 3$ matrices in 19 to 22 multiplications. While we do not yet resolve this or other open questions, our work opens up the potential for further enhancements to the CP formulation and search such as customized branching strategies and CP-based heuristics.

2 Fast Matrix Multiplication: Problem Statement

The multiplication of two matrices $A$ and $B$ of sizes $N \times M$ and $M \times P$, respectively, results in a product matrix $C$ of size $N \times P$. This operation can be represented by a binary third-order tensor $T_{NMP}$ ($T_N$ for square matrices $A$ and $B$ of size $N \times N$). An entry $T_{i,j,k}$ of this tensor is equal to 1 if and only if the $k^{th}$ entry in the output matrix $C$ uses the scalar product of the $i^{th}$ entry of $A$ and the $j^{th}$ entry of $B$. Here, $i, j$, and $k$ are indices of a matrix entry starting with 1 in the first row and column; and proceeding entry by entry, left to right, top
to bottom. For example, for $N = M = P = 2$, it must be that $T_{2,3,1} = 1$ since the first entry of $C$, $c_1$, is equal to $a_1b_1 + a_2b_3$. Similarly, $T_{1,2,1} = 0$ must hold since $a_1b_2$ is not part of $c_1$. Figures 1a and 1b show a complete example of the indexing and tensor representation.

The FMM problem for a given tensor $T_{NMP}$, rank $R \in \mathbb{Z}^+$, and field $\mathbb{F}$ (e.g., $\mathbb{F} = \{-1, 0, 1\}$) asks: can each entry $T_{i,j,k}$ of $T_{NMP}$ be expressed as the sum of exactly $R$ trilinear terms involving the factor matrices $U \in \mathbb{F}^{N \cdot M \times R}$, $V \in \mathbb{F}^{M \cdot P \times R}$, and $W \in \mathbb{F}^{N \cdot P \times R}$, as follows:

$$T_{i,j,k} = \sum_{r=1}^{R} U_{i,r} \cdot V_{j,r} \cdot W_{k,r} \quad \forall i \in \{1, \dots, N \cdot M\}, j \in \{1, \dots, M \cdot P\}, k \in \{1, \dots, N \cdot P\}$$

Note that we use the notation $\mathbb{F}^{L \times Q}$ to refer to the set of matrices of dimension $L \times Q$ and entries in $\mathbb{F}$. The CSP is to find factor matrices with entries in $\mathbb{F}$ that produce the tensor $T_{NMP}$ for a given rank $R$.

This decomposition is also referred to as the polyadic decomposition and its associated rank is the minimal $R$ needed. The rank can be interpreted as the number of multiplications required to compute the product. For example, for $2 \times 2$ matrices, the rank of the decomposition using Strassen's algorithm is 7. Figure 1 walks through an example of the low-rank decomposition of a $2 \times 2$ matrix multiplication using Strassen's algorithm. The matrix multiplication of the two $2 \times 2$ matrices can be seen in Figure 1a, its associated tensor representation $T_N$ in Figure 1b, the low-rank decomposition in Figure 1c, and the factor matrices $U$, $V$, and $W$ in Figure 1d.

# 3 Related Work

Since Strassen's discovery [16], there has been substantial research on finding faster algorithms for matrix multiplication. Mathematicians have discovered such algorithms manually over the years for a variety of matrix dimensions and ranks. In this section, however, we will focus on automated methods for discovering such algorithms and briefly discuss some of the existing methods. A recent survey on the topic can be found in [3].

## 3.1 Continuous Local Search Methods

The most common approach in the literature to compute the factor matrices $U$, $V$, and $W$ is to use (heuristic, continuous) local search methods for low-rank tensor decomposition. The state-of-the-art local method [13] uses alternating least squares with regularization. This method has been the most successful in finding fast algorithms whilst remaining computationally tractable and has been scaled up to $N = M = P = 4$, $R = 49^2$. However, this approach has limitations which include getting stuck at local minima, facing ill-conditioned linear least-squares problems, and solutions being only adequate up to machine precision. Additionally, these methods are not exhaustive and hence cannot be used to provide a proof of infeasibility for a given rank $R$.

## 3.2 AlphaTensor

More recently, DeepMind released AlphaTensor [7], a deep RL method that searches this large combinatorial space by playing a single-player game, the TensorGame, formulated as a Markov decision process (MDP). At every step $t$ of this MDP, the state is characterized

<footnote n="2">² Note that this particular result is not very useful as an $R = 49$ solution can be obtained by applying Strassen's $R = 7$ algorithm for $2 \times 2$ matrices on the four $2 \times 2$ blocks of the $4 \times 4$ matrices.</footnote>
14:4 Fast Matrix Multiplication Without Tears: A Constraint Programming Approach

$$ \begin{pmatrix} c_1 & c_2 \\ c_3 & c_4 \end{pmatrix} = \begin{pmatrix} a_1 & a_2 \\ a_3 & a_4 \end{pmatrix} \cdot \begin{pmatrix} b_1 & b_2 \\ b_3 & b_4 \end{pmatrix} $$

(a) Multiplication of two 2 × 2 matrices. We highlight the term $c_1 = a_1b_1 + a_2b_2$.

$$ T_{:,;,1} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad T_{:,;,2} = \begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad T_{:,;,3} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{pmatrix}, \quad T_{:,;,4} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} $$

(b) Tensor representation of the 2 × 2 matrix multiplication operation. $T_{:,;,1}$ represents $c_1$, the entry $T_{2,3,1}$ (in yellow) is set to 1 because the product $a_2b_3$ is required to compute $c_1$ (similarly for $T_{1,1,1}$ in red).

$$ m_1 = (a_1 + a_4)(b_1 + b_4) $$

$$ m_5 = (a_1 + a_2)(b_4) $$

$$ m_2 = (a_3 + a_4)(b_1) $$

$$ m_6 = (a_3 - a_1)(b_1 + b_2) $$

$$ m_3 = (a_1)(b_2 - b_4) $$

$$ m_7 = (a_2 - a_4)(b_3 + b_4) $$

$$ m_4 = (a_4)(b_3 - b_1) $$

$$ c_1 = m_1 + m_4 - m_5 + m_7 $$

$$ = (a_1 + a_4)(b_1 + b_4) + (a_4)(b_3 - b_1) - (a_1 + a_2)(b_4) + (a_2 - a_4)(b_3 + b_4) $$

$$ = a_1b_1 + a_1b_4 + a_4b_1 + a_4b_4 + a_4b_3 - a_4b_1 - a_1b_4 - a_2b_4 + a_2b_3 + a_2b_4 - a_4b_3 - a_4b_4 $$

$$ = a_1b_1 + a_2b_3 $$

$$ c_2 = m_3 + m_5 $$

$$ c_3 = m_2 + m_4 $$

$$ c_4 = m_1 - m_2 + m_3 + m_6 $$

(c) A low-rank decomposition of the 2 × 2 matrix multiplication using Strassen's algorithm. The $m$ terms are the multiplication terms and the $c$ terms represent the entries in the product matrix. Here $c_1 = m_1 + m_4 - m_5 + m_7$ gives $c_1 = a_1b_1 + a_2b_3$ after expansion.

$$ U = \begin{pmatrix} m_1 & m_2 & m_3 & m_4 & m_5 & m_6 & m_7 \\ 1 & 0 & 1 & 0 & 1 & -1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & -1 \end{pmatrix}, \quad \begin{array}{l} a_1 \\ a_2 \\ a_3 \\ a_4 \end{array} \\ \\ V = \begin{pmatrix} 1 & 1 & 0 & -1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 1 \\ 1 & 0 & -1 & 0 & 1 & 0 & 1 \end{pmatrix}, \quad \begin{array}{l} b_1 \\ b_2 \\ b_3 \\ b_4 \end{array} \\ \\ W = \begin{pmatrix} 1 & 0 & 0 & 1 & -1 & 0 & 1 \\ 0 & 0 & 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 1 & -1 & 1 & 0 & 0 & 1 & 0 \end{pmatrix}, \quad \begin{array}{l} c_1 \\ c_2 \\ c_3 \\ c_4 \end{array} $$

(d) The factor matrices $U$, $V$, and $W$ for Strassen's algorithm. The columns in $U$ and $V$ represent the coefficient of the $a$ and $b$ terms in each $m$. Each row in $W$ represents the coefficient of the $m$ terms in one $c$ term.

■ **Figure 1** A low-rank decomposition of a 2 × 2 matrix multiplication using Strassen's algorithm.
by a tensor $S_t$ which is initially set to the target multiplication tensor, i.e., $S_0 = T_N$. An action $a_t$ at iteration $t$ corresponds to the player selecting a triplet of vectors $(u^{(t)}, v^{(t)}, w^{(t)})$ which in turn will provide the next state $S_t = S_{t-1} - u^{(t)} \otimes v^{(t)} \otimes w^{(t)}$ where $\otimes$ denotes the outer tensor product. The goal of the player is to reach the zero tensor $S_t = \mathbf{0}$ in the fewest number of steps possible. This is done by providing a reward of $-1$ to the player after every non-terminal state whereas a large negative reward $-\gamma(S_{R_{\text{limit}}})$ is given to the player if the number of steps $R_{\text{limit}}$ is met, where $\gamma(S_{R_{\text{limit}}})$ upper bounds the rank of the tensor at iteration $R_{\text{limit}}$. If the agent successfully reaches the zero tensor, the sequence of actions taken constitutes a valid low-rank decomposition of $T_N$, and hence an FMM algorithm is found with the rank $R$ corresponding to the number of steps taken by the agent.

This approach is the first to directly incorporate learning into the search which resulted in the discovery of new minimal ranks for certain non-trivial cases. The largest case tackled by this method is $N = M = P = 5$, $R = 98$. The sole focus of this purely heuristic method is to find lower ranks than currently best-known ranks but it cannot prove the infeasibility of a given rank. Additionally, rather complex architectures and multiple training phases were required for successful learning. It is worth noting that AlphaTensor was trained for one week on 64 Tensor Processing Units (TPUs), Google's proprietary chip. The paper [7] does not provide any estimates of the amount of computation required to produce the reported results, namely how long the trained "agent" must be run to discover FMM algorithms. Our CP runs use much fewer resources while leveraging thread parallelism in the CP solver on readily-available CPU machines.

## 3.3 Integer Programming

The work that is the most related to our approach tackles this problem through a *mixed-integer linear program* (MILP) formulation in an unpublished technical report [14]. The goal of this methodology is to linearize the trilinear products in the low-rank decomposition of $T_N$ to a MILP that aims to 1) maximize the sparsity of the integer decision variables representing factor matrices $U, V$, and $W$ and 2) minimize the reconstruction loss (L1 norm) from the input $T_N$ and the multiplication tensor attained by the decision variables representing factor matrices. The report [14] focuses solely on presenting the MILP formulation for square matrices but does not include any computational experiments. However, the MILP formulations for $N \in \{2, 3\}$ are benchmark problems in MIPLIB 2017 [9]<fnref n="3" />. The linearization of the trilinear products likely leads to a weak linear programming relaxation as well as an explosion in the number of integer variables and constraints, which might explain why the MILP approach to FMM has not picked up significant interest. A CP formulation is more natural and compact, as we will show in this paper.

## 3.4 Classical AI Planning

Very recently, AI planning techniques were used for FMM [15]. They use a similar state space as AlphaTensor but use various planning tools (with and without exhaustive search) to solve this problem. They compared a number of heuristic and exact planning methods from the literature on matrices of size up to $3 \times 3$. However, the experiments show that planning approaches are severely limited, even failing to find Strassen's algorithm for the $2 \times 2$ case (see Table 1 in [15]). We will show that our CP approach is significantly more effective as we are able to attack the $3 \times 3$ case with $R = 23$, matching the known upper bound from the literature.

<footnote n="3">³ See https://miplib.zib.de/instance_details_fastxgemm-n3r21s3t6.html for example.</footnote>
## 3.5 Boolean Satisfiability-based Approaches

Boolean Satisfiability (SAT) formulations to solve the FMM problem were first proposed in 2011 [5]. More specifically, the authors studied the multiplication of 3 × 3 matrices and were able to find a new general algorithm with rank 23 with help from SAT solver tweaks and improvements in a few days with one CPU. More recently, Heule et al. improved on the SAT-based methods for the multiplication of 3 × 3 matrices and have successfully found many thousands new general algorithms with rank 23 [10, 11]. In their SAT-formulation, the authors transformed the multiplication to ‘and’ clauses and the addition to ‘xor’ clauses, then a Tseistin transformation is performed to formulate the algebraic FMM problem into a SAT problem.

The authors propose a random pairing method where value assignments are initialized based on observed results from previously known solutions with streamlining constraints to find new general solutions; then local search is used to find additional solutions. The SAT-based formulations of the FMM problem demonstrates to be very difficult for complete SAT solvers, thus making it difficult to provide proof of infeasibility.

# 4 Constraint Programming for Fast Matrix Multiplication

In the FMM problem, all variables have the same domain $F = \{-1, 0, 1\}^4$. Since the variable domains are small and this problem is highly structured, CP is a promising solution paradigm.

The base CP model for FMM is given in Equation (1). Let $\mathcal{U}$ denote the set $\{1, \dots, N \cdot M\}$, $\mathcal{V}$ denote the set $\{1, \dots, M \cdot P\}$, $\mathcal{W}$ denote the set $\{1, \dots, N \cdot P\}$, and $\mathcal{R}$ denote the set $\{1, \dots, R\}$. The CP model uses three sets of variables: $u_{i,r}$ where $i \in \mathcal{U}$, $v_{j,r}$ where $j \in \mathcal{V}$, and $w_{k,r}$ where $k \in \mathcal{W}$; $r \in \mathcal{R}$ in all three cases. Each variable $u_{i,r}$, $v_{j,r}$ and $w_{k,r}$ represents the value of the $i/j/k^{\text{th}}$ row and $r^{\text{th}}$ column of the matrices $U, V$, and $W$. The domain of all variables is $\{-1, 0, 1\}$. The set of constraints presented here requires that the decomposition algorithm’s output matches the original tensor multiplication $T_{NMP}$. Therefore the input to the CSP model is 4 integers: $(N, M, P)$ and $R$. The model then reads as:

$$ \sum_{r \in \mathcal{R}} (u_{i,r} \cdot v_{j,r} \cdot w_{k,r}) = T_{i,j,k}, \quad \forall i \in \mathcal{U}, j \in \mathcal{V}, k \in \mathcal{W} \\ u_{i,r}, v_{j,r}, w_{k,r} \in \{-1, 0, 1\}, \quad \forall i \in \mathcal{U}, j \in \mathcal{V}, k \in \mathcal{W}, r \in \mathcal{R} \quad (1) $$

The search space for this (NP-complete) problem grows very quickly with increasing matrix sizes $N, M, P$ and rank $R$. With only one set of equality constraints, a CP solver may struggle with constraint propagation, thus failing to scale with increasing $N, M, P$. To that end, we will introduce additional valid constraints to help CP prune and propagate more efficiently.

## 4.1 Symmetry Breaking

There are many symmetric solutions to the FMM problem. We can reduce the search space of our problem significantly by prohibiting symmetries.

<footnote n="4">⁴ One can consider bigger fields such as $\{-2, -1, 0, 1, 2\}$ but the bulk of the work in the literature has been with $\{-1, 0, 1\}$.</footnote>
### 4.1.1 Permutation Symmetry

Since addition is commutative, i.e., $(a_1 + a_2) = (a_2 + a_1)$, there are many equivalent solutions to the tensor decomposition problem. Therefore, any permutation of the columns of matrices $U$, $V$, and $W$ produces an equivalent solution. If we consider Strassen's solution for the $2 \times 2$ case, Figure 2 provides an example of two equivalent solutions.

$$
\begin{array}{l}
\text{sol1: } U = \begin{pmatrix} 1 & 0 & 1 & 0 & 1 & -1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & -1 \end{pmatrix} \quad V = \begin{pmatrix} 1 & 1 & 0 & -1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 1 \\ 1 & 0 & -1 & 0 & 1 & 0 & 1 \end{pmatrix} \quad W = \begin{pmatrix} 1 & 0 & 0 & 1 & -1 & 0 & 1 \\ 0 & 0 & 1 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & -1 & 1 & 0 & 0 & 1 & 0 \end{pmatrix} \\
\\
\text{sol2: } U = \begin{pmatrix} -1 & 0 & 0 & 0 & 1 & 1 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 1 & -1 & 0 & 1 & 0 \end{pmatrix} \quad V = \begin{pmatrix} 1 & -1 & 1 & 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 1 \end{pmatrix} \quad W = \begin{pmatrix} 0 & 1 & 0 & 1 & 0 & 1 & -1 \\ 0 & 0 & 0 & 0 & 1 & 1 & 1 \\ 0 & 1 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & -1 & 0 & 1 & 1 & 0 \end{pmatrix}
\end{array}
$$

■ **Figure 2** Two equivalent solutions for Strassen's solution of $2 \times 2$ matrix multiplication. sol2 is the lexicographic-strict presentation of this solution.

In order to break this symmetry, we introduce a lexicographic-strict<fnref n="5" /> constraint on the $u_{i,r}$ and $v_{j,r}$ variables. When applied to two variable arrays $x$ and $y$, the lexicographic ordering constraint enforces that $x$ is strictly less than $y$ in the defined lexicographic order. Because of the strictness, this also enforces that the two variable arrays must be different. This set of symmetry-breaking constraints is modelled as follows:

$$ \text{lexicographic-strict}([u_{:,r}; v_{:,r}], [u_{:,r+1}; v_{:,r+1}]), \quad \forall r \in \mathcal{R} $$

where $[u_{:,r}; v_{:,r}]$ represents the vector concatenating the $r$th column of the matrix $U$ and $V$. In Figure 2, sol2 satisfies the lexicographic-strict constraint.

### 4.1.2 Sign Symmetry

For the multiplicative $m_i$ terms, one can easily see that multiplying both sets of terms from $A$ and $B$ by $-1$ will result in the same solution. For example, $(a_1 + a_4)(b_1 + b_4) = (-a_1 - a_4)(-b_1 - b_4)$, where we could multiply any subset of columns of $U$ and $V$ by $-1$ to achieve the same solution. We call this symmetry the sign symmetry. In order to break it, we introduce the following constraints:

$$
\begin{array}{l}
u_{1,r} \le 0 \\
u_{i,r} \le \displaystyle\sum_{i'=1}^{i-1} |u_{i',r}| \qquad \forall r \in \mathcal{R}, i > 1, i \in \mathcal{U}
\end{array}
$$

The main idea of these constraints is to enforce that the first non-zero entry in a column of $U$ can only take on the value of $-1$, enforcing that the first entry of the columns is either $0$ or $-1$. The subsequent constraints ensure that for any column $r$, an entry in row $i > 1$ can only be $1$ if there has been an entry in the same column in an earlier row with value $-1$. This set of constraints applies to the concatenation of the columns in $U$ and $V$, however, in modelling, it only needs to be applied to the columns of the $U$ matrix as none of the columns can be zero, so the leading $-1$ must appear in the $U$ matrix. By applying these constraints, we make sure that $\{-u_{i,r}\}$ is infeasible for any feasible $\{u_{i,r}\}$. Employing this sign symmetry breaking constraint to sol2 from Figure 2, we arrive at sol3 shown in Figure 3.

<footnote n="5">⁵ https://www.ibm.com/docs/en/icos/22.1.0?topic=variables-lexicographic-constraint</footnote>