Graph Kernels, in depth: Weisfeiler-Lehman subtree features in code, the 1-WL ceiling shared with GNNs, shortest-path, random-walk and graphlet kernels, and kernel SVMs on graph data

What a graph kernel must be

A useful graph kernel needs three properties. It must be positive semidefinite: every Gram matrix it produces has no negative eigenvalues, which is what lets an SVM's dual problem stay convex. The easiest way to guarantee that is to define an explicit feature map phi and set k(G, G') = <phi(G), phi(G')>. It must be permutation invariant: renumbering vertices cannot change the value. And it must be computable in polynomial time, which rules out the ideal choice.

That ideal would be a complete kernel, where phi(G) = phi(G') only when G and G' are isomorphic. Gaertner, Flach and Wrobel showed in 2003 that computing any complete graph kernel is at least as hard as deciding graph isomorphism, and computing the inner product of full subgraph-count features is NP-hard. So every practical graph kernel is a deliberate compromise. It counts some substructure that is cheap to enumerate, accepts that some non-isomorphic graphs collide, and hopes the substructure it counts is the one that predicts the label.

The family map

Keep this picture in mind when choosing a kernel. The feature map decides which graphs look alike, and that decision is far more important than the classifier on top.

The Weisfeiler-Lehman subtree kernel

The WL kernel (Shervashidze et al., JMLR 2011) runs the colour-refinement step of the 1-dimensional Weisfeiler-Lehman isomorphism test. Every vertex starts with its label (or a constant if unlabelled). In each iteration, a vertex's new label is a compressed ID for the pair (own label, sorted multiset of neighbour labels). After h iterations, phi(G) is the histogram of every label seen in iterations 0 to h. A label at iteration i identifies a rooted subtree pattern of depth i, hence the name.

The compression dictionary must be shared across all graphs, or the same pattern gets different IDs in different graphs and the dot products become meaningless:

from collections import Counter

def wl_features(graphs, h=3):
    # graphs: list of (adj, labels); adj is a list of neighbour lists, labels a list of str
    table = {}                      # (own, neighbour multiset) -> compressed id, shared by all graphs
    feats = [Counter(lbl) for _, lbl in graphs]
    current = [list(lbl) for _, lbl in graphs]
    for it in range(h):
        nxt = []
        for gi, (adj, _) in enumerate(graphs):
            lab = current[gi]
            new = []
            for v in range(len(adj)):
                key = (lab[v], tuple(sorted(lab[u] for u in adj[v])))
                if key not in table:
                    table[key] = f"{it}:{len(table)}"
                new.append(table[key])
            feats[gi].update(new)
            nxt.append(new)
        current = nxt
    return feats

def kernel(f, g):
    return sum(cnt * g[k] for k, cnt in f.items() if k in g)

Cost: each iteration touches every edge once (plus a sort per vertex), so features cost O(h m) per graph up to log factors, and a Gram matrix over N graphs costs N² sparse dot products. That near-linear feature cost is why WL became the default baseline for graph classification.

Worked example: a cycle, two triangles and a path

Take three unlabelled graphs on six vertices: the cycle C6, two disjoint triangles (2C3), and the path P6. Run the code above with h = 3 and every vertex starting with label "1".

In C6 and in 2C3, every vertex has exactly two neighbours, and every neighbour has two neighbours, and so on. Every vertex therefore gets the same label at every iteration, and both graphs produce the same histogram: one label with count 6 at each of the four levels. The unnormalised kernel values are k(C6, C6) = k(2C3, 2C3) = k(C6, 2C3) = 144. After cosine normalisation, k(G, G') / sqrt(k(G, G) k(G', G')), the similarity is exactly 1.0, although one graph is connected and the other is not, and one has triangles and the other has none.

P6 is different because its two endpoints have degree 1. Refinement separates endpoints from interior vertices at iteration 1, and interior vertices next to endpoints from the central pair at iteration 2. We measured k(P6, P6) = 80 and k(P6, C6) = 72, so the normalised similarity is 72 / sqrt(144 x 80) = 0.67. The lesson is that WL sees degree sequences and local tree-unfoldings, never cycles as such. Any property the label depends on, such as ring membership in chemistry, must be either visible in those unfoldings or supplied as a vertex label.

The 1-WL ceiling, and why GNNs share it

Xu, Hu, Leskovec and Jegelka ("How Powerful are Graph Neural Networks?", ICLR 2019) and, independently, Morris et al. (AAAI 2019) proved that any message-passing GNN that aggregates neighbour states and updates a vertex state is at most as discriminative as 1-WL. GIN, whose aggregator is an injective sum followed by an MLP, reaches that bound. So the C6-versus-2C3 collision above is also a collision for GCN, GraphSAGE and GIN with constant input features. A WL kernel is therefore a strong baseline for a GNN, not a primitive one. When a GNN fails to beat it on small datasets, that is common in published benchmarks and not a bug.

Ways past the ceiling, roughly in order of cost:

• Add informative vertex and edge labels (atom type, bond order). Most real collisions disappear.
• Add structural features such as triangle counts or cycle membership as initial labels.
• Use higher-order k-WL variants that colour tuples of vertices. They are strictly more expressive, and the work grows as n^k.
• Combine kernels: a sum of PSD kernels is PSD, so WL plus a graphlet kernel covers both views.

Shortest-path, random-walk and graphlet kernels

Shortest-path kernel (Borgwardt and Kriegel, 2005). Compute all-pairs shortest paths (Floyd-Warshall, O(n³)), and describe the graph by the multiset of triples (label of u, label of v, distance). With discrete labels, phi is a histogram of those triples and the kernel is a dot product. This kernel sees global distances, so it separates C6 from 2C3 immediately: 2C3 has unreachable pairs and no distance 3. Cubic cost and dense distance matrices limit it to graphs of a few hundred vertices.

Random-walk kernel (Gaertner et al. 2003; Vishwanathan et al., JMLR 2010). Count pairs of matching walks in G and G' by walking on their direct product graph, whose adjacency for unlabelled graphs is the Kronecker product A_x = A kron A'. The geometric version sums lambda^k A_x^k over all lengths, giving 1^T(I - lambda A_x)^-11. That converges only if lambda is below 1 / rho(A_x), and for the Kronecker product rho is the product of the two spectral radii. Two problems. Naively the product graph has n n' vertices, so a direct solve costs O(n⁶) for equal sizes, and Vishwanathan et al. reduce it to roughly cubic with Sylvester-equation and conjugate-gradient solvers. Also, walks may step back and forth along one edge (tottering), which floods the count with uninformative walks and makes small lambda necessary, and that in turn makes the kernel almost a count of matching edges.

Graphlet kernel (Shervashidze et al., AISTATS 2009). Count occurrences of every non-isomorphic subgraph on 3, 4 or 5 vertices, usually by sampling, and normalise to a frequency vector. It captures triangles and small cycles that WL misses, but it ignores vertex labels in its basic form, and exhaustive counting of 5-vertex graphlets is expensive on dense graphs.

Training with a precomputed kernel

With features or a Gram matrix in hand, training is standard, with one twist: scikit-learn's SVC accepts kernel="precomputed", and at prediction time it needs the kernel between test and training graphs, shape (n_test, n_train).

import numpy as np
from sklearn.svm import SVC
from sklearn.model_selection import StratifiedKFold

def gram(fa, fb):
    return np.array([[kernel(a, b) for b in fb] for a in fa], dtype=float)

def cosine_normalise(K, da, db):
    return K / np.sqrt(np.outer(da, db))

feats = wl_features(graphs, h=3)          # graphs and y defined by your loader
diag = np.array([kernel(f, f) for f in feats])
for tr, te in StratifiedKFold(10, shuffle=True, random_state=0).split(feats, y):
    Ktr = cosine_normalise(gram([feats[i] for i in tr], [feats[i] for i in tr]), diag[tr], diag[tr])
    Kte = cosine_normalise(gram([feats[i] for i in te], [feats[i] for i in tr]), diag[te], diag[tr])
    clf = SVC(kernel="precomputed", C=1.0).fit(Ktr, y[tr])
    print((clf.predict(Kte) == y[te]).mean())

Select h and C inside each training fold (nested cross-validation). Selecting them on the test folds is the most common way graph-kernel papers overstate accuracy. Also note that the WL dictionary above was built over all graphs at once. That is harmless for features that are only counted, but if you later add feature selection or IDF-style weighting, fit those on training folds only.

Operational guidance

• Use explicit features when N is large. The Gram matrix is N²: 100,000 graphs means 10¹⁰ entries. WL features are sparse vectors, so put them in a sparse matrix and train a linear SVM or logistic regression directly, which is the same model in the primal, at linear cost in N.
• Normalise. Without normalisation, large graphs have large self-similarity and dominate. Cosine normalisation keeps the kernel PSD and puts every diagonal entry at 1.
• Watch diagonal dominance. As h grows, labels become so specific that each graph matches mostly itself. K approaches a scaled identity matrix and the SVM memorises. If training accuracy is 100% and validation accuracy collapses as h rises, lower h (2 to 5 is typical) or use a weighted sum over iterations that down-weights deep ones.
• Hash carefully. Implementations that hash label tuples instead of using a shared dictionary must use a stable hash across processes (not Python's salted hash() for strings), or parallel workers produce incompatible features.
• Check PSD numerically. Some popular similarity scores, such as optimal-assignment variants or graph edit distance turned into a kernel, are not guaranteed PSD. Inspect the smallest eigenvalue of a sample Gram matrix before trusting an SVM trained on it.

Failure modes

• Regular graphs collide under WL. As in the worked example, any two d-regular graphs of equal size with no vertex labels are indistinguishable at every h.
• Label leakage through the dictionary size. Rare WL labels that occur only in one class can act like molecule IDs. That shows up as near-perfect accuracy on duplicated or near-duplicated graphs. De-duplicate isomorphic graphs across train and test before evaluating.
• Continuous attributes. WL and shortest-path histograms need discrete labels. Binning continuous attributes throws information away. Use kernels designed for attributes (for example hashing or propagation-based ones) or a GNN.
• Random-walk kernel numerics. Choosing lambda above 1 / rho makes the series diverge. Choosing it far below makes the kernel degenerate to edge counts. Compute rho per dataset.
• Memory blow-up. All-pairs shortest paths on a 10,000-vertex graph is a 100-million-entry dense matrix per graph.

Trade-offs

Read the kernel trick and SVMs for why PSD matters and how the dual uses K, graph neural networks for message passing and its expressivity limits, graph attention networks for attention-based aggregation, and tree isomorphism for the canonical-labelling idea that WL relabelling generalises.

What to do next

1. Write down which substructure should predict your label (rings, distances, motifs), and choose the kernel family that counts it.
2. Implement or install a WL kernel with a shared label dictionary, and verify it on C6 versus P6 before trusting it.
3. De-duplicate isomorphic graphs across your train and test splits.
4. Run nested 10-fold cross-validation over h in 1..5 and C on a log grid, with cosine-normalised Gram matrices.
5. Inspect the Gram matrix: smallest eigenvalue (PSD check) and how close it is to the identity (diagonal dominance).
6. If N exceeds a few tens of thousands, switch to explicit sparse features and a linear model.
7. Compare against a GIN baseline. If neither beats WL, add structural vertex labels before reaching for bigger models.