\section{Introduction}

The majority of graph neural networks (GNNs) are based on the message passing paradigm \citep{Gilmer2017}, wherein node features are learned by means of a non-linear propagation on the graph. 
Multiple recent works have pointed to the limitations of the message passing approach. These include; limited expressive power~\citep{shervashidze2011weisfeiler,xu2018how,bouritsas2020improving,bodnar2021weisfeiler}, the related oversmoothing problem~\citep{NT2019, Oono2019} and the bottleneck phenomena \citep{Alon2020,wu2020graph}, which render such approaches inefficient, especially in deep GNNs. 
Multiple alternatives have been proposed, among which are higher-order methods \citep{maron2019provably,bodnar2021weisfeiler} and  decoupling the propagation and input graphs by modifying the topology, often referred to as {\em graph rewiring}. Topological modifications can take different forms such as graph sampling \citep{Hamilton2017}, $k$NN \citep{Klicpera2019}, using the complete graph \citep{Vaswani2017,Alon2020}, latent graph learning \citep{wang2019dynamic,kazi2020differentiable}, or multi-hop filters \citep{wu2019simplifying,rossi2020sign}. 
However, there is no agreement in the literature on when and how to modify the graph, and a single principled framework
for doing so. 



A somewhat underappreciated fact is that 
GNNs are intimately related to diffusion equations~\citep{Chamberlain2021a}, a connection that was exploited in the early work of Scarselli et al. \citep{Scarselli2008}. 
Diffusion PDEs have been historically important in computer graphics \citep{sun2009concise,bronstein2010scale,litman2013learning,patane2016star}, computer vision \citep{caselles1997geodesic,chan2001active,bertalmio2000image}, and image processing \citep{perona1990scale,sochen1998general,weickert1998anisotropic,tomasi1998bilateral,durand2002fast,buades2005non}, where they created an entire trend of variational and PDE-based approaches. 
In machine learning and data science, diffusion equations underpin such popular manifold learning methods as eigenmaps \citep{belkin2003laplacian} and diffusion maps \citep{coifman2006diffusion}, as well as the family of PageRank algorithms \citep{page1999pagerank,chakrabarti2007dynamic}. 
In deep learning, differential equations are used as models of neural networks
\citep{Chamberlain2021a, Chen2018a,Dupont2019,Xhonneux2020,Queiruga2020,Zhuang2020a} 
and for physics-informed learning \citep{Raissi2017a,cranmer2019learning,Sanchez-Gonzalez2019,cranmer2020lagrangian,shlomi2020graph,Li2020c}. 


\paragraph{Main contributions}
In this paper, we propose a novel class of GNNs  based on the discretised non-Euclidean diffusion PDE in joint positional and feature space, inspired by the {\em Beltrami flow} \citep{sochen1998general} used two decades ago in the image processing literature for edge-preserving image denoising. 
We show that the discretisation of the spatial component of the Beltrami flow offers a principled view on positional encoding and graph rewiring, whereas the discretisation of the temporal component can replace GNN layers with more flexible adaptive numerical schemes. 
Based on this model, we introduce Beltrami Neural Diffusion (BLEND) that generalises a broad range of GNN architectures and shows state-of-the-art performance on many popular benchmarks.   
In a broader perspective, our approach explores new tools from PDEs and differential geometry that are less well known in the graph ML community. 


\section{Background} \label{sec:background}



\paragraph*{Beltrami flow}
Kimmel et al. \citep{kimmel1997high,sochen1998general,kimmel2000images} considered images as 2-manifolds (parametric surfaces) $(\Sigma, g)$ embedded in some larger ambient space as $\mathbf{z}(\mathbf{u}) = (\mathbf{u},\alpha\mathbf{x}(\mathbf{u})) \subseteq \mathbb{R}^{d + 2}$ where $\alpha \geq 0$ is a scaling factor, $\mathbf{u} = (u_1, u_2)$ are the 2D {\em positional coordinates} of the pixels, and $\mathbf{x}$ are the $d$-dimensional {\em colour} or {\em feature coordinates} (with $d=1$ or $3$ for grayscale or RGB images, or $d=k^2$ when using $k\times k$ patches as features \citep{buades2005non}).
In these works, the image is evolved along the gradient flow of a functional $S[\mathbf{z},g]$ called the {\em Polyakov action
\citep{polyakov1981quantum}}, which roughly measures the smoothness of the embedding\footnote{Explicitly, by minimising the functional with respect to the embedding one finds Euler-Lagrange (EL) equations that can be used to dictate the evolution process of the embedding itself.}. For images embedded in Euclidean space with the functional $S$ minimised with respect to \emph{both} the embedding $\mathbf{z}$ and the metric $g$, one obtains the following PDE:
\vspace{-0.5mm}
\begin{equation}
\frac{\partial\mathbf{z}(\mathbf{u},t)}{\partial t} = \Delta_\mathbf{G} \mathbf{z}(\mathbf{u},t), \quad\quad \mathbf{z}(\mathbf{u},0) = \mathbf{z}(\mathbf{u}), \quad t\geq 0,  
\label{eqn:beltrami0}
\end{equation}
and boundary conditions as appropriate. 
Here $\Delta_\mathbf{G}$ is the {\em Laplace-Beltrami operator}, the Laplacian operator induced on $\Sigma$ by the Euclidean space we embed the image into. Namely, the embedding of the manifold allows us to \emph{pull-back} the Euclidean distance structure on the image: the distance between two nearby points $\mathbf{u}$ and $\mathbf{u}+\mathrm{d}\mathbf{u}$ is given by 
\begin{align}
\mathrm{d}\ell^2 = \mathrm{d}\mathbf{u}^\top \mathbf{G}(\mathbf{u}) \mathrm{d}\mathbf{u} = \mathrm{d}u_1^2 + \mathrm{d}u_2^2 + \alpha^2 \sum_{i=1}^d \mathrm{d}x_i^2,
\label{eq:metric}
\end{align}
where $\mathbf{G} = \mathbf{I} + \alpha^2 (\nabla_{\mathbf{u}}\mathbf{x}(\mathbf{u}))^\top \nabla_{\mathbf{u}}\mathbf{x}(\mathbf{u})$ is a $2\times 2$ matrix called the {\em Riemannian metric}.
The fact that the distance is a combination of the positional component (distance between pixels in the plane, $\| \mathbf{u} - \mathbf{u}' \|$) and the colour component (distance between the colours of the pixels, $\| \mathbf{x}(\mathbf{u}) - \mathbf{x}(\mathbf{u}') \|$) is crucial as it allows edge-preserving image diffusion.


When dealing with images, the evolution of the first two components of $(z_1, z_2) = \mathbf{u}$ is a nuisance amounting to the reparametrisation of the manifold and can be ignored. For grayscale images (the case when $d=1$ and $\mathbf{z} = (u_1, u_2, x)$), this is done by projection along the dimension $z_3$, in which case the Beltrami flow takes the form of an inhomogeneous diffusion equation of $x$,  
\begin{equation}
\frac{\partial x(\mathbf{u},t)}{\partial t} = \frac{1}{\sqrt{\mathrm{det}\mathbf{G}(\mathbf{u},t) }}\mathrm{div}\left( \frac{\nabla x(\mathbf{u},t) }{\sqrt{ \mathrm{det}\mathbf{G}(\mathbf{u},t) }} \right) \quad\quad t\geq 0.\vspace{-1mm}
\label{eqn:beltrami}
\end{equation}
\begin{wrapfigure}{r}{0.42\textwidth}\vspace{0mm}
\centering
    \includegraphics[width=0.35\textwidth]{beltrami1.png}\vspace{-1mm}
  \caption{\label{fig:beltrami}Two interpretations of the Beltrami flow: position-dependent bilateral kernel (top) and a Gaussian passed on the manifold (bottom). \vspace{-15mm}
  }
\end{wrapfigure}
The {\em diffusivity} \vspace{-2mm}
\begin{equation}
    a = \frac{1}{\sqrt{\mathrm{det}\mathbf{G}}} = \frac{1}{\sqrt{1 + \alpha^2 \| \nabla x \|^2}}
    \label{eqn:diffusivity}
    \end{equation}
    determining the speed of diffusion at each point, can be interpreted as an {\em edge indicator}:  diffusion is weak across edges where $\| \nabla x \| \gg 1$. 
The result is an adaptive diffusion \citep{perona1990scale} popular in image processing due to its ability to denoise images while preserving their edges.
For cases with $d>1$ (multiple colour channels), equation~(\ref{eqn:beltrami}) is applied to each channel separately; however, the metric $\mathbf{G}$ couples the channels, which results in their gradients becoming aligned \citep{kimmel2012numerical}.



\paragraph{Special cases} 
In the limit case $\alpha=0$, equation~(\ref{eqn:beltrami}) becomes the simple homogeneous isotropic diffusion $\tfrac{\partial }{\partial t}{x} = \mathrm{div}(\nabla x) = \Delta x$, where $\Delta = \tfrac{\partial^2}{\partial u_1^2} + \tfrac{\partial^2}{\partial u_2^2}$ is the standard Euclidean Laplacian operator. The solution is given in closed form as the convolution of the initial image and a Gaussian kernel with time-dependent variance,
\begin{align}
x(\mathbf{u},t) = x(\mathbf{u},0) \star \frac{1}{(4\pi t)^{d/2}} e^{-\| \mathbf{u}\|^2/4t} 
\end{align}
and can be considered a simple linear low-pass filtering.  In the limit $t\rightarrow \infty$, the image becomes constant and equal to the average colour.\footnote{Assuming appropriate boundary conditions.}


Another interpretation of the Beltrami flow is passing a Gaussian {\em on the manifold} (see Figure~\ref{fig:beltrami}, bottom), which can locally be expressed as non-linear filtering with the {\em bilateral kernel} \citep{tomasi1998bilateral} dependent on the joint positional and colour distance (Figure~\ref{fig:beltrami}, top),
\begin{equation}
x(\mathbf{u},t) = \frac{1}{(4\pi t)^{d/2}}\int_{\mathbb{R}^2} x(\mathbf{v},0) e^{-\|\mathbf{u}-\mathbf{v}\|^2/4t} e^{- \alpha^2 \|\mathbf{x}(\mathbf{u,0})-\mathbf{x}(\mathbf{v,0})\|^2/4t} \mathrm{d}\mathbf{v}. 
\label{eqn:bilateral}
\end{equation}
For $\alpha =0$, the bilateral filter~(\ref{eqn:bilateral}) reduces to a simple convolution with a time-dependent Gaussian.




\section{Discrete Beltrami flow on graphs}

We now develop the analogy of Beltrami flow for graphs. We consider a graph to be a {\em discretisation} of a continuous structure (manifold), and show that the evolution of the feature coordinates in time amounts to message passing layers in GNNs, whereas the evolution of the positional coordinates amounts to graph rewiring, which is used in some GNN architectures.  

\subsection{Graph Beltrami flow} 
Let $\mathcal{G}=(\mathcal{V}=\{1,\hdots, n\},\mathcal{E})$ be an undirected graph, where $\mathcal{V}$ and $\mathcal{E}$ denote node and edge sets, respectively. We further assume node-wise $d$-dimensional features
 $\mathbf{x}_i \in \mathbb{R}^d$ for $i=1,\hdots, n$.  
Denote by $\mathbf{z}_i = (\mathbf{u}_i, \alpha\mathbf{x}_i)$ the embedding of the graph in a joint space $\mathcal{C} \times \mathbb{R}^{d}$, where $\mathcal{C}$ is a $d'$-dimensional space with a metric $d_\mathcal{C}$ representing the node coordinates (for simplicity, we will assume $\mathcal{C} = \mathbb{R}^{d'}$ unless otherwise stated).  
We refer to $\mathbf{u}_i$ $\mathbf{x}_i$ and $\mathbf{z}_i$  as the {\em positional}, {\em feature} and {\em joint} coordinates of node $i$, respectively, and arrange them into the matrices $\mathbf{U}$, $\mathbf{X}$, and $\mathbf{Z}$,  of sizes $n\times d'$, $n\times d$, and $n\times (d'+d)$.  

For images, Beltrami flow amounts to evolving the embedding $\mathbf{z}$ along $\text{div}(a(\mathbf{z})\nabla \mathbf{z})$, with $a$ a diffusivity map.\footnote{Note that in equation \eqref{eqn:beltrami} the diffusivity function $a$ coincides with $(\text{det}(\mathbf{G}(\mathbf{u},t)))^{-\frac{1}{2}}$. Similarly to \cite[Section 4.2]{sochen1998general} we have neglected the extra term $1/a$ appearing in \eqref{eqn:beltrami}.}  Thus, we consider the {\em graph Beltrami flow} to be the discrete diffusion equation
\begin{equation}
    \frac{\partial \mathbf{z}_i(t) }{\partial t} = \hspace{-4mm}\sum_{j : (i,j) \in \mathcal{E}' } \hspace{-4mm} a(\mathbf{z}_i(t), \mathbf{z}_j(t) ) (\mathbf{z}_j(t) - \mathbf{z}_i(t) )\quad\quad \mathbf{z}_i(0) = \mathbf{z}_i;
    \quad i=1, \hdots, n;  \quad t\geq 0.
    \label{eqn:graph_beltrami}
\end{equation}
We motivate our definition as follows:
$\mathbf{g}_{ij} = \mathbf{z}_j - \mathbf{z}_i$  and $\mathbf{d}_i =  \sum_{j : (i,j) \in \mathcal{E} } \mathbf{g}_{ij} $ are the discrete analogies of the gradient $\nabla \mathbf{z}$ and divergence $\mathrm{div}(\mathbf{g})$, both with respect to a graph $(\mathcal{V},\mathcal{E}')$ that can be interpreted as the numerical stencil for the discretisation of the continuous Laplace-Beltrami operator in~(\ref{eqn:beltrami}). Note that $\mathcal{E}'$ can be different from the input $\mathcal{E}$ (referred to as `rewiring'). As discussed in Section~\ref{sec:gnns}, most GNNs use $\mathcal{E}'=\mathcal{E}$ (input graph is used for diffusion, no rewiring). Alternatively, 
the positional coordinates of the nodes can be used to define a new graph topology either with $\mathcal{E}(\mathbf{U}) = \{ (i,j) : d_\mathcal{C}(\mathbf{u}_i,\mathbf{u}_j) < r \}$, for some radius $r>0$, or using $k$ nearest neighbours.  
This new rewiring is precomputed using the input positional coordinates (i.e., $\mathcal{E}' = \mathcal{E}(\mathbf{U}(0))$) or updated throughout the diffusion (i.e., $\mathcal{E}'(t) = \mathcal{E}(\mathbf{U}(t))$). Therefore, \eqref{eqn:graph_beltrami} can be compactly rewritten as
\[
\frac{\partial \mathbf{z}_i(t) }{\partial t} = \text{div}\left( a(\mathbf{z}(t))\nabla\mathbf{z}_{i}(t)\right). 
\]
\noindent The function $a$ is the diffusivity controlling the diffusion strength between nodes $i$ and $j$ and is assumed to be normalised: $ \sum_{j:(i,j) \in \mathcal{E}'}  a(\mathbf{z}_i, \mathbf{z}_j) = 1$. The dependence of the diffusivity on the embedding $\mathbf{z}$ matches the smooth PDE analysed in e.g. \cite[Section 4.2]{sochen1998general} and is consistent with the form of attention mechanism used in e.g. \citep{Velickovic2018,Vaswani2017}. 
In matrix-form, we can also rewrite \eqref{eqn:graph_beltrami} as
\begin{align}
\begin{pmatrix}
\frac{\partial }{\partial t} \mathbf{U}(t), 
\frac{\partial  }{\partial t} \mathbf{X}(t)
\end{pmatrix}
&=
    (\mathbf{A}(\mathbf{U}(t), \mathbf{X}(t)) - \mathbf{I}) \begin{pmatrix}
\mathbf{U}(t), 
\mathbf{X}(t)
\label{eqn:graph_beltrami_m}
\end{pmatrix} \\
    \,\, \mathbf{U}(0) &=  \mathbf{U}; \,\,  \mathbf{X}(0) = \alpha \mathbf{X};  \,\, t\geq 0, \notag
\end{align}
where we emphasise the evolution of both the positional and feature  components, coupled through the matrix-valued function $\mathbf{A}$  
$$
a_{ij}(t) = \left\{ 
\begin{array}{ll}
 a((\mathbf{u}_i(t), \mathbf{x}_i(t)),(\mathbf{u}_j(t), \mathbf{x}_j(t))) & (i,j) \in \mathcal{E}(\mathbf{U}(t))  \\
0 & \text{else}.
\end{array}
\right.
$$
representing the diffusivity. 
The graph Beltrami flow produces an evolution of the joint positional and feature coordinates, $\mathbf{Z}(t) = (\mathbf{U}(t), \mathbf{X}(t))$. 
In Section~\ref{sec:gnns} we will show how the evolution of the feature coordinates $\mathbf{X}(t)$ results in feature diffusion or message passing on the graph, the core of GNNs.
As noted in Section~\ref{sec:background}, in the smooth case the Beltrami flow is obtained as gradient flow of an energy functional when minimised with respect to \emph{both} the embedding and the metric on the surface (an image). When the embedding takes values in the Euclidean space, this leads to equations of the form \eqref{eqn:beltrami} with no channel-mixing and an exact form of the diffusivity determined by the pull-back $\mathbf{G}$ of the Euclidean metric. To further motivate our approach, it is tempting to investigate whether a similar conclusion can be attained here. Although in the discrete case the operation of pull-back is not well-defined, we are able to derive that the gradient flow of a \emph{modified} graph Dirichlet energy gives rise to an equation of the form \eqref{eqn:graph_beltrami}. We note though that the gradient flow does not recover the exact form of the diffusivity implemented in this paper. This is not a limitation of the theory and should be expected: by requiring the gradient flow to avoid channel-mixing and imitate the image analogy in \citep{sochen1998general} and by inducing a \emph{discrete pull-back} condition, we are imposing constraints on the problem. We leave the theoretical implications for future work and refer to Appendix~\ref{app:theory} for a more thorough discussion, including definitions and proofs.

\begin{theorem}\label{thm:gradientflow}
Under structural assumptions on the diffusivity, graph Beltrami flow~\eqref{eqn:graph_beltrami} is the gradient flow of the discrete Polyakov functional. 
\end{theorem}

\subsection{Numerical solvers} 
\paragraph{Explicit vs implicit schemes}
Equation~(\ref{eqn:graph_beltrami}) is solved numerically, which in the simplest case is done by replacing the continuous time 
derivative $\tfrac{\partial}{\partial t}$ with forward time difference: 
\begin{equation}
\frac{\mathbf{z}^{(k+1)}_ i - \mathbf{z}^{(k)}_i}{\tau} = \hspace{-3mm}\sum_{j: (i,j) \in \mathcal{E}(\mathbf{U}^{(k)}) }
\hspace{-3mm}
a\left(\mathbf{z}^{(k)}_i, \mathbf{z}^{(k)}_j\right) (\mathbf{z}^{(k)}_j - \mathbf{z}^{(k)}_i).
\label{eq:g_diff}
\end{equation}
Here $k$ denotes the discrete time index (iteration) and $\tau$ is the time step (discretisation parameter).  
Rewriting~(\ref{eq:g_diff}) compactly in matrix-vector form with $\tau=1$
leads to the {\em explicit Euler scheme}:
\begin{equation}
\mathbf{Z}^{(k+1)}  = 
(\mathbf{A}^{(k)} - \mathbf{I}) \mathbf{Z}^{(k)} 
= \mathbf{Q}^{(k)} \mathbf{Z}^{(k)},
\label{eq:g_diff_m}
\end{equation}
where $a_{ij}^{(k)} = a(\mathbf{z}_i^{(k)},\mathbf{z}_j^{(k)})$ and the matrix $\mathbf{Q}^{(k)}$ (diffusion operator) is given by 
$$
q_{ij}^{(k)} = \left\{ 
\begin{array}{ll}
 1-\tau \hspace{-5pt}\displaystyle\sum_{l: (i,l)\in \mathcal{E}}\hspace{-5pt} a_{il}^{(k)} & i=j  \\
\tau a_{ij}^{(k)} & (i,j) \in \mathcal{E}(\mathbf{U}^{(k)})  \\
0 & \text{otherwise}
\end{array}
\right.
$$
The solution to the diffusion equation is computed by applying scheme~(\ref{eq:g_diff_m}) multiple times in sequence, starting from some initial 
$\mathbf{Z}^{(0)}$. 
It is `explicit' because the update $\mathbf{Z}^{(k+1)}$ is done  directly by the application of the diffusion operator $\mathbf{Q}^{(k)}$ on $\mathbf{Z}^{(k)}$ (as opposed to {\em implicit schemes} of the form $\mathbf{Z}^{(k)} =\mathbf{Q}^{(k)}\mathbf{Z}^{(k+1)}$ arising from backward time differences that require inversion of the diffusion operator \citep{Weickert1998}). 


\paragraph{Multi-step and adaptive schemes} 
Higher-order approximation of temporal derivatives amount to using intermediate fractional steps, which are then linearly combined.  Runge-Kutta (RK) \citep{runge1895numerische,kutta1901beitrag}, ubiquitously used in numerical analysis, is a classical family of explicit numerical schemes, including Euler as a particular case. 
The Dormand-Prince (DOPRI) \citep{dormand1980family} is an RK method based on fifth and fourth-order approximations, the difference between which is used as an error estimate guiding the time step size \citep{shampine1986some}. 

\paragraph{Adaptive spatial discretisation and rewiring} 
Many numerical PDE solvers also employ adaptive spatial discretisation. The choice of the stencil (mesh) for spatial derivatives is done based on the character of the solution at these points; in the simulation of phenomena such as shock waves it is often desired to use denser sampling in the respective regions of the domain, which can change in time. A class of techniques for adaptive rewiring of the spatial derivatives are known as Moving Mesh (MM) methods \citep{hawken1991review}. Interpreting the graph $\mathcal{E}'$ in~\eqref{eqn:graph_beltrami} as the numerical stencil for the discretisation of the continuous Laplace-Beltrami operator in~\eqref{eqn:beltrami}, we can regard rewiring
as a form of MM. 



\subsection{Relation to graph neural networks}
\label{sec:gnns}

Equation~(\ref{eq:g_diff}) has the structure of many GNN architectures of the `attentional' type \citep{gdlbook}, where the discrete time index $k$ corresponds to a (convolutional or attentional) layer of the GNN and multiple diffusion iterations amount to a deep GNN. 
In the diffusion formalism, the time parameter $t$ acts as a continuous analogy of the layers, in the spirit of neural differential equations \citep{Chen2018a}. Typical GNNs amount to explicit single-step (Euler) discretisation schemes, whereas our  continuous interpretation can exploit more efficient numerical schemes. 






\paragraph{GNNs as instances of graph Beltrami flow}
The graph Beltrami framework leads to a family of graph neural networks  that generalise many popular architectures (see Table~\ref{tab:comparison}). 
For example, GAT \citep{Velickovic2018} can be obtained as a particular setting of our framework where the input graph is fixed $(\mathcal{E}' = \mathcal{E})$ and only the feature coordinates $\mathbf{X}$ are evolved. Equation~(\ref{eq:g_diff_m}) in this case becomes
\begin{equation}
\mathbf{x}^{(k+1)}_ i = \mathbf{x}^{(k)}_i + \tau \hspace{-3mm}\sum_{j: (i,j) \in \mathcal{E} }
\hspace{-3mm}
a\left(\mathbf{x}^{(k)}_i, \mathbf{x}^{(k)}_j\right) (\mathbf{x}^{(k)}_j - \mathbf{x}^{(k)}_i)
\label{eq:g_diff_gat}
\end{equation}
and corresponds to the update formula of GAT with a residual connection and the assumption of no non-linearity between the layers.
The role of the diffusivity is played by a learnable parametric {\em attention} function, which is generally {\em time-dependent}: $a(\mathbf{z}_i^{(k)}, \mathbf{z}_j^{(k)},k)$. This results in separate attention parameters per layer $k$, which can be learned independently. 
Our intentionally simplistic choice of a {\em time-independent} attention function amounts to {\em parameter sharing} across layers. We will show in Section~\ref{sec:node_classification} that this leads to a smaller model that is less likely to overfit. 


Another popular architecture MoNet \citep{Monti2016} uses linear diffusion of the features with weights dependent on the structure of the graph expressed as `pseudo-coordinates', which can be cast as attention of the form $a(\mathbf{u}_i, \mathbf{u}_j)$.
Transformers \citep{vaswani2017attention} can be interpreted as feature diffusion on a fixed complete graph with $\mathcal{E}' = \mathcal{V}\times \mathcal{V}$ \citep{gdlbook}. Positional encoding (used in Transformers as well as in several recent GNN architectures \citep{bouritsas2020improving,dwivedi2020generalization}) amounts to attention dependent on both $\mathbf{X}$ and $\mathbf{U}$, which allows the diffusion to adapt to the local structure of the graph; importantly, the positional coordinates $\mathbf{U}$ are {\em precomputed}. 
Similarly, DeepSets \citep{zaheer2017deep} and PointNet \citep{qi2017pointnet} architectures can be interpreted as GNNs applied on a graph with no edges ($\mathcal{E}' = \emptyset$), where each node is treated independently of the others \citep{gdlbook}. 
DIGL \citep{Klicpera2019} performs graph rewiring as a pre-processing step 
using personalised page rank as positional coordinates, 
which are then fixed and not evolved. 
In the point cloud methods, DGCNN \citep{wang2019dynamic} and DGM \citep{kazi2020differentiable}, the graph is constructed based on the {\em feature coordinates} $\mathbf{X}$ and rewired adaptively (in our notation, $\mathcal{E}' = \mathcal{E}(\mathbf{X}(t))$). 







 \begin{table}[!h]
 \resizebox{0.95\textwidth}{!}{
 \begin{tabular}{l llll}
 \textbf{Method} & \textbf{Evolution}   & \textbf{Diffusivity}   & \textbf{Graph} $(\mathcal{V}, \mathcal{E}')$   & \textbf{Discretisation} \\
 \hline
 ChebNet   & Features $\mathbf{X}$    & Fixed $a_{ij}$   & Fixed $\mathcal{E}$   & Explicit fixed step \\
 GAT   & Features $\mathbf{X}$    &  $a(\mathbf{x}_i,\mathbf{x}_j)$   & Fixed $\mathcal{E}$   & Explicit fixed step \\
 %
MoNet   & Features $\mathbf{X}$    &  $a(\mathbf{u}_i,\mathbf{u}_j)$   & Fixed $\mathcal{E}$   & Explicit fixed step \\
 %
 Transformer & Features $\mathbf{X}$    & $a((\mathbf{u}_i, \mathbf{x}_i),(\mathbf{u}_j,\mathbf{x}_j))$   & Fixed $\mathcal{E}=\mathcal{V}\times \mathcal{V}$   & Explicit fixed step \\
 %
 DeepSet/PointNet & Features $\mathbf{X}$    & $a(\mathbf{x}_i)$   & Fixed $\mathcal{E}=\emptyset$   & Explicit fixed step \\ 
 %
  DIGL$^*$   &  Features $\mathbf{X}$    & $a(\mathbf{x}_i,\mathbf{x}_j)$   & Fixed $\mathcal{E}(\mathbf{U})$   & Explicit fixed step \\
 %
 DGCNN/DGM$^*$   &  Features $\mathbf{X}$    & $a(\mathbf{x}_i,\mathbf{x}_j)$   & Adaptive $\mathcal{E}(\mathbf{X})$   & Explicit fixed step \\
 %
 \hline
 {\bf Beltrami}   & Positions $\mathbf{U}$     & $a((\mathbf{u}_i, \mathbf{x}_i),(\mathbf{u}_j,\mathbf{x}_j))$   & Adaptive $\mathcal{E}(\mathbf{U})$   & Explicit adaptive step /
       \\
   & +Features $\mathbf{X}$    &    &    &  Implicit \\
 \hline
 \end{tabular}
 }
  \vspace{1mm}
 \caption{GNN architectures interpreted as particular instances of our framework. $^*$Attentional variant. }
 \vspace{-5mm}
 \label{tab:comparison}
 \end{table}




\paragraph{Graph rewiring}
Multiple authors have recently argued in favor of decoupling the input graph from the graph used for diffusion. Such rewiring can take the form of graph sampling \citep{Hamilton2017} to address scalability issues, data denoising \citep{Klicpera2019}, removal of information bottlenecks \citep{Alon2020}, or larger multi-hop filters \citep{rossi2020sign}. The graph construction can also be made differentiable and a task-specific rewiring can be learned \citep{wang2019dynamic,kazi2020differentiable}. 
The statement of Klicpera et al. \citep{Klicpera2019} that `diffusion improves graph learning', leading to the eponymous paradigm (DIGL), can be understood as a form of diffusion on the graph connectivity independent of the features. Specifically, the authors used as node positional encoding the Personalised PageRank (PPR), which can be interpreted as the steady-state of a diffusion process 
\begin{align}
\mathbf{U}_\mathrm{PPR} = \sum_{k \geq 0}(1-\beta)\beta^{k}\boldsymbol{\Delta}_{\mathrm{RW}}^{k}
=
(1-\beta)(\mathbf{I}-\beta\boldsymbol{\Delta}_{\mathrm{RW}} )^{-1},
    \quad 0<\beta<1,
    \label{eqn:PPR_steady}
\end{align}
where $\boldsymbol{\Delta}_{\mathrm{RW}}$ is the random walk graph Laplacian and $\beta \in (0,1)$ is a parameter such that $1-\beta$ represents the restart probability. The resulting positional encoding of dimension $d=n$ can be used to rewire the graph by $k$NN sampling, which corresponds to using $\mathcal{E}' = \mathcal{E}(\mathbf{U}_\mathrm{PPR})$ in our framework.  



\paragraph{Numerical schemes}
All the aforementioned GNN architectures can be seen as an explicit discretisation of equation~(\ref{eqn:graph_beltrami}) with fixed step size. On the other hand, our continuous diffusion framework offers an additional advantage of employing more efficient numerical schemes with adaptive step size. 
Graph rewiring of the form $\mathcal{E}' = \mathcal{E}(\mathbf{U}(t))$ can be interpreted as adaptive spatial discretisation (moving mesh method). 



\subsection{Extensions}
\begin{wrapfigure}{R}{0.4\textwidth
    \vspace{-19mm}
    \includegraphics[width=0.35\textwidth]{hyperbolic.pdf}\vspace{-2mm}
  \caption{\label{fig:hyperbolic}Graph Beltrami flow with hyperbolic positional coordinates. \vspace{-7mm}
  }
\end{wrapfigure}

\paragraph{Non-Euclidean geometry}
There are multiple theoretical and empirical arguments \citep{liu2019hyperbolic,chami2019hyperbolic} in favor of using hyperbolic spaces to represent real-life `small-world' graphs (in particular, scale-free networks can be obtained as $k$NN graphs in such spaces \citep{boguna2021network}). Our framework allows using a non-Euclidean metric $d_{\mathcal{C}}$ for the positional coordinates $\mathbf{U}$ (Figure~\ref{fig:hyperbolic}). In Section~\ref{sec:experiments} we show that hyperbolic positional encodings allow for a significant reduction in model size with only a marginal degradation of performance.  

\paragraph{Time-dependent diffusivity}
The diffusivity function $a$ which we assumed time-independent and which lead to parameter sharing across layers (i.e., updates of the form $\mathbf{Z}^{(k+1)} = \mathbf{Q}(\mathbf{Z}^{(k)}, \boldsymbol{\theta}) \mathbf{Z}^{(k)}$) can be made time-dependent of the form 
$\mathbf{Z}^{(k+1)} = \mathbf{Q}(\mathbf{Z}^{(k)}, \boldsymbol{\theta}^{(k)}) \mathbf{Z}^{(k)}$, where 
$\boldsymbol{\theta}$ and $\boldsymbol{\theta}^{(k)}$ denote shared and layer-dependent parameters, respectively. 


\paragraph{Onsager diffusion}
As we noted, the Beltrami flow diffuses each channel separately. A more general variant of diffusion allowing for feature mixing is the {\em Onsager diffusion} \citep{onsager1931reciprocal} of the form $\tfrac{\partial}{\partial t}\mathbf{Z}(t) = \mathbf{Q}(\mathbf{Z}(t))  \mathbf{Z}(t) \mathbf{W}(t)$, where 
the matrix-valued function $\mathbf{W}$ acts across the channels.  
GCN \citep{kipf2017} can be regarded a particular setting thereof, with update  of the form $\mathbf{X}^{(k+1)} = \mathbf{A}\mathbf{X}^{(k)}\mathbf{W}$. 


\paragraph{MPNNs}
Finally, we note that the Beltrami flow amounts to linear aggregation with non-linear coefficients, or the `attentional' flavor of GNNs \citep{gdlbook}. The more general message passing flavor \citep{Gilmer2017} is possible using a generic non-linear equation of the form $\tfrac{\partial}{\partial t}\mathbf{Z}(t) = \boldsymbol{\Psi}(\mathbf{Z}(t))$. 




\section{BLEND: Beltrami Neural Diffusion}

Beltrami Neural Diffusion (BLEND) is a novel class of graph neural network architectures derived from the graph Beltrami framework.  
We assume an input graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with $n$ nodes and $d$-dimensional node-wise features represented as a  matrix ${\mathbf{X}}_{\mathrm{in}}$. 
We further assume a $d'$-dimensional  positional encoding ${\mathbf{U}}_{\mathrm{in}}$ of the graph nodes.  
BLEND architectures implement 
a learnable joint diffusion process of $\mathbf{U}$ and $\mathbf{X}$ and runs it for time $T$, 
to produce an output node embeddings $\mathbf{Y}$,
\begin{align
    \mathbf{Z}(0) &= (\phi(\mathbf{U}_{\mathrm{in}}), \psi(\mathbf{X}_{\mathrm{in}}) ) \quad \quad
    \mathbf{Z}(T) = \mathbf{Z}(0) + \int_0^T \frac{\partial\mathbf{Z}(t)}{\partial t} \mathrm{d}t \quad\quad 
    \mathbf{Y} = \xi(\mathbf{Z}(T)),\notag 
\vspace{-6mm}
\end{align}

where 
$\phi$, $\psi$ are learnable positional and feature encoders and $\xi$ is a learnable decoder (possibly changing the output dimensions). Here the $\alpha$ in Equations~\eqref{eq:metric} and \eqref{eqn:graph_beltrami_m} is absorbed by $\psi$ and made learnable. 
$\frac{\partial\mathbf{Z}(t)}{\partial t} $ is given by the graph Beltrami flow equation~\eqref{eqn:graph_beltrami_m}, where the diffusivity function (attention) $a$ is also learnable. The choice of attention function depends on the geometry of the positional encoding and for Euclidean encodings we find the scaled dot product attention~\citep{Vaswani2017} performs well, in which case
\begin{align}
    a(\vec{z}_i,\vec{z}_j) = \softmax \left(\frac{(\mathbf{W}_K \vec{z}_{i} )^\top \mathbf{W}_Q \vec{z}_{j}}{d_k}\right)
    \label{eq:attention}
\end{align}
where $\mathbf{W}_K$ and $\mathbf{W}_Q$ are learned matrices, and $d_k$ is a hyperparameter.  


\section{Experimental results}
\label{sec:experiments}

In this section, we compare the proposed Beltrami framework to popular GNN architectures on standard node classification benchmarks and provide a detailed study of the choice of the positional encoding space. 
Implementation details, including runtimes and hyperparameter tuning are given in Section~\ref{app:implementation_details} and code is available at \url{https://github.com/twitter-research/graph-neural-pde}.




\paragraph{Datasets}
In our experiments, we use the following datasets: 
Cora~\citep{mccallum2000automating}, Citeseer~\citep{sen2008collective}, Pubmed~\citep{namata2012query},
CoauthorCS~\citep{shchur2018pitfalls}, Amazon, Computer, and  Photo~\citep{mcauley2015image}, and OGB-arxiv~\citep{hu2020open}. 
Since many works using the first three datasets rely on the 
Planetoid splits~\citep{yang2016revisiting},
we included them Table~\ref{tab:planetoid_results}, together with 
a more robust evaluation on 
100 random splits with 20 random initialisations ~\citep{shchur2018pitfalls}. In all cases the largest connected component is used and dataset statistics are given in Appendix~\ref{app:datasets}.

\paragraph{Baselines} We compare to the following GNN architectures: GCN~\citep{kipf2017}, GAT~\citep{Velickovic2018}, MoNet~\citep{Monti2016} and GraphSAGE~\citep{Hamilton2017}, and recent ODE-based GNN models: CGNN~\citep{Xhonneux2020}, GDE~\citep{Poli2019}, GODE~\citep{Zhuang2020a}, and two versions of LanczosNet~\citep{liao2019lanczosnet}.
We use two variants of our method: using fixed input graph (BLEND) and using $k$NN graph in the positional coordinates (BLEND-knn).
 



\subsection{Node Classification} \label{sec:node_classification}

\begin{wraptable}{R}{0.5\textwidth}\vspace{-10mm}
\resizebox{0.5\textwidth}{!}{
\begin{tabular}{llll}
\textbf{Method} & \textbf{CORA}       & \textbf{CiteSeer}   & \textbf{PubMed}     \\
\hline
\textit{GCN}             & 81.9$\pm$0.8 & 69.5$\pm$0.9 & 79.0$\pm$0.5 \\
\textit{GAT}             & 82.8$\pm$0.5 & 71.0$\pm$0.6 & 77.0$\pm$1.3 \\
\textit{MoNet}           & 82.2$\pm$0.7 & 70.0$\pm$0.6 & 77.7$\pm$0.6 \\
\textit{GS-maxpool}      & 77.4$\pm$1.0 & 67.0$\pm$1.0 & 76.6$\pm$0.8 \\
\textit{Lanczos}    & 79.5$\pm$1.8 & 66.2$\pm$1.9 & 78.3$\pm$0.3 \\
\textit{AdaLanczos}    & 80.4$\pm$1.1 & 68.7$\pm$1.0 & 78.1$\pm$0.4 \\
\textit{CGNN}     & 81.7$\pm$0.7 & 68.1$\pm$1.2 & 80.2$\pm$0.3 \\
\textit{GDE}    & \textcolor{blue}{\bf{83.8}}$\pm$0.5 & 72.5$\pm$0.5 & 79.9$\pm$0.3 \\
\textit{GODE}    & 83.3$\pm$0.3 & 72.4$\pm$0.6 & 80.1$\pm$0.3 \\
\hline
\textit{\textbf{ BLEND}}           & \textcolor{red}{\bf{84.2}}$\pm$0.6
& \textcolor{blue}{\bf{74.4}}$\pm$0.7 & \textcolor{blue}{\bf{80.7}}$\pm$ 0.7\\
\textit{\textbf{ BLEND-kNN
}}       & 83.1$\pm$0.8 
& \textcolor{red}{\bf{73.3}}$\pm$0.9 & 
\textcolor{red}{{\bf 81.5}}$\pm$0.5
\end{tabular}

}
\caption{Performance (test accuracy$\pm$std) of different GNN models using Planetoid splits.
}

\label{tab:planetoid_results}
\end{wraptable}
In these experiments, we followed the methodology of~\citep{shchur2018pitfalls} using 20 random weight initialisations for datasets with fixed Planetoid splits and 100 random splits otherwise. 
Where available, results from~\citep{shchur2018pitfalls} are reported. 
Hyperparameters with the highest validation accuracy were chosen and results are reported on a test set that is used only once. Hyperparameter search used Ray Tune~\citep{Liaw2018} with a thousand random trials using an asynchronous hyperband scheduler with a grace period and half life of ten epochs. 
The code to reproduce our results is included with the submission and will be released publicly following the review process. 
Experiments ran on  AWS p2.8xlarge machines, each with 8 Tesla V100-SXM2 GPUs.




\paragraph{Implementation details}

For all datasets excepting ogb-arxiv, adaptive explicit Dormand-Prince scheme  was used as the numerical solver; for ogb-arxiv, we used  the Runge-Kutta method. 
For the two smallest datasets (Cora and Citeseer) we performed direct backpropagation through each step of the numerical integrator. For the larger  datasets, to reduce memory complexity, we use Pontryagin's maximum principle to propagate gradients backwards in time~\citep{pontryagin2018mathematical}.
For the larger datasets, kinetic energy and Jacobian regularisation \citep{Finlay2020, Kelly2020} was employed. 
The regularisation ensures the learned dynamics is well-conditioned and easily solvable by a numeric solver, which reduced training time. 
We use constant initialisation for the attention weights, $\mathbf{W}_K,\mathbf{W}_Q$,
so training starts from a well-conditioned system that induces small regularisation penalty terms \citep{Finlay2020}.



The space complexity of BLEND is dominated by evaluating attention~\eqref{eq:attention} over edges and is $\mathcal{O}(|\mathcal{E}'|(d+d'))$ where $\mathcal{E}'$ is the edge set following rewiring and $d$ is dimension of features and $d'$ is the dimension of positional encoding.
The runtime complexity is $\mathcal{O}(|\mathcal{E}'|(d+d'))(E_b + E_f)$, split between the forward and backward pass and can be dominated by either depending on the number of function evaluations ($E_b$, $E_f$).




Tables~\ref{tab:planetoid_results}--\ref{tab:random_splits} summarise the results of our experiments. BLEND outperforms other GNNs in most of the experiments, showing state-of-the-art results on some datasets. 
Another important point to note is that compared GNNs use different sets of parameters per layer, whereas in BLEND, due to our choice of a time-independent attention, parameters are shared. This results in significantly fewer parameters: for comparison, the OGB versions of GCN, SAGE and GAT used in the ogb-arxiv experiment require 143K, 219K and 1.63M parameters respectively, compared to only 70K in BLEND.  















\begin{table*}[!tb]
\resizebox{\textwidth}{!}{
\begin{tabular}{lccccccc}
{\bf Method} & \textbf{CORA}       & \textbf{CiteSeer}   & \textbf{PubMed}  & \textbf{Coauthor CS}       & \textbf{Computer}   & \textbf{Photo} & \textbf{ogb-arxiv}   \\
\textit{GCN}           & 81.5$\pm$1.3 &  71.9 $\pm$1.9 & 77.8$\pm$2.9 & 91.1$\pm$0.5 & 82.6$\pm$2.4 & 91.2$\pm$1.2  & 71.74$\pm$0.29  \\
\textit{GAT}            & 81.8$\pm$1.3 & 71.4$\pm$1.9 &  78.7$\pm$2.3 & 90.5$\pm$0.6 & 78.0$\pm$19 & 85.7$\pm$20  & \textcolor{red}{\bf{73.01}}$\pm$0.19$^*$ \\
\textit{GAT-ppr}        & 81.6$\pm$0.3 & 68.5$\pm$0.2 & 76.7$\pm$0.3  & 91.3$\pm$0.1 & 85.4$\pm$0.3 & 90.9$\pm$0.3 & ---  \\
\textit{MoNet}           & 81.3$\pm$1.3 & 71.2$\pm$2.0 & 78.6$\pm$2.3  & 90.8$\pm$0.6 & 83.5$\pm$2.2 & 91.2$\pm$2.3 & --- \\
\textit{GS-mean}      & 79.2$\pm$7.7 & 71.6$\pm$1.9 & 77.4$\pm$2.2 & 91.3$\pm$2.8 & 82.4$\pm$1.8 & 91.4$\pm$1.3 & 71.49$\pm$0.27 \\
\textit{GS-maxpool}      & 76.6$\pm$1.9 & 67.5$\pm$2.3 & 76.1$\pm$2.3 & 85.0$\pm$1.1 & --- & 90.4$\pm$1.3 &  --- \\
\textit{CGNN}     & 81.4$\pm$1.6 & 66.9$\pm$1.8 & 66.6$\pm$4.4 & \textcolor{blue}{\bf{92.3}}$\pm$0.2 & 80.3$\pm$2.0 & 91.4$\pm$1.5 & 58.70$\pm$2.50 \\
\textit{GDE} & 78.7$\pm$2.2 & 71.8$\pm$1.1 & 73.9$\pm$3.7 & 91.6$\pm$0.1 & 82.9$\pm$0.6 & 92.4$\pm$2.0 & 56.66$\pm$10.9 \\
\hline
\textbf{\em BLEND}        & \textcolor{red}{\bf{84.8}}$\pm$0.9 & \textcolor{red}{\bf{75.9}}$\pm$1.3 & \textcolor{blue}{\bf{79.5}}$\pm$1.4 & \textcolor{red}{\bf{92.9}}$\pm$0.2 & \textcolor{red}{\bf{86.9}}$\pm$0.6 & 
\textcolor{blue}{{\bf 92.9}}$\pm$0.6 &  \textcolor{blue}{\bf{72.56}}$\pm$0.1 \\

\textbf{\em BLEND-kNN}         & \textcolor{blue}{\bf{82.5}}$\pm$0.9 & 
\textcolor{blue}{\bf{73.4}}$\pm$0.5 & 
\textcolor{red}{\bf{80.9}}$\pm$0.7 & \textcolor{blue}{{\bf 92.3}}$\pm$0.1 & \textcolor{blue}{\bf{86.7}}$\pm$0.6 & \textcolor{red}{\bf{93.5}}$\pm$0.3 & ---$^\dagger$ 
\end{tabular}
}
\caption{Performance (test accuracy$\pm$std) of different GNN models using random splits. 
$^*$\footnotesize{OGB GAT reference has 1.5M parameters vs ours 70K}.$\dagger$ \footnotesize{BLEND-kNN pre-processes the graph using the DIGL methodology~\citep{Klicpera2019} (Section~\ref{sec:gnns}), which constructs an n-dimensional representation of each node (an n-by-n matrix), then sparsifies into a kNN graph. The ogb-arxiv dataset has >150K nodes and goes OOM. This is not a limitation of BLEND, but that of DIGL. Other forms of initial rewiring are possible, but we chose to compare with DIGL (arguably the most popular graph rewiring) and so this result is missing}}
\label{tab:random_splits}
\vspace{-5mm}
\end{table*}



\subsection{Positional encoding}
\label{sec:positional_encoding}
\begin{wrapfigure}{r}{0.65\textwidth}
\vspace{-17mm}
\centering
  \includegraphics[width=\linewidth]{w_wo_pos_encs_abl.pdf}
  \vspace{-7mm}
\caption{An ablation study showing BLEND with and without positional encodings as well as GAT, the most similar conventional GNN with positional encodings added}
\vspace{-7mm}
\label{fig:pos_enc_ablation}
\end{wrapfigure}

In the second experiment, we investigated the impact of the positional encodings and vary the dimensionality and underlying geometry, in order to showcase the flexibility of our framework.  

Figure~\ref{fig:pos_enc_ablation} shows that for all datasets BLEND is superior to a Euclidean model where positional encodings are not used, which corresponds to  $Z=X$ (BLEND w/o positional in Figure~\ref{fig:pos_enc_ablation}) and a version of GAT where attention is over a concatenation of the same positional encodings used in BLEND and the features. The only exception is CoathorCS, where the performance without positional encodings is unchanged.

We experimented with three forms of positional encoding: {\em DIGL} PPR embeddings of dimension $n$ \citep{Klicpera2019}, {\em DeepWalk} embeddings~\citep{perozzi2014deepwalk} of dimensions $16$--$256$, and {\em hyperbolic} embeddings in the Poincare ball~\citep{chamberlain2017neural, nickel2018learning} of dimension $2$--$16$. 
Positional encodings are calculated as a preprocessing step and input as the $\mathbf{U}$ coordinates to BLEND. We calculated DeepWalk node embeddings using PyTorch Geometric's Node2Vec implementation with parameters $p=1$, $q=1$, context size 20, and 16 walks per node. We trained with walk length ranging between 40 and 120 and took the embeddings with the highest accuracy for the link prediction task. 
Shallow hyperbolic embeddings were generated using the HGCN~\citep{chami2019hyperbolic} implementation with the default parameters provided by the authors. 

Figure~\ref{fig:pos_encoding} (left) compares the performance of the best model using DIGL $n$-dimensional positional encodings against the best performing $d$-dimensional hyperbolic positional encodings with $d$ tuned over the range 2-16. The average performance with DIGL positional encodings is 85.48, compared to 85.28 for hyperbolic encodings. Only in one case the DIGL encodings
outperform the best hyperbolic encodings.
Figure~\ref{fig:pos_encoding} (right) show the change in performance with the dimension $d'$ of the positional encoding using a fixed hyperparameter configuration. As expected, we observe monotonic increase in performance as function of $d'$. Importantly most of the performance is captured by either 16 dimensional hyperbolic or positional encodings, with a small additional margin gained for going up to $d'=n$, which is impractical for larger graphs. 


\begin{figure}[t]
\centering
\begin{minipage}{.48\textwidth}
  \includegraphics[width=\linewidth]{GDC_HYP.pdf}
  \label{fig:test1}
\end{minipage}
\begin{minipage}{.5\textwidth}
  \includegraphics[width=\linewidth]{pos_emb_dim.pdf}
  \label{fig:test2}
\end{minipage}\vspace{-5mm}
\caption{Left: performance comparison between Euclidean and hyperbolic positional embeddings. Right: results of positional embeddings ablation. Hyperbolic or Euclidean embeddings with $d'=16$ allow to obtain performances comparable to euclidean embeddings with  $d'=n >> 16$.}
\label{fig:pos_encoding}
\end{figure}

\subsection{Additional ablations}

\begin{wrapfigure}{r}{0.55\textwidth}
\vspace{-18mm}
\centering
  \includegraphics[width=\linewidth]{step_size_abl.pdf}
\caption{step size against accurcy for explicit Euler compared to the adaptive dopri5.}
\label{fig:step_size_abl}
\vspace{-7mm}
\end{wrapfigure}

In addition to studying the affect of positional encodings , we performed ablation studies on the step size used in the integrator as well as different forms of attention. 

In Figure~\ref{fig:step_size_abl} we studied the affect of changing the step size of the integrator using the explicit Euler method with a fixed terminal time set to be the optimal terminal time. The left hand side of the figure shows the performance using the adaptive stepsize Dopri5 for comparison. Dopri5 gives the most consistent performance and is best if three of the six datasets.
Details of additional attention functions and their relative performance can be found in Appendix~\ref{app:attention}.  



\section{Related work}

\paragraph{Image processing, computer vision, and graphics.} 
Following the Perona-Malik scheme \citep{perona1990scale}, PDE-based approaches flourished in the 1990s with multiple versions of anisotropic \citep{weickert1998anisotropic} and non-Euclidean \citep{sochen1998general} diffusion used primarily for image denoising. 
The realisation of these ideas in the form of nonlinear filters \citep{tomasi1998bilateral,buades2005non} was adopted in the industry. 
PDE-based methods were also used for low-level vision tasks including  inpainting \citep{bertalmio2000image} and image segmentation \citep{caselles1997geodesic,chan2001active}.  
In computer graphics, fundamental solutions (`heat kernels') of diffusion equations were used as shape descriptors \citep{sun2009concise}.
The closed-form expression of such solutions using the Laplace-Beltrami operator served as inspiration for some of the early approaches for GNNs  \citep{HeBrLe2015,defferrard2016convolutional,kipf2017,LeMoBrBr2017}.


\paragraph{Neural differential equations} 
The interpretation of neural networks as discretised differential equations (`neural ODEs') \citep{Chen2018a} was an influential recent result with multiple follow-up works \citep{Dupont2019, Finlay2020,Liu2019, sdes2020}. 
In deep learning on graphs, this mindset has been applied to GNN architectures 
\citep{Avelar2019,Poli2019} and continuous message passing \citep{Xhonneux2020}. Continuous GNNs were also explored in~\citep{gu2020implicit} who, similarly to~\citep{Scarselli2008}, addressed the solutions of fixed point equations. Ordinary Differential Equations on Graph Networks (GODE)\citep{Zhuang2020a} approach the problem using the technique of invertible ResNets. Finally, \citep{Sanchez-Gonzalez2019} used graph-based ODEs to generate physics simulations. 

\paragraph{Physics-inspired learning}
Solving PDEs with deep learning has been explored by~\citep{Raissi2017a}. Neural networks appeared in~\citep{Li2020c} to accelerate PDE solvers with applications in the physical sciences. These have been applied to problems where the PDE can be described on a graph \citep{Li2020d}. \citep{belbute2020combining} consider the problem of predicting fluid flow and use a PDE inside a GNN. These approaches differ from ours in that they solve a given PDE, whereas we use the notion of discretising PDEs as a principle to understand and design GNNs. 

\paragraph{Neural ODEs on graphs}

The most similar work to this is GRAND~\citep{Chamberlain2021a} of which BLEND can be considered a non-Euclidean extension. In addition there are several other works that apply the neural ODE framework to graphs. In GDE~\citep{Poli2019}, GODE~\citep{Zhuang2020a} and CGNN~\citep{Xhonneux2020}, the goal is to adapt neural ordinary differential equations to graphs. In contrast, we consider non-Euclidean partial differential equations and regard GNNs as particular instances of their discretisation (both spatial and temporal). We can naturally use spaces with any metric, in particular, extending recent works on hyperbolic graph embeddings. None of the previous techniques explore the link to differential geometry. More specifically, GDE, GODE, and CGNN consider Neural ODEs of the canonical form $\tfrac{\partial x}{\partial t}=f(x,t,\theta)$ where $f$ is a graph neural network (GODE), the message passing component of a GNN (CGNN), or restricting $f$ to be layers of bijective functions on graphs (GDE). Furthermore, in CGNN only ODEs with closed form solutions are considered.
\citep{Sanchez-Gonzalez2019} on the other hand is quite distinct as they are not concerned with GNN design and instead use graph-based ODEs to generate physics simulations.








\section{Conclusion}


We developed a new class of graph neural networks based on the discretisation of a non-Euclidean diffusion PDE called Beltrami flow. 
We represent the graph structure and node features as a manifold in a joint continuous space, whose evolution by a parametric diffusion PDE (driven by the downstream learning task) gives rise to feature learning, positional encoding, and 
possibly also graph rewiring. 
Our experimental results show very good performance on popular benchmarks with a small fraction of parameters used by other GNN models.  
Perhaps most importantly, our framework establishes important links between GNNs, differential geometry, and PDEs -- fields with a rich history and profound results that in our opinion are still insufficiently explored in the graph ML community. 




\paragraph{Future directions} 
While we show that our framework generalises many popular GNN architectures we seek to use the graph as a numerical representation of an underlying continuous space. 
This view of the graph as an approximation of a continuous latent structure is a common paradigm of manifold learning \citep{tenenbaum2000global,belkin2003laplacian,coifman2006diffusion} and network geometry \citep{boguna2021network}. 
If adopted in graph ML, this mindset offers a rigorous mathematical framework for formalising and generalising some of the recent trends in the field, including the departure from the input graph as the basis for message passing \citep{Alon2020}, latent graph inference \citep{kipf2018neural,wang2019dynamic,franceschi2019learning,kazi2020differentiable} higher-order \citep{bodnar2021weisfeiler,maron2019provably} and directional \citep{beaini2020directional} message passing  (which can be expressed as anisotropic diffusion arising from additional structure of the underlying continuous space and different discretisation of the PDEs e.g. based on finite elements), and exploiting efficient numerical PDE solvers \citep{Chen2018a}. We leave these exciting directions for future research. 


\paragraph{Societal impact} GNNs have recently become increasingly utilized in industrial applications e.g. in recommender systems and social networks, and hence could potentially lead to a negative societal impact if used improperly. We would like to emphasize that our paper does not study such potential negative applications and the mathematical framework we develop could help to interpret and understand existing GNN models. We believe that better understanding of ML models is key to managing their potential societal implications and preventing their negative impact.




\paragraph{Limitations} 
The assumption that the graph can be modelled as a discretisation of some continuous space makes our framework applicable only to cases where the edge and node features are continuous in nature.  
Applications e.g. to knowledge graphs with categorical attributes could only be addressed by first embedding such attributes in a continuous space. Finally, the structural result presented in Theorem \ref{thm:gradientflow} that link our model to the Polyakov action have not been implemented. This will be addressed in future works.

\section{Acknowledgements}

We thank Nils Hammerla and Gabriele Corso for feedback on early version of this manuscript. MB and JR are supported in part by ERC Consolidator grant no 724228 (LEMAN).

