Graph Invariant Kernels

1 Introduction

Renewed interest has been manifested in recent years on graph kernels which can handle continuous (possibly high dimensional) attributes. A vast body of literature on Graph kernels is devoted to symbolic only structures. Which means that vertices (and possibly edges) are labeled by a number of discrete attributes. A popular approach to define graph kernels in this case is to count the number of common substructures (which we call patterns in this paper). Examples include tree-structured patterns [Ramon and G¨artner, 2003; Mah´e and Vert, 2009], paths [Kashima et al., 2003] or shortest paths [Borgwardt and Kriegel, 2005], cyclic and tree patterns [Horv´ath et al., 2004], neighborhoods [Costa and De Grave, 2010], or arbitrary frequent patterns [Deshpande et al., 2005].

Graphs with continuous attributes have been much less investigated, but one graph kernel that handles continuous attributes, is based on random walks (RW) labels [Kashima et al., 2003]. The RW kernel can be computed in O(V 2 ). Borgwardt and Kriegel (2005) proposed a variant based on shortest-paths (SP), which avoids the tottering problem of RW kernels. The SP kernel has a running time of O(V 4 ). This has been recently improved by Feragen et al. (2013) who introduced the GRAPHHOPPER kernel, also based on shortestpaths, reporting a running time of O(V 2 (E + log V + ∆2 ))

(where ∆ is the graph diameter). Kriege and Mutzel (2012) propose graph kernels in which subgraphs are matched with a score which computed as the product between a weight function and the kernels on vertex and edge attributes.

We upgraded existing graph kernels to continuous attributes by using graph and vertex invariants. Vertex invariants are functions that color vertices of a graph in a way that is not affected by isomorphism. They form the basis for several practical isomorphism checking algorithms [McKay and Piperno, 2014]. We consider the commonalities between graph kernels like the Weisfeiler-Lehman graph kernel (WLGK) [Shervashidze et al., 2011], the neighborhood subgraph pairwise distance kernel (NSPDK) [Costa and De Grave, 2010], the propagation kernels [Neumann et al., 2012] and GRAPHHOPPER [Feragen et al., 2013] and summarize them in a general formulation which we call Graph Invariant Kernels (GIK, pronounce "Geek"). GIKs decompose graphs into sets of vertices which are compared by a kernel that measures both their attribute and their structural similarity. The structural similarity indicates to which extent vertices play the same role in the graph they belong to. Our formulation allows arbitrary patterns (e.g. other than the shortest paths used by GRAPHHOPPER) and arbitrary graph and vertex invariants that can be obtained with color propagation schemas (e.g. Weisfeiler-Lehman, Propagation kernel). We also propose spectral coloring which exploits eigen-decompositions.

We show that the upgrade to continuous values provided by GIKs perform very well on a number of new and existing benchmarks. We experiment with different types of vertex invariants, including Weisfeiler-Lehman and spectral colors and compare the shortest paths used by GRAPHHOPPER with neighborhood subgraphs.

2 Notation and definitions

We consider undirected graphs G = (V,E) where V is the vertex set and E the edge set. Both vertices and edges may be labeled. ℓ : V ↦→ X is the vertex labeling function where X may incorporate both discrete and continuous attributes. When necessary, we distinguish continuous and discrete labels as ℓc(v) and ℓd(v), respectively. Similarly, we denote by l : E ↦→ Z the edge labeling function. When needed we represent the connectivity of the graph G with the adjacency matrix A. Two graphs G and G′ are isomorphic, written G ≈ G′ , if there exist a bijection f : V ↦→ V ′ (called an isomorphism) such that {u, v} ∈ E iff {f(u), f(v)} ∈ E′ . In the case of labeled graphs, we also require ℓ(v) = ℓ(f(v)) and l(u, v) = l(f(u), f(v)). An invariant I is a function on graphs such that G ≈ G′ implies I(G) = I(G′ ). If the reverse implication is also true, then I is called a complete invariant. A vertex invariant is a function L : V ↦→ C that assigns each vertex v of G a value L(v), called the color of v, that is preserved under isomorphism f, i.e. L(v) = L(f(v)). The symbol ⪯ denotes the lexicographical order on discrete colors.

3 Graph Invariant Kernels

GIKs measure the similarity between two attributed graphs G and G′ by comparing their vertices with a suitable kernel function kATTR(v, v′ ) between the continuous attributes ℓc and reweighing their similarity with a function w(v, v′ ):

We define w(v, v′ ) as a count on common graph invariants:

The weight w(v, v′ ) measures the structural similarity between vertices and can be designed combining an Rdecomposition relation [Haussler, 1999], a function δm(g, g′ ) and a kernel on vertex invariants kINV.

An R-decomposition relation is a binary relation R(G, g) which encodes that "g is part of G" and specifies a decomposition of G into its parts (patterns). We denote with R-1 (G) the multiset of all patterns in G.

According to Equation 2 the weight between a pair of vertices increases whenever the two vertices appear in the same pattern with the same structural role.

The function δm(g, g′ ) is used to determine whether the two patterns match, while the indicator function 1{v ∈ Vg ∧ v′ ∈ Vg′} is introduced to select only the subgraphs g and g′ in which the vertices v and v′ are involved respectively.

The kernel function kINV is used to measure the similarity between the colors produced by a vertex invariant L and encodes the extent to which the vertices play the same structural role in the pattern. A complete vertex invariant gives the most fine-grained matches and has the same effect as using an isomorphism map f, while weaker invariants induce spurious matches. Weaker invariants can be desirable because they allow to compare non isomorphic graphs. Vertex invariants can be weakened by exploiting the R-decomposition relation and computing the vertex invariants with respect to the patterns. We specialize the notation for kINV(v, v′ ) into local kg INV(v, v′ ) and global kG INV (v, v′ ) as to distinguish whether the vertex invariants in question were computed with respect to the patterns g ∈ R-1 (G) or the whole graph G respectively. We now explain how a suitable δm(g, g′ ) function can define a hard-match between subgraphs generated by some R-decomposition relation, while we refer the reader to Section 4 for the vertex invariants that can be used to define structural similarity kINV(v, v′ ) between vertices.

3.1 Hard-match kernels for discrete labels

If graphs are labeled by discrete symbols, a simple choice for δm(g, g′ ) is the hard match function

for a suitable definition of the equivalence relation ≡. By instantiating the definition of R and ≡, we can obtain a fairly large family of hard pattern-match graph kernels. For example, in the case of SP kernels [Borgwardt and Kriegel, 2005] R(G, g) iff g is a shortest-path between some pair of vertices u and v in G. If edges are unlabeled, then a special case of the kernel presented in [Borgwardt and Kriegel, 2005] can be interpreted as a hard pattern-match kernel by defining g ≡ g′ iff the two serializations (obtained by concatenating the vertex labels) of g and g′ are identical. A second example is the kernel described in [Horv´ath et al., 2004], where R(G, g) iff g is either a cycle in G, or a tree in the forest obtained by deleting all cycles from G. For this kernel, the equivalence g ≡ g′ is established by taking all possible serializations (for example all cyclic permutations in the case of cycles) and checking that the lexicographical minima are identical. As a last example, in the NSPDK [Costa and De Grave, 2010] the decomposition relation is parameterized by two natural numbers r and d and Rr,d(G, g) iff g is a pair of neighborhood subgraphs of G, g1 and g2, rooted at vertices v1 and v2, respectively, such that the two following conditions hold: (1) the shortestpath distance between v1 and v2 is d, and (2) for i = 1, 2, gi is the neighborhood subgraph rooted in vi. A neighborhood subgraph rooted in vi consists of the subgraph induced by all vertices of G whose shortest-path distance from vi is at most r. In [Costa and De Grave, 2010], ≡ is defined as the graph isomorphism relation, approximated by hashing a canonical representation of the patterns.

A weaker hard-match kernel can be obtained using a graph invariant I, i.e. g ≡ g′ if I(g) = I(g′ ). A simple graph invariant IV is the number of vertices of a graph. Vertex colors can be used to construct stronger graph invariants because they are preserved under graph isomorphism. A graph invariant can be obtained as the lexicographically sorted set of vertex colors or edge colors, where the color of an edge {u, v} ∈ E can be represented as (L(u), L(v), l(u, v)). The latter method was introduced by [Costa and De Grave, 2010] combined with distance coloring as in Section 4.3. Graph invariants are hashed to integers so that we can compute the subgraph match function δm(g, g′ ) in constant time. Collisions due to hashing can only weaken the graph invariant.

On the other hand graph invariants used by δm are always computed on patterns (i.e. at local level). Computing δm at the global level would not allow to capture meaningful correlations between a pair of graphs G and G′ leading to poor generalization performance.

GIKs can also be expressed as R-convolutional kernels on graphs [Haussler, 1999] so they are positive semidefinite (PSD).

4 Vertex invariants

In this section we review some possible vertex invariants to construct specific instances of the vertex coloring kernel.

4.1 Weisfeiler-Lehman coloring

The 1-dimensional Weisfeiler-Lehman (WL) refinement algorithm [Weisfeiler, 1976] finds a partition of the vertices of a graph in an iterative fashion:

• Initialization: All colors are initialized using vertex labels, i.e. ∀v ∈ V L(0) (v) .= ℓd(v)
• Recoloring step: ∀v ∈ V

where N(v) is the set of vertices adjacent to v and the function id returns the string of the vertex colors of the neighbors sorted in lexicographic order. Often a hash function is used to represent this string with an integer.

• Termination criterion: recoloring converges when the number of distinct colors stops increasing, which means that the vertices in the graph cannot be further partitioned. In practice the recoloring process is stopped after a predefined number h of iterations.

The 1-dimensional WL color of a vertex is finally obtained as the concatenation of the colors at all time steps:

For edge-labeled graphs, we propagate the vertex colors of the neighbors together with the edge labels:

The positive definite kernel on this kind of vertex invariant can be defined as:

4.2 Label updates in propagation kernels

When introducing propagation kernels (PROP), Neumann et al. (2012) proposed two WL-like label update approaches:

• Diffusion updates:

where T = D-1A is the transition matrix derived from the normalization of the rows of the adjacency matrix, and D is the diagonal degree matrix Dii = ∑Vk=0 aik

• Label propagation: as above except that before each iteration t + 1 the labels ℓd(vi) of the originally labeled vertices overwrite the value assigned to L(t) (vi) by the propagation process.

4.3 Distance coloring

This method was initially suggested in [Costa and De Grave, 2010] to compute the NSPDK and assumes connected and rooted patterns, i.e. g has a distinguished vertex r.

• Starting from r, a breadth-first visit is used to compute all-pairs distances d(v, w) between the vertices.
• Each vertex is colored as follows:

Information about the root of the subgraph can be included by adding to each vertex color information about its distance from the root vertex.

4.4 Spectral coloring

Eigenvalues and eigenvectors of a graph adjacency matrix are insensitive to vertex permutation. For this reason, starting from the seminal work of Umeyama (1988), spectral methods have been popular in pattern recognition as graph matching tool. Spectra can be also used for isomorphism testing if eigenvalues have bounded multiplicity [Babai et al., 1982].

A weighted adjacency matrix W can be defined by using vertex labels ℓc(v), assuming that X is a metric space endowed with a distance function d:

where aij are the entries in the adjacency matrix A of the graph and γ is a non-negative scalar hyperparameter of the heat kernel. In spectral coloring, a vertex invariant is obtained from the the Laplacian embedding of the graph as follows. Let L .= D -W denote the graph Laplacian matrix, with D the diagonal degree matrix. The embedding X ∈ RV×V is the solution of

This equation finds vertex coordinates x(v) (rows of X) whose euclidean distance preserves the graph connectivity and is equivalent to:

the solution can be obtained by solving an eigenvalue problem:

where the eigenvectors xi (columns of X) are sorted according to the corresponding eigenvalues, i.e. 0 = λ1 ≤ λ2 · · · ≤ λV . The row x(v) of X is the embedding coordinate of the vertex v in the graph. The Spectral vertex invariant is obtained by switching to zero the λ-eigenvectors whose eigenvalue λ has multiplicity µ > 1 and taking the absolute value of the components of the eigenvalues otherwise.

Some lifted inference approaches [Mladenov et al., 2012] perform color passing on factor graphs to cluster symmetric variables and speed up intractable inference problems. We solve a complementary problem, and propose the canonized Laplacian embedding to reveal the symmetries between the vertices of a graph. The graph Laplacian is preferred over the one of the adjacency matrix because its spectrum gives more insight on the connectivity of the graph. The smallest eigenvalue λ1 of the Laplacian spectrum of a connected graph is 0, has multiplicity 1 and its λ1-eigenvector is constant. All the other eigenvectors are orthogonal to the constant eigenvector and thus balanced. It is this appealing property that allows spectral coloring features to encode the symmetries of a graph. From an implementation point of view it is possible to limit the dimensionality of the embedding to the first h components. In case there are less than h vertices the embedding can be padded with zeros.

5 Algorithmic issues and running time

We limit our analysis to kernels that use an R-decomposition relation that generates connected subgraphs whose enumeration scales linearly with the number of vertices in the graph G. The time complexity of our method is V 2 (C1 + C2τV 2 g ) where V is the number of vertices, Vg is the number of vertices in a subgraph, C1 and C2 are two costs and τ is the subgraph matching count:

We assume that the cost of computing kINV or kATTR scales with the dimension dINV and dATTR of the corresponding features. According to Equation 2 we have C1 = dATTR and C2 = dgINV. If we use global instead of local graph invariants we can rewrite Equation 2 as:

and thus we obtain C1 = dATTR + dGINV and C2 = 1.

The worst case scenario for τ occurs when the subgraph matching function δm(g, g′ ) is always 1, in this case we have τ = V 2 . This can happen in two cases: 1) when the graph invariant I is constant and every pair of patterns matches and 2) when the patterns are all identical. In the former case we should resort to a more selective invariant, while in the latter case we should avoid computing w(v, v′ ) because its value is constant. This case in which all the patterns are identical is verified for 0-neighborhood subgraphs (vertices) or rneighborhood subgraphs in which r ≥ ∆2 , where ∆ is the diameter of the G. When the radius r is at least twice the diameter ∆ of the graph G, all the corresponding r-neighborhood subgraphs are the graph G itself. In both cases there is no structural contribution, these cases are easy to detect and being aware of this fact we can compute the kernel in O(V 2C1). Another option is to add structural information in order to avoid the extreme cases. We can incorporate discrete labels if present when r = 0 and root information when r ≥ ∆2 . Two rooted patterns must have both same graph invariant and same root vertex invariant in order to match . We obtain an upper bound on the number of vertices Vg of a r-neighborhood subgraph with maximum vertex degree d and diameter ∆ = 2r using the Moore bound [Miller and ˇSir´an, 2005]:

So Vg ∈ O(d2r ) and the complexity of our method is O(V 2 (C1 + C2τ d4r )).

6 Experimental evaluation

We answer the following experimental questions:

• Q1 How do GIKs compare to the state of the art?
• Q2 Do our kernels produce more accurate classifiers when continuous attributes are introduced?
• Q3 Can the graph representation of a sentence benefit from word vector attributes instead of discrete word attributes?
• Q4 What are the best graph invariants?

6.1 Datasets

PROTEINS and ENZYMESSYMM are sets of proteins from Dobson and Doig (2003) and from the BRENDA database [Schomburg et al., 2004], respectively. Vertices are secondary structure elements as in [Borgwardt et al., 2005].

SYNTHETICNEW is a dataset of 300 examples with two classes, it is based on a random graph G with 100 nodes and 196 edges, whose vertices were annotated with scalar attributes sampled fromN(0, 1). It was introduced in [Feragen et al., 2013].

FRANKENSTEIN is a dataset created by the fusion of the BURSI and MNIST datasets. BURSI is made by 4337 molecules with mutagenicity (AMES) classification, with 2401 mutagens and 1936 nonmutagens [Kazius et al., 2005]. Each molecule is represented as a graph whose vertices are labeled by the chemical atom symbol and edges by the bond type.

FRANKENSTEIN is a modified version of the BURSI dataset: we discarded bond type information and remapped the most frequent atom symbols (vertex labels) to MNIST digit images. The original atom symbols can only be recovered through the high dimensional MNIST vectors of pixel intensities, in this sense this is a challenging problem for a graph kernel that can handle continuous attributes.

QC and WEASEL are natural language processing datasets: QC is a dataset for question classification with fixed split (5452 train / 500 test) originally proposed in [Li and Roth, 2002]. We consider the classification task with six coarse labels. WEASEL is part of the CoNNL-2010 shared task dedicated to the detection of hedge cues, it is a binary classification task with fixed split (11111 train / 9634 test).

For both QC and WEASEL we represent a sentence as a graph whose vertices are tokens annotated with 300dimensional word-vector attributes, and whose edges represent the word adjacency and the typed dependency relations extracted with the Stanford dependency parser [De Marneffe et al., 2006]. Word vectors were obtained using the WORD2VEC software [Mikolov et al., 2013] and the whole Wikipedia corpus for training.

did the same for the γINV of the RBF kernel on the continuous vertex invariants derived from spectral graph coloring LSC. The number of dimensions dINV of LSC was set equal to the average number of vertices in a graph for GSGK and to the average number of vertices in a 3-neighborhood subgraph for LSGK. For FRANKENSTEIN we set γ = 0.0073, this value was taken from [Sevakula and Verma, 2012], while for QC we set γ = 1. This value was found during preliminary experiments in which we used the n-grams of word vectors (n = 1, 2, 3) matched with the RBF kernel.

6.3 Experiments

E1 In Table 1 we show the classification accuracy that we achieved without (kATTR = 1) and with (kATTR = RBF(γ)) kernel on continuous attributes. We use the bold font for the best results of each column. For FRANKENSTEIN we also measured the area under the roc curve (AUROC) which is 0.74 for GRAPHHOPPER with continuous attributes. With GWLWL we obtained 0.82 AUROC without continuous attributes and 0.86 with continuous attributes. For all datasets we used SVM-classifiers. Except for QC and WEASEL the accuracy was estimated by 10-times 10-fold cross-validation reporting means and standard deviations [Feragen et al., 2013]. QC and WEASEL have predefined train-test splits. The SVM regularization parameter was selected with an internal k-fold crossvalidation on the training data (k = 3 except k = 10 for QC and WEASEL).

E2 In Table 3 we report the accuracies obtained on the original implementations of WL and PROP using words instead of word vectors. The symbol T denotes the number of iterations and TV and w are parameters (see [Neumann et al., 2012]).

6.2 Kernels

We combined vertex invariants to construct some GIKs that can handle continuous attributes (see also Table 2).

NSK [Costa and De Grave, 2010] was instantiated using an R-decomposition relation that generates neighborhood subgraphs (NS) of increasing radius r = 0, . . . ,R and δm as an indicator function that matches two r-neighborhood subgraphs only if they have the same radius r and IWL subgraph invariant. Costa and De Grave (2010) originally used distance coloring IDC.

LWL and GWL are the local and global variant of the WL subtree kernel respectively. Strictly speaking, GWL is not equivalent to the WL subtree kernel. In fact instead of using an R-decomposition relation which generates vertex patterns we used r-neighborhood subgraphs.

LSGK and GSGK are local and global version of Spectral Graph Kernel in which we used our spectral coloring method.

Each GIK was instantiated in combination with two different kernels on subgraph invariants: IV and IWL. IWL is more selective and reduces the subgraph matching count τ. When necessary we use the subscript to specify which subgraph invariant was employed. We used an R-decomposition relation which generates neighborhood subgraphs of increasing radius r = 0, . . . ,R. Based on preliminary experiments, we fixed R = 3. The Gram matrices were normalized as proposed in [Costa and De Grave, 2010]: for each radius a different Gram matrix is extracted then normalized, we compute their sum and apply normalization again. For ENZYMESSYMM, PROTEINS, SYNTHETICNEW we chose the parameter γ of the RBF kernel as in [Feragen et al., 2013], to be 1/dATTR which is the inverse of the number of dimensions of the attributes. We These results can be compared with Table 1.

E3 In Table 4 we show the results on WEASEL for LWLV and GRAPHHOPPER. We used word vectors and no task specific knowledge. We compare to the state of the art which exploits task specific knowledge encoded in word lists. We set the radius of the r-neighborhood subgraphs to r = 0, 1 as done by Verbeke et al. (2012).

The experiments were run on a 16 cores machine (Intel Xeon CPU E5-2665@2.40GHZ and 96GB of RAM). The code of GIKs was written in the Python programming language, while for GRAPHHOPPER we used the MATLAB implementation from Feragen et al. The difference between the programming languages makes hard to perform fine-grained time measurements. We measured the wall-clock execution time.

In Table 5 we show the runtime of our most accurate and fast GIK selected from Table 1 compared to the runtime of GRAPHHOPPER. The experiment that had the highest runtime was on the dataset WEASEL: LWLV terminated in 73h 47' while GRAPHHOPPER terminated in 75h 43'. We also mention that the spectral kernels tend to be slower and in some cases (e.g. on FRANKENSTEIN) can have a slowdown of a factor of 2 with respect to GRAPHHOPPER.

6.4 Discussion

A1 Our experiments show that we obtained state-of-the-art results for all the datasets in Table 1. ENZYMESSYMM, PROTEINS, QC and WEASEL are are real-world datasets. On SYNTHETICNEW we verified that its random nature combined with the graph invariant IWL leads to diagonal kernel matrices and this explains the poor classification performance. On the other hand if we use the subgraph invariant IV, we outperform GRAPHHOPPER on its own artificial benchmark. In our experiments on FRANKENSTEIN (see Table 1), in some cases (GWLWL) we have an increase of 0.12 points of AUROC with respect to GRAPHHOPPER. This performance mismatch is also due to the nature of the dataset from which FRANKENSTEIN was originated. Indeed neighborhood subgraphs are known to perform well on BURSI [Costa and De Grave, 2010], and shortest paths are probably not the best pattern for this kind of dataset. We consider as upper bound the result obtained with NSPDK (R = 3, D = 0) on BURSI, which is 0.91 AUROC score. The best result obtained on FRANKENSTEIN with continuous attributes (0.86 AUROC score) is significantly higher than the best result obtained without (0.82 AUROC score). On QC we obtain the 94.6% (see Table 1) which is the state of the art. Our result can be compared to the 94.8% obtained by Croce et al. (2011), when they combine kernels on lexical centered trees LCT with a word vector representation obtained with latent semantic analysis. Nevertheless GIKs are a generic tool, not specialized for natural language as LCT and we did not rely on manually built word lists. The 58.3% f1-score (see Table 4) obtained on WEASEL with r = 1 is comparable with the top five results of the CoNNL2010 shared task. The state of the art for the task is 61.5% and was obtained by Verbeke et al. (2012) also exploiting task specific word lists not used by GIKS. The runtime measurements in Table 5 show that we could always instantiate GIKs obtaining comparable or higher accuracy with the same or inferior runtime except for SYNTHETICNEW and QC in which the higher runtime is compensated by higher accuracies.

A2 Comparing the results in Table 1 we see that it is always the case that continuous attributes increase the accuracy.

A3 On QC WL and PROP use words as discrete attributes and cannot obtain more than 89% of accuracy (see Table 3). GIKs successfully use word vectors, indeed LWLV achieves state-of-the-art results (see Table 1).

A4 If we consider the results in Table 1 we notice that LWLV generally tends to give the best results, indeed local versions favor spurious matches and also the subgraph invariant IV is weaker than IWL. Weaker invariants lower down the dimensionality of the kernel and thus make the learning system less prone to overfitting. This effect was probably beneficial due to the relatively small size of some of the datasets that we used. The spectral color invariant can obtain good results, but other GIKs yield equally or more accurate classifiers and have better runtime. What makes spectral color invariants interesting is the nice property of providing vertex invariants embedded in the euclidean space. This makes Spectral Graph Kernel amenable for future work on improving scalability using sketching techniques.

7 Conclusion

The GIK reformulation of well known graph kernels allows to obtain more insights in the exploration of graph kernels with continuous attributes. The underlying idea was to employ vertex invariants for soft subgraph matching. We contributed new insights into graph kernels and to upgrade existing ones for use with continuous attributes. Several graphkernel instances were then empirically evaluated on a number of new and existing benchmark datasets. The results showed that some combinations of graph and vertex invariants with continuous attributes lead to excellent performance. For future work we plan to improve runtime and scalability focusing on the local GIKs and spectral coloring which are more appealing for sketching.