Muhammet Balcilar

Pierre Héroux

Benoit Gaüzère

Pascal Vasseur

Sébastien Adam

Paul Honeine

Benoit Ga Üzère

Breaking the Limits of Message Passing Graph Neural Networks

come   L'archive ouverte pluridisciplinaire

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) [START_REF] Gilmer | Neural message passing from quantum chemistry[END_REF]. These algorithms spread each node features to the neighborhood nodes using train-able weights. These weights can be shared with respect to the distance between nodes (Chebnet GNN) [START_REF] Defferrard | Convolutional neural networks on graphs with fast localized spectral filtering[END_REF], to the connected nodes features (GAT for graph attention network) [START_REF] Veličković | Graph attention networks[END_REF] and/or to edge features [START_REF] Bresson | Residual gated graph convnets[END_REF]. 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 [START_REF] Dwivedi | Benchmarking graph neural networks[END_REF].

Despite these successes and these interesting computational properties, it has been shown that MPNNs are not powerful enough [START_REF] Xu | How powerful are graph neural networks[END_REF]. Considering two nonisomorphic 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 [START_REF] Maron | On the universality of invariant networks[END_REF][START_REF] Keriven | Universal invariant and equivariant graph neural networks[END_REF], 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 [START_REF] Morris | Weisfeiler and leman go neural: Higher-order graph neural networks[END_REF], 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 [START_REF] Abboud | The surprising power of graph neural networks with random node initialization[END_REF][START_REF] Sato | Random features strengthen graph neural networks[END_REF] or to add a unique label [START_REF] Murphy | Relational pooling for graph representations[END_REF] 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. [START_REF] Bouritsas | Improving graph neural network expressivity via subgraph isomorphism counting[END_REF][START_REF] Dasoulas | Coloring graph neural networks for node disambiguation[END_REF] 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 [START_REF] Battaglia | Relational inductive biases, deep learning, and graph networks[END_REF]. These may explain why theoretical powerful methods cannot outperform MPNNs on many downstream tasks, as reported in [START_REF] Dwivedi | Benchmarking graph neural networks[END_REF]. 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 [START_REF] Arvind | On weisfeiler-leman invariance: subgraph counts and related graph properties[END_REF][START_REF] Chen | Can graph neural networks count substructures?[END_REF][START_REF] Bouritsas | Improving graph neural network expressivity via subgraph isomorphism counting[END_REF][START_REF] Vignac | Building powerful and equivariant graph neural networks with structural message-passing[END_REF], which are informative for some social or chemical graphs. Finally, another important aspect mentioned by a recent paper [START_REF] Balcilar | Analyzing the expressive power of graph neural networks in a spectral perspective[END_REF] 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

Let G be a graph with n nodes and an arbitrary number of edges. Connectivity is given by the adjacency matrix A ∈ {0, 1} n×n and features are defined on nodes by X ∈ R n×f0 , with f 0 the length of feature vectors. For any matrix X, we used X i , X :j and X i,j to refer to its i-th column vector, j-th row vector and (i, j)-th entry, respectively. A graph Laplacian is given by

L = D -A (or L = I - D -1/2 AD -1/2
) where D ∈ R n×n is the diagonal degree matrix and I is the identity. Through an eigendecomposition, L can be written by L = U diag(λ)U T where each column of U ∈ R n×n is an eigenvector of L, λ ∈ R n gathers the eigenvalues of L and diag(•) creates a diagonal matrix whose diagonal elements are from a given vector. We use superscripts to refer to vectors or matrices evolving through iterations or layers. For instance, H (l) ∈ R n×f l refers to the node representation on layer l whose feature dimension is f l .

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 [START_REF] Wu | A comprehensive survey on graph neural networks[END_REF][START_REF] Chami | Machine learning on graphs: A model and comprehensive taxonomy[END_REF]. Spectral GNNs rely on the spectral graph theory [START_REF] Chung | Spectral graph theory[END_REF]. In this framework, signals on graphs are filtered using the eigendecomposition of the graph Laplacian [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. By transposing the convolution theorem to graphs, the spectral filtering in the frequency domain can be defined by x f lt = U diag(Ω(λ))U x, where Ω(.) is the desired filter function which needs to be learnt by backpropagation. On the other hand, spatial GNNs, such as GCN (graph convolutional network) [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF] and GraphSage [START_REF] Hamilton | Inductive representation learning on large graphs[END_REF], consider two operators, one that aggregates the connected nodes messages and one that updates the concerned node representation.

In a recent paper [START_REF] Balcilar | Analyzing the expressive power of graph neural networks in a spectral perspective[END_REF], it was explicitly shown that both spatial and spectral GNNs are MPNN, taking the general form

H (l+1) = σ s C (s) H (l) W (l,s) , (1)
where C (s) ∈ R n×n is the s-th convolution support that defines how the node features are propagated to the neighboring nodes and W (l,s) ∈ R f l ×f l+1 is the trainable matrix for the l-th layer and s-th support. Within this generalization, GNNs differ from each other by the design of the convolution supports C (s) . If the supports are designed in the spectral domain by Φ s (λ), the convolution support needs to be written as

C (s) = U diag(Φ s (λ))U .
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 nonisomorphic graphs is a universal approximator. Since the graph isomorphism problem is NP-intermediate [START_REF] Takapoui | Linear programming heuristics for the graph isomorphism problem[END_REF], 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 [START_REF] Arvind | On weisfeiler-leman invariance: subgraph counts and related graph properties[END_REF] that for k ≥ 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 ≥ 2, we have (k + 1)-WL ≈ (k)-FWL (Maron et al., 2019a).

In literature, some confusions occur among the two versions. Some papers use WL test order [START_REF] Morris | Weisfeiler and leman go neural: Higher-order graph neural networks[END_REF]Maron et al., 2019a), while others use FWL order under the name of WL such as in [START_REF] Abboud | The surprising power of graph neural networks with random node initialization[END_REF][START_REF] Arvind | On weisfeiler-leman invariance: subgraph counts and related graph properties[END_REF][START_REF] Takapoui | Linear programming heuristics for the graph isomorphism problem[END_REF]. 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 [START_REF] Immerman | Describing graphs: A firstorder approach to graph canonization[END_REF][START_REF] Barceló | The logical expressiveness of graph neural networks[END_REF]. 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 ≡ k-W L A H , if they are indistinguishable by a k-WL (or k-FWL) test.
Recently [START_REF] Brijder | On the expressive power of query languages for matrices[END_REF][START_REF] Geerts | On the expressive power of linear algebra on graphs[END_REF] 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. Definition 1. M L(L) is a matrix language with an allowed operation set L = {op 1 , . . . op n }, where op i ∈ {., +, , diag, tr, 1, , ×, 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. Definition 2. e(X) ∈ R is a sentence in M L(L) if it consists of any possible consecutive operations in L, operating on a given matrix X and resulting in a scalar value.

As an example, e(X) = 1 X 2 1 is a sentence of M L(L) with L = {., , 1}, computing the sum of all elements of square matrix X. In the following, we are interested in languages L 1 , L 2 and L 3 that have been used for characterizing the WL-test in [START_REF] Geerts | On the expressive power of linear algebra on graphs[END_REF]. These results are given next. Remark 1. Two adjacency matrices are indistinguishable by the 1-WL test if and only if e(A G ) = e(A H ) for all e ∈ L 1 with L 1 = {., , 1, diag}. Hence, all possible sentences in L 1 are the same for 1-WL equivalent adjacency matrices. Thus, [START_REF] Geerts | On the expressive power of linear algebra on graphs[END_REF])) Remark 4. Enriching the operation set to L + = L ∪ {+, ×, f } where L ∈ (L 1 , L 2 , L 3 ) does not improve the expressive power of the language. Thus, [START_REF] Geerts | On the expressive power of linear algebra on graphs[END_REF]))

A G ≡ 1-W L A H ↔ A G ≡ M L(L1) A H . (see Theorem 7.1 in (Geerts, 2020)) Remark 2. M L(L 2 ) with L 2 = {., , 1, diag,
= {., , 1, diag, tr, }. Thus, A G ≡ 3-W L A H ↔ A G ≡ M L(L3) A H . (see Theo- rem 9.2 in
A G ≡ M L(L) A H ↔ A G ≡ M L(L + ) A H . (see Proposition 7.5 in

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 [START_REF] Geerts | On the expressive power of linear algebra on graphs[END_REF] and more precisely the remarks of the preceding section. Proofs of the theorems are given in Appendix B. Theorem 1. MPNNs such as GCN, GAT, GraphSage, GIN (defined in Appendix H) cannot go further than operations in L + 1 . Thus, they are not more powerful than the 1-WL test.

This result has already been given in [START_REF] Xu | How powerful are graph neural networks[END_REF], which proposed GIN-(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.

Theorem 2. 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. Theorem 3. 3-star graphlets can be counted by sentences in L + 1 . Theorem 4. Triangle and 4-cycle graphlets can be counted by sentences in L + 2 . Theorem 5. Tailed triangle graphlets can be counted by sentences in 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. [START_REF] Dehmamy | Understanding the representation power of graph neural networks in learning graph topology[END_REF] 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 [START_REF] Geerts | On the expressive power of linear algebra on graphs[END_REF] to break the limits of 1-WL and reach 3-WL equivalence experimentally. GNNML1 relies on the node update schema given by :

H (l+1) =σ(H (l) W (l,1) +AH (l) W (l,2) +H (l) W (l,3) H (l) W (l,4) )
(2) 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: Theorem 6. GNNML1 can produce every possible sentences in M L(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 L + 1 language is needed. More precisely, according to Remarks 2 and 3, trace (tr) and element-wise multiplication ( ) 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 5 G ) = 0 but tr(A 5 H ) = 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 [START_REF] Geerts | On the expressive power of linear algebra on graphs[END_REF]).

In such cases, element-wise multiplication is useful. As an example, the sentence e(A) = 1 f ((A A 2 ) 2 1) where f (x) = x x for any vector x, gives e(A G ) = 6032 and e(A H ) = 5872 for the graphs of Figure 3. Thus, elementwise multiplication helps distinguishing these two graphs.

The calculation details can be found in Appendix D.

As shown by these examples, a model enriched by elementwise 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;[START_REF] Morris | Weisfeiler and leman go neural: Higher-order graph neural networks[END_REF], 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) = 1, a MPNN can calculate C 3 1 by 3 layered MPNN computing C(C(C1))
but not by (C 3 )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 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 A s ) l where l ∈ [0, l max ] is the layer number and s ∈ [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 GN-NML3 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.

Theorem 7. A convolution support given by

C (s) = U diag(Φ s (λ))U , (3)
where

Φ s (λ) = exp(-b(λ -f s ) 2 ), f s ∈ [λ min , λ 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 α s,i =

Φ (i) s (0) i! : C (s) = α s,0 L 0 + α s,1 L 1 + α s,2 L 2 + . . . . (4)
Since design parameters f s of each matrix are different, each C (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:

C = M mlp 4 (mlp 1 (C )|mlp 2 (C ) mlp 3 (C )) , (5)
where S) ] ∈ R n×n×S sparsify convolution support by defined receptive field mask M . The forward calculation of one layer MPNN becomes:

C = [C (1) | . . . |C (s) ] ∈ R n×n×S stacks initial convolution supports on the third dimension, mlp k (.) is a trainable model performed on third dimension of a given ten- sor, C = [C (1) | . . . |C (
H (l+1) = σ s (C (s) H (l) W (l,s) )|mlp5(H (l) ) mlp6(H (l) ) (6)
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 Φ s (λ) 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

B = I -D -1/2 AD -1/2 or B = A. Eigendecomposition: U diag(λ)U = B for s = 1 to S do C :s = sparse2vec M (U diag(Φs(λ))U ) end for Algorithm 2 GNNML3 Forward calculation Input: extracted edge features C ∈ R m×S , initial node feature H (0) ∈ R n×f 0 , receptive field mask M ∈ {0, 1} n×n , number of layers L, number of supports S Output: new node representation H (L) for l = 0 to L -1 do C=mlp l,4 (mlp l,1 (C )|mlp l,2 (C ) mlp l,3 (C )) for s = 1 to S do C (s) =vec2sparse( C:s,M )
end for

H (l+1) =σ( s (C (s) H (l) W (l,s) )|mlp l,5 (H (l) ) mlp l,6 (H (l) ))
end for (f s ) by given bandwidth (b). Therefore, our model will be able to learn properties depending on the spectrum of graph signal.

Algorithm 1 calculates the initial convolution supports (C ). Since the supports are computed for valid indices in the receptive field mask (where M i,j = 1), one can see C as extracted edge features where the edge indices are defined by M . In application, a function sparse2vec : R n×n → R m converts the sparse matrix to a vector by just keeping the components on valid indices of the mask. Algorithm 2 shows the forward calculation of the model for just one graph. To make the representation as simple as possible, we prefer to use tensor representation in Eq.( 5). However, implementation of Algorithm 2 just apply mlp k (.) to the valid indices defined by receptive field mask M . Thus, C, C have the dimension of R m×S , where m shows number of valid indices in M and mlp k (.) applies on columns of C . Beside, we use a function vec2sparse : R m → R n×n that converts the vector to the sparse convolution support according to a given mask M .

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 [START_REF] Arvind | On weisfeiler-leman invariance: subgraph counts and related graph properties[END_REF]. In Appendix E, a strongly regular graphs pair and the result of a sample sentence in 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 online1 . 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 bandpass 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 datasets2 . 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 nonconnected nodes share 6 common neighbors. Sr25 consists of 15 graphs, leading to 105 different pairs for comparison.

Moreover, we use the EXP dataset [START_REF] Abboud | The surprising power of graph neural networks with random node initialization[END_REF], 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. 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. λ 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

4-CYCLES MLP 1.0E-4 4.58E-1 3.13E-1 2.22E-1 1.73E-1 GCN 1.0E-4 3.22E-3 2.43E-1 1.42E-1 1.14E-1 GAT 1.0E-4 4.57E-3 2.47E-1 1.44E-1 1.12E-1 GIN 1.0E-4 1.47E-3 2.06E-1 1.18E-1 1.21E-1 CHEBNET 1.0E-4 7.68E-4 2.01E-1 1.15E-1 9.60E-2 PPGN 1.0E-4 9.19E-4 1.00E-4 2.61E-4 3.30E-4 GNNML1 1.0E-4 2.75E-4 2.45E-1 1.32E-1 1.14E-1 GNNML3 1.0E-4 7.24E-4 4.44E-4 3.18E-4 6.62E-4
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 ∈ L + 1 , e c (A) = 1 A diag(exp(-A 2 1))A1 with A the graph adjacency matrix. Since e c ∈ M L(L + 1 ), it may be learnable by 1-WL equivalent MPNNs. We used the RandomGraph dataset [START_REF] Chen | Can graph neural networks count substructures?[END_REF] 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 [START_REF] Balcilar | Analyzing the expressive power of graph neural networks in a spectral perspective[END_REF], 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 [START_REF] Balcilar | Analyzing the expressive power of graph neural networks in a spectral perspective[END_REF]. 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. 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 [START_REF] Bouritsas | Improving graph neural network expressivity via subgraph isomorphism counting[END_REF]. Obtained results confirm that models that are able to count substructures, such as PPGN and GNNML3,

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.

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 [START_REF] Takapoui | Linear programming heuristics for the graph isomorphism problem[END_REF]. 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) [START_REF] Arvind | On weisfeiler-leman invariance: subgraph counts and related graph properties[END_REF].

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

(0) v = 1
). Then, colors are updated by the following iteration:

H (t+1) v = σ H (t) v | H (t) u : u ∈ N (v) , (7)
where

H (t)
v is the color of vertex v at iteration t, N (v) is the set of neighbours of vertex v, | represents the concatenation operator and {.} is the order invariant multiset 3 . 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 σ(•) 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.

Higher order WL tests use the same algorithm while their color update schema is slightly different. The 2-WL test uses second order tuple of nodes (all ordered pairs of nodes), thus it needs H ∈ R n×n matrix, where the initial color set has two more colors than initial vertex colors as defined by:

H (0) v,u =    H (0) v if v = u edge if u ∈ N (v) nonedge if u ∈ N (v) (8)
Then, the iteration process is applied through the following 3 It is generally implemented by stacking all colors in the set and sorting them alphabetically schema where [n] is the set of node identifiers.

H (t+1) v,u = σ H (t) v,u | H (t) v,k : k ∈ [n] | H (t) k,u : k ∈ [n] , (9)
Although for k ≥ 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 ≥ 2, we have (k + 1)-WL ≈ (k)-FWL (Maron et al., 2019a). The iteration process of 2-FWL is given by the following equation;

H (t+1) v,u = σ H (t) v,u | H (t) v,k |H (t) k,u : k ∈ [n] , (10)
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) [START_REF] Morris | Weisfeiler and leman go neural: Higher-order graph neural networks[END_REF]Maron et al., 2019a) but some others such as [START_REF] Abboud | The surprising power of graph neural networks with random node initialization[END_REF][START_REF] Arvind | On weisfeiler-leman invariance: subgraph counts and related graph properties[END_REF][START_REF] Takapoui | Linear programming heuristics for the graph isomorphism problem[END_REF] 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.

B. Proofs of Theorems

B.1. Theorem.1

Proof. 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 L + 1 (Remark 4), Eq.(1) also can be explained by operations from L + 1 as well. Thus these methods cannot produce any sentence out of 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 L + 1 .

GCN uses C = (D + I) -0.5 (A + I)(D + I) -0.5 where D is the diagonal degree matrix [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF] in Eq.( 1). (D + I) -0.5 can be expressed as (D + I) -0.5 = diag(f (A1 + 1)), where f (x) = x -0.5 is element-wise operation on vector x. A + I can also be written A + diag(1). When we merge these equations, we get

C = diag(f (A1 + 1))(A + diag(1))diag(f (A1 + 1)).
The convolution support C is then written using operations from L + 1 . In the literature, GraphSage method was proposed to sample neighborhood and aggregate the neighborhood contribution by the mean operator or LSTM in [START_REF] Hamilton | Inductive representation learning on large graphs[END_REF]. 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 L + 1 , by observing that C (1) = diag(1) and C (2) = diag(f (A1))A, where f (x) = x -1 elementwise operation on vector x.

GIN [START_REF] Xu | How powerful are graph neural networks[END_REF] uses a convolution support C = A+I 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 L + 1 . C = A + × diag(1) and C mlp = diag(1). GAT [START_REF] Veličković | Graph attention networks[END_REF] can be expressed in Eq.( 1) by the convolution support designed by

C v,u = m(H v , H u )/ k∈ Ñ (v) m(H v , H k ), where Ñ (v) is the self-connection
(H v , H u ) = f 1 (H v ) + f 2 (H u ), we can define an intermediate ma- trix B = diag(f 1 (H))(A + I) + (A + I)diag(f 2 (H)).
Finally the GAT convolution support can be written by C = diag((B1) -1 )B using all operations included within the operation set L + 1 .

B.2. Theorem.2

Proof. Chebnet [START_REF] Defferrard | Convolutional neural networks on graphs with fast localized spectral filtering[END_REF] uses desired number k of convolution supports in Eq.( 1). As long as these convolutions can be written by operations in + 1 , we can conclude that Chebnet is no more powerful than 1-WL test. But if at least one convolution cannot be explained in L + 1 , we can say it is more powerful than 1-WL test.

Chebnet's convolution supports are -2) . The first support can always be written thanks to an operation from L 1 since C (1) = diag(1). Both normalized and combinatorial graph Laplacian can also be written as L = diag(A1)-A or L = diag(1) -diag(f (A1))Adiag(f (A1)) where f (x) = x -1/2 elementwise operation on vector x. If λ max for both graphs are the same, we can use a constant α = 2/λ max . The second convolution support can then be written as -2) are created by matrix multiplication and subtraction of previous supports which can all be expressed by mean of operations from L + 1 . Thus, if the maximum eigenvalues of tested graphs Laplacians are the same, Chebnet is not more powerful than 1-WL.

C (1) = I, C (2) = 2L/λ max -I, C (k) = 2C (2) C (k-1) -C (k
C (2) = α × L -diag(1). It is then expressed by means of operations from L + 1 . Other convolution supports C (k) = 2C (2) C (k-1) -C (k
However, if the maximum eigenvalues are not the same, C (2) cannot be expressed with the help of the constant value α. 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) andC (2) 1). If these two graphs are 1-WL equivalent, any sentence build on L + 1 applied on these graph is equivalent as well. For instance, we can use the sentences of e(X) = 1 X1 with operation in L + 1 . The output of the sentence should be same such e(L G ) = e(L H ) yields 1 L G 1 = 1 L H 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

G = α G × L G -diag(1)
H = α H × L H -diag(
(2)

G ) = e(C (2) H ) yields α G 1 L G 1 = α H 1 L H 1.
Last equation has contradiction to the previous one as long as the maximum eigenvalues are not same (i.e α G = α H ) and graphs are not regular (i.e 1 L G 1 > 0 and 1 L H 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 1 L G 1 ≥ 0 and 1 L H 1 ≥ 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 ← L + 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 ← 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 1 L1 = 0.

B.3. Theorem.3

Proof. The number of 3-star patterns can be determined by where d(v) is the degree of vertex v for undirected simple graphs [START_REF] Pinar | Escape: Efficiently counting all 5-vertex subgraphs[END_REF]. Using f (x) =

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 1 f (A1) using operations in L + 1 . According to the universal approximation theory of multi layer perceptron [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF], if we have enough layers, we can implement f (.) as an MLP in our model.

B.4. Theorem.4

Proof. The number of triangles can be determined by using trace operator as 1/6 × tr(A 3 ) [START_REF] Harary | On the number of cycles in a graph[END_REF] which can be written by means of operations from L + 2 . Number of 4-cycles is determined by 1/8 × (tr(A 4 ) + tr(A 2 ) -21 A 2 1) [START_REF] Harary | On the number of cycles in a graph[END_REF] which can be written by means of operations from L + 2 .

B.5. Theorem.5

Proof. 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 v t(v).(d(v) -2) for simple undirected graphs [START_REF] Pinar | Escape: Efficiently counting all 5-vertex subgraphs[END_REF]. Every node in a triangle has two closed walks of length 3. Thus, t(v) =

(A 3 )v,v 2
. It yields the number of tailed triangles can be found by 1 2 × 1 (A 3 diag(A1 -2))1. The computation of t(v) which involves the element-wise multiplication can be written with operations from L + 3 .

B.6. Theorem.6

Proof. Since the sentences in M L(L 1 ) produce a scalar value which can be reached in the graph readout layer as a sum thanks to 1 H (l end ) , we need to show that the MPNN can produce all possible vectors in L 1 on the last node representation layer. Since H (0) = 1, the output of the first layer consists of linear combination of [1, A1] because, in this case, the third term of the sum is just 1 • 1 = 1. On the second layer, the representation consists of a linear transformation of 4 different vectors [1, A1, A 2 1, A1•A1]. We can notice that these 4 vectors are the all possible vectors that L 1 can produce up to the second level. The diag operator can produce other outputs if we apply diag(A1).diag(A1)1 = A1 • A1. Because diag(1) = I cannot change anything if we use it any other expressions. Another selection would be A.diag(A1)1 = A 2 1 and last option gives diag(A1)A1 = A1 • A1. 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 , . . . h n ) in 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 , . . . h n . The second term produces Ah 1 , . . . Ah n . Finally, the term of the sum produces all pairs of elementwise multiplication such as

h 1 • h 1 , h 1 • h 2 , . . . h n • h n .
These are the all vectors that the language {., 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 [START_REF] Geerts | On the expressive power of linear algebra on graphs[END_REF] in Proposition 8.1.

B.7. Theorem.7

Proof. If the given function is Φ(λ), it can be written by power series using the Maclaurin expansion as follows:

Φ(λ) = Φ(0) 0! λ 0 + Φ (0) 1! λ 1 + Φ (2) (0) 2! λ 2 + . . . . (11)
Thus, the frequency response can be written by power series with coefficients

α i = Φ (i) (0) i!
. Using these coefficients, the convolution support can be formulated as

C = α 0 U IU +α 1 U diag(λ)U +α 2 U diag(λ) 2 U +. . . . (12) Since U IU = I = L 0 and U diag(λ) n U = L n ,
we can reach the final expression:

C = α 0 L 0 + α 1 L 1 + α 2 L 2 + . . . (13)
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 (Φ (n) (0) = 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.

C. L 1 Equivalent Graphs

A G =         
0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 1 0

         and A H =         
0 1 1 0 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 1 0

        
Their normalized Laplacian can be calculated by L = I -D -1/2 AD -1/2 and gives L G and L H as follows:

L G =            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            L H =            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           
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 = 1 C (2) G 1 = -9.9327 and y H = 1 C (2)
H 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 5 G ) = 0 whereas tr(A 5 H ) = 20. This observation is sufficient to claim that both graphs are not L 2 equivalent. According to these enumerations, their adjacency matrices are the following:

D. L 2 Equivalent Graphs

A G =         
0 1 0 1 0 1 0 1 0 0 1 0 1 1 1 0 0 0 0 0 0 1 0 0 1 0 1 0 0 1 1 1 0 0 0 1 0 1 0 0 0 1 1 0 0 0 1 0 0 1 1 0 0 1 0 0 0 0 1 1 0 0 1 0 1 0 0 1 1 0 1 0 0 1 0 0 1 0 1 0 0 0 0 0 0 1 1 1 0 1 0 0 1 0 1 1 0 0 1 0

         and A H =         
0 1 0 1 0 0 1 1 0 0 1 0 1 1 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 1 1 1 0 0 0 1 0 1 0 0 0 1 1 0 0 0 1 0 0 1 0 0 1 1 0 0 0 1 1 0 1 0 0 0 1 0 0 0 1 1 1 0 0 1 0 1 0 0 1 0 0 0 0 0 0 1 1 1 0 1 0 0 1 0 1 0 1 0 1 0

        
We have seen that their normalized Laplacian eigenvalues are λ G = λ 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: andtr(A 5 G ) = tr(A 5 H ) = 920). However, the sentence e(X) = 1 ((X X 2 ) 2 1) 2 which implements the element-wise multiplication from L 3 allows to distinguish both graphs. Indeed, the computation of this sentences on A G and A H gives 1 ((A G A 2 G ) 2 1) 2 = 6032 and 1 ((A H A 2 H ) 2 1) 2 = 5872. Thus, these two graphs are not L 3 equivalent (it means not 3-WL or 2-FWL equivalent as well) because the sample sentence can be explained in L 3 .

tr(A 2 G ) = tr(A 2 H ) = 40, tr(A 3 G ) = tr(A 3 H ) = 48, tr(A 4 G ) = tr(A 4 H ) = 360,

E. L 3 Equivalent Graphs

Strongly regular graphs are known to be 3-WL equivalent and L 3 equivalent as well. Figure 6 shows sample nonisomorphic graphs that are 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:

A G =                  
0 1 1 1 1 0 0 0 1 0 0 0 1 0 0 0 1 0 1 1 0 1 0 0 0 1 0 0 0 1 0 0 1 1 0 1 0 0 1 0 0 0 1 0 0 0 1 0 1 1 1 0 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 1 1 1 1 0 0 0 1 0 0 0 0 1 0 0 1 0 1 1 0 1 0 0 0 1 0 0 0 0 1 0 1 1 0 1 0 0 1 0 0 0 1 0 0 0 0 1 1 1 1 0 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 0 1 1 1 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 1 0 1 0 0 0 0 1 0 0 0 1 0 1 1 0 1 0 0 1 0 0 0 0 1 0 0 0 1 1 1 1 0 0 0 0 1 1 0 0 0 1 0 0 0 1 0 0 0 0 1 1 1 0 1 0 0 0 1 0 0 0 1 0 0 1 0 1 1 0 0 1 0 0 0 1 0 0 0 1 0 1 1 0 1 0 0 0 1 0 0 0 1 0 0 0 1 1 1 1 0

A H =                  
0 1 0 1 1 1 0 0 0 0 0 0 1 0 0 1 1 0 1 0 0 1 1 0 0 0 0 0 1 1 0 0 0 1 0 1 0 0 1 1 0 0 0 0 0 1 1 0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 1 1 1 0 0 1 0 1 0 1 1 1 0 0 0 0 0 0 1 1 0 0 1 0 1 0 0 1 1 0 0 0 0 0 0 1 1 0 0 1 0 1 0 0 1 1 0 0 0 0 0 0 1 1 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 1 1 1 0 0 0 0 0 0 1 1 0 0 1 0 1 0 0 1 1 0 0 0 0 0 0 1 1 0 0 1 0 1 0 0 1 1 0 0 0 0 0 0 1 1 1 0 1 0 1 0 0 1 1 1 0 0 0 0 0 0 1 0 0 1 0 1 0 1 0 1 1 0 0 0 0 0 1 1 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 1 1 0 0 1 0 1 1 0 0 1 0 0 0 0 0 0 1 1 1 0 1

0                  
The eigenvalues of the normalized Laplacian are equal (λ G = λ 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 L 3 can distinguish these graphs. For instance, we have seen that

1 ((A G A 2 G ) 2 1) 2 = 1 ((A H A 2 H ) 2 1) 2 = 331776.
In order to distinguish these two graphs, we need to mimic the 3-FWL (or 4-WL) test which needs a 3-order relationship between nodes. Thus, the adjacencies will be represented by A G , A H ∈ R 16×16×16 . For any 3 nodes there are 3 different pairs and thus 2 3 = 8 different states representing how these 3 nodes are connected or not. An additional state is used for the tensor diagonal. Thus, there is a total of 9 states. The node tuple is denoted by A i,j,k ∈ {0, . . . , 8}. 0 refers to the fact that none of three nodes are connected. 7 refers to the fact that all nodes are mutually connected (triangle). 8 is used for the tensor diagonal elements. We can then define an equivariant 3 dimensional tensor square operator by (A 2 ) i,j,k = s (A s,j,k .A i,s,k .A i,j,s ). By summing all elements of the 3-dimensional squared adjacency where the given adjacency is for instance 0, we can distinguish these two graphs. Indeed, (A 2 G (A G = 0)) = 205632 whereas (A 2 H (A H = 0)) = 208704. We can then conclude that these two graphs are not 3-FWL (or 4-WL) equivalent.

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 [START_REF] Xu | How powerful are graph neural networks[END_REF].

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 sources4 , EXP dataset from [START_REF] Abboud | The surprising power of graph neural networks with random node initialization[END_REF], Random graph dataset from [START_REF] Chen | Can graph neural networks count substructures?[END_REF], 2D-Grid and Band-Pass dataset from [START_REF] Balcilar | Analyzing the expressive power of graph neural networks in a spectral perspective[END_REF], Zinc12K from [START_REF] Dwivedi | Benchmarking graph neural networks[END_REF], Mnist-75 dataset from online source5 which was used in [START_REF] Balcilar | Analyzing the expressive power of graph neural networks in a spectral perspective[END_REF] with exactly the same procedure, PROTEINS, ENZYMES, MUTAG and PTC from TU dataset [START_REF] Morris | Tudataset: A collection of benchmark datasets for learning with graphs[END_REF] downloaded from resources of [START_REF] Xu | How powerful are graph neural networks[END_REF]. 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.

In GNNML3, we use the eigendecomposition of normalized Laplacian to calculate the initial edge feature for all problems, except for Zinc12K and substructure counting problems where the eigen decomposition was performed on the adjacency. Each initial convolution support is set such that Φ s (λ) = exp(-b(λ -f s ) 2 ), where the bandwidth parameter b is set to the value of 5. The spectrum has been uniformly sampled between minimum eigenvalue and the maximum eigenvalue with a selection of s n = [3, 5, 10] points in order to select the band specific parameter. Thus, band specific parameter of each frequency profile can be written f s = λ min + s-1 sn-1 (λ max -λ min ) for s ∈ {1, . . . , s n -1}. For the convolution support s = 0, we used all-pass filtering named identity matrix whose frequency response is Φ 0 (λ) = 1. Thus, we have a total of s n convolution supports. The 1-hop distance is always used for receptive field which corresponds to M = A + I. For the learning of convolution supports needed in Eq.( 5), we used a single layered MLP in each mlp k where mlp 1 , mlp 2 , mlp 3 : R S → R 2S with a sigmoid activation, and mlp 4 : R 4S → R S with ReLU activation as long as S is the number of initial convolutions extracted in the preprocessing step. In Eq.( 6), the size of the output of mlp 5 and mlp 6 is another hyperparameter where we used the same length with the first part of the Eq.( 6) defined by dimension of W (l,s) . 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 [START_REF] Xu | How powerful are graph neural networks[END_REF].

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 pytorchgeometric 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.

H. Summary of the Baseline Models

H.1. MPNN Baselines

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;

C = (D + I) -0.5 (A + I)(D + I) -0.5 , ( 14
)
where D is the diagonal degree matrix [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF] in Eq.( 1).

Chebnet relies on the approximation of a spectral graph analysis proposed in [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF], 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 [START_REF] Defferrard | Convolutional neural networks on graphs with fast localized spectral filtering[END_REF] as follows:

C (1) = I, C (2) = 2L/λ max -I,

C (k) = 2C (2) C (k-1) -C (k-2) , ∀k ≥ 2. ( 15
)
Graph Isomorphism Network (GIN) defined in [START_REF] Xu | How powerful are graph neural networks[END_REF]) has a single convolution support defined as follows:

C = A + (1 + )I, (16)
where is a parameter that makes the support trainable. Another version named GIN-0 is also defined in the same paper where = 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 [START_REF] Xu | How powerful are graph neural networks[END_REF].

Graph attention networks (GATs) in [START_REF] Veličković | Graph attention networks[END_REF] proposes to transpose the attention mechanism from [START_REF] Vaswani | Attention is all you need[END_REF] 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:

C (l,s) v,u = e v,u k∈ Ñ (v) e v,k , (17)
where e v,u = exp σ(a (l,s) [H

:v W (l,s) ||H (l)

:u W (l,s) ]) , and a (l,s) is another trainable weight. Convolution support will be calculated from node v to each element of Ñ (v), which shows the self-connection added neighborhood. In application of GAT, we use concatenation instead of sum in Eq.( 1) where the paper proposed both and there is slightly empirical advantage to use concatenation.

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 3dimensional input tensor where the adjacency, edge features (if it exists) and diagonalized node features are stacked on the 3rd dimension as:

H (0) = [A|E 1 | • • • |E e |diag(X 1 )| • • • |diag(X d )]. (18)
Here, X ∈ R n×d gathers node features and X i is its i-th column vector, E ∈ R n×n×e is edge features and E i ∈ R n×n is its i-th edge feature matrix, thus initial feature tensor is H (0) ∈ R n×n×(1+e+d) .

One layer forward calculation of PPNN would be:

H (l+1) = m 3 m 1 (H (l) ) • m 2 (H (l) )|H (l) (19)
where m 1 , m 2 : R n×n×dinp → R n×n×d mid and m 3 : R n×n×d mid +dinp → R n×n×dout are trainable models that can be implemented by a one layer MLP followed by nonlinearity. d inp is the feature length on the 3rd dimension. d mid , d out are the feature lengths which can be seen as hyperparameters of the layer. Multiplication (•) operates between matching features and means 2d matrix multiplication for each slice which has n × n dimensions. | operator is just the concatenation of two tensor on the 3rd dimension.

The output of the model would be: l) ) .

Y = l=1 mlp l diag(H (l) ) | offdiag(H (
(20)

We assign a function which selects the diagonal of each 2d slices of tensor as diag : R n×n×d → R n×1×d and function for selection the element out of the diagonal as offdiag : R n×n×d → R n×(n-1)×d . We use the sum operator which performs sum over the first 2 dimensions as : R d1×d2×d → R d and a trainable model that may be implemented by an MLP mlp l : R 2d → R dy , transforms the given vector into the targeted output representation length.

The one can see that in each layer, PPNN keeps H (l) ∈ R n×n×d l , thus its memory usage is in O(n 2 ). Since there is a matrix multiplication in Eq.( 19), its computation complexity is in O(n 3 ) when using the naive matrix multiplication operations. The PPNN paper mentioned that the computational complexity can be decreased by using effective matrix multiplication, but it is the same for all algorithms as well. For this reason, we think that taking the naive implementation into account makes more sense to do a fair comparison. In addition, again because of matrix multiplication, its update mechanism is not local. Because of calculation of the u, v node pairs representation in Eq.( 19), it needs to perform k H (l) u,k .H (l) k,v . That means that for each pair of nodes, k should be all nodes in the graph regardless how far away the node k from the concerned nodes u, v. In other words, very far away nodes feature affect the concerned node.

Even though PPNN (Maron et al., 2019a) is a very straight forward algorithm and has provable 3-WL power, the experimental results reported in the papers are not at the state of the art (Maron et al., 2019a;[START_REF] Dwivedi | Benchmarking graph neural networks[END_REF]. We believe that this can be at least partly explained by some implementation problems. Indeed, it was implemented by gathering same size graphs into batches in order to handle graphs of different size in a dataset. So the batches do not consist of randomly selected graphs in each epoch during the training phase. In our implementation, we first find the maximum size of the graph denoted as n max . Then, we create an initial tensor in Eq.( 18) in dimension of R nmax×nmax×1+e+d where left top n × n × 1 + e + d part of the tensor is valid, and the rest is zero. We also keep the valid part of the tensor diagonal and out of diagonal part mask in M 0 , M 1 ∈ {0, 1} nmax×nmax that shows which element is valid in the diagonal and which element is valid out of the diagonal of the representation tensor. Since some part of the tensor H (l) are not valid, we need to prevent to assign value after application of trainable model m k in Eq.( 19), because it affects the matrix multiplication result. One solution may be to mask the MLP result by M 0 + M 1 . Finally, we implement Eq.( 20) by selection diagonal and off-diagonal element by previously prepared mask matrices by diag(H (l) ) = M 0 H (l) and offdiag(H (l) ) = M 1 H (l) . By doing so, we can put any graph into same batch. These principles have been implemented as a class of the widely used open-source pytorch geometric library.

Figure 1 .

 1 Figure 1. Decalin (G) and Bicyclopentyl (H) graphs are L1 and also 1-WL equivalent, but Chebnet can distinguish them.

Figure 2 .

 2 Figure 2. Sample of patterns: 3-star, triangle, tailed triangle and 4-cycle graphlets used in our analysis.

Figure 3 .

 3 Figure 3. Cospectral and 4-regular graphs from (Van Dam & Haemers, 2003) are L1 and L2 equivalent.

Figure 4 .

 4 Figure 4. Decalin and Bicyclopentyl graphs are L1 equivalent and so 1-WL.

Figure 4 ,

 4 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

Figure 5

 5 Figure5shows two non-isomorphic but L 2 equivalent graphs, where vertices are enumerated.

Figure 5 .

 5 Figure 5. Cospectral and 4-regular graphs from (Van Dam & Haemers, 2003) are L2 equivalent.

Figure 6 .

 6 Figure 6. Strongly regular graph pair. 4 × 4-rook's graph and the Shrikhande graph from (Arvind et al., 2020) are L3 equivalent.

Table 1 .

 1 Number of undistinguished pairs of graphs in graph8c, sr25 and EXP datasets and binary classification accuracy on EXP dataset. An ideal method does not find any pair similar and classifies graphs with 100% accuracy. The number of pairs is 61M for graph8c, 105 pairs for sr25 and 600 for EXP.

	MODEL	GRAPH8C	SR25	EXP	EXP-CLASSIFY
	MLP	293K	105	600	50%
	GCN	4775	105	600	50%
	GAT	1828	105	600	50%
	GIN	386	105	600	50%
	CHEBNET	44	105	71	82%
	PPGN	0	105	0	100%
	GNNML1	333	105	600	50%
	GNNML3	0	105	0	100%

  Table1presents 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

Table 2 .

 2 Median of test set MSE error for graphlet counting problem on random graph dataset over 10 random runs.

	MODEL	3-STARS	CUSTOM	TRIANGLE	TAILED-TRI

Table 3 .

 3 Spectral expressive analysis results. R 2 for LowPass, HighPass and BandPass node regression tasks, accuracy on graph classification task. Results are median of 10 different runs.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.

	MODEL	LOW-PASS	HIGH-PASS	BAND-PASS	CLASSIFY
	MLP	0.9749	0.0167	0.0027	50.0%
	GCN	0.9858	0.0863	0.0051	77.9%
	GAT	0.9811	0.0879	0.0044	85.3%
	GIN	0.9824	0.2934	0.0629	87.6%
	CHEBNET	0.9995	0.9901	0.8217	98.2%
	PPGN	0.9991	0.9925	0.1041	91.2%
	GNNML1	0.9994	0.9833	0.3802	92.8%
	GNNML3	0.9995	0.9909	0.8189	97.8%

Table 4 .

 4 Results on Zinc12K and MNIST-75 datasets. The values are the MSE for Zinc12K and the accuracy for MNIST-75. Edge features are not used even if they are available in the datasets. For Zinc12K, all models use node labels. For MNIST-75, the model uses superpixel intensive values and node degree as node features.

	MODEL	ZINC12K	MNIST-75
	MLP	0.5869 ± 0.025	25.10% ± 0.12
	GCN	0.3322 ± 0.010	52.80% ± 0.31
	GAT	0.3977 ± 0.007	82.73% ± 0.21
	GIN	0.3044 ± 0.010	75.23% ± 0.41
	CHEBNET	0.3569 ± 0.012	92.08% ± 0.22
	PPGN	0.1589 ± 0.007	90.04% ± 0.54
	GNNML1	0.3140 ± 0.015	84.21% ± 1.75
	GNNML3	0.1612 ± 0.006	91.98% ± 0.18
	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.

  added neighborhood of v and m(.) is any trainable model. If we write the trainable model m(.) as a sum of each node such as m

Table 5 .

 5 Results on TU datasets. The values are the accuracy. Edge features are not used even if they are available in the datasets. The models use a one-hot encoding of node labels as node features, while the models also use extra 18 length continuous node features for ENZYMES-cont.

	MODEL	MUTAG	ENZYMES	ENZYMES-CONT	PROTEINS	PTC
	MLP	86.6% ± 4.95	30.8% ± 4.26	70.6% ± 5.22	74.3% ± 4.88	62.9% ± 5.89
	GCN	89.1% ± 5.81	49.0% ± 4.25	74.2% ± 3.26	75.2% ± 5.11	64.3% ± 8.35
	GAT	90.1% ± 5.84	54.1% ± 5.15	73.7% ± 4.47	75.9% ± 4.26	65.7% ± 7.97
	GIN	89.4% ± 5.60	55.8% ± 5.23	73.3% ± 4.48	76.1% ± 3.97	64.6% ± 7.00
	CHEBNET	89.7% ± 6.41	63.8% ± 7.92	75.3% ± 4.63	76.4% ± 5.34	65.5% ± 4.94
	PPGN	90.2% ± 6.62	55.2% ± 5.44	72.9% ± 4.18	77.2% ± 4.53	66.2% ± 6.54
	GNNML1	90.0% ± 0.42	54.9% ± 5.97	76.9% ± 5.14	75.8% ± 4.93	63.9% ± 6.37
	GNNML3	90.9% ± 5.46	63.6% ± 6.52	78.1% ± 5.05	76.4% ± 5.10	66.7% ± 6.49

Table 6 .

 6 Summary of the datasets used in our experiments.

		GRAPH8C	SR25	EXP	2D-GRID	RANDOM	BAND-PASS	PROTEINS	ENZYMES	MUTAG	PTC	MNIST-75	ZINC12K
	TASK	ISO	ISO	ISO&2CLASS	NREG	REG	2CLASS	2CLASS	6CLASS	2CLASS	2CLASS	10CLASS	REG
	GRAPHS	11117	15	1200	3	5K	5K	1113	600	188	344	70K	12K
	NODES	8.0	25.0	44.44	900.0	18.8	200.0	39.06	32.63	17.93	25.55	75.0	23.15
	EDGES	28.82	300.0	110.21	3480.0	62.67	1072.6	72.82	62.14	39.58	51.92	694.7	49.83
	FEATURE	MONO	MONO	MONO	1	MONO	1	3LABEL	3LABEL+18	7LABEL	19LABEL	1	21LABEL
	TRAIN	NA	NA	800	1	1500	3K	9-FOLD	9-FOLD	9-FOLD	9-FOLD	55K	10K
	VAL	NA	NA	200	1	1000	1K	1-FOLD	1-FOLD	1-FOLD	1-FOLD	5K	1K
	TEST	NA	NA	200	1	2500	1K	NA	NA	NA	NA	10K	1K

LITIS Lab, University of Rouen Normandy, France

InterDigital, France

LITIS Lab, INSA Rouen Normandy, France

MIS Lab, Université de Picardie Jules Verne, France. Correspondence to: Muhammet Balcilar <muhammetbalcilar@gmail.com>.

https://github.com/balcilar/gnn-matlang

http://users.cecs.anu.edu.au/∼bdm/data/graphs.html

http://users.cecs.anu.edu.au/∼bdm/data/graphs.html

https://graphics.cs.tu-dortmund.de/fileadmin/ls7www/misc/cvpr/mnist-superpixels.tar.gz

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.