# LAD2025, A constraint-based solver for the subgraph isomorphism problem

Christine Solnon

► To cite this version:

Christine Solnon. LAD2025, A constraint-based solver for the subgraph isomorphism problem. Artificial Intelligence (AIJ), 2026, 352, pp.104474. 10.1016/j.artint.2025.104474. hal-05464475

HAL Id: hal-05464475

https://hal.science/hal-05464475v1

Submitted on 19 Jan 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# LAD2025, a Constraint-Based Solver for the Subgraph Isomorphism Problem

Christine Solnon, CITI, Inria, INSA Lyon, F-69621 Villeurbanne
christine.solnon@insa-lyon.fr

December 21, 2025

## Abstract

The Subgraph Isomorphism Problem (SIP) is an $NP$-complete problem that aims at finding a copy of a pattern graph in a target graph. It may be modelled as a constraint satisfaction problem in a very straightforward way, and exact approaches for solving SIPs usually propagate constraints to reduce the search space. In particular, PathLAD is a solver introduced in 2016 that combines Locally All Different (LAD) constraints with path-based supplemental constraints. In this paper, we introduce LAD2025, which combines a complete refactoring of PathLAD with new features: new supplemental constraints, a weight-based variable ordering heuristic, random restarts with nogood recording, a new value ordering heuristic and a rule for selecting the level of filtering.

## 1 Introduction

The *Subgraph Isomorphism Problem* (SIP) aims at finding a copy of a pattern graph in a target graph. It has many applications as graphs are widely used to model objects such as molecules, biochemical reaction networks, images, 3D objects, or social networks, for example. It is also more general than classical $NP$-complete graph problems such as searching for a clique of order *k* or a Hamiltonian path.

The SIP may be modelled as a *Constraint Satisfaction Problem* (CSP) in a very straightforward way, and exact approaches for solving SIPs usually propagate constraints to reduce the search space. Different constraint propagation algorithms may be considered. Some of them are stronger than others, i.e., they prune more branches, but they also have higher time complexities. In [Sol19], four classical exact approaches are compared: VF2 [CFSV04], RI [BG17], Glasgow [MPT20], and PathLAD [KMS16]. This experimental study shows us that RI, which performs very fast and simple propagations, is the fastest approach on many easy instances, but it is not able to solve many other instances within a reasonable amount of time. PathLAD, which performs very strong and expensive propagations, is able to solve harder instances than RI but it is less efficient on very easy instances. Glasgow offers a good compromise between these two extremes, and it is the most successful approach when runtimes are limited to 1000 seconds. More recently, Wang et al [WJCL23, WJC24] have introduced PathLAD+ which improves PathLAD with new features, thus allowing it to solve more instances than Glasgow when runtimes are limited to one hour.

In this paper, we introduce LAD2025, which combines a complete refactoring of PathLAD with new heuristics, pruning rules, and search techniques, some of them being classical ingredients of *Constraint Programming* (CP) solvers, and some others being new. We show that LAD2025 is more successful than Glasgow and PathLAD+.
**Contributions and overview of the paper.** In Section 2, we define the SIP and describe the main existing exact approaches for solving it. In Section 3, we describe benchmarks and performance measures used in our experiments.

In Section 4, we describe how PathLAD has been refactored, mostly based on lessons learned from its analysis reported in [Knu25], and from the implementation of the Global All Different (GAD) constraint reported in [GMN08].

In Section 5, we show how to tighten domains by exploiting supplemental labels. These labels are constructed by enumerating paths, as proposed in [ALM$^+$14] and [MP15], and by computing maximum cliques in vertex neighbourhoods, as proposed in [KM21]. We also consider new labels obtained by enumerating all cliques up to a given order. These labels are used to filter domains before starting the search. They are also used to detect edge incompatibilities, thus allowing us to strengthen constraint propagation during the search.

In Section 6, we experimentally compare different variable ordering heuristics and we show that the weight-based ordering heuristic introduced in [BHLS04] strongly improves the solution process.

In [ADH$^+$19], it is shown that random restarts combined with nogood recording and a new randomised value ordering heuristic allow Glasgow to find solutions quicker. In Section 7, we introduce a new value ordering heuristic, that borrows features from phase savings used in SAT solvers [PD07], and that may be combined with the randomized value ordering heuristic introduced in [ADH$^+$19].

In Section 8, we introduce a simple rule for selecting the filtering level.

In Section 9, we experimentally compare LAD2025 with other exact solvers.

In Section 10, we summarise our contributions and discuss further work.

# 2 Background

## 2.1 Notations and preliminary definitions on graphs

The cardinality of a set $S$ is denoted $\#S$. Given two integer values $i$ and $j$, we denote $[i, j]$ the set of integer values ranging from $i$ to $j$ (this set is empty whenever $i > j$).

A graph is defined by an ordered pair $(V, E)$ such that $V$ is a set of vertices and $E \subseteq V \times V$ is a set of unordered pairs of vertices called edges. Two vertices $u$ and $v$ are said to be adjacent if $\{u, v\} \in E$. The neighbourhood of a vertex $u$, denoted $\mathcal{N}(u)$, is the set of its adjacent vertices, i.e., $\mathcal{N}(u) = \{v \in V \mid \{u, v\} \in E\}$. The degree of a vertex $u$, denoted $d^\circ(u)$, is the cardinality of its neighbourhood, i.e., $d^\circ(u) = \#\mathcal{N}(u)$. The density of a graph is $\frac{2\#E}{\#V(\#V-1)}$.

A clique in a graph $(V, E)$ is a set of vertices $S \subseteq V$ that are all pairwise adjacent, i.e., $\forall \{u, v\} \subseteq S$, $\{u, v\} \in E$. The order of a clique is its cardinality. A clique is maximum if there does not exist a clique in the graph with a larger order. We denote $\omega(G)$ the order of the maximum clique in a graph $G$, and $\omega(u)$ the order of the largest clique of $G$ that contains a vertex $u$ of $G$.

A $k$-path is a sequence of $k + 1$ vertices such that each vertex is adjacent to the next vertex in the sequence.

Throughout this paper, we consider two graphs: a pattern graph denoted $G = (V, E)$, and a target graph denoted $G' = (V', E')$. Without loss of generality, we assume that $V \cap V' = \emptyset$. Vertices in $V$ are called pattern vertices and are denoted $u$ or $v$ (possibly subscripted) whereas vertices in $V'$ are called target vertices and are denoted $u'$ or $v'$ (possibly subscripted).

We note $n$ the number of pattern vertices, $m$ the number of pattern edges, and $d$ the largest pattern vertex degree, i.e., $n = \#V$, $m = \#E$, and $d = \max_{u \in V} d^\circ(u)$. Similarly, $n'$, $m'$ and $d'$ denote the number of vertices, edges and largest degree in the target graph.
## 2.2 Definition of the SIP

Given $G = (V, E)$ and $G' = (V', E')$, the goal of the SIP is to find an injective mapping function $f: V \to V'$ that associates a different target vertex to each pattern vertex and such that each pattern edge is mapped to a target edge, i.e., $\forall \{u, v\} \in E, \{f(u), f(v)\} \in E'$. $f$ is said to be a sub-isomorphism function. Note that we consider the non-induced SIP such that two non-adjacent pattern vertices may be mapped to two adjacent target vertices (in other words, non-adjacency needs not to be preserved by the sub-isomorphism function $f$). We discuss the extension of LAD2025 to the induced SIP in Section 10.

## 2.3 Modelling the SIP as a CSP

The SIP may be modelled as a CSP in a straightforward way. For each pattern vertex $u \in V$, we introduce a variable $x_u$ which represents the target vertex that is mapped to $u$. The domain of $x_u$ is the set of target vertices that may be mapped to $u$ and it is denoted $D_u$.

Graphs may contain loops, i.e., edges whose endpoints are identical. In this case, a pattern vertex $u \in V$ with a loop must be mapped to a target vertex with a loop and, therefore, its domain $D_u$ is limited to target vertices with loops. Note that pattern vertices without loops may be mapped to target vertices with loops or not so that their domains are not limited to target vertices without loops.

The following two constraints ensure that the mapping is injective and that edges are preserved, respectively.

**Constraint 1 (injectivity)** $\forall \{u, v\} \subseteq V, x_u \neq x_v$

**Constraint 2 (edge preservation)** $\forall \{u, v\} \in E, \{x_u, x_v\} \in E'$

## 2.4 Constraint-based approaches for the SIP

In [Ull76], Ullman introduces an SIP algorithm that contains key ingredients of CP solvers even though CP was not yet known at this time. For each pattern vertex $u \in V$, the algorithm maintains a set of candidate target vertices that may be mapped to $u$ and this set actually corresponds to the domain $D_u$. These domains are filtered by a refinement procedure which ensures the Domain Consistency of Constraint 2. Let us first recall the definition of Domain Consistency (also known as Generalised Arc Consistency) [Knu25].

**Definition 1 (Domain Consistency (DC))** Let $c$ be a constraint defined over a tuple of variables $(x_1, \dots, x_k)$, and let $D_i$ be the domain of $x_i, \forall i \in [1, k]$. The constraint $c$ is DC if for each $i \in [1, k]$ and every value $v \in D_i$, there exists a supporting tuple, i.e., a tuple $(v_1, \dots, v_{i-1}, v, v_{i+1}, \dots, v_k)$ that satisfies $c$ and such that $v_j \in D_j, \forall j \in [1, k] \setminus \{i\}$.

Hence, the refinement procedure introduced by Ullman ensures that, for every pattern edge $\{u, v\} \in E$ and every target vertex $u' \in D_u$, there exists a target vertex $v' \in D_v$ such that $\{u', v'\} \in E'$.

This first CP-based algorithm for SIP may be improved by introducing stronger propagations. In particular, we may use the GAD propagation algorithm introduced in [Reg94] to ensure the DC of Constraint 1.

Also, we may filter domains by propagating invariant properties, i.e., properties that are satisfied by all sub-isomorphism functions. For example, for any sub-isomorphism function $f$ and any pattern vertex $u$, we have $d^\circ(u) \le d^\circ(f(u))$, because every neighbour of $u$ must be mapped to a different neighbour of $f(u)$. This invariant property is used to remove from the domain of $x_u$ every target vertex whose degree is smaller than the degree of $u$.
This degree-based invariant property may also be exploited to ensure stronger consistencies. More precisely, let $u \in V$ be a pattern vertex and let $u' \in D_u$ be a target vertex that may be mapped to $u$.

* In [LV02], the degree of $u$ is ensured to be smaller than or equal to the number of target vertices that are both adjacent to $u'$ and belong to domains of vertices adjacent to $u$, i.e., $d^\circ(u) \le \#(\mathcal{N}(u') \cap (\cup_{v \in \mathcal{N}(u)} D_v))$;

* Iterative Label Filtering (ILF) [ZDS10] ensures that the number of neighbours of $u$ with a given label $l$ is smaller than or equal to the number of neighbours of $v$ with a label compatible with $l$, where labels are initially defined by degrees and are iteratively refined by considering neighbour labels until reaching a fixed point or a given limit on the number of iterations;

* Locally All Different (LAD) filtering [Sol10] ensures that each neighbour $v \in \mathcal{N}(u)$ can be mapped to a different target vertex $v' \in \mathcal{N}(u') \cap D_v$ (this is explained with more details in Section 4).

We show in [Sol10] that LAD filtering is as strong as ILF when there is no limit on the number of iterations, whereas it is strictly stronger than the filtering of [LV02].

In [ALM$^{+}$14], LAD filtering is generalised in order to propagate other kinds of invariant properties among a subset of vertices, such as the number of $k$-paths between two vertices. This is defined by means of supplemental graphs in Glasgow [MP15], and by means of supplemental labels in [Knu25]. These supplemental labels are described in Section 5. In [KMS16], path-based supplemental labels are added to LAD, and the resulting solver is called PathLAD.

Glasgow is a state-of-the-art CP-based approach initially introduced in [MP15]. It considers a cheaper all-different propagator than the GAD algorithm of [Reg94] and it relies on specialised bit-parallel data structures. It also exploits path-based supplemental labels. Glasgow has been further improved in [ADH$^{+}$19] by introducing a new randomised value ordering heuristic which is combined with restarts and nogood recording (described in Section 7), and in [KM21] by exploiting new clique-based supplemental labels (described in Section 5).

In [WJCL23, WJC24], PathLAD has been improved by adding a new value ordering heuristic, an adaptative propagation method (described in Section 8), and a probing procedure used at the beginning of the search. This new approach, called PathLAD+, solves more instances than Glasgow within a time limit of one hour.

## 2.5 Lighter approaches for the SIP

CP-based approaches described in the previous section maintain a domain $D_u$ for each pattern vertex $u \in V$ and propagate constraints to filter these domains and detect inconsistencies whenever a domain becomes empty. VF2 [CFSV04] is a lighter approach that does not fully maintain domains: it maintains domains only for pattern vertices that are adjacent to vertices that are already mapped, and it uses simple propagation rules to filter these domains. RI [BG17] further improves this approach by introducing ordering heuristics.

In [Sol19], we report an experimental comparison of PathLAD, Glasgow, VF2, and RI on a wide benchmark described in the next section. On this benchmark, RI clearly outperforms VF2, and Glasgow clearly outperforms PathLAD. However, RI and Glasgow are rather complementary: RI is faster than Glasgow on many easy instances, but it is not able to solve many other instances that are easily solved by Glasgow, and it performs very poorly on hard instances. Hence, Glasgow is able to solve more instances than RI within a runtime limit of 1000 seconds.
2.6 Integer Linear Programming (ILP) Models for the SIP

An ILP model has been introduced in [LBLA+09]. This model introduces a binary variable $x_{uu'}$ for each couple of vertices $(u, u') \in V \times V'$, and a binary variable $y_{uvu'v'}$ for each couple of edges $(\{u, v\}, \{u', v'\}) \in E \times E'$. In [MdG23], this model is refined, and other models are also introduced. Experimental results reported in [MdG23] show us that none of these ILP models is competitive with Glasgow.

3 Experimental settings

3.1 Benchmarks

In this paper, we consider the same set of benchmarks as in [KMS16], and all these benchmarks
are available at http://perso.citi-lab.fr/csolnon/SIP.html.

LV is the benchmark introduced in [LV02], based on the Stanford GraphBase which contains various kinds of graphs described in [Knu93]. LV contains 3,767 instances obtained by considering all pairs $(G, G')$ such that $G$ (resp. $G'$) is one of the 50 (resp. 112) first graphs of the GraphBase (excluding the first graph which has no edge), and $G'$ has at least as many vertices and edges as $G$.

images contains 6,278 instances introduced in [DSdlH+11] and 24 instances introduced in [SDdlHJ15]. All graphs are built from segmented images: vertices correspond to segmented regions, and edges to adjacency relations between these regions.

meshes contains 3,018 instances introduced in [DSdlH+11]. Graphs are built from meshes that model 3D objects.

rand contains 1000 instances with randomly generated graphs: bounded valence graphs and 4D meshes introduced in [DSFSV03], and scale-free graphs introduced in [ZDS10]. All instances but 20 are feasible by construction: the target graph is randomly generated, and then a pattern is extracted from the target. The 20 infeasible instances contain scale-free graphs such that edges have been added to the pattern after its extraction.

randER contains 270 instances with randomly generated Erdös-Rényi graphs introduced in [DSFSV03]. All instances are feasible by construction. These instances are separated from the other *rand* instances because they are harder to solve.

phase contains 200 instances with randomly generated Erdös-Rényi graphs. However, unlike in *randER*, these instances are not feasible by construction, and they are expected to be really hard because they are close to the phase transition, as described in [MPST18].

In Figure 1, we display the number of vertices and edges of the pattern and target graphs for each instance. Some instances are very small: the smallest has 10 vertices in both the pattern and the target graph. Some others are rather large: the largest pattern graph (coming from *rand*) has 900 vertices, and the largest target graph (coming from *LV*) has 6,671 vertices.

Instances of LV have very different characteristics, as they are generated from graphs that have very different properties: some graphs are very sparse, whereas some others are very dense; some graphs have uniformly distributed edges whereas some others exhibit structures such as cliques, for example. All other benchmarks are rather homogeneous, i.e., instances within a same benchmark are quite similar (though they may have different sizes).
<figure><img src="image_2.png" /><figcaption>Figure 1: Number of vertices (left) and edges (right) of SIP instances. Each point $(x, y)$ corresponds to an instance such that $x = n$ and $y = n'$ on the left plot, and $x = m$ and $y = m'$ on the right plot. The colour of the point depends on the benchmark.</figcaption></figure>

In [Sol19], we report results of Glasgow, PathLAD, RI, and VF2 on these benchmarks, and we show that many instances are very easy. In particular, all instances of *images*, and also 900 instances of *rand*, 3,510 instances of *LV*, 2,914 instances of *meshes*, and 203 instances of *randER*, are solved by the four solvers within a time limit of one second. However, some other instances are very hard and cannot be solved by any of the four solvers within 1000 seconds: 120 instances of *LV*, 11 instances of *meshes*, and 83 instances of *phase*. Some of these unsolved instances are rather small. For example, an unsolved instance in *LV* (called g19-g20) has only 64 vertices in both the pattern and the target graphs, and there are 168 pattern edges and 728 target edges (corresponding to all valid knight and queen moves on an 8 × 8 chessboard, respectively [Knu25]). We refer the reader to [Sol19] for more details on these hard instances, and to [MPST18] for an extensive study of the relation between instance hardness and graph density for randomly generated graphs.

A challenge is to improve solvers on hard instances while not degrading performance on easy instances.

## 3.2 Performance measures

The source code of LAD2025 is available at https://perso.citi-lab.fr/csolnon/LAD.html.

All experiments are performed on a single core of a CPU Intel Xeon Gold 6444Y with 16 cores of 3.60 GHz. Given a solver *s* and an instance *i*, $t_i^s$ denotes the CPU time in seconds used by *s* to solve *i* (i.e., either find a solution or prove the inconsistency of *i*). Given a run time limit *l*, a solver *s*, and a benchmark *B*, $S_l^s = \{i \in B \mid t_i^s \le l\}$ (resp. $U_l^s = \{i \in B \mid t_i^s > l\}$) denotes the set of instances of *B* that are solved (resp. not solved) by *s* within the time limit *l*. We consider the following performance measures:

* the number of instances not solved by *s* within the time limit *l*, i.e., $\#U_l^s$;

* $t_l^s = \frac{\sum_{i \in S_l^s} t_i^s}{\#S_l^s}$, the average solving time for the instances solved by *s* within the time limit *l*.

* $\bar{t}_l^s = \frac{l * \#U_l^s + \sum_{i \in S_l^s} t_i^s}{\#B}$, a lower bound of the average solving time obtained by bounding the run time of each unsolved instance with *l*.

When there is no ambiguity, $\#U_l^s$ is noted $\#U$ and $t_l^s$ is noted $t$.
We consider two different time limits: $l = 100$ seconds during preliminary experiments used to evaluate the impact of the new features introduced to boost PathLAD; $l \in \{100, 2000\}$ seconds during the final comparison with other state-of-the-art solvers. To consider other time limits than 100s and 2000s, we also display cumulative plots, that display the evolution of the cumulative number of instances solved by a solver $s$ with respect to time, i.e., the function $f(l) = \#S_l^s$.

To evaluate the degree of dispersion and skewness of run times for each benchmark and each solver, we use box-plots to visualise minimum and maximum values as well as first, second, and third quartiles. We also display on these box-plots $t_l^s$ and $\bar{t}_l^s$.

To compare two solvers $a$ and $b$ for a given time limit $l$ such that $\#S_l^a \ge \#S_l^b = k$, we report the Aggregate Speed-Up (ASU) introduced in [HMN$^+$18] and defined by $ASU_l(\frac{b}{a}) = \frac{\text{argmin}_l \#S_l^b = k}{\text{argmin}_l \#S_l^a = k}$. In other words, $ASU_l(\frac{b}{a})$ measures the ratio between the smallest time limit that allows $b$ to solve $k$ instances and the smallest time limit that allows $a$ to also solve $k$ instances. Note that the sets of $k$ instances solved by each solver may be different.

To allow a per instance comparison of two solvers $a$ and $b$ for a given time limit $l$, we also display scatter plots that contain a point $(x, y)$ for each instance $i$ such that $x = t_i^a$ and $y = t_i^b$: every point below (resp. above) the line $x = y$ corresponds to an instance that is more efficiently solved by $b$ (resp. $a$). Whenever $t_i^a > l$ (resp. $t_i^b > l$), we set $x$ to $l$ (resp. $y$ to $l$). This allows us to visualise instances not solved by $a$ (resp. $b$) by looking at points on $x = l$ (resp. $y = l$).

# 4 Refactoring of PathLAD (V0)

In Section 4.1, we describe the LAD filtering algorithm introduced in [Sol10]. In [GMN08], Gent et al discuss implementation issues and report experimental evaluations of different key components of propagation algorithms for the allDifferent constraint. Also, the algorithm for ensuring the DC of LAD constraints is analysed in [Knu25], and Knuth has proposed some improvements for its implementation in personal communications. The most important improvements are described in Sections 4.2 to 4.5. In Section 4.6, we experimentally compare the old implementation of PathLAD with its refactoring (denoted V0).

## 4.1 Description of LAD filtering

For each pattern vertex $u$ and each target vertex $u' \in D_u$, the algorithm introduced in [Sol10] ensures the DC of a constraint denoted $LAD_{uu'}$ and defined as follows.

**Constraint 3 (LAD$_{uu'}$)** $x_u = u' \Rightarrow [\forall v \in \mathcal{N}(u), x_v \in \mathcal{N}(u')] \wedge \text{allDifferent}(\{x_v : v \in \mathcal{N}(u)\})$

To ensure the DC of $LAD_{uu'}$, we must remove $u'$ from $D_u$ whenever it is not possible to map every neighbour $v$ of $u$ to a different neighbour of $u'$ that belongs to the domain of $v$. This is implemented by searching for a covering matching in the bipartite graph $G_{uu'} = (V_{uu'}, E_{uu'})$ such that:

* $V_{uu'}$ associates a vertex with every neighbour of $u$ or $u'$, i.e., $V_{uu'} = \mathcal{N}(u) \cup \mathcal{N}(u')$;

* $E_{uu'}$ links every pattern vertex $v$ in the neighbourhood of $u$ with every target vertex $v'$ that both belongs to the neighbourhood of $u'$ and to the domain of $v$, i.e., $E_{uu'} = \{(v, v') \in \mathcal{N}(u) \times \mathcal{N}(u') | v' \in D_v\}$.

A covering matching of $G_{uu'}$ is a subset of edges $S \subseteq E_{uu'}$ such that (i) each pattern vertex $v \in \mathcal{N}(u)$ occurs in exactly one edge of $S$, and (ii) each target vertex $v' \in \mathcal{N}(u')$ occurs in at most one edge of $S$. Whenever such a covering matching does not exist, we remove $u'$ from the domain of $x_u$ to ensure the DC of $LAD_{uu'}$.
269

270

271

272

273

274

Each time a value $u'$ is removed from the domain of a variable $x_u$, we must propagate LAD$_{vv'}$ for every pattern vertex $v$ adjacent to $u$ and every target vertex $v' \in D_v$ adjacent to $u'$. Hence, we maintain a list $L$ of couples $(u, u')$ such that LAD$_{uu'}$ must be propagated and we iteratively remove a couple from $L$, ensure the DC of its associated LAD constraint, and add new constraints to $L$ whenever some domains have been reduced, until $L$ becomes empty (all LAD constraints are DC), or a domain has been wiped out (an inconsistency has been detected).

275

## 4.2 Computation of covering matchings

276

277

278

279

PathLAD propagates $O(n \cdot n')$ LAD constraints, and also one GAD constraint that ensures that the $x_u$ variables have all different values. In both cases, this implies to search for a covering matching in a bipartite graph (the bipartite graph for GAD associates an edge with every couple $(u, u') \in V \times V'$ such that $u' \in D_u$).

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

Covering matchings are incrementally computed. More precisely, for each bipartite graph (associated with a LAD$_{uu'}$ constraint or with the GAD constraint) we memorise the couples $(v, v')$ that belong to its last computed covering matching. Each time we have to search for a covering matching of a bipartite graph, we start the search from the initial matching that contains all memorised couples $(v, v')$ such that $v' \in D_v$, and we complete this initial matching by searching for augmenting paths in the residual graph. When using the algorithm of Ford and Fulkerson [FF56] with a breadth-first search to search for a covering matching (denoted BFS-FF), a single augmenting path is computed at each iteration and, therefore, the optimal matching is computed in $O(k \cdot e)$ where $k$ is the number of unmatched pattern vertices and $e$ is the number of edges in the bipartite graph ($e$ is upper bounded by $d \cdot d'$ for LAD constraints and by $n \cdot n'$ for the GAD constraint). When using the algorithm of Hopcroft and Karp [HK73] (denoted HK), several augmenting paths are computed in a single search and the optimal matching is computed in $O(\sqrt{k \cdot e})$. In [GMN08], it is experimentally shown that BFS-FF is never slower than HK, and it can be up to 30% faster than HK. Hence, we have replaced HK with BFS-FF in our new implementation. We have empirically noticed that $k$ is usually rather small, and this may explain why BFS-FF performs better than HK in practice.

296

## 4.3 Computation of strongly connected components

297

298

299

300

301

302

303

304

To ensure the DC of GAD, we must not only check that there exists a covering matching in a bipartite graph, but we must also remove every target vertex $u'$ from the domain of any variable $x_u$ whenever $x_u$ cannot be assigned to $u'$ without violating GAD. To do this, we must search for strongly connected components in the residual graph, as shown by Berge in [Ber73]. In our old implementation, we have used the algorithm of Kosaraju [CLRS09], which performs two depth-first searches. The algorithm of Tarjan introduced in [Tar72] has the same time complexity as Kosaraju's algorithm, but it is more efficient in practice as it performs only one depth-first search. Hence, we use Tarjan's algorithm in our new implementation.

305

## 4.4 Data structures for representing sets

306

307

308

309

310

We use the data structure described in [GMN08] to represent sets used to mark vertices (during BFS-FF, for example). The idea is to use an integer $c$ (called the certificate) and an array $a$ indexed by vertices: a vertex $v$ belongs to the set iff $a[v] = c$. Initially, all elements of $a$ are initialised to 0, and $c$ is initialised to 1, so that the set is empty. To add a vertex $v$ in the set we simply assign $c$ to $a[v]$. To remove all vertices from the set in constant time, we increment $c$.

311

312

313

Our old implementation already used sparse sets to represent domains: this allows us to restore a domain in constant time when backtracking [LCdSMSSL13]. In our new implementation, we also use a sparse set to maintain the current set of unassigned variables.
<figure><img src="image_3.png" /><figcaption>Figure 2: Scatter plot to compare PathLAD (x axis) and V0 (y axis), for a time limit of 100s.</figcaption></figure>

## 4.5 Bimodal search

For each pattern vertex $u \in V$ and each target vertex $u' \in D_u$, the bipartite graph $G_{uu'}$ used to ensure the DC of $LAD_{uu'}$ is not explicitly built, but it is inferred from $D_u$, $G$, and $G'$. More precisely, for each pattern vertex $v \in \mathcal{N}(u)$, we compute the set $S_v = D_v \cap \mathcal{N}(u')$ to obtain the set of edges that start from $v$ in $G_{uu'}$. To compute $S_v$, we take advantage of the fact that domains are represented with sparse sets and that edges of $G'$ are represented with both an adjacency matrix and adjacency lists. This allows us to (i) check in constant time if a value is present in $D_v$ or in $\mathcal{N}(u')$, and (ii) iterate on all values in $D_v$ or in $\mathcal{N}(u')$ in linear time. Hence, depending on whether $D_v$ contains more values than $\mathcal{N}(u')$ or not, we either iterate on $\mathcal{N}(u')$ and remove values that do not belong to $D_v$ or iterate on $D_v$ and remove values that do not belong to $\mathcal{N}(u')$. This is called bimodal search in [LBCBP$^{+}$25]. This bimodal search is already implemented in the first version of LAD and in PathLAD, and it is also implemented in LAD2025.

In [LBCBP$^{+}$25], bimodal search is also used for propagating GAD, both when searching for a covering matching and when computing strongly connected components. This implies using backtrackable doubly linked lists (called *dancing links* in [Knu22]) to manage the list of unvisited vertices. Experimental results reported in [LBCBP$^{+}$25] show us that this greatly speeds-up the propagation of GAD, especially for large sets of variables. We plan to investigate this technique in further work.

## 4.6 Experimental evaluation

Let V0 denote the new implementation of PathLAD as described in this section. We compare run times of the old implementation of PathLAD and V0 with a scatter plot displayed in Figure 2. Both implementations use the same ordering heuristics and ensure the same consistency at each node so that they build the same search tree. If V0 spends more time than PathLAD on some easy instances, V0 is faster for all instances for which PathLAD needs more than 0.1s, and there are 122 instances that are not solved by PathLAD within a time limit of 100s whereas they are solved by V0. The aggregate speed-up between PathLAD and V0 for a time limit of 100s is $ASU_{100}(\frac{PathLAD}{V0}) = 9.86$.
## 5 Supplemental labels (V1)

Supplemental labels are integer values which are associated with vertices (unary labels) or couples of vertices (binary labels).

### Definition 2 (Supplemental label)

* A unary supplemental label is a function $l: V \cup V' \to \mathbb{N}$ such that $l(u) \le l(f(u))$ for any sub-isomorphism function $f$ and any pattern vertex $u \in V$.

* A binary supplemental label is a function $l: V \times V \cup V' \times V' \to \mathbb{N}$ such that $l(u, v) \le l(f(u), f(v))$ for any sub-isomorphism function $f$ and any couple of pattern vertices $(u, v) \in V \times V$.

Supplemental labels are used to define incompatibility relationships which are used to prune the search space: a pattern vertex $u$ is incompatible with every target vertex $u'$ such that $l(u) > l(u')$ and a couple of pattern vertices $(u, v)$ is incompatible with every couple of target vertices $(u', v')$ such that $l(u, v) > l(u', v')$.

Supplemental labels may be computed by searching for occurrences of "motifs" in both $G$ and $G'$, and we define these motif-based labels in Section 5.1. In Section 5.2, we describe all supplemental labels used in LAD2025. In Section 5.3, we show how to use supplemental labels to tighten domains. In Section 5.4, we experimentally evaluate the interest of using supplemental labels.

### 5.1 Definition of motif-based supplemental labels

The authors of [ALM$^+$14] and [MP15] exploit the fact that sub-isomorphism functions preserve distances and paths to provide additional filtering procedures. This idea is generalised by means of supplemental graphs in [MPT20] and by means of motif-based supplemental labels in [Knu25]. Let us first define what is a motif.

**Definition 3 (motif)** A motif is a connected graph $m = (V_m, E_m)$ such that $V_m$ contains two special vertices, noted a and b, i.e., $\{a, b\} \subseteq V_m$.

Given a graph $G = (V, E)$ and a motif $m$, we define unary and binary labels as follows.

### Definition 4 (motif-based supplemental labels)

* For each vertex $u \in V$, the unary label $m(u)$ is the number of sub-isomorphism functions $g: V_m \to V$ such that $g(a) = u$.

* For each couple of vertices $(u, v) \in V \times V$, the binary label $m(u, v)$ is the number of sub-isomorphism functions $g: V_m \to V$ such that $g(a) = u$ and $g(b) = v$.

As edges are preserved by sub-isomorphism functions, we know that $m(u) \le m(f(u))$ and $m(u, v) \le m(f(u), f(v))$ for every pattern vertices $u, v \in V$, every motif $m$, and every sub-isomorphism function $f$ from $G$ to $G'$ (see [Knu25] for more details).

For example, the motif $P_2$ displayed in Fig. 3 is a 2-path whose endpoints are a and b. When computing labels associated with $P_2$ in the graph $G$, we have $P_2(1) = 5$ because there are five 2-paths starting from vertex 1. In the graph $G'$, we have $P_2(5') = 4$ and, therefore, vertex 1 of $G$ is not compatible with vertex 5' of $G'$. Also, we have $P_2(2, 4) = 2$ because there are two 2-paths from vertex 2 to vertex 4. As $P_2(2', 4') = 1$, the couple of vertices (2, 4) of $G$ is not compatible with the couple of vertices (2', 4') of $G'$.
<figure><img src="image_4.png" /><figcaption>Figure 3: A motif $P_2$, a pattern graph $G$ and a target graph $G'$.</figcaption></figure>

## 5.2 Supplemental labels considered in LAD2025

In LAD2025, we consider two different kinds of motif-based labels, $P_2$ and $K_i$, which are described in the next two paragraphs. We also consider a supplemental label which is not based on a motif, denoted $\omega$ and described in a third paragraph.

**$P_2$ labels.** The motif $P_2$ displayed in Fig. 3 is used in LAD2025 as it offers a good-tradeoff between the time needed to compute labels and the strength of the resulting pruning. The computation of $P_2$ labels is done in $O(p \cdot d)$ in $G$ (resp. $O(p' \cdot d')$ in $G'$) by visiting every neighbour of $v$ for every edge $\{u, v\} \in E$.

In Glasgow, the motif $P_2$ is also used to provide additional filtering, but this is done in a slightly different way. Before starting the search, two supplemental graphs $G_i$ and $G'_i$ are built for each $i \in [1, 4]$, where $G_i$ (resp. $G'_i$) contains the same set of vertices as $G$ (resp. $G'$), and an edge between two vertices $u$ and $v$ whenever $P_2(u, v) \ge i$ (resp. $u'$ and $v'$ whenever $P_2(u', v') \ge i$). Then, Glasgow searches for a mapping $f$ which is a sub-isomorphism function not only for $G$ and $G'$ but also for $G_i$ and $G'_i$, for all $i \in [1, 4]$ simultaneously. Limiting the number of couples of supplemental graphs $(G_i, G'_i)$ to $i = 4$ speeds-up the propagation, but it also filters fewer values in some cases.

**$K_i$ labels.** We also compute motif-based labels associated with cliques. More precisely, let $K_i$ be the motif which is a clique of order $i$ and where $a$ and $b$ are two different vertices of this clique. Given a parameter $C \ge 3$, we compute supplemental labels associated with every motif $K_i$ with $i \in [3, C]$. As far as we know, this is the first time that such supplemental labels are considered.

A key point is to enumerate all cliques efficiently. Let us first note that we cannot use existing algorithms for efficiently enumerating all maximal cliques, and then update labels of the nodes and edges of each maximal clique because a same non maximal clique may be included in several maximal cliques. Hence, given a graph $G$ and a maximum clique order $C$, $K_i$ labels are recursively computed as described in Algorithm 1.

<caption>**Algorithm 1:** computeLabels($G, C, S$)</caption>

**Input:** A graph $G = (V, E)$, the maximum clique order $C$, and a clique $S \subseteq V$ such that $\#S \le C$
1 &nbsp;&nbsp;if $\#S \ge 3$ then
2 &nbsp;&nbsp;&nbsp;&nbsp;for each vertex $u \in S$ do increment $K_{\#S}(u)$;
3 &nbsp;&nbsp;&nbsp;&nbsp;for each edge $\{u, v\} \subseteq S$ do increment $K_{\#S}(u, v)$;
4 &nbsp;&nbsp;if $\#S < C$ then
5 &nbsp;&nbsp;&nbsp;&nbsp;for each $u \in V$ such that $\forall v \in S, u \in \mathcal{N}(v) \land u > v$ do computeLabels($G, C, S \cup \{u\}$);

At the first call to computeLabels($G, C, S$), the clique $S$ is empty, and all vertex and edge labels have been initialised to 0. At each recursive call, the labels associated with the vertices and edges
of S are incremented (Lines 1-3). Then, if S has not reached the size limit, a vertex of G that is adjacent to every vertex in S is added to S and computeLabels is recursively called (Line 5). To break symmetries, we only add vertices that are greater than all vertices in S (with respect to a given total order on the set of vertices). We do not maintain a list of candidate vertices that can be added to S. Instead, we use a dedicated data structure to access in constant time to each neighbour of u that is larger than u, where u is the largest vertex in S (i.e., the last added one). To speed-up the computation, we upper bound C to the largest value k such that $\exists u \in V, K_k(u) > 0$ when enumerating cliques in the target graph $G'$: It is useless to enumerate cliques of order $i > k$ in $G'$ as $K_i$ labels of pattern and edge vertices are always equal to 0.

**$\omega$ labels.** Another kind of supplemental label is introduced in [KM21] and is based on maximum cliques in vertex neighbourhoods. More precisely, for each vertex $u \in V$ (resp. $u' \in V'$), the order of the largest clique that contains this vertex, i.e., $\omega(u)$ (resp. $\omega(u')$), is a supplemental label. Indeed, as edges are preserved by sub-isomorphism functions, we have $\omega(u) \le \omega(f(u))$ for every sub-isomorphism function $f$ and every pattern vertex $u \in V$.

To compute $\omega$ labels, we have to compute for each vertex the maximum clique that contains it. As proposed in [KM21], we do not search for cliques larger than $\omega(G)$ in the target graph $G'$ and set $\omega(u')$ to $\omega(G)$ whenever the order of the largest clique that contains $u'$ exceeds $\omega(G)$.

We could also consider similar labels for couples of vertices, i.e., compute for each edge $\{u, v\}$ the order of the largest clique that contains this edge. However, this is more expensive to compute. Indeed, a first possibility is to compute a maximum clique for each edge of $G$ and $G'$, but this is prohibitively long for very dense graphs. Another possibility is to compute the order of the largest clique that contains $\{u, v\}$ when computing $\omega(u)$ and $\omega(v)$, but this prevents us from using bounding functions (that are used to prune parts of the search space that cannot contain cliques of order greater than the largest clique found so far), as the largest clique that contains both $u$ and $v$ may be smaller than $\omega(u)$ and $\omega(v)$.

## 5.3 Using supplemental labels to tighten domains

$\omega$ labels are unary, whereas the motifs $P_2$ and $K_i$ are used to generate both unary and binary labels. Unary labels are used to filter domains before starting the search: for each pattern vertex $u \in V$, we remove from $D_u$ every target vertex $u'$ such that $l(u) > l(u')$ for some label $l \in \{\omega, P_2, K_i | i \in [3, C]\}$.

Binary labels may be used to define the following binary constraint, for each $(u, v) \in V \times V$ such that $P_2(u, v) \ge 1$ or $K_i(u, v) \ge 1$ for some $i \in [3, C]$.

**Constraint 4** $(x_u, x_v) \in \{(u', v') \in V' \times V' | P_2(u, v) \le P_2(u', v') \wedge \forall i \in [3, C], K_i(u, v) \le K_i(u', v')\}$ This constraint is propagated by Glasgow for $P_2$ motifs, for example (using supplemental graphs and up to $i=4$ occurrences of $P_2$ as explained in Section 5.2).

The authors of [ALM+14] define a stronger global constraint called SND. This constraint is parameterised by a neighbourhood function which associates a subset of vertices with every vertex. The SND constraint associated with a neighbourhood function $N$, a pattern vertex $u \in V$, a target vertex $u' \in D_u$, and a binary label $l$ is defined as follows.

**Constraint 5 (SND<sub>uu'</sub><sup>l,N</sup>)**
$x_u = u' \Rightarrow [\forall v \in N(u), x_v \in N(u')] \wedge [l(u, v) \le l(u', x_v)] \wedge \text{allDifferent}(\{x_v | v \in N(u)\})$

This constraint is a generalisation of LAD<sub>uu'</sub>: LAD<sub>uu'</sub> is equivalent to SND<sub>uu'</sub><sup>l,N</sup> when $N$ is defined by vertex adjacency (i.e., $N = \mathcal{N}$) and all pairs of vertices have the same supplemental label (i.e., $\forall u, v \in V, \forall u', v' \in V', l(u, v) = l(u', v')$). To ensure the DC of SND<sub>uu'</sub><sup>l,N</sup>, we have to search for a covering matching in a bipartite graph that contains an edge $\{v, v'\}$ for every pattern vertex $v \in N(u)$ and every target vertex $v' \in N(u') \cap D_v$ such that $l(u, v) \le l(u', v')$.

Two specific neighbourhoods are considered in [ALM+14]:
451 • $\mathcal{N}$ is the neighbourhood which associates to $u$ all its adjacent vertices;

452 • all is the neighbourhood which associates to $u$ all vertices but $u$, i.e., $\text{all}(u) = V \setminus \{u\}$.

453 When supplemental labels are computed by enumerating a motif $m$, the choice of $N$ in $\{\mathcal{N}, \text{all}\}$ highly depends on whether the two special vertices $a$ and $b$ are connected by an edge or not in $m$.

455 Neighbourhood function considered when $\{a, b\}$ is an edge of $m$. Let us first consider the case where $a$ and $b$ are connected by an edge, which is the case when $m = K_i$, for each $i \in [3, C]$. In this case, the neighbourhood all does not allow us to filter more values than the neighbourhood $\mathcal{N}$ because, for any vertices $u, v \in V$ such that $v$ is not adjacent to $u$, $K_i(u, v) = 0$ and, therefore, $(u, v)$ is compatible with any couple $(u', v') \in V' \times V'$ according to $K_i$.

460 Hence, in this case we set the neighbourhood function $N$ to $\mathcal{N}$ because (i) $\text{SND}_{uu'}^{K_i, \mathcal{N}}$ is as strong as $\text{SND}_{uu'}^{K_i, \text{all}}$, and (ii) the DC of $\text{SND}_{uu'}^{K_i, \mathcal{N}}$ may be ensured while ensuring the DC of $\text{LAD}_{uu'}$. Indeed, as both constraints consider the same neighbourhood function $\mathcal{N}$, we simply modify the definition of the set of edges $E_{uu'}$ of the bipartite graph $G_{uu'}$ used to ensure the DC of LAD as follows:

$$E_{uu'} = \{(v, v') \in \mathcal{N}(u) \times \mathcal{N}(u') | v' \in D_v \wedge K_i(u, v) \le K_i(u', v'), \forall i \in [3, C]\}$$

465 In other words, we simply have to remove from the bipartite graph every edge which is incompatible with respect to $K_i$, and then we can ensure the DC of both $\text{LAD}_{uu'}$ and $\text{SND}_{uu'}^{K_i, \mathcal{N}}$, $\forall i \in [3, C]$ by searching for a covering matching in the resulting bipartite graph.

468 Neighbourhood function considered when $\{a, b\}$ is not an edge of $m$. Let us now consider the case where $a$ and $b$ are not connected by an edge, which is the case when $m = P_2$. In this case, it may be interesting to consider the neighbourhood function all. For example, let us consider the vertices 1 and 7' in the graphs $G$ and $G'$ displayed in Fig. 3:

472 • if $N = \mathcal{N}$, then $N(1) = \{2, 4\}$ and $N(7') = \{1', 2'\}$ and, as $P_2(1, 2) = P_2(1, 4) = P_2(7', 1') = P_2(7', 2') = 1$, the constraint $\text{SND}_{17'}^{P_2, \mathcal{N}}$ is DC and 7' is not removed from $D_1$;

474 • if $N = \text{all}$, then $N(1) = \{2, 3, 4, 5\}$ and $N(7') = \{1', 2', 3', 4', 5', 6'\}$ and, as $P_2(1, 5) = 2$ and $P_2(7', u') = 1, \forall u' \in N(7')$, we must remove 7' from $D_1$ to ensure the DC of $\text{SND}_{17'}^{P_2, \text{all}}$.

476 If $\text{SND}_{uu'}^{P_2, \text{all}}$ may allow us to filter more values than $\text{SND}_{uu'}^{P_2, \mathcal{N}}$, it is also more expensive as we need to search for covering matchings of new bipartite graphs that may be very large. A tradeoff has to be found between the cost of the propagation and its strength. Preliminary experiments have shown us that for some instances with dense graphs (such as those of *randER* and *phase*, for example) ensuring the DC of $\text{SND}_{uu'}^{P_2, \text{all}}$ constraints is unreasonably long.

481 Hence, we ensure a weaker consistency than DC for $\text{SND}_{uu'}^{P_2, \text{all}}$. More precisely, we ensure that for each vertex $v \in \text{all}(u)$ there exists a different vertex $v' \in \text{all}(u')$ such that $P_2(u, v) \le P_2(u', v')$. In other words, we ensure that there exists an injective mapping $m: \text{all}(u) \to \text{all}(u')$ such that $\forall v \in \text{all}(u), P_2(u, v) \le P_2(u', m(v))$. This is weaker than DC because we do not enforce $m(v)$ to belong to $D_v$.

486 Efficient propagation of $\text{SND}_{uu'}^{P_2, \text{all}}$. Let us first describe a naive procedure for deciding whether there exists an injective mapping $m: \text{all}(u) \to \text{all}(u')$ such that $\forall v \in \text{all}(u), P_2(u, v) \le P_2(u', m(v))$. First, we define two arrays $A_u$ and $A'_{u'}$, indexed by the vertices in $\text{all}(u)$ and $\text{all}(u')$, respectively, and such that $\forall v \in \text{all}(u), A_u[v] = P_2(u, v)$, and $\forall v' \in \text{all}(u'), A'_{u'}[v'] = P_2(u', v')$. Then, we sort $A_u$ and $A'_{u'}$ by decreasing values and check that $A_u[i] \le A'_{u'}[i]$ for each index $i$ of
A. For example, let us consider the SIP instance displayed in Fig. 3, and let us consider the case where $u = 2$ and $u' = 3'$. We build the arrays:

$$A_2 = \begin{bmatrix} 1 & 1 & 2 & 2 \\ 1 & 3 & 4 & 5 \end{bmatrix} \quad A'_{3'} = \begin{bmatrix} 1 & 2 & 1 & 1 & 1 & 1 \\ 1' & 2' & 4' & 5' & 6' & 7' \end{bmatrix}$$

In other words, there is one 2-path from 2 to each vertex in $\{1,3\}$, and there are two 2-paths from 2 to each vertex in $\{4,5\}$. Similarly, there is one 2-path from $3'$ to each vertex in $\{1',4',5',6',7'\}$, and there are two 2-paths from $3'$ to $2'$. Once sorted, these arrays become:

$$A_2 = \begin{bmatrix} 2 & 2 & 1 & 1 \end{bmatrix} \quad A'_{3'} = \begin{bmatrix} 2 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}$$

As the second element of $A_2$ is greater than the second element of $A'_{3'}$, we conclude that $3'$ must be removed from $D_2$.

However, sorting $A_u$ and $A'_{u'}$ is expensive as these arrays have $n-1$ and $n'-1$ elements, respectively. To lower the time complexity, we take inspiration from the first step of the counting sort [CLRS09], as displayed in Algorithm 2.

<caption>**Algorithm 2: Propagation of SND$_{uu'}^{P_2,all}$**</caption>

**Input:** A pattern vertex $u \in V$ and a target vertex $u' \in D_u$
1 &nbsp;&nbsp;**for** each $i \in [1, d]$ **do** $B_u[i] \leftarrow \#\{v \in \text{all}(u) \mid P_2(u, v) = i\}$;
2 &nbsp;&nbsp;**for** each $i \in [1, d-1]$ **do** $B'_{u'}[i] \leftarrow \#\{v' \in \text{all}(u') \mid P_2(u', v') = i\}$;
3 &nbsp;&nbsp;**for** each $i \in [1, d]$ **do** $B'_{u'}[i] \leftarrow \#\{v' \in \text{all}(u') \mid P_2(u', v') \ge d\}$;
4 &nbsp;&nbsp;**if** $\exists i \in [1, d]$, $\sum_{j=i}^d B_u[j] > \sum_{j=i}^d B'_{u'}[j]$ **then** remove $u'$ from $D_u$;

Indeed, the number of 2-paths between $u$ and another vertex $v \in \text{all}(u)$ cannot exceed $d$ and therefore all values in $A_u$ belong to $[0, d]$. Similarly, all values in $A'_{u'}$ belong to $[0, d']$ but, when a value in $A'_{u'}$ is greater than $d$, we can replace it with $d$ without changing the final result. Hence, we can define two arrays $B_u$ and $B'_{u'}$ indexed by $[1, d]$; for each $i \in [1, d]$, $B_u[i]$ gives the number of vertices $v$ such that $A_u[v] = i$ (Line 1) and $B'_{u'}[i]$ gives the number of vertices $v'$ such that $A'_{u'}[v'] = i$ (Line 2), except when $i = d$ in which case $B'_{u'}[d]$ gives the number of vertices $v'$ such that $A'_{u'}[v'] \ge d$ (Line 3). Given $B_u$ and $B'_{u'}$, we can build the sorted arrays $A_u$ and $A'_{u'}$ in linear time (except that values larger than $d$ in $A'_{u'}$ are replaced with $d$), just like in the counting sort: the $B_u[0]$ last elements of the sorted array are equal to 0, the $B_u[1]$ previous elements are equal to 1, etc. However, it is not necessary to explicitly build the sorted arrays $A$ and $A'$, and we can simply compare cumulated sums of elements of $B$ and $B'$, starting from the end of the two arrays. More precisely, if for some $i \in [1, d]$ we have $\sum_{j=i}^d B_u[i] > \sum_{j=i}^d B'_{u'}[j]$, we can remove $u'$ from $D_u$ (line 4). On our example, we build the arrays:

$$B_2 = \begin{bmatrix} 0 & 2 & 2 & 0 & 0 \\ 0 & 1 & 2 & 3 & 4 \end{bmatrix} \quad B'_{3'} = \begin{bmatrix} 0 & 5 & 1 & 0 & 0 \\ 0 & 1 & 2 & 3 & 4 \end{bmatrix}$$

In other words, $A_2$ has no element equal to 0, 3, or 4, whereas it has two elements equal to 1 and two elements equal to 2. Similarly, $A'_{3'}$ has no element equal to 0, 3, or 4, whereas it has five elements equal to 1 and one element equal to 2. Finally, as $\sum_{j=2}^4 B_2[j] = 2$ whereas $\sum_{j=2}^4 B'_{3'}[j] = 1 < 2$, we can remove $3'$ from $D_2$.

The computation of arrays $B_u$ and $B'_{u'}$ is done incrementally when enumerating 2-paths to compute $P_2$ labels: when $P_2(u, v)$ is set to 1 for some $(u, v) \in V \times V$, we increment $B_u[1]$; then each time $P_2(u, v)$ is incremented from $k$ to $k+1$ with $k > 1$, we decrement $B_u[k]$ and increment $B_u[k+1]$. In this case, the time complexity of Algorithm 2 is $O(d)$ (provided that sums of values are incrementally computed at Line 4).
<table><caption>Table 1: Impact of supplemental labels. The first column specifies the considered set $L_i$ of supplemental labels. For each set $L_i$ and each benchmark, $\#S_{pre}$ is the number of instances solved before starting search, $t_{pre}$ is the average time needed to compute all labels and filter domains before starting search, $rm_V$ is the percentage of target vertices removed from domains before starting search, and $rm_E$ is the percentage of edges removed from bipartite graphs $G_{uu'}$ ($rm_V$ and $rm_E$ are computed on average on instances that are not solved before starting the search).</caption><thead><tr><th>Supp. labels</th><th></th><th>#S<sub>pre</sub></th><th>t<sub>pre</sub></th><th>rm<sub>V</sub></th><th>rm<sub>E</sub></th><th></th><th>#S<sub>pre</sub></th><th>t<sub>pre</sub></th><th>rm<sub>V</sub></th><th>rm<sub>E</sub></th></tr></thead><tbody><tr><td>L<sub>0</sub> = Ø</td><td rowspan="6">LV</td><td>2179</td><td>0.007</td><td>8.45</td><td>-</td><td rowspan="6">randER</td><td>0</td><td>0.062</td><td>4.14</td><td>-</td></tr><tr><td>L<sub>1</sub> = {ω}</td><td>2623</td><td>0.012</td><td>12.36</td><td>-</td><td>0</td><td>0.063</td><td>6.78</td><td>-</td></tr><tr><td>L<sub>2</sub> = L<sub>1</sub> ∪ {P<sub>2</sub>}</td><td>2743</td><td>0.015</td><td>14.56</td><td>-</td><td>0</td><td>0.069</td><td>13.12</td><td>-</td></tr><tr><td>L<sub>3</sub> = L<sub>2</sub> ∪ {K<sub>3</sub>}</td><td>2762</td><td>0.012</td><td>15.25</td><td>1.30</td><td>0</td><td>0.143</td><td>21.76</td><td>4.36</td></tr><tr><td>L<sub>4</sub> = L<sub>3</sub> ∪ {K<sub>4</sub>}</td><td>2764</td><td>0.013</td><td>15.28</td><td>1.34</td><td>0</td><td>0.188</td><td>28.47</td><td>4.50</td></tr><tr><td>L<sub>5</sub> = L<sub>4</sub> ∪ {K<sub>5</sub>}</td><td>2765</td><td>0.014</td><td>15.31</td><td>1.37</td><td>0</td><td>0.212</td><td>31.82</td><td>4.33</td></tr><tr><td>L<sub>6</sub> = L<sub>5</sub> ∪ {K<sub>6</sub>}</td><td></td><td>2766</td><td>0.019</td><td>15.32</td><td>1.38</td><td></td><td>0</td><td>0.232</td><td>32.40</td><td>4.29</td></tr><tr><td>L<sub>0</sub> = Ø</td><td rowspan="6">rand</td><td>20</td><td>0.010</td><td>65.13</td><td>-</td><td rowspan="6">phase</td><td>0</td><td>0.003</td><td>0.00</td><td>-</td></tr><tr><td>L<sub>1</sub> = {ω}</td><td>20</td><td>0.010</td><td>65.21</td><td>-</td><td>0</td><td>0.003</td><td>0.00</td><td>-</td></tr><tr><td>L<sub>2</sub> = L<sub>1</sub> ∪ {P<sub>2</sub>}</td><td>20</td><td>0.010</td><td>66.53</td><td>-</td><td>0</td><td>0.006</td><td>0.00</td><td>-</td></tr><tr><td>L<sub>3</sub> = L<sub>2</sub> ∪ {K<sub>3</sub>}</td><td>20</td><td>0.015</td><td>67.44</td><td>6.51</td><td>0</td><td>0.009</td><td>0.00</td><td>0.00</td></tr><tr><td>L<sub>4</sub> = L<sub>3</sub> ∪ {K<sub>4</sub>}</td><td>20</td><td>0.015</td><td>67.44</td><td>6.51</td><td>0</td><td>0.048</td><td>0.00</td><td>0.00</td></tr><tr><td>L<sub>5</sub> = L<sub>4</sub> ∪ {K<sub>5</sub>}</td><td>20</td><td>0.016</td><td>67.44</td><td>6.51</td><td>0</td><td>0.441</td><td>0.00</td><td>0.00</td></tr><tr><td>L<sub>6</sub> = L<sub>5</sub> ∪ {K<sub>6</sub>}</td><td></td><td>20</td><td>0.017</td><td>67.44</td><td>6.51</td><td></td><td>0</td><td>2.961</td><td>0.00</td><td>0.00</td></tr><tr><td>L<sub>0</sub> = Ø</td><td rowspan="4">meshes</td><td>2884</td><td>0.007</td><td>27.73</td><td>-</td><td rowspan="4">images</td><td>5795</td><td>0.029</td><td>76.81</td><td>-</td></tr><tr><td>L<sub>1</sub> = {ω}</td><td>2884</td><td>0.007</td><td>27.73</td><td>-</td><td>5795</td><td>0.029</td><td>82.54</td><td>-</td></tr><tr><td>L<sub>2</sub> = L<sub>1</sub> ∪ {P<sub>2</sub>}</td><td>2885</td><td>0.010</td><td>27.80</td><td>-</td><td>6216</td><td>0.016</td><td>95.54</td><td>-</td></tr><tr><td>L<sub>3</sub> = L<sub>2</sub> ∪ {K<sub>3</sub>}</td><td>2893</td><td>0.010</td><td>27.97</td><td>0.00</td><td>6249</td><td>0.015</td><td>93.12</td><td>0.27</td></tr></tbody></table>

Algorithm 2 is called only once for each pattern vertex $u \in V$ and each target vertex $u' \in D_u$ before starting the search. It is useless to call it during the search as domains are not considered. During the search, we propagate Constraint 4 with a simple forward-checking: whenever $D_u$ is reduced to a singleton $\{u'\}$, we remove from $D_v$ any value $v'$ such that $P_2(u, v') > P_2(u', v')$.

An alternative is to ensure the DC of Constraint 4, using one of the numerous algorithms that have been designed for propagating table constraints [YXW20]: this may filter more values than forward-checking, but this is also more expensive. Hence, we choose not to ensure DC on this constraint.

## 5.4 Experimental results

Before starting search, we compute all supplemental labels, filter domains with respect to unary labels, ensure the DC of $LAD_{uu'}$ and $SND_{uu'}^{K_i, N}$ constraints (as described in the previous section), the DC of the GAD constraint, and propagate $SND_{uu'}^{P_2, all}$ constraints (as described in the previous section). In some cases, this allows us to detect an inconsistency, in which case the instance is solved at the root of the search tree.

In Table 1, we report the number of instances solved during this preprocessing step, for each benchmark, and for different sets of supplemental labels. When this set is empty ($L_0$), we only ensure the DC of LAD and GAD, and this allows us to solve 2179 LV instances, 20 rand instances, 2884 meshes instances, and 5795 images instances. When considering all labels ($L_6$), the number of instances solved during the preprocessing step is increased by 587 for LV, 9 for meshes, and 454 for images. Instances of meshes and images never contain cliques of order larger than 3, and that is why we do not report results with $K_i$ when $i > 3$ as they are the same as with $K_3$. All instances
which are solved during the preprocessing step are infeasible (inconsistency is detected when a domain is wiped out). Note that *randER* only contains feasible instances whereas *rand* contains only 20 infeasible instances. This explains why no instance of *randER* and only 20 instances of *rand* are solved during the preprocessing step. However, when an instance is not solved, variable domains may be tightened. For example, for *randER*, only 4.14% of the values are removed without using any label; $\omega$ labels slightly increase this percentage whereas $P_2$ and $K_i$ allow us to filter much more values, up to 32.4% when considering all labels in $L_6$. In some cases, this percentage may decrease when considering more labels. This is the case for *images*: it decreases from 95.54% with $L_2$ to 93.12% with $L_3$. This comes from the fact that these percentages are computed for unsolved instances only, and there are fewer unsolved instances with $L_3$ than with $L_2$. $K_i$ labels are used to remove some edges from the bipartite graphs used to ensure the DC of SND and LAD constraints, as explained in the previous section. The percentage of removed edges is rather low on average for all instances, because it is close to 0 for many instances. However, on some instances, this percentage may be much larger (e.g., up to 41% for LV) and, in this case, the time needed to compute covering matchings is strongly decreased. The benchmark *phase* is the only benchmark for which domains are not tightened and no edge is removed from bipartite graphs. For this benchmark, it is useless to compute labels.

The average preprocessing time when considering no label ($L_0$) ranges from 3ms for *phase* to 62ms for *randER*. When considering more labels, this time usually increases, but it may decrease in some cases because inconsistencies are detected earlier when considering more labels, and we stop the computation of labels as soon as an inconsistency has been detected. For example, some graphs of *LV* contain large cliques (e.g., g21, g40, and g50 contain cliques of order 32, 36, and 42, respectively), and computing $\omega$ labels for these graphs is expensive. As $K_3$ labels are quickly computed, we first compute $K_3$ labels and if an inconsistency is detected thanks to them we do not compute $\omega$ labels. For all benchmarks but *randER* and *phase*, the average time of the preprocessing step never exceeds 19ms for the largest set of labels $L_6$, and this time is usually negligible compared to the solving time. However, for *randER* and *phase*, which contain dense graphs, the preprocessing step needs more time (especially for *phase* instances for which it is close to 3s). For *randER*, it is usually worth spending time for computing labels as this allows us to tighten domains and remove some edges from bipartite graphs. However, for *phase*, labels are useless. For this benchmark, the parameter $C$ should be set to 0 in order to avoid computing $K_i$ labels. However, the default value for this parameter is $C = 6$ as this allows us to solve more instances and remove more values and edges without significantly increasing CPU time for most benchmarks.

PathLAD and its new implementation V0 only consider $P_2$ labels: unary $P_2$ labels are used to remove from $D_u$ every target vertex $u'$ such that $P_2(u) > P_2(u')$, and binary $P_2$ labels are used to filter edges of every bipartite graph $G_{uu'}$ (every edge $(v, v')$ such that $P_2(u, v) > P_2(u', v')$ is removed from $E_{uu'}$). Let V1 denote the solver obtained when considering labels $\omega, P_2, K_3, K_4, K_5$, and $K_6$. We compare CPU times of V0 and V1 on a per instance basis in Fig. 4. It shows us that, if for some instances V1 needs more time than V0 (points in the upper left triangle), for a majority of instances V1 is more efficient (points in the lower right triangle). In particular, V1 is able to solve 37 more instances than V0 when considering a time limit of 100s (mostly *LV* and *randER* instances). The aggregate speed-up between V0 and V1 for a time limit of 100s is $ASU_{100}(\frac{V0}{V1}) = 2.66$.

## 6 Variable Ordering Heuristics (V2)

*Variable Ordering Heuristics* (VrOHs) are used at each node of the search tree to choose the next variable to assign. In our context, variables are associated with pattern vertices so that the goal is to choose a pattern vertex which is not yet mapped to a target vertex. In PathLAD (and in V0 and V1 too), we use the VrOH called *minDom*, i.e., we select the variable $x_u$ such that $\#D_u$ is minimal. In case of ties, we use the same rule as the one used in Glasgow to break ties, i.e., we
<figure><img src="image_5.png" /><figcaption>Figure 4: Scatter plot to compare V0, that considers $P_2$ labels (x axis) with V1, that considers labels $\omega, P_2, K_3, K_4, K_5$, and $K_6$ (y axis), for a time limit of 100 seconds.</figcaption></figure>

select the variable $x_u$ such that $d^\circ(u)$ is maximal.

Many other generic VrOHs have been proposed. A well-known VrOH is the weighted degree heuristic defined in [BHLS04], denoted $wdeg$. More precisely, a weight $w(c)$ is associated with every constraint $c$. This weight is initialised to 1 when starting the search, and it is incremented each time the propagation of $c$ wipes out a domain. $wdeg$ selects the variable that maximises the sum of the weights over all constraints involving this variable and at least another non assigned variable. $wdeg$ is usually combined with minDom by selecting the variable that minimises the ratio between the size of the domain and the weighted sum. This VrOH is denoted $\frac{dom}{wdeg}$.

In [WLPT19], variants of $\frac{dom}{wdeg}$ are introduced, and it is shown that the variant denoted ca.cd obtains the best results. The idea is to associate weights with constraint/variable couples, and to increment the weight associated with a couple $(c, x)$ whenever the propagation of $c$ wipes out a domain, and $x$ is an unassigned variable involved in $c$. Also, instead of incrementing weights of 1, they are incremented of $\frac{1}{\#fut(c) \times (1+\#D_x)}$ where $fut(c)$ is the set of variables of $c$ that are not yet assigned and $D_x$ is the domain of $x$.

In [LYL21], four other VrOHs are introduced, i.e., FRB, FLB, FRBA, and FLBA, and we have compared them on our SIP benchmark. The best performing one is FRBA, and this is consistent with experiments reported in [LYL21]. FRBA selects the variable $x_u$ that maximises

$$ \frac{failNum(x_u)/assignNum(x_u) + 1/(totalFailure - lastFailure(x_u) + 1)}{\#D_u} $$

where $failNum(x_u)$ is the total number of failures caused by $x_u$, $assignNum(x_u)$ is the total number of times $x_u$ has been assigned since the beginning of search, $totalFailure$ is the total number of failures detected since the beginning of search, and $lastFailure(x_u)$ is the totalFailure value of the last failure caused by $x_u$.

In Table 2, we report experimental results obtained with minDom (corresponding to V1), ca.cd, FRBA, and $\frac{dom}{wdeg}$. Both ca.cd, FRBA, and $\frac{dom}{wdeg}$ perform much better than minDom. The two benchmarks that are the most positively impacted by these weight-based VrOHs are LV and meshes, for which the number of unsolved instances drops down from 142 to 109, 108, or 103 for LV and from 45 to 5, 3, or 1 for meshes. For randER and phase, the number of solved instances is slightly improved for $\frac{dom}{wdeg}$. For images and rand that contain only very easy instances, VrOHs have no impact at all. When considering all benchmarks, the best performing VrOH is $\frac{dom}{wdeg}$, and we call
<table><caption>Table 2: Comparison of VrOHs for a time limit of 100 seconds. #U is the number of unsolved instances and t is the average solving time of the solved instances.</caption><thead><tr><th rowspan="2"></th><th colspan="2">LV</th><th colspan="2">images</th><th colspan="2">meshes</th><th colspan="2">rand</th><th colspan="2">randER</th><th colspan="2">phase</th><th colspan="2">All</th></tr><tr><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th></tr></thead><tbody><tr><td>minDom</td><td>131</td><td>0.35</td><td>0</td><td>0.02</td><td>38</td><td>0.19</td><td>0</td><td>0.03</td><td>32</td><td>6.43</td><td>166</td><td>17.55</td><td>367</td><td>0.29</td></tr><tr><td>ca.cd</td><td>109</td><td>0.32</td><td>0</td><td>0.02</td><td>5</td><td>0.15</td><td>0</td><td>0.03</td><td>29</td><td>6.51</td><td>170</td><td>14.62</td><td>313</td><td>0.28</td></tr><tr><td>FRBA</td><td>108</td><td>0.37</td><td>0</td><td>0.02</td><td>3</td><td>0.16</td><td>0</td><td>0.03</td><td>35</td><td>6.34</td><td>164</td><td>27.58</td><td>310</td><td>0.31</td></tr><tr><td>dom/wdeg</td><td>103</td><td>0.34</td><td>0</td><td>0.02</td><td>1</td><td>0.11</td><td>0</td><td>0.03</td><td>30</td><td>6.86</td><td>163</td><td>21.08</td><td>297</td><td>0.29</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: Scatter plot to compare V1 (x-axis), that uses the *minDom* VrOH, with V2 (y-axis), that uses the *dom/wdeg* VrOH, for a time limit of 100 seconds.</figcaption></figure>

V2 the variant obtained from V1 by replacing *minDom* with *dom/wdeg*.

In Fig. 5, we compare the CPU times of V1 and V2 on a per instance basis. V2 solves 70 more instances than V1 within a time limit of 100s, and the aggregate speed-up between V1 and V2 for this time limit is $ASU_{100}(\frac{V1}{V2}) = 5.88$.

# 7 Restarts and value ordering heuristics (V3)

Combining randomised search with restarts is a classical approach for avoiding heavy tails when exploring search spaces [GSK98]. To avoid exploring several times the same part of the search space, random restarts are usually combined with nogood recording: before restarting the search, nogood constraints are added to forbid the exploration of previously seen nodes. In [LSTV07], Lecoutre et al introduce negative last decision nogoods that are extracted at each restart of a backtrack search algorithm. In [LSZ16], it is shown that these nogoods are increasing: it means that it suffices to record the sequence of decisions in the rightmost branch exactly as it is instead of extracting a set of nogoods. Authors of [LSZ16] also introduce a global constraint to ensure the DC of increasing nogoods in linear time with respect to the sum of all domains, i.e., in $O(nn')$ when considering the SIP.

Increasing nogoods are defined for binary search, i.e., when each node of the search tree has at most two children: one with a positive decision ($x = v$ where $x$ is a variable and $v$ is a value in the domain of $x$), and one with a negative decision ($x \neq v$). PathLAD and LAD2025 perform n-ary searches such that each node of a search tree is associated with a variable and has one child
<table><caption>Table 3: Impact of restarts with nogood recording and of VIOHs for a time limit of 100 seconds. V2 is the variant without random restarts nor VIOH. RR uses random restarts and randomly chooses values. RRbia uses random restarts and the degree-biased VIOH. RRsav uses random restarts and the degree-biased VIOH with saving. #U is the number of unsolved instances and t is the average solving time of the solved instances.</caption><thead><tr><th rowspan="2"></th><th colspan="2">LV</th><th colspan="2">images</th><th colspan="2">meshes</th><th colspan="2">rand</th><th colspan="2">randER</th><th colspan="2">phase</th><th colspan="2">All</th></tr><tr><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th></tr></thead><tbody><tr><td>V2</td><td>103</td><td>0.34</td><td>0</td><td>0.01</td><td>1</td><td>0.11</td><td>0</td><td>0.03</td><td>30</td><td>6.86</td><td>163</td><td>21.08</td><td>297</td><td>0.29</td></tr><tr><td>RR</td><td>92</td><td>0.42</td><td>0</td><td>0.02</td><td>2</td><td>0.04</td><td>0</td><td>0.03</td><td>32</td><td>6.01</td><td>151</td><td>18.70</td><td>277</td><td>0.29</td></tr><tr><td>RRbia</td><td>101</td><td>0.43</td><td>0</td><td>0.02</td><td>1</td><td>0.05</td><td>0</td><td>0.03</td><td>29</td><td>5.41</td><td>119</td><td>14.49</td><td>250</td><td>0.30</td></tr><tr><td>RRsav</td><td>90</td><td>0.47</td><td>0</td><td>0.01</td><td>2</td><td>0.05</td><td>0</td><td>0.03</td><td>28</td><td>5.74</td><td>109</td><td>12.40</td><td>229</td><td>0.32</td></tr></tbody></table>

for each value in the domain of this variable. We have made experiments with binary search, and noticed that n-ary search is often more efficient than binary search when solving SIPs. Hence, we have adapted increasing nogoods to n-ary search: when search must be restarted, we build the sequence of positive decisions associated with the rightmost branch of the n-ary tree together with the negation of each decision associated with a child of any node of this rightmost branch, and we use this sequence of positive and negative decisions to define an increasing nogood as proposed in [LSZ16].

Search is restarted when the number of failures reaches a given cutoff limit, and we consider a geometric progression for this cutoff limit: it is initialised to 10, and it is multiplied by 1.1 after each restart.

Restarts must be combined with randomised search, so that decisions after a restart are expected to be different from decisions before the restart. Randomisation is usually achieved through VrOHs. In particular, weight-based VrOHs ensure that the selected variables are modified after each failure as weights are updated. In [ADH$^+$19], randomisation is achieved through *Value Ordering Heuristics* (VIOHs), and a new randomised VIOH is introduced: given a pattern vertex $u$, they branch on the target vertex $u' \in D_u$ with probability

$$ \frac{2^{d^{\circ}(u')}}{\sum_{v' \in D_u} 2^{d^{\circ}(v')}}, $$

In other words, if $D_u$ contains two target vertices $u'$ and $v'$ such that $d^{\circ}(u') = d^{\circ}(v') + 1$, then $u'$ is twice more likely to be selected than $v'$. This VIOH is said to be degree-biased.

We propose a new VIOH which combines this degree-biased VIOH with saving. Phase saving has been introduced for SAT solvers in [PD07]. The idea is to reuse the last value of a Boolean variable (its phase) when it is reconsidered. This VIOH has been extended to constrained optimisation problems in [VP17, DCS18]: the idea is to give priority to the value the variable is assigned to in the last solution found, thus concentrating the search around the current best solution, just as in Large Neighbourhood Search (LNS). For satisfaction problems, this is meaningless as the search terminates once a solution is found. Hence, instead of giving priority to values in the last solution, we propose to give priority to values in the largest consistent assignment built since the beginning of the search. More precisely, each time a domain is wiped out, we consider the assignment $A$ before the decision that triggered the failure, and if the number of assigned variables in $A$ is larger than the largest memorised assignment, then we memorise $A$. When choosing a value for a pattern vertex $u$, if $x_u$ is assigned to $u'$ in the last memorised assignment and if $u' \in D_u$, then the first value selected for $x_u$ is $u'$ and the next values are selected according to the degree-biased VIOH; otherwise all values for $x_u$ are selected according to the degree-biased VIOH. This new VIOH is called degree-biased VIOH with saving.
<figure><img src="image_7.png" /><figcaption>Figure 6: Scatter plot to compare V2 (x-axis), that does not perform random restarts and has no V1OH, with V3 (y-axis), that performs random restarts with nogood recording and uses the degree-biased V1OH with saving.</figcaption></figure>

In Table 3, we report experimental results obtained without random restarts and when no V1OH is used (in which case, values are selected in the order defined by the domain), with results obtained when using random restarts with nogood recording and with different V1OHs: random choice, degree-biased V1OH, and degree-biased V1OH with saving. Random restarts improve results for *LV* and *phase*, but they slightly degrade results for *meshes* and *randER*. The degree-biased V1OH strongly improves results for *phase*, and the degree-biased V1OH with saving further improves these results. For *LV*, the degree-biased V1OH degrades results, while saving improves them.

This new variant of LAD with random restarts, nogood recording, and degree-biased V1OH with saving is denoted V3. In Fig. 6, we compare the CPU times of V2 and V3 on a per instance basis. V3 solves 68 more instances than V2 within a time limit of 100s, and the aggregate speed-up between V2 and V3 for this time limit is $ASU_{100}(\frac{V^2}{V_3}) = 4.97$.

# 8 Per instance filtering selection (LAD2025)

LAD filtering is expensive on dense graphs as each bipartite graph $G_{uv'}$ has $O(dd')$ edges. In [WJC24], the filtering level is dynamically adapted during the search. More precisely, Constraints 1 and 2 are first propagated with a simple forward-checking, i.e., when the domain of a variable $x_u$ is reduced to a singleton $\{u'\}$, $u'$ is removed from the domains of all other variables and, for each pattern edge $\{u, v\}$, every target vertex $v' \in D_v$ which is not adjacent to $u'$ is removed from $D_v$. Then, depending on whether the target graph is sparse or not, and also on whether *nbNodes* and *nbFail/nbNodes* exceed some given thresholds or not, the DC of LAD and/or GAD constraints may be ensured (where *nbNodes* is the number of explored nodes, *nbFail* is the number of failures, and the target graph is defined to be sparse when its median vertex degree is smaller than a parameter *deg*<sub>m</sub> which is set to 20 in [WJC24]).

In [KMS16], another possibility is explored to take advantage of solver complementarity: a per instance selection model is learned on a training set of instances and, each time a new instance must be solved, the learned model is used to select the expected best solver for this instance.

In this section, we consider a very simple rule to select the filtering level: when the density of both *G* and *G'* exceeds 0.15, we propagate Constraint (2) with a simple forward-checking instead of ensuring the DC of LAD constraints (in all cases, we ensure the DC of the GAD constraint).
<table><caption>Table 4: Number of unsolved instances (#U) and average solving time of the solved instances in seconds (t) when the time limit is 100 seconds for the different variants from PathLAD to LAD2025.</caption><thead><tr><th rowspan="2"></th><th colspan="2">LV</th><th colspan="2">images</th><th colspan="2">meshes</th><th colspan="2">rand</th><th colspan="2">randER</th><th colspan="2">phase</th><th colspan="2">All</th></tr><tr><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th><th>#U</th><th>t</th></tr></thead><tbody><tr><td>PathLAD</td><td>203</td><td>1.01</td><td>0</td><td>0.86</td><td>45</td><td>0.26</td><td>0</td><td>0.26</td><td>97</td><td>8.74</td><td>181</td><td>27.13</td><td>526</td><td>0.86</td></tr><tr><td>V0</td><td>151</td><td>0.56</td><td>0</td><td>0.03</td><td>38</td><td>0.20</td><td>0</td><td>0.02</td><td>49</td><td>9.04</td><td>166</td><td>13.76</td><td>404</td><td>0.37</td></tr><tr><td>V1</td><td>131</td><td>0.35</td><td>0</td><td>0.02</td><td>38</td><td>0.19</td><td>0</td><td>0.03</td><td>32</td><td>6.43</td><td>166</td><td>17.55</td><td>367</td><td>0.29</td></tr><tr><td>V2</td><td>103</td><td>0.34</td><td>0</td><td>0.02</td><td>1</td><td>0.11</td><td>0</td><td>0.03</td><td>30</td><td>6.86</td><td>163</td><td>21.08</td><td>297</td><td>0.29</td></tr><tr><td>V3</td><td>90</td><td>0.47</td><td>0</td><td>0.02</td><td>2</td><td>0.05</td><td>0</td><td>0.03</td><td>28</td><td>5.74</td><td>109</td><td>12.40</td><td>229</td><td>0.32</td></tr><tr><td>LAD2025</td><td>90</td><td>0.48</td><td>0</td><td>0.02</td><td>2</td><td>0.05</td><td>0</td><td>0.03</td><td>17</td><td>4.33</td><td>84</td><td>8.25</td><td>193</td><td>0.29</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 7: Scatter plot to compare V3 (x-axis), that does not use the selection rule, with LAD2025 (y-axis), that uses the selection rule.</figcaption></figure>

This choice is done at the beginning of the search and it is not changed during the search, in the spirit of a per-instance selector though we have not used machine learning to learn the selection rule. Our goal is to show that the solution process is improved when replacing the LAD filtering with a simple forward-checking on dense graphs. We believe better results could be obtained by using an automatic per-instance selection tools such as Llama [Kot13], and by considering more variants (e.g., we could also automatically set the parameter C used to upper bound the order of enumerated cliques when computing supplemental labels). However, this kind of study goes beyond the scope of this paper.

We denote LAD2025 the new variant obtained from V3 by using our filtering level selection rule. In Table 4, we report results obtained with all variants introduced in this paper, from PathLAD to LAD2025. The filtering level selection rule does not change performance for images, meshes, and rand as these benchmarks only contain sparse graphs. For LV, it does not change the number of solved instances. For randER and phase, it allows us to solve 11 and 25 more instances as these benchmarks contain many instances with very dense graphs on which forward-checking is more efficient than LAD filtering.

In Fig. 7, we compare the CPU times of V3 and LAD2025 on a per instance basis. The aggregate speed-up between V3 and LAD2025 for a limit of 100s is $ASU_{100}(\frac{V3}{LAD2025}) = 3.07$.

The number of unsolved instances (for a time limit of 100 seconds) has decreased from 526 with PathLAD to 193 with LAD2025. On easy benchmarks (images and rand), the improvement is mostly coming from refactoring. For LV, which is composed of instances with very different
<figure><img src="image_9.png" /><figcaption>Figure 8: Evolution of the cumulative number of instances solved by each variant with respect to time (with a logscale). The right figure is a zoom where $x \ge 0.05s$ and $y \ge 13400$ instances.</figcaption></figure>

<table><caption>Table 5: Number of unsolved instances within a time limit of 100s (#U<sub>100</sub>) and 2000s (#U<sub>2000</sub>).</caption><thead><tr><th rowspan="2"></th><th colspan="2">LV</th><th colspan="2">meshes</th><th colspan="2">rand</th><th colspan="2">randER</th><th colspan="2">phase</th><th colspan="2">All</th></tr><tr><th>#U<sub>100</sub></th><th>#U<sub>2000</sub></th><th>#U<sub>100</sub></th><th>#U<sub>100</sub></th><th>#U<sub>100</sub></th><th>#U<sub>2000</sub></th><th>#U<sub>100</sub></th><th>#U<sub>2000</sub></th><th>#U<sub>100</sub></th><th>#U<sub>2000</sub></th><th>#U<sub>100</sub></th><th>#U<sub>2000</sub></th></tr></thead><tbody><tr><td>RI</td><td>485</td><td>348</td><td>0</td><td>330</td><td>323</td><td>24</td><td>18</td><td>32</td><td>4</td><td>186</td><td>169</td><td>1057</td><td>862</td></tr><tr><td>VF3</td><td>439</td><td>328</td><td>0</td><td>329</td><td>323</td><td>27</td><td>19</td><td>31</td><td>6</td><td>184</td><td>161</td><td>1010</td><td>837</td></tr><tr><td>PathLAD</td><td>203</td><td>137</td><td>0</td><td>45</td><td>30</td><td>0</td><td>0</td><td>97</td><td>55</td><td>181</td><td>162</td><td>526</td><td>384</td></tr><tr><td>Glasgow</td><td>165</td><td>141</td><td>0</td><td>36</td><td>29</td><td>0</td><td>0</td><td>7</td><td>0</td><td>96</td><td>74</td><td>304</td><td>244</td></tr><tr><td>PathLAD+23</td><td>159</td><td>106</td><td>0</td><td>18</td><td>10</td><td>0</td><td>0</td><td>90</td><td>2</td><td>154</td><td>88</td><td>421</td><td>206</td></tr><tr><td>PathLAD+24</td><td>140</td><td>100</td><td>0</td><td>14</td><td>6</td><td>0</td><td>0</td><td>84</td><td>0</td><td>119</td><td>65</td><td>357</td><td>171</td></tr><tr><td>LAD2025</td><td>90</td><td>68</td><td>0</td><td>2</td><td>0</td><td>0</td><td>0</td><td>17</td><td>0</td><td>84</td><td>54</td><td>193</td><td>122</td></tr></tbody></table>

characteristics, each new variant (except per instance filtering selection) allows us to solve more instances. For *meshes*, the greatest improvement comes from the VrOH $\frac{dom}{wdeg}$ (V2) which allows us to solve all instances but one. For *randER*, each variant improves results. For *phase*, each variant but V1 improves results, the greatest improvement coming from random restarts and the VIOH (V3). In Fig. 8, we display the evolution of the cumulative number of instances solved by each variant with respect to time.

## 9 Experimental comparison with other exact approaches

Let us now compare LAD2025 with other exact approaches described in Section 2, i.e., PathLAD, RI, Glasgow, and PathLAD+ (we denote PathLAD+23 the solver described in [WJCL23] and PathLAD+24 the solver described in [WJC24]). We also include results of VF3 [CFSV18] which has been initially designed to solve the induced SIP, and has been recently extended to solve the non-induced SIP (called edge-induced SIP by authors of VF3).

In Table 5, we display the number of unsolved instances for each benchmark, each solver and two different time limits, i.e., 100s and 2000s. In Figure 9, we display box-plots of solving times for each benchmark and each solver.
<figure><img src="image_10.png" /><figcaption>Figure 9: Box-plots of solving times for each benchmark and each solver. Black squares correspond to $t_l^s$ (average solving time of the solved instances) and black dotted lines to $\bar{t}_l^s$ (lower bound of the average solving time obtained by bounding the run time of each unsolved instance with 2000).</figcaption></figure>

<table><caption>Table 6: Aggregate speed-ups between all couples of solvers for a time limit of 100s ($l = 100$) and a time limit of 2000s ($l = 2000$). Each cell at line $i$ and column $j$ displays $ASU_l(j)$ for $l \in \{100, 2000\}$ if $i$ solves more instances than $j$ within a time limit of $l$ seconds, and '-' otherwise.</caption><thead><tr><th rowspan="2"></th><th colspan="2">RI</th><th colspan="2">VF3</th><th colspan="2">PathLAD</th><th colspan="2">Glasgow</th><th colspan="2">PathLAD+23</th><th colspan="2">PathLAD+24</th></tr><tr><th>l=100</th><th>l=2000</th><th>l=100</th><th>l=2000</th><th>l=100</th><th>l=2000</th><th>l=100</th><th>l=2000</th><th>l=100</th><th>l=2000</th><th>l=100</th><th>l=2000</th></tr></thead><tbody><tr><td>VF3</td><td>2</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>PathLAD</td><td>31</td><td>493</td><td>30</td><td>468</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Glasgow</td><td>482</td><td>7,503</td><td>482</td><td>7,010</td><td>66</td><td>109</td><td>-</td><td>-</td><td>13</td><td>-</td><td>3</td><td>-</td></tr><tr><td>PathLAD+23</td><td>30</td><td>493</td><td>30</td><td>471</td><td>3</td><td>12</td><td>-</td><td>4</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>PathLAD+24</td><td>612</td><td>4,658</td><td>554</td><td>4,172</td><td>10</td><td>25</td><td>-</td><td>8</td><td>2</td><td>4</td><td>-</td><td>-</td></tr><tr><td>LAD2025</td><td>1,542</td><td>17,288</td><td>1,335</td><td>15,807</td><td>165</td><td>646</td><td>12</td><td>89</td><td>50</td><td>33</td><td>22</td><td>11</td></tr></tbody></table>
**Results on LV.** LAD2025 is the most successful solver on this benchmark, and it solves 50 (resp. 32) more instances than PathLAD+24, the second most successful solver, when the time limit is 100s (resp. 2000s). It also has the smallest average solving time, and the smallest third quartile. It has the second smallest median solving time, slightly larger than the one of RI (19ms for LAD2025 instead of 15ms for RI), but RI has a much higher average solving time (21.5s for RI instead of 3.4s for LAD2025). This comes from the fact that LV contains many easy instances that are very quickly solved by RI, whereas LAD2025 (as well as all constraint-based solvers) may need more time on these instances because it has to instantiate data structures (such as sparse sets used to manage domains, for example). These data structures are useless for solving easy instances, but they pay-off on many other instances and allow LAD2025 to solve much more instances than RI for longer time limits.

There are 28 instances that are solved only by LAD2025. In particular, instance g19-g20 (described in Section 3.1) is solved by LAD2025 in 65s whereas it is not solved by any of the other solvers within 2000s. There are only two instances (g29-g50 and g31-g50) that are not solved by LAD2025 whereas they are solved by another solver. More precisely, Glasgow is the only solver able to solve these two instances (that are infeasible), and it solves them without backtracking. This may come from the fact that Glasgow is the only solver that fully ensures the DC of edge constraints of supplemental graphs associated with $P_2$ motifs (as explained in Section 5.2). PathLAD+ and PathLAD only propagate unary $P_2$ labels and binary $P_2$ labels associated with pattern edges (and not $P_2$ labels associated with couples of non-adjacent pattern vertices). LAD2025 first ensures a weaker consistency than DC for $SND_{uu'}^{P_2,all}$ (as explained in Section 5.3), and then propagates Constraint 4 with a simple forward checking which is weaker than DC.

**Results on images.** This benchmarks only contains instances that are quickly solved by all solvers. RI is the fastest approach and its average solving time is 3ms. LAD2025 is the second fastest approach: its average solving time is 17ms, whereas the average solving times of PathLAD+24, Glasgow, PathLAD, PathLAD+23, and VF3 are 35ms, 85ms, 860ms, 901ms, and 989ms, respectively. Also, all quartiles of LAD2025 run times are smaller than those of all other approaches but RI.

On this benchmark, the run time of VF3 highly depends on whether the instance is feasible or not: there are 52 feasible instances and, for these instances, VF3 never needs more than 0.015s; there are 6250 infeasible instances and, for these instances, VF3 needs between 0.946s and 1.019s.

**Results on meshes.** This benchmark contains many easy instances, and a few hard instances that cannot be solved by some solvers within 2000s. All these hard instances have the same pattern graph (i.e., pattern3), and the degree of vertex 2 in this graph is equal to 131, whereas all other vertices for all patterns have much smaller degrees (never exceeding 10, except for two vertices that have a degree equal to 21). VrOHs that dynamically identify critical variables, such as $\frac{dom}{wdeg}$, greatly help solving instances that involve this pattern, and allow LAD2025 to solve all instances within a time limit of 2000s whereas the second most successful approach, PathLAD+24, fails at solving 6 instances within the same limit. The average solving time of LAD2025 is also much smaller.

This benchmark contains 2929 infeasible instances and 89 feasible instances. Neither VF3 nor RI are able to solve any of the feasible instances. If RI is able to solve many infeasible instances in a few milliseconds, VF3 needs at least 0.961s to solve these instances.

**Results on rand.** This benchmark only contains easy instances, that are quickly solved by all constraint-based solvers. However, if RI and VF3 perform very well on many instances of this benchmark, they fail at solving most of the infeasible scale-free instances (18 for RI and 19 for VF3) within a time limit of 2000s. Hence, if they have the smallest values for the first, second
<figure><img src="image_11.png" /><figcaption>Figure 10: Evolution of the cumulative number of instances solved by each solver with respect to time (with a logscale). The right figure is a zoom when $x \ge 0.05s$ and $y \ge 13400$ instances.</figcaption></figure>

and third quartiles, they also have the largest maximum values and average solving times. On this benchmark, LAD2025 has the smallest average solving time and maximum solving time.

**Results on randER.** On this benchmark, PathLAD has very poor performance compared to all other approaches: when the time limit is 2000s, it fails at solving 55 instances whereas LAD2025, PathLAD+24, and Glasgow are able to solve all instances, and PathLAD+23, RI, and VF3 fail at solving a few instances (from 2 to 6). Some instances contain rather dense graphs, and FC is much more effective than LAD on these instances.

Among the three solvers that are able to solve all instances, Glasgow has the smallest average solving time and quartile values. LAD2025 has the second best average solving times and the second best first, second, and third quartiles. However, the maximum solving time of LAD2025 is larger than the maximum solving time of PathLAD+24 (1706s for LAD2025 and 1321s for PathLAD+24).

**Results on phase.** This benchmark contains randomly generated instances that are expected to be hard because they are close to the phase transition region. Indeed, when considering a time limit of 2000s, more than 25% of these instances cannot be solved by any of the considered approaches, and three solvers (RI, VF3, and PathLAD) solve less than 25% of these instances. On this benchmark, LAD2025 is the most successful approach and it is able to solve 11 (resp. 35) more instances than PathLAD+24, the second most successful approach within 100s (2000s). LAD2025 also has the smallest average and median run times, as well as the smallest first quartile. However, it does not have the smallest minimum solving time (the smallest solving time of LAD2025 is 9ms whereas it is equal to 4ms for Glasgow and 6ms for PathLAD+24). This may come from the fact that the computation of labels is time consuming and useless on these instances.

**Impact of the time limit on the total number of solved instances.** When the time limit is set to 100s, the most successful solver is LAD2025, the second most successful solver is Glasgow (which solves 111 fewer instances than LAD2025), and the third most successful solver is PathLAD+24 (which solves 53 fewer instances than Glasgow). When the time limit is set to 2000s, LAD2025 still is the most successful solver, but the second most successful solver is PathLAD+24
807

(which solves 49 fewer instances than LAD2025), and the third most successful solver is Path-LAD+23 (which solves 35 fewer instances than PathLAD+23).

808

809

Aggregate speed-ups are displayed in Table 6, and we observe that, if LAD2025 is more efficient than all other approaches for the two considered time limits, ASUs strongly depend on the considered time limit. More precisely, when the time limit is 100s, LAD2025 is aggregately 12 times as fast as Glasgow, which is aggregately 3 times as fast as PathLAD+2024, which is aggregately twice as fast as PathLAD+23, which is aggregately 3 times as fast as PathLAD, which is aggregately 30 times as fast as VF3, which is aggregately twice as fast as RI. However, when the time limit is set to 2000s, LAD2025 is aggregately 11 times as fast as PathLAD+24, which is aggregately 4 times as fast as PathLAD+23, which is aggregately 4 times as fast as Glasgow, which is aggregately 109 times as fast as PathLAD, which is aggregately 468 times as fast as VF3, which is aggregately twice as fast as RI.

810

811

812

813

814

815

816

817

818

819

To allow us to compare solvers for other time limits than 100s and 2000s, we display the evolution of the cumulative number of solved instances with respect to time in Fig. 10. We do not display results of PathLAD and VF3 to avoid overloading the figure. For time limits smaller than 0.02s, RI is the most successful solver: it is able to solve more than 12,000 instances within a time limit of 0.02s. However, for larger time limits, it is strongly outperformed by LAD2025 which is the most successful solver for all time limits between 0.02s and 2000s. PathLAD+24 is more successful than PathLAD+23, for all time limits. However, the relation between PathLAD+24 and Glasgow depends on the time limit: Glasgow is more successful than PathLAD+24 for very short time limits (smaller than 4ms), and when the time limit is between 0.2s and 200s; PathLAD+24 is more successful than Glasgow when the time limit is between 4ms and 0.2s, and when the time limit is larger than 200s.

820

821

822

823

824

825

826

827

828

829

830

Finally, in Fig. 11, we compare the CPU times of LAD2025, PathLAD+24, Glasgow and RI on a per instance basis. This shows the complementarity of these approaches. LAD2025 is very often the most efficient solver on *LV* and *meshes* instances, and also on many hard *phase* instances. However, some *phase* instances are solved in much less than 10s by PathLAD+24 and Glasgow whereas LAD2025 needs 10s or so (due to the computation of labels before starting the search, which is time consuming and useless on these instances). Glasgow is very often the most efficient solver on *randER* instances. RI is very often the most efficient solver on *images* instances.

831

832

833

834

835

836

837

# 10 Discussion

838

In this paper, we have seen that implementation does matter: PathLAD and V0 build the same search tree, but V0 is 10 times as fast as PathLAD when looking at the aggregate speed-up. However, a careful implementation is not enough to obtain state-of-the-art results. Some improvements described in this paper are well-known techniques used in CP solvers, i.e., weight-based variable ordering heuristics and random restarts with nogood recording. The fact that these techniques improve our solver demonstrates their robustness. Some other improvements are coming from other SIP solvers: supplemental labels based on 2-paths and on maximum cliques in vertex neighbourhoods, and degree-biased value ordering heuristic. Finally, two improvements are new: supplemental labels based on clique enumeration and the fact that we use the largest consistent assignment as a value ordering heuristic. Also, even though using supplemental labels based on 2-path enumeration is not new, we have introduced a new algorithm to efficiently propagate SND<sub>uu'</sub><sup>P<sub>2</sub>,all</sup> constraints before starting the search. All these improvements allow LAD2025 to be the most successful approach, for all time limits ranging from 0.02s to 2000s. However, LAD2025 is complementary to Glasgow, which outperforms it on *randER* instances, and RI, which outperforms it on *images* instances. As shown in [KMS16], we could consider a portfolio approach to take advantage of this complementarity.

839

840

841

842

843

844

845

846

847

848

849

850

851

852

853

854

LAD2025 may be easily extended to the case of labelled graphs, such that labels are associated
<figure><img src="image_12.png" /><figcaption>Figure 11: Scatter plots to compare LAD2025, PathLAD+24, Glasgow, and RI.</figcaption></figure>

with vertices and edges, and the sub-isomorphism function must map vertices and edges with compatible labels: labels associated with vertices are used to filter domains before starting the search whereas labels associated with edges are used to remove incompatible edges from bipartite graphs associated with LAD constraints, just like what is done with supplemental labels. LAD2025 may also be easily extended to the case of directed graphs, as directed edges may be handled by adding labels to edges (where labels indicate whether the edge is incoming, outgoing, or both incoming and outgoing).

In this paper, we have considered the non-induced SIP. Our approach could be extended to the
induced case, such that two non-adjacent pattern vertices must be mapped to non-adjacent target vertices. In this case, we may also define LAD constraints associated with non-adjacent vertices: given a pattern vertex $u$ and a target vertex $u' \in D_u$, we may not only ensure that each vertex $v$ adjacent to $u$ can be mapped to a different vertex that is adjacent to $u'$ and belongs to $D_v$, like in the non-induced case, but also ensure that each non-adjacent vertex $v \in V \setminus \mathcal{N}(u)$ can be mapped to a different vertex that is not adjacent to $u'$ and belongs to $D_v$. However, a trade-off has to be found between the strength of the filtering and the number of LAD constraints to propagate, depending on the size of both $\mathcal{N}(u)$ and $V \setminus \mathcal{N}(u)$. We plan to investigate this in future work.

LAD2025 has been extensively tested by comparing its results with those of other exact solvers on all our benchmarks. However, this does not ensure that there is no bug in the code and that LAD2025 actually outputs correct results for instances that do not belong to our test set. Hence, we plan to extend LAD2025 with proof logging in order to produce certificates of correctness that can be checked by independent programs, as proposed in [GMM$^+$24], for example.

**Acknowledgements:** The author would like to express her gratitude to Donald E. Knuth for his careful analysis of LAD that helped her to improve it and also greatly motivated her to write this paper.

This work was partially supported by the PEPR IA HOLIGRAIL project (ANR-23-PEIA-0010).

## References

[ADH$^+$19] Blair Archibald, Fraser Dunlop, Ruth Hoffmann, Ciaran McCreesh, Patrick Prosser, and James Trimble. Sequential and parallel solution-biased search for subgraph algorithms. In Louis-Martin Rousseau and Kostas Stergiou, editors, *16th International Conference on Integration of Constraint Programming, Artificial Intelligence, and Operations Research, CPAIOR 2019*, volume 11494 of *Lecture Notes in Computer Science*, pages 20–38. Springer, 2019.

[ALM$^+$14] G. Audemard, C. Lecoutre, M. Samy Modeliar, G. Goncalves, and D. Porumbel. Scoring-based neighborhood dominance for the subgraph isomorphism problem. In *20th Int. Conf. on Principles and Practice of Constraint Programming CP*, pages 125–141, 2014.

[Ber73] Claude Berge. *Graphs and Hypergraphs*. 1973.

[BG17] V. Bonnici and R. Giugno. On the variable ordering in subgraph isomorphism algorithms. *IEEE/ACM Trans. Comput. Biology Bioinform.*, 14(1):193–203, 2017.

[BHLS04] Frédéric Boussemart, Fred Hemery, Christophe Lecoutre, and Lakhdar Sais. Boosting systematic search by weighting constraints. In Ramón López de Mántaras and Lorenza Saitta, editors, *Proceedings of the 16th European Conference on Artificial Intelligence, ECAI'2004, including Prestigious Applicants of Intelligent Systems, PAIS 2004, Valencia, Spain, August 22-27, 2004*, pages 146–150. IOS Press, 2004.

[CFSV04] L. P. Cordella, P. Foggia, C. Sansone, and M. Vento. A (sub)graph isomorphism algorithm for matching large graphs. *IEEE Trans. Pattern Anal. Mach. Intell.*, 26(10):1367–1372, 2004.

[CFSV18] V. Carletti, P. Foggia, A. Saggese, and M. Vento. Challenging the time complexity of exact subgraph isomorphism for huge and dense graphs with VF3. *IEEE Trans. Pattern Anal. Mach. Intell.*, 40(4):804–818, 2018.
[CLRS09] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. *Introduction to Algorithms*. The MIT Press, 2009.

[DCS18] Emir Demirovic, Geoffrey Chu, and Peter J. Stuckey. Solution-based phase saving for CP: A value-selection heuristic to simulate local search behavior in complete solvers. In *Principles and Practice of Constraint Programming - 24th International Conference, CP 2018, Proceedings*, volume 11008 of *Lecture Notes in Computer Science*, pages 99-108. Springer, 2018.

[DSdlH$^{+}$11] G. Damiand, C. Solnon, C. de la Higuera, J.-C. Janodet, and E. Samuel. Polynomial algorithms for subisomorphism of nD open combinatorial maps. *Computer Vision and Image Understanding (CVIU)*, 115(7):996-1010, 2011.

[DSFSV03] M. De Santo, P. Foggia, C. Sansone, and M. Vento. A large database of graphs and its use for benchmarking graph isomorphism algorithms. *Pattern Recogn. Lett.*, 24(8):10671079, 2003.

[FF56] L. R. Ford and D. R. Fulkerson. Maximal flow through a network. *Canadian Journal of Mathematics*, 8:399-404, 1956.

[GMM$^{+}$24] Stephan Gocht, Ciaran McCreesh, Magnus O. Myreen, Jakob Nordström, Andy Oertel, and Yong Kiam Tan. End-to-end verification for subgraph solving. In Michael J. Wooldridge, Jennifer G. Dy, and Sriraam Natarajan, editors, *Thirty-Eighth AAAI Conference on Artificial Intelligence, AAAI 2024, Thirty-Sixth Conference on Innovative Applications of Artificial Intelligence, IAAI 2024, Fourteenth Symposium on Educational Advances in Artificial Intelligence, EAAI 2014, February 20-27, 2024, Vancouver, Canada*, pages 8038-8047. AAAI Press, 2024.

[GMN08] Ian P. Gent, Ian Miguel, and Peter Nightingale. Generalised arc consistency for the all different constraint: An empirical survey. *Artif. Intell.*, 172(18):1973-2000, 2008.

[GSK98] Carla P. Gomes, Bart Selman, and Henry A. Kautz. Boosting combinatorial search through randomization. In Jack Mostow and Chuck Rich, editors, *Proceedings of the Fifteenth National Conference on Artificial Intelligence and Tenth Innovative Applications of Artificial Intelligence Conference, AAAI 98, IAAI 98, July 26-30, 1998, Madison, Wisconsin, USA*, pages 431-437. AAAI Press / The MIT Press, 1998.

[HK73] John E. Hopcroft and Richard M. Karp. An $n^{5/2}$ algorithm for maximum matchings in bipartite graphs. *SIAM Journal on Computing*, 2(4):225-231, 1973.

[HMN$^{+}$18] R. Hoffmann, C. McCreesh, S. N. Ndiaye, P. Prosser, C. Reilly, C. Solnon, and J. Trimble. Observations from parallelising three maximum common (connected) subgraph algorithms. In *Integration of Constraint Programming, Artificial Intelligence, and Operations Research - 15th International Conference, CPAIOR*, volume 10848 of LNCS, pages 298-315. Springer, 2018.

[KM21] Sonja Kraiczy and Ciaran McCreesh. Solving graph homomorphism and subgraph isomorphism problems faster through clique neighbourhood constraints. In *Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21*, pages 1396-1402. International Joint Conferences on Artificial Intelligence Organization, 2021.
[KMS16] L. Kotthoff, C. McCreesh, and C. Solnon. Portfolios of subgraph isomorphism algorithms. In *10th Int. Conf. on Learning and Intelligent Optimization LION*, volume 10079 of LNCS, pages 107-122. Springer, 2016.

[Knu93] D. E. Knuth. *The Stanford GraphBase - a platform for combinatorial computing*. ACM, 1993.

[Knu22] Donald Ervin Knuth. *The art of computer programming, Volume 4B, Combinatorial Algorithms, Part 2*. Addison-Wesley, 2022.

[Knu25] Donald Ervin Knuth. *The art of computer programming, Volume 4, Fascicle 7A: Constraint Satisfaction*. Addison-Wesley, 2025.

[Kot13] Lars Kotthoff. Llama: leveraging learning to automatically manage algorithms. *ArXiv preprint arXiv:1306.1031*, 2013.

[LBCBP$^+$25] Sulian Le Bozec-Chiffoleau, Nicolas Beldiceanu, Charles Prud'homme, Gilles Simonin, and Xavier Lorca. Bimodal depth-first search for scalable gac for alldifferent. In James Kwok, editor, *Proceedings of the Thirty-Fourth International Joint Conference on Artificial Intelligence, IJCAI-25*, pages 2610-2618. International Joint Conferences on Artificial Intelligence Organization, 2025.

[LBLA$^+$09] Pierre Le Bodic, Hervé Locteau, Sébastien Adam, Pierre Héroux, Yves Lecourtier, and Arnaud Knippel. Symbol Detection Using Region Adjacency Graphs and Integer Linear Programming. In International Association for Pattern Recognition TC-10 and TC-11, editors, *International Conference on Document Analysis and Recognition*. Institute of Electrical and Electronics Engineers, 2009.

[LCdSMSSL13] Vianney Le Clément de Saint-Marcq, Pierre Schaus, Christine Solnon, and Christophe Lecoutre. Sparse-Sets for Domain Implementation. In *CP workshop on Techniques for Implementing Constraint programming Systems (TRICS)*, pages 1-10, 2013.

[LSTV07] Christophe Lecoutre, Lakhdar Sais, Sébastien Tabary, and Vincent Vidal. Nogood recording from restarts. In *IJCAI 2007, Proceedings of the 20th International Joint Conference on Artificial Intelligence*, pages 131-136, 2007.

[LSZ16] Jimmy H. M. Lee, Christian Schulte, and Zichen Zhu. Increasing nogoods in restart-based search. In Dale Schuurmans and Michael P. Wellman, editors, *Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence*, pages 3426-3433. AAAI Press, 2016.

[LV02] J. Larrosa and G. Valiente. Constraint satisfaction algorithms for graph pattern matching. *Mathematical. Structures in Comp. Sci.*, 12(4):403-422, 2002.

[LYL21] Hongbo Li, Minghao Yin, and Zhanshan Li. Failure Based Variable Ordering Heuristics for Solving CSPs. In *27th International Conference on Principles and Practice of Constraint Programming (CP 2021)*, volume 210 of Leibniz International Proceedings in Informatics (LIPIcs), pages 9:1-9:10. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021.

[MdG23] Etienne Macé de Gastines. *Modélisation, approche polyédrale et approche arborescente pour des problèmes de comparaison de graphes*. Theses, Normandie Université, 2023.
<table>
  <tbody>
    <tr>
      <td>990</td>
      <td>[MP15]</td>
      <td>C. McCreesh and P. Prosser. A parallel, backjumping subgraph isomorphism algorithm using supplemental graphs. In Principles and Practice of Constraint Programming, volume 9255 of LNCS, pages 295–312. Springer, 2015.</td>
    </tr>
    <tr>
      <td>991</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>992</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>993</td>
      <td>[MPST18]</td>
      <td>C. McCreesh, P. Prosser, C. Solnon, and J. Trimble. When Subgraph Isomorphism is Really Hard, and Why This Matters for Graph Databases. Journal of Artificial Intelligence Research, 61:723 – 759, 2018.</td>
    </tr>
    <tr>
      <td>994</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>995</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>996</td>
      <td>[MPT20]</td>
      <td>Ciaran McCreesh, Patrick Prosser, and James Trimble. The glasgow subgraph solver: Using constraint programming to tackle hard subgraph isomorphism problem variants. In Graph Transformation - 13th International Conference, ICGT 2020, Proceedings, volume 12150 of Lecture Notes in Computer Science, pages 316–324. Springer, 2020.</td>
    </tr>
    <tr>
      <td>997</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>998</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>999</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1000</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1001</td>
      <td>[PD07]</td>
      <td>Knot Pipatsrisawat and Adnan Darwiche. A lightweight component caching scheme for satisfiability solvers. In Theory and Applications of Satisfiability Testing - SAT 2007, 10th International Conference, Proceedings, volume 4501 of Lecture Notes in Computer Science, pages 294–299. Springer, 2007.</td>
    </tr>
    <tr>
      <td>1002</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1003</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1004</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1005</td>
      <td>[Reg94]</td>
      <td>J.-C. Regin. A filtering algorithm for constraints of difference in CSPs. In Proc. 12th Conf. American Assoc. Artificial Intelligence, volume 1, pages 362–367. Amer. Assoc. Artificial Intelligence, 1994.</td>
    </tr>
    <tr>
      <td>1006</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1007</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1008</td>
      <td>[SDdlHJ15]</td>
      <td>C. Solnon, G. Damiand, C. de la Higuera, and J.-C. Janodet. On the complexity of submap isomorphism and maximum common submap problems. Pattern Recognition, 48(2):302–316, 2015.</td>
    </tr>
    <tr>
      <td>1009</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1010</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1011</td>
      <td>[Sol10]</td>
      <td>C. Solnon. Alldifferent-based filtering for subgraph isomorphism. Artif. Intell., 174(12-13):850–864, 2010.</td>
    </tr>
    <tr>
      <td>1012</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1013</td>
      <td>[Sol19]</td>
      <td>Christine Solnon. Experimental evaluation of subgraph isomorphism solvers. In Graph-Based Representations in Pattern Recognition - 12th IAPR-TC-15 International Workshop, GbRPR 2019, Proceedings, volume 11510 of Lecture Notes in Computer Science, pages 1–13. Springer, 2019.</td>
    </tr>
    <tr>
      <td>1014</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1015</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1016</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1017</td>
      <td>[Tar72]</td>
      <td>Robert Tarjan. Depth-first search and linear graph algorithms. SIAM Journal on Computing, 1(2):146–160, 1972.</td>
    </tr>
    <tr>
      <td>1018</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1019</td>
      <td>[Ull76]</td>
      <td>J. R. Ullmann. An algorithm for subgraph isomorphism. J. ACM, 23(1):31–42, 1976.</td>
    </tr>
    <tr>
      <td>1020</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1021</td>
      <td>[VP17]</td>
      <td>Julien Vion and Sylvain Piechowiak. Une simple heuristique pour rapprocher DFS et LNS pour les COP. In Journées Francophones de Programmation par Con-traintes (JFPC), 2017, Proceedings, pages 39–44, 2017.</td>
    </tr>
    <tr>
      <td>1022</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1023</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1024</td>
      <td>[WJC24]</td>
      <td>Yiyuan Wang, Chenghou Jin, and Shaowei Cai. Pathlad+: Towards effective exact methods for subgraph isomorphism problem. Artificial Intelligence, 337:104219, 2024.</td>
    </tr>
    <tr>
      <td>1025</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1026</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1027</td>
      <td>[WJCL23]</td>
      <td>Yiyuan Wang, Chenghou Jin, Shaowei Cai, and Qingwei Lin. Pathlad+: An improved exact algorithm for subgraph isomorphism problem. In Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence, IJCAI, pages 5639–5647. ijcai.org, 2023.</td>
    </tr>
    <tr>
      <td>1028</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1029</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>1030</td>
      <td></td>
      <td></td>
    </tr>
  </tbody>
</table>
1031

[WLPT19]

1032

1033

1034

Hugues Wattez, Christophe Lecoutre, Anastasia Paparrizou, and Sébastien Tabary.
Refining constraint weighting. In *31st IEEE International Conference on Tools with Artificial Intelligence, ICTAI 2019, Portland, OR, USA, November 4-6, 2019*, pages 71–77. IEEE, 2019.

1035

[YXW20]

1036

1037

1038

1039

1040

Roland H. C. Yap, Wei Xia, and Ruiwei Wang. Generalized arc consistency algorithms for table constraints: A summary of algorithmic ideas. In *The Thirty-Fourth AAAI Conference on Artificial Intelligence, AAAI 2020, The Thirty-Second Innovative Applications of Artificial Intelligence Conference, IAAI 2020, The Tenth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2020, New York, NY, USA, February 7-12, 2020*, pages 13590–13597. AAAI Press, 2020.

1041

[ZDS10]

1042

S. Zampelli, Y. Deville, and C. Solnon. Solving subgraph isomorphism problems with constraint programming. *Constraints*, 15(3):327–353, 2010.