1 Introduction
Many realworld domains incorporate large and complex networks of interconnected units. Examples include social networks, the Internet, or biological and chemical systems. These networks raise interesting questions regarding their structure. One of those questions asks whether a given network contains a particular pattern, which typically represents a specific behaviour of interest [1, 7, 15]. The problem of locating a particular pattern in the given network can be restated as a problem of locating a subgraph isomorphic to the given pattern graph in the network graph.
Formally, the Subgraph Isomorphism (SubIso) problem is, given two undirected graphs and , to decide whether there is a (not necessarily induced) subgraph of isomorphic to . Or, in other words, whether there is an adjacencypreserving injective mapping from vertices of to vertices of . Since we do not require the subgraph to be induced (or the mapping to preserve nonadjacencies), some authors call this variant of the problem Subgraph Monomorphism.
For many applications it is not enough to just learn that the pattern does occur in the network, but it is necessary to actually obtain the location of an occurrence of the pattern or rather of all occurrences of the pattern [21, 27]. Because of that, we aim to solve the problem of subgraph enumeration, in which it is required to output all subgraphs of the network graph isomorphic to the pattern graph. In Subgraph Enumeration (SubEnum), given again two graphs and , the goal is to enumerate all subgraphs of isomorphic to . Note, that SubEnum is at least as hard as SubIso. We call the variants, where the problem is required to be induced IndSubIso and IndSubEnum, respectively.
As Clique, one of the problems on the Karp’s original list of 21 NPcomplete problems [18], is a special case of SubIso, the problem is NPcomplete. Nevertheless, there are many heuristic algorithms for SubEnum, many of them based on ideas from constraint programming (see Section 1.1), which give results in reasonable time for most instances. However, for each of them there are rather small instances which they find genuinely hard, as pointed out by McCreesh et al. [25]. Therefore, developing an alternative approach that could possibly cope with these hard instances would be of interest.
In this paper we focus on the well known randomized color coding approach [2], which presumably has almost optimal worst case time complexity. Indeed, its time complexity is with memory requirements of , where and denote the number of vertices in the network graph and the pattern graph , respectively, and is the treewidth of graph —a measure of treelikeness (see Section 1.2 for exact definitions). Moreover, we presumably cannot avoid the factor exponential in treewidth in the worst case running time, as Marx [23] presented an ETH^{1}^{1}1Exponential Time Hypothesis [17]based lower bound for Partitioned Subgraph Isomorphism problem.
Proposition 1 (Marx [23])
If there is a recursively enumerable class of graphs with unbounded treewidth, an algorithm , and an arbitrary function such that correctly decides every instance of Partitioned Subgraph Isomorphism with the smaller graph in in time , then ETH fails.
As the memory requirements of the color coding approach grow exponentially with treewidth of the pattern graph, existing implementations for subgraph enumeration based on this principle restrict themselves to paths [16] or trees [27], both having treewidth 1. As the real world applications might employ networks of possibly tens to hundreds of thousands of vertices and also pattern graphs with structure more complicated than trees, we need to significantly reduce the memory usage of the algorithm.
Using the principle of inclusionexclusion, Amini et al. [3, Theorem 15] suggested a modification of the color coding algorithm, which can decide whether the pattern occurs in the graph in expected time with memory requirements reduced to .^{2}^{2}2While the formulation of Theorem 15 in [3] might suggest that the algorithm actually outputs a witnessing occurrence, the algorithm merely decides whether the number of occurrences is nonzero (see the proof of the theorem). While single witnessing occurrence can be found by means of selfreduction (which is complicated in case of randomized algorithm), the inclusionexclusion nature of the algorithm does not allow to find all occurrences of pattern in the graph, which is our main goal.
Therefore, our approach rather follows the paradigm of generating only those parts of a dynamic programming table that correspond to subproblems with a positive answer, recently called “positive instance driven” approach [30]. This further prohibits the use of the inclusionexclusion approach of Amini et al. [3], since the inclusionexclusion approach tends to use most of the table and the term is itself prohibitive in the memory requirements for .
Because of the time and memory requirements of the algorithm, for practical purposes we restrict ourselves to pattern graphs with at most vertices.
Altogether, our main contribution is twofold:

We provide a practical implementation of the color coding algorithm of Alon, Yuster, and Zwick [2] capable of processing large networks and (possibly disconnected) pattern graphs of small, yet not a priory bounded, treewidth.

We supply a routine to extract the occurrences of the subgraphs found from a run of the algorithm.
It is important to note that all the modifications only improve the practical memory requirements and running time. The theoretical worst case time and space complexity remain the same as for the original color coding algorithm and the algorithm achieves these, e.g., if the network graph is complete. Also, in such a case, there are occurrences of the pattern graph in the network implying a lower bound on the running time of the enumeration part.
In Section 2 we describe our modifications to the algorithm and necessary tools used in the process. Then, in Section 5, we benchmark our algorithm on synthetic and realistic data and compare its performance with available existing implementations of algorithms for subgraph isomorphism and discuss the results obtained. Section 6 presents future research directions.
1.1 Related work
There are several algorithms tackling SubIso and its related variants. Some of them only solve the variant of subgraph counting, our main focus is however on algorithms actually solving SubEnum. Following Carletti et al. [10] and Kimmig et al. [19], we categorize the algorithms by the approach they use (see also Kotthoff et al. [20] for more detailed description of the algorithms). Many of the approaches can be used both for induced and noninduced variants of the problem, while some algorithms are applicable only for one of them.
Vast majority of known algorithms for the subgraph enumeration problem is based on the approach of representing the problem as a searching process. Usually, the state space is modelled as a tree and its nodes represent a state of a partial mapping. Finding a solution then typically resorts to the usage of DFS in order to find a path of mappings in the state space tree which is compliant with isomorphism requirements. The efficiency of those algorithms is largely based on early pruning of unprofitable paths in the state space. Indeed, McCreesh et al. [25] even measure the efficiency in the number of generated search tree nodes. The most prominent algorithms based on this idea are Ullmann’s algorithm [31], VF algorithm and its variants [8, 10, 12, 13] (the latest VF3 [8] only applies to IndSubEnum) and RI algorithm [7]. The differences between these algorithms are based both on employed pruning strategies and on the order in which the vertices of pattern graph are processed (i.e. in the shape of the state space tree).
Another approach is based on constraint programming, in which the problem is modelled as a set of variables (with respective domains) and constraints restricting simultaneous variable assignments. The solution is an assignment of values to variables in a way such that no constraint remains unsatisfied. In subgraph isomorphism, variables represent pattern graph vertices, their domain consists of target graph vertices to which they may be mapped and constraints ensure that the properties of isomorphism remain satisfied. Also in this approach, a state space of assignments is represented by a search tree, in which nonprofitable branches are to be filtered. Typical algorithms in this category are LAD algorithm [28], Ullmann’s bitvector algorithm [32], and Glasgow algorithm [24]. These algorithms differ in the constraints they use, the way they propagate constraints, and in the way they filter state space tree.
There are already some implementations based on the color coding paradigm, where the idea is to randomly color the input graph and search only for its subgraphs, isomorphic to the pattern graph, that are colored in distinct colors (see Section 2.1 for more detailed description). This approach is used in subgraph counting algorithms, e.g., in ParSE [33], FASCIA [26], and in [1], or in algorithms for path enumeration described in [27] or in [16]. Each of these algorithms, after the color coding step, tries to exploit the benefits offered by this technique in its own way; although usually a dynamic programming sees its use. Counting algorithms as ParSE and FASCIA make use of specifically partitioned pattern graphs, which allow to use combinatorial computation. Weighted path enumeration algorithms [27, 16] describe a dynamic programming approach and try to optimize it in various ways. However, to the best of our knowledge there is no color coding algorithm capable of enumerating patterns of treewidth larger than 1.
Our aim is to make step towards competitive implementation of color coding based algorithm for SubEnum, in order to see, where this approach can be potentially beneficial against the existing algorithms. To this end, we extend the comparisons of SubEnum algorithms [9, 20, 25] to color coding based algorithms, including the one proposed in this paper.
1.2 Basic definitions
All graphs in this paper are undirected and simple. For a graph we denote its vertex set, the size of this set, its edge set, and the size of its edge set.
As already said, we use the color coding algorithm. The algorithm is based on a dynamic programming on a nice tree decomposition of the pattern graph. We first define a tree decomposition and then its nice counterpart.
Definition 1
A tree decomposition of a graph is a triple , where is a tree rooted at node and is a mapping satisfying: (i) ; (ii) , such that ; (iii) the nodes form a connected subtree of .
We shall denote bag as . The width of tree decomposition is . Treewidth of graph is the minimal width of a tree decomposition of over all such decompositions.
Definition 2
A tree decomposition of a graph is nice if , , and each node is of one of the following four types:

Leaf node— has no children and ;

Introduce node— has exactly one child and for some ;

Forget node— has exactly one child and for some ;

Join node— has exactly two children and .
Note that for practical purposes, we use a slightly modified definition of nice tree decomposition in this paper. As the algorithm starts the computation in a leaf node, using the standard definition with empty bags of leaves [14] would imply that the tables for leaves would be somewhat meaningless and redundant. Therefore, we make bags of leaf nodes contain a single vertex.
Definition 3
For a tree decomposition , we denote by the set of vertices in and in for all descendants of in . Formally .
Note that, by Definition 3, for the root of we have and .
2 Algorithm Description
In this section we first briefly describe the idea of the original color coding algorithm [2], show, how to alter the computation in order to reduce its time and memory requirements, and describe implementation details and further optimizations of the algorithm.
2.1 Idea of the Algorithm
The critical idea of color coding is to reduce the problem to its colorful version. For a graph and a pattern graph , we color the vertices of with exactly colors. We use the randomized version, i.e., we create a random coloring . After the coloring, the algorithm considers as valid only subgraphs of that are colorful copies of as follows.
Definition 4
Subgraph of a graph is a colorful copy of with respect to coloring , if is isomorphic to and all of its vertices are colored by distinct colors in .
As the output of the algorithm heavily depends on the chosen random coloring of
, in order to reach some predefined success rate of the algorithm, we need to repeat the process of coloring several times. The probability of a particular occurrence of pattern graph
becoming colorful with respect to the random coloring is , which tends to for large . Therefore, by running the algorithm times, each time with a random coloring , the probability that an existing occurrence of the pattern will be revealed in none of the runs is at most . While using more colors can reduce the number of iterations needed, it also significantly increases the memory requirements. Hence, we stick to colors. Even though it is possible to derandomize such algorithms, e.g., by the approach shown in [14], in practice the randomized approach usually yields the results much quicker, as discussed in [27]. Moreover, we are not aware of any actual implementation of the derandomization methods.The main computational part of the algorithm is a dynamic programming. The target is to create a graph isomorphism . We do so by traversing the nice tree decomposition of the pattern graph and at each node of the tree decomposition, we construct possible partial mappings with regard to required colorfulness of the copy. Combination of partial mappings consistent in colorings then forms a desired resulting mapping.
The semantics of the dynamic programming table is as follows. For any tree decomposition node , any partial mapping and any color subset , we define if there is an isomorphism of to a subgraph of such that:

for all ;

is a colorful copy of using exactly the colors in , that is, and is injective on .
If there is no such isomorphism, then we let . We denote all configurations for which as nonzero configurations.
The original version of the algorithm is based on topdown dynamic programming approach with memoization of already computed results. That immediately implies a big disadvantage of this approach—it requires the underlying dynamic programming table (which is used for memoization) to be fully available throughout the whole run of the algorithm. To avoid this inefficiency in our modification we aim to store only nonzero configurations, similarly to the recent “positive instance driven” dynamic programming approach [30].
3 Obtaining a Nice Tree Decomposition
The algorithm requires a nice tree decomposition of the pattern graph for its work. The running time of the algorithm actually does not depend on the treewidth of the graph, but rather on the width of the tree decomposition supplied. On one hand, for a given graph and an integer , the problem of determining whether the treewidth of is at most , is NPcomplete [4]. On the other hand, as the main algorithm is exponential in the size of the pattern anyway, we can afford to use exponential time algorithms in order to obtain a tree decomposition of minimum possible width. We employ known technique based on graph triangulation, chordal graphs and elimination orderings, combining the ideas of Bodlaender et al. [5, 6]. Time complexity of this algorithm is .
After obtaining tree decomposition of the pattern graph, we construct a nice tree decomposition as described by Cygan et al. [14]—any tree decomposition of a graph that consists of nodes with width at most , can be, in time, transformed to a nice tree decomposition of with nodes and width bounded by . Our further modification to the nice tree decomposition does not affect the running time.
3.1 Initial Algorithm Modification
In our implementation, we aim to store only nonzero configurations, therefore we need to be able to construct nonzero configurations of a parent node just from the list of nonzero configurations in its child/children.
We divide the dynamic programming table into lists of nonzero configurations, where each nice tree decomposition node has a list of its own. Formally, for every node , let us denote by a list of all mappings with a list of their corresponding color sets , for which . The list for all is, in terms of contained information, equivalent to maintaining the whole table —all configurations not present in the lists can be considered as configurations with a result equal to zero.
Dynamic Programming Description
We now describe how to compute the lists for each type of a nice tree decomposition node.
For a leaf node , there is only a single vertex in to consider. We can thus map to all possible vertices of , and we obtain a list with partial mappings , in which the color list for each mapping contains a single color set .
For an introduce node and its child in , we denote by the vertex being introduced in , i.e., . For all nonzero combinations of a partial mapping and a color set in the list , we try to extend by all possible mappings of the vertex to the vertices of . We denote one such a mapping as . We can consider mapping as correct, if (i) the new mapping of the vertex extends the previous colorset , that is, , and (ii) is edge consistent, that is, for all edges between currently mapped vertices, i.e., in our case , there must be an edge . However, because was by construction already edge consistent, it suffices to check the edge consistency only for all edges in with as one of their endpoints, i.e., for all edges with . After checking those two conditions, we can add to .
For a forget node and its child in , there is not much work to do, as in the bottomup construction, we directly obtain the parent list of nonzero configurations. We denote by , the vertex being forgotten in , i.e., . In this case, for all partial mappings in the list , we create a new mapping that excludes the mapping for vertex , i.e., . Color sets corresponding to are carried over to , as they represent colors already used in the construction. After this step, we might need to merge color lists of previously different mappings, as after the removal of the mapping , they might become the same mappings.
For a join node , we denote by and its children in . We traverse the children lists and and look for partial mappings and , for which holds. Such mappings are the only ones to potentially form a new nonzero configuration in the parent list, as due to the fact that , we construct the new partial mapping as . However, for each such mapping , we must also construct the new list of color sets, which would afterwards be corresponding to the mapping in the parent list. We do that by traversing color lists corresponding to mappings and in and , respectively, and for particular sets and from the color lists of and construct a new color set . We must check, whether the intersection of and contains exactly the colors to color the vertices in . That is, for mapping , we add to a color set , if .
Because we build the result from the leaves of the nice tree decomposition, we employ a recursive procedure on its root, in which we perform the computations in a way of a postorder traversal of a tree. From each visited node, we obtain a bottomup dynamic programming list of nonzero configurations. After the whole nice tree decomposition is traversed, we obtain a list of configurations, that were valid in its root. Such configurations thus represent solutions found during the algorithm, from which we afterwards reconstruct results. Note that as we prepend a root with no vertices in its bag to the nice tree decomposition, there is a nonzero number of solutions if and only if, at the end of the algorithm, the list contains a single empty mapping using all colors.
3.2 Further Implementation Optimizations
Representation of Mappings
For mapping representation, we suppose that the content of all bags of the nice tree decomposition stays in the same order during the whole algorithm. This natural and easily satisfied condition allows us to represent a mapping in a nice tree decomposition node simply by an ordered tuple of vertices from . From this, we can easily determine which vertex from is mapped to which vertex in . Also, for a mapping in an introduce or a forget node, we can describe a position in the mapping, on which the process of introducing/forgetting takes place.
Representation of Color Sets
We represent color sets as bitmasks, where the th bit states whether color is contained in the set or not. For optimization purposes, we represent bitmasks with an integer number. As we use colors in the algorithm and restricted ourselves to pattern graphs with at most vertices, we represent a color set with a bit number.
Compressing the Lists
Because we process the dynamic programming lists one mapping at a time, we store these lists in a compressed way and decompress them only on a mapping retrieval basis. We serialize dynamic programming lists into a simple buffer of bytes. We store a dynamic programming list as a continuous group of records, each of which represents one partial mapping and its corresponding list of color sets. Each record contains:

a mapping in the form of ordered tuple of vertices,

the number of color sets included,

and color sets corresponding to the mapping in the form of nonnegative integer numbers.
While deserializing, the number of vertices in a mapping can be easily determined by the size of the bag of a particular node.
As there is no requirement on the order of the color sets, we store these sets sorted in the increasing order of their number representation. This allows us to effectively use delta compression. Moreover, we use variable length encoding to store the numbers into buffer.
For several routines we use the LibUCW library,^{3}^{3}3LibUCW is downloadable from http://www.ucw.cz/libucw/. a C language library highly optimized for performance. The employed routines include data structures like growing buffers, hash tables, redblack trees, fast sorters, fast buffered input/output, and also the efficient variable length encoding of integers.
Masking Unprofitable Mappings
Our implementation supports an extended format of input graphs where one can specify for each vertex of the network, which vertices of the pattern can be mapped to it. This immediately yields a simple degreebased optimization. Before the run of the main algorithm, we perform a linear time preprocessing of input graphs and only allow a vertex to be mapped to a vertex if .
Mapping Expansion Optimizations
The main “bruteforce” work of the algorithm is performed in two types of nodes—leaf and introduce nodes, as we need to try all possible mappings of a particular vertex in a leaf node or all possible mappings of an introduced vertex in a introduce node to a vertex from . We describe ways to optimize the work in introduce nodes in this paragraph.
Let be an introduce node, the vertex introduced and a mapping from a nonzero configuration for the child of . We always need to check whether the new mapping of is edge consistent with the mapping of the remaining vertices for the corresponding bag, i.e., whether all edges of incident on would be realized by an edge in . Therefore, if has any neighbors in , then a vertex of is a candidate for the mapping of only if it is a neighbor of all vertices in the set , i.e., the vertices of , where the neighbors of in are mapped. Hence, we limit the number of candidates by using the adjacency lists of the already mapped vertices.
In the case we have to use different approach. The pattern graphs tend to be smaller than the input graphs by several orders of magnitude. Hence, if the introduced vertex is in the same connected component of as some vertex already present in the bag, a partial mapping processed in an introduce node anchors the possible resulting component to a certain position in . During the construction of possible mapping of , it is useless to try mapping to vertices in that are very distant to the current position and, therefore, could by no means form a resulting subgraph isomorphic to the component of . We obtain the maximal possible distance to be considered in as a minimal distance of to a vertex in . I.e., we determine . Due to the limit on pattern graph size, we precompute distances between every pair of vertices by using BFS on each vertex before the start of the algorithm. Then it suffices to try vertices from that are in distance at most from in ; again, by a simple BFS usage.
Only if there is no vertex in the bag sharing a connected component of with , we have to fall back to trying all possible mappings.
4 Result Reconstruction
It is usual in dynamic programming to store a witness or all witnesses for each reasonable value in the table. However, this would completely neglect the effect of compression of the lists and significantly increase the memory requirements.
Hence, to enumerate all results, we recursively traverse the computed dynamic programming lists starting with the root of underlying nice tree decomposition. The main idea is to ask the child (or children ) of a node, how was a certain partial mapping obtained during the computation. In other words, for each partial mapping in of interest, we will extract partial mappings in (or possibly also partial mappings in ) that lead to the addition of to . To preserve the colorfulness (and thus validity) of glued partial mappings, with a partial mapping we also recursively pass a set of colors that can still be used in the choice of corresponding partial mappings to glue with . By applying this approach recursively and by gluing the partial information together in mappings, we obtain all possible ways an empty mapping in could have been obtained—which is exactly the same as all possible ways a pattern graph can be mapped to the original graph. For the efficiency of the computation, we process each node only once and we thus recursively pass all partial mappings of interest and their remaining colorsets at once.
Formally, we will consider a recursive call of function , where , is a set of pairs , is a partial mapping and is a color subset . This function returns a list containing, for each pair in , a list of mappings that lead to the appearance of in .
In detail, the function proceeds as follows. If node has a child , we first scan for pairs , in which and and colorsets and are in some sense consistent. Pairs satisfying these conditions are then recursively passed to in a list , from which we obtain a list . In case of a join node, we perform this process for both children. After that, we construct the resulting list from recursively computed list(s) by an addition of a mapping of the introduced vertex (in the case of a introduce node) or by all possible combinations of mappings (in the case of a join node). In the end, we sometimes need to flatten the resulting list, so all resulting partial mappings corresponding to a pair are in a single list . We now describe the process and consistency conditions formally for each of the node types.
For a leaf node , we simply return a list of single element lists for each pair . Each such list contains a partial mapping .
For an introduce node and its child in , we denote by the vertex being introduced in , i.e., . We create consisting of by setting, for each pair , and . Note that must be in , as otherwise would not be present in (or equivalently in ). Then we obtain . Let be constructed from and be the corresponding list of . We construct the list corresponding to by extending each by .
For a forget node and its child in , we denote by the vertex being forgotten in , i.e., . We add to for each all from such that and . Then we obtain . We construct from as follows. We do not need to modify any in any , as the computation in forget node adds no new information to the mapping. However, as there may have been multiple (say ) pairs constructed from a single , we flatten all lists obtained from such pairs into a single list corresponding to .
For a join node , we denote by and its children in . In this case, we create and consisting of and , respectively. For each , we find all and , such that , and and add them to and , respectively. Then we obtain and . Let and , …, and be the pairs constructed from single and let and be the corresponding lists in and , respectively. We obtain corresponding to as a union of lists , where contains the union of and for each and each .
5 Experimental Results
The testing was performed on a 64bit linux system with Intel Xeon CPU E31245v6@3.70GHz and 32GB 1333MHz DDR3 SDRAM memory.
The module was compiled with gcc
compiler (version 7.3.1)
with O3
optimizations enabled.
Implementation and instances utilized in the testing are available at http://users.fit.cvut.cz/malikjo1/subiso/.
All results are an average of 5 independent measurements.
We evaluated our implementation in several ways. Firstly, we compare available implementations on two different real world source graphs and a set of moreorless standard target graph patterns. Secondly, we compare available implementations on instances from ICPR2014 Contest on Graph Matching Algorithms for Pattern Search in Biological Databases [11] with suitably small patterns. We also adapt the idea of testing the algorithms on ErdősRényi random graphs [25].
5.1 Algorithm Properties and Performance
In the first two subsection we used two different graphs of various properties as target graph . The first instance, Images, is built from an segmented image, and is a courtesy of [29]. It consists of vertices and edges. The second instance, Trans, is a graph of transfers on bank accounts. It is a very sparse network, which consists of vertices and undirected edges.
For the pattern graphs, we first use a standard set of basic graph patterns, as the treewidth of such graphs is well known and allows a clear interpretation of the results. In particular, we use paths, stars, cycles, an complete graphs on vertices, denoted , , , and with treewidth 1, 1, 2, and , respectively. We further used grids on vertices, with treewidth . Secondly, we use a special set of pattern graphs in order to demonstrate performance on various patterns. Patterns , , , and have , , , and vertices, , , , edges, and treewidth , , , and , respectively. Patterns , , and appear in both dataset, pattern in neither and pattern is disconnected. Description of these pattern graphs is shown in Table 1 and their illustration is shown in Figure 5.
Pattern  Description  

path on vertices  
star on vertices  
cycle on vertices  
complete graph on vertices  
grid graph on vertices  
Pattern  subgraph of both datasets with vertices, edges  
Pattern  subgraph of both datasets with vertices, edges  
Pattern  graph with vertices, edges  
Pattern  disconnected graph with vertices, edges 
Due to randomization, in order to achieve some preselected constant error rate, we need to repeat the computation more than once. The number of found results thus depends not only on the quality of the algorithm, but also on the choice of the number of its repetitions. Hence, it is logical to measure performance of the single run of the algorithm. Results from such a testing, however, should be still taken as a rough average, because the running time of a single run of the algorithm depends on many factors.
Therefore, we first present measurements, where we average the results of many single runs of the algorithm. We average not only the time and space needed, but also the number of found subgraphs. To obtain the expected time needed to run the whole algorithm, it suffices to sum the time needed to create a nice tree decomposition and times the time required for a single run, if there are runs in total.
Pattern  Comp. time [ms]  Comp. memory [MB]  Occurrences 

240  12.73  3488.21  
160  8.52  732.46  
90  10.54  76.18  
4  5.37  114.72  
20  7.24  239.17  
70  9.34  26.64  
5  6.46  0  
90  13.42  0  
Pattern  80  9.14  292.48 
Pattern  10  7.17  6.85 
Pattern  10  5.30  0 
Pattern  40  10.14  426.76 
Pattern  Comp. time [s]  Comp. memory [MB]  Occurrences 

6.15  33.82  54572.94  
0.32  34.40  562.54  
0.17  34.83  91.49  
0.11  30.11  15.32  
0.24  33.84  0  
0.09  27.91  0.72  
2.16  48.68  0  
Pattern  54.42  109.88  94145.82 
Pattern  0.18  34.12  1127.11 
Pattern  0.29  35.16  0 
Pattern  0.25  47.54  178.29 
5.2 Comparison on Real World Graphs and Fixed Graph Patterns
We compare our implementation to three other tools for subgraph enumeration: RI algorithm [7] (as implemented in [22]), LAD algorithm [28] and color coding algorithm for weighted path enumeration [16] (by setting, for comparison purposes, all weights of edges to be equal). The comparison is done on the instances from previous subsection and only on pattern graphs which occur at least once in a particular target graph.
In comparison, note the following specifics of measured algorithms. The RI algorithm does not support outputting solutions, which might positively affect its performance. LAD algorithm uses adjacency matrix to store input graphs, and thus yields potentially limited use for graphs of larger scale. Neither of RI or LAD algorithms supports enumeration of disconnected patterns.^{4}^{4}4 When dealing with disconnected patterns, one could find the components of the pattern one by one, omitting the vertices of the host graph used by the previous component. However, this would basically raise the running time of the algorithm to the power equal to the number of components of the pattern graph. Also we did not measure the running time of the weighted path algorithm on nonpath queries and also on Trans dataset, as its implementation is limited to graph sizes of at most .
We run our algorithm repeatedly to achieve an error rate of . In order to be able to measure the computation for larger networks with many occurrences of the pattern, we measure only the time required to retrieve no more than first solutions and we also consider running time greater than 10 minutes (600 seconds) as a timeout. Since we study noninduced occurrences (and due to automorphisms) there might be several ways to map the pattern to the same set of vertices. Other measured algorithms do count all of them. Our algorithm can behave also like this, or can be switched to count only occurrences that differ in vertex sets. For the sake of equal measurement, we use the former version of our algorithm.
Pattern  Our algorithm  RI algorithm  LAD algorithm  Weighted path 

31.12  0.11  28.86  362.41  
53.17  1.25  13.63  
104.30  3.7  8.18  
0.94  0.07  0.43  –  
4.98  0.14  35.18  –  
151.25  3.44  174.27  –  
Pattern  43.11  0.82  36.60  – 
Pattern  91.93  0.41  0.83  – 
Pattern  23.54  –  –  – 
Pattern  Our algorithm  RI algorithm  LAD algorithm  Weighted path 

11.64  2.53  59.57  –  
44.72  4.77  34.00  –  
295.11  24.11  28.98  –  
19.25  0.56  24.58  –  
4.38  0.70  2.24  –  
Pattern  61.36  11.85  52.77  – 
Pattern  23.91  1.63  31.67  – 
Pattern  481.25  –  –  – 
From Table 4 and Table 5, we can see that RI algorithm outperforms all other measured algorithms. We can also say our algorithm is on par with LAD algorithm, as the results of comparison of running times are similar, but vary instance from instance. Our algorithm nevertheless clearly outperforms another color coding algorithm, which on one hand solves more complicated problem of weighted paths, but on the another, is still limited only to paths. Also, our algorithm is the only algorithm capable of enumerating disconnected patterns.
The weak point of the color coding approach (or possibly only of our implementation) appears to be the search for a pattern of larger size with very few (or possibly zero) occurrences. To achieve the desired error rate, we need to repeatedly run the algorithm many times. Therefore our algorithm takes longer time to run on some instances (especially close to zerooccurrence ones), which are easily solved by the other algorithms.
5.3 ICPR2014 Contest Graphs
To fully benchmark our algorithm without limitations on time or number of occurrences found, we perform a test on ICPR2014 Contest on Graph Matching Algorithms for Pattern Search in Biological Databases [11].
In particular, we focus our attention on a Molecules dataset, containing 10,000 (target) graphs representing the chemical structures of different small organic compounds and on a Proteins dataset, which contains 300 (target) graphs representing the chemical structures of proteins and protein backbones. Target graphs in both datasets are sparse and up to 99 vertices or up 10,081 vertices for Molecules and Proteins, respectively.
In order to benchmark our algorithm without limiting its number of iterations, we focus on pattern graphs of small sizes, which offer reasonable number of iterations for an error rate of . Both datasets contain 10 patterns for each of considered sizes constructed by randomly choosing connected subgraphs of the target graphs. We obtained an average matching time of all pattern graphs of a given size to all target graphs in a particular dataset.
Targets  Pattern size  Our algorithm  LAD algorithm  RI algorithm 

Molecules  4  0.01  0.01  0.01 
Molecules  8  0.67  0.14  0.01 
Proteins  8  19.45  8.83  0.51 
From the results in Table 6, we can see our algorithm being on par with LAD algorithm, while being outperformed by RI algorithm. However, we mainly include these results as a proof of versatility of our algorithm. As discussed in [25], benchmarks created by constructing subgraphs of target graphs do not necessarily encompass the possible hardness of some instances and might even present a distorted view on algorithms’ general performance. Thus, in the following benchmark we opt to theoretically analyze our algorithm.
5.4 ErdősRényi Graph Setup
In order to precisely analyze the strong and weak points of our algorithm we measure its performance is a setting where both the pattern and the target are taken as an ErdősRényi random graph of fixed size with varying edge density and compare the performance of our algorithm with the analysis of McCreesh et al. [25], which focused on algorithms Glasgow, LAD, and VF2.
An ErdősRényi graph is a random graph on vertices where each edge is included in the graph independently at random with probability
. We measure the performance on target graph of 150 vertices and pattern graph of 10 vertices with variable edge probabilities. As our algorithm cannot be classified in terms of search nodes used (as in
[25]), we measure the time needed to complete 10 iterations of our algorithm.From Fig. 8
we can see our algorithm indeed follows a well observed phase transition (transition between instances without occurrence of the pattern and with many occurrences of the pattern). If we compare our results from Fig.
6 to the results of [25], we can see that hard instances for our algorithm start to occur later (in terms of edge probabilities). However, due to the almost linear dependency of treewidth on edge probabilities (see Fig. 7), hard instances for our algorithm concentrate in the “upper right corner” of the diagram, which contains dense graphs with naturally large treewidth.6 Conclusion
We described an efficient implementation of the well known color coding algorithm for the subgraph isomorphism problem. Our implementation is the first colorcoding based algorithm capable of enumerating all occurrences of patterns of treewidth larger than one. Moreover, we have shown that our implementation is competitive with existing stateoftheart solutions in the setting of locating small pattern graphs. As it exhibits significantly different behaviour than other solutions, it can be an interesting contribution to the portfolio of known algorithms [20, 25].
As an obvious next step, the algorithm could be made to run in parallel. We also wonder whether the algorithm could be significantly optimized even further, possibly using some of the approaches based on constraint programming.
References
 [1] Alon, N., Dao, P., Hajirasouliha, I., Hormozdiari, F., Sahinalp, S.: Biomolecular network motif counting and discovery by color coding. Bioinformatics 24, 241–249 (2008)
 [2] Alon, N., Yuster, R., Zwick, U.: Colorcoding. J. ACM 42(4), 844–856 (1995)
 [3] Amini, O., Fomin, F.V., Saurabh, S.: Counting subgraphs via homomorphisms. SIAM J. Discrete Math. 26(2), 695–717 (2012)
 [4] Arnborg, S., Corneil, D.G., Proskurowski, A.: Complexity of finding embeddings in a ktree. SIAM J. Algebraic Discrete Methods 8(2), 277–284 (1987)
 [5] Bodlaender, H.L., Fomin, F.V., Koster, A.M.C.A., Kratsch, D., Thilikos, D.M.: On exact algorithms for treewidth. ACM Trans. Algorithms 9(1), 12:1–12:23 (2012)
 [6] Bodlaender, H.L., Koster, A.M.C.A.: Treewidth computations I. Upper bounds. Inf. Comput. 208(3), 259–275 (2010)
 [7] Bonnici, V., Giugno, R., Pulvirenti, A., Shasha, D., Ferro, A.: A subgraph isomorphism algorithm and its application to biochemical data. BMC Bioinformatics 14, 1–13 (2013)
 [8] Carletti, V., Foggia, P., Saggese, A., Vento, M.: Introducing VF3: A new algorithm for subgraph isomorphism. In: GbRPR 2017. LNCS, vol. 10310, pp. 128–139 (2017)
 [9] Carletti, V., Foggia, P., Vento, M.: Performance comparison of five exact graph matching algorithms on biological databases. In: ICIAP 2013. LNCS, vol. 8158, pp. 409–417. Springer (2013)
 [10] Carletti, V., Foggia, P., Vento, M.: VF2 plus: An improved version of VF2 for biological graphs. In: GbRPR 2015. LNCS, vol. 9069, pp. 168–177. Springer (2015)
 [11] Carletti, V., Foggia, P., Vento, M., Jiang, X.: Report on the first contest on graph matching algorithms for pattern search in biological databases. In: GbRPR 2015. LNCS, vol. 9069, pp. 178–187. Springer (2015)
 [12] Cordella, L.P., Foggia, P., Sansone, C., Vento, M.: Performance evaluation of the VF graph matching algorithm. In: 10th International Conference on Image Analysis and Processing, ICIAP 1999. pp. 1172–1177. IEEE Computer Society (1999)
 [13] Cordella, L.P., Foggia, P., Sansone, C., Vento, M.: A (sub)graph isomorphism algorithm for matching large graphs. IEEE Trans. Pattern Anal. Mach. Intell. 26(10), 1367–1372 (2004)
 [14] Cygan, M., Fomin, F.V., Kowalik, L., Lokshtanov, D., Marx, D., Pilipczuk, M., Pilipczuk, M., Saurabh, S.: Parameterized Algorithms. Springer (2015)

[15]
Dahm, N., Bunke, H., Caelli, T., Gao, Y.: Efficient subgraph matching using topological node feature constraints. Pattern Recognition
48(2), 317 – 330 (2015)  [16] Hüffner, F., Wernicke, S., Zichner, T.: Algorithm engineering for colorcoding with applications to signaling pathway detection. Algorithmica 52(2), 114–132 (2008)
 [17] Impagliazzo, R., Paturi, R.: On the complexity of kSAT. J. Comput. Syst. Sci. 62(2), 367–375 (2001)
 [18] Karp, R.M.: Reducibility among combinatorial problems. In: Symposium on the Complexity of Computer Computations, COCO 1972. pp. 85–103. The IBM Research Symposia Series, Plenum Press, New York (1972)
 [19] Kimmig, R., Meyerhenke, H., Strash, D.: Shared memory parallel subgraph enumeration. In: 2017 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW). pp. 519–529. IEEE Computer Society (2017)
 [20] Kotthoff, L., McCreesh, C., Solnon, C.: Portfolios of subgraph isomorphism algorithms. In: LION 10. LNCS, vol. 10079, pp. 107–122. Springer (2016)
 [21] Kuramochi, M., Karypis, G.: Frequent subgraph discovery. In: 2001 IEEE International Conference on Data Mining. pp. 313–320. IEEE Computer Society (2001)
 [22] Leskovec, J., Sosič, R.: Snap: A generalpurpose network analysis and graphmining library. ACM Trans. on Intel. Systems and Technology (TIST) 8(1), 1 (2016)

[23]
Marx, D.: Can you beat treewidth? Theory of Computing
6(1), 85–112 (2010)  [24] McCreesh, C., Prosser, P.: A parallel, backjumping subgraph isomorphism algorithm using supplemental graphs. In: CP 2015. LNCS, vol. 9255, pp. 295–312. Springer (2015)
 [25] McCreesh, C., Prosser, P., Solnon, C., Trimble, J.: When subgraph isomorphism is really hard, and why this matters for graph databases. J. Artif. Intell. Res. 61, 723–759 (2018)
 [26] Slota, G.M., Madduri, K.: Fast approximate subgraph counting and enumeration. In: ICPP 2013. pp. 210–219. IEEE Computer Society (2013)
 [27] Slota, G.M., Madduri, K.: Parallel colorcoding. Parallel Computing 47, 51 – 69 (2015)
 [28] Solnon, C.: AllDifferentbased filtering for subgraph isomorphism. Artif. Intell. 174(1213), 850–864 (2010)
 [29] Solnon, C., Damiand, G., de la Higuera, C., Janodet, J.C.: On the complexity of submap isomorphism and maximum common submap problems. Pattern Recognition 48(2), 302 – 316 (2015)
 [30] Tamaki, H.: Positiveinstance driven dynamic programming for treewidth. In: ESA 2017. LIPIcs, vol. 87, pp. 68:1–68:13. Schloss Dagstuhl (2017)
 [31] Ullmann, J.R.: An algorithm for subgraph isomorphism. J. ACM 23(1), 31–42 (1976)

[32]
Ullmann, J.R.: Bitvector algorithms for binary constraint satisfaction and subgraph isomorphism. J. Exp. Algorithmics
15, 1.6:1.1–1.6:1.64 (2011)  [33] Zhao, Z., Khan, M., Kumar, V.S.A., Marathe, M.V.: Subgraph enumeration in large social contact networks using parallel color coding and streaming. In: ICPP 2010. pp. 594–603. IEEE Computer Society (2010)
Comments
There are no comments yet.