Improving Graph Neural Network Expressivity via Subgraph Isomorphism Counting

Introduction

The field of graph representation learning has undergone a rapid growth in the past few years. In particular, Graph Neural Networks (GNNs), a family of neural architectures designed for irregularly structured data, have been successfully applied to problems ranging from social networks and recommender systems to bioinformatics , chemistry and physics , to name a few. Most GNN architectures are based on message passing , where the representation of each node is iteratively updated by aggregating information from its neighbours.

A crucial difference from traditional neural networks operating on grid-structured data is the absence of canonical ordering of the nodes in a graph. To address this, the aggregation function is constructed to be invariant to neighbourhood permutations and, as a consequence, to graph isomorphism. This kind of symmetry is not always desirable and thus different inductive biases that disambiguate the neighbours have been proposed. For instance, in geometric graphs, such as 3D molecular graphs and meshes, directional biases are usually employed in order to model the positional information of the nodes ; for proteins, ordering information is used to disambiguate amino-acids at different positions in the sequence ; in multi-relational knowledge graphs, a different aggregation is performed for each relation type .

The structure of the graph itself does not usually explicitly take part in the aggregation function. In fact, most models rely on multiple message passing steps as a means for each node to discover the global structure of the graph. However, since message-passing GNNs are at most as powerful as the Weisfeiler Leman test (WL) , they are limited in their abilities to adequately exploit the graph structure, e.g. by counting substructures . This uncovers a crucial limitation of GNNs, as substructures have been widely recognised as important in the study of complex networks. For example, in molecular chemistry, functional groups and rings are related to a plethora of chemical properties, while cliques are related to protein complexes in Protein-Protein Interaction networks and community structure in social networks, respectively .

Therefore, three major questions arise when designing GNN architectures: (a) How to go beyond isotropic, i.e., locally symmetric, aggregation functions ? (b) How to ensure that GNNs are aware of the structural chatacteristics of the graph? (c) How to achieve the above two without sacrificing invariance to isomorphism and hence the ability of GNNs to generalise?

In this work we attempt to simultaneously provide an answer to the above. We propose to break local symmetries by introducing structural information in the aggregation function, hence addressing (a) and (b). In particular, the contribution of each neighbour (message) is transformed differently depending on its structural relationship with the central node. This relationship is expressed by counting the appearance of certain substructures. Since substructure counts are vertex invariants, i.e. they are invariant to vertex permutations, it is easy to see that the resulting GNN will be invariant to isomorphism, hence also addressing (c). Moreover, by choosing the substructures, one can provide the model with different inductive biases, based on the graph distribution at hand.

We characterise the expressivity of our message-passing scheme, coined as Graph Substructure Network (GSN), showing that GSN is strictly more expressive than traditional GNNs for the vast majority of substructures, while retaining the locality of message passing, as opposed to higher-order methods that follow the WL hierarchy (see Section 2). In the limit, our model can yield a unique representation for every isomorphism class and is thus universal. We provide an extensive experimental evaluation on hard instances of graph isomorphism testing (strongly regular graphs), as well as on real-world networks from the social and biological domains, including the recently introduced large-scale benchmarks . We observe that when choosing the structural inductive biases based on domain-specific knowledge, GSN achieves state-of-the-art results.

Preliminaries

Let $G=(\mathcal{V}_{G},\mathcal{E}_{G})$ be a graph with vertex set $\mathcal{V}_{G}$ and edge set $\mathcal{E}_{G}$, directed or undirected. A subgraph ${G_{S}=(\mathcal{V}_{G_{S}},\mathcal{E}_{G_{S}})}$ of $G$ is any graph with $\mathcal{V}_{G_{S}}\subseteq\mathcal{V}_{G},\ \mathcal{E}_{G_{S}}\subseteq\mathcal{E}_{G}$. When $\mathcal{E}_{G_{S}}$ includes all the edges of $G$ with endpoints in $\mathcal{V}_{G_{S}}$, i.e., {\mathcal{E}_{G_{S}}=\mathcal{E}_{G}\cap\big{(}\mathcal{V}_{G_{S}}\times\mathcal{V}_{G_{S}}\big{)}}, the subgraph is said to be induced.

Analogously, we define edge structural roles via edge automorphisms, i.e., bijective mappings from the edge set onto itself, that preserve edge adjacency (two edges are adjacent if they share a common endpoint). In particular, every vertex automorphism $g$ induces an edge automorphism by mapping each edge $(u,v)$ to $(g(u),g(v))$. Note that the edge automorphism group is larger than that of induced automorphisms, but strictly larger only for 3 trivial cases . However, induced automorphisms provide a more natural way to express edge structural roles. In the same way as before, we construct the edge automorphism group, from which we deduce the partition of the edge set in edge orbits $\{O_{H,1}^{E},O_{H,2}^{E},\dots,O_{H,d_{H}}^{E}\}$.

2 Weisfeiler-Leman tests

The Weisfeiler-Leman graph-isomorphism test , also known as naive vertex refinement, 1-WL, or just WL), is a fast heuristic to decide if two graphs are isomorphic or not. The WL test proceeds as follows: every vertex $v$ is initially assigned a colour $c^{0}_{v}$ that is later iteratively refined by aggregating neighbouring information:

where $\Lbag\cdot\Rbag$ denotes a multiset (a set that allows element repetitions) and $\mathcal{N}(v)$ is the neighbourhood of $v$. The WL algorithm terminates when the colours stop changing, and outputs a histogram of colours. Two graphs with different histograms are non-isomorphic; if the histograms are identical, the graphs are possibly, but not necessarily, isomorphic. Note that the neighbour aggregation in the WL test is a form of message passing, and GNNs are the learnable analogue.

A series of works on improving GNN expressivity mimic the higher-order generalisations of WL, known as $k$-WL and $k$-Folklore WL (WL hierarchy) and operate on $k$-tuples of nodes (see Appendix B.1). The $(k+1)$-FWL is strictly stronger than $k$-FWL, $k$-FWL is as strong as $(k+1)$-WL and 2-FWL is strictly stronger than the simple 1-WL test.

Graph Substructure Networks

Graphs consist of nodes (or edges) with repeated structural roles. Thus, it is natural for a neural network to treat them in a similar manner, akin to weight sharing between local patches in CNNs for images or positional encodings in language models for sequential data . Nevertheless, GNNs are usually unaware of the nodes’ different structural roles, since all nodes are treated equally when performing local operations. Despite the initial intuition that the neural network would be able to discover these roles by constructing deeper architectures, it has been shown that GNNs are ill-suited for this purpose and are blind to the existence of structural properties, e.g. triangles or larger cycles .

To this end, we propose to explicitly encode structural roles as part of message passing, in order to capture richer topological properties. Our method draws inspiration from , where it was shown that GNNs become universal when the nodes in the graph are uniquely identified, i.e when they are equipped with different features. However, it is not clear how to choose these identifiers in a permutation equivariant way. Structural roles, when treated as identifiers, although not necessarily unique, are not only permutation equivariant, but also more amenable to generalisation due to their repetition across different graphs. Thus, they can constitute a trade-off between uniqueness and generalisation.

Similarly, we can define edge structural features $\mathbf{x}_{H,i}^{E}(u,v)$ by counting occurrences of edge automorphism orbits:

and the combined edge features ${\mathbf{x}^{E}_{u,v}=[\mathbf{x}_{H_{1}}^{E}(u,v),\ldots,\mathbf{x}_{H_{K}}^{E}(u,v)]}$. An example of vertex and edge structural features is illustrated in Figure 1.

2 Structure-aware message passing

The key building block of our architecture is the graph substructure layer, defined in a general manner as a Message Passing Neural Network (MPNN) , where now the messages from the neighbouring nodes also contain the structural information. In particular, each node $v$ updates its state $\mathbf{h}^{t}_{v}$ by combining its previous state with the aggregated messages:

It is important to note here that contrary to identifier-based GNNs that obtain universality at the expense of permutation equivariance (since the identifiers are arbitrarily chosen with the sole requirement of being unique), GSNs retain this property, hence they are by construction ivariant to isomorphism. This stems from the fact that the process generating our structural identifiers (i.e. subgraph isomorphism) is permutation equivariant itself (proof provided in the Appendix A.1).

3 How powerful are GSNs?

We now turn to the expressive power of GSNs in comparison to MPNNs and the WL tests, a key tool for the theoretical analysis of the expressivity of graph neural networks so far. Since GSN is a generalisation of MPNNs, it is easy to see that it is at least as powerful. Importantly, GSNs have the capacity to learn functions that traditional MPNNs cannot learn. The following observation derives directly from the analysis of the counting abilities of the 1-WL test and its extension to MPNNs (for proofs, see Appendices A.2-A.4).

GSN is strictly more powerful than MPNN and the 1-WL test when one of the following holds:

$H$ is any graph except for star graphs of any size, and structural features are inferred by subgraph matching, i.e. we count all subgraphs $G_{S}\cong H$ for which it holds that $\mathcal{E}_{G_{S}}\subseteq\mathcal{E}_{G}$. Or,

$H$ is any graph except for single edges and single nodes, and structural features are inferred by induced subgraph matching, i.e. we count all subgraphs $G_{S}\cong H$ for which it holds that \mathcal{E}_{G_{S}}=\mathcal{E}_{G}\cap\big{(}\mathcal{V}_{G_{S}}\times\mathcal{V}_{G_{S}}\big{)}.

It is easy to see that GSN model class contains MPNNs, and is thus at least as expressive. We can also show that GSN is at least as expressive as the 1-WL test by repurposing the proof of Theorem 3 in (see Appendix A.2)

Given the first part of the proposition, in order to show that GSNs are strictly more expressive than the 1-WL test, it suffices to show that GSN can distinguish a pair of graphs that 1-WL deems isomorphic. showed that 1-WL, and consequently MPNNs, can count only forests of stars. Thus, if the subgraphs are required to be connected, then they can only be star graphs of any size (note that this contains single nodes and single edges). In addition, showed that 1-WL, and consequently MPNNs, cannot count any connected induced subgraph with 3 or more nodes, i.e. any connected subgraph apart from single nodes and single edges.

If $H$ is a substructure that 1-WL cannot learn to count, i.e. the ones mentioned above, then there is at least one pair of graphs with different number of counts of $H$, that 1-WL deems isomorphic. Thus, by assigning counting features to the nodes/edges of the two graphs based on appearances of $H$, a GSN can obtain different representations for $G_{1}$ and $G_{2}$ by summing up the features. Hence, $G_{1}$, $G_{2}$ are deemed non-isomorphic. An example is depicted in Figure 2 (left), where the two non-isomorphic graphs are distinguishable by GSN via e.g. cycle counting, but not by 1-WL. ∎

Universality. A natural question that emerges is what are the sufficient conditions under which GSN can solve graph isomorphism. This would entail that GSN is a universal approximator of functions defined on graphs . To address this, we can examine whether there exists a specific substructure collection that can completely characterise each graph. As of today, we are not aware of any results in graph theory that can guarantee the reconstruction of a graph from a smaller collection of its subgraphs. However, the Reconstruction Conjecture , states that a graph with size $n\geq 3$ can be reconstructed from its vertex-deleted subgraphs (proven for $n\leq 11$ ). Consequently, (proof in the Appendix A.3):

If the Reconstruction Conjecture holds and the substructure collection $\mathcal{H}$ contains all graphs of size $k=n-1$, then GSN can distinguish all non-isomorphic graphs of size $n$ and is therefore universal.

GSN-v vs GSN-e. We can also examine the expressive power of the two proposed variants. A crucial observation that we make is that for each graph $H$ in the collection, the vertex structural identifiers can be reconstructed by the corresponding edge identifiers. Thus, we can show that for every GSN-v there exists a GSN-e that can simulate the behaviour of the former (proof provided in the Appendix A.4).

For a given subgraph collection $\mathcal{H}$, let $C^{V}$ the set of functions that can be expressed by a GSN-v with arbitrary depth and with, and $C^{E}$ the set of functions that can be expressed by a GSN-e with the same properties. Then, it holds that $C^{E}\supseteq C^{V}$, or in other words GSN-e is at least as expressive as GSN-v.

Comparison with higher-order WL tests. Finally, the expressive power of GSN can be compared to higher-order versions of the WL test. In particular, for each $k$-th order Folklore WL test in the hierarchy, it is known that there exists a family of graphs that will make the test fail. These are known in the literature as $k$-isoregular graphs , and the most well-known example is the Strongly Regular (SR) graph family, for $k=2$ (more details can be found in Appendix B).

Hence, if we can find a substructure collection that allows GSN to distinguish certain pairs from these families, then this guarantees that the corresponding $k$-FWL test is no stronger than GSN. In this work, we identify such counterexamples for the 2-FWL test. Formally:

There exist substructure families with $\mathcal{O}(1)$ size, i.e., independent of the size of the graph $n$, such that 2-FWL is no stronger than GSN.

We provide numerous counterexamples that prove this claim. Figure 2 (right) provides a typical pair of SR graphs that can be distinguished with a 4-clique, while in section 5.1 this is extended to a large-scale study, where other constant size substructures (paths, cycles and cliques) can achieve similar results.

Although it is not clear if there exists a certain substructure collection that results in GSNs that align with the WL hierarchy, we stress that this is not a necessary condition in order to design more powerful GNNs. In particular, despite the increase in expressivity, k-WL tests are not only more computationally involved, but they also process the graph in a non-local fashion. However, locality is presumed to be a strong inductive bias of GNNs and key to their excellent performance in real-world scenarios.

4 How to choose the substructures?

Expressivity. The Reconstruction Conjecture provides a sufficient, albeit impractical condition for universality. This motivates us to analyse the constant size case $k=\mathcal{O}(1)$ for practical scenarios, similar to the argument put forward for hard instances of graph isomorphism (Proposition 3.4).

In particular, one can count only the most discriminative subgraphs, i.e. the ones that can achieve the maximum possible vertex disambiguation, similarly to identifier-based approaches. Whenever these subgraph counts can provide a unique identification of the vertices, then universality will also hold (Corollary 3.1. in ).

We conjecture, that in real-world scenarios the number of subgraphs needed for unique, or near-unique identification, are far fewer than those dictated by Corollary 3.2. This is consistent with our experimental findings, where we observed that certain small substructures such as paths and trees, significantly improve vertex disambiguation, compared to the initial vertex features (see Figure 5 (left) and Table 2 in the appendix). As expected this allows for better fitting of the training data, which validates our claim that GNN expressivity improves.

Generalisation. However, none of the above claims can guarantee good generalisation in unseen data. For example, in Figure 5, we observe that the test set performance does not follow the same trend with train performance when choosing substructures with strong vertex disambiguation. Aiming at better generalisation, it is desirable to make use of substructures for which there is prior knowledge of their importance in certain network distributions and have been observed to be intimately related to various properties. For example, small substructures (graphlets) have been extensively analysed in protein-protein interaction networks , triangles and cliques characterise the structure of ego-nets and social networks in general , simple cycles (rings) are central in molecular distributions, directed and temporal motifs have been shown to explain the working mechanisms of gene regulatory networks, biological neural networks, transportation networks and food webs .

In Figure 5 (right), we showcase the importance of these inductive biases: a cycle-based GSN predicting molecular properties achieves smaller generalisation gap compared to a traditional MPNN, while at the same time generalising better with less training data. Choosing the best substructure collection is still an open problem that does not admit a straightforward solution due to its combinatorial nature. Alternatively, various heuristics can be used, e.g., motif frequencies or feature selection strategies. Answering this question is left for future work.

5 Complexity

The complexity of GSN comprises two parts: precomputation (substructure counting) and training/testing. The key appealing property is that training and inference are linear w.r.t the number of edges, $\mathcal{O}(|\mathcal{E}|)$, as opposed to higher-order methods with $\mathcal{O}(n^{k})$, and with $\mathcal{O}(n^{2})$ training complexity and relational pooling with $\mathcal{O}(n!)$ training complexity in absence of approximations.

The worst-case complexity of subgraph isomorphism of fixed size $k$ is $\mathcal{O}(n^{k})$, by examining all the possible $k$-tuples in the graph. However, for specific types of subgraphs, such as paths and cycles, the problem can be solved even faster (see e.g. ). Approximate counting algorithms are also widely used, especially for counting frequent network motifs , and can provide a considerable speed-up. Furthermore, recent neural approaches provide fast approximate counting.

In our experiments, we performed exact counting using the common isomorphism algorithm VF2 . Although its worst case complexity is $\mathcal{O}(n^{k})$, it scales better in practice, for instance when the candidate subgraph is infrequently matched or when the graphs are sparse, and is also trivially parallelisable. In Figure 3, we show a quantitative analysis of the empirical runtime of the counting algorithm against the worst case, for three different graph distributions: molecules, protein contact maps, social networks. It is easy to see that when the graphs are sparse (for the first two cases) and the number of matches is small, the algorithm is significantly faster than the worst case, while it scales better with the size of the graph $n$. Even, in the case of social networks, where several examples are near-complete graphs, both the runtime and the growth w.r.t both $n$ and $k$ are better than the worst case. Overall, the preprocessing computational burden in most of the cases remains negligible for relatively small and sparse graphs, as it is the case of molecules.

Related Work

WL hierarchy. The seminal results in the theoretical analysis of the expressivity of GNNs and k-GNNs established that traditional message passing-based GNNs are at most as powerful as the 1-WL test. showed that graph isomorphism is equivalent to universal invariant function approximation. Higher-order Invariant Graph Networks (IGNs) have been studied in a series of works , establishing connections with the WL hierarchy, similarly to . The main drawbacks of these methods are the training and inference time complexity and memory requirements of $\mathcal{O}(n^{k})$ and the super-exponential number of parameters (for linear IGNs) making them impractical, as well as their non-local nature making them more prone to overfitting. Finally, also analysed the expressive power of MPNNs and other more powerful variants and provided generalisation bounds.

Unique identifiers. From a different perspective, and showed the connections between GNNs and distributed local algorithms and suggested more powerful alternatives based on either local orderings or unique global identifiers (in the form of random features in ) that make GNNs universal. Similarly, propose to use random colorings in order to uniquely identify the nodes. However, these methods lack a principled permutation equivariant way to choose orderings/identifiers. To date this is an open problem in graph theory called graph canonisation and it is at least as hard as solving graph isomorphism itself. A possible workaround is proposed in , where the authors take into account all possible vertex permutations. However, obviously this quickly becomes intractable ($\mathcal{O}(n!)$) even when considering small-sized graphs.

More expressive permutation equivariant GNNs. Concurrently with our work, other more expressive GNNs have been proposed using equivariant message passing. In , the authors propose to linearly transform each message with a different kernel based on the local isomorphism class of the corresponding edge (similar to our definition of structural roles). However, as also noted by the authors, taking into account all possible local isomorphism classes leads to insufficient weight sharing and hence to overfitting. In contrast, in GSN, usually the substructure collection is small ($\mathtt{\sim}$5-10 graphs) and the substructures are repetitive in the graph distribution, and as a result generalisation improves. Vignac et al. propose a message passing scheme where matrices of order equal to the size of the graph are propagated instead of vectors. This can be perceived as a practical unique identification scheme, but the neural network complexity becomes quadratic in the number of nodes. Finally, and enhance the aggregation function with distance encodings (a strategy more relevant for vertex-level tasks) and graph eigenvectors respectively as alternative symmetry breaking mechanisms. In the experimental section, GSN is compared against these methods in real-world scenarios.

Quantifying expressivity. Solely quantifying the expressive power of GNNs in terms of their ability to distinguish non-isomorphic graphs does not provide the necessary granularity: even the 1-WL test can distinguish almost all (in the probabilistic sense) non-isomorphic graphs . As a result, there have been several efforts to analyse the power of $k$-WL tests in comparison to other graph invariants , while recently approached GNN expressivity by studying their ability to count substructures.

2 Substructures in Complex Networks.

The idea of analysing complex networks based on small-scale topological characteristics dates back to the 1970’s and the notion of triad census for directed graphs . The seminal paper of coined the term network motifs as over-represented subgraph patterns that were shown to characterise certain functional properties of complex networks in systems biology. The closely related concept of graphlets , different from motifs in being induced subgraphs, has been used to analyse the distribution of real-world networks and as a topological signature for network similarity. Our work is similar in spirit with the graphlet degree vector (GDV) , a node-wise descriptor based on graphlet counting.

Substructures have been also used in the context of ML. In particular, subgraph patterns have been used to define Graph Kernels (GKs) , with the most prominent being the graphlet kernel . Motif-based node embeddings and diffusion operators that employ adjacency matrices weighted according to motif occurrences, have recently been proposed for graph representation learning. Our formulation provides a unifying framework for these methods and it is the first to analyse their expressive power. Finally, GNNs that operate in larger induced neighbourhoods or higher-order paths have prohibitive complexity since the size of these neighbourhoods typically grows exponentially.

Experimental Evaluation

In the following section we evaluate GSN in comparison to the state-of-the-art in a variety of datasets from different domains. We are interested in practical scenarios where the collection of subgraphs, as well as their size, are kept small. Depending on the dataset domain we experimented with typical substructure families (cycles, paths and cliques) and maximum substructure size $k$ (note that for each setting, our substructure collection consists of all the substructures of the family with size $\leq k$). We also experimented with both graphlets and motifs and observed similar performance in most cases. To showcase that structural features can be used as an off-the-shelf strategy to boost GNN performance, we usually choose a base message passing architecture and minimally modify it into a GSN. Unless otherwise stated, the base architecture is a general-purpose MPNN with MLPs used in the message and update functions. Additional implementation details can be found in the Appendix C.

We tested the ability of GSNs to decide if two graphs are non-isomorphic on a collection of Strongly Regular graphs of size up to 35 nodes, attempting to disambiguate pairs with the same number of nodes (for different sizes the problem becomes trivial). As we are only interested in the bias of the architecture itself, we use GSN with random weights to compute graph representations. Two graphs are deemed isomorphic if the Euclidean distance of their representations is smaller than a predefined threshold $\epsilon$. Figure 4 shows the failure percentage of our isomorphism test when using different graphlet substructures (cycles, paths, and cliques) of varying size $k$. Interestingly, the number of failure cases of GSN decreases rapidly as we increase $k$; cycles and paths of maximum length $k=6$ are enough to tell apart all the graphs in the dataset. Note that the performance of cliques saturates, possibly because the largest clique in our dataset has 5 nodes. Observe also the discrepancy between GSN-v and GSN-e. In particular, vertex-wise counts do not manage to distinguish all graphs, although missing only a few instances, which is in accordance with Theorem 3.3. Finally, 1-WL and 2-FWL equivalent models demonstrate 100% failure, as expected from theory.

2 TUD Graph Classification Benchmarks

We evaluate GSN on datasets from the classical TUD benchmarks. We use seven datasets from the domains of bioinformatics and computational social science and compare against various GNNs and Graph Kernels. The base architecture that we used is GIN . We follow the same evaluation protocol of , performing 10-fold cross-validation and then reporting the performance at the epoch with the best average accuracy across the 10 folds. Table 1 lists all the methods evaluated with the split of . We select our model by tuning architecture and optimisation hyperparameters and substructure related parameters, that is: (i) $k$, (ii) motifs against graphlets. Following domain evidence we choose the following substructure families: cycles for molecules, cliques for social networks. Best performing substructures both for GSN-e and GSN-v are reported. As can be seen, our model obtains state-of-the-art performance in most of the datasets, with a considerable margin against the main GNN baselines in some cases.

3 ZINC Molecular graphs

We evaluate GSN on the task of regressing the “penalized water-octanol partition coefficient - logP” (see for details) of molecules from the ZINC database . We use structural features obtained with $k$-cycle counting and report the result of the best performing substructure w.r.t. the validation set.

As dictated by the evaluation protocol of , the total number of parameters of the model is approximately 100K, which is achieved by selecting an appropriate network width.A larger version of GSN using 500K parameters attains 0.101 ± 0.010 test MAE. The data split is obtained from and the evaluation metric is the Mean Absolute Error (MAE). We compare against a variety of baselines, ranging from traditional message passing NNs to recent more expressive architectures and a molecular-specific one which is based on the junction tree molecular decomposition . Wherever possible, we compare two variants, one that does not take edge features into account and one that does. In both cases, GSN achieves state-of-the-art results outperforming all the baseline architectures.

4 OGB-MOLHIV

We use ogbg-molhiv from the Open Graph Benchmark - OGB - as a graph-level binary classification task, where the aim is to predict if a molecule inhibits HIV replication or not. We obtain baseline results from a plethora of different methodsthe most representative ones from the OGB public leaderboard: https://ogb.stanford.edu/docs/leader_graphprop/##ogbg-molhiv, ranging from algorithms specific to molecular graphs to general purpose architectures. Following the same rationale as in the previous experiments, we choose a base architecture and modify it into a GSN variant by introducing structural features in the aggregation function (cycle counts, similar to other molecular datasets). Here we use the following two base architectures: (a) GIN-VN, a variation of GIN that allows for edge features and is extended with a virtual node, i.e. a node connected to every node in the graph. (b) Directional Graph Networks (DGN), a GNN that propagates messages in an anisotropic manner, based on a predefined graph vector field. Observe that the vector field is an alternative way to break local symmetries. The authors of DGN use vector fields defined by the eigenvectors of the graph, while in our case the vector field is defined by graph substructures. More information can be found in the supplementary material.

Using the evaluator provided by the authors, we report the ROC-AUC metric at the epoch with the best validation performance (substructures are also chosen based on the validation set). By examining the results in Table 3 the following observations can be made, (a) general purpose GNNs benefit from symmetry breaking mechanisms, either in the form of eigenvectors or in the form of substructures. (b) Cyclical substructures are a good inductive bias when learning on molecules (e.g. P-WL is a graph kernel based on topological features that contain information of graph cycles, similar to GSN). (c) Further evidence for that is provided by observing the performance of molecular fingerprints methods. In specific, a method based on the extended-connectivity fingerprints , which mainly focuses on the structure of the molecule reports 0.8060 ± 0.0010 test performance, while one that additionally uses the MACCS fingerprints , which mainly encode the presence of certain functional groups (i.e. both structure and attributes), reports 0.8232 ± 0.0047 test performance. These methods although not directly comparable to ours, currently achieve the best results in this dataset, thus this further motivates research on making GNNs structure-aware.

5 Ablation Studies

In Figure 5 (left), we compare the training and test error for different substructure families (cycles, paths and non-isomorphic trees – for each experiment we use all the substructures of size $\leq k$ in the family). Additionally, we measure the “uniqueness” of the identifiers each substructure yields as follows: for each graph $G$ in the dataset, we measure the number of unique vertex features $u_{G}$ (input vertex features concatenated with vertex structural identifiers for GSN-v). Then, we sum them up over the entire training set and divide by the total number of nodes, yielding the disambiguation score $\delta=\frac{\sum_{G}u_{G}}{\sum_{G}|V_{G}|}$. The disambiguation scores for the different types of substructures are illustrated as horizontal bars in Figure 5 (the exact values can be found in Appendix C.2, Table 2).

A first thing to notice is that the training error is tightly related to the disambiguation score. As identifiers become more discriminative, the model gains expressive power. On the other hand, the test error is not guaranteed to decrease when the identifiers become more discriminative. For example, although cycles have smaller disambiguation scores, they manage to generalise much better than the other substructures, the performance of which is similar to the baseline architecture (MPNN with MLPs). This is also observed when comparing against (MPNN-r method in Table 3), where, akin to unique identifiers, random features are used to strengthen the expressivity of GNN architectures. This approach also fails to improve the baseline architecture in terms of the performance in the test set. This validates our intuition that unique identifiers can be hard to generalise when chosen in a non-permutation equivariant way and motivates once more the importance of choosing the identifiers not only based on their discriminative power, but also in a way that allows incorporating the appropriate inductive biases. Finally, we observe a substantial jump in performance when using GSN with cycles of size $k\geq 6$. This is not surprising, as cyclical patterns of such sizes (e.g. aromatic rings) are very common in organic molecules.

5.2 Generalisation

We repeat the experimental evaluation on ZINC using different fractions of the training set and compare the vanilla MPNN model against GSN. In Figure 5 (right), we plot the training and test errors of both methods. Regarding the training error, GSN consistently performs better, following our theoretical analysis on its expressive power. More importantly, GSN manages to generalise much better even with a small fraction of the training dataset. Observe that GSN requires only 20% of the samples to achieve approximately the same test error that MPNN achieves when trained on the entire training set.

5.3 Structural Features & Message Passing:

We perform an ablation study on the abilities of the structural features to predict the task at hand, when given as input to a graph-agnostic network. In particular, we compare our best performing GSN with a DeepSets model that treats the input features and the structural identifiers as a set. For fairness of evaluation the same hyperparameter search is performed for both models (see Appendix C.5. Interestingly, as we show in Table 4, our baseline attains particularly strong performance across a variety of datasets and often outperforms other traditional message passing baselines. This demonstrates the benefits of these additional features and motivates their introduction in GNNs, which are unable to compute them. As expected, we observe that applying message passing on top of these features, brings performance improvements in the vast majority of the cases, sometimes considerably, as in the ZINC dataset.

Conclusion

In this paper, we propose a novel way to design structure-aware graph neural networks. Motivated by the limitations of traditional GNNs to capture important topological properties of the graph, we formulate a message passing scheme enhanced with structural features that are extracted by subgraph isomorphism. We show both theoretically and empirically that our construction leads to improved expressive power and attains state-of-the-art performance in real-world scenarios. In future work, we will further explore the expressivity of GSNs as an alternative to the $k$-WL tests, as well as their generalisation capabilities. Another important direction is to infer prominent substructures directly from the data and explore the ability of graph neural networks to compose substructures.

Acknowledgements

This research was partially supported by the ERC Consolidator Grant No. 724228 - LEMAN (GB and MB). The work of GB is partially funded by a PhD scholarship from the Department of Computing, Imperial College London. SZ acknowledges support from the EPSRC Fellowship DEFORM: Large Scale Shape Analysis of Deformable Models of Humans (EP/S010203/1) and a Google Faculty award. MB acknowledges support from Google Faculty awards and the Royal Society Wolfson Research Merit award.

References

Appendix A Deferred Proofs

It is obvious that subgraph isomorphism is invariant to the ordering of the vertices, i.e. we will always obtain the same matching between subgraphs $G_{S}$ and graphs $H$ in the subgraph collection. Thus, each vertex $v$ (edge $(v,u)$) in the graph will be assigned the same structural identifiers $\mathbf{x}^{V}_{v}$ ($\mathbf{x}^{E}_{v,u}$) regarless of the vertex ordering, and $S^{V}$ and $S^{E}$ are permutation equivariant, i.e., for any permutation matrix $\mathbf{P}\in\{0,1\}^{n\times n}$ it holds that

where in the last step we used the permutation equivariant property of MPNNs. Similarly, we can show that a GSN-e layer is permutation equivariant, since it can be expressed as a traditional MPNN layer by replacing the edge features with the concatenation of the original edge features and the edge structural identifiers $f(\mathbf{A},\mathbf{H},\mathbf{E})=g(\mathbf{A},\mathbf{H},[\mathbf{E};S^{E}(A)])$.

Overall, a GSN network is permutation equivariant as composition of permutation equivariant functions, or permutation invariant when composed with a permutation invariant layer at the end, i.e. the READOUT function.

A.2 Proof of Theorem 3.1: GSN is at least as powerful as the 1-WL test

To show that GSN it is at least as expressive as the 1-WL test, we will repurpose the proof of Theorem 3 in and demand the injectivity of the update function (w.r.t. both the hidden state $\mathbf{h}^{t}_{v}$ and the message $\mathbf{m}^{t+1}_{v}$), and the injectivity of the aggregation w.r.t. the multiset of the hidden states of the neighbours $\Lbag\mathbf{h}^{t}_{u}\Rbag_{u\in\mathcal{N}(v)}$. It suffices then to show that if injectivity is preserved then GSNs are at least as powerful as the 1-WL.

We will show the above statement for vertex-labelled graphs, since traditionally the 1-WL test does not take into account edge labels.if one considers a simple 1-WL extension that concatenates edge labels to neighbour colours, then the same proof applies. We can rephrase the statement as follows: If GSN deems two graphs $G_{1}$, $G_{2}$ as isomorphic, then also 1-WL deems them isomorphic. Given that the graph-level representation is extracted by a readout function that receives the multiset of the vertex colours in its input (i.e. the graph-level representation is the vertex colour histogram at some iteration $t$), then it suffices to show that if for the two graphs the multiset of the vertex colours that GSN infers is the same, then also 1-WL will infer the same multiset for the two graphs.

Consider the case where the two multisets that GSN extracts are the same: i.e. ${\Lbag\mathbf{h}^{t}_{v}\Rbag_{v\in\mathcal{V}_{G_{1}}}=\Lbag\mathbf{h}^{t}_{u}\Rbag_{u\in\mathcal{V}_{G_{2}}}}$. Then both multisets contain the same distinct colour/hidden representations with the exact same multiplicity. Thus, it further suffices to show that if two vertices $v,u$ (that may belong to the same or to different graphs) have the same GSN hidden representations $\mathbf{h}^{t}_{v}=\mathbf{h}^{t}_{u}$ at any iteration $t$, then they will also have the same colours $c^{t}_{v}=c^{t}_{u}$, extracted by 1-WL. Intuitively, this means that GSN creates a partition of the vertices of each graph that is at least as fine-grained as the one created by 1-WL. We prove by induction (similarly to ) that GSN model class contains a model where this holds (w.l.o.g. we show that for GSN-v; same proof applies to GSN-e).

$\mathbf{h}^{t-1}_{v}=\mathbf{h}^{t-1}_{u}$, which from the induction hypothesis implies that $c^{t-1}_{v}=c^{t-1}_{u}.$

$\mathbf{m}^{t}_{v}=\mathbf{m}^{t}_{u}$, where the message function is defined as in Eq. (5) of the main paper. Additionally here we require $M^{t}$ to be injective w.r.t. the multiset of the hidden representations of the neighbours:Lemma 5 from states that such a function always exists assuming that the elements of the multiset originate from a countable domain

From the induction hypothesis we know that ${\mathbf{h}^{t-1}_{w}=\mathbf{h}^{t-1}_{z}}$ implies that ${c^{t-1}_{w}=c^{t-1}_{z}}$ for any ${w\in\mathcal{N}_{v},z\in\mathcal{N}_{u}}$, thus ${\Lbag c^{t-1}_{w}\Rbag_{w\in\mathcal{N}_{v}}=\Lbag c^{t-1}_{z}\Rbag_{z\in\mathcal{N}_{u}}}$.

A.3 Proof of Corollary 3.2

In order to prove the universality of GSN, we will show that when the substructure collection contains all graphs of size $n-1$, then there exists a parametrisation of GSN that can infer the isomorphism classes of all vertex-deleted subgraphs of the graph $G$ (the deck of $G$).

The reconstruction conjectures states that two graphs with at least three vertices are isomorphic if and only if they have the same deck. Thus, the deck is sufficient to distinguish all non-isomorphic graphs. The deck can be defined as follows:

Let $\mathcal{H}=\{H_{1},H_{2},\cdots,H_{K}\}$ the set of all possible graphs of size $n-1$. The vertex-deleted subgraphs of $G$ are by definition all the induced subgraphs of $G$ with size $n-1$, which we denote as

Then, the deck $D(G)$ can be defined as a vector of size $|\mathcal{H}|$, where at the $j$-th dimension

where $\mathds{1}[\cdot]$ the indicator function. The structural feature $x_{H_{j},i}^{V}(v)$ for each substructure $H_{j}$ and orbit $i$ are computed as follows:

where $f_{G_{S}}(\cdot)=\emptyset$ if $G_{S}\not\simeq H_{j}$, otherwise it is the bijective mapping function. The deck can be inferred as follows:

where we used that $\sum_{i=1}^{d_{H_{j}}}\mathds{1}[f_{G_{S}}(v)\in O_{H_{j},i}^{V}]$ = 1, since each vertex can be mapped to a single orbit only. Thus, the deck can be inferred by a simple GSN-v parametrisation (a linear layer with depth equal to $|\mathcal{H}|$ that performs orbit-wise summation and division by the constant $n-1$ for each vertex separately, followed by a sum readout). Since GSN-v can be inferred by GSN-e (Theorem 3.3), then GSN-e is also universal. ∎

A.4 Proof of Theorem 3.3

Without loss of generality we will show Theorem 3.3 for the case of a single substructure $H$. In order to show that GSN-e can express GSN-v, we will first prove the following: the vertex identifier of a vertex $v$ can be inferred by the edge identifiers of its incident edges.

For brevity we will use the following notation: vertex orbits are indexed as follows $O^{V}_{1},O^{V}_{2},\ldots O^{V}_{d_{V}}$ and edge orbits $O_{1,1}^{E},O_{1,2}^{E},\ldots O_{d_{V},d_{V}}^{E}$ with O_{i,j}^{E}=\big{(}O^{V}_{i},O^{V}_{j}\big{)}. Structural features are denoted accordingly: $x_{i}^{V}(v)$ and $x_{ij}^{E}(v,u)$

Let’s assume that there exists only one matched subgraph $G_{S}\simeq H$ and the bijection between $\mathcal{V}_{G_{S}}$ and $\mathcal{H}$ is denoted as $f$. Then, for an abitrary vertex $v$ and a vertex orbit $O^{V}_{i}$, one of the following holds:

$v\not\in\mathcal{V}_{G_{S}}$. Then $x_{i}^{V}(v)=0$ and $x_{ij}^{E}(v,u)=0,\ \forall u\in\mathcal{N}(v)$,

Note that since $M$ is a universal multiset function approximator, then there exist parameters of $M$ with which the above function can be computed. The next $L$ layers of GSN-e can implement a traditional MPNN where now the input vertex features are $[\mathbf{h}^{t}_{v};\mathbf{x}^{V}_{v}]$ (which is exactly the formulation of GSN-v) and this concludes the proof.

Appendix B Comparison with higher-order Weisfeiler-Leman tests

Following the terminology introduced in , we describe the so-called Folklore WL family ($k$-FWL). Note that, in the majority of papers on GNN expressivity another family of WL tests is discussed, under the terminology $k$-WL with expressive power equal to $(k-1)$-FWL. In contrast, in most graph theory papers on graph isomorphism the $k$-WL term is used to describe the algorithms referred to as $k$-FWL in GNN papers. Here, we follow the $k$-FWL convention to align with the work mostly related to ours.

The $k$-FWL operates on $k$-tuples of vertices ${\mathbf{v}=(v_{1},v_{2},\dots,v_{k})}$ to which an initial colour $c^{0}_{\mathbf{v}}$ is assigned based on their isomorphism types (see section B.2), which can loosely be thought of as a generalisation of isomorphism that also preserves the ordering of the vertices in the tuple. Then, at each iteration the colour is refined as follows:

where $\mathbf{v}_{u,j}=(v_{1},v_{2},\dots,v_{j-1},u,v_{j+1},\dots,v_{k})$.

The multiset \Lbag\big{(}c^{t}_{\mathbf{v}_{u,1}},c^{t}_{\mathbf{v}_{u,2}},\dots,c^{t}_{\mathbf{v}_{u,k}}\big{)}\Rbag_{u\in\mathcal{V}} can be perceived as a form of generalised neighbourhood. Observe that all possible tuples in the graph store information necessary for the updates, thus each $k$-tuple receives information from the entire graph, contrary to the local nature of the 1-WL test.

B.2 Why does 2-FWL fail on strongly regular graphs?

We first formally define what an isomorphism type is and what are the properties of the SR family:

Two k-tuples $\mathbf{v}^{a}=\{v_{1}^{a},v_{2}^{a},\dots,v_{k}^{a}\}$, $\mathbf{v}^{b}=\{v_{1}^{b},v_{2}^{b},\dots,v_{k}^{b}\}$ will have the same isomorphism type iff:

$\forall\,i,j\in\{0,1,\dots,k\},\quad v_{i}^{a}=v_{j}^{a}\Leftrightarrow v_{i}^{b}=v_{j}^{b}$

$\forall\,i,j\in\{0,1,\dots,k\},\quad v_{i}^{a}\sim v_{j}^{a}\Leftrightarrow v_{i}^{b}\sim v_{j}^{b}$, where $\sim$ means that the vertices are adjacent.

Note that this is a stronger condition than isomorphism, since the mapping between the vertices of the two tuples needs to preserve order. In case the graph is employed with edge and vertex features, these need to be preserved as well (see ) for the extended case).

A SR($n$,$d$,$\lambda$,$\mu$)-graph is a regular graph with $n$ vertices and degree $d$, where every two adjacent vertices have always $\lambda$ mutual neighbours, while every two non-adjacent vertices have always $\mu$ mutual neighbours.

Now we can proceed to the details of the proof. For the 2-FWL test, when working with simple undirected graphs without self-loops, we have the following 2-tuple isomorphism types:

$\mathbf{v}=\{v_{1},v_{1}\}$: vertex type. Mapped to the colour $c^{0}=c_{\alpha}$

$\mathbf{v}=\{v_{1},v_{2}\}$ and $v_{1}\not\sim v_{2}$: non-edge type. Mapped to the colour $c^{0}=c_{\beta}$

$\mathbf{v}=\{v_{1},v_{2}\}$ and $v_{1}\sim v_{2}$: edge type. Mapped to the colour $c^{0}=c_{\gamma}$

For each 2-tuple $\mathbf{v}=\{v_{1},v_{2}\}$, a generalised “neighbour” is the following tuple: {(\mathbf{v}_{u,1},\mathbf{v}_{u,2})=\big{(}(u,v_{2}),(v_{1},u)\big{)}}, where $u$ is an arbitrary vertex in the graph.

Now, let us consider a strongly regular graph SR($n$,$d$,$\lambda$,$\mu$). We have the following cases:

generalised neighbour of a vertex type tuple: {(\mathbf{v}_{u,1},\mathbf{v}_{u,2})=\big{(}(u,v_{1}),(v_{1},u)\big{)}}. The corresponding neighbour colour tuples are:

$(c_{\beta},c_{\beta})$ if $v_{1}\not\sim u$ ,

$(c_{\gamma},c_{\gamma})$ if $v_{1}\sim u$.

generalised neighbour of a non-edge type tuple: {(\mathbf{v}_{u,1},\mathbf{v}_{u,2})=\big{(}(u,v_{2}),(v_{1},u)\big{)}}. The corresponding neighbour colour tuples are:

$(c_{\gamma},c_{\beta})$ if $v_{2}\sim u$ and $v_{1}\not\sim u$,

$(c_{\beta},c_{\gamma})$ if $v_{1}\sim u$ and $v_{2}\not\sim u$,

$(c_{\beta},c_{\beta})$ if $v_{1}\not\sim u$ and $v_{2}\not\sim u$,

$(c_{\gamma},c_{\gamma})$ if $v_{1}\sim u$ and $v_{2}\sim u$.

generalised neighbour of an edge type tuple:

$(c_{\gamma},c_{\beta})$ if $v_{2}\sim u$ and $v_{1}\not\sim u$,

$(c_{\beta},c_{\gamma})$ if $v_{1}\sim u$ and $v_{2}\not\sim u$,

$(c_{\beta},c_{\beta})$ if $v_{1}\not\sim u$ and $v_{2}\not\sim u$,

$(c_{\gamma},c_{\gamma})$ if $v_{1}\sim u$ and $v_{2}\sim u$.

From the analysis above, it is clear that all 2-tuples in the graph of the same initial type are assigned the same colour in the 1st iteration of 2-FWL. In other words, the vertices cannot be further partitioned, so the algorithm terminates. Therefore, if two SR graphs have the same parameters $n$,$d$,$\lambda$,$\mu$ then 2-FWL will yield the same colour distribution and thus the graphs will be deemed isomorphic.

Appendix C Experimental Settings - Additional Details

In this section, we provide additional implementation details of our experiments. All experiments were performed on a server equipped with 8 Tesla V100 16 GB GPUs, except for the Collab dataset where a Tesla V100 GPU with 32 GB RAM was used due to larger memory requirements. Experimental tracking and hyperparameter optimisation were done via the Weights & Biases platform (wandb) . Our implementation is based on native PyTorch sparse operations in order to ensure complete reproducibility of the results. PyTorch Geometric was used for additional operations (such as preprocessing and data loading).

In each one of the different experiments we aim to show that structural identifiers can be used off-the-shelf and independently of the architecture. At the same time we aim to suppress the effect of other confounding factors in the model performance, thus wherever possible we build our model on top of a baseline architecture. For more details, please see the relevant subsections. Interestingly, we observed that in most of the cases it was sufficient to replace only the first layer of the baseline architecture with a GSN layer, in order to obtain a boost in performance.

Throughout the experimental evaluation the structural identifiers $\mathbf{x}^{V}_{v}$ and $\mathbf{x}^{E}_{u,v}$ are one-hot encoded, by taking into account the unique count values present in the dataset. Other more sophisticated methods can be used, e.g. transformation to continuous features via a normalisation scheme or binning. However, we found that the number of unique values in our datasets were usually relatively small (which is a good indication of recurrent structural roles) and thus such methods were not necessary.

For the Strongly Regular graphs dataset (available from http://users.cecs.anu.edu.au/~bdm/data/graphs.html) we use all the available families of graphs with size of at most 35 vertices:

The total number of non-isomorphic pairs of the same size is $\approx 7*10^{7}$. We used a simple 2-layer architecture with width 64. The message aggregation was performed as in the general formulation of Eq. (5) of the main paper, where the update and the message functions are MLPs. The prediction is inferred by applying a sum readout function in the last layer and then passing the output through a MLP. Regarding the substructures, we use graphlet counting, as certain motifs (e.g. cycles of length up to 7) are known to be unable to distinguish strongly regular graphs (since they can be counted by the 2-FWL ).

Given the adversities that strongly regular graphs pose in graph isomorphism testing, it would be interesting to see how this method can perform in other categories of hard instances, such as the classical CFI counter-examples for k-WL proposed in , and explore further its expressive power and combinatorial properties. We leave this direction to future work.

C.2 TUD Graph Classification Benchmarks

For this family of experiments, due to the usually small size of the datasets, we choose a parameter-efficient architecture, in order to reduce the risk of overfitting. In particular, we follow the simple GIN architecture and we concatenate structural identifiers to vertex or edge features depending on the variant. Then for GSN-v, the hidden representation is updated as follows:

where $\mathbf{x}^{E}_{v,v}$ is a dummy variable (also one-hot encoded) used to distinguish self-loops from edges. Empirically, we did not find training the $\epsilon$ parameter used in GIN to make a difference.

We implement an architecture similar to GIN , i.e. message passing layers: 4 , jumping knowledge from all the layers (including the input), transformation of each intermediate graph-level representation: linear layer, readout: sum for biological and mean for social networks. Vertex features are one-hot encodings of the categorical vertex labels. Similarly to the baseline, the hyperparameters search space is the following: batch size in {32, 128} (except for Collab where only 32 was searched due to GPU memory limits), dropout in {0,0.5}, network width in {16,32} for biological networks, 64 for social networks, learning rate in {0.01, 0.001}, decay rate in {0.5,0.9} and decay steps in {10,50} (number of epochs after which the learning rate is reduced by multiplying with the decay rate). For social networks, since they are not attributed graphs, we also experimented with using the degree as a vertex feature, but in most cases the structural identifiers were sufficient.

Model selection is done in two stages. First, we choose a substructure that we perceive as promising based on indications from the specific domain: triangles for social networks and Proteins, and 6-cycles (motifs) for molecules. Under this setting we tune model hyperparameters for a GSN-e model. Then, we extend our search to the parameters related to the substructure collection: i.e. the maximum size $k$ and motifs vs graphlets. In all the molecular datasets we search cycles with $k=3,\ldots,12$, except for NCI1, where we also consider larger sizes due to the presence of large rings in the dataset (macrocycles ). For social networks, we searched cliques with $k=3,4,5$. In Table 5 we report the hyperparameters chosen by our model selection procedure, including the best performing substructures.

The seven datasetsmore details on the description of the datasets and the corresponding tasks can be found at . we chose are the intersection of the datasets used by the authors of our main baselines: the Graph Isomorphism Network (GIN) , a simple, yet powerful GNN with expressive power equal to the 1-WL test, and the Provably Powerful Graph Network (PPGN) , a polynomial alternative to the Invariant Graph Network , that increases its expressive power to match the 2-FWL. We also compare our results to other GNNs as well as Graph Kernel approaches. Our main baseline from the GK family is the Graph Neural Tangent Kernel (GNTK) , which is a kernel obtained from a GNN of infinite width. This operates in the Neural Tangent Kernel regime .

C.3 Graph Regression on ZINC

Experimental Details: The ZINC dataset includes 12k molecular graphs of which 10k form the training set and the remaining 2k are equally split between validation and test (splits obtained from https://github.com/graphdeeplearning/benchmarking-gnns). Molecule sizes range from 9 to 37 vertices/atoms. Vertex features encode the type of atoms and edge features the chemical bonds between them. Again, here vertex and edge features are one-hot encoded.

Our MPNN baseline model updates vertex representations as follows: $\mathbf{h}^{t+1}(v)=\text{MLP}^{t+1}(\mathbf{h}^{t}_{v},\mathbf{m}^{t+1}_{v})$, \mathbf{m}^{t+1}_{v}=\sum_{u\in\mathcal{N}(v)}\text{MLP}^{t}\big{(}\mathbf{h}^{t}_{v},\mathbf{h}^{t}_{u},\mathbf{e}_{u,v}\big{)}. Our instantiation of GSN is a simple extension where structural identifiers are also given as input to the message MLP.

Following the same rationale as before, the network configuration is minimally modified w.r.t. the baselines provided in , while here no hyperparameter tuning is done and we use the default ones provided by the authors. In particular, the parameters are the following: message passing layers: 4, transformation of the output of the last layer: MLP, readout: sum, batch size: 128, dropout: 0.0, network width: 128, learning rate: 0.001. The learning rate is reduced by 0.5 (decay rate) after 5 epochs (decay rate patience) without improvement in the validation loss. Training is stopped when the learning rate reaches the minimum learning rate value of $10^{-5}$. Validation and test metrics are inferred using the model at the last training epoch.

We select our best performing substructure related parameters based on the performance in the validation set in the last epoch. We search cycles with $k$ $=3,\ldots,10$, graphlets vs motifs, and GSN-v vs GSN-e. The chosen hyperparameters for GSN are: GSN-e, cycle graphlets of 10 vertices and for GSN-EF: GSN-v, cycle motifs of 8 vertices. Once the model is chosen, we repeat the experiment 10 times with different seeds and report the mean and standard deviation of the test MAE in the last epoch.

Disambiguation Scores: In Table 6, we provide the disambiguation scores $\delta$ as defined in section 5.5.1 of the main paper for different types of substructures. These are computed based on vertex structural identifiers (GSN-v).

C.4 Graph Classification on ogbg-molhiv

The ogbg-molhiv dataset contains $\approx$ 41K graphs, with 25.5 vertices and 27.5 edges on average. As most molecular graphs, the average degree is small (2.2) and they exhibit a tree-like structure (average clustering coefficient 0.002). The average diameter is 12 (more details in ). Below we describe how we extend the two base architectures:

GIN+VN. We follow the design choices of the authors of and extend their architectures to include structural identifiers. Initial vertex and edge features are multi-hot encodings passed through linear layers that project them in the same embedding space, i.e. ${\mathbf{h}^{0}_{v}=\mathbf{W}^{0}_{h}\cdot\mathbf{h}^{in}_{v}}$, ${\mathbf{e}^{t}_{v,u}=\mathbf{W}^{t}_{e}\cdot\mathbf{e}^{in}_{u,v}}$. The baseline model is a modification of GIN that allows for edge features: for each neighbour, the hidden representation is added to an embedding of its associated edge feature. Then the result is passed through a ReLU non-linearity which produces the neighbour’s message. Formally, the aggregation is as follows:

In order to allow global information to be broadcasted to the vertices, a virtual node takes part in the message passing. The virtual node representation, denoted as $\mathbf{G}^{t}$, is initialised as a zero vector $\mathbf{G}^{0}$ and then Message Passing becomes:

We modify this model, as follows: first the substructure counts are embedded into the same embedding space as the rest of the features. Then, for GSN-v, they are added to the corresponding vertex embeddings: ${\mathbf{\acute{h}}^{t}_{v}=\mathbf{h}^{t}_{v}+\mathbf{W}^{t}_{V}\cdot\mathbf{x}^{V}_{v}}$, or for GSN-e, they are added to the edge embeddings $\mathbf{\acute{e}}^{t}_{v,u}=\mathbf{e}^{t}_{v,u}+\mathbf{W}^{t}_{E}\cdot\mathbf{x}^{E}_{u,v}$.

DGN + substructures. We use the directional average operator as defined in :

where $\alpha^{i}_{v,u}$ are weighting average coefficients. In our case, each orbit $i$ induces a separate set of averaging coefficients. For example, for GSN-e $\alpha_{v,u}=\frac{|x^{E}_{v,u}|}{\epsilon\ +\ \sum_{u\in\mathcal{N}(v)}|x^{E}_{v,u}|}$, where $x^{E}_{v,u}$ denotes edge-wise substructure counts (the index of the orbit $i$ was dropped to simplify notation). Similarly, for GSN-v, $\alpha_{v,u}=\frac{|x^{V}_{v}-x^{V}_{u}|}{\epsilon\ +\ \sum_{u\in\mathcal{N}(v)}|x^{V}_{v}-x^{V}_{u}|}$. Subsequently, the vertex representation is updated as follows: ${\mathbf{h}^{t+1}_{v}=\text{MLP}^{t+1}(\mathbf{m}^{t+1}_{v})}$. Observe that this model is simpler than the aforementioned, in terms of both its parameter count and its expressive power. Since the MOLHIV dataset poses a significant challenge w.r.t. generalisation (the data splits reflect different molecular distributions), architectures biased towards simpler solutions usually perform better, sincey the mitigate the risk of overfitting.

In both cases we use the the same hyperparameters as the ones provided by the authors, and only select the substructure related parameters based on the highest validation ROC-AUC (choosing the best scoring epoch as in ). We search cycles with $k=3,\ldots,12$, graphlets vs motifs, and GSN-v vs GSN-e. The chosen hyperparameters are: GSN-e, cycle graphlets of 6 vertices. We repeat the experiment 10 times with different seeds and report the mean and standard deviation of the train, validation and test ROC-AUC, again by choosing the best scoring epoch w.r.t the validation set.

C.5 Structural Features & Message Passing

The baseline architecture treats the input vertex and edge features, along with the structural identifiers, as a set. In particular, we consider each graph as a set of independent edges $(v,u)$ endowed with the features of the endpoint vertices $\mathbf{h}_{v},\mathbf{h}_{u}$, the structural identifiers $\mathbf{x}^{V}_{v},\mathbf{x}^{V}_{u}$ and the edge features $\mathbf{e}(v,u)$, and we implement a DeepSets universal set function approximator to learn a prediction function:

with $\mathcal{E}_{G}$ the edge set of the graph and $\rho,\phi$ MLPs. This baseline is naturally extended to the case where we consider edge structural identifiers by replacing $(\mathbf{x}^{V}_{v},\mathbf{x}^{V}_{u})$ with $\mathbf{x}^{E}_{v,u}$. For fairness of evaluation, we follow the exact same parameter tuning procedure as the one we followed for our GSN models for each benchmark, i.e. for the TUD datasets we first tune network and optimisation hyperaparameters (network width was set to be either equal to the ones we tuned for GSN, or such that the absolute number of learnable parameters was equal to those used by GSN; depth of the MLPs was set to 2) and subsequently we choose the substructure related parameters based on the evaluation protocol of . For ZINC and ogbg-molhiv we perform only substructure selection, based on the performance on the validation set. Using the same widths as in GSN leads to smaller baseline models w.r.t the absolute number of parameters, and we interestingly observed this to lead to particularly strong performance in some cases, especially Proteins and MUTAG, where our DeepSets implementation attains state-of-art results. This finding motivated us to explore ‘smaller’ GSNs (with either reduced layer width or a single message passing layer). These GSN variants exhibited a similar trend, i.e. to perform better than their ‘larger’ counterparts over these two datasets. We hypothesise this phenomenon to be mostly due to the small size of these datasets, which encourages overfitting when using architectures with larger capacity. In Table 4 in the main paper, we report the result for the best performing architectures, along with the number of learnable parameters.