PDE-GCN: Novel Architectures for Graph Neural Networks Motivated by Partial Differential Equations
Year: 2021

Abstract:
Graph neural networks are increasingly becoming the go-to approach in various fields such as computer vision, computational biology and chemistry, where data are naturally explained by graphs. However, unlike traditional convolutional neural networks, deep graph networks do not necessarily yield better performance than shallow graph networks. This behavior usually stems from the over-smoothing phenomenon. In this work, we propose a family of architectures to control this behavior by design. Our networks are motivated by numerical methods for solving Partial Differential Equations (PDEs) on manifolds, and as such, their behavior can be explained by similar analysis. Moreover, as we demonstrate using an extensive set of experiments, our PDE-motivated networks can generalize and be effective for various types of problems from different fields. Our architectures obtain better or on par with the current state-of-the-art results for problems that are typically approached using different architectures.

heading: 1 Introduction
text: In recent years, Graph Convolutional Networks (GCNs) [1, 2, 3] have drawn the attention of researchers and practitioners in a variety of domains and applications, ranging from computer vision and graphics [4, 5, 6, 7, 8] to computational biology [9, 10, 11], recommendation systems [12] and social network analysis [13, 14]. However, GCNs still suffer from two main problems. First, they are usually shallow as opposed to the concept of deep convolutional neural networks (CNNs) [15, 16] due to the over-smoothing phenomenon [17, 18, 19], where the node feature vectors become almost identical, such that they are indistinguishable, which yields non-optimal performance. Furthermore, GCNs are typically customized to a specific domain and application. That is, as we demonstrate in Sec. 4.1, a successful point-cloud classification network [5] can perform poorly on a citation graph node-classification problem [19], and vice-versa. Furthermore, because many GCNs lack theoretical guarantees, it is difficult to reason about their success on one problem and lack on others. These observations motivate us to develop a profound understanding of graph networks and their dynamics.
To this end, we suggest a novel, universal approach to the design of GCN architectures. Our inspiration stems from the similarities and equivalence between Partial Differential Equations (PDEs) and deep networks explored in [20, 21, 22]. Furthermore, as GCNs can be thought of as a generalization of
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
CNNs, and a standard convolution can be represented as a combination of differential operators on a structured grid [22], we adopt this interpretation to explore versions of GCNs as PDEs on graphs or manifolds. We therefore call our network architectures PDE-GCN, and demonstrate that our approach is general with respect to the given task. That is, our architectures behave similarly for different problems, and their performance is on par or better than other domain-specific GCNs. Furthermore, our family of architectures are backed by theoretical guarantees that allow us to explain the behavior of the GCNs in some of the results that we present. To be more specific, our contribution is as follows:
• We introduce and implement general graph convolution operators, based on graph gradient and divergence. This abstraction of the spatial operation on the graph leads to a more general and flexible approach for architecture design.
• We propose treating a variety of graph related problems as discretized PDEs, and formulate the dynamics that match different problems such as node-classification and dense shape-correspondence. This is in direct effort to propose a family of networks that can solve multiple problems, instead of GCNs which are tailored for a specific application.
• Our method allows constructing a deep GCN without over-smoothing, with theoretical guarantees. • We validate our method by conducting numerous experiments on multiple datasets and applications,
achieving significantly better or similar accuracy compared to state-of-the-art models.

heading: 2 Related work
text: Graph Convolutional Networks: GCNs are typically divided into spectral [1, 3, 2] and spatial [23, 24, 25, 4] categories. Most of those can be implemented using the Message-Passing Neural Network paradigm [25], where each node aggregates features (messages) from its neighbors, according to some scheme. The works [3, 2] use polynomials of the graph Laplacian to parameterize the convolution operator. DGCNN [5] constructs a k-nearest-neighbors graph from point-clouds and dynamically updates it. MoNet [4] learns a Gaussian mixture model to weight the edges of meshes for shape analysis tasks. Works like GraphSAGE [26] and GAT [27] propose methods for inductive and transductive learning on non-geometrical graphs.
Several of the methods above suffer from over-smoothing [17, 19], leading to undesired node features similarity for deep networks. To overcome this problem, some approaches rely on imposing regularization and augmentation. For example, PairNorm [17] introduces a novel normalization layer, and DropEdge [28] randomly removes edges to decrease the degree of nodes. Other methods prevent over-smoothing by dedicated construction. For instance, JKNet [29] combines all intermediate representations at each stage of the network. APPNP [30] replaces learnt convolutional layers with a pre-defined kernel based on PageRank, yielding a shallow network which preserves locality, making it robust to over-smoothing. GCNII [19] proposes to add an identity residual where the initial features of the networks are added to the features of the l-th layer, scaled by some coefficient.
Another approach is to construct a network that inherently does not over-smooth, as we suggest in this work. Our network is based on discretized PDEs, hence we are able to motivate our choices by well studied theory and numerical experiments [31]. On a similar note, the recent DiffGCN [8] also makes use of discretized operators to approximate the graph gradient and Laplacian. However, DiffGCN is specifically tailored for geometric tasks since it projects its components on the x, y, z axes, and is applied using a ResNet [16] (diffusion) structure only. The recent GRAND [32] applies attention mechanism with diffusive dynamics, using several integration schemes. Here we propose a network for both geometric and non-geometric tasks, and also utilize both the diffusion or hyperbolic layer dynamics, including a mixture of the two.
PDEs and CNNs: In a recent series of works, the connection between PDEs and CNNs was studied [20, 21, 22, 33, 34, 35, 36, 37]. It was shown that it is possible to treat a deep neural network as a dynamical system driven by some PDE, where each convolution layer is considered a time step of a discretized PDE. The connection between PDEs and CNNs was also used to reduce the computational burden [38]. Besides the interpretation of CNNs as a PDE solver, it was shown that it is also possible to learn a symbolic representation of PDEs in a task-and data driven approach in [39, 40]. In the context of GCNs and PDEs, it was recently shown [41, 42] that GCNs can be utilized to enhance the solution of PDEs for problems like fluid flow dynamics. However, in this work, we harness PDE concepts to design and construct GCNs for a variety of applications.

heading: 3 Methods

heading: 3.1 Partial Differential Equations on manifolds
text: We show now that GCNs can be viewed as discretizations of PDEs on manifolds, similarly to CNNs that are viewed as discretized PDEs on regular grids [22]. Consider a general manifoldM where a vector function f resides (also dubbed the features function), along with its continuous differential operators such as the gradient ∇, divergence ∇· and the Laplacian ∆ that reside on the manifoldM. Given these differential operators, one can model different processes onM. In particular, we consider two PDEs – the non-linear diffusion and the non-linear hyperbolic equations
ft = ∇ ·K∗σ(K∇f), f(t = 0) = f0, t ∈ [0, T ], (1) ftt = ∇ ·K∗σ(K∇f), f(t = 0) = f0, ft(t = 0) = 0, t ∈ [0, T ], (2)
respectively, equipped with appropriate boundary conditions. Here K is a coefficient matrix that can change in time and represents the propagation over the manifoldM, K∗ is its conjugate transpose and σ(·) is a non-linear activation function. Eq. (1)-(2) define a non-linear operator that takes initial feature vectors f0 at time 0 and propagates them to time T , yielding fT where they can be used for different tasks. We now provide two theorems that characterize the behavior of Eq. (1)-(2), based on ideas from [22]1. Theorem 1. If the activation function σ(·) is monotonically non-decreasing and sign-preserving, then the forward propagation through the diffusive PDE in (1) for t ∈ [0,∞) yields a non-increasing feature norm, that is,
∂ ∂t ‖f‖2 ≤ 0.
Theorem 2. Assume that the activation function σ(·) is monotonically non-decreasing, signpreserving and satisfies |σ(x)| ≤ |x|, and define energy
Enet = ‖ft‖2 + (K∇f, σ(K∇f)) , 1See proofs in Appendix A.
then the forward propagation through the hyperbolic PDE in (2) satisfies Enet ≤ cK , where cK is a constant that depends on K but independent of time.
The outcome of those theorems is that the dynamics described in Eq. (1) is smoothing, while the one in Eq. (2) is bounded by a conserving mapping. An illustration of this behavior is presented in Fig. 1.
In the physical world, diffusion and hyperbolic equations are used for different applications. Similarly, many computational models for image segmentation [43], denoising [44] and deblurring are based on anisotropic diffusion which are similar to the model in Eq. (1). On the other hand, applications that require conservation such as volume/distance preservation as in the dense shape correspondence task [45] and protein folding [11], are typically better treated using a hyperbolic equation as in Eq. (2). Those insights motivate us to construct two types of layers according to Eq. (1)-(2) using discretized differential operators on graphs.

heading: 3.2 Discretized differential operators on graphs
text: The models (1) and (2) reside in a continuous manifoldM, on which a continuous function vector f is defined. A graph can be thought of as a discretization of that manifold to a finite space. Assume we are given an undirected graph G = (V, E) where V ∈ M is the set of n vertices of the graph and E is the set of m edges of the graph. Let us denote by f i ∈ Rc the value of the discrete version of f , on the i-th node of G. c is the number of channels, which is the width of the neural network. We define G, the discrete gradient operator on the graph, also known as the incidence matrix, as follows:
(Gf)ij = Wij(fi − fj), (3)
where nodes i and j are connected via the (i, j)-th edge, Wij is an edge weight matrix which can be learnt, and fi and fj are the features on the i-th and j-th nodes, respectively. The gradient operator can be thought of as a weighted directional derivative of the function f in the direction defined by the nodes i and j. Furthermore, if we choose the scaled identity matrix Wij = d−1ij I, where dij is the distance between the two nodes, then the discrete gradient is a second order approximation to the true gradient of the function f on the edges of the graph. In this work, we use Wij = γijI, where the scale γij is the geometric mean of the degree of the nodes i, j. Note, that the gradient operator is a mapping from the vertex space to the edge space.
Given the gradient matrix, it is possible to define the divergence matrix [46], which is an approach that is used extensively in mimetic discretizations of PDEs. To this end, we define the inner product between an edge feature vector q and the gradient of a node feature vector f as
(q,Gf) = q>Gf = f>G>q. (4)
The divergence is naturally defined as the operator that maps edge operator q to the node space, that is ∇· ≈ −G>. As usual, the graph Laplacian operator can be obtained by taking the divergence of the gradient. In graph theory it is defined as a positive matrix that is, ∆ ≈ G>G 2. We also define the weighted line integral over an edge. Similarly to Eq. (3)-(4), we define
(Af)ij = 1
2 Wij(fi + fj), (q,Af) = q
>Af = f>A>q. (5)
The operator A approximates the mass operator on the edges. The right equation in Eq. (5) suggests that an appropriate averaging operator for edge features is the transpose of the nodal edge average.
The advantage of defining such operators is that we are able to design networks with architectures that mimic the continuous operators (1) and (2) on the discrete level, as we show in the next section.

heading: 3.3 PDE-GCN: Graph Convolutional Networks by Partial Differential Equations
text: In order to use the computational models in Eq. (1)-(2), we form their discrete versions:
f (l+1) = f (l) − hG>K>l σ(KlGf (l)), (6) f (l+1) = 2f (l) − f (l−1) − h2G>K>l σ(KlGf (l)). (7)
2In the field of numerical PDEs, the Laplacian is defined as −G>G, but the combinatorial Laplacian is defined as G>G.
Here, in Eq. (6) we use the forward Euler to discretize Eq. (1), and in Eq. (7) we discretize the second order time derivative in Eq. (2), using the leapfrog method. In both cases, f (l) are the node features and Kl is a 1 × 1 trainable convolution of the l-th layer. The similarity between (6) and ResNet is well documented in the context of CNNs [22]. The hyper-parameter h is the step-size, and it is chosen such that the stability of the discretization is kept. We use σ = tanh for the activation function as it yields slightly better results in our experiments, although other functions such as ReLU can also be used, as reported in Sec. 4.7. Each of Eq. (6)-(7) defines a PDE-GCN block. We denote the former (diffusive equation) by PDE-GCND and the latter (hyperbolic equation) by PDE-GCNH.
To complete the description of our network, a few more details are required, as follows:
The (opening) embedding layer. The input vertex features uV are fed through an embedding (1× 1 convolution) layer Ko to obtain the initial features f0 of our PDE-GCN network: f0 = KouV . In cases where input edge attributes (features) uE are available (as in the experiment in Sec. 4.6), we transform them to the vertex space by taking their divergence and average, and concatenate them to the input of the embedding layer Ko as follows: f0 = Ko(uV ⊕A>uE ⊕G>uE).
The (closing) embedding layer. Given the final vertex features f (L) of our PDE-GCN with L layers, we obtain the output of the network by performing: uout = Kcf (L). Here Kc is a 1× 1 convolution layer mapping the hidden feature space to the output shape.
Initialization. The 1 × 1 convolutions Kl in Eq. (6)-(7) are initialized as identity. This way, the network begins from a diffusion/hyperbolic equation, which serves as a prior and guides the network to initially behave like classical methods [43, 45] and to further improve during training.
The choice of dynamics. For some applications, anisotropic diffusion is appropriate, while for others conservation is more important. However, in some applications this may not be clear. To this end, it is possible to combine Eq. (1)-(2) to obtain the continuous process
αftt + (1− α)ft = ∇ ·K∗σ(K∇f), f(t = 0) = f0, ft(t = 0) = 0 t ∈ [0, T ], (8)
where α = sigmoid(β), meaning 0 ≤ α ≤ 1, and β is a single trainable parameter. The discretization of this PDE leads to the following network dynamics:
α(f (l+1) − 2f (l) + f (l−1)) + h(1− α)(f (l+1) − f (l)) = −h2G>K>l σ(KlGf (l)), (9)
where f (l+1) is updated by the known f (l) and f (l−1). We denote a layer that is governed by Eq. (9) by PDE-GCNM (where M stands for mixture). Note, that it is also possible to learn a combination coefficient αi per layer, although we did not read a benefit from such scheme.
As we show in our numerical experiments, learning α yields results that are consistent with our understanding of the problem, that is, graph node-classification gravitates towards no second order derivatives while applications that require conservation gravitate towards the hyperbolic equation.

heading: 4 Experiments
text: In this section we demonstrate our approach on various problems from different fields and applications ranging from 3D shape-classification [47] to protein-protein interaction [26] and node-classification [48]. The experiments also vary in their output type. Node classification is similar to segmentation problems that are typically solved by anisotropic diffusion while the dense shape correspondence problem is conservative and therefore can be thought of as a problem that does not require smoothing.
In all experiments, we use the suitable PDE-GCN (D, H or M) block as described in Sec. 3, with various depths (number of layers) and widths (number of channels), as well as the appropriate final convolution steps, depending on the task at hand. A detailed description of the architectures used in our experiments is given in Appendix B. We use the Adam [49] optimizer in all experiments, and perform grid search over the hyper-parameters of our network. The selected hyper-parameters are reported in Appendix C. Our objective function in all experiments is the cross-entropy loss, besides inductive learning on PPI where we use the binary cross-entropy loss. Our code is implemented using PyTorch [50], trained on an Nvidia Titan RTX GPU.
We show that for all the considered tasks and datasets, our method is either remarkably better or on par with state-of-the-art models.

heading: 4.1 Generalization of GCNs to different applications
text: To gain deeper understanding about the effectiveness and generalization of various GCN methods to different tasks, we start by picking two datasets - ModelNet-10 [47] for 3D shape-classification, and Cora [48] for semi-supervised node classification. Those datasets are not only different in terms of application (global versus local classification), but also stem from different domains. While in ModelNet-10 the data has geometrical meaning, the data in Cora has no obvious geometrical interpretation. Therefore, we suggest that success in both applications should be obtained from a generalizable GCN. We compare our PDE-GCND with two recent and popular networks - DGCNN [5] and GCNII [19]. For ModelNet-10 shape-classification, we randomly sample 1,024 points from each shape to form a point cloud, and connect its points using k-nearest-neighbors (k-nn) algorithm with k = 10 to obtain a graph and follow the training scheme of [5]. On Cora, we follow the same procedure as in [19]. We evaluate all models with 4 layers, as well as 2 layers for DGCNN on Cora.
Our results, reported in Tab. 1 suggest that while each of the considered methods obtains high accuracy on the dataset it originally was tested on (Cora for GCNII and ModelNet-10 for DGCNN), obtaining a similar measure of success on a different dataset was not possible when using the very same networks. Additionally, while our attempts to add a diffusion-equation dynamics to DGCNN (i.e., updating features as in Eq. (6)) showed an increase in performance – a large gap to state-of-theart model still exists. On top of that, we also see that DGCNN suffers from over-smoothing, as its accuracy significantly decreases when adding more layers. Last but not least, we observe that our PDE-GCND obtains high accuracy on both datasets, similar or better than state-of-the-art models.

heading: 4.2 Learning PDE network dynamics
text: In this experiment, we delve on the ability to learn the appropriate PDE that better models a given problem. To this end, we use the mixture model from Eq. (9) so that the resulting PDE is a combination of the diffusion and hyperbolic dynamics. We use a 8 layer mixed PDE-GCN, starting with α = 0.5, such that it is balanced between a PDE-GCND and PDE-GCNH. By learning the parameter α in (9), we allow to choose a mixed PDE between a purely conservative network and a
diffusive one. We consider two problems: semi-supervised node classification on Cora, and dense shape correspondence on FAUST [51].
Our results, reported in Fig. 2 suggest that just as in classical works [43, 52], problems like nodeclassification obtain better performance with an anisotropic diffusion like in Eq. (6), and for problems involving dense-correspondences like in [45, 11] that tend to conserve the energy of the underlying problem, a hyperbolic equation type of PDE as in Eq. (7) is more appropriate.

heading: 4.3 Semi-supervised node classification
text: In this set of experiments we use three datasets – Cora, Citeseer and PubMed [48]. For all datasets we use the standard training/validation/testing split as in [53], with 20 nodes per class for training, 500 validation nodes and 1,000 testing nodes and follow the training scheme of [19]. The statistics of the datasets are reported in Tab. 2. We compare our results using PDE-GCND with recent and popular models like GCN [3], GAT [27], APPNP [30], JKNet [29] and DropEdge [28]. We note that our network does not over-smooth, as an increase in the number of layers does not cause performance degradation. For example, on CiteSeer, we obtain 75.6% accuracy with 32 layers, compared to 74.6% with two layers3. Overall, our results in Tab. 3 show that our PDE-GCND achieves similar or better accuracy than the considered methods.

heading: 4.4 Fully-supervised node classification
text: We follow [54] and use 7 datasets: Cora, CiteSeer, PubMed, Chameleon, Cornell, Texas and Wisconsin. We also use the same train/validation/test splits of 60%, 20%, 20%, respectively. In addition, we report the average performance over 10 random splits from [54]. We fix the number of channels to 64 and perform grid search to determine the hyper-parameters parameters, which are reported in Appendix C. We compare our network with GCN, GAT, three variants of Geom-GCN [54], APPNP, JKNet, Incep and GCNII in Tab. 4. Our experiments read either similar or better than the state-of-the-art on Cora, CiteSeer and PubMed datasets. On Chameleon [55], Cornell, Texas and Wisconsin datasets, we improve state-of-the-art accuracy by a significant margin. For example, we obtain 93.24% accuracy on Texas with our PDE-GCNM, compared to 77.84% with GCNII*. Similar improvements hold for Cornell and Wisconsin datasets. The common factor for these datasets is their small size, as depicted from Tab. 2. We argue that the success of our network stems from its capability of apriori extracting features from graphs, due to its utilization of discretized differential operators and PDE guided construction. On Chameleon [55], using PDE-GCNM, we improve the current state-of-the-art accuracy of GCNII* from 62.48% to 66.01%. Also, we note that unlike in the semi-supervised case, where some of the labels are missing, here it is possible to obtain meaningful results with the hyperbolic equation based PDE-GCNH as we do not have unknown nodes in the fully-supervised case, which would be otherwise preserved using the hyperbolic equation dynamics.

heading: 4.5 Inductive Learning
text: We follow [19] and employ the PPI dataset [26] for the inductive learning task. We use a 8 layer PDE-GCND network, without dropout or weight-decay, and a learning rate of 0.001. We compare our results with methods like GraphSAGE, GAT, JKNet, GeniePath, GCNII and others. As reported in Tab. 5, our PDE-GCN D achieves 99.07 Micro-averaged F1 score, superior to methods like GAT, JKNet and GeniePath, also close to state-of-the-art GCNII* with a score of 99.58.
3Note that this result is also a new state-of-the-art accuracy.

heading: 4.6 Dense shape correspondence
text: Finding dense correspondences between shapes is a classical experiment for hyperbolic dynamics, as we are interested in modeling local motion dynamics. In essence, learning to find correspondences between shapes is similar to learning a transformation from one shape to the other. To this end, we use the FAUST dataset [51] containing 10 scanned human shapes in 10 different poses with 6,890 nodes each. We follow the train and test split from [4], where the first 80 subjects are used for training and the remaining 20 subjects for testing. Our metric is the correspondence percentage with zero geodesic error, i.e., the percentage of perfectly matched shapes from all our test cases. We follow the pre-processing and training scheme of [62] where we use Cartesian coordinates to describe distances between nodes, with initial features of a constant 1 ∈ Rn where n is the number of nodes. We use a 8 layer PDE-GCND and PDE-GCNH variants, both with constant learning rate of 0.001 with no weight-decay or dropout, and compare to recent and popular methods like ACNN [60], MoNet [4], FMNet [61], and SplineCNN [62]. As expected from the discussion in Sec. 3.1, and reported in Tab. 6, the hyperbolic equation proves to be a better fit for this kind of problem. Also, our PDE-GCNH achieves a promising 99.9% correspondence rate with zero geodesic error, outperforming the rest of the considered methods.

heading: 4.7 Ablation study
text: Our method has two main dynamics – diffusion and hyperbolic. To verify our proposal in Sec. 3.1, we examine the performance of the hyperbolic PDE-GCNH on semi-supervised node-classification on Cora and CiteSeer datasets. Our results in Tab. 7 show that indeed for problems where we wish to obtain a piecewise-constant prediction, the diffusive formulation of our network, PDE-GCND, is
Table 5: Protein-protein interaction (PPI). Results are reported in micro-averaged F1 score.
Method Micro-averaged F1
GraphSAGE [26] 61.20 VR-GCN [56] 97.80 GaAN [57] 98.71 GAT [27] 97.30 JKNet [29] 97.60 GeniePath [58] 98.50 Cluster-GCN [59] 99.36 GCNII [19] 99.54 GCNII* [19] 99.58 PDE-GCND (Ours) 99.07 PDE-GCNM (Ours) 99.18
Table 6: Dense shape correspondence (%) with zero geodesic error
Method Faust
ACNN [60] 63.8 MoNet [4] 89.1 FMNet [61] 98.2 SplineCNN [62] 99.2
PDE-GCND (Ours) 64.2 PDE-GCNH (Ours) 99.9
more suitable. Furthermore, we study the importance of the positive-semi definiteness of our learnt operator as described in Eq. (6)-(7), by removing the KTl term from the dynamics equations. This yields a non-symmetric operator that does not guarantee positive-semi-definiteness. We note that the enforcement of the latter is important to obtain higher accuracy which is improved by up to 3% with the introduction of a positive semi-definite operator. Also, we report on the use of σ = ReLU , from the discussion in Sec. 3.3, where we favor tanh as an activation function, and the ability of our PDE-GCNM to reproduce the results of PDE-GCND in the case of semi-supervised learning.

heading: 5 Summary
text: In this paper we explored new architectures for graph neural networks. Our motivation stems from the similarities between graph networks and time dependent partial differential equations that are discretized on manifolds and graphs. By adopting an appropriate PDE, and embedding the finite graph in an infinite manifold, we are able to define networks that are either diffusive, conservative, or a combination of both.
Not all natural phenomena are solved using the same PDE and we should not expect that all graph problems should be solved by the same network dynamics. To this end we allow the data to choose which type of network is appropriate for the solution of the problem (diffusive or hyperbolic). Indeed, numerical experiments show that the network gravitates towards a hyperbolic one for problems where conservation is required, and towards a diffusive one when anisotropic diffusion is favorable.
Finally, we showed that the proposed networks can be made deep without over-smoothing and, can deliver the state-of-the-art performance or improve it for virtually every problem we worked with. In particular, our network dramatically improved the state-of-the-art for problems that are data-poor. We believe that for such problems the structure imposed by our dynamics and operators regularizes the network and therefore yields implicit regularization.
Acknowledgments and Disclosure of Funding
The research reported in this paper was supported by grant no. 2018209 from the United States - Israel Binational Science Foundation (BSF), Jerusalem, Israel. ME is supported by Kreitman High-tech scholarship.

Review summaries:
Motivated by the connections between PDEs and conventional convolutional networks, this paper proposes the view of graph convolutional networks (GCNs) as discretizations of PDEs on manifolds. Using this view, the authors propose a a family of graph convolutional architectures that utilize discrete graph operators as their basis and argue that this approach can address existing issues with graph convolutional networks, such as over-smoothing. Moreover, the authors also argue that this approach can have a more even performance accorss a diverse range of tasks, instead of being specialized to only a specific type of problem.

The authors propose a new class of graph neural networks motivated by numerical methods for solving PDEs on manifolds. The idea is that there is a correspondence between specific GNNs (e.g., GNNs that learn the weight of edges such as GATs) and discretized PDEs on manifolds. Here, the discretization is spatially. For time (i.e., the depth of the GNN) one can use numerical solvers for PDEs similar to what was proposed by the original Neural ODE paper. The main motivation is that there are types of PDEs that do not suffer from the oversmoothing problem and, therefore, that their discretized versions might be a better fit for some problem domains where oversmoothing is a problem.