Details of our implementation

Consider the graph $G$ and the path $T$ given here:

               0
               A
          _.-'
G:     1 B
         |'.         E 3
         |  `.      ,
       2 C.._ `.  ,'
             ``-D'
                2

T:       a -> b -> c -> d

Each vertex in $G$ has a color in $\{0,1,2,3\}$ which is indicated next to the vertex. Note that we in the following associate with each vertex a value that is a vertex identifier, and hence unique. We model attributes of vertices (and edges) by accessor functions taking the identity as input and yielding the value as output. We associate each subgraph of $G$ with its colorset, i.e., the set of colors assigned to its vertices. We say that subgraph of $G$ is colorful if there are no two vertices in the subgraph that share a color, i.e., the cardinality of its colorset is equal to its order. We now want to find a colorful match for the directed path $T$ of length $m=4$ in $G$, in other words we want a vertex assignment from $T$ to a colorful subgraph of $G$. By inspection we see that there are two possibilities, the colorful path ABDE in $G$, or the reverse of this, starting at E. These are the only possibilities because C and D have the same color 2. Now, consider the case where we have matched the $k$ length path $(a,b,c)$ and wish to extend this with a match for $d$. The following table shows the possibilities.

If we are only interested in answering whether there exists a match, we only have to store the colors and the prefix endpoint (the vertex to which the next vertex will be connected) of any matched $k$ length subpath in $G$.

Now consider each extension of a match having a cost depending on which edge in $G$ is used (independently from previous edges used), and instead of wanting to determine whether a match exists, we want the cheapest match. Again it is enough to store the colorset and the prefix endpoint, but with the caveat that we also need to store information that lets us recover a cheapest path corresponding with a given colorset and endpoint. One way to do this is to also store at each step the value of the cheapest prefix as well as the prefix itself. Alternatively, Huffner [Huffner] proposes storing the order of the colors instead of the prefix itself, and then recover the sequence at the end by backtracking. Storing each entry in the prefix costs $\log(n)$ bits, while storing each entry in the color order costs $\log(m)$ bits, representing savings in storage space at the cost of the recovery.

The above approach of path extension can be thought of as processing of the ordered list of vertices in $T$ much like we would compute the sum of integers in a list using a binary operation. At each step a vertex from $T$ is combined with a table to produce a table. The result is then read out of the final table. Considering that a list is a recursive data structure, this pattern of recursion is called a fold or "reduction" of the list.

In the path matching approach above we extended a prefix of length $i$ with endpoint $u$ matched to $T$ vertex $q_i$ by considering all edges $(u,v)$ in $G$ as a match for edge $(q_i, q_{i+1})$ in $T$. If we think in terms of edges instead of vertices, the fold can equivalently be expressed as a fold over the edges of $T$ in pre-order. For the path problem, we know that each edge to be added originates (is incident to) at the last vertex added. This is not so in general for trees. The next edge origin might be any already visited vertex. Consider the tree $G_1$ with vertices A,B,C,D,E below:

                           ,-"""--.
                          /       |  A root
                         /   A    |
                       ,'   ,'    /
                      /    /     /
                    ,'   ,'     /
B prefix endpoint  /   ,'      /
                 ,'   B._     |
               ,'  ,-  \ `--._|
             ,' ,-'     \     |`-..
             | C     _.  D    /     E    Next child
             L     ,'  `.    |
              `._,'      `--'


                      D last visited

Assume we have already visited vertices A, B, C, and D in pre-order, indicated by the drawn envelope of this subtree. The vertex D is the last visited vertex, A is the root, while the next parent is B representing the prefix endpoint. All the candidate colorful trees in $G_2$ are extended at the point corresponding to B adding an edge to a point which will be put in correspondence with E. Consequently, all trees sharing the same colorset and prefix endpoint also share the same possible extensions, which means that we only need to store information on one tree for each possible (colorset, endpoint) pair in order to determine whether there exists a tree match or compute the minimum cost match. In the following we will use a mapping data structure "table" that maps a (colorset, vertex) pair to stored information.

In the following we turn to details of how we implement the above ideas. We write function applications and definitions much like it is done in the functional programming languages OCaml and Haskell. Let (h::t) denote a list that has h as a first element and the list t containing the rest of the elements. The infix list construction operator :: is often called "cons" and we could have equivalently defined (h::t) in terms of a prefix application like (::) h t when we use the convention that we can convert an infix operator into a prefix function by enclosing it in parentheses. Also let [] denote the empty list. The list containing elements 1,2, and 3 in that order can then be written as (1::2::3::[]) or (::) 1 ((::) 2 ((::) 3 [])). Now recall that (+) x y = x + y, and let us substitute (::) in the above list with (+), and [] with z = 0. Then we get (+) 1 ((+) 2 ((+) 3 z)) = (1 + (2 + (3 + 0))). This substitution of function calls (and z) for type constructors in recursively defined data types is called a fold. For lists as we defined them, we can define a fold function as follows:

fold g z []    = z
fold g z (h::t) = g h (fold g z t)

We can get the above fold example by letting g = (+) and z = 0, i.e., fold (+) 0 (1::2::3::[]).

Using a fold, we can define standard list processing functions in term of fold as:

map f list = fold (x t -> (f x)::t) [] list
filter p list = fold (x t -> if p x then x::t else t) [] list
concat l1 l2 = fold (::) l2 l1

where (arguments -> expression) means the unnamed function that when given arguments arguments returns the value of expression. The function map applies the function f to all elements in list in turn and returns the list of results, filter takes a predicate p and a list list and returns the sub-list of elements in list for which p is true, while concat concatenates two lists. Concatenation of lists is for conveniance often denoted by an infix operator. We will use ++, i.e., a ++ b = concat a b.

Of particular interest for us later is the foldl function:

foldl f z []     = z
foldl f z (h::t) = foldl (f z h) t

This function is called a left fold in that it instead of processing a list in a right to left manner, processes it in a left to right manner instead. So foldl (+) 0 (1::2::3::[]) is ((0 + 1) + 2) + 3. We can define foldl in terms of fold as:

foldl f z l = fold (x g -> (a -> g (f a x))) (x -> x) l z

As Hutton [Hutton] points out, while programs written using fold can be less readable, they can be constructed systematically and can more easily be analyzed and transformed. This, and the expressive power of using folds, are reasons why we in the following express color coding in terms of folds.

Returning to the path matching example, we now see that if we have an operation g that takes as input a table and a vertex and returns a table, we can compute the final table t as:

finaltable vertices = foldl g emptytable vertices

For path matching, we know that the prefix endpoint is the latest seen vertex, hence vertices can be a simple ordered list of vertices. In the tree case, as we have discussed above, this assumption does not hold and vertices must contain information about which vertex is the endpoint for each prefix.

Now consider a recursively defined data type Tree that consists of vertices Vertex v s where v is the vertex label and s is a list of subtree root vertices. Note that this data type has two constructors, Vertex for creating the (label, subtree) tuple, and "cons" or (::) for creating the subtree list. Similarly to fold for lists, we can define a fold on trees by substituting functions for the type constructors (and again an initial element z for [] in the subtree list signifying that Vertex label [] is a leaf):

treefold  f g z (Vertex label subtrees) =
  f label (treefolds f g z subtrees)
treefolds f g z []     = z
treefolds f g z (h::t) = g (treefold f g z h) (treefolds f g z t)

The function f takes the place of the Vertex constructor, the function g takes the place of "cons" in the subtree list, and z the place of []. The auxiliary function subtrees is essentially a fold over the subtree list. We can define the following tree processing functions using treefold:

preorder tree = treefold (::) (++) [] tree
order tree = treefold (y x -> x + 1) (+) 0 tree

Here, preorder returns the list of vertex labels in preorder, while order returns the number of vertices in the tree.

An ordered tree, which our Tree data type represents, is isomorphic to its pre-order edge list. We can obtain this list using treefold as:

f label s = (label, fold ((i, l) res -> (label, i)::(l ++ res)) [] s)
preedges tree = snd (treefold f (::) [] tree)

where snd (a,b) = b. As a path is isomorphic to its preorder vertex list, and we could fold over this in the path matching example above, we can now analogously fold over the pre-order edge list when matching trees:

finaltable tree = foldl g emptytable (preedges tree)

We now turn to the question of how to implement the binary operation g we want to use in the fold over pre-order edge list. To this end let a (sub-)tree in $G$ be stored in a structure that has the information retrieval functions colorset, endpoint, score defined on it. Also, let addvertex tree edge v add the vertex v to the structure instance tree and update information stored in it accordingly using that the current edge in the origin tree is edge. Also let init v produce a structure containing only the root v. Of the above functions, the addvertex function needs to have access to information from which it can compute the new colorset and endpoint of the tree, as well as the score. We will not address this function in further detail. Also let empty denote an empty tree structure.

As discussed above, we need to be able to store trees indexed by their colorsets and endpoints. Let table be a generic map type from (colorset, endpoint) keys to trees. Associated with the map type are functions items which yields a sequence (e.g., a list) of elements (trees) stored in the map instance, update which, given a tree, stores the tree at its (colorset, endpoint) entry either if there is no tree there already, or if the already stored tree has a lower score, and new () which creates a new empty store. Also, let the function max extract a maximum score tree from the store.

We can now define a function that takes a tree and an augmented edge and produces all colorful extensions of this tree as:

extend edge tree =
  [addvertex tree edge v | v <- successors G (endpoint tree) and
                       (color v) not in (colorset tree) ]

The notation [expr | v <- gen and pred] is a list comprehension and means the list of all values of expression expr for all elements v generated by expression gen for which expression pred evaluates to true. Functions successor G u yields all v such that there is an edge (u, v) in graph G, and color v yields the color associated with vertex v in G. We can now define the operation g that takes as input a table instance and an augmented edge and returns a table instance as:

g table edge =
  let op table' tree = foldl update table' (extend edge tree) in
  foldl op (new ()) (items table)

The function g works by iteratively extending each tree in the given table into a list of trees from which an initially empty table is populated. We can then now encode a trial of color coding as:

trial tree =
  let table = foldl update (new ()) [init u | u <- (vertices G)] in
  max (foldl g table (preedges tree))

The trial function first creates an initial table containing only single vertex trees, one for each vertex in $G$. Then it folds g over the augmented edge list computed from tree creating a final table from which a maximum score tree is extracted and returned. Once enough trials have been performed, the best matching tree can be translated into an assignment.

In fact, any subtree produced during the execution of trial can be translated into a sub-assignment. In fact, we can consider sub-trees as sub-assignments and vice versa. Now, given a sub-assignment $f'$ of with value $\sigma'$ we could decide if $f'$ could be extended to an optimum assignment $f$ with value $\sigma^*$ if we knew the largest possible value $\sigma''$ associated with this extension. If we knew this, we could at any stage in the algorithm, and in the computation of g in particular, prune away subassignments that would never be extendable to an optimum assignment, i.e., prune $f'$ if $\sigma' + \sigma'' < \sigma^*$. The time and space savings of this would be significant! Note that while $\sigma''$ and $\sigma^*$ are in general not available, any approximation $\sigma''_{\sim}$ and $\sigma^*_{\sim}$ such that $\sigma''_{\sim} \geq \sigma''$ and $\sigma^*_{\sim} \leq \sigma^*$ can be used. For $\sigma^*_{\sim}$ we can at any time use the best encountered value so far. For the subassignment $f'$ of size $k>0$ we know that the extension to $f$ involves $m-k$ edges and $m-k$ vertices. We can thus find an upper bound $\eta_i + \alpha_i$ for $\sigma''$ by computing $w(e, e'')$ for $(e,e'') \in E(G_1) \times E(G_2)$ and computing $\eta_i$ as the sum of the $i$ largest of these, and computing $v(a,b)$ for $(a,b) \in V(G_1) \times V(G_2)$ and letting $\alpha_i$ be the sum of the $i$ largest of these. We then perform pruning by setting $\sigma''_{\sim} = \eta_{m-k} + \alpha_{m-k}$ when considering $f'$ of size $k$. Alternatively, as Huffner et al. [Huffner] suggest for their their path matching algorithm, we can exhaustively compute $\sigma_i$ for $i \leq l$ for small $l$, and compute $\sigma_i$ for $i > l$ as the smallest sum of $\sigma_j$ such that $j \leq l$ and $\sum_j j = i$. Given a function keep tree that returns false if the sub-assignment associated with tree cannot be extended to full assignment with a score that exceeds the currently best known $\sigma''_{\sim}$, we can perform the filtering by inserting a call to keep into extend as:

extend edge tree =
  filter keep [addvertex tree edge v
               | v <- successors G (endpoint tree)
                 and ((color v) not in (colorset tree)]

Trial Time Complexity Analysis

In the above algorithm, when processing edge $i$ (1-indexed) we are processing the trees in a map data structure that maps colorsets of length $i$ and a prefix endpoint to trees. Consequently, we are processing at most $n{m \choose i}$ trees for edge $i$. Now consider a particular tree with colorset $C$ and endpoint $u$. This particular prefix can only be extended with vertices that are adjacent to $u$ in $G$. This means that we for all trees in the map with colorset $C$, we only need to consider edges in $(u,v) \in E(G)$ for which $u$ is an actually occurring prefix endpoint. In short, for each colorset $C$ we need to consider at most $| E(G) | $ possible extensions, making the total number of extensions we need to process $| E(G) | {m \choose i}$ for edge $i$. Since we are extending as many times as there are edges in the tree, each time with a colorset size that has increased by one, we get that we are processing at most $$ \sum_{i=1}^{m-1} | E(G) | {m \choose i} \in O( | E(G) | 2^m) $$ extensions. We now assume that computing the score of an extension takes constant time, as does accessing the list of successors of a vertex in $G$. The two remaining steps of processing an extension is updating the structure (addvertex), and the update of the map data structure. Updating the tree structure entails adding a color to a colorset, storing the new score, and adding a vertex to a list. For colorsets smaller than the length of machine registers, this update can be considered constant time. This is not a real limitation as the processing of $2^k$ trees becomes impractical for $k$ significantly smaller than modern machine register lengths. There are two subtle points associated with the implementation of the tree structure. It would be easy for addvertex to return a complete copy of the original tree together with updates making up the extension. This would make the running time of addvertex a function of the input. This can be avoided by having all extensions share the common part of the structure, for example by storing the vertices encountered in reverse order in a linked list. The other subtle point is keeping track of the next endpoint. If we keep track of the index of this endpoint and keep the tree structure in a list as proposed, the endpoint function has to traverse a list of $O(i)$ every time for origin edge number $i$. However, since the structure of the tree is known a priori, we can let the tree representation have a stack that can in constant time be used to compute the next endpoint. Like the tree vertex list, extensions of a tree share the common part of the stack. Consequently, the extensions of a tree can each be done in constant time with this scheme at the cost of storing an additional $O(log(m))$ size stack for each tree. If we allocate a $n2^m$ table of empty tree references up front, new does nothing and update takes constant time. Using combinatorial number systems [Knuth] and color sets encoded in register size bitvectors, we can produce all colorsets of a given size $i$ in ${m \choose i}$ time. This means that items for the lookup-table takes at most $n{m \choose i}$ time (it can also disregard empty entries in the table). Since items is called once for each edge, we get that the resulting running time is $$ \sum_{i=1}^{m-1} | E(G) | {m \choose i} + n{m \choose i} \in O( | E(G) | 2^m) $$ for any connected graph. This approach also has a best case running time in $\Omega(| E(G) | 2^m)$, and takes $O(n2^m)$ and $\Omega(n2^m)$ space. At the cost of an additional $O(n2^m)$ in processing time, we can instead of allocating the table up front, allocate a table of size $n{m \choose i}$ at step $i$. Space consumption is then only $O(n { m \choose \lceil m/2 \rceil })$., but the asymptotic running times stay the same. Another option is to use a dynamic map structure that stores only what is needed. Allocating a new structure at step $i$ with new then becomes constant, but update becomes $O(\log (n{m \choose i}))$ resulting in a running time $$ O(\sum_{i=1}^{m-1} | E(G) | {m \choose i}\log(n{m \choose i})) \subseteq O(m | E(G) | 2^m) . $$ This bound is no longer tight, and using a dynamic table can in practice be faster than the $O(| E(G) | 2^m)$ approach when pruning is effective. Also the $O(n { m \choose \lceil m/2 \rceil })$ space bound is not tight, as only actually constructed colorful extensions for each origin edge are stored and processed.

Classification of problem instances

Borrowing the term "deceptive problem instances" from the evolutionary algorithms literature we define deceptive matching problem instances as follows. First, for a problem instance $(G_1, G_2)$, and an assignment $f : V(G_1) \to V(G_2)$, let $\sigma_{(G_1, G_2)}(f)$ be the score of the assignment for this problem instance. Now, let $F_{(G_1, G_2)}$ be the set of optimal assignments all sharing the same optimal score $\sigma^*_{(G_1, G_2)} = \sigma(f)_{(G_1, G_2)}$ for $f\in F_{(G_1, G_2)}$. Also, let $\mathcal{S}_k(G)$ denote the set of order $k$ subgraphs of $G$. For an assignment $f : V(G_1) \to V(G_2)$ we define the set of $k$-deceptive subgraphs of $G_1$ as the set $D^k_{(G_1, G_2)}(f)$ of $k$ order subgraphs for which the restriction of $f$ is not optimal. Formally, $$ D^k_{(G_1, G_2)}(f) = \{G \in \mathcal{S}_k(G_1) | \sigma^*_{(G, G_2)} < \sigma(f|_{V(G)})_{(G, G_2)}\}. $$ For the purpose of this discussion, we will call a problem instance $(G_1, G_2)$ $k$-deceptive if there exists $f \in F_{(G_1, G_2)}$ such that $D^k_{(G_1, G_2)}(f) \neq \emptyset$. Such an $f$ will also be called $k$-deceptive. We could use $D^k_{(G_1, G_2)}$ to build hierarchies of deceptiveness, levels could be defined by requiring that all $f \in F_{(G_1, G_2)}$ be $k$-deceptive, then by intersections of $D^k_{(G_1, G_2)}(f)$ for $f \in F_{(G_1, G_2)}$, being deceptive for $l \leq k$, and so forth. However, an investigation of such lies outside the scope of this discussion. However, we can say that the degree of $k$-deceptiveness of a problem instance $(G_1, G_2)$ is related to the minimum size of $D^k_{(G_1, G_2)}(f)$ taken over all $f \in F_{(G_1, G_2)}$.

We now define quasi-deceptiveness. Formally, let $$ Q^k_{(G_1, G_2)}(f) = \{G \in \mathcal{S}_k(G_1) | \sigma^*_{(G, G_2)} \leq \sigma(f|_{V(G)})_{(G, G_2)}\}. $$ We now have that $D^k_{(G_1, G_2)}(f) \subseteq Q^k_{(G_1, G_2)}(f)$. Analogous to k-deceptiveness we define $(G_1, G_2)$ quasi $k$-deceptive if there exists $f \in F_{(G_1, G_2)}$ such that $Q^k_{(G_1, G_2)}(f) \neq \emptyset$. This $f$ will also be referred to as being quasi $k$-deceptive. If an instance or assignment is quasi $k$-deceptive but not $k$-deceptive we call it strictly quasi $k$-deceptive. Again, we say that the degree of (strictly) quasi $k$-deceptiveness of $(G_1, G_2)$ is related to the minimum size of ($Q^k_{(G_1, G_2)}(f) - D^k_{(G_1, G_2)}(f)$) $Q^k_{(G_1, G_2)}(f)$ taken over all $f \in F_{(G_1, G_2)}$.

Our subsampling strategy rests on the assumption that we can treat deceptiveness as noise that can be addressed by sub-sampling. We also note that any non-weighted subgraph monomorphism problem is at most strictly quasi deceptive when a solution exists.

Ranking performance

The overall average rank/correct ratio is 0.85. This means that if we sort the individual correspondences decreasingly according to their strength, a prefix of length $0.85 * c$ consists of correct correspondences. As such, this ratio serves as an indicator of the usefulness of the strength measure computed for each individual correspondence in an assignment. Note that the normalization achieved by dividing by $c$ makes this independent from the quality of the assignment.