Breaking the Limits of Message Passing Graph Neural Networks

Introduction

In the past few years, finding the best inductive bias for relational data represented as graphs has gained a lot of interest in the machine learning community. Node-based message passing mechanisms relying on the graph structure have given rise to the first generation of Graph Neural Networks (GNNs) called Message Passing Neural Networks (MPNNs) (Gilmer et al., 2017). These algorithms spread each node features to the neighborhood nodes using trainable weights. These weights can be shared with respect to the distance between nodes (Chebnet GNN) (Defferrard et al., 2016), to the connected nodes features (GAT for graph attention network) (Veličković et al., 2018) and/or to edge features (Bresson & Laurent, 2018). When considering sparse graphs, the memory and computational complexity of such approaches are linear with respect to the number of nodes. As a consequence, these algorithms are feasible for large sparse graphs and thus have been applied with success on many downstream tasks (Dwivedi et al., 2020).

Despite these successes and these interesting computational properties, it has been shown that MPNNs are not powerful enough (Xu et al., 2019). Considering two non-isomorphic graphs that are not distinguishable by the first order Weisfeiler-Lehman test (known as the 1-WL test), existing maximum powerful MPNNs embed them to the same point. Thus, from a theoretical expressive power point of view, these algorithms are not more powerful than the 1-WL test. Beyond the graph isomorphism issue, it has also been shown that many other combinatorial problems on graph cannot be solved by MPNNs (Sato et al., 2019).

In (Maron et al., 2019b; Keriven & Peyré, 2019), it has been proven that in order to reach universal approximation, higher order relations are required. In this context, some powerful models that are equivalent to the 3-WL test were proposed. For instance, (Maron et al., 2019a) proposed the model PPGN (Provably Powerful Graph Network) that mimics the second order Folklore WL test (2-FWL), which is equivalent to the 3-WL test. In (Morris et al., 2019), they proposed to use message passing between 1, 2 and 3 order node tuples hierarchically, thus reaching the 3-WL expressive power. However, using such relations makes both memory usage and computational complexities grown exponentially. Thus, it is not feasible to have universal approximation models in practice.

In order to increase the theoretical expressive power of MPNNs by keeping the linear complexity mentioned above, some researchers proposed to partly randomize node features (Abboud et al., 2020; Sato et al., 2020) or to add a unique label (Murphy et al., 2019) in order to have the ability to distinguish two non-isomorphic graphs that are not distinguished by the 1-WL test. These solutions need massively training samples and involve slow convergence. (Bouritsas et al., 2020; Dasoulas et al., 2020) proposed to use a preprocessing step to extract some features that cannot be extracted by MPNNs. Thus, the expressive power of their GNN is improved. However, these handcrafted features need domain expertise and a feature selection process among an infinite number of possibilities.

All these studies target more theoretically powerful models, closer to universal approximation. However, this does not always induce a better generalization ability. Since most of the realistic problems are given with many node/edge features (which can be either continuous or discrete), there is almost no pair of graphs that are not distinguishable by the 1-WL test in practice. In addition, theoretically more powerful methods use non-local updates, breaking one of the most important inductive bias in Euclidean learning named locality principle (Battaglia et al., 2018). These may explain why theoretical powerful methods cannot outperform MPNNs on many downstream tasks, as reported in (Dwivedi et al., 2020). On the other hand, it is obvious that 1-WL equivalent GNNs are not expressive enough since they are not able to count some simple structural features such as cycles or triangles (Arvind et al., 2020; Chen et al., 2020; Bouritsas et al., 2020; Vignac et al., 2020), which are informative for some social or chemical graphs. Finally, another important aspect mentioned by a recent paper (Balcilar et al., 2021) concerns the spectral ability of GNN models. It is shown that a vast majority of the MPNNs actually work as low-pass filters, thus reducing their expressive power.

In this paper, we propose to design graph convolution in the spectral domain with custom non-linear functions of eigenvalues and by masking the convolution support with desired length of receptive field. In this way, we have (i) a spatially local updates process, (ii) linear memory and computational complexities (except the eigendecomposition in preprocessing step), (iii) enough spectral ability and (iv) a model that is theoretically more powerful than the 1-WL test, and experimentally as powerful as PPGN. Experiments show that the proposed model can distinguish pairs of graphs that cannot be distinguished by 1-WL equivalent MPNNs. It is also able to count some substructures that 1-WL equivalent MPNNs cannot. Its spectral ability enables to produce various kind of spectral components in the output, while the vast majority of the GNNs including higher order WL equivalent models do not. Finally, thanks to the sparse matrix multiplication, it has linear time complexity except the eigendecomposition in preprocessing step.

The paper is structured as follows. In Section 2, we set the notations and the general framework used in the following. Section 3 is dedicated to the characterization of WL test, which is the backbone of our theoretical analysis. It is followed by our findings in Section 4 on analysing the expressive power of MPNNs and our solutions to improve expressive power of MPNNs in Section 5. The experimental results and conclusion are the last two section of this paper.

Generalization of Spectral and Spatial MPNN

GNN models rely on a set of layers where each layer takes the node representation of the previous layer $H^{(l-1)}$ as input and produces a new representation $H^{(l)}$, with $H^{(0)}=X$. According to the domain which is considered to design the layer computations, GNNs are generally classified as either spectral or spatial (Wu et al., 2019; Chami et al., 2020). Spectral GNNs rely on the spectral graph theory (Chung, 1997). In this framework, signals on graphs are filtered using the eigendecomposition of the graph Laplacian (Shuman et al., 2013). By transposing the convolution theorem to graphs, the spectral filtering in the frequency domain can be defined by $x_{flt}=U{diag}(\Omega(\boldsymbol{\lambda}))U^{\top}x$, where $\Omega(.)$ is the desired filter function which needs to be learnt by back-propagation. On the other hand, spatial GNNs, such as GCN (graph convolutional network) (Kipf & Welling, 2017) and GraphSage (Hamilton et al., 2017), consider two operators, one that aggregates the connected nodes messages and one that updates the concerned node representation.

In a recent paper (Balcilar et al., 2021), it was explicitly shown that both spatial and spectral GNNs are MPNN, taking the general form

One can see that as long as $C^{(s)}$ matrices are sparse (number of edges is defined by some constant multiplied by the number of nodes), MPNN in Eq.1 has linear memory and computational complexities with respect to the number of nodes. Because, the valid entries in $C^{(s)}$ that we need to keep is linear with respect to the number of nodes and thank to the sparse matrix multiplication $C^{(s)}H^{(l)}$ takes linear time with respect to the number of edges thus nodes as well.

Characterization of Weisfeiler-Lehman

The universality of a GNN is based on its ability to embed two non-isomorphic graphs to distinct points in the target feature space. A model that can distinguish all pairs of non-isomorphic graphs is a universal approximator. Since the graph isomorphism problem is NP-intermediate (Takapoui & Boyd, 2016), the Weisfeiler-Lehman Test (abbreviated WL-test), which gives sufficient but not enough evidence of graph isomorphism, is frequently used for characterizing GNN expressive power. The classical vertex coloring WL test can be extended by taking into account higher order of node tuple within the iterative process. These extensions are denoted as $k$-WL test, where $k$ is equals to the order of the tuple. These tests are described in Appendix A.

It is shown in (Arvind et al., 2020) that for $k\geq 2$, $(k+1)$-WL $>$ $(k)$-WL, i.e., higher order of tuple leads to a better ability to distinguish two non-isomorphic graphs. For $k=1$, this statement is not true, and 2-WL is not more powerful than 1-WL (Maron et al., 2019a). To clarify this point, the Folkore WL (FWL) test has been defined such that 1-WL=1-FWL, but for $k\geq 2$, we have $(k+1)$-WL $\approx$ $(k)$-FWL (Maron et al., 2019a).

In literature, some confusions occur among the two versions. Some papers use WL test order (Morris et al., 2019; Maron et al., 2019a), while others use FWL order under the name of WL such as in (Abboud et al., 2020; Arvind et al., 2020; Takapoui & Boyd, 2016). In this paper, we explicitly mention both WL and FWL equivalent.

In order to better understand the capability of WL tests, some papers attempt to characterize these tests using a first order logic (Immerman & Lander, 1990; Barceló et al., 2019). Consider two unlabeled and undirected graphs represented by their adjacency matrices $A_{G}$ and $A_{H}$. These two graphs are said $k$-WL (or $k$-FWL) equivalent, and denoted $A_{G}\equiv_{k-WL}A_{H}$, if they are indistinguishable by a $k$-WL (or $k$-FWL) test.

Recently (Brijder et al., 2019; Geerts, 2020) proposed a new Matrix Language called MATLANG. This language includes different operations on matrices and makes some explicit connections between specific dictionaries of operations and the 1-WL and 3-WL tests. Expressive power varies with the operations included in each dictionnary.

$ML(\mathcal{L})$ is a matrix language with an allowed operation set $\mathcal{L}=\{op_{1},\dots op_{n}\}$, where $op_{i}\in\{.,+,^{\top},diag,tr,{\bf 1},\odot,\times,f\}$. The possible operations are matrices multiplication and addition, matrix transpose, vector diagonalization, matrix trace computation, column vector full of 1, element-wise matrix multiplication, matrix/scalar multiplication and element-wise custom function operating on scalars or vectors.

As an example, $e(X)={\bf 1}^{\top}X^{2}{\bf 1}$ is a sentence of $ML(\mathcal{L})$ with $\mathcal{L}=\{.,^{\top},{\bf 1}\}$, computing the sum of all elements of square matrix $X$. In the following, we are interested in languages $\mathcal{L}_{1},\mathcal{L}_{2}$ and $\mathcal{L}_{3}$ that have been used for characterizing the WL-test in (Geerts, 2020). These results are given next.

Two adjacency matrices are indistinguishable by the 1-WL test if and only if $e(A_{G})=e(A_{H})$ for all $e\in\mathcal{L}_{1}$ with $\mathcal{L}_{1}=\{.,^{\top},{\bf 1},diag\}$. Hence, all possible sentences in $\mathcal{L}_{1}$ are the same for 1-WL equivalent adjacency matrices. Thus, $A_{G}\equiv_{1-WL}A_{H}\leftrightarrow A_{G}\equiv_{ML(\mathcal{L}_{1})}A_{H}$. (see Theorem 7.1 in (Geerts, 2020))

$ML(\mathcal{L}_{2})$ with $\mathcal{L}_{2}=\{.,^{\top},{\bf 1},diag,tr\}$ is strictly more powerful than $\mathcal{L}_{1}$, i.e., than the 1-WL test, but less powerful than the 3-WL test. (see Theorem 7.2 and Example 7.3 in (Geerts, 2020))

Two adjacency matrices are indistinguishable by the 3-WL test if and only if they are indistinguishable by any sentence in $ML(\mathcal{L}_{3})$ with $\mathcal{L}_{3}=\{.,^{\top},{\bf 1},diag,tr,\odot\}$. Thus, $A_{G}\equiv_{3-WL}A_{H}\leftrightarrow A_{G}\equiv_{ML(\mathcal{L}_{3})}A_{H}$. (see Theorem 9.2 in (Geerts, 2020))

Enriching the operation set to $\mathcal{L}^{+}=\mathcal{L}\cup\{+,\times,f\}$ where $\mathcal{L}\in(\mathcal{L}_{1},\mathcal{L}_{2},\mathcal{L}_{3}$) does not improve the expressive power of the language. Thus, $A_{G}\equiv_{ML(\mathcal{L})}A_{H}\leftrightarrow A_{G}\equiv_{ML(\mathcal{L}^{+})}A_{H}$. (see Proposition 7.5 in (Geerts, 2020))

How Powerful are MPNNs?

This section presents some results about the theoretical expressive power of state-of-the-art MPNNs. Those results are derived using the MATLANG language (Geerts, 2020) and more precisely the remarks of the preceding section. Proofs of the theorems are given in Appendix B.

MPNNs such as GCN, GAT, GraphSage, GIN (defined in Appendix H) cannot go further than operations in $\mathcal{L}_{1}^{+}$. Thus, they are not more powerful than the 1-WL test.

This result has already been given in (Xu et al., 2019), which proposed GIN-$\epsilon$ (GIN for Graph Isomorphism Network) and showed that it is the unique MPNN which is provably exact the same powerful with the 1-WL test, while the rest of MPNNs are known to be less powerful than 1-WL test.

Chebnet is also known to be not more powerful than the 1-WL test. However, the next theorem states that it is true if the maximum eigenvalues are the same for both graphs. For a pair of graphs whose maximum eigenvalues are not equal, Chebnet is strictly more powerful than the 1-WL test.

Chebnet is more powerful than the 1-WL test if the Laplacian maximum eigenvalues of the non-regular graphs to be compared are not the same. Otherwise Chebnet is not more powerful than 1-WL.

Figure 1 shows two graphs that are 1-WL equivalent and are generally used to show how MPNNs fail. However, their normalized Laplacian’s maximum eigenvalues are not the same. Thus, Chebnet can project these two graphs to different points in feature space. Details can be found in Appendix C.

As stated in the introduction, comparison with the WL-test is not the only way to characterize the expressive power of GNNs. Powerful GNNs are also expected to be able to count relevant substructures in a given graph for specific problems. The following theorems describe the matrix language required to be able to count the graphlets illustrated in Figure 2, which are called 3-star, triangle, tailed triangle and 4-cycle.

3-star graphlets can be counted by sentences in $\mathcal{L}_{1}^{+}$.

Triangle and 4-cycle graphlets can be counted by sentences in $\mathcal{L}_{2}^{+}$.

Tailed triangle graphlets can be counted by sentences in $\mathcal{L}_{3}^{+}$.

These theorems show that 1-WL equivalent MPNNs can only count 3-star patterns, while 3-WL equivalent MPNNs can count all graphlets shown in Figure 2.

(Dehmamy et al., 2019) has shown that a MPNN is not able to learn node degrees if the MPNN has not an appropriate convolution support (e.g. $A$). Therefore, to achieve a fair comparison, we assume that node degrees are included as a node feature. Note however, that the number of 3-star graphlets centered on a node can be directly derived from its degrees (see Appendix B.3). Therefore, any graph agnostic MLP can count the number of 3-star graphlets given the node degree.

MPNN Beyond 1-WL

In this section, we present two new MPNN models. The first one, called GNNML1 is shown to be as powerful as the 1-WL test. The second one, called GNNML3 exploits the theoretical results of (Geerts, 2020) to break the limits of 1-WL and reach 3-WL equivalence experimentally. GNNML1 relies on the node update schema given by :

where $W^{(l,s)}$ are trainable parameters. Using this model, the new representation of a node consists of a sum of three terms : (i) a linear transformation of the previous layer representation of the node, (ii) a linear transformation of the sum of the previous layer representations of its connected nodes and (iii) the element-wise multiplication of two different linear transformations of the previous layer representation of the node.

The expressive power of GNNML1 is defined by the following theorem. Its proof is given in Appendix B:

GNNML1 can produce every possible sentences in $ML(\mathcal{L}_{1})$ for undirected graph adjacency $A$ with monochromatic edges and nodes. Thus, GNNML1 is exactly as powerful as the 1-WL test.

Hence, this model has the same ability as the 1-WL test to distinguish two non-isomorphic graphs, i.e., the same as GIN. This is explained by the third term in the sum of Eq.(2) since it can produce feature-wise multiplication on each layer. Since node representation is richer, we also assume that it would be more powerful for counting substructures. This assumption is validated by experiments in Section 6.

To reach more powerful models than 1-WL, theoretical results (see Remarks 1, 2 and 3 in Section 3) show that a model that can produce different outputs than $\mathcal{L}_{1}^{+}$ language is needed. More precisely, according to Remarks 2 and 3, trace ($tr$) and element-wise multiplication ($\odot$) operations are required to go further than 1-WL.

In order to illustrate the impact of the trace operation, one can use 1-WL equivalent Decalin and Bicyclopentyl graphs in Figure 1. It is easy to show that $tr(A_{G}^{5})=0$ but $tr(A_{H}^{5})=20$, $tr(A^{5})$ giving the number of 5-length closed walks. Thus, if a model can apply a trace operator over some power of adjacency, it can easily distinguish these two graphs. Computational details concerning this example are given in Appendix C.

Despite this interesting property of the trace operator, it is not sufficient to distinguish cospectral graphs, since cospectral graphs (see Figure 3) have the same number of closed walks of any length (see Proposition 5.1 in (Geerts, 2020)). In such cases, element-wise multiplication is useful. As an example, the sentence $e(A)={\bf 1}^{\top}f((A\odot A^{2})^{2}{\bf 1})$ where $f(x)=x\odot x$ for any vector $x$, gives $e(A_{G})=6032$ and $e(A_{H})=5872$ for the graphs of Figure 3. Thus, element-wise multiplication helps distinguishing these two graphs. The calculation details can be found in Appendix D.

As shown by these examples, a model enriched by element-wise multiplication and trace operator can go further than the 1-WL test. However, these operations need to keep the power of the adjacency matrix explicitly and to multiply these dense matrices to each other by matrix or element-wise multiplication. Such a strategy is actually used by higher order GNNs such as (Maron et al., 2019a; Morris et al., 2019), which are provably more powerful than existing MPNNs.

However, MPNNs cannot calculate the power of a given adjacency explicitly. Indeed, a MPNN layer multiplies the previous representation of the nodes by sparse adjacency matrix or more generally sparse convolution supports $C$ in Eq.(1). More precisely, if the given node features are $H^{(0)}={\bf 1}$, a MPNN can calculate $C^{3}{\bf 1}$ by 3 layered MPNN computing $C(C(C{\bf 1}))$ but not by $(C^{3}){\bf 1}$. Since a MPNN does not keep $C^{3}$ explicitly, it cannot take its trace or multiply element-wise to another power of support. This is a major disadvantage of MPNNs, but it explains why MPNNs need just linear time and memory complexity, making them useful in practice.

A solution to the problem mentioned above is to design graph convolution supports by the element-wise multiplication of the $s$-power of the adjacency matrix and a given receptive field, i.e., by $C^{(s)}=M\odot A^{s}$ where $M$ masks the components of the powered matrix and keeps the convolution support sparse. $M=A+I$ is an example of mask that gives a maximum 1-length receptive field. This model cannot calculate all possible element-wise multiplications between all possible matrices, but it can produce any sentence in a form of $(M\odot A^{s})^{l}$ where $l\in[0,l_{max}]$ is the layer number and $s\in[0,s_{max}]$ is the pre-computed power of convolution supports. In this proposition, the receptive field mask and the number of power of adjacency should be computed in a pre-processing step. However, we cannot initially know which power of adjacency matrix is necessary for a given problem. One solution is to tune it as an hyperparameter of the model. Another problem of this approach is that using powers of adjacency makes the convolution supports filled with high values that have to be normalized.

To overcome these problems, we propose through our GNNML3 model to design convolution supports in the spectral domain as functions of eigenvalues of the normalized Laplacian matrix or of the adjacency matrix. The following theorem, with proof given in Appendix B, shows that such supports can be written as power series of the graph Laplacian or the adjacency matrix.

where $\Phi_{s}(\lambda)=exp(-b(\lambda-f_{s})^{2})$, $f_{s}\in[\lambda_{min},\lambda_{max}]$ is a scalar design parameter of each convolution support and $b>0$ is a general scalar design parameter, can be expressed as a linear combination of all powers of graph Laplacian (or adjacency) as follows, with $\alpha_{s,i}=\frac{\Phi_{s}^{(i)}(0)}{i!}$:

Since design parameters $f_{s}$ of each matrix are different, each $C^{\prime(s)}$ in Eq.(4) consists of different linear combinations of power series of the graph Laplacian (or adjacency). Thus, necessary powers of the graph Laplacian (or adjacency) and its diagonal part (for trace operation) can be learned and their element-wise multiplication can be produced by:

where we concatenate MPNN representation under learned convolution with element-wise product of node representations as in GNNML1.

There is an infinite number of selections of $\Phi_{s}(\lambda)$ that make the convolution support written by power series of graph Laplacian (or adjacency). However, we can design each convolution support to be sensitive on each band of spectrum ($f_{s}$) by given bandwidth ($b$). Therefore, our model will be able to learn properties depending on the spectrum of graph signal.

The limit of the proposed method is similar to the limit of 3-WL (or 2-FWL) test. For instance, it fails to distinguish strongly regular graphs, that can be defined by 3 parameters: the degree of the nodes, the number of common neighbours of adjacent node pairs, and the number of common neighbours of non-adjacent node pairs. Such graphs are provably known to be 3-WL equivalent (Arvind et al., 2020). In Appendix E, a strongly regular graphs pair and the result of a sample sentence in $\mathcal{L}_{3}$ are presented.

Experimental Results

This section presents the experimental results obtained by the proposed models GNNML1 and GNNML3. All codes and datasets are available online https://github.com/balcilar/gnn-matlang. We use GCN, GAT, GIN and Chebnet as 1-WL MPNN baselines and PPGN as 3-WL baseline (see Appendix H). Experiments aim to answer four questions:

Q1: How many pairs of non-isomorphic simple graphs that are either 1-WL or 3-WL equivalent are not distinguished by the models? Q2: Can the models generalize the counting of some substructures in a given graph? Q3: Can the models learn low-pass, high-pass and band-pass filtering effects and generalize the classification problem according to the frequency of the signal? Q4: Can the models generalize downstream graph classification and regression tasks?

In order to perform experimental expressive power tests, we use graph8c and sr25 datasetshttp://users.cecs.anu.edu.au/$\sim$bdm/data/graphs.html. Graph8c is composed of all the $11\,117$ possible connected non-isomorphic simple graphs with 8 nodes. We compare all possible pairs of graphs of this dataset, leading to more than 61M comparisons. According to our test, we found that 312 pairs out of 61M are 1-WL equivalent and none of the pairs are 3-WL equivalent. The sr25 dataset contains strongly regular graphs where each graph has 25 nodes, each node’s degree is 12, connected nodes share 5 common neighbours and non-connected nodes share 6 common neighbors. Sr25 consists of 15 graphs, leading to 105 different pairs for comparison.

Moreover, we use the EXP dataset (Abboud et al., 2020), having 600 pairs of 1-WL equivalent graphs. This dataset also includes a binary classification task. Depending on graph features, each graph of a pair of 1-WL equivalent graphs is assigned to two different classes. We split the dataset into 400, 100, and 100 pairs for train, validation and test sets respectively. The test set is used to measure the generalization ability: a model that fails to distinguish 1-WL equivalent graphs inevitably fails to learn this task.

We use 3-layer graph convolution followed by sum readout layer, and then a linear layer to convert the readout layer representation into a 10-length feature vector. We keep the parameter budget around 30K for all methods. For graph8c, sr25 and EXP tasks, there is no learning. Model weights are randomly initialized and 10-length graph representations are compared by the Manhattan distance. If the distance is less than $10^{-3}$ in all 100 independent runs, we assume the pairs are similar. For EXP-classification task, we train the model and pick the best one according to validation set performance and report its performance on test set.

Table 1 presents the obtained results. One can see that 99.5% of the graphs in graph8c dataset can be distinguished even by graph agnostic method MLP (293K out of 61M is not separable by MLP). This can be explained by the fact that the node degrees has been added as node features. Hence, all methods initially know the result of first iteration of 1-WL test. Thus, MLP (and also first iteration of 1-WL test) can distinguish pairs of graphs when multiset of node degrees are not same. GNNML1 and GIN’s result is very closed to the theoretical limit of 1-WL test which is 312 pairs for graph8c dataset. The difference can be explained by threshold value to make decision if the two representations are equal and/or the number of layers in the model. It is possible that 1-WL test may need more than 3 iteration to distinguish some pairs. Due to having less expressive power of GCN and GAT compare to the 1-WL test, their performances are worse than 1-WL test. Since graph8c dataset has 1-WL equivalent non-regular graph pairs that have different maximum eigenvalue, Chebnet could detect these pairs and reaches better performance than theoretical limit of 1-WL test as stated by Theorem 2.

On EXP dataset, composed of 1-WL equivalent graph pairs, MPNNs cannot distinguish any pair of graphs, except Chebnet which is able to distinguish all the pairs with different maximum eigenvalues. In EXP there is no regular graphs and only 71 graph pairs have similar maximum eigenvalues. Chebnet fails on these pairs but distinguishes the others, as stated by Theorem 2. One can note that using a fixed value for maximum eigenvalue (e.g. $\lambda_{max}=2$ as it is usually done in practice) reduces Chebnet performance to those of MPNNs.

Similarly to results on EXP, 1-WL equivalent MPNNs except Chebnet fail to predict of EXP classification task and do not perform better than random prediction. On the contrary, PPGN and GNNML3 have perfect results on graph8c, EXP and EXP-classify tasks thanks to their 3-WL equivalence. However, since strongly regular graphs are 3-WL equivalent, no model less or as powerful as 3-WL test can distinguish the pairs in sr25 dataset. To obtain a better result on this dataset, we need to go further than 3-WL (see Appendix E). These experiments reply to Q1.

To bring an answer to Q2, we propose to count 3-star, triangle, tailed-triangle and 4-cycle substructures (Fig. 2). In addition to these 4 graphlets, we also create another task (noted as CUSTOM in Table 2) that aims to approximate a custom sentence $e_{c}\in\mathcal{L}_{1}^{+}$ , $e_{c}(A)={\bf 1}^{\top}A\,diag(exp(-A^{2}{\bf 1}))A{\bf 1}$ with $A$ the graph adjacency matrix. Since $e_{c}\in ML(\mathcal{L}_{1}^{+})$, it may be learnable by 1-WL equivalent MPNNs. We used the RandomGraph dataset (Chen et al., 2020) with same partitioning: 1500, 1000 and 2500 graphs for train, validation and test respectively. To create the ground truth of number of graphlets, we count them according to theorem proofs in Appendix B.3, B.4, B.5 and normalized the number to a unitary standard deviations, to keep the errors in the same scale as in Table 2. We use 4 convolution layers, a graph readout layer computing a sum and followed by 2 fully connected layers. All methods parameter budget is around 30K. We keep the maximum number of iterations to 200 and we stop the algorithm if the error goes below ${10}^{-4}$.

The results in Table 2 are consistent with Theorems 3, 4, 5. 3-WL models are able to count graphlets and approximate our custom function (result $<$ ${10}^{-3}$), while 1-WL equivalent models can only count the 3-stars graphlet, as stated in Theorem 3. Custom function approximation results also show that GNNML1 and Chebnet provide better approximation of the target other MPNNs, which is again consistent with our analysis.

Question Q3 concerns the spectral expressive power of models. Such an analysis is important when input-output relations depend on the spectral properties of the graph signal such as in image/signal processing applications. As shown in (Balcilar et al., 2021), the vast majority of existing MPNNs operate as low-pass filters which limits their capacity. To lead this analysis, we use the datasets presented in (Balcilar et al., 2021). First, we evaluate if the models can learn low-pass, high-pass and band-pass filtering effects, through a node regression problem. Model performances are thus reported $R^{2}$ using mean square error (MSE) loss. The original data consists in a 2-d grid graph of size 100x100. Since the PPGN’s memory and computational complexity is prohibitive with a reasonable computer, we select 3 different 30x30 regions of the original 2-d grid graph as training, validation and test sets. A second dataset consists of 5K planar graphs, split into 3K, 1K and 1K sets for train, validation and test. They are used to evaluate if the models can classify graphs into binary classes where the ground truth labels were determined according to the frequency of the signal on the graph. Since the problem is binary graph classification we use binary cross entropy loss.

The results of spectral expressive power analysis are presented in Table 3. Node regression results show that 1-WL equivalent existing MPNNs can mostly learn low-pass effects. By applying different weights to self node and neighbourhood, GNNML1 can learn high pass effect relatively well. PPGN also learns high-pass effect better than 1-WL equivalent methods. Band-pass can be generalized by Chebnet and GNNML3 thanks to the convolutions designed in spectral domain. The reason why the band-pass regression results are worse than the low and high-pass results is that the ground truth band-pass effect is created by very stiff frequency function and Chebnet also GNNML3 need more convolution supports to learn it. Because of non-local process in PPGN, it cannot learn the band-pass effect and provide no better result than 1-WL MPNNs in graph classification problem. Thus, Chebnet and GNNML3 give the best results on all spectral ability test, thanks to their spectral convolutions process.

For answering the last question Q4, we apply the different models on some common benchmark tasks and datasets. Table 4 and Table 5 present the performance of both baseline models and the proposed ones on these benchmark datasets. The results on Zinc12K and MNIST-75 datasets are very interesting because of the nature of these two problems. The solution of the Zinc12K dataset mostly depends on structural features of the graph. For instance, a recent study reaches 0.14 MAE by using handcrafted features, which cannot be extracted by a 3-WL equivalent model (Bouritsas et al., 2020). Obtained results confirm that models that are able to count substructures, such as PPGN and GNNML3, perform better than others with a large margin. On the other hand, since MNIST-75 dataset is based on image analysis, it needs a model with a higher spectral ability. Therefore, Chebnet and GNNML3 perform significantly better than other models on this task. Our proposal GNNML3 gives comparable results on other TU datasets in (Morris et al., 2020) such as MUTAG, ENZYMES, PROTEINS and PTC presented in Appendix F.

Conclusion

Despite a computational and memory efficiency, MPNN is known to have an expressive power limited to 1-WL test. MPNN is then unable to distinguish 1-WL equivalent graphs and cannot count some substructures of the graph. In this paper, we have presented new models, by translating the insights of MATLANG to the GNN world. This solution gives access to a new MPNN that is theoretically more powerful than the 1-WL test, and experimentally as powerful as 3-WL existing models for distinguishing non-isomorphic graphs and for counting substructures without feature engineering nor node permutations in the training phase. The proposed MPNN is also powerful in terms of spectral expressive ability, going beyond low-pass filtering, which is another expressive perspective of GNNs. Experimental results confirm the theorems stated in the paper. The proposed method has a big advantage over all studied MPNN on graph isomorphism and substructure counting tasks. With respect to the 3-WL equivalent baseline PPGN, the biggest advantage of our proposal is its complexity. Proposed GNNML3 needs linear memory and time complexity with respect to the number of nodes, while PPGN needs quadratic memory and cubic time complexity, making the model infeasible for large graphs. The second advantage over PPGN is that since it is created in the spectral domain, its convolution process takes care of signal frequencies, making it more efficient in terms of output signal frequency profile.

Acknowledgments

This work was partially supported by the Normandy Region (grant RiderNet), the French Agence National de Recherche (grant APi, ANR-18-CE23-0014) and the PAUSE Program of Collège de France.

References

Appendix A Weisfeiler-Lehman Test

The universality of a GNN is based on its ability to embed two non-isomorphic graphs to distinct points in the target feature space. A model which can distinguish all pairs of non-isomorphic graphs is a universal approximator. Since it is not known if the graph isomorphism problem can be solved in polynomial time or not, this problem is neither NP-complete nor P, but NP-intermediate (Takapoui & Boyd, 2016). One of the oldest but prominent polynomial approach is the Weisfeiler-Lehman Test (abbreviated WL-test) which gives sufficient but not enough evidence. WL test can be extended by taking into account higher order of node tuple within the iterative process. These extensions are denoted as $k$-WL test, where $k$ is equal to the order of the tuple. It is important to mention that an higher order of tuple leads to a better ability to distinguish two non-isomorphic graphs (with the exception for $k=2$) (Arvind et al., 2020).

The 1-WL test, known as vertex coloring, starts with the given initial color of nodes if available. Otherwise all nodes are colored with the same color ($H_{v}^{(0)}=1$). Then, colors are updated by the following iteration:

where $H_{v}^{(t)}$ is the color of vertex $v$ at iteration $t$, $\mathcal{N}(v)$ is the set of neighbours of vertex $v$, $|$ represents the concatenation operator and $\{.\}$ is the order invariant multisetIt is generally implemented by stacking all colors in the set and sorting them alphabetically. In order to avoid the new color of vertex become bigger after each iteration due to the concatenation operation and to keep the color description simple, the recoloring $\sigma(\cdot)$ function is applied after each iteration. It assigns a new simple color identifier to the any newly created color. The test is performed in parallel for two graphs. The iterative process is stopped when the color histograms are kept unchanged between two consecutive iterations. The color histograms associated to the compared graphs are examined. If in any iteration the histograms are different, we can conclude that the graphs are not isomorphic. However, the opposite conclusion can not be drawn if color histograms are equal as two same histograms may be computed even for non-isomorphic graphs.

Then, the iteration process is applied through the following schema where $[n]$ is the set of node identifiers.

Although for $k\geq 2$, $(k+1)$-WL is more powerful than $(k)$-WL, it is not true for $k=1$, thus 2-WL (Eq.(9)) is no more powerful than 1-WL (Eq.(7)) (Maron et al., 2019a). To clarify this point, Folkore WL (FWL) test is defined such that 1-WL=1-FWL, but for $k\geq 2$, we have $(k+1)$-WL $\approx$ $(k)$-FWL (Maron et al., 2019a). The iteration process of 2-FWL is given by the following equation;

In the literature, there are different interpretations of the order of the WL test. Some papers use WL test order to denote the iteration given by Eq.(7) and Eq.(9) (Morris et al., 2019; Maron et al., 2019a) but some others such as (Abboud et al., 2020; Arvind et al., 2020; Takapoui & Boyd, 2016) use FWL order under the name of WL. In this paper, we explicitly mention both WL and FWL equivalent such as 3-WL (or 2-FWL) to alleviate ambiguities.

Appendix B Proofs of Theorems

All these methods can be written in Eq.(1) by different convolution matrices $C$. The main idea of the proof is that as long as convolution matrices $C$ can be explained by operations from the enriched set $\mathcal{L}_{1}^{+}$ (Remark 4), Eq.(1) also can be explained by operations from $\mathcal{L}_{1}^{+}$ as well. Thus these methods cannot produce any sentence out of $\mathcal{L}_{1}^{+}$. As a consequence, their expressive power is not more than 1-WL test. To provide a proof, the mentioned methods’ convolution matrices have to be expressed using operations from $\mathcal{L}_{1}^{+}$.

GCN uses $C=(D+I)^{-0.5}(A+I)(D+I)^{-0.5}$ where $D$ is the diagonal degree matrix (Kipf & Welling, 2017) in Eq.(1). $(D+I)^{-0.5}$ can be expressed as $(D+I)^{-0.5}=diag(f(A{\bf 1}+{\bf 1}))$, where $f(x)=x^{-0.5}$ is element-wise operation on vector $x$. $A+I$ can also be written $A+diag({\bf 1})$. When we merge these equations, we get $C=diag(f(A{\bf 1}+{\bf 1}))(A+diag({\bf 1}))diag(f(A{\bf 1}+{\bf 1}))$. The convolution support $C$ is then written using operations from $\mathcal{L}_{1}^{+}$.

In the literature, GraphSage method was proposed to sample neighborhood and aggregate the neighborhood contribution by the mean operator or LSTM in (Hamilton et al., 2017). Since we restrict the method using full sampling and mean aggregator, we can define GraphSage by the general framework given by Eq.(1) with two convolution supports which are the identity matrix $C^{(1)}=I$ and the row normalized adjacency matrix $C^{(2)}=D^{-1}A$. These convolution supports can also be expressed by operations from $\mathcal{L}_{1}^{+}$, by observing that $C^{(1)}=diag({\bf 1})$ and $C^{(2)}=diag(f(A{\bf 1}))A$, where $f(x)=x^{-1}$ elementwise operation on vector $x$.

GIN (Xu et al., 2019) uses a convolution support $C=A+I\epsilon$ in Eq.(1) which is followed by a custom number of MLP layers. Each of these layers correspond to a convolution support that can by expressed as $C_{mlp}=I$ in Eq.(1). Finally, these convolution supports can be written thanks to operations from $\mathcal{L}_{1}^{+}$. $C=A+\epsilon\times diag({\bf 1})$ and $C_{mlp}=diag({\bf 1})$.

B.2 Theorem.2

Chebnet (Defferrard et al., 2016) uses desired number $k$ of convolution supports in Eq.(1). As long as these convolutions can be written by operations in $\mathcal{L}_{1}^{+}$, we can conclude that Chebnet is no more powerful than 1-WL test. But if at least one convolution cannot be explained in $\mathcal{L}_{1}^{+}$, we can say it is more powerful than 1-WL test.

Chebnet’s convolution supports are $C^{(1)}=I,~{}~{}C^{(2)}=2L/\lambda_{\max}-I,~{}~{}C^{(k)}=2C^{(2)}C^{(k-1)}-C^{(k-2)}$. The first support can always be written thanks to an operation from $\mathcal{L}_{1}$ since $C^{(1)}=diag({\bf 1})$. Both normalized and combinatorial graph Laplacian can also be written as $L=diag(A{\bf 1})-A$ or $L=diag({\bf 1})-diag(f(A{\bf 1}))Adiag(f(A{\bf 1}))$ where $f(x)=x^{-1/2}$ elementwise operation on vector $x$. If $\lambda_{max}$ for both graphs are the same, we can use a constant $\alpha=2/\lambda_{max}$. The second convolution support can then be written as $C^{(2)}=\alpha\times L-diag({\bf 1})$. It is then expressed by means of operations from $\mathcal{L}_{1}^{+}$. Other convolution supports $C^{(k)}=2C^{(2)}C^{(k-1)}-C^{(k-2)}$ are created by matrix multiplication and subtraction of previous supports which can all be expressed by mean of operations from $\mathcal{L}_{1}^{+}$. Thus, if the maximum eigenvalues of tested graphs Laplacians are the same, Chebnet is not more powerful than 1-WL.

However, if the maximum eigenvalues are not the same, $C^{(2)}$ cannot be expressed with the help of the constant value $\alpha$. It means that different coefficients should be used for each graph. For two tested graphs $G$ and $H$, we can write second kernel of Chebnet as $C^{(2)}_{G}=\alpha_{G}\times L_{G}-diag({\bf 1})$ and $C^{(2)}_{H}=\alpha_{H}\times L_{H}-diag({\bf 1})$. If these two graphs are 1-WL equivalent, any sentence build on $\mathcal{L}_{1}^{+}$ applied on these graph is equivalent as well. For instance, we can use the sentences of $e(X)={\bf 1}^{\top}X{\bf 1}$ with operation in $\mathcal{L}_{1}^{+}$. The output of the sentence should be same such $e(L_{G})=e(L_{H})$ yields ${\bf 1}^{\top}L_{G}{\bf 1}={\bf 1}^{\top}L_{H}{\bf 1}$. If we assume that Chebnet cannot separate these two graphs, we can calculate one layer ChebNet’s output by second support with the same sentence and they should be the same such $e(C^{(2)}_{G})=e(C^{(2)}_{H})$ yields $\alpha_{G}{\bf 1}^{\top}L_{G}{\bf 1}=\alpha_{H}{\bf 1}^{\top}L_{H}{\bf 1}$. Last equation has contradiction to the previous one as long as the maximum eigenvalues are not same (i.e $\alpha_{G}\neq\alpha_{H}$) and graphs are not regular (i.e ${\bf 1}^{\top}L_{G}{\bf 1}>0$ and ${\bf 1}^{\top}L_{H}{\bf 1}>0$ for normalized laplacian). This contradiction says that assumption is wrong, so one layer Chebnet’s second support can distinguish 1-WL equivalent graphs whose maximum eigenvalues are not same and graphs are not regular with the same degree. ∎

Since the graph laplacians are positive semi-definite, it always yields ${\bf 1}^{\top}L_{G}{\bf 1}\geq 0$ and ${\bf 1}^{\top}L_{H}{\bf 1}\geq 0$ and they are zero as long as the graphs are regular with the same degree. Thus, if we add smallest positive scalar value on the diagonal of the laplacian such $L\leftarrow L+\epsilon I$, we get rid of the necessity that graphs must be non-regular. So Chebnet become more powerful and will be able to distinguish all 1-WL equivalent regular graphs whose maximum eigenvalues are different. Considering the graph8c task, we have seen that classic ChebNet could not distinguish 44 pairs where there are 312 1-WL equivalent pairs. If we use $L\leftarrow L+0.01I$, the number of undistinguished pairs of graph decreased from 44 to 19, where 19 undistinguished pairs are all 1-WL equivalent and have exact the same maximum eigenvalues. On the other hand, original Chebnet was not able to distinguish 44-19=25 graphs pairs whose maximum eigenvalues are different but all of them are regular thus ${\bf 1}^{\top}L{\bf 1}=0$.

B.3 Theorem.3

The number of 3-star patterns can be determined by $\sum_{v}\binom{d(v)}{3}$ where $d(v)$ is the degree of vertex $v$ for undirected simple graphs (Pinar et al., 2017). Using $f(x)=\frac{x!}{(x-3)!3!}$ as a function that operates on each element of a given vector $x$, we can calculate the number of 3-star patterns in a given adjacency matrix $A$ by ${\bf 1}^{\top}f(A{\bf 1})$ using operations in $\mathcal{L}_{1}^{+}$. According to the universal approximation theory of multi layer perceptron (Hornik et al., 1989), if we have enough layers, we can implement $f(.)$ as an MLP in our model. ∎

B.4 Theorem.4

The number of triangles can be determined by using trace operator as $1/6\times tr(A^{3})$ (Harary & Manvel, 1971) which can be written by means of operations from $\mathcal{L}_{2}^{+}$.

Number of 4-cycles is determined by $1/8\times(tr(A^{4})+tr(A^{2})-2{\bf 1}^{\top}A^{2}{\bf 1})$ (Harary & Manvel, 1971) which can be written by means of operations from $\mathcal{L}_{2}^{+}$. ∎

B.5 Theorem.5

If $t(v)$ denotes the number triangles including vertex $v$ and $d(v)$ denotes the degree of vertex $v$, the number of tailed triangles can be found by $\sum_{v}t(v).(d(v)-2)$ for simple undirected graphs (Pinar et al., 2017). Every node in a triangle has two closed walks of length 3. Thus, $t(v)=\frac{(A^{3})_{v,v}}{2}$. It yields the number of tailed triangles can be found by $\frac{1}{2}\times{\bf 1}^{\top}(A^{3}\odot diag(A{\bf 1}-2)){\bf 1}$. The computation of $t(v)$ which involves the element-wise multiplication can be written with operations from $\mathcal{L}_{3}^{+}$. ∎

B.6 Theorem.6

Since the sentences in $ML(\mathcal{L}_{1})$ produce a scalar value which can be reached in the graph readout layer as a sum thanks to ${\bf 1}^{\top}H^{(l_{end})}$, we need to show that the MPNN can produce all possible vectors in $\mathcal{L}_{1}$ on the last node representation layer. Since $H^{(0)}={\bf 1}$, the output of the first layer consists of linear combination of $[{\bf 1},A{\bf 1}]$ because, in this case, the third term of the sum is just ${\bf 1}\circ{\bf 1}={\bf 1}$. On the second layer, the representation consists of a linear transformation of 4 different vectors $[{\bf 1},A{\bf 1},A^{2}{\bf 1},A{\bf 1}\circ A{\bf 1}]$. We can notice that these 4 vectors are the all possible vectors that $\mathcal{L}_{1}$ can produce up to the second level. The $diag$ operator can produce other outputs if we apply $diag(A{\bf 1}).diag(A{\bf 1}){\bf 1}=A{\bf 1}\circ A{\bf 1}$. Because $diag({\bf 1})=I$ cannot change anything if we use it any other expressions. Another selection would be $A.diag(A{\bf 1}){\bf 1}=A^{2}{\bf 1}$ and last option gives $diag(A{\bf 1})A{\bf 1}=A{\bf 1}\circ A{\bf 1}$. So up to $l=2$ the proof is true. Then, we follow an inductive reasoning and assume that in the $k$-th layer, Eq.(2) produces all possible vectors ($h_{1},\dots h_{n}$) in $\mathcal{L}_{1}$ and we show that it is true for $k+1$-th layer as well. In the $k+1$-th layer, the first term of the sum keeps $h_{1},\dots h_{n}$. The second term produces $Ah_{1},\dots Ah_{n}$. Finally, the term of the sum produces all pairs of element-wise multiplication such as $h_{1}\circ h_{1},h_{1}\circ h_{2},\dots h_{n}\circ h_{n}$. These are the all vectors that the language $\{.,{\bf 1},diag\}$ can produce using one extra $A$ and/or $diag$ operator. The transpose operator is neglected because the adjacency matrix is symmetric. Furthermore, since at the readout layer these vectors are to be summed up, their order or the fact that they are transposed or not does not matter.

Beside, it was also shown that $diag(.)$ operator can be implemented by element-wise multiplication of vectors in (Geerts, 2020) in Proposition 8.1. ∎

B.7 Theorem.7

If the given function is $\Phi(\lambda)$, it can be written by power series using the Maclaurin expansion as follows:

Thus, the frequency response can be written by power series with coefficients $\alpha_{i}=\frac{\Phi^{(i)}(0)}{i!}$. Using these coefficients, the convolution support can be formulated as

Since $UIU^{\top}=I=L^{0}$ and $Udiag(\lambda)^{n}U^{\top}=L^{n}$, we can reach the final expression:

The convolution support $C$ is expressed as power series of graph laplacian $L$ as long as all order derivation of frequency response is not zero ($\Phi^{(n)}(0)\neq 0$). Since the selection of the function is based on $exp(.)$ and its derivation is never null, we can conclude that designed convolution support can be written by power series of graph Laplacian. ∎

Figure 4, shows Decalin and Bicyclopentyl graphs, with a proposed node enumeration. According to these enumerations, their adjacency matrices are $A_{G}$ and $A_{H}$, respectively

Their normalized Laplacian can be calculated by $L=I-D^{-1/2}AD^{-1/2}$ and gives $L_{G}$ and $L_{H}$ as follows:

$\scriptstyle{L_{G}}=\left(\begin{smallmatrix}1&-0.33&-0.41&0&0&0&-0.41&0&0&0\\ -0.33&1&0&0&0&-0.41&0&0&0&-0.41\\ -0.41&0&1&-0.5&0&0&0&0&0&0\\ 0&0&-0.5&1&-0.5&0&0&0&0&0\\ 0&0&0&-0.5&1&-0.5&0&0&0&0\\ 0&-0.41&0&0&-0.5&1&0&0&0&0\\ -0.41&0&0&0&0&0&1&-0.5&0&0\\ 0&0&0&0&0&0&-0.5&1&-0.5&0\\ 0&0&0&0&0&0&0&-0.5&1&-0.5\\ 0&-0.41&0&0&0&0&0&0&-0.5&1\\ \end{smallmatrix}\right)$

$\scriptstyle{L_{H}}=\left(\begin{smallmatrix}1&-0.33&-0.41&0&0&-0.41&0&0&0&0\\ -0.33&1&0&0&0&0&-0.41&0&0&-0.41\\ -0.41&0&1&-0.5&0&0&0&0&0&0\\ 0&0&-0.5&1&-0.5&0&0&0&0&0\\ 0&0&0&-0.5&1&-0.5&0&0&0&0\\ -0.41&0&0&0&-0.5&1&0&0&0&0\\ 0&-0.41&0&0&0&0&1&-0.5&0&0\\ 0&0&0&0&0&0&-0.50&1&-0.5&0\\ 0&0&0&0&0&0&0&-0.50&1&-0.5\\ 0&-0.41&0&0&0&0&0&0&-0.5&1\\ \end{smallmatrix}\right)$

Their second Chebnet convolution supports are $C^{(2)}_{G}=2/2L_{G}-I$ and $C^{(2)}_{H}=2/1.8418L_{H}-I$ because their maximum eigenvalues are 2.0 and 1.8418 respectively. Finally, when computing the output of the first layer by linear activation function without any learning parameters, we obtain $y_{G}={\bf 1}^{\top}C^{(2)}_{G}{\bf 1}=-9.9327$ and $y_{H}={\bf 1}^{\top}C^{(2)}_{H}{\bf 1}=-9.9269$. We observe a slight difference between these two values, which means that Chebnet can project both graphs to the different points, thus it is able to distinguish them.

Since the maximum eigenvalues of graphs Laplacians are different, they are not cospectral as well. It means that they can also be distinguished on the basis of the number closed walks for some lengths which can be determined by trace operator. Indeed, even if up to 4th power of the adjacency matrix, the trace operator gives the same values for both graphs, we can observe that $tr(A_{G}^{5})=0$ whereas $tr(A_{H}^{5})=20$. This observation is sufficient to claim that both graphs are not $\mathcal{L}_{2}$ equivalent.

Figure 5 shows two non-isomorphic but $\mathcal{L}_{2}$ equivalent graphs, where vertices are enumerated.

According to these enumerations, their adjacency matrices are the following:

We have seen that their normalized Laplacian eigenvalues are $\lambda_{G}=\lambda_{H}=[0,0.44,0.61,0.75,1.25,1.25,1.25,1.25,1.56,1.64]$. Thus, they are cospectral. Considering that for cospectral graphs, the trace of any power of the adjacency matrix which gives the number of closed walks, is the same, we conclude that the trace operator does not help to distinguish these two graphs.

For instance, it can be verified that the trace of the adjacency matrix up to its 5th power is equal: $tr(A_{G}^{2})=tr(A_{H}^{2})=40$, $tr(A_{G}^{3})=tr(A_{H}^{3})=48$, $tr(A_{G}^{4})=tr(A_{H}^{4})=360$, and $tr(A_{G}^{5})=tr(A_{H}^{5})=920$).

However, the sentence $e(X)={\bf 1}^{\top}((X\odot X^{2})^{2}{\bf 1})^{2}$ which implements the element-wise multiplication from $\mathcal{L}_{3}$ allows to distinguish both graphs. Indeed, the computation of this sentences on $A_{G}$ and $A_{H}$ gives ${\bf 1}^{\top}((A_{G}\odot A_{G}^{2})^{2}{\bf 1})^{2}=6032$ and ${\bf 1}^{\top}((A_{H}\odot A_{H}^{2})^{2}{\bf 1})^{2}=5872$. Thus, these two graphs are not $\mathcal{L}_{3}$ equivalent (it means not 3-WL or 2-FWL equivalent as well) because the sample sentence can be explained in $\mathcal{L}_{3}$.

Strongly regular graphs are known to be 3-WL equivalent and $\mathcal{L}_{3}$ equivalent as well. Figure 6 shows sample non-isomorphic graphs that are $\mathcal{L}_{3}$ equivalent.

When we enumerate the nodes from the top-left to the bottom-right according to their locations in the Figure 6, their adjacency matrices are the following:

The eigenvalues of the normalized Laplacian are equal ($\lambda_{G}=\lambda_{H}$). Both normalized Laplacians have 3 distinct eigenvalues which are 0, 0.667 and 1.33 with the respective multiplicity of 1, 6 and 9. Thus the graphs are cospectral. Since they are 3-WL equivalent, none of the sentences in $\mathcal{L}_{3}$ can distinguish these graphs. For instance, we have seen that ${\bf 1}^{\top}((A_{G}\odot A_{G}^{2})^{2}{\bf 1})^{2}={\bf 1}^{\top}((A_{H}\odot A_{H}^{2})^{2}{\bf 1})^{2}=331776$.

Appendix F Result of TU Datasets

Table 5 shows the results of 10-fold cross validation over studied datasets named MUTAG, ENZYMES, PROTEINS and PTC. All these datasets consist of chemical molecules where nodes refer to atoms while edges refer to atomic bonds. For these molecular datasets, node features is a one hot coding of atom types and none of the model use any edge feature even if it exists for MUTAG. In addition to these results, we also provide results on the ENZYMES dataset using extra 18-length continuous features on atoms. Using these continuous features, graph agnostic method MLP performance increases drastically from 30.8% to 70.6%, showing that these continuous features contain at least a part of the structural information. Models were ran for a fixed number of epochs on each fold and we select the epoch where the general accuracy is maximum on the validation set. The test procedure and train/validation split was taken from (Xu et al., 2019).

Appendix G Datasets and Application Details

Table 6 shows the summary of the dataset used in experimental evaluation. The evaluation has been performed on four differents tasks depending on the dataset. These are graph isomorphism (Iso), graph regression (Reg), node regression (NReg) and $n$-class graph classification task (#-Class). We did not use any edge features even if some were available. All features were defined on nodes. These features were discrete node labels coded by one-hot vectors (#Label) and/or continuous features referred by numbers in Tab. 6. We can notice that some graphs have no feature on nodes.

We get the Graph8c and Sr25 dataset from online sourceshttp://users.cecs.anu.edu.au/$\sim$bdm/data/graphs.html, EXP dataset from (Abboud et al., 2020), Random graph dataset from (Chen et al., 2020), 2D-Grid and Band-Pass dataset from (Balcilar et al., 2021), Zinc12K from (Dwivedi et al., 2020), Mnist-75 dataset from online sourcehttps://graphics.cs.tu-dortmund.de/fileadmin/ls7-www/misc/cvpr/mnist-superpixels.tar.gz which was used in (Balcilar et al., 2021) with exactly the same procedure, PROTEINS, ENZYMES, MUTAG and PTC from TU dataset (Morris et al., 2020) downloaded from resources of (Xu et al., 2019). All dataset except for EXP, Random and 2-D grid graph were used on a single task. We used EXP for graph isomorphism test and binary classification task. 2D-Grid graph was used for three different node regression tasks respectively on low-pass, band-pass and high-pass filtering effect prediction. Finally, Random graph is used on five different substructure counting tasks.

In all cases, we used roughly 30K trainable parameters for all problems and all models. We tuned the number of layers from 2 to 5 and the number of convolution kernels in Chebnet from 3 to 5. We used Adam optimization with learning rate in $[10^{-2},10^{-3}]$ and a weight decay in $[10^{-3},10^{-4},10^{-5}]$. We also used dropout layer before all graph convolution layers under selection of $[0,0.1,0.2]$ dropout rate. We used ReLU as non-linearity operation in all layers if it is not mentioned explicitly for any specific model. For classification problems, the loss function was implemented through cross-entropy. For regression problems, mean squared error was used as the loss function except on Zinc12K dataset where the loss function was mean absolute error. Unless otherwise specified, we used both sum and max readout layer after last layer of graph convolution. It is then followed by a fully connected layer which ended up with output layer.

Mentioned hyperparamters are optimzed for concerned model according to validation set performance if it is available. For TU dataset, since the validation and test set is not available in public split, we first created a hyperparameter tuning task by dividing the dataset one time into pre-training (80%) and pre-validation (20%). The optimal value of the parameters is searched on the basis of the performance on the pre-validation set. Then, these hyperparameter values for the general test procedure as defined in (Xu et al., 2019).

Our tests were conducted with implementations of Chebnet, GCN, GIN and GAT layer provided by pytorch-geometric (Fey & Lenssen, 2019). Besides, PPGN, GNNML1 and GNNML3 layer were implemented as a class of pytorch-geometric and our models were tested on the basis of these implementation. By doing so, we integrate the PPGN into the widely used graph library pytorch-geometric and make it publicly available beside our own proposals.

Appendix H Summary of the Baseline Models

In this section of the appendix, we present the baseline methods which are GCN, GIN, Chebnet and GAT thanks to the general framework given by Eq.(1). Each model differs from others by selection of their convolution support $C$.

GCN uses a single convolution support given by;

where $D$ is the diagonal degree matrix (Kipf & Welling, 2017) in Eq.(1).

Chebnet relies on the approximation of a spectral graph analysis proposed in (Hammond et al., 2011), based on the Chebyshev polynomial expansion of the scaled graph Laplacian. The number of convolution supports $C^{(k)}$ can be chosen. They are defined by (Defferrard et al., 2016) as follows:

Graph Isomorphism Network (GIN) defined in (Xu et al., 2019) has a single convolution support defined as follows:

where $\epsilon$ is a parameter that makes the support trainable. Another version named GIN-0 is also defined in the same paper where $\epsilon=0$, which makes $C=A+I$. GIN proposes to use a desired number of MLP after each graph convolution. In our implementation, we use one MLP ($C=I$) after each GIN graph convolution as described in (Xu et al., 2019).

Graph attention networks (GATs) in (Veličković et al., 2018) proposes to transpose the attention mechanism from (Vaswani et al., 2017) into the graph world by the way of sparse attention instead of full attention in transformers. GAT convolution support can be seen as weighted, self loop added adjacency. It can be represented in Eq.(1) by defining its trainable convolution supports as follows:

All MPNN baselines start with a given node features $H^{(0)}$ and provide the node representation of the next layer by Eq.(1). After the last layer, we apply a graph readout function which summarizes the learned node representation. Graph readout layer is followed by a desired number of fully connected layers ended with a number of neuron defined by targeted number of classes.

H.2 PPGN Baseline

PPGN (Maron et al., 2019a) starts the process with a 3-dimensional input tensor where the adjacency, edge features (if it exists) and diagonalized node features are stacked on the 3rd dimension as:

One layer forward calculation of PPNN would be: