From Continuous Dynamics to Graph Neural Networks: Neural Diffusion and
  Beyond
Graph neural networks (GNNs) have demonstrated significant promise in
modelling relational data and have been widely applied in various fields of
interest. The key mechanism behind GNNs is the so-called message passing where
information is being iteratively aggregated to central nodes from their
neighbourhood. Such a scheme has been found to be intrinsically linked to a
physical process known as heat diffusion, where the propagation of GNNs
naturally corresponds to the evolution of heat density. Analogizing the process
of message passing to the heat dynamics allows to fundamentally understand the
power and pitfalls of GNNs and consequently informs better model design.
Recently, there emerges a plethora of works that proposes GNNs inspired from
the continuous dynamics formulation, in an attempt to mitigate the known
limitations of GNNs, such as oversmoothing and oversquashing. In this survey,
we provide the first systematic and comprehensive review of studies that
leverage the continuous perspective of GNNs. To this end, we introduce
foundational ingredients for adapting continuous dynamics to GNNs, along with a
general framework for the design of graph neural dynamics. We then review and
categorize existing works based on their driven mechanisms and underlying
dynamics. We also summarize how the limitations of classic GNNs can be
addressed under the continuous framework. We conclude by identifying multiple
open research directions.
Andi Han, Dai Shi, Lequan Lin, Junbin Gao
# From Continuous Dynamics to Graph Neural Networks:

###### Abstract

Graph neural networks (GNNs) have demonstrated significant promise in modelling relational data and have been widely applied in various fields of interest. The key mechanism behind GNNs is the so-called message passing where information is being iteratively aggregated to central nodes from their neighbourhood. Such a scheme has been found to be intrinsically linked to a physical process known as heat diffusion, where the propagation of GNNs naturally corresponds to the evolution of heat density. Analogizing the process of message passing to the heat dynamics allows to fundamentally understand the power and pitfalls of GNNs and consequently informs better model design. Recently, there emerges a plethora of works that proposes GNNs inspired from the continuous dynamics formulation, in an attempt to mitigate the known limitations of GNNs, such as oversmoothing and oversquashing. In this survey, we provide the first systematic and comprehensive review of studies that leverage the continuous perspective of GNNs. To this end, we introduce foundational ingredients for adapting continuous dynamics to GNNs, along with a general framework for the design of graph neural dynamics. We then review and categorize existing works based on their driven mechanisms and underlying dynamics. We also summarize how the limitations of classic GNNs can be addressed under the continuous framework. We conclude by identifying multiple open research directions.

## 1 Introduction

Graph neural networks (GNNs) [58, 91, 38] have emerged as one of the most popular choices for processing and analyzing relational data. The main goal of GNNs is to acquire expressive representations at node, edge or graph level, primarily for tasks including but not limited to node classification, link prediction and graph property prediction. Such a learning process requires not only utilizing the node features but also the underlying graph topology. GNNs have been successful in various fields of application where relational structures are critical for learning quality, including recommendation systems [31], transportation networks [54], particle system [82], molecule and protein designs [50, 95], neuroscience [8] and material science [100], among many others.

The effectiveness of GNNs primarily roots in the message passing paradigm, in which information is iteratively aggregated among the neighbouring nodes to update representation of the center node. One prominent example is the graph convolution network (GCN) [58] where each layer updates node representation by taking a degree-weighted average of its neighbours' representation. Although being effective in capturing dependencies among the nodes, the number of layers of message passing requires to be carefully chosen so as to avoid performance degradation. This is unlike the classic (feedforward) neural networks where increasing
depth generally leads to improved predictive performance. In particular, due to the nature of message passing, deeper GNNs have a tendency to over-smooth the features, leading to un-informative node representations. On the other hand, shallow GNNs are likely to suffer from information bottleneck where signals from distant nodes (with regards to graph topology) exert little influence on the centered node.

In order to analyze and address the aforementioned issues, many recent works have cast GNNs as discretization of certain continuous dynamical systems, by analogizing propagation through the layers to time evolution of dynamics. Indeed, the idea of viewing neural networks in the continuous limits is not new and has been explored for classic neural networks [27, 45, 65, 61, 80]. The continuous formulation provides a unified framework for understanding and designing the propagation of neural networks, leveraging various mathematical tools from control theory and differential equations. One seminal work [19] proposes neural ordinary differential equation (Neural ODE), which can be viewed as a continuous-depth residual network [48]. Another work explores the connections between the convolutional networks with partial differential equations (PDEs) [80].

In GNNs, in particular, the message passing mechanism can be viewed as realization of a class of PDEs, namely the heat diffusion equation [15]. Such a novel perspective allows to better dissect the behaviours of GNNs. As an example, because heat diffusion is known to dissipate energy and converge to a steady state of heat equilibrium, the phenomenon of oversmoothing corresponds to the equilibrium state of diffusion where heat distributes uniformly across spatial locations. Apart from offering theoretical insights into the evolution of GNNs, the continuous dynamics perspective also allows easy integration of structural constraints and desired properties into the dynamical system, such as energy conservation, (ir)reversibility, and boundary conditions. Meanwhile advanced numerical integrators, like high-order and implicit schemes, can be employed to enhance the efficiency and stability of discretization [13]. Finally, the variety of continuous dynamics grounded with physical substance can inform better designs of GNNs to enhance the representation power and overcome the performance limitations [15, 86, 29, 78, 68, 103, 21, 30, 42, 57].

The theory of differential equations and dynamical systems is well-developed with a long-standing history [92, 72], while graph neural network is a comparatively nascent field, garnering attention within the past decade. This offers great potential in harnessing the rich foundations of the dynamical system theory for enhancing the understanding and functionality of GNNs. In this work, we provide a comprehensive review of existing developments in continuous dynamics inspired GNNs, which we collectively refer to as graph neural dynamics. To the best of our knowledge, this is the first survey on the connections between continuous dynamical systems and graph neural networks. Nevertheless, we would like to highlight a related line of research that utilizes deep learning to solve differential equations [52, 44, 60, 59]. In particular, graph neural operators [3, 62] employ GNNs to solve PDEs by discretizing the domains and constructing graph structure based on spatial proximity. In contrast, the idealization that we focus in this work is the reverse where PDEs inform the designs of graph neural networks.

Organization.We organize the rest of the paper as follows. In Section 2, we introduce diffusion equations from first principles and then review discrete operators on graphs such as gradient and divergence, which are crucial for formulations of continuous GNNs. This section concludes with the connection between diffusion equation with the famous graph convolutional network. Section 3 then presents the framework for designing GNNs from continuous dynamics and summarizes the existing works based on the underlying physical processes. We then in Section 4 discuss the various numerical schemes for propagating and learning the continuous GNN dynamics. Finally, in Section 5, we explain how the various dynamics help to tackle the shortcomings of existing GNNs, in terms of oversmoothing, oversquashing, poor performance on heterophilic graphs (where neighbouring nodes do not share similar information), as well as adversarial robustness and
training stability.

## 2 From diffusion equations to graph neural networks

A fundamental building block of graph neural networks is message passing where information flows between neighbouring nodes. Message passing has natural connection to diffusion equations, which describe how certain quantities of interest, such as mass or heat disperse spatially, as a function of time.

Diffusion equation in continuous domains.The two laws governing any diffusion equation are the _Fick's law_ and _mass conservation law_. The former states the diffusion happens from the regions of higher concentration to regions of lower concentration and the rate of diffusion is proportional to the mass difference (the gradient). The latter states the mass cannot be created nor destroyed through the process of diffusion. Formally, let \(x:\Omega\times\mathbb{R}\rightarrow\mathbb{R}\) represent the mass distribution, i.e., \(x(u,t)\) over the spatial locations \(u\in\Omega\) and time \(t\in\mathbb{R}\). Denote \(\nabla x:=\frac{\partial x}{\partial u}\) as the gradient of mass across the space. The Fick's law states that a flux \(J\) exists towards lower mass concentrated regions, i.e., \(J=-D\nabla x\), where \(D\) is the diffusivity coefficient that potentially depends on both space and time. When \(\Omega\subseteq\mathbb{R}^{d}\), we can write \(\nabla x=[\frac{\partial x}{\partial u_{1}},...\frac{\partial x}{\partial u_{ d}}]\) and \(D\in\mathbb{R}^{d\times d}\) and \(J\) corresponds to a flux in \(d\) directions. The mass conservation law then leads to the following continuity equation \(\frac{\partial x}{\partial t}=-\mathrm{div}J\), where \(\mathrm{div}\) is the divergence operator that computes the mass changes at certain location and time. Combining the two equations yields the fundamental (heat) diffusion equation

\[\frac{\partial x}{\partial t}=\mathrm{div}(D\nabla x). \tag{1}\]

The diffusion process is called _homogeneous_ when the diffusivity coefficient \(D\) is space independent, and is called _isotropic_ when \(D\) depends only on the location, and called _anisotropic_ when \(D\) depends on both the location and the direction of the gradient. In the case of homogeneous diffusion, we can write the diffusion equation in terms of the Laplacian operator \(\Delta\coloneqq\mathrm{div}\cdot\nabla\), i.e., \(\frac{\partial x}{\partial t}=D\Delta x\).

Graphs and discrete differential operators.In order to properly characterize diffusion dynamics over graphs, it is necessary to generalize the concepts from differential geometry to the discrete space, such as gradient and divergence.

A graph can be represented by a tuple \(\mathcal{G}=(\mathcal{V},\mathcal{E})\) where \(\mathcal{V},\mathcal{E}\) denote the set of nodes and edges respectively. In this work, we focus on the undirected graphs, i.e., if \((i,j)\in\mathcal{E}\), then \((j,i)\in\mathcal{E}\). Leveraging tools from differential geometry, let \(L^{2}(\mathcal{V})\) and \(L^{2}(\mathcal{E})\) be Hilbert spaces for real-valued functions on \(\mathcal{V}\) and \(\mathcal{E}\) respectively with the inner product given by

\[\langle f,g\rangle_{L^{2}(\mathcal{V})}=\sum_{i\in\mathcal{V}}f_{i}g_{i}, \qquad\langle F,G\rangle_{L^{2}(\mathcal{E})}=\sum_{(i,j)\in\mathcal{E}}F_{ij} G_{ij}\]

for \(f,g:\mathcal{V}\rightarrow\mathbb{R}\) and \(F,G:\mathcal{E}\rightarrow\mathbb{R}\). This allows to generalize the definitions of gradient and divergence to graph domains [12]. Formally, graph gradient is defined as \(\nabla:L^{2}(\mathcal{V})\to L^{2}(\mathcal{E})\) such that \((\nabla f)_{ij}=f_{j}-f_{i}\). Here we assume the edges are anti-symmetric, namely \((\nabla f)_{ij}=-(\nabla f)_{ji}\). Graph divergence \(\mathrm{div}:L^{2}(\mathcal{E})\to L^{2}(\mathcal{V})\), is defined as the converse of graph gradient such that \((\mathrm{div}F)_{i}=\sum_{j:(i,j)\in\mathcal{E}}F_{ij}\).

By the discrete nature, graphs and functions/signals on a graph can be compactly represented via vectors and matrices. In a graph with \(n\) nodes (\(|\mathcal{V}|=n\)), each \(f\in L^{2}(\mathcal{V})\) can be written as \(n\)-dimensional vector \(\mathbf{f}\) and \(F\in L^{2}(\mathcal{E})\) can be written as a matrix \(\mathbf{F}\) of size \(n\times n\) (with nonzero \(i,j\)-th entry \(\mathbf{F}_{i,j}\) only when
\((i,j)\in\mathcal{E}\)) An adjacency matrix is \(\mathbf{A}\in\{0,1\}^{n\times n}\) such that \(\mathbf{A}_{i,j}=1\) if \((i,j)\in\mathcal{E}\) and \(0\) otherwise. The degree matrix \(\mathbf{D}\) is a diagonal degree matrix with \(i\)-th diagonal entry \(\mathbf{D}_{i,i}=\deg_{i}=\sum_{j}\mathbf{A}_{i,j}\). Graph Laplacian is then defined as \(\mathbf{L}=\mathbf{D}-\mathbf{A}\) and one can verify

\[\operatorname{div}(\nabla\mathbf{f})=\sum_{j:(i,j)\in\mathcal{E}}(f_{j}-f_{i}) =-\mathbf{L}\mathbf{f}.\]

Further, graph gradient can be represented with the incidence matrix \(\mathbf{G}\in\mathbb{R}^{e\times n}\). Specifically \(\mathbf{G}_{k,i}=1\) if edge \(k\) enters node \(i\), \(-1\) if edge \(k\) leaves node \(i\) and \(0\) otherwise. Graph divergence is thus given by \(-\mathbf{G}^{\top}\), which is the negative adjoint of the gradient. The edge direction in \(\mathbf{G}\) can be arbitrarily chosen for undirected graphs because the Laplacian is indifferent to the choice of direction as \(\mathbf{L}=\mathbf{G}^{\top}\mathbf{G}\).

Graph diffusion and graph neural networks.In this work, we consider \(\mathbf{x}:\mathcal{V}\rightarrow\mathbb{R}^{c}\) as a multi-channel signal (or feature) over the node set. We denote \(\mathbf{x}_{i}\in\mathbb{R}^{c}\) as the signal over node \(i\) and let \(\mathbf{X}\in\mathbb{R}^{n\times c}\) collects the signals over all the nodes. The previously defined discrete gradient and divergence operators lead to the following graph diffusion process, which can be seen as a discrete version of heat diffusion equation in (1):

\[\frac{\partial\mathbf{x}_{i}}{\partial t}=\operatorname{div}(D\nabla\mathbf{X })_{i}=\sum_{j:(i,j)\in\mathcal{E}}D(\mathbf{x}_{i},\mathbf{x}_{j},t)(\mathbf{ x}_{j}-\mathbf{x}_{i}),\]

where the diffusivity \(D(\mathbf{x}_{i},\mathbf{x}_{j},t)\) is often scalar-valued and positive, which applies channel-wise. In the special case of homogeneous diffusion, i.e., \(D(\mathbf{x}_{i},\mathbf{x}_{j},t)=1\), the diffusion process can be written as \(\frac{\partial\mathbf{x}_{i}}{\partial t}=\operatorname{div}(\nabla\mathbf{X })_{i}=-(\mathbf{L}\mathbf{X})_{i}\). This process is known as the Laplacian smoothing where the signals become progressively smooth by within-neighbourhood averaging. The solution to the diffusion equation is given by the heat kernel, i.e., \(\mathbf{X}(t)=\exp(-t\mathbf{L})\mathbf{X}(0)\). In fact, the graph heat equation can be also derived as the gradient flow of the so-called Dirichlet energy, which measures the local variations: \(\mathcal{E}_{\operatorname{dir}}(\mathbf{X})=\sum_{(i,j)\in\mathcal{E}}\| \mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}=\operatorname{tr}(\mathbf{X}^{\top} \mathbf{L}\mathbf{X})\).

The graph diffusion process is related to the famous graph convolutional network (GCN) [58], where the latter can be viewed as taking the Euler discretization of the former with a unit stepsize. In addition, GCN defines diffusion with (symmetrically) normalized adjacency \(\widehat{\mathbf{A}}=\mathbf{D}^{-1/2}\mathbf{A}\mathbf{D}^{-1/2}\) and Laplacian \(\widehat{\mathbf{L}}=\mathbf{I}-\widehat{\mathbf{A}}\). The use of normalized Laplacian \(\widehat{\mathbf{L}}\) in place of the combinatorial Laplacian \(\mathbf{L}\) switches the dynamics from being homogeneous to isotropic (where diffusion is weighted by node degrees). More precisely, the discretized dynamics gives the update \(\mathbf{X}^{\ell+1}=\mathbf{X}^{\ell}-\widehat{\mathbf{L}}\mathbf{X}^{\ell}= \widehat{\mathbf{A}}\mathbf{X}^{\ell}\). It can be readily verified that GCN corresponds to the gradient flow of a normalized Dirichlet energy \(\mathcal{E}_{\operatorname{dir}}(\mathbf{X})=\sum_{(i,j)\in\mathcal{E}}\| \mathbf{x}_{i}/\sqrt{\deg_{i}}-\mathbf{x}_{j}/\sqrt{\deg_{j}}\|^{2}=\operatorname {tr}(\mathbf{X}^{\top}\widehat{\mathbf{L}}\mathbf{X})\). Additional channel mixing \(\mathbf{W}\) and nonlinear activation \(\sigma(\cdot)\) are added, leading to a single GCN layer, \(\mathbf{X}^{\ell+1}=\sigma\big{(}\widehat{\mathbf{A}}\mathbf{X}^{\ell}\mathbf{ W}^{\ell}\big{)}\).

## 3 A general framework of continuous dynamics informed GNNs beyond diffusion

Diffusion dynamics has been shown to underpin the design of message passing and graph convolutional networks. While showing promise in many applications, vanilla diffusion dynamics may suffer from performance degradation due to overly smoothing graph signals (oversmoothing), message passing bottlenecks (oversquashing) and graph heterophily. This has motivated the consideration of alternative dynamics other than isotropic diffusion, including anisotropic diffusion, diffusion with source term, geometric diffusion,
oscillation, convection, advection and reaction. Many of them are directly adapted from the existing physical processes.

This section summarizes existing developments on graph neural dynamics under a general framework and categories the literature by the underlying continuous system. We follow the same notations in Section 2 where a graph is represented as \(\mathcal{G}=(\mathcal{V},\mathcal{E})\) with \(|\mathcal{V}|=n\) and \(|\mathcal{E}|=e\). The graph signals or features are encoded in a matrix \(\mathbf{X}\), which is in general time-dependent. Unless mentioned otherwise, we omit the time dependence and treat \(\mathbf{X}\) as \(\mathbf{X}(t)\) for notation clarity. We also denote \(\mathbf{X}^{0}=\mathbf{X}(0)\) as the initial conditions for the system.

The general framework for designing continuous graph dynamics is given as follows, which relates the time derivatives of the signals with spatial derivatives on graphs.

\[\Big{[}\frac{\partial\mathbf{X}}{\partial t},\frac{\partial^{2}\mathbf{X}}{ \partial t^{2}}\Big{]}=\mathcal{F}_{\mathcal{G}}(\mathbf{X},\nabla\mathbf{X}), \tag{2}\]

where \(\mathcal{F}_{\mathcal{G}}\) is a spatial coupling function that is usually parameterized by neural networks. The initial condition of the system is generally the input graph signals (after an encoder). We have summarized and compared existing works discussed in the survey in Table 1 in terms of the driving mechanisms and problems addressed. As we shall see, most of the existing works consider the first-order time derivative while some works explore the second-order time derivative to encode oscillatory systems.

### Anisotropic diffusion

The first class of dynamics, _anisotropic diffusion_, generalizes the isotropic diffusion in GCN, offering great flexibility in controlling the local diffusivity patterns. In image processing, anisotropic diffusion has been extensively applied for low-level tasks, such as image denoising, restoration and segmentation [96]. In particular, the Perona-Malik model [73] sets the diffusivity coefficient \(D\propto|\nabla x|^{-1}\), which is often called the edge indicator as it preserves the sharpness of signals by slowing down diffusion in regions of high variations.

The idea of using anisotropic diffusion to define continuous GNN dynamics has been firstly considered by [75] and formalized by [15], where a class of graph neural diffusion (GRAND) dynamics is proposed. The anisotropic diffusion dynamics is formally given by

\[\text{GRAND}:\quad\frac{\partial\mathbf{X}}{\partial t}=\mathrm{div}\big{(} \mathbf{A}(\mathbf{X})\cdot\nabla\mathbf{X}\big{)},\]

where \(\mathbf{A}(\mathbf{X})\in\mathbb{R}^{n\times n}\) encodes the anisotropic diffusivity along the edges. Here \(\nabla\mathbf{X}\in\mathbb{R}^{n\times n\times c}\) collects the gradient along edges and we use \(\cdot\) to represent the elementwise multiplication of diffusivity coefficients broadcast across the channel dimension. The coefficient is determined by the feature similarity, i.e., \(\mathbf{A}(\mathbf{X})=[a(\mathbf{x}_{i},\mathbf{x}_{j})]_{(i,j)\in\mathcal{E}}\) where \(a(\mathbf{x}_{i},\mathbf{x}_{j})=(\mathbf{W}_{K}\mathbf{x}_{i})^{\top}\mathbf{ W}_{Q}\mathbf{x}_{j}\) computes the dot product attention with learnable parameters \(\mathbf{W}_{K},\mathbf{W}_{Q}\). Softmax normalization is performed on \(\mathbf{A}(\mathbf{X})\) to ensure it is right-stochastic (row sums up to one). Notably, the explicit-Euler discretization of GRAND corresponds to the propagation of graph attention network (GAT) [91]. Several versions of GRAND are proposed to render the dynamics more adaptable compared to the discretized version in [91]. In particular, the diffusivity matrix \(\mathbf{A}(\mathbf{X})\) can be fixed as \(\mathbf{A}(\mathbf{X}^{0})\) that only depends on the initial features. This leads to a GAT propagation with shared attention weights. In addition, \(\mathbf{A}(\mathbf{X})\) can vary according to a dynamically rewired edge set based on the attention score, i.e., \(\mathcal{E}\leftarrow\{(i,j):(i,j)\in\mathcal{E},a(\mathbf{x}_{i},\mathbf{x}_ {j})>\rho\}\) for some threshold \(\rho>0\). Instructively, this interprets the graph structure as discrete realization of certain underlying domains where graph rewiring dynamically changes the spatial discretization.
Building on the formulation of GRAND, BLEND further [16] augments the input signals \(\mathbf{x}_{i}\) with position

\begin{table}
\begin{tabular}{l|l l|l l l} \hline \hline  & \multirow{2}{*}{Methods} & \multirow{2}{*}{Mechanism} & \multirow{2}{*}{OSM} & \multicolumn{3}{c}{Problems Tackled} \\  & & & & OSM & OSQ & HETERO & STAB \\ \hline \multirow{8}{*}{Anisotropic diffusion} & GRAND [15] & Attention diffusivity \& rewiring & & \(\boldsymbol{\bigvee}^{*}\) & \(\boldsymbol{\bigvee}\) \\  & BLEND [16] & Position augmentation & & \(\boldsymbol{\bigvee}^{*}\) & & \(\boldsymbol{\bigvee}\) \\  & Mean Curvature [83] & \multirow{2}{*}{Non-smooth edge indicators} & & & \(\boldsymbol{\bigvee}\) \\  & Beltrami [83] & & & & \\  & \(p\)-Laplacian [35] & \(p\)-Laplacian regularization & \(\boldsymbol{\bigvee}^{\dagger}\) & & \\  & DIFFormer [98] & Full graph transformer & & \(\boldsymbol{\bigvee}\) & & \\  & DIGNN [34] & Parameterized Laplacian & \(\boldsymbol{\bigvee}^{\dagger}\) & & \\  & GRAND++ [86] & Source term & & & \\  & PDE-GCN\({}_{\text{D}}\)[29] & Nonlinear diffusion & & \(\boldsymbol{\bigvee}\) & & \\  & GIND [18] & Implicit nonlinear diffusion & & \(\boldsymbol{\bigvee}\) & \(\boldsymbol{\bigvee}\) \\ \hline \multirow{2}{*}{Oscillation} & PDE-GCN\({}_{\text{H}}\)[29] & Wave equation & & \(\boldsymbol{\bigvee}\) & \\  & GraphCON [78] & Damped coupled oscillation & & \(\boldsymbol{\bigvee}\) & \(\boldsymbol{\bigvee}\) \\ \hline \multirow{2}{*}{Non-local dynamics} & FLODE [68] & Fractional Laplacian & \(\boldsymbol{\bigvee}\) & \(\boldsymbol{\bigvee}\) & \(\boldsymbol{\bigvee}\) \\  & QDC [67] & Quantum diffusion & & \(\boldsymbol{\bigvee}\) & \(\boldsymbol{\bigvee}\) \\  & TIDE [7] & Learnable heat kernel & & \(\boldsymbol{\bigvee}\) & \\  & G2TN [89] & Hypo-elliptic diffusion & & \(\boldsymbol{\bigvee}\) & \\ \hline \multirow{8}{*}{Diffusion with external forces} & CDE [103] & Convection diffusion & & \(\boldsymbol{\bigvee}\) & \\  & GREAD [21] & Reaction diffusion & & \(\boldsymbol{\bigvee}\) & \\  & ACMP [94] & Allen-Chan retraction & & \(\boldsymbol{\bigvee}\) & \\  & & with negative diffusivity & & & \\  & ODNet [66] & Diffusion with confidence & & \(\boldsymbol{\bigvee}\) & \(\boldsymbol{\bigvee}\) \\  & ADR-GNN [30] & Advection reaction diffusion & & \(\boldsymbol{\bigvee}\) & \\  & G\({}^{2}\)[79] & Diffusion gating & & \(\boldsymbol{\bigvee}\) & \(\boldsymbol{\bigvee}\) \\  & MHKG [81] & Reverse diffusion & & \(\boldsymbol{\bigvee}\) & \(\boldsymbol{\bigvee}\) \\  & A-DGN [41] & Anti-symmetric weight & & \(\boldsymbol{\bigvee}\) & \(\boldsymbol{\bigvee}\) \\ \hline \multirow{8}{*}{Geometry underpinned} & NSD [10] & Sheaf diffusion & & \(\boldsymbol{\bigvee}\) & \\  & Hamiltonian\({}_{G}\), etc. [42] & Bracket dynamics with & & \\  & HamGNN [57] & higher-order cliques & & & \\  & HamGNN [102] & Learnable Hamiltonian dynamics & & & \(\boldsymbol{\bigvee}\) \\ \hline Gradient flow & GRAFF [39] & Parameterized gradient flow & & \(\boldsymbol{\bigvee}\) & \\ \hline Multi-scale diffusion & GradFUFG [46] & 
\begin{tabular}{l} Separated diffusion for low-pass \\ and high-pass at different scales \\ \end{tabular} & & \(\boldsymbol{\bigvee}\) & \\ \hline \hline \end{tabular}

* The ability of the methods to mitigate oversquashing is through graph rewiring. \({}^{\dagger}p\)-Laplacian and DIGNN avoids oversmoothing provided the input dependent regularization (i.e., a source term) is added.

\end{table}
Table 1: Summary of continuous dynamics informed graph neural networks, including the driving mechanism and the problems addressed, including oversmoothing (OSM), oversquashing (OSQ), graph heterophily (HETERO), and stability (STAB) with respect to perturbation or with training, such as gradient vanishing or explosion associated with increased depth. Note we show the problems addressed either theoretically or empirically in the paper (unless proven otherwise in subsequent literature). More detailed discussions are in Section 5.
coordinates \(\mathbf{u}_{i}\) for each node. The diffusion process then becomes

\[\text{BLEND}:\quad\frac{\partial[\mathbf{X},\mathbf{U}]}{\partial t}=\mathrm{div} (\mathbf{A}([\mathbf{X},\mathbf{U}])\cdot\nabla[\mathbf{X},\mathbf{U}]).\]

The joint diffusion over both positions and signals is motivated by diffusion on Riemannian manifolds with the Laplace-Beltrami operator, in which the gradient, diffusivity and divergence all depend on the Riemannian metric (that varies according to the position). In a similar vein, augmenting the position information while diffusing on graphs produces joint evolution of features as well as topology, which further allows graph rewiring to improve the information flow. In particular, based on the evolved positional information, the graph is dynamically rewired as \(\mathcal{E}\leftarrow\{(i,j):d_{\mathcal{C}}(\mathbf{u}_{i},\mathbf{u}_{j})<r\}\) or by k-nearest neighbour graph. The positional information can be pre-computed by personalized PageRank [36], deepwalk [74] or even learned hyperbolic embeddings [17].

In BLEND, the diffusivity is given by the attention score over both the positional features and signals, which is in fact the core idea of transformers [90] that augments samples with positional encoding. [98] develop a transformer-based diffusion process (called DIFFormer) where attention diffusivity coefficients are computed on a fully connected graph \(\mathcal{V}\times\mathcal{V}\). The input graph (represented via input adjacency matrix \(\mathbf{A}^{0}\)) serves as a geometric prior that augments the learned attention matrix.

\[\text{DIFFormer}:\quad\frac{\partial\mathbf{X}}{\partial t}=\mathrm{div}_{ \mathcal{V}\times\mathcal{V}}\big{(}(\mathbf{A}^{0}+\mathbf{A}(\mathbf{X})) \cdot\nabla\mathbf{X}\big{)},\]

where \(\mathrm{div}_{\mathcal{V}\times\mathcal{V}}(\mathbf{F})_{i}=\sum_{j\in \mathcal{V}}\mathbf{F}_{i,j}\) denotes the divergence operator on the complete graph and \(\mathbf{A}(\mathbf{X})\) computes the diffusivity by a transformer block [90], different to the ones used by GRAND and BLEND. A recent work [97] has shown the evolution via both local message passing through \(\mathbf{A}^{0}\) and global attention diffusion through \(\mathbf{A}(\mathbf{X})\) improves the generalization under topological distribution shift, i.e., when training and test graph topology differs.

Apart from its connection to transformer-based dynamics, BLEND is in fact motivated by a geometric evolution known as Beltrami flow where the diffusion over a non-Euclidean domains also depends on the varying metric. In [83], the Beltrami flow, along with mean-curvature flow, is generalized to graph domains by explicitly factorizing out the edge indicator \(\|\delta\mathbf{x}_{i}\|:=\sqrt{\sum_{j:(i,j)\in\mathcal{E}}\|\mathbf{x}_{j}- \mathbf{x}_{i}\|^{2}}\). Formally, let \(\mathbf{S}_{\mathrm{mc}}(\mathbf{X})\), \(\mathbf{S}_{\mathrm{bel}}(\mathbf{X})\) be the diffusivity coefficients of mean curvature and Beltrami diffusion, with their elements defined as \([\mathbf{S}_{\mathrm{mc}}(\mathbf{X})]_{i,j}=\frac{1}{\|\delta\mathbf{x}_{i} \|}+\frac{1}{\|\delta\mathbf{x}_{j}\|}\) and \([\mathbf{S}_{\mathrm{bel}}(\mathbf{X})]_{i,j}=\frac{1}{\|\delta\mathbf{x}_{i} \|^{2}}+\frac{1}{\|\delta\mathbf{x}_{i}\|\|\delta\mathbf{x}_{j}\|}\). The diffusion processes [83] propose are

\[\text{Mean Curvature}:\quad\frac{\partial\mathbf{X}}{\partial t} =\mathrm{div}\Big{(}\big{(}\mathbf{A}(\mathbf{X})\odot\mathbf{S}_{ \mathrm{mc}}(\mathbf{X})\big{)}\cdot\nabla\mathbf{X}\Big{)},\] \[\text{Beltrami}:\quad\frac{\partial\mathbf{X}}{\partial t} =\mathrm{div}\Big{(}\big{(}\mathbf{A}(\mathbf{X})\odot\mathbf{S} _{\mathrm{bel}}(\mathbf{X})\big{)}\cdot\nabla\mathbf{X}\Big{)},\]

where \(\mathbf{A}(\mathbf{X})\) is computed from the dot product attention following the previous works [15, 16]. For both dynamics, non-smooth signals are preserved by slowing down diffusion where signal abruptly changes. As commented by the paper, positional information can be added in a similar way as BLEND [16]. It should be noticed that although BLEND originates from the Beltrami flow, the dynamics turns out to be the same as GRAND, augmented with positional embeddings.

A more general \(p\)-Laplacian based graph neural network is proposed by [35] where the diffusion is derived from a \(p\)-Laplacian regularization framework.

\[p\text{-Laplacian}:\quad\frac{\partial\mathbf{X}}{\partial t}=\mathrm{div} \big{(}\|\nabla\mathbf{X}\|^{p-2}\cdot\nabla\mathbf{X}\big{)}-\mu(\mathbf{X}- \mathbf{S}),\]
with \(\mathbf{S}\) as a source term and \(\mu>0\) controlling the regularization strength. The diffusivity \(\|\nabla\mathbf{X}\|^{p-2}\) is an \(n\times n\) matrix with elements \([\|\nabla\mathbf{X}\|^{p-2}]_{i,j}=\|\mathbf{x}_{j}-\mathbf{x}_{i}\|^{p-2}\) if \((i,j)\in\mathcal{E}.\) The injection of source information can be understood physically as the heat exchange from the system to the outside. In the paper the source term is simply the input feature matrix, i.e., \(\mathbf{S}=\mathbf{X}^{0}\). When \(p=2\), the diffusion reduces to the heat diffusion with classic Laplacian. When \(p=1\), the dynamics recovers the mean curvature flow (although the definition slightly differs from the one in [83]). As a result, with properly selected \(p\), the dynamics is flexible in adapting to different types of graphs and able to perverse the boundaries without oversmoothing the signals. It is worth mentioning that unlike previous works that use discretization to update \(\mathbf{X}\), the paper directly solves for the equilibrium state by setting \(\frac{\partial\mathbf{X}}{\partial t}=0\), which leads to an implicit graph diffusion layer given by \(p\)-Laplacian message passing.

The \(p\)-Laplacian diffusion corresponds to the gradient flow of a \(p\)-Dirichlet energy given with a regularization term \(\|\mathbf{X}-\mathbf{S}\|^{2}\). A similar idea has been considered in [22] where the Dirichlet energy is replaced with the total variation on graph gradients, which leverages the \(L_{1}\) norm (which is different to the case of \(p=1\) in the \(p\)-Laplacian diffusion [83]). A dual-optimization scheme is introduced due to the non-differentiablity of the objective at zero.

A recent paper [34] parameterizes the graph gradient and divergence instead of the diffusivity coefficients as in previous works, and defines a parameterized graph Laplacian for diffusion process. In particular, [34] consider weighted inner products for both the \(L^{2}(\mathcal{V})\) and \(L^{2}(\mathcal{E})\), i.e., \(\langle f,g\rangle_{L^{2}(\mathcal{V})}=\sum_{i\in\mathcal{V}}\chi_{i}f_{i}g_ {i}\) and \(\langle F,G\rangle_{L^{2}(\mathcal{E})}=\phi_{i,j}F_{i,j}G_{i,j}.\) The graph gradient is defined as \((\nabla_{\Theta}f)_{i,j}\coloneqq\psi_{i,j}(f_{j}-f_{i})\), which also leads to a notion of graph divergence \((\mathrm{div}_{\Theta}F)_{i}=\frac{1}{2\chi_{i}}\sum_{j\in\mathcal{N}_{i}} \psi_{i,j}\phi_{i,j}(F_{i,j}-F_{j,i})\), where \(\chi_{i},\phi_{i,j},\psi_{i,j}\) are strictly positive real-valued functions on nodes and edges respectively. Here we denote \(\nabla_{\Theta},\mathrm{div}_{\Theta}\) to emphasize that the gradient and divergence operators are parameterized. Because the graph gradient is parameterized and may not be anti-symmetric, i.e., \((\nabla_{\Theta}f)_{i,j}\neq-(\nabla_{\Theta}f)_{j,i}\), the divergence encodes directional information. The paper also parameterizes the weighting functions to be node-dependent, i.e., \(\chi_{i}=\chi_{i}(\mathbf{x}_{i}),\phi_{i,j}=\phi_{i,j}(\mathbf{x}_{i},\mathbf{ x}_{j})\) and \(\psi_{i,j}=\psi_{i,j}(\mathbf{x}_{i},\mathbf{x}_{j})\), which involves learnable parameters. The diffusion process the paper considers is thus

\[\text{DIGNN}:\quad\frac{\partial\mathbf{X}}{\partial t}=\mathrm{div}_{\Theta} (\nabla_{\Theta}\mathbf{X})-\mu(\mathbf{X}-\mathbf{S}),\]

where a regularization term \(\|\mathbf{X}-\mathbf{S}\|^{2}\) is added similarly in [35].

The idea of adding an energy source has also been considered in GRAND++ [86] based on the framework of anisotropic diffusion of GRAND:

\[\text{GRAND++}:\quad\frac{\partial\mathbf{X}}{\partial t}=\mathrm{div}( \mathbf{A}(\mathbf{X})\cdot\nabla\mathbf{X})+\mathbf{S}\]

where \(\mathbf{S}\in\mathbb{R}^{n\times c}\) is a source term. The paper proposes a random walk viewpoint to show that, without the source term, the dynamics reduces to GRAND and is guaranteed to converge to a stationary distribution independent of the initial conditions \(\mathbf{X}^{0}\). Each row of the source term \(\mathbf{S}\) is defined as \(\mathbf{s}_{i}=\sum_{j\in\mathcal{I}}\delta_{ij}(\mathbf{x}_{j}^{0}-\bar{ \mathbf{x}}^{0})\) with \(\mathcal{I}\subseteq\mathcal{V}\) a selected node subset used as the source term and \(\bar{\mathbf{x}}^{0}=\frac{1}{|\mathcal{I}|}\sum_{j\in\mathcal{I}}\mathbf{x}_{j }^{0}\) is the average signal. \(\delta_{ij}\) denotes the initial transition probability from node \(j\) to \(i\). By construction, the limiting signal distribution is close to an interpolation of the source signals in the selected subset \(\mathcal{I}\). Similar idea of source term injection has also been considered in earlier work [99], which can be seen as the homogeneous diffusion with both a source term and a residual term. The proposed model (called CGNN) follows the dynamics \(\frac{\partial\mathbf{X}}{\partial t}=-\mathbf{L}\mathbf{X}+\mathbf{X}\mathbf{ W}+\mathbf{X}^{0}\), for some learnable channel mixing matrix \(\mathbf{W}\).

Anisotropic diffusion is generally nonlinear in the sense that diffusivity depends nonlinearly on the mass along the evolution. This is the case in edge-preserving dynamics as the diffusivity explicitly depends
nonlinearly on the gradient. GRAND-based dynamics is also nonlinear as long as the attention coefficients is computed for each timestep. In [29, 18], additional nonlinearity is further incorporated by factoring the diffusivity \(D\) as the composition of a linear operator \(\mathcal{K}\) and its adjoint \(\mathcal{K}^{*}\), i.e., \(D=\mathcal{K}^{*}\mathcal{K}\). Then a pointwise nonlinearity function \(\sigma(\cdot)\) is added, leading to

\[\frac{\partial\mathbf{X}}{\partial t}=\mathrm{div}\big{(}\mathcal{K}^{*} \sigma(\mathcal{K}\nabla\mathbf{X})\big{)}.\]

In the case when \(\sigma\) is the identity map, the dynamics recovers the anisotropic diffusion. Such a nonlinear system has been firstly proposed by [80] for defining convolutional residual networks for images. In [29], the idea is generalized to graphs by defining \(\mathcal{K}\) as learnable pointwise convolution. Specifically, the dynamics, called PDE-\(\text{GCN}_{\text{D}}\), can be written in terms of the gradient and divergence operator as follows.

\[\text{PDE-GCN}_{\text{D}}:\quad\frac{\partial\mathbf{X}}{\partial t}=-\mathbf{ G}^{\top}\mathbf{K}^{\top}\sigma(\mathbf{K}\mathbf{G}\mathbf{X}),\]

where \(\mathbf{K}\) is a learnable parameter and \(\mathbf{G}\) is the gradient operator defined in Section 2. It can be readily observed that when \(\sigma\) is identity and \(\mathbf{K}=\mathbf{I}\), the dynamics reduces to the heat diffusion implemented by GCN.

In the follow-up work [18], the linear operator \(\mathcal{K}\) (parameterized by \(\mathbf{K}\)) is applied over the channel space instead of the edge space, i.e.,

\[\text{GIND}:\quad\frac{\partial\mathbf{X}}{\partial t}=-\mathbf{G}^{\top} \sigma(\mathbf{G}\mathbf{X}\mathbf{K}^{\top})\mathbf{K}.\]

Motivated by [43, 18] consider an implicit propagation of GIND as \(\mathbf{Z}=-\mathbf{G}^{\top}\sigma\big{(}\mathbf{G}(\mathbf{Z}+b_{\mathbf{ \Omega}}(\mathbf{X}^{0}))\mathbf{K}^{\top}\big{)}\mathbf{K}\), where \(b_{\mathbf{\Omega}}(\cdot)\) is an affine transformation with parameter \(\mathbf{\Omega}\). The model corresponds to a refinement process for the flux \(\mathbf{Z}\) and the output is given by a decoder over \(\mathbf{Z}+\mathbf{X}^{0}\). It has been shown the equilibrium state of the implicit diffusion corresponds to the minimizer of a convex objective function provided the nonlinearity is monotone and Lipschitz and \(\mathbf{K}\otimes\mathbf{G}\) is upper bounded in norm. This result guarantees the convergence of the dynamics and allows structural constraints to be embedded to the dynamics by explicitly modifying the objective.

### Oscillations

The phenomenon of oscillation has been widely found in physics, which primarily features the repetitive motion. Unlike diffusion that dissipates energy, oscillatory system often preserves energy and is thus reversible. Oscillatory processes are often modelled as second-order ordinary/partial differential equations. One simple example of oscillatory process is characterized by the wave equation \(\frac{\partial^{2}x}{\partial t^{2}}=c\Delta x\), which is a hyperbolic PDE. The wave equation has been considered in [29] for defining dynamics on graphs, which follows the nonlinear formalism in [80]:

\[\text{PDE-GCN}_{\text{H}}:\quad\frac{\partial^{2}\mathbf{X}}{\partial t^{2}} =\mathrm{div}\big{(}\mathcal{K}^{*}\sigma(\mathcal{K}\nabla\mathbf{X})\big{)}= -\mathbf{G}^{\top}\mathbf{K}^{\top}\sigma(\mathbf{K}\mathbf{G}\mathbf{X}).\]

In addition, GraphCON [78] considers a more general oscillatory dynamics which combines a damped oscillating process with a coupling function. That is,

\[\text{GraphCON}:\quad\frac{\partial^{2}\mathbf{X}}{\partial t^{2}}=\sigma(F_{ \mathcal{G}}(\mathbf{X}))-\gamma\mathbf{X}-\alpha\frac{\partial\mathbf{X}}{ \partial t},\]
where \(F_{\mathcal{G}}(\mathbf{X})_{i}=F_{\mathcal{G}}(\mathbf{x}_{i},\{\mathbf{x}_{j}\}_ {j\in\mathcal{N}_{i}})\) is the coupling function, \(\alpha\geq 0\) and \(\sigma(\cdot)\) is some nonlinear activation. When \(\sigma(F_{\mathcal{G}}(\mathbf{X}))=0\), \(\alpha=0\), the system reduces to the classic harmonic oscillation for each node independently. Adding the damping term \(-\frac{\partial\mathbf{X}}{\partial t}\) mimics the frictional force that diminishes the oscillation. Finally, due to the presence of interdependence between the nodes, the coupling function is required to model the interactions between the nodes. The paper mainly considers two choices of coupling function, with isotropic and anisotropic diffusion (which leads to GCN and GAT respectively). Formally, \(F_{\mathcal{G}}(\mathbf{X})_{i}=\sum_{j:(i,j)\in\mathcal{E}}\mathbf{A}(\mathbf{ X})_{i,j}\mathbf{x}_{j}\) represents the message passing with normalized adjacency or learned attention scores. When \(\gamma=1\), \(\sigma\) is identity, the GraphCON can be rewritten as \(\frac{\partial^{2}\mathbf{X}}{\partial t}=\mathrm{div}(\mathbf{A}(\mathbf{X}) \cdot\nabla\mathbf{X})-\alpha\frac{\partial\mathbf{X}}{\partial t}\), which is effectively the wave equation with a damping term. GraphCON is flexible in that the coupling term can accommodate arbitrary message passing scheme. In addition, it possesses greater expressive power by showing the GNN induced by the coupling function approaches the steady state of GraphCON, while the latter explores the entire trajectory.

### Non-local dynamics

We have so far focused on local dynamics, in the sense that the diffusion or oscillation happens locally within the neighbourhood. Thus it usually requires sufficiently large timestep for one node's influence to reach a distant node (with respect to graph topology). Graph rewiring employed in GRAND and BLEND can be utilized to enable long-range diffusion by modifying the graph topology. This section explores various dynamics-based formulations that transcend the local community when propagating information, resulting in non-localized dynamics.

Fractional Laplacian.Fractional Laplacian [76, 64] has been effective to represent complex anomalous processes, such as fluids dynamics in porous medium [14]. A recent work [68] utilizes the fractional graph Laplacian/adjacency matrix \(-\widehat{\mathbf{A}}^{\alpha}\) (for some \(\alpha\in\mathbb{R}\)) in order to define a non-local diffusion process as

\[\text{FLODE}:\quad\frac{\partial\mathbf{X}}{\partial t}=-\widehat{\mathbf{A}} ^{\alpha}\mathbf{X}.\]

where \(\widehat{\mathbf{A}}\) is the symmetric normalized adjacency matrix. A critical difference compared to \(p\)-Laplacian in terms of the order is that here \(\alpha\) can be fractional, instead of being restricted to integers. The fractional Laplacian is often dense, and thus the corresponding diffusion is non-local where long-range interactions are captured. When coupled with a symmetric channel mixing matrix \(\mathbf{W}\), i.e., \(\frac{\partial\mathbf{X}}{\partial t}=-\widehat{\mathbf{A}}^{\alpha}\mathbf{ X}\mathbf{W}\), the flexibility in the choice of \(\alpha\) allows dynamics to accommodate both smoothing and sharpening effects, which avoids oversmoothing and is suited for heterophilic graphs. The work also extends the formulation to directed graphs and correspondingly defines the notions of oversmoothing and Dirichlet energy with the asymmetric Laplacian. On directed graphs, the fractional Laplacian is defined through singular value decomposition, i.e., \(\widehat{\mathbf{A}}^{\alpha}\coloneqq\mathbf{U}\Sigma^{\alpha}\mathbf{V}^{H}\), where \(\mathbf{U},\mathbf{V}\in\mathbb{C}^{n\times n}\) are unitary singular vectors and \(\mathbf{\Sigma}\in\mathbb{R}^{n\times n}\) contains singular values on the diagonal. Finally, the paper also considers a Schrodinger equation based diffusion as \(\frac{\partial\mathbf{X}}{\partial t}=\mathrm{i}\widehat{\mathbf{A}}^{\alpha} \mathbf{X}\), where \(\mathrm{i}=\sqrt{-1}\) represents the imaginary unit.

Quantum diffusion kernel.The Schrodinger's equation \(\mathrm{i}\frac{\partial\psi(u,t)}{\partial t}=\mathcal{H}\psi(u,t)\) has also been considered in [67], where \(\mathcal{H}\) is the Hamiltonian operator (composed of a kinetic and potential energy operator). \(\psi(u,t)\) denotes the (complex-valued) quantum wave function at position \(u\) at time \(t\), and the square modulus of the wave function \(|\psi(u,t)|^{2}\) indicates the probability density of a particle. With the quantum system defined by \(\psi(u,t)\), a quantum state \(|\psi(t)\rangle\) refers to the superposition of all the states, i.e., \(|\psi(t)\rangle=\int\psi(u,t)|u\rangle du\)
where \(|u\rangle\) is the position state (a basis vector for representing the quantum state). One can recover the components of the quantum state by the inner product between quantum state vector \(|\psi(t)\rangle\) and position vector \(|u\rangle\), i.e., \(\psi(u,t)=\langle u|\psi(t)\rangle\). In a simple quantum system where the Hamiltonian \(\mathcal{H}\) is time-independent, the solution to the Schrodinger's equation can be written in terms of quantum state as \(|\psi(t)\rangle=e^{-\mathrm{i}\mathcal{H}t}|\psi(0)\rangle\).

On a graph with \(n\) nodes, \(|\psi(t)\rangle\in\mathbb{C}^{n}\) can be interpreted as the state vector of particles across all nodes, and the position state refers to the node set. In a system without a potential, the Hamiltonian reduces to the negative Laplacian, and the Schrodinger's equation reduces to \(\mathrm{i}\frac{\partial|\psi(t)\rangle}{\partial t}=-\mathbf{L}|\psi(t)\rangle\). The eigenvectors of \(\mathbf{L}\), denoted as \(\{|\mathbf{\phi}_{i}\rangle\}_{i=1}^{n}\) provides a complete orthogonal basis, which can be treated as the position basis so that one can express the solution as \(|\psi(t)\rangle=\sum_{k=1}^{n}c_{k}e^{\mathrm{i}\lambda_{k}t}|\mathbf{\phi}_{k}\rangle\), where \(\lambda_{k}\) is the \(k\)-th eigenvalue and \(c_{k}=\langle\psi(0)|\mathbf{\phi}_{k}\rangle\). Instead of working with the solution which involves the complex unit, [67] define a kernel that models the the average overlap between any two nodes \(i,j\) as the inner product between \(|\psi(t)\rangle_{i}\) and \(|\psi(t)\rangle_{j}\), which is \(\sum_{k,l=1}^{n}c_{k}^{\star}c_{l}|\mathbf{\phi}_{k}\rangle_{i}^{\ast}|\mathbf{\phi}_{ l}\rangle_{j}\). To engineer observation operators as a spectral filter, the paper further leverages a Gaussian filter \(\mathcal{P}=\sum_{k=1}^{n}e^{-(\lambda_{k}-\mu)^{2}/2\sigma^{2}}\), which leads to the proposed quantum diffusion kernel (QDC) \(\mathbf{Q}\in\mathbb{R}^{n\times n}\), where

\[\text{QDC : }\quad\mathbf{Q}_{i,j}=\sum_{k=1}^{n}e^{-(\lambda_{k}-\mu)^{2}/2 \sigma^{2}}|\mathbf{\phi}_{k}\rangle_{i}^{\ast}|\mathbf{\phi}_{k}\rangle_{j}\]

where the initial quantum state is assumed to be equally delocalized. The kernel matrix \(\mathbf{Q}\) is interpreted as the transition matrix between the nodes and hence can be supplemented for any message-passing-based graph neural networks. Hence the resulting quantum diffusion corresponds to the anisotropic diffusion where the diffusivity is given by the quantum diffusion kernel. The kernel matrix can be further sparsified either using a threshold or KNN. The kernel allows non-local message passing due to the quantum interference across all the position states. Further a multi-scale variant of quantum diffusion is proposed that combines the propagation from quantum diffusion and standard graph diffusion.

Time derivative diffusion.Another work [7] introduces a non-local message passing scheme by combining local message passing with a learnable-timestep heat kernel. Recall the heat diffusion follows \(\frac{\partial\mathbf{X}}{\partial t}=-\mathbf{L}\mathbf{X}\), where its solution is given by the heat kernel as \(\mathbf{X}(t)=\exp(-t\mathbf{L})\mathbf{X}(0)\). One can generalize the heat kernel to capture the transition between any two states in time, i.e., \(\mathbf{X}(t)=\exp(-(t-s)\mathbf{L})\mathbf{X}(s)\). Instead of setting a pre-defined \(t,s\), the paper parameterizes the heat kernel as \(\exp(-t_{\theta}\mathbf{L})\) where \(t_{\theta}\) is the learnable timestep in order to adapt the diffusion range to different types of a dataset. Further, to simultaneously account for local message passing, the paper combines the adjacency matrix \(\mathbf{A}\) with learned heat kernel, which leads to the proposed dynamics as \(\mathbf{X}(t)=\mathbf{A}\exp(-t_{\theta}\mathbf{L})\mathbf{X}(s)\). This corresponds to the following continuous dynamics (up to some normalizing constants):

\[\text{TIDE : }\quad\frac{\partial\mathbf{X}}{\partial t}=\mathrm{div}( \mathbf{A}\exp(-t_{\theta}\mathbf{L})\cdot\nabla\mathbf{X}).\]

It is clear that when \(t_{\theta}=0\), the model reduces to the classic GCN. The learnability of \(t_{\theta}\) ensures the model is flexible to capture both local and multi-hop communication.

Hypo-elliptic diffusion.In [89], non-local and higher-order information is captured via the so-called _hypo-elliptic diffusion_, which is based on a tensor-valued Laplacian that aggregates the entire trajectories of random walks on graphs. The sequential nature of the path can be characterized with the free (associative) algebra, which lifts the sequence injectively to a vector space with non-commutative multiplication (i.e., an
algebra). More formally, an algebra \(H\) over \(\mathbb{R}^{c}\) can be realized as a sequence of tensors with increasing order, i.e., \(H\coloneqq\{\mathbf{v}=(\mathbf{v}^{0},\mathbf{v}^{1},\mathbf{v}^{2},...)\colon \mathbf{v}^{m}\in(\mathbb{R}^{c})^{\otimes m}\}\), where \((\mathbb{R}^{c})^{\otimes m}\) denotes the space of \(m\)-order tensors. For example, \((\mathbb{R}^{c})^{\otimes 0}\equiv\mathbb{R}\), \((\mathbb{R}^{c})^{\otimes 1}\equiv\mathbb{R}^{c}\). The scalar multiplication and vector addition of \(H\) is defined according to \(\lambda\mathbf{v}\coloneqq(\lambda\mathbf{v}^{m})_{m\geq 0}\) for \(\lambda\in\mathbb{R}\) and \(\mathbf{v}+\mathbf{w}\coloneqq(\mathbf{v}^{m}+\mathbf{w}^{m})_{m\geq 0}\). The algebra multiplication of \(H\) is given by \(\mathbf{v}\cdot\mathbf{w}\coloneqq\big{(}\sum_{k=0}^{m}\mathbf{v}^{k}\otimes \mathbf{w}^{m-k}\big{)}_{m\geq 0}\). \(H\) can be further made into a Hilbert space with the chosen inner product \(\langle\mathbf{v},\mathbf{w}\rangle\coloneqq\sum_{m\geq 0}\langle\mathbf{v}^{m}, \mathbf{w}^{m}\rangle^{m}\) where \(\langle\cdot,\cdot\rangle^{m}\) denotes the classic inner product on the tensor space \((\mathbb{R}^{c})^{\otimes m}\).

With the properly defined algebra, one can lift a sequence to such space, thus summarizing the full details of its trajectory. Specifically, denote the space of sequences in \(\mathbb{R}^{c}\) as \(\mathrm{Seq}(\mathbb{R}^{c})\coloneqq\bigcup_{k=0}^{\infty}(\mathbb{R}^{c})^{ k+1}\), where \((\mathbb{R}^{c})^{k+1}\) denotes the \(k+1\)-product space of \(\mathbb{R}^{c}\). A sequence, denoted as \(\mathbf{x}=(\mathbf{x}^{0},\mathbf{x}^{1},...,\mathbf{x}^{k})\in(\mathbb{R}^{ c})^{k+1}\) is an element of the sequence space \(\mathrm{Seq}(\mathbb{R}^{c})\). Let \(\varphi:\mathbb{R}^{c}\to H\) be an algebra lifting, which allows to define a sequence feature map \(\tilde{\varphi}(\mathbf{x})=\varphi(\mathbf{x}^{0})\cdot\varphi(\mathbf{x}^{1 }-\mathbf{x}^{0})\cdots\varphi(\mathbf{x}^{k}-\mathbf{x}^{k-1})\in H\). One example of such injective map is the tensor exponential given by \(\varphi(\mathbf{u})=\exp_{\otimes}(\mathbf{u})\coloneqq\big{(}\frac{\mathbf{ u}^{\otimes m}}{m!}\big{)}_{m\geq 0}\) for \(\mathbf{u}\in\mathbb{R}^{c}\). Such a feature map is able to summarize the entire sequence path up to step \(k\).

On a graph, let \(\mathbf{x}_{i}^{k}\in\mathbb{R}^{c}\) denote the signals at node \(i\in\mathcal{V}\) at diffusion step \(k\geq 0\). Instead of focusing on the signals at a particular timestep \(k\), the paper leverages the sequence map to capture the entire past trajectory of the diffusion process through \(\tilde{\varphi}(\mathbf{x}_{i})\in H\) where \(\mathbf{x}_{i}\coloneqq(\mathbf{x}_{i}^{0},\mathbf{x}_{i}^{1},...,\mathbf{x}_{ i}^{k})\). The corresponding diffusion process requires a tensor adjacency matrix, \(\widetilde{\mathbf{A}}\in H^{n\times n}\) with the entries \(\widetilde{\mathbf{A}}_{i,j}=\varphi(\mathbf{x}_{j}^{0}-\mathbf{x}_{i}^{0})\in H\) if \((i,j)\in\mathcal{E}\) and \(0\) otherwise. The associated Laplacian \(\widetilde{\mathbf{L}}\) can be defined accordingly. For example, the random walk Laplacian has entries \(\widetilde{\mathbf{L}}_{i,i}=1\) and \(\widetilde{\mathbf{L}}_{i,j}=-\frac{1}{\deg_{i}}\varphi(\mathbf{x}_{j}^{0}- \mathbf{x}_{i}^{0})\) if \((i,j)\in\mathcal{E}\) and \(0\) otherwise. The classic graph heat diffusion is generalized to hypo-elliptic graph diffusion as

\[\text{G2TN}:\quad\frac{\partial\tilde{\varphi}(\mathbf{X})}{\partial t}=- \widetilde{\mathbf{L}}\tilde{\varphi}(\mathbf{X}),\]

where we let \(\tilde{\varphi}(\mathbf{X})\coloneqq[\tilde{\varphi}(\mathbf{x}_{1}),...,\tilde {\varphi}(\mathbf{x}_{n})]\in H^{n}\). The multiplication of \(\widetilde{\mathbf{L}}\) and \(\tilde{\varphi}(\mathbf{X})\) is defined over the space of algebra \(H\) (similarly as how classic matrix multiplication works). In the case of random walk Laplacian, the paper verifies that the solution of the hypo-elliptic diffusion summarizes the entire random walk histories of each node, in contrast to the snapshot state at each timestep given by the classic diffusion equation. Notice here instead of working with node signals \(\mathbf{x}\) directly, the diffusion concerns all the node signals along the trajectory, i.e., \(\varphi(\mathbf{x})\). Thus \(\tilde{\varphi}(\mathbf{x})\) is more expressive compared to \(\mathbf{x}\). The paper further adapts the attention mechanism to define a weighted hypo-elliptic adjacency as \(\widetilde{\mathbf{A}}_{i,j}=a(\mathbf{x}_{i}^{0},\mathbf{x}_{j}^{0})\varphi( \mathbf{x}_{j}^{0}-\mathbf{x}_{i}^{0})\), which correspondingly defines an anisotropic hypo-elliptic diffusion.

### Diffusion with external forces

Most of the aforementioned dynamics is only controlled by a single mechanism, either diffusion or oscillation. This section discusses systems that impart external forces to the diffusion dynamics, such as _convection_, _advection_, and _reaction_. In particular, convection and advection are widely known in physical sciences that describe how the mass transports as a result of the movement of the underlying fluid. Such a process is characterized by \(\frac{\partial x}{\partial t}=-\mathrm{div}(\mathbf{v}x)\) where \(\mathbf{v}\) represents the velocity field of the fluid motion. Reaction process is more general and often found in chemistry where chemical substance interacts with each other and leads to a change of mass and substance, i.e., \(\frac{\partial x}{\partial t}=r(x)\), for some reaction function \(r(\cdot)\). Other mechanisms such as gating, reverse diffusion and anti-symmetry have also been exploited in literature to modulate and control the diffusion dynamics.
Convetion-diffusion.Convection-diffusion dynamics combines the convection with diffusion process in which mass not only transports but also disperse in space. [103] generalize the convection-diffusion equation (CDE) to graphs as

\[\text{CDE}:\quad\frac{\partial\mathbf{X}}{\partial t}=\operatorname{div}( \mathbf{A}(\mathbf{X})\cdot\nabla\mathbf{X})+\operatorname{div}(\mathbf{V} \circ\mathbf{X})\]

where \(\mathbf{V}\) denotes a velocity field and \((\mathbf{V}\circ\mathbf{X})_{i,j}:=\mathbf{V}_{i,j}\odot\mathbf{x}_{j}\) with \(\odot\) representing the elementwise product. In particular, [103] define \(\mathbf{V}_{i,j}=\sigma(\mathbf{W}(\mathbf{x}_{j}-\mathbf{x}_{i}))\in\mathbb{ R}^{c}\) for some nonlinear activation \(\sigma(\cdot)\) and learnable matrix \(\mathbf{W}\). Such a choice is motivated from the heterophilic graphs where neighbouring nodes exhibit diverse features. Hence \(\mathbf{V}_{i,j}\) captures the dissimilarity between the nodes, which ultimately guides the diffusion process. Accordingly, \(\operatorname{div}(\mathbf{V}\circ\mathbf{X})_{i}=\sum_{j:(i,j)\in\mathcal{E} }\mathbf{V}_{i,j}\odot\mathbf{x}_{j}\) measures the flow of density at node \(i\) in the direction of \(\mathbf{V}_{i,j}\). The nonlinear parameterization of the gradient further enhances the dynamics to adapt to different graphs with varying homophily levels.

Reaction-diffusion.A more general reaction-diffusion dynamics is considered in [21]:

\[\text{GREAD}:\quad\frac{\partial\mathbf{X}}{\partial t}=\operatorname{div}( \mathbf{A}(\mathbf{X})\cdot\nabla\mathbf{X})+R(\mathbf{X}),\]

where \(R(\mathbf{X})\) is the reaction term, and proper choice of \(R(\cdot)\) recovers many existing works as special cases. For example, the Fisher reaction [33] is given by \(R(\mathbf{X})=\kappa\mathbf{X}\odot(\mathbf{1}-\mathbf{X})\), which can be used to model the spread of biological populations where \(\kappa\) represents the intrinsic growth rate. Other reaction processes include Allen-Cahn [1]\(R(\mathbf{X})=\mathbf{X}\odot(1-\mathbf{X}\odot\mathbf{X})\) and Zeldovich [37]\(R(\mathbf{X})=\mathbf{X}\odot(\mathbf{X}-\mathbf{X}\odot\mathbf{X})\). Apart from the physics informed choices of reaction term, [21] also consider \(R(\mathbf{X})=\mathbf{X}^{0}\), which follows GRAND++, CGNN to incorporate a source term, and also proposes several high-pass reaction term based on graph structure, aiming to induce a sharpening effect. For example, the blurring-sharpening reaction is defined as \(R(\mathbf{X})=(\mathbf{A}(\mathbf{X})-\mathbf{A}(\mathbf{X})^{2})\mathbf{X}\), which corresponds to performing a low-pass filter \(\mathbf{A}(\mathbf{X})-\mathbf{I}\) followed by a high-pass filter \(\mathbf{I}-\mathbf{A}(\mathbf{X})\). The paper also considers two learnable coefficients \(\alpha,\beta\) to control the emphasis on diffusion and reaction terms respectively, i.e., \(\alpha\operatorname{div}(\mathbf{A}(\mathbf{X})\cdot\nabla\mathbf{X})+\beta R( \mathbf{X})\).

A closely related work [94] adopts the Allen-Cahn reaction term for the reaction-diffusion process. Further, the paper allows negative diffusion coefficients, which is able to induce a repulsive force between the nodes:

\[\text{ACMP}:\quad\frac{\partial\mathbf{X}}{\partial t}=\operatorname{div}(( \mathbf{A}(\mathbf{X})-\mathbf{B})\cdot\nabla\mathbf{X})+\mathbf{X}\odot( \mathbf{1}-\mathbf{X}\odot\mathbf{X})\]

where \(\mathbf{B}>0\) is a bias term controlling the strength and direction of message passing. This allows \(\mathbf{A}(\mathbf{X})-\mathbf{B}\) to model the interactive forces and can become negative. Further, the first term of ACMP corresponds to the gradient flow of a pseudo Dirichlet energy given by \(\sum_{(i,j)\in\mathcal{E}}(\mathbf{A}_{i,j}-\mathbf{B}_{i,j})\|\mathbf{x}_{i}- \mathbf{x}_{j}\|^{2}\) where we ignore the dependence of \(\mathbf{A}\) on \(\mathbf{X}\) for the time being. This suggests, when \(\mathbf{A}_{i,j}-\mathbf{B}_{i,j}>0\), node \(i\) is attracted by node \(j\) by minimizing the difference between \(\mathbf{x}_{i}\) and \(\mathbf{x}_{j}\) and when \(\mathbf{A}_{i,j}-\mathbf{B}_{i,j}<0\), node \(i\) is repelled by node \(j\). It is noticed that the presence of negative weights can cause the energy to be unbounded and thus the dynamics may not be convergent. To resolve the issue, the paper considers an external potential, namely the double-well potential \((\delta/4)\sum_{i\in\mathcal{V}}(1-\|\mathbf{x}_{i}\|^{2})^{2}\). ACMP is indeed derived as the gradient flow of the pseudo Dirichlet energy combined with the double-well potential. Theoretically, the Dirichlet energy \(\mathcal{E}_{\operatorname{dir}}\) of ACMP evolution is upper bounded due to the Allen-Cahn reaction term as well as lower bounded due to the repulsive forces. Hence, the system remains stable while avoiding oversmoothing. For practical implementation, ACMP sets \(\mathbf{B}=\beta>0\), a tunable hyperparameter for simplicity of optimization. Two channel-wise coefficient vectors are added to balance the diffusion and reaction similarly in [21].
The idea of incorporating repulsive force in the message passing has also been explored in [66]. The work proposes to view the message passing mechanism on graphs in the framework of opinion dynamics. The work explores the notion of _bounded confidence_ from the Hegselmann-Krause (HK) model [49] where only similar opinions (up to some cut-off threshold) are exchanged. This motivates the following dynamics on graphs that separates messages according to the similarity level of graph signals:

\[\text{ODNet}:\quad\frac{\partial\mathbf{X}}{\partial t}=\operatorname{div}_{ \mathcal{V}\times\mathcal{V}}(\Phi(\mathbf{A}(\mathbf{X}))\cdot\nabla\mathbf{ X})+R(\mathbf{X}),\]

where \(\operatorname{div}_{\mathcal{V}\times\mathcal{V}}\) defines message aggregation over the complete graph. \(\Phi(\mathbf{A}(\mathbf{X}))\) is an elementwise scalar-valued function (called influence function) on diffusivity and is required to be non-decreasing. In [66], \(\Phi(\cdot)\) is chosen to be piecewise linear as follows.

\[\Phi(s)=\begin{cases}\mu s,&\text{if }s>\epsilon_{2},\\ s,&\text{if }\epsilon_{1}\leq s\leq\epsilon_{2},\\ \nu(1-s),&\text{otherwise}\end{cases}\]

where \(\mu>0\) and \(\nu\leq 0\) are hyperparameters. In addition, \(\epsilon_{1},\epsilon_{2}\) defines the influence regions (which resembles bounded confidence in HK model). Because \(\nu\) can be negative, it is able to induce repulsive forces by separating the node representations. Empirically, the work chooses \(\nu=0\) for homophilic graphs and \(\nu<0\) for heterophilic graphs. It is noticed that when \(\nu<0\), the message can propagate even for unconnected nodes in the case of \(\nu<0\).

Advection-diffusion-reaction.ADR-GNN [30] further adds an explicit advection term on top of the reaction-diffusion process considered in [21].

\[\text{ADR-GNN}:\quad\frac{\partial\mathbf{X}}{\partial t}=\operatorname{div}( \mathbf{A}(\mathbf{X})\cdot\nabla\mathbf{X})+\operatorname{div}(\mathbf{V} \circ\mathbf{X})+R(\mathbf{X}).\]

The work considers homogeneous diffusion with channel scaling, i.e., \(\operatorname{div}(\mathbf{A}(\mathbf{X})\cdot\nabla\mathbf{X})=-\mathbf{L} \mathbf{X}\mathrm{diag}(\boldsymbol{\theta})\), where \(\boldsymbol{\theta}\in\mathbb{R}^{c}\) is the channel-wise scaling factor. In contrast to CDE [103], the advection term concerns two directional velocity fields \(\mathbf{V}_{i,j}\), \(\mathbf{V}_{j,i}\in\mathbb{R}^{c}\), which measures both in-flow and out-flow of density, with \((\mathbf{V}\circ\mathbf{X})_{i,j}\coloneqq\mathbf{V}_{j,i}\odot\mathbf{x}_{j }-\mathbf{V}_{i,j}\odot\mathbf{x}_{i}\). By further ensuring channel-wise row stochasticity of both \(\mathbf{V}_{i,j}\) and \(\mathbf{V}_{j,i}\), \(\mathbf{V}_{i,j}\odot\mathbf{x}_{i}\) is interpreted as the mass of node \(j\) to be transported to node \(i\). Thus the advection term, given by \(\operatorname{div}(\mathbf{V}\circ\mathbf{X})_{i}=\sum_{j:(i,j)\in\mathcal{E }}\mathbf{V}_{j,i}\odot\mathbf{x}_{j}-\mathbf{x}_{i}\), quantifies the net flow of density at node \(i\). The reaction term \(R(\mathbf{X})\) is parameterized by additive and multiplicative MLPs with a source term, \(R(\mathbf{X})=\sigma(\mathbf{X}\mathbf{W}_{1}+\tanh(\mathbf{X}\mathbf{W}_{2}) \odot\mathbf{X}+\mathbf{X}^{0}\mathbf{W}_{3})\). Different to previous works, the paper considers operator splitting for discretizing the continuous dynamics, i.e., by separating the propagation of the three processes. Such a scheme allows separate treatment and analysis of each process. Particularly, [30] demonstrate the mass preserving property and stability of the advection operator. Empirically, ADR-GNN has shown promising performance for modelling not only the static graphs but also spatial temporal graphs where advection-diffusion-reaction process has been successful [32].

Gating.It has been shown that controlling the speed of diffusion through convection/advection term is able to counteract the smoothing process. [79] adopt a different strategy by explicitly modelling a gating function on the diffusion:

\[\text{G}^{2}:\quad\frac{\partial\mathbf{X}}{\partial t}=\mathbf{T}(\mathbf{X} )\odot\operatorname{div}(\mathbf{A}(\mathbf{X})\cdot\nabla\mathbf{X})\]
where \(\mathbf{T}(\mathbf{X})\in[0,1]^{n\times c}\) collects the rate of speed for each node and across each channel. Specifically, the rates depend on the graph gradient as \(T(\mathbf{X})_{i,k}=\tanh(\sum_{j\in\mathcal{N}_{i}}|\hat{\mathbf{X}}_{j,k}- \hat{\mathbf{X}}_{i,k}|^{p}),p>0\) where \(\hat{\mathbf{x}}_{i}=\sum_{j\in\mathcal{N}_{i}}\widehat{\mathbf{A}}(\mathbf{X })_{i,j}\mathbf{x}_{j}\) and \(\widehat{\mathbf{A}}(\mathbf{X})\) is another message aggregation. Conceptually, the gating rates \(T(\mathbf{X})_{i,:}\) for node \(i\) depend on the channel-wise graph gradients to all its neighbours. The use of \(\tanh(\cdot)\) ensures when \(\sum_{j\in\mathcal{N}_{i}}|\hat{\mathbf{X}}_{j,k}-\hat{\mathbf{X}}_{i,k}|^{p}\to 0\), the rate \(T(\mathbf{X})_{i,k}\) vanishes at a faster rate. This correspondingly shuts down the update for node \(i\) and thus avoid oversmoothing. The paper also considers more general choices of coupling functions in place of \(\operatorname{div}(\mathbf{A}(\mathbf{X})\cdot\nabla\mathbf{X})\) where nonlinearity is added.

The idea of gating has been similarly explored in DeepGRAND [71], which utilizes a channel-wise scaling factor \(\langle\mathbf{X}\rangle^{p}\in\mathbb{R}^{n\times d}\) in place of \(T(\mathbf{X})\) where \(\langle\mathbf{X}\rangle^{p}_{:,k}=\|\mathbf{X}_{:,k}\|^{p}\mathbf{1}_{n}\). The dynamics also incorporates a perturbation to the diffusivity as \(\mathbf{A}(\mathbf{X})-(1+\epsilon)\mathbf{I}\). The scaling factor and perturbation help regulate the convergence of node features so that the node features neither explodes nor converges too fast to the steady state.

Reverse diffusion.Similar to the idea in [21] that simultaneously accounts for low-pass and high-pass filters, [81] introduce a reverse diffusion process based on the heat kernel. When coupled with heat diffusion, it leads to a process that accommodates both smoothing and sharpening effects:

\[\text{MHKG}:\quad\frac{\partial\mathbf{X}}{\partial t}=\big{(}\mathrm{diag}( \boldsymbol{\theta}_{1})\exp(f(\widehat{\mathbf{L}}))+\mathrm{diag}(\boldsymbol {\theta}_{2})\exp(g(\widehat{\mathbf{L}}))-\mathbf{I}\big{)}\mathbf{X},\]

where \(\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}\) are learnable filters and \(f,g\) are scalar-valued functions defined over the eigenvalues of \(\widehat{\mathbf{L}}\), e.g., \(f(\widehat{\mathbf{L}})\coloneqq\mathbf{U}f(\boldsymbol{\Lambda})\mathbf{U}^ {\top}\) where \(f(\boldsymbol{\Lambda})=\mathrm{diag}([f(\lambda_{i})]_{i=1}^{n})\), with \(\boldsymbol{\Lambda}\) being the diagonal matrix collecting the eigenvalues and \(\mathbf{U}\) collecting the eigenvectors of the normalized Laplacian \(\widehat{\mathbf{L}}\). In particular, \(f,g\) are assumed to be opposite in terms of monotonicity. In the simplest case, suppose \(f(\widehat{\mathbf{L}})=-\widehat{\mathbf{L}}\) and \(g(\widehat{\mathbf{L}})=\widehat{\mathbf{L}}\), then the two terms in MHKG correspond to the heat kernel and its reverse. The filtering coefficients \(\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}\) controls the relative dominance of the two terms, where the former smooths while the latter sharpens the signals.

Anti-symmetry.In [41], a stable and non-dissipative system is proposed by imposing the additional anti-symmetric structure for channel mixing matrix as

\[\text{A-DGN}:\quad\frac{\partial\mathbf{X}}{\partial t}=\mathbf{X}(\mathbf{W}- \mathbf{W}^{\top})+F_{\mathcal{G}}(\mathbf{X}),\]

where we omit the nonlinearity and a bias term to show only the driving factors. Here, \(F_{\mathcal{G}}(\mathbf{X})\) is a coupling function similar as in [78], such as a simple homogeneous message passing \(F_{\mathcal{G}}(\mathbf{X})_{i}=\sum_{j\in\mathcal{N}_{i}}\mathbf{V}\mathbf{x }_{j}\) for some weight matrix \(\mathbf{V}\) or the one with attention mechanism [91]. The incorporation of anti-symmetry constraint for the channel mixing renders the system to be stable and non-dissipative, both due to the fact that the Jacobian of dynamics has pure imaginary eigenvalues, i.e., all real parts of the eigenvalues are zero. This suggests the solutions to the system remain bounded under perturbation of initial conditions, which concludes the stability of the evolution. In addition, the zero real part of the eigenvalues suggests the sensitivity of the node signals to its initial values, i.e., the magnitude of \(\frac{\partial\mathbf{x}_{i}(t)}{\partial\mathbf{x}_{i}(0)},\forall i,t\) stays constant throughout the dynamics. This result infers that oversmoothing in the limit does not occur as the final state still depends on the initial conditions as \(\lim_{t\to\infty}\frac{\partial\mathbf{x}_{i}(t)}{\partial\mathbf{x}_{i}(0)}\neq 0\). Further, this also suggests the magnitude of \(\frac{\partial L_{\mathrm{loss}}}{\partial\mathbf{x}_{i}(0)}\) remains unchanged over time and hence gradient vanishing or explosion is avoided during backpropagation. This allows the dynamics to be propagating to the limit and capture long-range interactions without facing the issue of oversmoothing or training instability.
### Geometry-underpinned dynamics

Previous sections have viewed graphs from its trivial topology, and standard dynamics on graphs amounts to propagating information from node to edge space and back, only utilizing the connectivity between nodes. In fact, graphs usually possess more intricate geometries and viewed as discrete realization of general topological spaces. In this section, we show how dynamics on graphs can be underpinned with additional geometric structure, such as sheaves and stalks in [10], and cliques and cochains in [42].

Sheaf diffusion.[47, 4] and [10] leverage _cellular sheave theory_ to endow a geometric structure for graphs. Specifically, each node \(i\in\mathcal{V}\) and edge \(e_{ij}=\{i,j\}\in\mathcal{E}\) (undirected) is equipped with a vector space structure \(\mathcal{F}(i),\mathcal{F}(e_{ij})\), with a linear map \(\boldsymbol{\mathcal{F}}_{i\,\triangle\,e_{ij}}:\mathcal{F}(i)\to\mathcal{F}( e_{ij})\) (called restriction map) that connects the node to edge spaces. Its adjoint operator \(\boldsymbol{\mathcal{F}}_{i\,\triangle\,e_{ij}}^{\top}:\mathcal{F}(e_{ij}) \to\mathcal{F}(i)\) does the reverse. The direct sum of all vector spaces of nodes is called the space of \(0\)-cochains, denoted by \(C^{0}(\mathcal{G};\mathcal{F})\coloneqq\bigoplus_{i\in\mathcal{V}}\mathcal{F} (i)\). Suppose, without loss of generality, that all vector spaces \(\mathcal{F}(i),\mathcal{F}(e_{ij})\) are \(d\)-dimensional, we can represent \(\boldsymbol{\mathcal{F}}_{i\,\triangle\,e_{ij}}\) as a \(d\times d\) matrix. Further, suppose \(\mathbf{x}\in C^{0}(\mathcal{G};\mathcal{F})\), then each \(\mathbf{x}_{i}\in\mathbb{R}^{d}\) and \(\mathbf{x}\in\mathbb{R}^{nd}\) by stacking the feature vectors across all nodes.

Under the construction of vector spaces on nodes and edges, the concepts of graph gradient and divergence require to utilize the restriction maps because the vector spaces are not directly comparable. That is, the graph gradient (also known as the co-boundary map) is defined as \((\nabla_{\mathcal{F}}\mathbf{x})_{i,j}\coloneqq\boldsymbol{\mathcal{F}}_{j\, \triangle\,e_{ij}}\mathbf{x}_{j}-\boldsymbol{\mathcal{F}}_{i\,\triangle\,e_{ ij}}\mathbf{x}_{i}\), where the restriction maps transport the signals to a common disclosure space. The graph divergence is thus similarly defined as \((\operatorname{div}_{\mathcal{F}}\mathbf{G})_{i}=\sum_{j\in\mathcal{N}_{i}} \boldsymbol{\mathcal{F}}_{i\,\triangle\,e_{ij}}^{\top}\mathbf{G}_{e_{ij}}\), for \(\mathbf{G}_{e_{ij}}\in\mathcal{F}(e_{ij})\). This leads to the definition of sheaf Laplacian as \(\mathbf{L}_{\mathcal{F}}(\mathbf{x})_{i}=\operatorname{div}_{\mathcal{F}}( \nabla_{\mathcal{F}}\mathbf{x})_{i}=\sum_{j\in\mathcal{N}_{i}}\boldsymbol{ \mathcal{F}}_{i\,\triangle\,e_{ij}}^{\top}(\boldsymbol{\mathcal{F}}_{i\, \triangle\,e_{ij}}\mathbf{x}_{i}-\boldsymbol{\mathcal{F}}_{j\,\triangle\,e_{ ij}}\mathbf{x}_{j})\). The sheaf Laplacian is an \(nd\times nd\) block positive semi-definite matrix, with the diagonal blocks \((\mathbf{L}_{\mathcal{F}})_{i,i}=\sum_{j\in\mathcal{N}_{i}}\boldsymbol{ \mathcal{F}}_{i\,\triangle\,e_{ij}}^{\top}\boldsymbol{\mathcal{F}}_{i\, \triangle\,e_{ij}}\) and off-diagonal blocks \((\mathbf{L}_{\mathcal{F}})_{i,j}=-\boldsymbol{\mathcal{F}}_{i\,\triangle\,e_{ ij}}^{\top}\boldsymbol{\mathcal{F}}_{j\,\triangle\,e_{ij}}\). A symmetrically normalized sheaf Laplacian can be similarly computed by \(\mathbf{D}_{\mathcal{F}}^{-1/2}\mathbf{L}_{\mathcal{F}}\mathbf{D}_{\mathcal{F} }^{-1/2}\) with \(\mathbf{D}_{\mathcal{F}}\) is the block diagonal of \(\mathbf{L}_{\mathcal{F}}\). The corresponding sheaf diffusion process is

\[\frac{\partial\mathbf{X}}{\partial t}=\operatorname{div}_{\mathcal{F}}( \nabla_{\mathcal{F}}\mathbf{X})=-\mathbf{L}_{\mathcal{F}}\mathbf{X}\]

where here \(\mathbf{X}\in\mathbb{R}^{(nd)\times c}\), where \(c\) is the feature channels and \(\operatorname{div}_{\mathcal{F}},\nabla_{\mathcal{F}},\mathbf{L}_{\mathcal{F}}\) are applied channel-wise. The Sheaf diffusion turns out to be the gradient flow of the sheaf Dirichlet energy \(\operatorname{tr}(\mathbf{X}^{\top}\mathbf{L}_{\mathcal{F}}\mathbf{X})\), which measures the smoothness of signals in the disclosure space. For practical settings where sheaf structure is unavailable, one can construct such a feature through input embedding. It is worth highlighting that, when \(d=1\), the sheaf Laplacian reduces to the classic graph Laplacian and the sheaf diffusion becomes the heat diffusion. In [10], a variety of restriction maps are constructed, which leads to dynamics flexible enough to handle different types of graphs and avoid oversmoothing. The paper also develops a general framework for learning the sheaf Laplacian from the features and include channel mixing and nonlinearity to increase the expressive power:

\[\text{NSD}:\quad\frac{\partial\mathbf{X}}{\partial t}=-\sigma(\mathbf{L}_{ \mathcal{F}(\mathbf{X})}(\mathbf{I}\otimes\mathbf{W}_{1})\mathbf{X}\mathbf{W} _{2})\]

where \(\mathbf{W}_{1}\in\mathbb{R}^{d\times d}\) transforms the feature vectors and \(\mathbf{W}_{2}\in\mathbb{R}^{c\times c^{\prime}}\) mixes the channels and \(\otimes\) denotes the Kronecker product. \(\mathbf{L}_{\mathcal{F}(\mathbf{X})}\) is parameterized via a matrix-valued function on the current feature values.

In [10], the sheaf Laplacian is parameterized by \(\boldsymbol{\mathcal{F}}_{i\,\triangle\,e_{ij}}=\sigma(\mathbf{V}[\mathbf{x}_{ i}\|\mathbf{x}_{j}])\) where \(\cdot\|\cdot\) denotes concatenation. A follow-up work [5] devises a sheaf attention mechanism to further enhance the diffusion process. Let \(\mathbf{A}(\mathbf{X})\in\mathbb{R}^{n\times n}\) be a matrix of learnable attention coefficients (the same as in GAT [91]), and let
\(\mathbf{A}(\mathbf{X})\otimes\mathbf{1}_{d\times d}\) that assigns uniform attention coefficients for each feature dimension within a vector space. The attentive sheaf diffusion (SheafAN) is introduced as \(\frac{\partial\mathbf{X}}{\partial t}=(\widehat{\mathbf{A}}(\mathbf{X})\odot \mathbf{A}_{\mathcal{F}(\mathbf{X})}-\mathbf{I})\mathbf{X}\), where \((\mathbf{A}_{\mathcal{F}(\mathbf{X})})_{i,j}=\boldsymbol{\mathcal{F}}_{i\, \circ\,e_{ij}}^{\top}\boldsymbol{\mathcal{F}}_{j\,\circ\,e_{ij}}\) is the sheaf adjacency matrix with self-loop, i.e., \(e_{ii}\in\mathcal{E}\). A second-order sheaf PDE (NSP) is proposed in [84] using the wave equation as \(\frac{\partial^{2}\mathbf{X}}{\partial t^{2}}=\mathrm{div}_{\mathcal{F}}( \nabla_{\mathcal{F}}\mathbf{X})\).

Bracket dynamics.[42] propose to use _geometric brackets_ that implicitly parameterize dynamics on graphs that satisfy certain properties while equipping graphs with higher-order structures. The formulation requires concepts from structure-preserving bracket-based dynamics and exterior calculus. In general, for a state variable \(x\), its dynamics can be given by some combination of reversible and irreversible brackets. The _reversible_ bracket (also known as Poisson bracket) is denoted by \(\{A,E\}\coloneqq\langle\frac{\partial A}{\partial x},\tilde{L}\frac{\partial E }{\partial x}\rangle\) for some skew-symmetric operator \(\tilde{L}\)1 and some inner product \(\langle\cdot,\cdot\rangle\). The reversibility is a result of energy conservation. The _irreversible_ bracket is defined by \([A,E]\coloneqq\langle\frac{\partial A}{\partial x},M\frac{\partial E}{ \partial x}\rangle\) for some (either positive or negative) semi-definite operator \(M\). The irreversibility describes the loss of energy from the system due to friction or dissipation. The double bracket \(\{\{A,E\}\}\coloneqq\langle\frac{\partial A}{\partial x},\tilde{L}^{2}\frac{ \partial E}{\partial x}\rangle\) is an irreversible bracket.

Footnote 1: Here, in order not to be confused with the Laplacian \(L\) used in previous discussions, we use \(\tilde{L}\).

For simplicity, the paper considers \(A=x\) and one can simplify the brackets as \([x,E]=M\frac{\partial E}{\partial x}\) and \(\{x,E\}=\tilde{L}\frac{\partial E}{\partial x}\). The paper considers four different types of dynamics leveraging both reversible and irreversible brackets:

\[\begin{array}{ll}\text{Hamiltonian}:&\frac{\partial x}{\partial t}=\{x,E\} ;\\ \text{Gradient}:&\frac{\partial x}{\partial t}=-[x,E];\\ \text{Double bracket}:&\frac{\partial x}{\partial t}=\{x,E\}+\{\{x,E\}\};\\ \text{Metriplectic}:&\frac{\partial x}{\partial t}=\{x,E\}+[x,S],\end{array} \tag{3}\]

where \(E(x)\) is referred to as the energy of the state and \(S(x)\) is the entropy. The dynamics of each process captures fundamentally different systems and has natural substance in physics. The Hamiltonian dynamics leads to a complete, isolated system in the sense that no energy is lost to the external environment. In contrast, both the gradient and double bracket dynamics are incomplete where the energy is lost through the process. Metriplectic dynamics is complete by further requiring the degeneracy conditions \(\tilde{L}\frac{\partial S}{\partial x}=M\frac{\partial E}{\partial x}=0\). These conditions ensure the conservation of energy, i.e., \(\frac{\partial E}{\partial t}=0\) and the entropy inequality, i.e., \(\frac{\partial S}{\partial t}\geq 0\) in an isolated system.

In order to properly generalize the dynamics to discrete domains like graphs, one requires to identify the state variable, an inner product structure as well as energy and entropy. Rather than only considering node features as the state variable, [42] extend the framework to higher-order clique cochains, including edges and cycles. Formally, let \(\Omega_{k}\) be the set of \(k\)-cliques on a graph \(\mathcal{G}\), which contains ordered, complete, subgraphs generated by \((k+1)\)-nodes. For example, nodes, edges and triangles, correspond to the \(0\)-clique, \(1\)-clique and \(2\)-clique respectively. The exterior derivative operator is denoted as \(d_{k}:\mathfrak{F}(\Omega_{k})\rightarrow\mathfrak{F}(\Omega_{k+1})\) where \(\mathfrak{F}(\Omega)\) represents a function space over the domain \(\omega\). The specific structure of \(\mathfrak{F}\) depends on the chosen inner product \(\langle\cdot,\cdot\rangle\). The dual derivative \(d_{k}^{*}:\mathfrak{F}(\Omega_{k+1})\rightarrow\mathfrak{F}(\Omega_{k})\) is given as the adjoint of the \(d_{k}\) that satisfies \(\langle d_{k}f,G\rangle_{k+1}=\langle f,d_{k}^{*}G\rangle_{k}\) for any \(f\in\mathfrak{F}(\Omega_{k}),G\in\mathfrak{F}(\Omega_{k+1})\). One common choice of \(\mathfrak{F}\) is the \(L^{2}\) space, where one can derive \(d_{k}^{*}=d_{k}^{\top}\). For example, \(d_{0}\), is the graph gradient on nodes and \(d_{0}^{*}\) becomes the graph divergence on edges. The classic graph Laplacian can be computed as \(d_{0}^{*}d_{0}=d_{0}^{\top}d_{0}=\mathbf{G}^{\top}\mathbf{G}\) as we have shown in Section 2.

Instead, the paper pursues an inner product parameterized by positive definite matrices \(\mathbf{A}_{0},\mathbf{A}_{1},...,\mathbf{A}_{k}\) up to \(k\)-cliques. For example on node space \(\Omega_{0}\), where the \(L^{2}\) space has inner product \(\mathbf{f}^{\top}\mathbf{g}\) for \(\mathbf{f},\mathbf{g}\in\mathfrak{F}(\Omega_{0})\)
the generalized inner product is given by \(\mathbf{f}^{\top}\mathbf{A}_{0}\mathbf{g}\). Under such choice, one can show \(d_{k}^{\star}=\mathbf{A}_{k}^{-1}d_{k}^{\top}\mathbf{A}_{k+1}\).

The state variable is set to be a node-edge feature pairs, denoted by \(\mathbf{x}=(\mathbf{q},\mathbf{p})\), which can be treated as the position and momentum of a phase space. Further, the following operators are utilized to extend the dynamics in (3) to graphs,

\[\tilde{\mathbf{L}}=\begin{pmatrix}0&-d_{0}^{\ast}\\ d_{0}&0\end{pmatrix},\qquad\tilde{\mathbf{G}}=\begin{pmatrix}d_{0}^{\ast}d_{0} &0\\ 0&d_{1}^{\ast}d_{1}+d_{0}d_{0}^{\ast}\end{pmatrix},\quad\text{and}\quad\tilde{ \mathbf{M}}=\begin{pmatrix}0&0\\ 0&\mathbf{A}_{1}d_{1}^{\ast}d_{1}\mathbf{A}_{1}\end{pmatrix}.\]

It can be verified that \(\tilde{\mathbf{L}}\) is skew-symmetric and \(\tilde{\mathbf{G}},\tilde{\mathbf{M}}\) are symmetric positive definite with respect to the block-diagonal inner product parameterized by \(\mathbf{A}=\mathrm{diag}(\mathbf{A}_{0},\mathbf{A}_{1})\). Furthermore, let \(\mathbf{X}=(\mathbf{Q},\mathbf{P})\) denote the tuple of node and edge feature matrices, the energy considered is the total kinetic energy on both node and edge spaces, i.e., \(E(\mathbf{X})=\frac{1}{2}(\|\mathbf{Q}\|^{2}+\|\mathbf{P}\|^{2})\). The gradient with respect to the generalized inner product (called \(\mathbf{A}\)-gradient) can be computed as \(\nabla_{\mathbf{A}}E(\mathbf{X})=[\mathbf{A}_{0}^{-1}\frac{\partial E}{ \partial\mathbf{Q}},\mathbf{A}_{1}^{-1}\frac{\partial E}{\partial\mathbf{P}} ]^{\top}=[\mathbf{A}_{0}^{-1}\mathbf{Q},\mathbf{A}_{1}^{-1}\mathbf{P}]^{\top}\). For the Metriplectic dynamics, it is in general non-trivial to identify an entropy such that the degeneracy conditions hold. Hence the paper constructs a separate energy and entropy function pair as \(E_{m}(\mathbf{X})=f_{E}(\mathbf{Q})+g_{E}(d_{0}d_{0}^{\top}\mathbf{P})\) and \(S_{m}(\mathbf{X})=g_{S}(d_{1}^{\top}d_{1}\mathbf{P})\), for some node function \(f_{E}\) and edge functions \(g_{E},g_{S}\) applied channel-wise. The \(\mathbf{A}\)-gradient is derived as

\[\nabla_{\mathbf{A}}E_{m}(\mathbf{X})=\begin{pmatrix}\mathbf{A}_{0}^{-1} \mathbf{1}\otimes\nabla_{\mathbf{A}}f_{E}(\mathbf{Q})\\ d_{0}d_{0}^{\top}\mathbf{1}\otimes\nabla_{\mathbf{A}}g_{E}(d_{0}d_{0}^{\top} \mathbf{P})\end{pmatrix},\qquad\qquad\nabla_{\mathbf{A}}S_{m}(\mathbf{X})= \begin{pmatrix}0\\ \mathbf{A}_{1}^{-1}d_{1}^{\top}d_{1}\mathbf{1}\otimes\nabla_{\mathbf{A}}g_{S}( d_{1}d_{1}^{\top}\mathbf{P})\end{pmatrix}.\]

Importantly, the degeneracy conditions \(\tilde{\mathbf{L}}\nabla_{\mathbf{A}}S=\tilde{\mathbf{M}}\nabla_{\mathbf{A}} E=0\) are satisfied by construction.

Finally the generalized dynamics from (3) to graphs are

\[\text{Hamiltonian}_{G}: \frac{\partial\mathbf{X}}{\partial t}=\tilde{\mathbf{L}}(\mathbf{ X})\nabla_{\mathbf{A}}E(\mathbf{X});\] \[\text{Gradient}_{G}: \frac{\partial\mathbf{X}}{\partial t}=-\tilde{\mathbf{G}}(\mathbf{ X})\nabla_{\mathbf{A}}E(\mathbf{X});\] \[\text{Double bracket}_{G}: \frac{\partial\mathbf{X}}{\partial t}=\tilde{\mathbf{L}}(\mathbf{ X})\nabla_{\mathbf{A}}E(\mathbf{X})+\tilde{\mathbf{L}}^{2}(\mathbf{X})\nabla_{ \mathbf{A}}E(\mathbf{X});\] \[\text{Metrplectic}_{G}: \frac{\partial\mathbf{X}}{\partial t}=\tilde{\mathbf{L}}(\mathbf{ X})\nabla_{\mathbf{A}}E_{m}(\mathbf{X})+\tilde{\mathbf{M}}(\mathbf{X})\nabla_{ \mathbf{A}}S_{m}(\mathbf{X}),\]

where the operators \(\tilde{\mathbf{L}}(\mathbf{X}),\tilde{\mathbf{G}}(\mathbf{X}),\tilde{\mathbf{ M}}(\mathbf{X})\) are state-dependent through attention mechanism to construct the metric tensor \(\mathbf{A}_{0},\mathbf{A}_{1}\), which thus parameterizes the exterior derivative operators.

**Remark 1** (Connection to GCN and GAT/GRAND).: GCN can be seen as the Gradient\({}_{G}\) dynamics with \(\mathbf{A}_{0},\mathbf{A}_{1}\) parameterized by node degrees. Setting \(\mathbf{A}_{0}\) as a diagonal matrix (on nodes) with diagonal entries \(a_{0,ii}=\sqrt{\deg_{i}}\) and \(\mathbf{A}_{1}\) as the identity matrix (over edges), we can verify \((d_{0}^{\ast}d_{0}\mathbf{A}_{0}^{-1}\mathbf{Q})_{i}=\sum_{j\in\mathcal{N}_{i }}(\mathbf{q}_{i}/\deg_{i}-\mathbf{q}_{j}/\sqrt{\deg_{i}\deg_{j}})=(\tilde{ \mathbf{L}}\mathbf{Q})_{i}\). In this case, Gradient\({}_{G}\) recovers the heat equation (with normalized Laplacian) when \(\mathbf{P}=0\) and thus leads to GCN under discretization.

Similarly, GAT can be seen as the same Gradient\({}_{G}\) dynamics while learning a metric structure \(\mathbf{A}_{0},\mathbf{A}_{1}\). That is, choose \(a_{0,ii}=\sqrt{\sum_{j\in\mathcal{N}_{i}}\exp(\mathrm{attn}(\mathbf{q}_{i}, \mathbf{q}_{j}))}\) and \(a_{1,e_{ij}}=\exp(\mathrm{attn}(\mathbf{q}_{i},\mathbf{q}_{j}))\). Let the attention coefficient be \(a(\mathbf{q}_{i},\mathbf{q}_{j})=a_{1,e_{ij}}/a_{0,ii}\). Then one can show Gradient\({}_{G}\) dynamics with \(\mathbf{P}=0\) corresponds to a (symmetrically) normalized version of GRAND.

This result demonstrates the irreversible nature of GCN or GAT/GRAND dynamics where energy dissipates, while other dynamics are either conservative or partially dissipative.
Hamiltonian mechanics.The _Hamiltonian mechanics_ has also been considered in [57] where instead of using the edge features as the momentum, the work parameterizes the momentum with a neural network from the node features. This leads to distinction compared to the previous works in the evolution of the node features, which is decoupled from the graph structure. Let \(\mathbf{Q}\) denote the node features and the momentum is computed as \(\mathbf{P}=\mathrm{MLP}_{\Theta}(\mathbf{Q})\). The Hamiltonian mechanics is determined by a Hamiltonian function \(\mathcal{H}(\mathbf{Q},\mathbf{P})\), which characterizes the total energy of the system. The dynamics is governed by the Hamiltonian equation

\[\frac{\partial\mathbf{Q}}{\partial t}=\frac{\partial\mathcal{H}}{\partial \mathbf{P}},\qquad\frac{\partial\mathbf{P}}{\partial t}=-\frac{\partial \mathcal{H}}{\partial\mathbf{Q}}.\]

The paper motivates a variety of parameterization for the Hamiltonian \(\mathcal{H}\). One specific choice of Hamiltonian is \(\mathcal{H}(\mathbf{Q},\mathbf{P})=\mathrm{tr}(\mathbf{P}^{\top}\mathbf{M}( \mathbf{Q})\mathbf{P})\) where \(\mathbf{M}(\mathbf{Q})\) represents the inverse metric tensor at \(\mathbf{Q}\) (also learnable in the local neighbourhood). Its solution \(\mathbf{Q}(t)\) recovers the geodesic (a locally shortest curve) on the manifold with metric parameterized by \(\mathbf{M}(\mathbf{Q})^{-1}\) at \(\mathbf{Q}\). Hence, node is effectively embedded to a (implicit) manifold space where the metric is learnable. Unlike previous methods, where the dynamics of the features depend on a coupling function regulated by the graph topology, here the evolution of \(\mathbf{Q}\) is independent across the nodes. To further incorporate the graph structure, message passing by neighbourhood aggregation is performed after the feature evolution via Hamiltonian equation. Multiple layers of Hamiltonian dynamics and message passing are stacked to model complex geometries and node embeddings.

A follow-up work [102] extends the formulation by considering general graph-coupled Hamiltonian function \(\mathcal{H}_{\mathcal{G}}(\mathbf{Q},\mathbf{P})\). For example, the paper considers a Hamiltonian function defined as the norm of the output from a two-layer GCN, where \(\frac{\partial\mathcal{H}_{\mathcal{G}}}{\partial\mathbf{Q}},\frac{\partial \mathcal{H}_{\mathcal{G}}}{\partial\mathbf{P}}\) are computed from auto-differentiation. Further, the paper studies various notions of stability on graphs from the theory of dynamical system, including BIBO, Lyapunov, structural and conservative stability. The work conducts a systematic analysis and comparison of proposed Hamiltonian-based dynamics with existing graph neural dynamics, such as GRAND, BLEND, Mean Curvature and Beltrami. It is found that the conservative Hamiltonian dynamics has shown improved robustness against adversarial attacks.

### Dynamics as gradient flow

Most of the aforementioned designs of GNNs are inspired by evolution of some underlying dynamics. The learnable parameters, such as channel mixing, are usually added after discretization to increase the expressive power. In [39], the dynamics is instead given as the _gradient flow_ of some learnable energy. The framework is general and includes many of the existing works as special cases (as long as the channel mixing matrix is symmetric)2. The parameterized energy takes the following form:

Footnote 2: Although many existing works can be written as gradient flow of some energy, their motivation comes mostly from the dynamics, not from the energy.

\[\mathcal{E}_{\theta}(\mathbf{X})=\frac{1}{2}\mathrm{tr}(\mathbf{X}^{\top} \mathbf{X}\mathbf{\Omega})-\frac{1}{2}\mathrm{tr}(\mathbf{X}^{\top}\mathbf{A} \mathbf{X}\mathbf{W})+\varphi^{0}(\mathbf{X},\mathbf{X}^{0}),\]

where \(\mathbf{A}\) is the (normalized) adjacency matrix and \(\mathbf{\Omega},\mathbf{W}\in\mathbb{R}^{c\times c}\) are assumed to be symmetric3. The first term determines the external forces exerted upon the system and the second term reflects the pairwise interactions while the last term quantifies the energy preserved by the source term \(\mathbf{X}^{0}\). Although \(\varphi^{0}\) can be general, the paper considers a form \(\varphi^{0}(\mathbf{X},\mathbf{X}^{0})=\mathrm{tr}(\mathbf{X}^{\top}\mathbf{X }^{0}\tilde{\mathbf{W}})\). The gradient flow of \(\mathcal{E}_{\theta}(\mathbf{X})\) yields the
dynamics of the following general form,

\[\text{GRAFF}:\quad\frac{\partial\mathbf{X}}{\partial t}=-\nabla\mathcal{E}_{ \theta}(\mathbf{X})=-\mathbf{X}\mathbf{\Omega}+\mathbf{A}\mathbf{X}\mathbf{W}- \mathbf{X}^{0}\tilde{\mathbf{W}}.\]

This formulation includes many of the existing dynamics-motivated GNNs. When \(\mathbf{\Omega}=\mathbf{W},\tilde{\mathbf{W}}=0\), this corresponds to the evolution of (residual) GCN or GAT/GRAND [15] if \(\mathbf{A}\) is constructed by attention mechanism. If \(\tilde{\mathbf{W}}\neq 0\), this corresponds to the GRAND++ [86] and thus also the CGNN [99]. The decrease of the general energy does not necessarily lead to a decrease in the Dirichlet energy (which is a special case of \(\mathcal{E}_{\theta}(\mathbf{X})\) with \(\mathbf{\Omega}=\mathbf{W}=\mathbf{I},\varphi^{0}=0\)). This is mainly due to the occurrence of both attractive and repulsive effects along the positive and negative eigen-directions of \(\mathbf{W}\). More formally, one decompose \(\mathbf{W}=\mathbf{\Theta}_{+}^{\top}\mathbf{\Theta}_{+}-\mathbf{\Theta}_{- }^{\top}\mathbf{\Theta}_{-}\) and rewrite the energy (without \(\varphi^{0}\)) as

\[\mathcal{E}_{\theta}(\mathbf{X})=\frac{1}{2}\sum_{i\in\mathcal{V}}\langle \mathbf{x}_{i},(\mathbf{\Omega}-\mathbf{W})\mathbf{x}_{i}\rangle+\frac{1}{4} \sum_{(i,j)\in\mathcal{E}}\|\mathbf{\Theta}_{+}(\nabla\mathbf{X})_{i,j}\|^{2 }-\frac{1}{2}\sum_{(i,j)\in\mathcal{E}}\|\mathbf{\Theta}_{-}(\nabla\mathbf{F} )_{i,j}\|^{2}.\]

It has been shown that the gradient flow minimizes the gradient along the positive eigen-directions (which leads to smoothing effect) while maximizing the gradient along the negative eigen-directions (which leads to sharpening effect). In contrast, minimizing the Dirichlet energy always induces smoothing effect because \(\mathbf{W}=\mathbf{I}\) with only positive eigen-directions. This allows GRAFF to avoid oversmoothing and produce sharpening effects as long as the \(\mathbf{W}\) has sufficiently large negative spectrum.

### Multi-scale diffusion

Previous sections have mostly focused on dynamics with local diffusion. In other words, the signal/density at certain node changes depending on its immediate neighbourhood. The communication with distant nodes only happens when diffusion time is sufficiently large. Section 3.3 discusses dynamics that leverages non-local diffusion, but nonetheless mostly restricted to a single scale at each diffusion step, e.g., single \(\alpha\) in fractional diffusion [68] and \(t_{\theta}\) for time derivative diffusion [7]. Multi-scale graph neural networks, such as ChebyNet [23], LanczosNet [63] and Framelet GNN [104] are capable of capturing multi-scale graph properties through spectral filtering on the graph Laplacian. The eigen-pairs of graph Laplacian encode structural information at different levels of granularity and separate processing of each resolution provides insights into both local and global patterns.

Several recent works have adapted the continuous dynamics formulation to multi-scale GNN. [46] introduce a multi-scale diffusion process via _graph framelet_. Apart from the multi-scale properties shared with other spectral GNNs, graph framelet further separates the low-pass from high-pass filters. Let \(\mathcal{W}_{0,J}\in\mathbb{R}^{n\times n}\) denote the low-pass framelet transform and \(\mathcal{W}_{r,j}\in\mathbb{R}^{n\times n}\), \(r=1,...,R\), \(j=1,...,J\) denote the high-pass framelet transforms, where \(R\) is the number of high-pass filter banks and \(J\) is the scale level. For a multi-channel graph signal \(\mathbf{X}\in\mathbb{R}^{n\times c}\), \(\mathcal{W}_{0,J}\mathbf{X}\) and \(\mathcal{W}_{r,j}\mathbf{X}\), \(r=1,...,R\), \(j=1,...,J\) represent low-pass and high-pass framelet coefficients. For notation clarity let \(\mathcal{I}=\{(0,J)\}\cup\{(r,j)\}_{1\leq r\leq R,1\leq j\leq J}\) be the framelet index set. Due to the reconstruction property of framelets, it satisfies that \(\sum_{(r,j)\in\mathcal{I}}\mathcal{W}_{r,j}^{\top}\mathcal{W}_{r,j}\mathbf{X}= \mathbf{X}\).

The spatial framelet diffusion proposed in [46] is given as

\[\text{GradFUFG}:\quad\frac{\partial\mathbf{X}}{\partial t}=-\sum_{(r,j)\in \mathcal{I}}\big{(}\mathcal{W}_{r,j}^{\top}\mathcal{W}_{r,j}\mathbf{X} \mathbf{\Omega}_{r,j}-\mathcal{W}_{r,j}^{\top}\mathbf{\widehat{A}}\mathcal{W}_ {r,j}\mathbf{X}\mathbf{W}_{r,j}\big{)},\]

where \(\mathbf{\Omega}_{r,j},\mathbf{W}_{r,j}\) are symmetric channel mixing matrices. When \(\mathbf{\Omega}_{r,j}=\mathbf{W}_{r,j}=\mathbf{I}\), by the tightness of framelet transform, GradFUFG reduces to the heat equation as \(\frac{\partial\mathbf{X}}{\partial t}=-\mathbf{\widehat{L}}\mathbf{X}\). The spatial-based framelet
diffusion can be seen as gradient flow of a generalized Dirichlet energy, which also includes the parameterized energy considered by GRAFF [39] as a special case. It can be seen, due to the separate modelling of low-pass and high-pass as well as different scale levels, the dynamics can adapt to datasets of different homophily levels as well as has the potential of avoiding oversmoothing.

## 4 Numerical solvers for differential equations on graphs

Most of the works based on continuous GNN dynamics consider the following architecture, which is firstly introduced in GRAND [15].

\[\mathbf{X}(0)=\mathrm{Emb}(\mathbf{X}_{\mathrm{in}}),\qquad\mathbf{X}(T)= \mathbf{X}(0)+\int_{0}^{T}\frac{\partial\mathbf{X}(t)}{\partial t}dt,\qquad \mathbf{Y}=\mathrm{Out}(\mathbf{X}(T)),\]

where \(\mathrm{Emb}(\cdot),\mathrm{Out}(\cdot)\) are respectively the input embedding and output decoding functions, which are usually learnable. For second-order dynamics like in PDE-GCN\({}_{\mathrm{H}}\) and GraphCON, the same formulation applies if the second-order ODE is rewritten as a system of first-order ODEs.

There exist a variety of numerical solvers for differential equations, once the continuous dynamics is formulated on graphs. In fact, with the discrete differential operators defined on graphs, the PDE reduces to an ODE, which can be solved with standard numerical integrators. One natural strategy for discretizing a continuous dynamics is through _finite differences_. Particularly, the forward Euler and leapfrog methods are commonly employed for first-order and second-order ODE respectively. For a first-order ODE given by \(\frac{\partial\mathbf{X}}{\partial t}=F(\mathbf{X},t)\), the forward Euler discretization leverages forward finite time difference, which gives the update \(\mathbf{X}(t+1)=\mathbf{X}(t)+\tau F(\mathbf{X}(t),t)\) for some stepsize \(\tau>0\). The forward Euler is an explicit method in that \(\mathbf{X}(t+1)\) depends on \(\mathbf{X}(t)\) explicitly. In contrast, backward Euler method is implicit by discretization \(\mathbf{X}(t+1)=\mathbf{X}(t)+\tau F(\mathbf{X}(t+1),t+1)\), which involves solving a (nonlinear) system of equations to obtain \(\mathbf{X}(t+1)\).

For a second-order ODE \(\frac{\partial^{2}\mathbf{X}}{\partial t}=F(\mathbf{X},t)\), the leapfrog method that leverages centered finite differences yields the update \(\mathbf{X}(t+1)=2\mathbf{X}(t)-\mathbf{X}(t-1)+\tau^{2}F(\mathbf{X}(t),t)\), which has been considered for PDE-GCN\({}_{H}\)[29]. The leapfrog method can be equivalently derived by rewriting the the second-order ODE into a system of first-order ODEs (concerning both the position and velocity) and then use forward Euler method to alternatively update the two. In GraphCON [78], similar idea has been applied for the more complex second-order dynamics.

Higher-order methods employ high-order approximations for solving both first- or second-order ODEs and can be either explicit or implicit. Common choices include Runge-Kutta, which is a class of single-step methods that requires evaluating \(F\) at multiple extrapolated state. One popular variant of Runge-Kutta method is Runge-Kutta 4 that balances the accuracy and efficiency. It computes \(\mathbf{X}(t+1)=\mathbf{X}+\frac{\tau}{\mathfrak{e}}(\mathbf{K}_{1}+2\mathbf{ K}_{2}+2\mathbf{K}_{3}+\mathbf{K}_{4})\) where \(\mathbf{k}_{1}=F(\mathbf{X}(t),t),\mathbf{K}_{2}=F(\mathbf{X}(t)+\tau\mathbf{K }_{1}/2,t+\tau/2)\), \(\mathbf{K}_{3}=F(\mathbf{X}(t)+\tau\mathbf{K}_{2}/2,t+\tau/2)\), \(\mathbf{K}_{4}=F(\mathbf{X}(t)+\tau\mathbf{K}_{3},t+\tau)\). Multi-step methods, such as Adams-Bashford, store multiple previous steps in order to approximate high-order derivatives. More advanced solvers set the adaptive stepsize based on the error estimated at each iteration. Dormand-Prince is a Runge-Kutta method with step-size control and has been widely adopted as the default solver for ODEs [26].

For the backpropgation, the adjoint sensitivity method introduced in [19] can be used, which requires to solve an adjoint ODE for the gradient.
On oversmoothing, oversquashing, heterophily and robustness of GNN dynamics

As we have briefly mentioned, the main hurdles for successful designs of GNN include oversmoothing, oversquashing, graph heterophily and adversarial perturbations. The framework of continuous GNNs allows principled analysis and treatment for the issues identified. In fact, many of the previously introduced dynamics are motivated to overcome the limitations of existing GNN designs. Nonetheless, there exist possible trade-offs such as between oversmoothing and oversquashing [40, 81], which makes it challenging for GNN dynamics to avoid both phenomenon at the same time. Conceptually, common strategies for mitigating oversquashing relies on graph rewiring [88] that increases propagation strength, which likely enhances the smoothing effects. Apart from the trade-off, new problems may emerge as a result of resolving the existing ones. One example is the training difficulty (from gradient vanishing or explosion) associated with complex and long-range dynamics.

This section provides a holistic overview on these undesired behaviours of GNNs through the lens of continuous dynamics, along with a discussion on how existing approaches address the issues.

Oversmoothing.Oversmoothing refers to a phenomenon that node signals become progressively similar as a result of message passing, and converge to a state independent of the input signals. Although there are many different characterizations for oversmoothing, most if not all rely on the (normalized) Dirichlet energy that measures node similarity. Recall from Section 2, the Dirichlet energy is defined as \(\mathcal{E}_{\mathrm{dir}}(\mathbf{X})=\frac{1}{2}\sum_{(i,j)\in\mathcal{E}} \|\mathbf{x}_{i}/\sqrt{\mathrm{deg}_{i}}-\mathbf{x}_{j}/\sqrt{\mathrm{deg}_{ j}}\|^{2}=\mathrm{tr}(\mathbf{X}^{\top}\widehat{\mathbf{L}}\mathbf{X})\), where \(\widehat{\mathbf{L}}\) is the symmetrically normalized Laplacian. Because oversmoothing claims a loss of distinguishability across the nodes, a natural quantification of oversmoothing is \(\mathcal{E}_{\mathrm{dir}}(\mathbf{X}(t))\to 0\) as \(t\to\infty\). From the definition of the Dirichlet energy, in the limiting state, nodes with same degree collapse into a single representation. A stronger notion of oversmoothing has been considered in [78, 77] requiring an exponential decay of the Dirichlet energy for oversmoothing to occur, i.e., \(\mathcal{E}_{\mathrm{dir}}(\mathbf{X}(t))\leq C_{1}e^{-\mathcal{C}_{2}t}, \forall t>0\) and for some constants \(C_{1},C_{2}>0\). This is equivalent to claiming that the system has a node-wise constant, exponentially stable equilibrium state. In this work, we focus on a notion of oversmoothing called low-frequency dominance (LFD) put forward by [39]. A dynamics is said to be low-frequency dominant if \(\mathcal{E}_{\mathrm{dir}}(\mathbf{X}(t)/\|\mathbf{X}(t)\|)\to 0\). It can be readily verifies that when \(\mathcal{E}_{\mathrm{dir}}(\mathbf{X}(t))\to 0\), then \(\mathcal{E}_{\mathrm{dir}}(\mathbf{X}(t)/\|\mathbf{X}(t)\|)\to 0\), while the reverse argument is false (see Appendix B.2 of [39] more discussions). The normalization by \(\|\mathbf{X}(t)\|\) reveals the dependence of asymptotic behaviour to the dominant spectrum of the dynamics, where low-frequency is related to the smoothing effect.

In the sense of LFD, it can be shown that the dynamics relying on (anisotropic) diffusion would incur oversmoothing in the limit, such as the linear versions of GRAND [15], BLEND [16], and DIFFormer [98]. See for example [39, Theorem B.4] for a formal proof. Dynamics that involves (non-smooth) edge indicators, like Mean Curvature, Beltrami flows [83], \(p\)-Laplacian [35] and DIGNN [34] also correspond to LFD dynamics although the convergence is slower compared to (smooth) diffusion dynamics. Such a statement has been formalized in [34].

We summarize existing remedies for oversmoothing as follows.

1. _Source term_: One remedy that steers the dynamics away from the low-frequency dominance is to include a source term. For example, in \(p\)-Laplacian [35] and DIGNN [34], the source term is implicitly added in terms of input dependent regularization, while GRAND++ [86] and CGNN [99] explicitly inject the source information to the dynamics. The presence of a source term ensures the Dirichlet energy is lower bounded above zero where the limiting state is non-constant.
2. _Non-smoothing dynamics_: Choosing a non-smoothing dynamics can help avoid oversmoothing, such as an oscillatory process as in PDE-GCN\({}_{\mathrm{H}}\)[29] and GraphCON [78]. Oscillatory system usually converses rather than dissipates energy. In particular, it has been proved that the constant equilibrium state of the dynamics is not exponentially stable because any perturbation remains due to the energy conservation. The reversible dynamics, such as Hamiltonian and metriplectic in [42] is also non-smoothing by construction.
3. _Diffusion modulators_: Imposing external forces to counteract or modulate the diffusion mechanism is beneficial for mitigating oversmoothing. Examples include the high-pass filters in GREAD [21] and GradFUFG [46], repulsive diffusion coefficients in ACMP [94] and ODNet [66], negative spectrum in GRAFF [39], reverse heat kernel in MHKG [81], anti-symmetry in A-DGN [41], diffusion gating in \(\mathrm{G}^{2}\)[79] and DeepGRAND [71]. In essence, oversmoothing can be avoided in the above dynamics by diminishing smoothing effects in the limit of evolution.
4. _Higher-order geometries_: [10] demonstrate the cause of oversmoothing from the perspective of choosing a trivial geometry. In contrast, by enriching the graph with a higher-order sheaf topology, the asymptotic behaviour of the diffusion process can be better controlled. Indeed, the sheaf Dirichlet energy can be increased with a suitable choice of non-symmetric sheaf structure on graphs, thus avoiding oversmoothing under this regime [10, Proposition 17].

Heterophily.Unlike homophilic graphs where neighbouring nodes tend to share similar features and labels, nodes in heterophilic graphs usually exhibit dissimilar patterns compared to the neighbours. GNN that is dominated by smoothing effect is generally less preferred for such graphs. Following [39], the ability of a dynamics to adapt for heterophilic graphs is measured in terms of high-frequency dominance (HFD). A dynamics is said to be high-frequency dominant if \(\mathcal{E}_{\mathrm{dir}}(\mathbf{X}(t)/\|\mathbf{X}(t)\|)\to\rho_{L}\), where \(\rho_{L}\leq 2\) is the largest eigenvalue of \(\widehat{\mathbf{L}}\). Intuitively, HFD dynamics is dominated by a sharpening effect and eventually converge to the highest-frequency eigenvectors of the Laplacian where node information gets separated. In this work, we equate the ability of a dynamics handling heterophily with its ability to become dominated by the highest frequency. Although HFD dynamics is also undesired due to information loss related to the low frequencies, here we focus on the ability, meaning that a system should have the flexibility to converge to the eigenspace associated with the highest frequency, which allows separability of bipartite graphs. Under such notion, dynamics that is purely smoothing cannot be HFD. The incorporation of a source term and (positive) graph re-weighting alone are not sufficient to turn a smoothing dynamics into a sharpening one even though empirically these strategies boost performance on heterophilic datasets. This is because the driving mechanism of the resulting system is still diffusion.

Existing dynamics that provably accommodates separating effects can be classified as follows.

1. _Sharpening induced forces_: One scheme for a system to be HFD is to explicitly introduce sharpening induced forces. Most of the diffusion modulators identified for tackling oversmoothing amounts to incorporate such anti-smoothing forces, including GREAD [21], GradFUFG [46], ACMP [94], ODNet [66], GRAFF [39], and MHKG [81]. Conceptually, the dynamics is dominated by the high frequency if the magnitude of the sharpening effect surpasses the smoothing effect.
2. _Higher-order geometries_: Once the graph is equipped with a sheaf topology, [10] have shown that choosing a higher-dimensional stalks and non-symmetric restriction maps allows sheaf diffusion to gain linear separability in the limit for heterophilic graphs.
We highlight that other schemes, such as oscillatory dynamics in GraphCON [78] and convection/advection process in [103, 30] also demonstrate empirical success under heterophilic settings. Nevertheless, theoretical guarantee for theses schemes in adapting to heterophilic graphs is currently unavailable.

Oversquashing.The phenomenon of oversquashing has firstly been empirically observed by [2] and later formally studied by [88, 25]. Such a phenomenon has been identified as a key factor hindering expressive power of GNNs [25]. Oversquashing concerns a scenario where long-range communication between distant nodes matters for the downstream tasks whereas message passing squashes information due to the exponentially growing size of the receptive field. A closely related concept is under-reaching, which refers to shallow GNN not being able to fully explore long-range interactions [6]. A critical difference between under-reaching and oversquashing is that the former can often be alleviated by increasing the depth of a GNN while the latter can occur even with deep GNNs. One major cause of such a phenomenon is the local diffusion process according to the given graph structure, where messages are exchanged only within the neighbourhood. This can be formally characterized by the sensitivity across nodes, measured through the Jacobian. More formally, consider a (discretized) dynamics \(\mathbf{X}^{\ell+1}=\mathcal{A}_{\mathcal{G}}(\mathbf{X}^{\ell})\), where \(\mathcal{A}_{\mathcal{G}}\) is a coupling function that aggregates information based on input graph \(\mathcal{G}\). The sensitivity between nodes \(i,j\in\mathcal{V}\) after \(\ell\) iterations of updates can be computed as \(\big{\|}\frac{\partial\mathbf{x}_{i}^{\ell}}{\partial\mathbf{x}_{j}^{0}} \big{\|}_{2}\), where a smaller value indicates lower sensitivity and thus higher degree of oversquashing. In the simple case of isotropic diffusion \(\mathcal{A}_{\mathcal{G}}(\mathbf{X}^{\ell})=\sigma(\mathbf{A}\mathbf{X}^{ \ell}\mathbf{W})\), one can show \(\big{\|}\frac{\partial\mathbf{x}_{i}^{\ell}}{\partial\mathbf{x}_{j}^{0}} \big{\|}_{2}\leq c^{\ell}(\mathbf{A}^{\ell})_{i,j}\), where \(\mathbf{A}^{\ell}\) denotes the \(\ell\)-th matrix power of \(\mathbf{A}\) and \(c>0\) is a Lipschitz constant that depends on the activation function and \(\mathbf{W}\). It can be readily noticed that the graph structure encoded in \(\mathbf{A}\) critically determines the severity of oversquashing.

Below are strategies employed in GNN system that avoids oversquashing.

1. _Graph topology modification_: One natural strategy in addressing the oversquashing is to enhance communications by rewiring the graph structure. GRAND [15] and BLEND [16] dynamically modify the connectivity according to the updated node features, which amounts to changing the spatial discretization over time. DIFFormer [98], on the other hand, adopts a fully connected graph and thus avoids oversquashing by allowing message passing between all pairs of nodes.
2. _Non-local dynamics_: Alternatively, non-local dynamics constructs dense graph communication function and thus mitigates oversquashing by design. This can be achieved for example by taking the fractional power of message passing matrix as in FLODE [68], leveraging the dense quantum diffusion kernel in QDC [67], learning the timestep of heat kernel in TIDE [7], and capturing the entire path of diffusion as in G2TN [89].
3. _Multi-scale and spectral filtering_: Multi-scale diffusion in GradFUFG [46] has the potential to mitigate oversquashing by simultaneously accounting for information at different frequencies. Further, properly controlling the magnitude of high-pass versus low-pass filtering coefficients as in MHKG [81] (in the HFD setting) can provably improve the upper bound of the node sensitivity and thus oversquashing is alleviated.

Despite the success of aforementioned techniques in addressing oversquashing, a potential downside is the increased complexity of propagating information at each update, which is often associated with the more dense message passing matrix. We also remark that some methods, including GIND [18] and A-DGN [41] are able to capture long-range dependencies without suffering from oversmoothing. In particular, the driving mechanisms for such purpose are the implicit diffusion in GIND and anti-symmetric weight in A-DGN where
the former can be viewed as an infinite-depth GNN and the latter allows to build deep GNNs with constant sensitivity to avoid oversmoothing.

Stability, gradient vanishing and explosion.Lastly, we highlight several other potential pitfalls of GNN dynamics, that are often less explored in the literature compared to the previous issues, such as robustness and stability, and gradient vanishing and exploding during training.

Stability and robustness against adversarial attacks is a critical factor in assessing the performance of GNNs. For continuous dynamics, [83] verify that the graph neural diffusion, like GRAND [15] and BLEND [16], are empirically more stable to graph topology perturbation compared to other discrete variants, which follows from the derived theoretical results based on heat kernel. [83] also prove that, due to the row-stochastic normalization of the diffusivity matrix in graph neural diffusion, the stability against node perturbation is guaranteed. In [41], the robustness to the change in initial conditions is explicitly enforced with the required anti-symmetry (due to the zero real parts of the Jacobian eigenvalues). In [102], the conservative stability offered by Hamiltonian dynamics has shown enhanced robustness to adversarial attacks. In addition, [97] improve the generalization ability of GNNs with respect to topological distribution shifts, via local diffusion and global attention.

Another flight of deep neural networks is gradient vanishing and exploding, which refers to the situation when gradient exponentially converges to zero or infinity as a result of increasing depth. Because classic graph neural networks are often designed to be shallow as in GCN [58] and GAT [91], in order to avoid oversmoothing, gradient vanishing or explosion have received less attention in the GNN community. However, in the continuous regimes, especially when depth increases, these problems can emerge and severely escalate the training difficulty. In [78], it has been verified that GraphCON is able to mitigate the the vanishing and exploding gradient problems because the gradient of GraphCON is upper bounded at most quadratically in depth and is shown to be independent of the depth in terms of decaying behaviour. Due to the anti-symmetric channel mixing in A-DGN [41], the magnitude of the gradient stays constant during backpropagation and hence avoids the problem.

## 6 Open research questions

The recent success of graph neural dynamics has presented numerous opportunities for future explorations. This section summarizes several exciting research directions that remain open.

Dynamics with high-order graph structures and spatial derivatives.Most existing continuous GNNs are limited to evolution over nodes, except for [42] that considers higher-order cliques. Meanwhile graphs often encode more intricate topology where higher-order substructures, such as paths, triangles and cycles are crucial for downstream tasks [85]. For example, aromatic rings are cycle structures commonly appearing in molecule, which determines various chemical properties, such as solubility and acidity. Thus designing a dynamics aware of such local substructures is beneficial. In addition, current dynamics often rely on a coupling function that involves only the first-order spatial derivative (i.e., the gradient). It is thus interesting to investigate how to properly define higher-order spatial derivatives (e.g., the Hessian) and incorporate into the dynamics formulation.

Explore the potential of graph neural dynamics for other applications.Existing studies proposing continuous GNNs mostly focus on the node-level prediction tasks while the continuous formalism presents a
general framework for modelling not limited to the evolution of nodes, but also edges, communities and even the entire graph. It is thus rewarding to adapt graph neural dynamics for other types of applications, such as link and graph-level prediction, anomaly detection, time series forecasting. For example, [28] leverage both linear and anisotropic diffusion for molecule property prediction; [9] utilize graph heat and wave equation for graph topology recovery; [30] demonstrate the promise of using advection process for spatial temporal graph learning.

Expressivity of graph neural dynamics.The unprecedented success of GNNs have propelled researchers to study their expressive power in distinguishing non-isomorphic graphs [101, 70], counting subgraph structures [20], etc. While the continuous formalism present a framework for interpreting and designing GNNs, there have been few works that characterize the expressivity of graph neural dynamics, especially on graph and subgraph levels. There also lacks theoretical understanding on how the choice of numerical integrators affects the performance of the neural dynamics. In addition, an equally fruitful direction is to explore whether the theory of dynamical systems, such as energy conservation and reversibility can be leveraged to design theoretically more powerful graph neural networks.

Continuous formulation for other graph types.Most of the continuous GNNs focus primarily on dynamics over static, undirected, homogeneous graphs. Nevertheless, other graph types have also witnessed wide applicability, including _signed_ or _directed graphs_ (where edges can be negative or directed) [24, 51, 69, 87], _heterogeneous graphs_ (where nodes and edges have multiple types) [93, 53], _geometric graphs_ (where nodes or edges respect geometric constraints) [11], _spatial-temporal graphs_ (where graph topology also evolves in time) [56, 55]. The continuous formulation of GNNs for these more complex graph types could be beneficial for both enhancing the understanding and representation power of existing GNNs. However, Such generalization requires nontrivial efforts. For example, the dynamics should preserve symmetries for geometric graphs and account for both evolution of graph topology and signals for spatial-temporal graphs.

## 7 Conclusion

In this work, we conduct a comprehensive review of recent developments on continuous dynamics informed GNNs, an increasingly growing field that marries the theory of dynamical system and differential equation with graph representation learning. In particular, we provide mathematical formulations for a diverse range of graph neural dynamics, and show how the common flights of GNNs can be alleviated through the continuous formulation. We highlight several open challenges and fruitful research directions that warrant further exploration. We hope this survey brings attention to the potential of the continuous formalism of GNNs, and severs as a starting point for future endeavors in harnessing classic theories of continuous dynamics for enhancing the explainability and expressivity of GNN designs.
