PathLAD+: An Improved Exact Algorithm for Subgraph Isomorphism Problem

1 Introduction

The (non-induced) subgraph isomorphism problem (SIP), which is also known as the subgraph matching problem, involves deciding if there exists a copy of a pattern graph in a target graph. As one of the basic concepts of graph theory, the SIP can be seen as a generalization of both the maximum clique problem and the problem of testing whether a graph contains a Hamiltonian cycle. Recently, the SIP has been used in various domains, such as symbol recognition [Llad´os et al., 2001], social networks [Snijders et al., 2006], computer vision [Damiand et al., 2011], biochemical data [Bonnici et al., 2013], RDF query processing [Kim et al., 2015] and graph databases [Wang et al., 2022]. For example, the SIP has also been used in the field of cheminformatics to search for some similarities between chemical compounds from their structural formula [Ohlrich et al., 1993].

∗Corresponding author

It is well known that the SIP is NP-complete in the general case [Johnson and Garey, 1979]. For the optimized version of the SIP, i.e., the maximum common induced subgraph (MCS) problem, many methods have been presented to deal with the MCS problem [McCreesh et al., 2016a; McCreesh et al., 2017; Liu et al., 2020; Gocht et al., 2020; Zhou et al., 2022; Liu et al., 2022]. The MCS approaches can be directly used into solving the SIP, but they usually have poor performance practically due to the characteristics of the SIP as the decision problem. Thus, up to now, the SIP has been still considered as a challenging problem.

In the last decade, lots of researchers focused on designing several exact methods to address the SIP [Zampelli et al., 2010; Solnon, 2010; Bonnici et al., 2013; Audemard et al., 2014; McCreesh and Prosser, 2015; Carletti et al., 2017; McCreesh et al., 2018; Archibald et al., 2019; Solnon, 2019; McCreesh et al., 2020]. We list some representative solvers for the SIP as below. An early algorithm for the SIP named VF2 was proposed, which used a state space representation of the matching process and introduced a set of five feasibility rules for pruning the search tree [Cordella et al., 2004]. Bonnici et al. [2013] developed a new search strategy called RI based on the pattern graph topology, which significantly reduced the search space without using any complex pruning rules or reduction procedures. Solnon [2010] introduced a new filtering algorithm called LAD based on local alldifferent constraints. The LAD algorithm was further improved by combining the local all-different constraints with the exploitation of path length properties, resulting in the PathLAD algorithm [Kotthoff et al., 2016]. Very recently, Kraiczy and McCreesh [2021] improved the Glasgow [McCreesh et al., 2020] by using a new form of filtering based upon clique-finding and designed a new algorithm called Glasgow+Clq. According to the literature, the current best algorithm for the SIP is Glasgow+Clq [Kraiczy and McCreesh, 2021].

1.1 Our Contribution

Motivated to contribute to further improving the performance of SIP, in this work, we choose the PathLAD algorithm as a baseline algorithm. Our proposed algorithm is divided into two parts, including a probing search procedure and a main search procedure. Below are three main novel ideas in our proposed algorithm.

First, we propose a probing search procedure in which the algorithm tries several times to judge whether the pattern graph is isomorphic to a subgraph of the target graph quickly. It has two main purposes. On the one hand, if this procedure can successfully solve a given instance, we can obtain an outcome within a short time. On the other hand, if this procedure cannot get any outcomes (i.e., reaching cutoff time), instead of using some traditional restart mechanisms, we can still extract the information from this search procedure to guide a main search procedure. That is, the algorithm learns about the useful search information between pattern vertices and target vertices and then uses this information in our proposed main procedure.

Second, we design a new matching ordering method to decide which target vertex from the domain of the corresponding pattern vertex is selected. Recently, several matching ordering methods were proposed [Archibald et al., 2019; Wang et al., 2022]. For example, Archibald et al. [2019] found that it is effective to preferentially select vertices with high degree values when selecting a matching target vertex. We follow this line of research by attempting to apply the degree information of target vertices in the matching process. At the same time, we also use the useful search information generated from the probing search procedure as another matching criterion. Thus, our proposed matching ordering method considers the above two principles, resulting in a novel scoring function denoted as oscore used in the matching process.

Third, we present an adaptive propagation method to dynamically use different strong propagation methods for the SIP. Previous algorithms have always applied strong propagation methods to remove some unnecessary vertices from the corresponding domains, but these propagation methods need to cost lots of run time in practice, which reduces the performance of these algorithms. In some cases, instead of strong propagation methods, some weak propagation methods can make backtrack quickly or reduce the size of the corresponding domain effectively on some branches. Conversely, the performance of these algorithms would be also bad if they don't use any strong propagation methods because the algorithms fail to backtrack some branches immediately. Based on the above considerations, we analyze the properties of strong propagation methods and combine the search information generated from the main search procedure to dynamically employ different strong propagation methods during the search. To our best knowledge, it is the first time that different propagation methods are dynamically used to accelerate the search procedure for addressing SIP.

By incorporating these ideas, we develop an improved exact algorithm for the SIP called PathLAD+. Extensive experiments are carried out to evaluate PathLAD+ on the benchmarks used in the literature. Experimental results show that PathLAD+ outperforms four state-of-the-art SIP algorithms for all the benchmarks. In addition, our experimental analyses report that the proposed strategies play important roles in the outstanding performance of our proposed algorithm.

In the next section, we introduce some necessary background knowledge. After that, we present our proposed algorithm and its components. Experimental results are shown in Section 4. Finally, we make conclusions.

2 Preliminaries

2.1 Basic Definitions and Notations

Let G = (V, E) be an undirected graph where a vertex set is V = {v1, v2, . . . , vn} and an edge set E = {e1, e2, . . . , em}. Each edge e is a 2-element subset of V , i.e., e = (v, u). For an edge e = (v, u), we say vertices v and u are the endpoints of edge e. For a vertex v ∈ V , the neighborhood of vertex v is denoted as NG(v) = {u | (v, u) ∈ E} and its degree is denoted as degG(v) = |NG(v)|. A finite walk is a sequence of edges (e′1 , e′2 , . . . , e′q-1 ) for which there is a sequence of vertices (v′1 , v′2 , . . . , v′q ) such that e′i = (v′i , v′i+1 ) for i ∈ [1, q - 1]. A path is a finite walk in which all vertices and all edges are distinct, denoted as ζ G = (v′1 , v′2 , . . . , v′q ). The length of ζ G is denoted as |ζ G| = q.

Given a pattern graph Gp = (Vp,Ep) and a target graph Gt = (Vt, Et), the SIP is to decide whether Gp is isomorphic to some subgraph of Gt. Formally, the aim of the SIP is to obtain an injective matching f : Vp → Vt that associates a different target vertex to each pattern vertex, and preserves pattern edges, i.e., (f(v), f(u)) ∈ Et for ∀(v, u) ∈ Ep. It is noted that the subgraph is not necessarily induced, which means that two pattern vertices that are not linked by an edge may be matched to two target vertices that are linked by an edge. During the search procedure, the current matched list of pattern and target pairs is denoted as D = {{vp 1 , vt1 }, . . . , {vp r , vtr }}. For a pattern vertex vp ∈ Vp, the domain of vertex vp is defined as the set of target vertices that may be matched to vp , i.e., Dom(vp ) = {vt1 , vt2 , . . . , vt l }, and the size of its domain is |Dom(vp )| = l.

2.2 Some PropagationMethods for the SIP

Recently, three filtering propositions [Zampelli et al., 2010; McCreesh and Prosser, 2015] have been used in Glasgow+Clq [Kraiczy and McCreesh, 2021]. We first introduce three propositions that are used to judge whether the pattern vertices can be matched to the corresponding target vertices.

Proposition 1. Given a pattern graph Gp = (Vp,Ep) and a target graph Gt = (Vt,Et), if vp ∈ Vp can be matched to vt ∈ Vt (i.e., f(vp ) = vt ), it must satisfy degGp (vp ) ≤ degGt (vt ).

Proposition 2. Given a pattern graph Gp = (Vp, Ep) and a target graph Gt = (Vt,Et), if vp ∈ Vp can be matched to vt ∈ Vt (i.e., f(vp ) = vt ), it must satisfy the i-th value of ND(vt ) is not less than the same position of ND(vp ) where ND(vp ) = {degGp (up ) | up ∈ NGp (vp )}, ND(vt ) = {degGt (ut ) | ut ∈ NGt (vt )} and the positions of elements in ND(vp ) and ND(vt ) both are arranged in a descending order of the degree values.

Proposition 3. Given a pattern graph Gp = (Vp,Ep) and a target graph Gt = (Vt,Et), if vp and up in Vp can be matched to vt and ut in Vt (i.e., f(vp ) = vt and f(up ) = ut ) as well as Pathp and Patht are not empty, it must satisfy |Pathp | ≤ |Patht | where Pathp = {ζ p | ζ p = (vp , . . . , up ), |ζ p | = 3} and Patht = {ζ t | ζ t = (vt , . . . , ut ), |ζ t | = 3} .

Some propagation methods of difference constraints [Solnon, 2010] are used in the PathLAD [Kotthoff et al., 2016], which are shown as below.

• Vertex constraint denoted as FC(Diff): whenever a pattern vertex vp is matched to a target vertex vt , FC(Diff) removes vt from the domains of all nonmatched pattern vertices. The time complexity of FC(Diff) is O(|Vp|).
• Edge constraint denoted as FC(Edges): whenever a pattern vertex vp is matched to a target vertex vt , FC(Edges) removes any target vertex not adjacent to vt from the domain of every pattern vertex adjacent to vp . The time complexity of FC(Edges) is O(degGp (vp ) · |Vt|).
• Global neighborhood constraint denoted as GAC(allDiff): It ensures that all pattern vertices can be assigned to different target vertices. In detail, if a set of k pattern vertices can be found with only k target vertices among the domains of their pattern vertices, then those target vertices can be removed from the domains of other pattern vertices. The time complexity of GAC(allDiff) is O(|Vp|2 · |Vt|2 ).
• Filtering method denoted as LAD-filtering: for vt ∈ Dom(vp ), a bipartite graph is defined as G(vp ,vt ) = (NGp (vp ),NGt (vt ),E(vp ,vt )) where E(vp ,vt ) = {(v′ , u′ )|v′ ∈ NGp (vp ), u′ ∈ NGt (vt ), u′ ∈ Dom(v′ )}. If there does not exist a matching of the bipartite graph G(vp ,vt ) that covers NGp (vp ), the pattern vertices adjacent to vp cannot be matched to all different target vertices and thus vt can be removed from Dom(vp ). The time complexity of LAD-filtering is O(|Vp| · |Vt| · deg2Gp (vp ) · deg2Gt (vt )).

Note that two strong propagation methods LAD-filtering and GAC(allDiff) are implemented by the matching algorithm Hopcroft and Karp, and more details can be seen [Solnon, 2010].

3 The PathLAD+ Algorithm

This section describes the proposed PathLAD+ algorithm in Algorithm 1. Details of important functions in PathLAD+ will be presented in the following subsections. We use switchL and switchA to control whether the proposed algorithm uses LAD-filtering and GAC(allDiff), respectively. Meanwhile, nbnodes records the sum of call times and backtrack times of SearchSIP, nbfail records the number of backtrack times of SearchSIP, and Nb is used in our proposed adaptive propagation method. The output value st has three values: true means that the algorithm can return a successful matched list; false means the pattern graph is not isomorphic

Algorithm 1 PathLAD+

Input: Pattern graph Gp, target graph Gt and the cutoff time Output: outcome st

• 1: reduce the domain of pattern vertices based on Propositions 1 and 2;
• 2: if some domains become empty then
• 3: return false;
• 4: end if
• 5: nbnodes := nbfail := 0 and Nb := +∞;
• 6: switchL := switchA := 1;
• 7: st := unknown;
• 8: ProSearch(Gp,Gt);
• 9: return SearchSIP(Gp,Gt, ∅, cutoff);

to any subgraph of the target graph; unknown means that the algorithm cannot solve a given instance within a cutoff time. In the beginning, the algorithm reduces the domain of pattern vertices in a given pattern graph according to Propositions 1 and 2. If any domain becomes an empty set, the algorithm returns false. Otherwise, six variables are initialized accordingly (Lines 5-7). Then, the proposed algorithm can be divided into two procedures, including a probing search procedure (ProSearch in Line 8) and a main search procedure (SearchSIP in Line 9).

3.1 The Search Framework for SIP

The main function SearchSIP is shown in Algorithm 2, which is a recursive function. The input variable D is denoted as an already-matched list. If all vertices in the pattern graph are matched, which means that the algorithm has found a matched list for all pattern vertices, the algorithm returns true (Lines 1-2). Otherwise, if the time limit is reached, the algorithm returns unknown (Lines 3-4). The value of nbnodes is increased by 1 (Line 6). The algorithm chooses a non-matched pattern vertex vp i with the smallest domain size, breaking ties by picking the one with the biggest degree value (Line 7). Afterward, the algorithm arranges the positions of target vertices in Dom(vp i ) based on a novel matching ordering method (i.e., oscore), which will be introduced in Section 3.3 (Line 8). In Lines 9-23, the algorithm tries to match each target vertex in Dom(vp i ) orderly. Before executing Line 9, the algorithm will store the domain of all pattern vertices. In Line 14, the algorithm restores the domain of all pattern vertices to their previous saved domain in Line 9. The algorithm orderly tries to match a target vertex vt i in theDom(vp i ) to the selected pattern vertex vp i (Line 10). In each time, the algorithm reduces the domain of each non-matched pattern vertex based on an adaptive propagation method APM, which will be mentioned in Section 3.4 (Line 11). If the domain of some pattern vertex becomes empty, which means that vt i cannot be matched to vp i , the algorithm will restore and then continue to select the next target vertex in Dom(vp i ) (Lines 1216). The corresponding values of nbfail and nbnodes will be increased by 1 (Line 13). If the algorithm doesn't obtain any empty domains, it will search for the next pattern vertex (Line 17). st stores the backtracking result of SearchSIP. If st equals false, the algorithm needs to restore the related domains and then continue to select the next target vertex (Lines

Algorithm 2 SearchSIP

• Input: Pattern graph Gp, target graph Gt, an alreadymatched list of pattern and target pairs D = {{vp 1 , vt1 }, . . . , {vp i-1 , vt i-1 }} and the cutoff time

Output: outcome st

• 1: if all the pattern vertices have been matched to respective target vertices then
• 2: return true;
• 3: else if elapsed time > cutoff then
• 4: return unknown;
• 5: end if
• 6: nbnodes++;
• 7: select a vertex vp i with the smallest domain size |Dom(vp i )| from all non-matched pattern vertices, breaking ties by picking the one with the biggest degree value; /∗ recording info values in ProSearch, see Sec. 3.2 ∗/
• 8: sort the position of vertices in Dom(vp i ) based on the descending order of oscore values; /∗ see Sec. 3.3 ∗/
• 9: for each target vertex vt i ∈ Dom(vp i ) satisfying Proposition 3 do
• 10: match vt i to vp i ;
• j ) for each non-matched pat-
• 11: Dom(vp j ) := APM(Gt, vp tern vertex vp j ;
• 12: if some domains become empty then
• 13: nbfail++ and nbnodes++;
• 14: restore the domain of some pattern vertices;
• 15: continue;
• 16: end if
• 17: st := SearchSIP(Gp,Gt,D∪ {vp i , vt i }, cutoff);
• 18: if st == false then
• 19: restore the domain of some pattern vertices;
• 20: else
• 21: return st;
• 22: end if
• 23: end for
• 24: return false;

18-19).

3.2 The Probing Search Procedure for SIP

Before calling a main search procedure, a probing search procedure ProSearch plans to use less time to try to successfully solve an instance. If ProSearch can successfully solve an instance, the algorithm can return a matched list of pattern and target pairs quickly or can determine that the pattern graph is not isomorphic to a subgraph of the target graph. Otherwise, the algorithm can also grasp some useful information, denoted as info in our algorithm, from this search procedure, i.e., recording which target vertices are included in the domain of the corresponding pattern vertex during this procedure. It means that the information obtained from ProSearch can reflect which vertex pair has more potential.

The specific way of updating info values is presented as follows. At first, the info value of each pair of pattern and target vertices is initialized to 0. In the search procedure, assuming that we select a pattern vertex vp , we scan all target vertices in the domain of vp (Line 7 in Algorithm 2). For each target vertex vt i in Dom(vp ), the info({vp , vt }) is increased by 1.

The proposed probing search procedure ProSearch works as follows. During the search procedure of ProSearch, the algorithm uses four propagation methods and the third proposition which have already been introduced in Section 2.2 to reduce the domain of pattern vertices. ProSearch has two search modes. In the first mode, to explore the vertices in the deep depth of the search tree, the algorithm runs SearchSIP with a cutoff time of 10 seconds without sorting the domains in any way (i.e., the default lexicographical order). Because some SIP instances can be found with only a small number of conflicts, the heavy commitment to early branching choices made by backtracking search can be extremely costly for these instances [Archibald et al., 2019]. Based on the above consideration, in the second mode, the algorithm runs SearchSIP 20 times with a cutoff time of 1 second each time. To diversify early branch selections, the second mode sorts the position of target vertices in the corresponding domains randomly each time.

According to preliminary experiments, we found updating info during the main search procedure caused some vertex pairs with high info values and led to the poor performance. Thus, we restrict updating info only in ProSearch.

3.3 ANovel Matching OrderingMethod

In the search procedure of SearchSIP, among non-matched pattern vertices, we select a pattern vertex with the smallest domain size. After choosing a pattern vertex, the next key step is how to select a target vertex from the domain of the selected pattern vertex. Whatever matching ordering method is used, the method will only affect the performance for some instances that have a successful matched list, whereas it has no influence on some instances where a given pattern graph is not isomorphic to any subgraph of the target graph because a complete search must be performed.

During the probing search procedure, we use info to collect useful information on the relationship between pattern and target vertices. After some pattern vertices have already matched to different target vertices, we assume that a target vertex vt i is often included in the domain of a pattern vertex Dom(vp ), i.e., info({vp , vt i }) with a high value. We believe that vt i has more potential to match vp compared to other target vertices because a matched pair {vp , vt i } would bring few conflicts. This means that some other pattern vertices are more likely to successfully find the corresponding target vertices in the following search when vt i matches vp . At the same time, we consider the structure information of the target graph such as the degree value in our proposed matching ordering method.

As a result, we have the notion of a novel ordering score, which is formally defined as follows.

Definition 1. For a pattern graph Gp = (Vp, Ep) and a target graph Gt = (Vt,Et), the ordering score function, denoted as oscore is a function on vp ∈ Vp and vt ∈ Dom(vp ) such that

In our proposed matching ordering method, when the algorithm chooses a pattern vertex vp , the positions of target vertices in the domain of a selected pattern are arranged in a descending order of the oscore values (Line 8 in Algorithm 2). The proposed matching ordering method depends on the search information of ProSearch. In the experimental section, we will show that this method has outstanding performance compared to several other sorting methods.

3.4 Adaptive PropagationMethod

For the SIP, Glasgow [McCreesh et al., 2020] and PathLAD [Kotthoff et al., 2016] have outstanding performance, but they are completely different search strategies. Especially, Glasgow combines a weak propagation method with a fast restart mechanism. According to our preliminary experiments, Glasgow can make at least 104 recursive calls per second for some instances. On the contrary, PathLAD uses a strong propagation method at each stage of the search, and thus it sometimes makes less than one recursive call per second when dealing with some large target graphs. Based on our observations, no current algorithms for the SIP use different strengths of propagation methods at different stages of the search. Thus, our motivation is to design a method that can flexibly use some propagation methods in the search.

In the case of ProSearch procedure, the algorithm always uses LAD-filtering and GAC(allDiff). Both of them have high time complexity. Therefore, in the main search procedure, we design an adaptive propagation method to guide the use of strong propagation methods GAC(allDiff) and LAD-filtering. The pseudo-code of APM is shown in Algorithm 3.

Let us consider LAD-filtering first. We define a target graph to be sparse if the median of its vertex degrees is less than degm. In our work, degm is set to 20. When a pattern vertex vp is matched to a target vertex vt i ∈ Dom(vp ), LAD-filtering ensures that every pattern vertex in NGp (vp ) can match different target vertices in NGt (vt i ). Its execution time is based on the degree values of vp and vt i . Because degGt (vt i ) must be larger than or be equal to degGp (vp ), we just need to focus on the degree of target vertex vt i . If the target graph is sparse, the execution time of LAD-filtering is reasonable and we think that using it at every stage of the main search procedure is feasible.

In other cases, if the algorithm often backtracks due to lots of conflicts, the algorithm can actually turn to use some weak propagation methods including FC(Diff) and FC(Edges) instead of strong propagation methods. Although using weak propagation methods may result in searching deeper on the wrong branch compared to strong propagation methods, the algorithm can backtrack faster because the complexity of these weak methods is quite low. For such cases, calling LAD-filtering multiple times during the main search procedure will waste a lot of computation time. In our work, we analyze whether backtracking often occurs in the main search procedure by observing the values of nbnodes and nbfail. Meanwhile, we use a parameter max tries as the upper bound of nbnodes. In detail, on the one hand, If nbnodes is smaller than max tries, it may occur in the early stage of the search procedure. Because the backtracking for branch selection is costly, we want to explore more conflicts by us-

Algorithm 3 APM

Input: Target graph Gt and a non-matched pattern vertex vp Output: The reduced domain D(vp ) of vp 1: reduce D(vp ) based on FC(Diff) and FC(Edges); 2: if Gt is not a sparse graph then 3: if nbnodes > max tries && nbfail/nbnodes > β1 then 4: switchL := 0; 5: if switchL == 0 at the first time then 6: Nb := nbnodes; 7: end if 8: end if 9: if switchL == 0 && nbnodes > 2Nb && nbfail/nbnodes > β2 then switchA := 0; reduce D(vp ) based on LAD-filtering; reduce D(vp ) based on GAC(allDiff);

10: 11: end if 12: end if 13: if switchL == 1 then 14: 15: end if 16: if switchA == 1 then 17: 18: end if 19: return D(vp );

ing LAD-filtering on the wrong branches as early as possible. On the other hand, if nbnodes is larger than max tries, the algorithm has already explored some parts of the whole search space. For this case, we think the relationship between nbfail and nbnodes can provide some useful information for a given instance. If the number of failed vertices in the search procedure is large (i.e., nbfail/nbnodes > β1 where β1 is a parameter), it means that the algorithm has already backtracked a lot and thus the algorithm no longer uses LAD-filtering (Lines 3-4).

In the following, we consider another strong propagation method GAC(allDiff). Although this constraint has high time complexity in theory, it is actually faster than LAD-filtering in most cases. We will explain this reason as below. GAC(allDiff) constructs a bipartite graph between pattern vertices and target vertices. If a pattern vertex vp is matched to a target vertex vt i , GAC(allDiff) will remove vt i from the domain of some other pattern vertices and ensures that all pattern vertices can still match different target vertices. Removing a selected target vertex from the generated bipartite graph only needs to find the next free target vertices for some pattern vertices by looking for an augmenting path [Derigs, 1981]. In fact, the size of a given target graph is usually larger than that of a corresponding pattern graph. Thus, when the sizes of the two graphs are quite different, GAC(allDiff) is likely to be run in linear time.

In the main search procedure, after disabling the LAD-filtering, the algorithm begins to consider whether to disable the GAC(allDiff). When LAD-filtering is forbidden for the first time, we use variable Nb to record the current value of nbnodes (Lines 5-6). GAC(allDiff) will continue to be used until nbnodes has been increased twofold since LAD-filtering is disabled at the first time, i.e., nbnodes > 2Nb. At this time, if the algorithm still backtracks frequently (i.e., nbfail/nbnodes > β2 where β2 is a parameter), we disable GAC(allDiff) (Lines 9-10). In the subsequent search procedure, the algorithm doesn't employ any strong propagation methods including LAD-filtering and GAC(allDiff).

Remark that, in our work, the switch of propagation methods is one-way. The strength of weak propagation methods increases significantly with search depth, so there is no need to switch back to using strong propagation methods when the search depth reaches a certain point. Based on our preliminary experiments, we found that one-way switching was both straightforward and effective, whereas two-way switching exhibited poor performance on some instances. Recently, researchers have developed dynamic choice methods for several well-known problems, such as CSP [Stergiou, 2021]. One crucial step in algorithm design is to dynamically combine various methods for a particular problem. It's worth noting that our method is the first to use a dynamic choice approach to select propagation methods for the SIP.

Here, we will summarize the impact of the parameters β1 and β2 on the propagator choices as below. Parameters β1 and β2 are two thresholds that define whether a given instance is easy to backtrack due to numerous conflicts. When the conflict ratio is larger than these two parameters, we turn to using some simple propagation methods to make backtrack fast. Specifically, a higher value of β1 indicates a greater tolerance for conflicts, allowing us to use all propagation methods. However, if the conflict ratio surpasses β1, we discard the LAD-filtering method. On the other hand, a larger value of β2 implies a higher tolerance for conflicts to solely rely on the strong propagation method GAC(allDiff). When the conflict ratio exceeds β2, we also abandon the GAC(allDiff).

4 Experimental Evaluation

In this section, we carry out experiments to evaluate PathLAD+ on a broad range of various benchmarks, compared against the state-of-the-art algorithms for the SIP.

4.1 Benchmarks

For our experiments, we select all used instances from [Kraiczy and McCreesh, 2021; Liu et al., 2022] which can also download from the website 1 . In total, we choose 15396 instances, which can be grouped into 8 benchmarks.

• images-CVIU11 (6278 instances): This benchmark includes 43 pattern graphs and 146 target graphs, which have been generated from segmented images [Damiand et al., 2011]. In the benchmark, pattern graphs have between 22 and 151 vertices, whereas target graphs have between 1072 and 5972 vertices.
• meshes-CVIU11 (3018 instances): It is composed of 6 pattern graphs and 503 target graphs, which have been generated from meshes modeling 3D object [Damiand et al., 2011]. The number of vertices for pattern graphs is from 40 to 199, while the number of vertices is from 208 to 5873.

1 http://liris.cnrs.fr/csolnon/SIP.html

• images-PR15 (24 instances): There are 24 pattern graphs that have between 4 and 170 vertices and 1 target graph that has 4838 vertices. All the graphs have been derived from segmented images [Solnon et al., 2015].
• scalefree (100 instances): Each instance contains a target graph whose vertices are between 200 and 1000 and a pattern graph whose vertices are 90% of the vertices of the corresponding target graph. All the instances in the benchmark have been randomly generated using a power law distribution of degrees [Solnon, 2010].
• si (1170 instances): Each instance is composed of a target graph (between 200 and 1296 vertices) and a pattern graph (between 20% and 60% of the vertices of the corresponding target graph). This benchmark is from bounded valence graphs, modified bounded valence graphs, 4D meshes, and random generated graphs [Solnon, 2010].
• phase-transition (200 instances): These random instances are chosen to be close to the satisfiableunsatisfiable phase transition. Pattern graphs have 30 vertices, while target graphs have 150 vertices [McCreesh et al., 2016b].
• LV (1176 instances): The selected 49 graphs whose vertices are between 10 and 128 are considered as pattern and target graphs, and this benchmark has already been used as the tested benchmark [Liu et al., 2020]. These graphs have different properties [Solnon, 2010], such as connected, biconnected, triconnected, etc.
• LargerLV (3430 instances): From the above 49 LV graphs as the pattern graph and the other 70 graphs as the target graph whose vertices are between 138 and 6671. More details of the target graphs can be seen on the website 2 .

4.2 Experiment Setup

We compare PathLAD+ with four state-of-the-art SIP algorithms, including Glasgow+Clq [Kraiczy and McCreesh, 2021], PathLAD [Kotthoff et al., 2016], RI [Bonnici et al., 2013] and VF2 [Cordella et al., 2004]. The codes of these competitors are kindly provided by the authors. Our source code is publicly available at github 3 . Our proposed algorithm and four competitors are all implemented in C++ and compiled by g++ with '-O3' option. All the algorithms are run on Intel(R) Xeon(R) Platinum 8260 CPU @ 2.40GHz 512GB RAM under CentOS 7.9. The cutoff time is 3600 seconds for each instance. According to our preliminary experiments, parameters max tries, β1, and β2 are set to 1000, 0.85, and 0.8, respectively.

For each algorithm, we report the number of instances for each benchmark (#inst) and the number of successful solved instances (#solved). The bold values in the tables indicate the best solution among all the algorithms.

2 https://github.com/ciaranm/cpaior2021-finding-subgraphswith-side-constraints/tree/main/instances/largerGraphs

3 https://github.com/yiyuanwang1988/PathLAD-Plus

4.3 Experiment Results

We show the experiment results of our proposed algorithm and all competitors in Table 1. As observed from the results of Table 1, PathLAD+ performs much better than our baseline algorithm PathLAD on all the benchmarks. Overall, the performance of PathLAD+ totally dominates Glasgow+Clq, PathLAD, RI, and VF2. Because all algorithms can solve simple instances very well, we mainly focus on some hard instances. We can find that the performance of PathLAD+ is significantly better than all competitors on some hard benchmarks, especially in meshes-CVIU11. In this benchmark, all competitors have at least more than 30 unsolvable instances, whereas PathLAD+ only has 10 unsolvable instances within a cutoff time. Among the selected 15396 instances, PathLAD+ can solve 15197 instances within a cutoff time whereas the current best algorithm Glasgow+Clq can only solve 15141 instances. Furthermore, to intuitively display the performance of each algorithm, we report detailed results in Figure 1, which verifies the effectiveness of our proposed algorithm.

4.4 Analysis of Proposed Strategies

To confirm the effectiveness of our proposed matching ordering method, we evaluate different matching ordering methods on our baseline algorithm PathLAD, including 1) PathLAD-our uses our proposed matching ordering method; 2) PathLAD-degree selects a target vertex with the biggest degree value from the given domain; 3) PathLAD-random chooses a random target vertex from the given domain; 4) PathLAD-anti picks a target vertex with the smallest degree value from the given domain. Since different matching ordering methods only affect some isomorphic satisfiable instances [Archibald et al., 2019], we have shown the performance of different matching ordering methods in these instances in Figure 2. Results show that our proposed matching ordering method performs better than other methods. Moreover, the proposed sorting method effectively utilizes the useful information generated from the probing search procedure, and it clearly improves the performance of SIP.

We compare PathLAD with one alternative algorithm PathLAD-1 that uses the adaptive propagation method. PathLAD-1 and PathLAD don't use any matching ordering methods, and the effectiveness of the adaptive propagation method can be clearly observed in Figure 3. The different colored points show the instance from the different benchmarks. Figures 2 and 3 intuitively show that the proposed two strategies make an important role in our proposed algorithm. Besides, because PathLAD is a baseline algorithm of our proposed algorithm, we compare PathLAD+ with PathLAD in terms of run time in Figure 4. Once again, the results show the superiority of PathLAD+.

5 Conclusion

In this paper, we propose a probing search procedure, a novel matching ordering method, and an adaptive propagation method for the SIP. Based on the above strategies, we develop an efficient algorithm called PathLAD+. Experiments show PathLAD+ significantly outperforms the state-of-theart SIP algorithms.

As for future work, the proposed adaptive propagation method can be considered as a general idea to solve some other NP-hard problems [Chen et al., 2023].

Acknowledgements

This work was supported by CAS Project for Young Scientists in Basic Research (Grant No.YSBR-040), NSFC (61806050), Jilin Science and Technology Association QT202005, and Science and Technology Development Program of Jilin Provice (YDZJ202201ZYTS412 and 20230101060JC). We would like to thank the anonymous referees for their helpful comments.