Year: 2020

Abstract:
In order to overcome the expressive limitations of graph neural networks (GNNs), we propose the first method that exploits vector flows over graphs to develop globally consistent directional and asymmetric aggregation functions. We show that our directional graph networks (DGNs) generalize convolutional neural networks (CNNs) when applied on a grid. Whereas recent theoretical works focus on understanding local neighbourhoods, local structures and local isomorphism with no global information flow, our novel theoretical framework allows directional convolutional kernels in any graph. First, by defining a vector field in the graph, we develop a method of applying directional derivatives and smoothing by projecting node-specific messages into the field. Then we propose the use of the Laplacian eigenvectors as such vector field, and we show that the method generalizes CNNs on an n-dimensional grid, and is provably more discriminative than standard GNNs regarding the Weisfeiler-Lehman 1-WL test. Finally, we bring the power of CNN data augmentation to graphs by providing a means of doing reflection, rotation and distortion on the underlying directional field. We evaluate our method on different standard benchmarks and see a relative error reduction of 8% on the CIFAR10 graph dataset and 11% to 32% on the molecular ZINC dataset. An important outcome of this work is that it enables to translate any physical or biological problems with intrinsic directional axes into a graph network formalism with an embedded directional field.

heading: 1 INTRODUCTION
text: One of the most important distinctions between convolutional neural networks (CNNs) and graph neural networks (GNNs) is that CNNs allow for any convolutional kernel, while most GNN methods are limited to symmetric kernels (also called isotropic kernels in the literature) (Kipf & Welling, 2016; Xu et al., 2018a; Gilmer et al., 2017). There are some implementation of asymmetric kernels using gated mechanisms (Bresson & Laurent, 2017; Veličković et al., 2017), motif attention (Peng et al., 2019), edge features (Gilmer et al., 2017) or by using the 3D structure of molecules for message passing (Klicpera et al., 2019).
However, to the best of our knowledge, there are currently no methods that allow asymmetric graph kernels that are dependent on the full graph structure or directional flows. They either depend on local structures or local features. This is in opposition to images which exhibit canonical directions: the horizontal and vertical axes. The absence of an analogous concept in graphs makes it difficult to define directional message passing and to produce an analogue of the directional frequency filters (or Gabor filters) widely present in image processing (Olah et al., 2020).
We propose a novel idea for GNNs: use vector fields in the graph to define directions for the propagation of information, with an overview of the paper presented in 1. Hence, the aggregation or message passing will be projected onto these directions so that the contribution of each neighbouring node nv will be weighted by its alignment with the vector fields at the receiving node nu. This enables our method to propagate information via directional derivatives or smoothing of the features.
We also explore using the gradients of the low-frequency eigenvectors of the Laplacian of the graph φk, since they exhibit interesting properties (Bronstein et al., 2017; Chung et al., 1997). In particular, they can be used to define optimal partitions of the nodes in a graph, to give a natural ordering (Levy, 2006), and to find the dominant directions of the graph diffusion process (Chung & Yau, 2000). Further, we show that they generalize the horizontal and vertical directional flows in a grid (see
figure 2), allowing them to guide the aggregation and mimic the asymmetric and directional kernels present in computer vision. In fact, we demonstrate mathematically that our work generalizes CNNs by reproducing all convolutional kernels of radius R in an n-dimensional grid, while also bringing the powerful data augmentation capabilities of reflection, rotation or distortion of the directions.
We further show that our directional graph network (DGN) model theoretically and empirically allows for efficient message passing across distant communities, which reduces the well-known problem of over-smoothing, and aligns well with the need of independent aggregation rules (Corso et al., 2020). Alternative methods reduce the impact of over-smoothing by using skip connections (Luan et al., 2019), global pooling (Alon & Yahav, 2020), or randomly dropping edges during training time (Rong et al., 2020), but without solving the underlying problem. In fact, we also prove that DGN is more discriminative than standard GNNs in regards to the Weisfeiler-Lehman 1-WL test, showing that the reduction of over-smoothing is accompanied by an increase of expressiveness.
Our method distinguishes itself from other spectral GNNs since the literature usually uses the low frequencies to estimate local Fourier transforms in the graph (Levie et al., 2018; Xu et al., 2019). Instead, we do not try to approximate the Fourier transform, but only to define a directional flow at each node and guide the aggregation.

heading: 2 THEORETICAL DEVELOPMENT

heading: 2.1 INTUITIVE OVERVIEW
text: One of the biggest limitations of current GNN methods compared to CNNs is the inability to do message passing in a specific direction such as the horizontal one in a grid graph. In fact, it is difficult to define directions or coordinates based solely on the shape of the graph.
The lack of directions strongly limits the discriminative abilities of GNNs to understand local structures and simple feature transformations. Most GNNs are invariant to the permutation of the neighbours’ features, so the nodes’ received signal is not influenced by swapping the features of 2 neighbours. Therefore, several layers in a deep network will be employed to understand these simple changes instead of being used for higher level features, thus over-squashing the message sent between 2 distant nodes (Alon & Yahav, 2020).
In this work, one of the main contributions is the realisation that low-frequency eigenvectors of the Laplacian can overcome this limitation by providing a variety of intuitive directional flows. As a first example, taking a grid-shaped graph of sizeN×M with N2 < M < N , we find that the eigenvector
associated to the smallest non-zero eigenvalue increases in the direction of the width N and the second one increases in the direction of the height M . This property generalizes to n-dimensional grids and motivated the use of gradients of eigenvectors as preferred directions for general graphs.
We validated this intuition by looking at the flow of the gradient of the eigenvectors for a variety of graphs, as shown in figure 2. For example, in the Minnesota map, the first 3 non-constant eigenvectors produce logical directions, namely South/North, suburb/city, and West/East.
Another important contribution also noted in figure 2 is the ability to define any kind of direction based on prior knowledge of the problem. Hence, instead of relying on eigenvectors to find directions in a map, we can simply use the cardinal directions or the rush-hour traffic flow.

heading: 2.2 VECTOR FIELDS IN A GRAPH
text: Based on a recent review from Bronstein et al. (2017), this section presents the ideas of differential geometry applied to graphs, with the goal of finding proper definitions of scalar products, gradients and directional derivatives.
Let G = (V,E) be a graph with V the set of vertices and E ⊂ V × V the set of edges. The graph is undirected meaning that (i, j) ∈ E iff (j, i) ∈ E. Define the vector spaces L2(V ) and L2(E) as the set of maps V → R and E → R with x,y ∈ L2(V ) and F ,H ∈ L2(E) and scalar products
〈x,y〉L2(V ) := ∑ i∈V xiyi , 〈F ,H〉L2(E) := ∑ (i,j)∈E F(i,j)H(i,j) (1)
Think of E as the “tangent space” to V and of L2(E) as the set of “vector fields” on the space V with each row Fi,: representing a vector at the i-th node. Define the pointwise scalar product as the map L2(E)×L2(E)→ L2(V ) taking 2 vector fields and returning their inner product at each point of V , at the node i is defined by equation 2.
〈F ,H〉i := ∑
j:(i,j)∈E
Fi,jHi,j (2)
In equation 3, we define the gradient ∇ as a mapping L2(V ) → L2(E) and the divergence div as a mapping L2(E)→ L2(V ), thus leading to an analogue of the directional derivative in equation 4.
(∇x)(i,j) := x(j)− x(i) , (divF )i := ∑
j:(i,j)∈E
F(i,j) (3)
Definition 1. The directional derivative of the function x on the graph G in the direction of the vector field F̂ where each vector is of unit-norm is
DF̂x(i) := 〈∇x, F̂ 〉i = ∑
j:(i,j)∈E
(x(j)− x(i))F̂i,j (4)
|F | will denote the absolute value of F and ||Fi,:||Lp the Lp-norm of the i-th row of F . We also define the forward/backward directions as the positive/negative parts of the field F±.

heading: 2.3 DIRECTIONAL SMOOTHING AND DERIVATIVES
text: Next, we show how the vector field F is used to guide the graph aggregation by projecting the incoming messages. Specifically, we define the weighted aggregation matrices Bav and Bdx that allow to compute the directional smoothing and directional derivative of the node features.
The directional average matrixBav is the weighted aggregation matrix such that all weights are positives and all rows have an L1-norm equal to 1, as shown in equation 5 and theorem 2.1, with a proof in the appendix C.1.
Bav(F )i,: = |Fi,:|
||Fi,:||L1 + (5)
The variable is an arbitrarily small positive number used to avoid floating-point errors. The L1norm denominator is a local row-wise normalization. The aggregator works by assigning a large weight to the elements in the forward or backward direction of the field, while assigning a small weight to the other elements, with a total weight of 1. Theorem 2.1 (Directional smoothing). The operation y = Bavx is the directional average of x, in the sense that yu is the mean of xv , weighted by the direction and amplitude of F .
The directional derivative matrix Bdx is defined in (6) and theorem 2.2, with the proof in appendix C.2. Again, the denominator is a local row-wise normalization but can be replaced by a global normalization. diag(a) is a square, diagonal matrix with diagonal entries given by a. The aggregator works by subtracting the projected forward message by the backward message (similar to a center derivative), with an additional diagonal term to balance both directions.
Bdx(F )i,: = F̂i,: − diag (∑
j
F̂:,j ) i,:
F̂i,: =
( Fi,:
||Fi,:||L1 +
) (6)
Theorem 2.2 (Directional derivative). Suppose F̂ have rows of unit L1 norm. The operation y = Bdx(F̂ )x is the centered directional derivative of x in the direction of F , in the sense of equation 4, i.e.
y = DF̂x = ( F̂ − diag (∑ j F̂:,j )) x
These aggregators are directional, interpretable and complementary, making them ideal choices for GNNs. We discuss the choice of aggregators in more details in appendix A, while also providing alternative aggregation matrices such as the center-balanced smoothing, the forward-copy, the phantom zero-padding, and the hardening of the aggregators using softmax/argmax on the field. We further provide a visual interpretation of theBav andBdx aggregators in figure 3. Interestingly, we also note in appendix A.1 thatBav andBdx yield respectively the mean and Laplacian aggregations when F is a vector field such that all entries are constant Fij = ±C. 𝑭𝑣,𝑢3 𝑭𝑣,𝑢1 𝑭𝑣,𝑢2 Direc�onal smoothing aggrega�on 𝑩𝑎𝑣 𝑭 𝒙 Direc�onal deriva�ve aggrega�on 𝑩𝑑𝑥 𝑭 𝒙Graph features focused on the neighbourhood of 𝒏𝒗
𝑣: Node receiving the message𝑢1,2,3: Neighbouring node𝒙𝑢 : Feature at node 𝑢𝑭𝑣,𝑢 : Direc�onal vector field between the node 𝑣 and 𝑢 Weighted forward deriva�ve with 𝑢1 Weighted backwardderiva�ve with 𝑢2 Weighted backwardderiva�ve with 𝑢3+ +Sum of the absolute weights
Figure 3: Illustration of how the directional aggregation works at a node nv , with the arrows representing the direction and intensity of the field F .

heading: 2.4 GRADIENT OF THE EIGENVECTORS AS INTERPRETABLE VECTOR FIELDS
text: In this section we give theoretical support for the choice of gradients of the eigenfunctions of the Laplacian as sensible vectors along which to do directional message passing since they are interpretable and allow to reduce the over-smoothing.
As usual the combinatorial, degree-normalized and symmetric normalized Laplacian are defined as
L = D −A , Lnorm = D−1L , Lsym = D− 1 2LD− 1 2 (7)
The problems of over-smoothing and over-squashing are critical issues in GNNs (Alon & Yahav, 2020; Hamilton, 2020). In most GNN models, node representations become over-smoothed after several rounds of message passing (i.e., convolutions), as the representations tend to reach a meanfield equilibrium equivalent to the stationary distribution of a random walk (Hamilton, 2020). Oversmoothing is also related to the problem of over-squashing, which reflects the inability for GNNs to propagate informative signals between distant nodes (Alon & Yahav, 2020) and is a major bottleneck to training deep GNN models (Xu et al., 2019). Both problems are related to the fact that the influence of one node’s input on the final representation of another node in a GNN is determined by the likelihood of the two nodes co-occurring on a truncated random walk (Xu et al., 2018b).
We show in theorem 2.3 (proved in appendix C.3) that by passing information in the direction of φ1, the eigenvector associated to the lowest non-trivial frequency of Lnorm, DGNs can efficiently share information between the farthest nodes of the graph, when using the K-walk distance to measure the difficulty of passing information. Thus, DGNs provide a natural way to address both the oversmoothing and over-squashing problems: they can efficiently propagate messages between distant nodes and in a direction that counteracts over-smoothing. Definition 2 (K-walk distance). The K-walk distance dK(vi, vj) on a graph is the average number of times vi is hit in a K step random walk starting from vj . Theorem 2.3 (K-Gradient of the low-frequency eigenvectors). Let λi and φi be the eigenvalues and eigenvectors of the normalized Laplacian of a connected graph Lnorm and let a, b = arg max1≤i,j≤n{dK(vi, vj)} be the nodes that have highest K-walk distance. Let m = arg min1≤i≤n(φ1)i and M = arg max1≤i≤n(φ1)i, then dK(vm, vM ) − dK(va, vb) has order O(1− λ2).
As another point of view into the problem of oversmoothing, consider the hitting time Q(x, y) defined as the expected number of steps in a random walk starting from node x ending in node y with the probability transition P (x, y) = 1dx . In appendix C.4 we give an informal argument supporting the following conjecture. Definition 3 (Gradient step). Suppose the two neighboring nodes x and z are such that φ(z)−φ(x) is maximal among the neighbors of x, then we will say z is obtained from x by taking a step in the direction of the gradient∇φ. Conjecture 2.4 (Gradient steps reduce expected hitting time). Suppose that x, y are uniformly distributed random nodes such that φi(x) < φi(y). Let z be the node obtained from x by taking one step in the direction of∇φi, then the expected hitting time is decreased proportionally to λ−1i and
Ex,y[Q(z, y)] ≤ Ex,y[Q(x, y)]
The next two corollaries follow from theorem 2.3 (and also conjecture 2.4 if it is true). Corollary 2.5 (Reduces over-squashing). Following the direction of ∇φ1 is an efficient way of passing information between the farthest nodes of the graph (in terms of the K-walk distance). Corollary 2.6 (Reduces over-smoothing). Following the direction of∇φ1 allows the influence distribution between node representations to be decorrelated from random-walk hitting times (assuming the definition of influence introduced in Xu et al. (2018b)).
Our method also aligns perfectly with a recent proof that multiple independent aggregators are needed to distinguish neighbourhoods of nodes with continuous features (Corso et al., 2020).
When using eigenvectors of the Laplacian φi to define directions in a graph, we need to keep in mind that there is never a single eigenvector associated to an eigenvalue, but a whole eigenspace.
For instance, a pair of eigenvalues can have a multiplicity of 2 meaning that they can be generated by different pairs of orthogonal eigenvectors. For an eigenvalue of multiplicity 1, there are always two unit norm eigenvectors of opposite sign, which poses a problem during the directional aggregation. We can make a choice of sign and later take the absolute value (i.e. Bav in equation 5). An alternative is to take a sample of orthonormal basis of the eigenspace and use each choice to augment the training (see section 2.8). Although multiplicities higher than one do happen for low-frequencies (square grids have a multiplicity 2 for λ1) this is not common in “real-world graphs”; we found no λ1 multiplicity greater than 1 on the ZINC and PATTERN datasets (see appendix B.4). Further, although all φ are orthogonal, their gradients, used to define directions, are not always locally orthogonal (e.g. there are many horizontal flows in the grid). This last concern is left to be addressed in future work.

heading: 2.5 GENERALIZATION OF THE CONVOLUTION ON A GRID
text: In this section we show that our method generalizes CNNs by allowing to define any radius-R convolutional kernels in grid-shaped graphs. The radius-R kernel at node u is a convolutional kernel that takes the weighted sum of all nodes v at a distance d(u, v) ≤ R. Consider the lattice graph Γ of size N1 × N2 × ... × Nn where each vertices are connected to their direct non-diagonal neighbour. We know from Lemma C.1 that, for each dimension, there is an eigenvector that is only a function of this specific dimension. For example, the lowest frequency eigenvectorφ1 always flows in the direction of the longest length. Hence, the Laplacian eigenvectors of the grid can play a role analogous to the axes in Euclidean space, as shown in figure 2.
With this knowledge, we show in theorem 2.7 (proven in C.7), that we can generalize all convolutional kernels in an n-dimensional grid. This is a strong result since it demonstrates that our DGN framework generalizes CNNs when applied on a grid, thus closing the gap between GNNs and the highly successful CNNs on image tasks. Theorem 2.7 (Generalization radius-R convolutional kernel in a lattice). For an n-dimensional lattice, any convolutional kernel of radius R can be realized by a linear combination of directional aggregation matrices and their compositions.
As an example, figure 4 shows how a linear combination of the first and m-th aggregators B(∇φ1,m) realize a kernel on an N ×M grid, where m = dN/Me and N > M .

heading: 2.6 EXTENDING THE RADIUS OF THE AGGREGATION KERNEL
text: Having aggregation kernels for neighbours of distance 2 or 3 is important to improve the expressiveness of GNNs, their ability to understand patterns, and to reduce the number of layers required. However, the lack of directions in GNNs strongly limits the radius of the kernels since, given a graph of regular degree d, a mean/sum aggregation at a radius-R will result in a heavy over-squashing of dR messages. Using the directional fields, we can enumerate different paths, thus assigning a different weight for different R-distant neighbours. This method, proposed in appendix A.7, avoids the over-squashing, but empirical results are left for future work.

heading: 2.7 COMPARISON WITH WEISFEILER-LEHMAN (WL) TEST
text: We also compare the expressiveness of the Directional Graph Networks with the classical WL graph isomorphism test which is often used to classify the expressivity of graph neural networks (Xu et al., 2018a). In theorem 2.8 (proven in appendix C.8) we prove that DGNs are capable of distinguishing pairs of graphs that the 1-WL test (and so ordinary GNNs) cannot differentiate.
Theorem 2.8 (Comparison with 1-WL test). DGNs using the mean aggregator, any directional aggregator of the first eigenvector and injective degree-scalers are strictly more powerful than the 1-WL test.

heading: 2.8 DATA AUGMENTATION
text: Another important result is that the directions in the graph allow to replicate some of the most common data augmentation techniques used in computer vision, namely reflection, rotation and distortion. The main difference is that, instead of modifying the image (such as a 5◦ rotation), the proposed transformation is applied on the vector field defining the aggregation kernel (thus rotating the kernel by−5◦ without changing the image). This offers the advantage of avoiding to pre-process the data since the augmentation is done directly on the kernel at each iteration of the training.
The simplest augmentation is the vector field flipping, which is done changing the sign of the field F , as stated in definition 4. This changes the sign ofBdx, but leavesBav unchanged. Definition 4 (Reflection of the vector field). For a vector field F , the reflected field is −F .
Let F1,F2 be vector fields in a graph, with F̂1 and F̂2 being the field normalized such that each row has a unitary L2-norm. Define the angle vector α by 〈(F̂1)i,:, (F̂2)i,:〉 = cos(αi). The vector field F̂⊥2 is the normalized component of F̂2 perpendicular to F̂1. The equation below defines F̂ ⊥ 2 . The next equation defines the angle
(F̂⊥2 )i,: = (F̂2 − 〈F̂1, F̂2〉F̂1)i,: ||(F̂2 − 〈F̂1, F̂2〉F̂1)i,:||
Notice that we then have the decomposition (F̂2)i,: = cos(αi)(F̂1)i,: + sin(αi)(F̂⊥2 )i,:.
Definition 5 (Rotation of the vector fields). For F̂1 and F̂2 non-colinear vector fields with each vector of unitary length, their rotation by the angle θ in the plane formed by {F̂1, F̂2} is
F̂ θ1 = F̂1diag(cos θ)+ F̂ ⊥ 2 diag(sin θ) , F̂ θ 2 = F̂1diag(cos(θ+α))+ F̂ ⊥ 2 diag(sin(θ+α)) (8)
Finally, the following augmentation has a similar effect to a wave distortion applied on images. Definition 6 (Random distortion of the vector field). For vector field F and anti-symmetric random noise matrixR, its randomly distorted field is F ′ = F +R ◦A.

heading: 3 IMPLEMENTATION
text: We implemented the models using the DGL and PyTorch libraries and we provide the code at the address https://anonymous.4open.science/r/a752e2b1-22e3-40ce-851c-a564073e1fca/. We test our method on standard benchmarks from Dwivedi et al. (2020) and Hu et al. (2020), namely ZINC, CIFAR10, PATTERN and MolHIV with more details on the datasets and how we enforce a fair comparison in appendix B.1.
For the empirical experiments we inserted our proposed aggregation method in two different type of message passing architecture used in the literature: a simple one similar to the one present in GCN (equation 9a) (Kipf & Welling, 2016) and a more complex and general one typical of MPNN (9b) (Gilmer et al., 2017) with or without edge features eji. Hence, the time complexity O(Em) is identical to the PNA (Corso et al., 2020), where E is the number of edges and m the number of aggregators, with an additional O(Ek) to pre-compute the k-first eigenvectors, as explained in the appendix B.2.
X (t+1) i = U ( ⊕ (j,i)∈E X (t) j ) (9a) X(t+1)i = U ( X (t) i , ⊕ (j,i)∈E M ( X (t) i , X (t) j , eji︸︷︷︸
optional
)) (9b)
where ⊕
is an operator which concatenates the results of multiple aggregators, X is the node features, M is a linear transformation and U a multiple layer perceptron.
We tested the directional aggregators across the datasets using the gradient of the first k eigenvectors ∇φ1,...,k as the underlying vector fields. Here, k is a hyperparameter, usually 1 or 2, but could be bigger for high-dimensional graphs. To deal with the arbitrary sign of the eigenvectors, we take the absolute value of the result of equation 6, making it invariant to a reflection of the field. In case of a disconnected graph, φi is the i-th eigenvector of each connected component. Despite the numerous aggregators proposed in appendix A, onlyBdx andBav are tested empirically.

heading: 4 RESULTS AND DISCUSSION
text: Directional aggregation Using the benchmarks introduced in section 3, we present in figure 5 a fair comparison of various aggregation strategies using the same parameter budget and hyperparameters. We see a consistent boost in the performance for simple, complex and complex with edges models using directional aggregators compared to the mean-aggregator baseline.
In particular, we see a significant improvement in ZINC and MolHIV using the directional aggregators. We believe this is due to the capacity to move efficiently messages across opposite parts of the molecule and to better understand the role of atom pairs. Further, the thesis that DGNs can bridge the gap between CNNs and GNNs is supported by the clear improvements on CIFAR10 over the baselines. This contrasts with the positional encoding which showed no clear improvement.
With our theoretical analysis in mind, we expected to perform well on PATTERN since the flow of the first eigenvectors are meaningful directions in a stochastic block model and passing messages using those directions allows the network to efficiently detect the two communities. The results match our expectations, outperforming all the previous models.
Comparison to the literature In order to compare our model with the literature, we fine-tuned it on the various datasets and we report its performance in figure 6. We observe that DGN provides significant improvement across all benchmarks, highlighting the importance of anisotropic kernels. In the work by Dwivedi et al. (2020), they proposed the use of positional encoding of the eigenvectors in node features, but these bring significant improvement when many eigenvectors and high network depths are used. Our results outperform them with fewer parameters, less depth, and only 1-2 eigenvectors, further motivating their use as directional flows instead of positional encoding.
Data augmentation To evaluate the effectiveness of the proposed augmentation, we trained the models on a reduced version of the CIFAR10 dataset. The results in figure 7 show clearly a higher expressive power of the dx aggregator, enabling it to fit well the training data. For a small dataset, this comes at the cost of overfitting and a reduced test-set performance, but we observe that randomly rotating or distorting the kernels counteracts the overfitting and improves the generalization.
As expected, the performance decreases when the rotation or distortion is too high since the augmented graph changes too much. In computer vision images similar to CIFAR10 are usually rotated by less than 30◦ (Shorten & Khoshgoftaar; O’Gara & McGuinness, 2019). Further, due to the constant number of parameters across models, less parameters are attributed to the mean aggregation in
the directional models, thus it cannot fit well the data when the rotation/distortion is too strong since the directions are less informative. We expect large models to perform better at high angles.

heading: 5 CONCLUSION
text: The proposed DGN method allows to solve many problems of GNNs, including the lack of anisotropy, the low expressiveness, the over-smoothing and over-squashing. For the first time in graph networks, we generalize the directional properties of CNNs and their data augmentation capabilities. Based on an intuitive idea and backed by a set of strong theoretical and empirical results, we believe this work will give rise to a new family of directional GNNs. Future work can focus on the implementation of radius-R kernels and improving the choice of multiple orthogonal directions.
Broader Impact This work will extend the usability of graph networks to all problems with physically defined directions, thus making GNN a new laboratory for physics, material science and biology. In fact, the anisotropy present in a wide variety of systems could be expressed as vector fields (spinor, tensor) compatible with the DGN framework, without the need of eigenvectors. One example is magnetic anisotropicity in metals, alloys and also in molecules such as benzene ring, alkene, carbonyl, alkyne that are easier or harder to magnetise depending on the directions or which way the object is rotated. Other examples are the response of materials to high electromagnetic fields (e.g. to study material responses at terahertz frequency); all kind of field propagation in crystals lattices (vibrations, heat, shear and frictional force, young modulus, light refraction, birefringence); multi-body or liquid motion; traffic modelling; and design of novel materials and constrained structures. This also enables GNNs to be used for virtual prototyping systems since the added directional constraints could improve the analysis of a product’s functionality, manufacturing and behavior.

heading: AUTHOR CONTRIBUTIONS
text: Anonymous

heading: ACKNOWLEDGMENTS
text: Anonymous

heading: A APPENDIX - CHOICES OF DIRECTIONAL AGGREGATORS
text: This appendix helps understand the choice of Bav and Bdx in section 2.3 and presents different directional aggregators that can be used as an alternative to the ones proposed.
A simple alternative to the directional smoothing and directional derivative operator is to simply take the forward/backward values according to the underlying positive/negative parts of the field F , since it can effectively replicate them. However, there are many advantage of using Bav,dx. First, one can decide to use either of them and still have an interpretable aggregation with half the parameters. Then, we also notice that Bav,dx regularize the parameter by forcing the network to take both forward and backward neighbours into account at each time, and avoids one of the neighbours becoming too important. Lastly, they are robust to a change of sign of the eigenvectors since Bav is sign invariant and Bdx will only change the sign of the results, which is not the case for forward/backward aggregations.

heading: A.1 RETRIEVING THE MEAN AND LAPLACIAN AGGREGATIONS
text: It is interesting to note that we can recover simple aggregators from the aggregation matrices Bav(F ) and Bdx(F ). Let F be a vector field such that all edges are equally weighted Fij = ±C for all edges (i, j). Then, the aggregatorBav is equivalent to a mean aggregation:
Bav(F )x = D −1Ax
Under the condition Fij = C, the differential aggregator is equivalent to a Laplacian operator L normalized using the degreeD
Bdx(CA)x = D −1(A−D)x = −D−1Lx

heading: A.2 GLOBAL FIELD NORMALIZATION
text: The proposed aggregators are defined with a row-wise normalized field
F̂i,: = Fi,:
||Fi,:||LP
meaning that all the vectors are of unit-norm and the aggregation/message passing is done only according to the direction of the vectors, not their amplitude. However, it is also possible to do a global normalization of the field F by taking a matrix-norm instead of a vector-norm. Doing so will modulate the aggregation by the amplitude of the field at each node. One needs to be careful since a global normalization might be very sensitive to the number of nodes in the graph.

heading: A.3 CENTER-BALANCED AGGREGATORS
text: A problem arises in the aggregators Bdx and Bav proposed in equations 5 and 6 when there is an imbalance between the positive and negative terms of F±. In that case, one of the directions overtakes the other in terms of associated weights.
An alternative is also to normalize the forward and backward directions separately, to avoid having either the backward or forward direction dominating the message.
Bav−center(F )i,: = F ′+i,: + F ′− i,:
||F ′+i,j + F ′− i,j ||L1
, F ′±i,: = |F±i,: |
||F±i,: ||L1 + (10)
The same idea can be applied to the derivative aggregator equation 11 where the positive and negative parts of the field F± are normalized separately to allow to project both the forward and backward messages into a vector field of unit-norm. F+ is the out-going field at each node and is used for the forward direction, while F− is the in-going field used for the backward direction. By averaging the forward and backward derivatives, the proposed matrix Bdx-center represents the centered derivative matrix.
Bdx-center(F )i,: = F ′ i,: − diag ∑ j F ′:,j  i,: , F ′i,: = 1 2  F + i,:
||F+i,:||L1 + ︸ ︷︷ ︸ forward field
+ F−i,:
||F−i,: ||L1 + ︸ ︷︷ ︸ backward field
 (11)

heading: A.4 HARDENING THE AGGREGATORS
text: The aggregation matrices that we proposed, mainly Bdx and Bav depend on a smooth vector field F . At any given node, the aggregation will take a weighted sum of the neighbours in relation to the direction of F . Hence, if the field Fv at a node v is diagonal in the sense that it gives a non-zero weight to many neighbours, then the aggregator will compute a weighted average of the neighbours.
Although there are clearly good reasons to have this weighted-average behaviour, it is not necessarily desired in every problem. For example, if we want to move a single node across the graph, this behaviour will smooth the node at every step. Instead, we propose below to soften and harden the aggregations by forcing the field into making a decision on the direction it takes.
Soft hardening the aggregation is possible by using a softmax with a temperature T on each row to obtain the field Fsofthard.
(Fsofthard)i,: = sign(Fi,:)softmax(T |Fi,:|) (12)
Hardening the aggregation is possible by using an infinite temperature, which changes the softmax functions into argmax. In this specific case, the node with the highest component of the field will be copied, while all other nodes will be ignored.
(Fhard)i,: = sign(Fi,:)argmax(|Fi,:|) (13)
An alternative to the aggregators above is to take the softmin/argmin of the negative part and the softmax/argmax of the positive part.

heading: A.5 FORWARD AND BACKWARD COPY
text: The aggregation matrices Bav and Bdx have the nice property that if the field is flipped (change of sign), the aggregation gives the same result, except for the sign of Bdx. However, there are cases where we want to propagate information in the forward direction of the field, without smoothing it with the backward direction. In this case, we can define the strictly forward and strictly backward fields below, and use them directly with the aggregation matrices.
Fforward = F + , Fbackward = F − (14)
Further, we can use the hardened fields in order to define a forward copy and backward copy, which will simply copy the node in the direction of the highest field component.
Fforward copy = F + hard , Fbackward copy = F − hard (15)

heading: A.6 PHANTOM ZERO-PADDING
text: Some recent work in computer vision has shown the importance of zero-padding to improve CNNs by allowing the network to understand it’s position relative to the border (Islam et al., 2020). In contrast, using boundary conditions or reflection padding makes the network completely blind to positional information. In this section, we show that we can mimic the zero-padding in the direction of the field F for both aggregation matricesBav andBdx.
Starting with theBav matrix, in the case of a missing neighbour in the forward/backward direction, the matrix will compensate by adding more weights to the other direction, due to the denominator which performs a normalization. Instead, we would need the matrix to consider both directions separately so that a missing direction would result in zero padding. Hence, we define Bav,0pad below, where either the F+ or F− will be 0 on a boundary with strictly in-going/out-going field.
(Bav,0pad)i,: = 1
2
( |F+i,:|
||F+i,:||L1 + + |F−i,: | ||F−i,: ||L1 +
) (16)
Following the same argument, we define Bdx,0pad below, where either the forward or backward term is ignored. The diagonal term is also removed at the boundary so that the result is a center derivative equal to the subtraction of the forward term with the 0-term on the back (or vice-versa), instead of a forward derivative.
Bdx−0pad(F )i,: =  F ′+i,: if ∑ j F ′− i,j = 0 F ′−i,: if ∑ j F ′+ i,j = 0
1 2 ( F ′+i,: + F ′− i,: − diag (∑ j F ′+ :,j + F ′− :,j ) i,: ) , otherwise
F ′+i,: = F+i,:
||F+i,:||L1 + F ′−i,: =
F−i,:
||F−i,: ||L1 +
(17)

heading: A.7 EXTENDING THE RADIUS OF THE AGGREGATION KERNEL
text: We aim at providing a general radius-R kernelBR that assigns different weights to different subsets of nodes nu at a distance R from the center node nv .
First, we decompose the matrix B(F ) into positive and negative parts B±(F ) representing the forward and backward steps aggregation in the field F .
B(F ) = B+(F )−B−(F ) (18)
Thus, defining B±fb(F )i,: = F±i,:
||Fi,:||Lp , we can find different aggregation matrices by using different combinations of walks of radius R. First demonstrated for a grid in theorem 2.7, we generalize it in equation 19 for any graph G. Definition 7 (General radius R n-directional kernel). Let Sn be the group of permutations over n elements with a set of directional fields Fi.
BR := ∑
V={v1,v2,...,vn}∈Nn ||V ||L1≤R, −R≤vi≤R︸ ︷︷ ︸
Any choice of walk V with at mostR steps using all combinations of v1, v2, ..., vn
∑ σ∈Sn︸︷︷︸ optional
permutations
aV N∏ j=1 (B sgn(vσ(j)) fb (Fσ(j))) |vσ(j)|
︸ ︷︷ ︸ Aggregator following the steps V , permuted by Sn
(19)
In this equation, n is the number of directional fields and R is the desired radius. V represents all the choices of walk {v1, v2, ..., vn} in the direction of the fields {F1,F2, ...,Fn}. For example, V = {3, 1, 0,−2} has a radius R = 6, with 3 steps forward of F1, 1 step forward of F2, and 2 steps backward of F4. The sign of eachB±fb is dependant to the sign of vσ(j), and the power |vσ(j)| is the
number of aggregation steps in the directional field Fσ(j). The full equation is thus the combination of all possible choices of paths across the set of fields Fi, with all possible permutations. Note that we are restricting the sum to vi having only a possible sign; although matrices don’t commute, we avoid choosing different signs since it will likely self-intersect a lower radius walk. The permutations σ are required since, for example, the path up→ left is different (in a general graph) than the path left→ up.
This matrix BR has a total of ∑R r=0(2n) r = (2n) R+1−1
2n−1 parameters, with a high redundancy since some permutations might be very similar, e.g. for a grid graph we have that up → left is identical to left → up. Hence, we can replace the permutation Sn by a reverse ordering, meaning that ∏N j Bj = BN ...B2B1. Doing so does not perfectly generalize the radius-R kernel for all graphs, but it generalizes it on a grid and significantly reduces the number of parameters to∑R r=0 ∑min(n,r) l=1 2 r ( n l )( r−1 l−1 ) .

heading: B APPENDIX - IMPLEMENTATION DETAILS

heading: B.1 BENCHMARKS AND DATASETS
text: We use a variety of benchmarks proposed by Dwivedi et al. (2020) and Hu et al. (2020) to test the empirical performance of our proposed methods. In particular, to have a wide variety of graphs and tasks we chose:
1. ZINC, a graph regression dataset from molecular chemistry. The task is to predict a score that is a subtraction of computed properties logP − SA, with logP being the computed octanol-water partition coefficient, and SA being the synthetic accessibility score (Jin et al., 2018).
2. CIFAR10, a graph classification dataset from computer vision (Krizhevsky, 2009). The task is to classify the images into 10 different classes, with a total of 5000 training image per class and 1000 test image per class. Each image has 32× 32 pixels, but the pixels have been clustered into a graph of ∼ 100 super-pixels. Each super-pixel becomes a node in an almost grid-shaped graph, with 8 edges per node. The clustering uses the code from Knyazev et al. (2019), and results in a different number of super-pixels per graph.
3. PATTERN, a node classification synthetic benchmark generated with Stochastic Block Models, which are widely used to model communities in social networks. The task is to classify the nodes into 2 communities and it tests the fundamental ability of recognizing specific predetermined subgraphs.
4. MolHIV, a graph classification benchmark from molecular chemistry. 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.
Our goal is to provide a fair comparison to demonstrate the capacity of our proposed aggregators. Therefore, we compare the various methods on both types of architectures using the same hyperparameters tuned in previous works (Corso et al., 2020) for similar networks. The models vary exclusively in the aggregation method and the width of the architectures to keep a set parameter budget.
In CIFAR10 it is impossible to numerically compute a deterministic vector field with eigenvectors due to the multiplicity of λ1 being greater than 1. This is caused by the symmetry of the square image, and is extremely rare in real-world graphs. Therefore, we used as underlying vector field the gradient of the coordinates of the image. Note that these directions are provided in the nodes’ features in the dataset and available to all models, that they are co-linear to the eigenvectors of the grid as per lemma C.1, and that they mimic the inductive bias in CNNs.
B.2 IMPLEMENTATION AND COMPUTATIONAL COMPLEXITY
Unlike several more expressive graph networks (Kondor et al., 2018; Maron et al., 2018), our method does not require a computational complexity superlinear with the size of the graph. The calculation
of the first k eigenvectors during pretraining, done using Lanczos method (Lanczos, 1950) and the sparse module of Scipy, has a time complexity of O(Ek) where E is the number of edges. During training the complexity is equivalent to a m-aggregator GNN O(Em) (Corso et al., 2020) for the aggregation and O(Nm) for the MLP.
To all the architectures we added residual connections (He et al., 2016), batch normalization (Ioffe & Szegedy, 2015) and graph size normalization (Dwivedi et al., 2020).
For all the datasets with non-regular graphs, we combine the various aggregators with logarithmic degree-scalers as in Corso et al. (2020).
An important thing to note is that, for dynamic graphs, the eigenvectors need to be re-computed dynamically with the changing edges. Fortunately, there are random walk based algorithms that can estimate φ1 quickly, especially for small changes to the graph (Doshi & Eun, 2000). In the current empirical results, we do not work with dynamic graphs.

heading: B.3 RUNNING TIME
text: The precomputation of the first four eigenvectors for all the graphs in the datasets takes 38s for ZINC, 96s for PATTERN and 120s for MolHIV on CPU. Table 1 shows the average running time on GPU for all the various model from figure 5. On average, the epoch running time is 16% slower for the DGN compared to the mean aggregation, but a faster convergence for DGN means that the total training time is on average 8% faster for DGN.

heading: B.4 EIGENVECTOR MULTIPLICITY
text: The possibility to define equivariant directions using the low-frequency Laplacian eigenvectors is subject to the uniqueness of those vectors. When the dimension of the eigenspaces associated with the lowest eigenvalues is 1, the eigenvectors are defined up to a constant factor. In section 2.4, we propose the use of unit vector normalization and an absolute value to eliminate the scale and sign ambiguity. When the dimension of those eigenspaces is greater than 1, it is not possible to define equivariant directions using the eigenvectors.
Fortunately, it is very rare for the Laplacian matrix to have repeated eigenvalues in real-world datasets. We validate this claim by looking at ZINC and PATTERN datasets where we found no graphs with repeated Fiedler vector and only one graph out of 26k with multiplicity of the second eigenvector greater than 1.
When facing a graph that presents repeated Laplacian eigenvalues, we propose to randomly shuffle, during training time, different eigenvectors randomly sampled in the eigenspace. This technique will act as a data augmentation of the graph during training time allowing the network to train with multiple directions at the same time.

heading: C APPENDIX - MATHEMATICAL PROOFS

heading: C.1 PROOF FOR THEOREM 2.1 (DIRECTIONAL SMOOTHING)
text: The operation y = Bavx is the directional average of x, in the sense that yu is the mean of xv , weighted by the direction and amplitude of F .
Proof. This should be a simple proof, that if we want a weighted average of our neighbours, we simply need to multiply the weights by each neighbour, and divide by the sum of the weights. Of course, the weights should be positive.

heading: C.2 PROOF FOR THEOREM 2.2 (DIRECTIONAL DERIVATIVE)
text: Suppose F̂ have rows of unit L1 norm. The operation y = Bdx(F̂ )x is the centered directional derivative of x in the direction of F , in the sense of equation 4, i.e.
y = DF̂x = ( F̂ − diag (∑ j F̂:,j )) x
Proof. Since F rows have unit L1 norm, F̂ = F . The i-th coordinate of the vector( F − diag (∑ j F:,j )) x is
Fx− diag ∑
j
F x  i = ∑ j Fi,jx(j)− ∑ j Fi,j x(i) =
∑ j:(i,j)∈E (x(j)− x(i))Fi,j
= DF x(i)

heading: C.3 PROOF FOR THEOREM 2.3 (K-GRADIENT OF THE LOW-FREQUENCY EIGENVECTORS)
text: Let λi and φi be the eigenvalues and eigenvectors of the normalized Laplacian of a connected graph Lnorm and let a, b = arg max1≤i,j≤n{dK(vi, vj)} be the nodes that have highest K-walk distance. Let m = arg min1≤i≤n(φ1)i and M = arg max1≤i≤n(φ1)i, then dK(vm, vM ) − dK(va, vb) has order O(1− λ2).
Proof. For this theorem, we use the indices i = 0, ..., (N − 1), sorted such that λi ≤ λi+1. Hence, λ0 = 0 and λ1 is the first non-trivial eigenvalue.
First we need the following proposition:
Proposition 1 (K-walk distance matrix). The K-walk distance matrix P associated with a graph is the matrix such that (P )i,j = dK(vi, vj) can be written as ∑K p=1W
p, where W = D−1A is the random walk matrix.
Let’s defineW = D−1A the random walk matrix of the graph.
First, we are going to show that W is jointly diagonalizable with Lnorm and we are going to relate its eigenvectors φ′i and its eigenvalues λ ′ i with the ones ofW .
Indeed,Lsym is a symmetric real matrix which is semi-positive definite diagonalizable by the spectral theorem. Since the matrix Lnorm is similar toD 1 2LnormD − 12 =D− 1 2LD− 1 2 = Lsym and the matrix of similarity isD 1 2 , a positive definite matrix, Lnorm is diagonalizable and semi-positive definite.
By
Lnorm = D −1L = D−1(L+D −D) = I +D−1(L−D) = I −D−1A = I −W
the random walk matrix is jointly diagonalizable with the random walk Laplacian. Also their eigenvalues and eigenvectors are related to each other by φi = φ′n−1−i and λ ′ i = 1− λn−1−i
Moreover, the constant eigenvector associated with eigenvalue 0 of the Random walk Laplacian, is the eigenvector associated with the highest eigenvalue of the Random walk matrix and by the formula obtained, λ′n−1 = 1− λ0 = 1 Now, we are going to approximate the K-walk distance matrix P using the 2 eigenvectors of the Random walk matrix associated with the highest eigenvalues.
By Proposition 1 we have that P = ∑K p=1W p, which can be written as
K∑ p=1 ( n−1∑ i=0 φ′iφ ′T i (λ ′ i)) p = K∑ p=1 n−1∑ i=0 φ′iφ ′T i (λ ′ i) p
by eigen-decomposition.
Since λn−1−i = 1− λ′i and λ2 λ1, we have that λ′n−2 λ′n−3, hence we can approximate
P = K∑ p=1 ( n−1∑ i=0 φ′iφ ′ i(λ ′ i) p) ≈ K∑ p=1 ( n−1∑ i=n−2 φ′iφ ′T i (λ ′ i) p) +O(λ′n−3) =
= K∑ p=1 ( 1∑ i=0 φiφ T i (1− λi)p) +O(1− λ2) = K∑ p=1 (φ0φ T 0 + φ1φ T 1 (1− λ1)p) +O(1− λ2)
= Kφ0φ T 0 + κφ1φ T 1 +O(1− λ2) where κ = ∑K p=1(1− λ1)p is a positive constant.
Now we are going to show that the farthest nodes with respect to the K-walk distance are the ones associated with the highest and lowest value of φ1.
Indeed if we want to choose i, j to be at the farthest distance we need to minimise
(P )i,j = (Kφ0φ T 0 + κφ1φ T 1 )i,j =
K n + κφ1(i)φ1(j)
which is minimum when φ1(i)φ1(j) is minimum.
We are going to show that exist p, q such that φ1(p) < 0,φ1(q) > 0. Since the eigenvector is nonzero, without loss of generality assume φ1(0) 6= 0. Since φ0 and φ1 are eigenvectors associated with different eigenvalues of a real symmetric matrix, they are orthogonal:
n−1∑ i=0 φ0(i) · φ1(i) = 0
and since φ0 is constant the previous equation leads to
n−1∑ i=0 φ1(i) = 0⇐⇒ φ1(0) = − n−1∑ i=1 φ1(i)
If such p, q didn’t exist then we would get that ∀i, j φ1(i) ·φ1(j) ≥ 0, hence multiplying both sides of the previous equation by φ1(0) we get
φ1(0) 2 = − n−1∑ i=1 φ1(i) · φ1(0)⇒ φ1(0)2 ≤ 0
Which is a contradiction since by assumption φ1(0) > 0; hence exist p, q such that φ1(p) < 0,φ1(q) > 0.
Since φ1 attains both positive and negative values, the quantity φ1(i)φ1(j) is minimised when it has negative sign and highest absolute value, hence when i, j are associated with the negative and positive values with the highest absolute value: the lowest and the highest value of φ1. Hence, dK(vM , vm)− dK(va, vb) = O(1− λ2)
C.4 INFORMAL ARGUMENT IN SUPPORT OF CONJECTURE 2.4 (GRADIENT STEPS REDUCE EXPECTED HITTING TIME)
Suppose that x, y are uniformly distributed random nodes such that φi(x) < φi(y). Let z be the node obtained from x by taking one step in the direction of ∇φi, then the expected hitting time is decreased proportionally to λ−1i and
Ex,y[Q(z, y)] ≤ Ex,y[Q(x, y)]
As a reminder, the definition of a gradient step is given in the definition 3, copied below.
Suppose the two neighboring nodes x and z are such that φ(z) − φ(x) is maximal among the neighbors of x, then we will say z is obtained from x by taking a step in the direction of the gradient ∇φ. In (Chung & S.T.Yau, 2000), it is shown the hitting time Q(x, y) is given by the equation
Q(x, y) = vol
( G(y, y)
dy − G(x, y) dx ) With λk and φk being the k-th eigenvalues and eigenvectors of the symmetric normalized Laplacian Lsym, vol the sum of the degrees of all nodes, dx the degree of node x and G Green’s function for the graph
G(x, y) = d 1 2 x d −1 2 y ∑ k>0 1 λk φk(x)φk(y)
Since the sign of the eigenvector is not deterministic, the choice φi(x) < φi(y) is used to simplify the argument without having to consider the change in sign.
Supposing λ1 λ2, the first term of the sum ofG has much more weight than the following terms. With z obtained from x by taking a step in the direction of the gradient of φ1 we have
φ1(z)− φ1(x) > 0
We want to show that the following inequality holds
Ex,y(Q(z, y)) < Ex,y(Q(x, y))
this is equivalent to the following inequality
Ex,y[G(z, y)] > Ex,y[G(x, y)]
By the hypothesis λ1 λ2, we can approximate G(x, y) ∼ d 1 2 x d −1 2 y
1 λ1 φ1(x)φ1(y) so the last
inequality is equivalent to Ex,y [ d 1 2 z d −1 2 y 1
λ1 φ1(z)φ1(y)
] > Ex,y [ d 1 2 x d −1 2 y 1
λ1 φ1(x)φ1(y) ] Removing all equal terms from both sides, the inequality is equivalent to
Ex,y [ d 1 2 z φ1(z) ] > Ex,y [ d 1 2 xφ1(x) ] But showing this last inequality is not easy. We know that φ1(z) > φ1(x) and from the choice of z being a step in the direction of∇φ1, we know it is less likely to be on the border of the graph so we believe E(dz) ≥ E(dx). Thus we also believe that the conjecture should hold in general. We believe this should be true even without the assumption on λ1 and λ2 and for more eigenvectors than φ1.

heading: C.5 PROOF FOR LEMMA C.1 (COSINE EIGENVECTORS)
text: Consider the lattice graph Γ of size N1 ×N2 × ...×Nn, that has vertices ∏ i=1,...,n{1, ..., Ni} and the vertices (xi)i=1,...,n and (yi)i=1,...,n are connected by an edge iff |xi − yi| = 1 for one index i and 0 for all other indices. Note that there are no diagonal edges in the lattice. The eigenvector of the Laplacian of the grid L(Γ) are given by φj .
Lemma C.1 (Cosine eigenvectors). The Laplacian of Γ has an eigenvalue 2− 2 cos ( π Ni ) with the associated eigenvector φj that depends only the variable in the i-th dimension and is constant in all others, with φj = 1N1 ⊗ 1N2 ⊗ ...⊗ x1,Ni ⊗ ...⊗ 1Nn , and x1,Ni(j) = cos ( πj n − π 2n
) Proof. First, recall the well known result that the path graph on N vertices PN has eigenvalues
λk = 2− 2 cos ( πk
n ) with associated eigenvector xk with i-th coordinate
xk(i) = cos
( πki
n + πk 2n
)
The Cartesian product of two graphs G = (VG, EG) and H = (VH , EH) is defined as G × H = (VG×H , EG×H) with VG×H = VG × VH and ((u1, u2), ((v1, v2)) ∈ EG×H iff either u1 = v1 and (u2, v2) ∈ EH or (u1, v1) ∈ VG and u2 = v2. It is shown in (Fiedler, 1973) that if (µi)i=1,...,m and (λj)j=1,...,n are the eigenvalues of G and H respectively, then the eigenvalues of the Cartesian product graphG×H are µi+λj for all possible eigenvalues µi and λj . Also, the eigenvectors associated to the eigenvalue µi + λj are ui ⊗ vj with ui an eigenvector of the Laplacian of G associated to the eigenvalue µi and vj an eigenvector of the Laplacian of H associated to the eigenvalue λj .
Finally, noticing that a lattice of shape N1 × N2 × ... × Nn is really the Cartesian product of path graphs of length N1 up to Nn, we conclude that there are eigenvalues 2 − 2 cos ( π Ni ) . Denoting by 1Nj the vector in R Nj with only ones as coordinates, then the eigenvector associated to the
eigenvalue 2− 2 cos ( π Ni ) is
1N1 ⊗ 1N2 ⊗ ...⊗ x1,Ni ⊗ ...⊗ 1Nn
where x1,Ni is the eigenvector of the Laplacian of PNi associated to its first non-zero eigenvalue. 2− 2 cos ( π Ni ) .

heading: C.6 RADIUS 1 CONVOLUTION KERNELS IN A GRID
text: In this section we show any radius 1 convolution kernel can be obtained as a linear combination of the Bdx(∇φi) and Bav(∇φi) matrices for the right choice of Laplacian eigenvectors φi. First we show this can be done for 1-d convolution kernels. Theorem C.2. On a path graph, any 1D convolution kernel of size 3 k is a linear combination of the aggregatorsBav,Bdx and the identity I .
Proof. Recall from the previous proof that the first non zero eigenvalue of the path graph PN has associated eigenvector φ1(i) = cos(πiN − π 2N ). Since this is a monotone decreasing function in i, the i-th row of ∇φ1 will be (0, ..., 0, si−1, 0,−si+1, 0, ..., 0)
with si−1 and si+1 > 0. We are trying to solve
(aBav + bBdx + cId)i,: = (0, ..., 0, x, y, z, 0, ..., 0)
with x, y, z, in positions i− 1, i and i+ 1. This simplifies to solving
a 1 ‖s‖L1 |s|+ b 1 ‖s‖L2 s+ c(0, 1, 0) = (x, y, z)
with s = (si−1, 0,−si+1), which always has a solution because si−1, si+1 > 0.
Theorem C.3 (Generalization radius-1 convolutional kernel in a grid). Let Γ be the n-dimensional lattice as above and let φj be the eigenvectors of the Laplacian of the lattice as in theorem C.1. Then any radius 1 kernel k on Γ is a linear combination of the aggregators Bav(φi),Bdx(φi) and I .
Proof. This is a direct consequence of C.2 obtained by adding n 1-dimensional kernels, with each kernel being in a different axis of the grid as per Lemma C.1. See figure 4 for a visual example in 2D.
C.7 PROOF FOR THEOREM 2.7 (GENERALIZATION RADIUS-R CONVOLUTIONAL KERNEL IN A LATTICE)
For an n-dimensional lattice, any convolutional kernel of radius R can be realized by a linear combination of directional aggregation matrices and their compositions.
Proof. For clarity, we first do the 2 dimensional case for a radius 2, then extended to the general case. Let k be the radius 2 kernel on a grid represented by the matrix
a5×5 =  0 0 a−2,0 0 0 0 a−1,−1 a−1,0 a−1,1 0
a0,−2 a0,−1 a0,0 a0,1 a0,2 0 a1,−1 a1,0 a1,1 0 0 0 a2,0 0 0  since we supposed the N1 × N2 grid was such that N1 > N2, by theorem C.1, we have that φ1 is depending only in the first variable x1 and is monotone in x1. Recall from C.1 that
φ1(i) = cos
( πi
N1 +
π
2N1 ) The vector N1π ∇ arccos(φ1) will be denoted by F1 in the rest. Notice all entries of F1 are 0 or ±1. Denote by F2 the gradient vector N2π ∇ arccos(φk) where φk is the eigenvector given by theorem C.1 that is depending only in the second variable x2 and is monotone in x1 and recall
φk(i) = cos
( πi
N2 +
π
2N2
)
For a matrix B, let B± the positive/negative parts of B, ie matrices with positive entries such that B = B+ −B−. LetBr1 be a matrix representing the radius 1 kernel with weights
a3×3 =
( 0 a−1,0 0
a0,−1 a0,0 a0,1 0 a1,0 0
)
The matrix Br1 can be obtained by theorem C.3. Then the radius 2 kernel k is defined by all the possible combinations of 2 positive/negative steps, plus the initial radius-1 kernel.
Br2 = ∑
−2≤i,j≤2 |i|+|j|=2
( ai,j(F sgn(i) 1 ) |i|(F sgn(j) 2 ) |j| )
︸ ︷︷ ︸ Any combination of 2 steps + Br1︸︷︷︸ all possible single-steps
with sgn the sign function sgn(i) = + if i ≥ 0 and − if i < 0. The matrix Br2 then realises the kernel a5×5.
We can further extend the above construction to N dimension grids and radius R kernels k
∑ V={v1,v2,...,vN}∈Nn
||V ||L1≤R −R≤vi≤R︸ ︷︷ ︸
Any choice of walk V with at mostR-steps
aV
N∏ j=1 (F sgn(vj) j ) |vj |
︸ ︷︷ ︸ Aggregator following the steps defined in V
with Fj = Nj π ∇ arccosφj ,φj the eigenvector with lowest eigenvalue only dependent on the j-th variable and given in theorem C.1 and ∏
is the matrix multiplication. V represents all the choices of walk {v1, v2, ..., vn} in the direction of the fields {F1,F2, ...,Fn}. For example, V = {3, 1, 0,−2} has a radius R = 6, with 3 steps forward of F1, 1 step forward of F2, and 2 steps backward of F4.

heading: C.8 PROOF FOR THEOREM 2.8 (COMPARISON WITH 1-WL TEST)
text: DGNs using the mean aggregator, any directional aggregator of the first eigenvector and injective degree-scalers are strictly more powerful than the 1-WL test.
Proof. We will show that (1) DGNs are at least as powerful as the 1-WL test and (2) there is a pair of graphs which are not distinguishable by the 1-WL test which DGNs can discriminate.
Since the DGNs include the mean aggregator combined with at least an injective degree-scaler, Corso et al. (2020) show that the resulting architecture is at least as powerful as the 1-WL test.
Then, to show that the DGNs are strictly more powerful than the 1-WL test it suffices to provide an example of a pair of graphs that DGNs can differentiate and 1-WL cannot. Such a pair of graphs is illustrated in figure 8.
The 1-WL test (as any MPNN with, for example, sum aggregator) will always have the same features for all the nodes labelled with a and for all the nodes labelled with b and, therefore, will classify the graphs as isomorphic. DGNs, via the directional smoothing or directional derivative aggregators based on the first eigenvector of the Laplacian matrix, will update the features of the a nodes differently in the two graphs (figure 8 presents also the aggregation functions) and will, therefore, be capable of distinguishing them.

Review summaries:
The authors propose a convolution as a message passing of node features over edges where messages are aggregated weighted by a "direction" edge field. Furthermore, the authors propose to use the gradients of Laplace eigenfunctions as direction fields. Presumably, the aggregation is done with different direction fields derived from the Laplace eigenfunctions with lowest eigenvalues, which are then linearly combined with learnable parameters. Doing so allows their graph network to behave more like a conventional CNN, in which the kernels have different parameters for signals from different directions. The authors achieve good results on several benchmarks. Furthermore, the authors prove that their method reduces to a conventional CNN on a rectangular grid and have theoretical results that suggest that their method suffers less from the "over-smoothing" and "over-squashing" problems.