GEOMETRIC RANDOM WALK GRAPH NEURAL NET-
Year: 2021

Abstract:
Graph neural networks have recently attracted a lot of attention and have been applied with great success to several important graph problems. The Random Walk Graph Neural Network model was recently proposed as a more intuitive alternative to the well-studied family of message passing neural networks. This model compares each input graph against a set of latent “hidden graphs” using a kernel that counts common random walks up to some length. In this paper, we propose a new architecture, called Geometric Random Walk Graph Neural Network (GRWNN), that generalizes the above model such that it can count common walks of infinite length in two graphs. The proposed model retains the transparency of Random Walk Graph Neural Networks since its first layer also consists of a number of trainable “hidden graphs” which are compared against the input graphs using the geometric random walk kernel. To compute the kernel, we employ a fixed-point iteration approach involving implicitly defined operations. Then, we capitalize on implicit differentiation to derive an efficient training scheme which requires only constant memory, regardless of the number of fixed-point iterations. The employed random walk kernel is differentiable, and therefore, the proposed model is end-to-end trainable. Experiments on standard graph classification datasets demonstrate the effectiveness of the proposed approach in comparison with stateof-the-art methods.

heading: 1 INTRODUCTION
text: Recent years have witnessed an enormous growth in the amount of data represented as graphs. Indeed, graphs emerge naturally in several domains, including social networks, bioinformatics, and neuroscience, just to name a few. Besides the increase in the amount of graph-structured data, there is also a growing interest in applying machine learning techniques to data modeled as graphs. Among others, the graph classification and graph regression tasks have attracted a great deal of attention in the past years. These tasks have served as the fundamental building block within applications that deal with problems ranging from drug design Kearnes et al. (2016) to session-based recommendation Wu et al. (2019).
Graph Neural Networks (GNNs) provide a powerful tool for machine learning on graphs, So far, the field of GNNs has been largely dominated by message passing architectures. Indeed, most of them share the same basic idea, and can be reformulated into a single common framework, socalled message passing neural networks (MPNNs) Gilmer et al. (2017). These models employ a message passing procedure to aggregate local information of vertices. For graph-related tasks, MPNNs usually apply some permutation invariant readout function to the vertex representations to produce a representation for the entire graph. The family of MPNNs has been heavily studied in the past few years, and there are now available very expressive models which have achieved state-ofthe-art results in several tasks Xu et al. (2019); Morris et al. (2019). Although the family of MPNNs is perhaps the most successful story in the field of graph representation learning, there exist models that follow different design paradigms and do not fall into this family. An example of such a model is the recently proposed Random Walk Graph Neural Network (RWNN) Nikolentzos & Vazirgiannis (2020). This model contains a number of trainable “hidden graphs”, and it compares the input graphs against these graphs using a random walk kernel which counts the number of common walks in two graphs. The emerging kernel values are fed into a fully-connected neural network which acts as the classifier or regressor. The employed random walk kernel is differentiable, and thus RWNN is end-
to-end trainable. However, this kernel considers only random walks of a small length. Such local patterns may fail to capture the overall large-scale shape of the graphs, while several interesting properties of graphs depend on the graph’s global structure. Furthermore, increasing the length of the walks has a direct impact on the model’s computational complexity.
In this paper, we propose a novel approach to tackle these challenges. Specifically, we propose a new architecture, called Geometric Random Walk Graph Neural Network (GRWNN), that generalizes the RWNN model such that it can count common walks of infinite length in two graphs. The model contains a number of trainable “hidden graphs”, and it compares the input graphs against these graphs using the geometric random walk kernel. Thus, instead of walks of small length, the proposed model considers walks of infinite length. To compute the kernel, GRWNN uses a fixedpoint iteration approach. The kernel values are then passed on to a fully-connected neural network which produces the output. The proposed neural network is end-to-end trainable since we can directly differentiate through the fixed-point equations via implicit differentation, which leads to a very efficient implementation in terms of memory requirements. Hence, we can still update the “hidden graphs” during training with backpropagation. We compare the performance of the proposed model to state-of-the-art graph kernels and recently-proposed neural architectures on several graph classification datasets. Results show that in most cases, the GRWNN model matches or outperforms competing methods. Our main contributions are summarized as follows:
• We propose a novel neural network model, Geometric Random Walk Graph Neural Network, which employs the geometric random walk kernel to produce graph representations. The model counts common walks of infinite length in the input graph and a set of randomly initialized “hidden graphs”.
• We employ an efficient scheme to compute the random walk graph kernel using fixed-point iterations. We show that we can directly differentiate through the fixed-point equations via implicit differentation, which leads to an efficient implementation.
• We evaluate the model’s performance on several standard graph classification datasets and show that it achieves results similar and in some cases superior to those obtained by recent GNNs and graph kernels.

heading: 2 RELATED WORK
text: Graph kernels have a long history in the field of graph representation learning Kriege et al. (2020). A graph kernel is a kernel function between graphs, i. e., a symmetric positive semidefinite function defined on the space of graphs. These methods generate implicitly (or explicitly) graph representations and enable the application of kernel methods such as the SVM classifier to graphs. Most graph kernels are instances of the R-convolution framework Haussler (1999), and they compare substructures extracted from the graphs to each other. Such substructures include shortest paths Borgwardt & Kriegel (2005), random walks Gärtner et al. (2003); Kashima et al. (2003), small subgraphs Shervashidze et al. (2009), and others. Our work is related to random walk kernels, i. e., kernels that compare random walks to each other. The first such kernels were proposed by Gärtner et al. Gärtner et al. (2003) and by Kashima et al. Kashima et al. (2003). The work of Kashima et al. was later refined by Mahé et al. Mahé et al. (2004). Vishwanathan et al. Vishwanathan et al. (2010) and Kang et al. Kang et al. (2012) proposed new algorithms for efficiently computing random walk kernels. These algorithms improve the time complexity of kernel computation. Sugiyama and Borgwardt studied the problem of halting (i. e., longer walks are downweighted so much that the kernel value is completely dominated by the comparison of walks of length 1) that occurs in random walk kernels, and showed that its extent depends on properties of the graphs being compared Sugiyama & Borgwardt (2015). Zhang et al. defined a different kernel which does not compare random walks to each other, but instead, compares the return probabilities of random walks Zhang et al. (2018b). Finally, Kalofolias et al. proposed a variant of the random walk kernel where structurally dissimilar vertices are not just down-weighed, but are not allowed to be visited during the simultaneous walk Kalofolias et al. (2021).
Although the first GNNs were proposed several years ago Sperduti & Starita (1997); Scarselli et al. (2009); Micheli (2009), until recently, these models had attracted limited attention. In recent years, with the rise of deep learning, a lot of models started to emerge Bruna et al. (2014); Li et al. (2015); Duvenaud et al. (2015); Atwood & Towsley (2016); Defferrard et al. (2016); Lei et al. (2017). Most
models update the representation of each vertex by aggregating the feature vectors of its neighbors. This update procedure can be viewed as a form of message passing algorithm and thus, these models are known as message passing neural networks (MPNNs) Gilmer et al. (2017). To compute a feature vector for the entire graph, MPNNs apply some permutation invariant readout function to all the vertices of the graph. The family of MPNNs has been heavily studied in the past few years, and there are now available several sophisticated models which can produce expressive graph representations Xu et al. (2019); Morris et al. (2019); Dehmamy et al. (2019); Morris et al. (2020). Despite the general recent focus on MPNNs, some works have proposed architectures that are not variants of this family of models Niepert et al. (2016); Maron et al. (2019b;a); Nikolentzos & Vazirgiannis (2020). The work closest to ours is the one reported in Nikolentzos & Vazirgiannis (2020) which presents the Random Walk Graph Neural Network (RWNN) model. In fact, in this paper, we generalize the RWNN model to compare random walks of infinite length in two graphs. Recently, another method that uses random walks to extract features which are then processed by a standard convolutional neural network was proposed Toenshoff et al. (2021). However, the proposed approach decouples data representation from learning since random walks are sampled in a preprocessing stage. Our work is also related to implicit models which have been applied successfully to many problems de Avila Belbute-Peres et al. (2018); Chen et al. (2018); Amos et al. (2018); Bai et al. (2019). The outputs of these models are determined implicitly by a solution of some underlying sub-problem. Implicit models have also been defined in the context of graph representation learning. For instance, Gu et al. proposed IGNN, a model that seeks the fixed-point of some equation which is equivalent to running an infinite number of message passing iterations Gu et al. (2020). Thus, the final representation potentially contains information from all neighbors in the graph capturing long-range dependencies. Gallicchio and Micheli proposed a similar model which generates graph representations based on the fixed point of a recursive/dynamical system, but is actually only partially trained Gallicchio & Micheli (2020). In contrast to these approaches whose objective is to apply a large (or infinite) number of message passing layers implicitly, in our setting, we employ a fixed-point iteration approach to compute the random walk kernel and then we directly differentiate through the fixed point equations via implicit differentation.

heading: 3 PRELIMINARIES

heading: 3.1 NOTATION
text: Let [n] = {1, . . . , n} ⊂ N for n ≥ 1. Let G = (V,E) be an undirected graph, where V is the vertex set and E is the edge set. We will denote by n the number of vertices and by m the number of edges. The adjacency matrix A ∈ Rn×n of a graph G is a symmetric (typically sparse) matrix used to encode edge information in the graph. The element of the ith row and jth column is equal to the weight of the edge between vertices vi and vj if such an edge exists, and 0 otherwise. The degree d(v) of a vertex v is equal to the sum of the weights of the edges that are adjacent to the vertex. For vertex-attributed graphs, every vertex in the graph is associated with a feature vector. We use X ∈ Rn×d to denote the vertex features where d is the feature dimensionality. The feature of a given vertex vi corresponds to the ith row of X.
The direct (tensor) productG× = (V×, E×) of two graphsG = (V,E) andG′ = (V ′, E′) is defined as follows:
V× = {(v, v′) ∈ V × V ′} E× = {( (v, v′), (u, u′) ) ∈ V× × V× | (v, u) ∈ E, and (v′, u′) ∈ E′ } We denote by A× the adjacency matrix ofG×, and denote by ∆× and d̄× the maximum and average of the vertex degrees of G×, respectively. Thus, d̄× = 1/n ∑ v∈V× d(v). A walk in a graph is a sequence of vertices such that consecutive vertices are linked by an edge. Performing a random walk on the direct product G× of two graphs G and G′ is equivalent to performing a simultaneous random walk on the two graphs G and G′.
We use ⊗ to represent the Kronecker product, and use to represent elementwise multiplication between two matrices or vectors of the same dimension. For a p × q matrix V, vec(V) ∈ Rpq represents the vectorized form of V, obtained by stacking its columns. Let also vec−1 denote the inverse vectorization operator which transforms a vector into a matrix, i. e., for a pq vector v, V = vec−1(v) where V ∈ Rp×q (see the appendix for the exact definition of the vec and vec−1 operators).

heading: 3.2 RANDOM WALK KERNEL
text: Given two graphs G and G′, the random walk kernel counts all pairs of matching walks on G and G′ Gärtner et al. (2003). There are different variants of the kernel. For instance, the p-step random walk kernel (where p ∈ N) counts all pairs of matching walks up to length p on two graphs. The number of matching walks can be obtained through the adjacency matrix A× of the product graph G× Vishwanathan et al. (2010) since a random walk on G× is equivalent to a simultaneous random walk on the two graphs. Assuming a uniform distribution for the starting and stopping probabilities over the vertices of two graphs, the p-step random walk kernel is defined as:
κp(G,G′) = |V×|∑ i=1 |V×|∑ i=j [ p∑ l=0 λlA l × ] ij
where λ0, λ1, λ2, . . . , λp are positive, real-valued weights, and A0× is the identity matrix, i. e., A 0 × = I. For p→∞, we obtain κ∞(G,G′) which is known as the random walk kernel. It turns out that if the sequence of weights λ0, λ1, λ2, . . . coresponds to the geometric sequence defined as λl = λl, then the limit κ∞(G,G′) can be computed analytically as follows:
k∞(G,G′) = |V×|∑ i=1 |V×|∑ i=j [ ∞∑ l=0 λlAl× ] ij = |V×|∑ i=1 |V×|∑ i=j [ (I− λA×)−1 ] ij = 1>(I− λA×)−11 (1)
It is well-known that the geometric series of matrices I+λA×+(λA×)2 + . . . converges only if the the largest-magnitude eigenvalue of A× (which is also the maximum eigenvalue if G× is a graph with non-negative edge weights), denoted by µmax× , is strictly smaller than 1/λ. Therefore, the geometric random walk kernel k∞ is well-defined only if λ < 1/µmax× . Interestingly, the maximum eigenvalue of A× is sandwiched between the average and the maximum of the vertex degrees ofG× Brouwer & Haemers (2011). We thus have that d̄× ≤ µmax× ≤ ∆×, and by setting λ < 1/∆×, the geometric series of matrices is guaranteed to converge.
By defining initial and stopping probability distributions over the vertices of G and G′, we can obtain a probabilistic variant of the geometric random walk kernel. Let p and p′ be two vectors that represent the initial probability distributions over the vertices of G and G′. Likewise, let q and q′ denote stopping probability distributions over the vertices ofG andG′. For uniform distributions for the initial and stopping probabilities over the vertices of the two graphs, we have pi = qi = 1/|V | and p′i = q ′ i = 1/|V
′|. Then, p× = pp′> and q× = qq′>, and the variant of the geometric random walk kernel can be computed as k∞(G,G′) = vec(q×)>(I− λA×)−1vec(p×).

heading: 4 GEOMETRIC RANDOM WALK GRAPH NEURAL NETWORKS
text: The proposed GRWNN model maps input graphs to vectors by comparing them against a number of “hidden graphs”, i. e., graphs whose adjacency and attribute matrices are trainable. The function that we employ to compare the input graphs against the “hidden graphs” is the geometric random walk graph kernel, one of the most well-studied kernels between graphs Gärtner et al. (2003); Mahé et al. (2004); Vishwanathan et al. (2010). The proposed GRWNN model contains N “hidden graphs” in total. The graphs may differ from each other in terms of size (i. e., number of vertices). Furthermore, the vertices and/or edges of those graphs can be annotated with continuous multi-dimensional features. As mentioned above, both the structure and the vertex attributes (if any) of these “hidden graphs” are trainable. Thus, the adjacency matrix of a “hidden graph” Gi of size n is described by a trainable matrix Wi ∈ Rn×n, while the vertex attributes are contained in the rows of another trainable matrix Qi ∈ Rn×d. Note that the “hidden graphs” correspond to weighted graphs, which can be directed or undirected graphs with or without self-loops. In our implementation, we constraint them to be undirected graphs without self-loops (n(n−1)/2 trainable parameters in total).
To compare an input graph G against a “hidden graph” Gi, the model uses the geometric random walk kernel that was introduced in the previous section:
k∞(G,Gi) = |V×|∑ i=1 |V×|∑ i=j [ ∞∑ l=0 λlAl× ] ij = |V×|∑ i=1 |V×|∑ i=j [ (I− λA×)−1 ] ij = 1>(I− λA×)−11 (2)
where A× = A⊗Ai and Ai is the adjacency matrix of “hidden graph”Gi obtained as Ai = f(W). Here, f(·) is a function whose output is non-negative and potentially bounded, i. e., f(Wi) = ReLU(Wi) or f(Wi) = σ(Wi) where σ(·) denotes the sigmoid activation function. Then, given the set Gh = {G1, G2, . . . , GN} where G1, G2, . . . , GN denote the N “hidden graphs”, we can compute N kernel values in total. These kernel values can be thought of as features of the input graph, and can be concatenated to form a vector representation of the input graph. This vector can then be fed into a fully-connected neural network to produce the output.
Following Vishwanathan et al. (2010), to compute the geometric random walk graph kernel shown in Equation equation 2 above, we employ a two-step approach . We first need to solve the following linear system for z: (I− λA×) z = 1 Then, given z, we can compute the kernel value as k∞(G,Gi) = 1>z. To solve the above linear system, we capitalize on fixed point methods. We first rewrite the above system as:
z = 1 + λA× z (3)
Now, solving for z is equivalent to finding a fixed point of Equation equation 3 Nocedal & Wright (2006). Such a fixed point can be obtained by simply iterating the first part of the forward pass. Letting z(t) denote the value of z at iteration t, we set z(0) = 1, and then compute the following:
z(t+1) = 1 + λA× z (t)
repeatedly until ||z(t+1)− z(t)|| < , where || · || denotes the Euclidean norm and some predefined tolerance or until a specific number of iterations has been reached. As mentioned in the previous section, the above problem is guaranteed to converge if the maximum eigenvalue of A× is strictly smaller than 1/λ, thus if all the eigenvalues of λA× lie inside the unit disk. If the values of the elements of Ai are bounded, we can compute an upper bound on the maximum degree of G× and set the parameter λ to some value smaller than the inverse of the upper bound.
Efficient implementation. If the input graph G consists of n vertices and a “hidden graph” Gi consists of c vertices, then A× is an nc × nc matrix. Thus, multiplying A× by some vector inside the fixed-point algorithm requires O(n2c2) operations in total. Fortunately, to compute the kernel, it is not necessary to explicitly compute matrix A×. Specifically, the Kronecker product and vec operator are linked by the well-known property Bernstein (2009):
vec(ABC) = (C> ⊗A)vec(B) (4)
Then, let Z ∈ Rn×c be a matrix such that Z = vec−1(z). Recall also that A× = A⊗Ai. Based on the above and on Equation equation 4, we can write:
A× z = (A⊗Ai)vec(Z) = vec(Ai ZA>) = vec(Ai vec−1(z)A>) (5)
The above matrix-vector product can be computed in O(n2c) time in case n > c. If A is sparse, then it can be computed yet more efficiently. Furthermore, we do not need to compute and store matrix A× which might not be feasible due to high memory requirements. Then, instead of solving the system of Equation equation 3, we solve the following equivalent system:
z = 1 + λ vec(Ai vec−1(z)A>) (6)
Node attributes. In many real-world problems, vertices of the input graphs are annotated with real-valued multi-dimensional vertex attributes. We next generalize the proposed model to graphs that contain such vertex attributes. Let X ∈ Rn×d denote the matrix that contains the vertex attributes of the input graphG. As already mentioned, we also associate a trainable matrix Qi ∈ Rc×d to each “hidden graph”Gi, where c is the number of vertices ofGi. Then, let S = σ(XQ>i ) ∈ Rc×n where σ(·) denotes the sigmoid function. The (j, k)th element of matrix S is equal to the inner product (followed by a sigmoid) between the attributes of the jth vertex of the input graph G and the kth vertex of the “hidden graph” Gi. Roughly speaking, this matrix encodes the similarity between the attributes of the vertices of the two graphs. Note that instead of directly using matrix X, we can first transform it into matrix X̃ using a single- or a multi-layer perceptron. Let s = vec(S) where s ∈ Rnc. Each element of s corresponds to a vertex of G× and quantifies the similarity between
the attributes of the pair of vertices (i. e., one from G and one from Gi) it represents. Then, we can compute the geometric random walk kernel as follows:
k∞(G,G′) = |V×|∑ i=1 |V×|∑ i=j [ ∞∑ l=0 λl ( (s s>) A× )l] ij
= |V×|∑ i=1 |V×|∑ i=j [( I− λ(s s>) A× )−1] ij = 1> ( I− λ(s s>) A× )−1 1
(7)
Note that since the elements of s take values between 0 and 1, the same applies to the elements of the output of the outer product s s>. Therefore, the maximum degree of the vertices of the graph derived from the matrix s s> A× is not greater than that of the graph derived from matrix A×, and we do not thus need to set λ to a new value. Then, to compute the kernel, we first need to solve the following system: z = 1 + λ(s s> A×)z (8) Again, naively computing the right part of the above Equation is expensive and requires O(n2c2) operations in total. The following result shows that in fact we can compute the above in a more time and memory efficient manner. Proposition 1. Let A1 ∈ Rn×n and A2 ∈ Rm×m be two real matrices. Let also s,y ∈ Rnm be two real vectors. Then, we have that:(
s s> (A1 ⊗A2) ) y = s vec ( A2 vec−1(y s)A>1 ) Based on the above result (the proof is left to the appendix), the system that needs to be solved is:
z = 1 + λ ( s vec ( Ai vec−1(z s)A> )) Since we store matrix A as a sparse matrix, if there are O(m) non-zero entries in A, then computing one iteration of the above equation for all N “hidden graphs” takes O ( Nc(n(d + c) + m)
) computational time where d is the size of the vertex attributes.
Implicit differentiation. Clearly, iteratively computing Equation equation 3 or Equation equation 8 to find the fixed point corresponds to a differentiable module. However, to train the model, we need to backpropagate the error through the fixed point solver in the backward pass. That would require storing all the intermediate terms, which could be prohibitive in practice. Fortunately, thanks to recent advances in implicit layers and equilibrium models Bai et al. (2019), this can be performed in a simple and efficient manner which requires constant memory, and assumes no knowledge of the fixed point solver. Specifically, based on ideas from Bai et al. (2019), we derive the form of implicit backpropagation specific to the employed fixed point iteration layer. Theorem 1. Let fθ be the system of Equation equation 3 or Equation equation 8, and z? ∈ Rnc be a solution to that linear system. Let also gθ(z?;A,X) = fθ(z?;A,X)− z?. Since z? is a fixed point, we have that gθ(z?;A,X) → 0 and z? is thus the root of gθ. Let y ∈ R denote the ground-truth target of the input sample, h : R → R be any differentiable function and let L : R × R → R be a loss function that computes:
` = L ( h(1>z?), y ) = L ( h ( 1>FindRoot(gθ;A,X) ) , y )
(9)
Then, the gradient of the loss w.r.t. (·) (e. g., θ, A or X) is:
∂` ∂(·) = − ∂` ∂z? ( J−1gθ ∣∣ z? )∂fθ(z?;A,X) ∂(·) = − ∂` ∂h ∂h ∂z? ( J−1gθ ∣∣ z? )∂fθ(z?;A,X) ∂(·)
(10)
where J−1gθ ∣∣ z? is the inverse Jacobian of gθ evaluated at z?.
The above formula gives a form for the necessary Jacobian without needing to backpropagate through the method used to obtain the fixed point. Thus, as mentioned above, we only need to find the fixed point, and we can compute the necessary Jacobians at this specific point using the above analytical form. No intermediate terms of the iterative method used to compute the fixed point need to be stored in memory, while there is also no need to unroll the forward computations
within an automatic differentiation layer. Still, to compute the analytical backward gradient at the solution of the fixed point equation, it is necessary to first compute the exact inverse Jacobian J−1gθ which has a cubic cost. As shown in Bai et al. (2019), we can instead compute the − ∂`∂z? ( J−1gθ ∣∣ z? ) term by solving the following linear system:
x =
( ∂fθ(z ?;A,X)
∂z?
)> x + ( ∂`
∂z? )> which in fact is also a fixed point equation and can be solved via some iterative procedure. Note that the first term of the above Equation is a vector-Jacobian product which can be efficiently computed via autograd packages (e. g., PyTorch Paszke et al. (2017)) for any x, without explicitly writing out the Jacobian matrix. Finally, we can compute ∂`∂(·) as follows:
∂`
∂(·) =
( ∂fθ(z ?;A,X)
∂(·)
)> x
where again this product is itself a vector-Jacobian product, computable via normal automatic differentiation packages.

heading: 5 EXPERIMENTAL EVALUATION
text: We next evaluate the proposed GRWNN model on standard graph classification datasets.
Datasets. We evaluate the proposed model on 10 publicly available graph classification datasets including 5 bio/chemo-informatics datasets: MUTAG, D&D, NCI1, PROTEINS, ENZYMES, and 5 social interaction datasets: IMDB-BINARY, IMDB-MULTI, REDDIT-BINARY, REDDIT-MULTI5K, COLLAB Kersting et al. (2016). To show that the proposed model also scales to larger datasets, we additionally use two Open Graph Benchmark (OGB) datasets Hu et al. (2020). Specifically, we use a molecular property prediction dataset: ogbg-molhiv, and a code summarization dataset: ogbg-code2. More details about the datasets are given in the appendix.
Experimental Setup. In the case of the 10 standard benchmark datasets, we compare the proposed model against the following three graph kernels: (1) graphlet kernel (GR) Shervashidze et al. (2009), (2) shortest path kernel (SP) Borgwardt & Kriegel (2005), and (3) Weisfeiler-Lehman subtree kernel (WL) Shervashidze et al. (2011), and against the following six neural network models: (1) DGCNN Zhang et al. (2018a), (2) DiffPool Ying et al. (2018), (3) ECC Simonovsky & Komodakis (2017), (4) GIN Xu et al. (2019), (5) GraphSAGE Hamilton et al. (2017), and (6) RWNN Nikolentzos & Vazirgiannis (2020). We also compare the proposed model against GRWNN-fixed, a variant of the model whose “hidden graphs” are randomly initialized and kept fixed during training. To evaluate the proposed model, we employ the experimental protocol proposed in Errica et al. (2020). Therefore, we perform 10-fold cross-validation to obtain an estimate of the generalization performance of each method, while within each fold a model is selected based on a 90%/10% split of the training set. We use exactly the same splits as in Errica et al. (2020) and in Nikolentzos & Vazirgiannis (2020), hence, for the different datasets, we use the results reported in these two papers.
For all datasets, we set the batch size to 64 and the number of epochs to 300. We use the Adam optimizer with initial learning rate 0.001 and applied an adaptive learning rate decay based on validation results. We use a 1-layer perceptron to transform the vertex attributes. We apply layer normalization Ba et al. (2016) on the generated graph representations (i. e., vector consisting of kernel values). The hyper-parameters we tune for each dataset are: (1) the number of “hidden graphs” ∈ {32, 64}, (2) the number of vertices of the “hidden graphs” ∈ {5, 10}, (3) the hidden-dimension size of the vertex features ∈ {32, 64} for the bio/chemo-informatics datasets and ∈ {8, 16} for the social interaction datasets, and (4) the dropout ratio ∈ {0, 0.1}. For both OGB datasets, we used the available predefined splits. We compare the proposed model against the following neural network models: GCN Kipf & Welling (2017), GIN Xu et al. (2019), GCN-FLAG Kong et al. (2020), GIN-FLAG Kong et al. (2020), PNA Corso et al. (2020), GSN Bouritsas et al. (2020), HIMP Fey et al. (2020), and DGN Beaini et al. (2020). For all models, we use the results that are reported in the respective papers. For ogbg-code2, we did not add the inverse edges to the graphs. All reported results are averaged over 10 runs.
For both OGB datasets, we set the batch size to 128. For the ogb-molhiv dataset, we set the number of epochs to 300, the number of “hidden graphs” to 200, the number of vertices of the “hidden graphs” to 5, the hidden-dimension size of the vertex features to 128 and the dropout ratio to 0.1. Furthermore, we employ the probabilistic variant of the geometric random walk kernel and use uniform distributions for the initial and stopping probabilities over the vertices of the two compared graphs. For the ogb-code2 dataset, we set the number of epochs to 100, the number of “hidden graphs” to 200, the number of vertices of the “hidden graphs” to 5, the hidden-dimension size of the vertex features to 128 and the dropout ratio to 0.2. For both datasets, we apply layer normalization Ba et al. (2016) on the generated graph representations.
Implementation Details. To set the value of parameter λ, we assume a transductive setting, where we are given a collection of graphs beforehand. Therefore, we can find the vertex of highest degree across all graphs and set the value of λ accordingly. In the inductive learning setting, since we do not know a priori target graphs that the model may receive in the future, λ should be small enough so that λ < 1/µmax× for any pair of an unseen graph and a “hidden graph”. This is a limitation of the proposed model since in case the model receives at test time a graph whose largest eigenvalue is higher than expected, we need to set λ to a smaller value and retrain the model.
To compute the fixed point of Equation equation 3 or equation 8, we followed the naive approach where we simply performed multiple times the forward iteration. In practice, there are more efficient fixed point iteration methods, such as Anderson Acceleration Walker & Ni (2011), that converge faster than the naive forward iteration at the cost of some additional memory complexity. However, as shown next, we found that in our setting, the naive forward iteration converges in a small number of steps, while the additional cost introduced by more efficient methods associated with the generation and manipulation of new tensors made them overall slower than the naive forward iteration even though they required fewer iterations to converge.
The model was implemented with PyTorch Paszke et al. (2019), and all experiments were run on a single machine equipped with an NVidia Titan Xp GPU card.
Results. Table 1 illustrates average prediction accuracies and standard deviations for the 10 standard graph classification datasets. We observe that the proposed GRWNN model is the best-performing method on 2 out of the 10 datasets, while it provides the second best and third best accuracy on 3
and 1 out of the remaining 8 datasets, respectively. The most successful method is the WL kernel which performs best on 4 of the 10 datasets, while it outperforms the other approaches with quite wide margins in most cases. Among the neural network models, the proposed GRWNN model outperforms the baseline models on 4 out of the 10 datasets. On the remaining 6 datasets, GIN is the best-performing model on half of them, and RWNN on the other half. On the ENZYMES and IMDB-BINARY datasets, our model offers respective absolute improvements of 3.1%, and 1.6% in accuracy over GIN. Overall, the model exhibits highly competitive performance on the graph classification datasets, while the achieved accuracies follow different patterns from all the baseline methods. Furthermore, the proposed model outperforms GRWNN-fixed on all datasets, demonstrating that the set of trainable “hidden graphs” is an indispensable component of the model.
The Table shown in Figure 1a illustrates the performance on the two OGB datasets. Note that the proposed model does not utilize the edge features that are provided for the different datasets. Still, we can see that it managed to outperform several of the baselines on the ogbg-molhiv dataset, where it achieved the fourth best ROC-AUC. On the ogbg-code2 dataset, GRWNN outperformed GIN, while it achieved an F1-score similar to that of GCN. However, all these three models achieved a much smaller F1-score than the one achieved by PNA which is the best-performing model.
As already discussed, the running time of the model depends on the number of fixed point iterations that need to be performed until convergence. Figure 1b (top) illustrates the average number of iterations (across all batches) for the forward and backward pass for different values of λ and for each epoch. The model was trained on a single split of the ENZYMES dataset. The maximum eigenvalue of all graphs of the dataset is equal to 5.47, while the highest degree is equal to 9. The number of nodes of the “hidden graphs” was set to 5. If the elements of the adjacency matrices of the “hidden graphs” take values no greater than one, then no vertex of G× can have a degree greater than 9∗4 = 36. Thus, setting λ < 1/36 guarantees convergence. In practice, as shown in the Figure, we found that even if λ takes larger values, we only need a small number of iterations. For λ = 1/5, we can see that the fixed point equation fails to converge since the average number of iterations is close to 100 (which is the upper limit we have set). For λ = 1/10 and for smaller values of λ, the system converges in a small number of iterations. In terms of performance, as shown in Figure 1b (bottom), the model achieves the highest levels of validation accuracy for λ = 1/20 and λ = 1/30, while for λ = 1/5, the model yields much worse performance compared to the other values of λ. Similar behavior was observed on the other datasets.

heading: 6 CONCLUSION
text: In this paper, we introduced the GRWNN model, a new architecture which generates graph representations by comparing the input graphs against “hidden graphs” using the geometric random walk kernel. To compute the kernel, the proposed model uses a fixed point iteration algorithm, and to update the “hidden graphs” during backpropagation, the model capitalizes on implicit differentation.The model was evaluated on several graph classification datasets where it achieved competitive performance.

heading: A APPENDIX
text: The appendix is organized as follows. In section B, we define some basic concepts from linear algebra. In section C, we prove the Proposition 1. In section D, we give more details about the direct differentiation through the fixed point, while section E provides a detailed description of the graph classification datasets. Section F gives details about the parameter λ, while section G about the running time of the model. Finally, in section H, we perform a sensitivity analysis.

heading: B LINEAR ALGEBRA CONCEPTS
text: In this Section, we provide definitions for concepts of linear algebra, namely the vectorization operator, the inverse vectorization operator and the Kronecker product, which we use heavily in the main paper. Definition 1. Given a real matrix A ∈ Rm×n, the vectorization operator vec : Rm×n → Rmn is defined as:
vec(A) =  A:1 A:2 ...
A:n  where A:i is the ith column of A.
Definition 2. Given a real vector b ∈ Rmn, the inverse vectorization operator vec−1 : Rnm → Rn×m is defined as:
vec−1(b) =  b1 bn+1 . . . bn(m−1)+1 b2 bn+2 . . . bn(m−1)+2 ... ... ...
... bn b2n . . . bnm  Definition 3. Given real matrices A ∈ Rn×m and B ∈ Rp×q , the Kronecker product A ⊗ B ∈ Rnp×mq defined as:
A⊗B =  A11 B A12 B . . . A1mB A21 B A22 B . . . A2mB ... ... ... ...
An1 B An2 B . . . AnmB 

heading: C PROOF OF PROPOSITION 1
text: For convenience, we restate the Proposition below. Proposition 2. Let A1 ∈ Rn×n and A2 ∈ Rm×m be two real matrices. Let also s,y ∈ Rnm be two real vectors. Then, we have that:(
s s> (A1 ⊗A2) ) y = s vec ( A2 vec−1(s y)A>1 ) Proof. Let Ds denote a diagonal matrix with the vector s as its main diagonal. Then, we have:(
s s> (A1 ⊗A2) ) y = ( Ds (A1 ⊗A2)Ds ) y
The Hadamard product of two vectors s and y is the same as matrix multiplication of vector y by the corresponding diagonal matrix Ds of vector s, i. e., Ds y = s y. Thus, it follows that:(
Ds (A1 ⊗A2)Ds ) y = Ds (A1 ⊗A2) (s y)
Note that the Kronecker product and vec operator are linked by the well-known property Bernstein (2009)(Proposition 7.1.9): vec(ABC) = (C> ⊗A)vec(B) Therefore, we have that:(
Ds (A1 ⊗A2) (s y) = Ds vec ( A2 vec−1(s y)A>1 ) = s vec(A2 vec−1(s y)A>2 )
which concludes the proof.
D IMPLICIT DIFFERENTIATION
Clearly, iteratively computing Equation (3) or Equation (8) (main paper) to find the fixed point corresponds to a differentiable module. However, to train the model, we need to backpropagate the error through the fixed point solver in the backward pass. That would require storing all the intermediate terms, which could be prohibitive in practice. Fortunately, thanks to recent advances in implicit layers and equilibrium models Bai et al. (2019), this can be performed in a simple and efficient manner which requires constant memory, and assumes no knowledge of the fixed point solver. Specifically, based on ideas from Bai et al. (2019), we derive the form of implicit backpropagation specific to the employed fixed point iteration layer. Theorem 2. Let fθ be the system of Equation (3) or Equation (8) (main paper), and z? ∈ Rnc be a solution to that linear system. Let also gθ(z?;A,X) = fθ(z?;A,X)− z?. Since z? is a fixed point, we have that gθ(z?;A,X) → 0 and z? is thus the root of gθ. Let y ∈ R denote the ground-truth target of the input sample, h : R → R be any differentiable function and let L : R × R → R be a loss function that computes:
` = L ( h(1>z?), y ) = L ( h ( 1>FindRoot(gθ;A,X) ) , y )
(11)
Then, the gradient of the loss w.r.t. (·) (e. g., θ, A or X) is: ∂`
∂(·) = − ∂` ∂z? ( J−1gθ ∣∣ z? )∂fθ(z?;A,X) ∂(·) = − ∂` ∂h ∂h ∂z? ( J−1gθ ∣∣ z? )∂fθ(z?;A,X) ∂(·)
(12)
where J−1gθ ∣∣ z? is the inverse Jacobian of gθ evaluated at z?.
The above Theorem gives a form for the necessary Jacobian without needing to backpropagate through the method used to obtain the fixed point. We can thus treat the fixed point algorithm as a black box, and we do not need to store intermediate terms associated with the fixed point algorithm into memory. We only need to apply some algorithm that will produce a solution to the system (i. e., it will compute the fixed point).
Following Bai et al. (2019), we differentiate the two sides of the fixed point equation z? = fθ(z
?;A,X) wih respect to (·): dz?
d(·) =
dfθ(z ?;A,X)
d(·) = ∂fθ(z
?;A,X)
∂(·) + ∂fθ(z
?;A,X) ∂z? dz? d(·)
=⇒ ( I− ∂fθ(z ?;A,X)
∂z?
) dz?
d(·) = ∂fθ(z
?;A,X)
∂(·) Since gθ(z?) = fθ(z?;A,X)− z?, we have:
Jgθ ∣∣ z? = − ( I− ∂fθ(z ?;A,X)
∂z? ) which implies the following:
∂`
∂(·) =
∂` ∂z? dz? d(·) = − ∂` ∂z? ( J−1gθ ∣∣ z? )∂fθ(z?;A,X) ∂(·) .
Unfortunately, computing the exact inverse Jacobian J−1gθ has a cubic cost. As shown in Bai et al. (2019), we can instead compute the − ∂`∂z? ( J−1gθ ∣∣ z? ) term of the gradient (which contains the Jacobian) by solving the following linear system:
x> = − ∂` ∂z? ( J−1gθ ∣∣ z? ) = ∂` ∂z? ( I− ∂fθ(z ?;A,X) ∂z? )−1 =⇒ x = (( I− ∂fθ(z ?;A,X)
∂z?
)>)−1( ∂`
∂z? )> =⇒ ( I− ∂fθ(z ?;A,X)
∂z?
)> x = ( ∂`
∂z? )> =⇒ x = ( ∂fθ(z ?;A,X)
∂z?
)> x + ( ∂`
∂z? )> which in fact is also a fixed point equation and can be solved via some iterative procedure. Note that the first term of the above Equation is a vector-Jacobian product which can be efficiently computed via autograd packages (e. g., PyTorch Paszke et al. (2017)) for any x, without explicitly writing out the Jacobian matrix. Finally, we can compute ∂`∂(·) as follows:
∂`
∂(·) =
( ∂fθ(z ?;A,X)
∂(·)
)> x
where again this product is itself a vector-Jacobian product, computable via normal automatic differentiation packages.

heading: E DATASETS
text: We evaluated the proposed model on 10 publicly available graph classification datasets including 5 bio/chemo-informatics datasets: MUTAG, D&D, NCI1, PROTEINS and ENZYMES, as well as 5
social interaction datasets: IMDB-BINARY, IMDB-MULTI, REDDIT-BINARY, REDDIT-MULTI5K and COLLAB Kersting et al. (2016). A summary of the 10 datasets is given in Table 2. MUTAG consists of 188 mutagenic aromatic and heteroaromatic nitro compounds. The task is to predict whether or not each chemical compound has mutagenic effect on the Gram-negative bacterium Salmonella typhimurium Debnath et al. (1991). ENZYMES contains 600 protein tertiary structures represented as graphs obtained from the BRENDA enzyme database. Each enzyme is a member of one of the Enzyme Commission top level enzyme classes (EC classes) and the task is to correctly assign the enzymes to their classes Borgwardt et al. (2005). NCI1 contains more than four thousand chemical compounds screened for activity against non-small cell lung cancer and ovarian cancer cell lines Wale et al. (2008). PROTEINS contains proteins represented as graphs where vertices are secondary structure elements and there is an edge between two vertices if they are neighbors in the amino-acid sequence or in 3D space. The task is to classify proteins into enzymes and non-enzymes Borgwardt et al. (2005). D&D contains over a thousand protein structures. Each protein is a graph whose nodes correspond to amino acids and a pair of amino acids are connected by an edge if they are less than 6 Ångstroms apart. The task is to predict if a protein is an enzyme or not Dobson & Doig (2003). IMDB-BINARY and IMDB-MULTI were created from IMDb, an online database of information related to movies and television programs. The graphs contained in the two datasets correspond to movie collaborations. The vertices of each graph represent actors/actresses and two vertices are connected by an edge if the corresponding actors/actresses appear in the same movie. Each graph is the ego-network of an actor/actress, and the task is to predict which genre an ego-network belongs to Yanardag & Vishwanathan (2015). REDDIT-BINARY and REDDIT-MULTI-5K contain graphs that model the social interactions between users of Reddit. Each graph represents an online discussion thread. Specifically, each vertex corresponds to a user, and two users are connected by an edge if one of them responded to at least one of the other’s comments. The task is to classify graphs into either communities or subreddits Yanardag & Vishwanathan (2015). COLLAB is a scientific collaboration dataset that consists of the ego-networks of several researchers from three subfields of Physics (High Energy Physics, Condensed Matter Physics and Astro Physics). The task is to determine the subfield of Physics to which the ego-network of each researcher belongs Yanardag & Vishwanathan (2015).
We also evaluated the proposed model on two datasets from the Open Graph Benchmark (OGB) Hu et al. (2020), a collection of large-scale and diverse benchmark datasets for machine learning on graphs. A summary of the two datasets is given in Table 3. The ogbg-molhiv dataset is a molecular property prediction dataset that is adopted from the MoleculeNet Wu et al. (2018). The dataset consists of 41, 127 molecules and corresponds to a binary classification dataset where the task is to predict whether a molecule inhibits HIV virus replication or not. The molecules in the training, validation and test sets are divided using a scaffold splitting procedure that splits the molecules based on their two-dimensional structural frameworks. The scaffold splitting attempts to separate structurally different molecules into different subsets. The ogbg-code2 dataset contains a large number of Abstract Syntax Trees (ASTs) that are extracted from approximately 450, 000 Python method definitions. For each method, the AST edges, the AST nodes, and the tokenized method name are retrieved. Given the body of a method represented by the AST and its node features, the task (which is known as code summarization) is to predict the sub-tokens forming the name of the method. The ASTs for the training set are obtained from GitHub projects that do not appear in the validation and test sets. We refer the reader to Hu et al. (2020) for more details about the OGB datasets.
F PARAMETER λ
Given an input graphG and a “hidden graph”Gi, since the “hidden graph” is trainable, the maximum vertex degree of the product graph G× is not known beforehand. However, in case the weights of the edges of the “hidden graph” are bounded, we can compute an upper bound to that. Let ∆ denote the maximum vertex degree of G, c denote the number of vertices of the “hidden graph” Gi, and b the maximum edge weight of the “hidden graph”. Then, we have that ∆× ≤ ∆cb, and therefore, to guarantee convergence, we need to set λ ≤ 1/∆cb. In practice, we empirically found that even if λ takes higher values, the geometric series converges within a small number of iterations. Table 4 shows the value of λ that we used for each dataset.

heading: G RUNTIME ANALYSIS
text: The proposed model is indeed computationally more expensive than the RWNN model due to the fixed point iteration which is not parallelizable. However, as already discussed, we empirically observed that the forward iteration converges in a small number of steps, thus incurring a relatively small overhead in the model’s running time. We have computed the average running time per epoch of the proposed model, and 3 of the baselines (2-RWNN, 3-RWNN and GIN) on the 10 graph classification datasets. We use the same values for the common hyperparameters (e. g., number and size of hidden graphs for GRWNN and RWNN, and hidden dimension size, batch size, etc for all 3 models). The results are shown in Table 5 (in seconds). As we can see, the proposed model is not much more expensive than the baselines. In fact, in most cases, its average running time per epoch is 1 − 3 times higher than that of the baselines, which is by no means prohibitive for real-world scenarios.

heading: H SENSITIVITY ANALYSIS
text: The proposed GRWNN model involves two main parameters: (1) the number of “hidden graphs”, and (2) the number of vertices of “hidden graphs”. We next investigate how these two parameters influence the performance of the GRWNN model. Specifically, in Figure 2, we examine how the different values of these parameters affect the performance of GRWNN on the ENZYMES dataset. We observe that the accuracy on the test set increases as the number of “hidden graphs” increases. The number of “hidden graphs” seems to have a significant impact on the performance of the model. When the number of graphs is set equal to 5, the model achieves an accuracy smaller than 50%, while when the number of graphs is set equal to 100, it yields an accuracy greater than 65%. On the other hand, the number of vertices of the “hidden graphs” does not affect that much the performance of the model.

Review summaries:
The paper discusses random walk graph kernels for graph neural network feature generation in other works. They work by generating features for graphs by comparing random walks present in each graph compared to a set of trainable hidden graphs. It is stated in this paper that an advantage of random walk kernels is that they are end to end differentiable, but current implementations have the weakness of only being efficient for walks of very short length. They introduce a way to efficiently extend walks to ‘infinite length’, theoretically making features generated by the kernel more informative. They compare their results with top GNNs and with other graph kernel approaches such as the mentioned random walk kernels on graph classification datasets. Out of 10 datasets they achieve best performance on 2 (1 of the 2 is tied for best), second best on 3, and third best on 1, It should be noted that many placements are well within the margin of error. 

The paper deals with learning with graphs, specifically graph-level classification. Given a graph G, the main idea of the paper is to compare the graph G to a set of N graphs using the (infinite-length) random walk kernel (Vishwanathan et al. 2010), resulting in a vector with N components. The resulting vector is fed into an MLP performing the classification tasks. Crucially, the N graphs are learned using end-to-end learning.  Further, the authors show how to compute the RW kernel efficiently, building on previous work  (Vishwanathan et al. 2010), and how to do differentiate through the architecture efficiently, building on (Bai et al. 2016). The paper can be seen as an extension of the RWNN model (Nikolentzos et al. 2020), which uses a finite-length random kernel.   The methodology is complemented with an experimental study showing promising performance on standard benchmarks datasets.      

In this paper, the authors proposed an interesting geometric random walk neural network for graph representation. The proposed method can be viewed as an improvement of the random walk neural network in [Nikolentzos & Vazirgiannis, 2020], which learns graph representations by modeling the geometric series of the “product” between observed graphs and some hidden and learnable graphs.  As a result, the graph representations are learned by finding fixed points of linear systems.