We introduce the FRactional-Order graph Neural Dynamical network (FROND), a new continuous graph neural network (GNN) framework. Unlike traditional continuous GNNs that rely on integer-order differential equations, FROND employs the Caputo fractional derivative to leverage the non-local properties of fractional calculus. This approach enables the capture of long-term dependencies in feature updates, moving beyond the Markovian update mechanisms in conventional integer-order models and offering enhanced capabilities in graph representation learning.
We offer an interpretation of the node feature updating process in FROND from a non-Markovian random walk perspective when the feature updating is particularly governed by a diffusion process.
We demonstrate analytically that oversmoothing can be mitigated in this setting.
Experimentally, we validate the FROND framework by comparing the fractional adaptations of various established integer-order continuous GNNs, demonstrating their consistently improved performance and underscoring the framework's potential as an effective extension to enhance traditional continuous GNNs.
The code is available at <https://github.com/zknus/ICLR2024-FROND>.

§ INTRODUCTION

Graph Neural Networks (GNNs) have excelled in diverse domains, e.g., chemistry <cit.>, finance <cit.>, and social media <cit.>.  The message passing scheme <cit.>, where features are aggregated along edges and iteratively propagated through layers, is crucial for the success of GNNs.
Over the past few years, numerous types of GNNs have been proposed, including Graph Convolutional Networks (GCN) <cit.>, Graph Attention Networks (GAT) <cit.>, and GraphSAGE <cit.>. Recent works, such as <cit.>, have incorporated various continuous dynamical processes to propagate information over graph nodes, giving rise to a class of continuous GNNs based on integer-order differential equations.
These continuous models have demonstrated notable performance, for instance, in enhancing robustness and addressing heterophilic graphs <cit.>.

Within these integer-order continuous GNNs, the differential operator $\ud^\beta/\ud t^\beta$ has been constrained to integer values of $\beta$, primarily 1 or 2. However, over recent decades, the wider scientific community has explored fractional-order differential operators, where $\beta$ can be any real number. These expansions have proven pivotal in various applications characterized by non-local and memory-dependent behaviors, with examples including viscoelastic materials <cit.>, anomalous transport mechanisms <cit.>, and fractal media <cit.>. Unlike conventional integer-order derivatives that measure the function's instantaneous rate of change and focus on the local vicinity, fractional-order derivatives <cit.> consider the entire historical trajectory of the function.

We introduce the FRactional-Order graph Neural Dynamical network (FROND) framework, a new approach that broadens the capabilities of traditional integer-order continuous GNNs by incorporating fractional calculus.
It naturally generalizes the integer-order derivative $\ud^\beta/\ud t^\beta$ in these GNNs to accommodate any positive real number $\beta$.
This modification gives FROND the ability to incorporate memory-dependent dynamics for information propagation and feature updating, enabling refined graph representations and improved performance potentially.
Importantly, this technique assures at least equivalent performance to integer-order models, as setting $\beta$ to integer values reverts the models to their traditional integer-order forms.

Several works like <cit.> have combined fractional graph shift operators with integer-order ordinary differential equations (ODEs). These studies are distinct from our research, wherein we focus on incorporating time-fractional derivatives for updating graph node features, modeled as a memory-inclusive dynamical process. Other works like <cit.> have used fractional calculus in gradient propagation for the training process, which is different from leveraging fractional differential equations (FDEs) in modeling the node feature updating.  We provide a detailed discussion of the differences between FROND and these works in <ref>.

Many real-world graph datasets, such as the World Wide Web, the Internet, and various biological and social networks, are known to exhibit scale-free hierarchical structures. These structures suggest a pervasive self-similarity across different scales, hinting at an underlying fractal behavior <cit.>. It has been well-established that dynamical processes with self-similarity on such fractal media are more accurately described using FDEs. For instance, the dispersion of heat or mass over these structures is best modeled using fractional diffusion equations <cit.>. Further investigations have revealed a direct connection between the fractal dimension of these structures and the order $\beta$ in fractional derivatives $\ud^\beta/\ud t^\beta$ <cit.>. This revelation births a compelling insight: the optimal $\beta$ in our models, which may differ from integers, can pave the way for enhanced node classification and potentially unearth insights into the inherent “fractality” of graph datasets.

Main contributions.
Our objective in this paper is to formulate a generalized fractional-order continuous GNN framework. Our key contributions are summarized as follows:

* We propose a novel, generalized continuous GNN framework that incorporates non-local fractional derivatives $\ud^\beta/\ud t^\beta$. This framework generalizes the prior class of integer-order continuous GNNs, subsuming them as special instances with $\beta$ setting as integers.
This approach also lays the groundwork for a diverse new class of GNNs that can accommodate a broad array of learnable memory-dependent feature-updating processes.

* We provide an interpretation from the perspective of a non-Markovian graph random walk when the feature-updating dynamics are inspired by the fractional heat diffusion process. Contrasting with the Markovian random walk implicit in traditional integer-order graph neural diffusion models whose convergence to the stationary equilibrium is exponentially swift, we establish that in FROND, convergence follows a slow algebraic rate. This characteristic enhances FROND's ability to mitigate oversmoothing, as verified by our experimental results.

* We underscore the compatibility of FROND, emphasizing its capability to be seamlessly integrated to augment the performance of existing integer-order continuous GNNs across diverse datasets. Our exhaustive experiments, encompassing the fractional differential extension of <cit.>, substantiate this claim.
Through detailed ablation studies, we provide insights into the choice of numerical schemes and parameters.

§ PRELIMINARIES

In this section, we briefly introduce fractional calculus and integer-order continuous GNNs. For a comprehensive review of fractional calculus, readers are referred to <ref>.

§.§ Caputo Fractional Derivative

The literature offers various fractional derivative definitions, notably by Riemann, Liouville, Chapman, and Caputo <cit.>. Our study leverages the Caputo fractional derivative, due to the reasons listed in <ref>.
The traditional first-order derivative of a scalar function $f(t)$ represents the local rate of change of the function at a point, defined as: $\frac{\ud f(t)}{\ud t}=\lim _{\Delta t \rightarrow 0} \frac{f(t+\Delta t)-f(t)}{\Delta t}$.
Let $F(s)$ denote the Laplace transform of $f(t)$, assumed to exist on $\left[s_0, \infty\right)$ for some $s_0 \in \mathbb{R}$. Under certain conditions <cit.>, the Laplace transform of $\frac{\ud f(t)}{\ud t}$ is given by:

\begin{align}
\mathcal{L}\left\{\frac{\ud f(t)}{\ud t}\right\}=s F(s)-f(0) \label{eq.first_lap}
\end{align}

The Caputo fractional derivative of order $\beta \in (0,1]$ for a function $f(t)$ is defined as follows:

\begin{align} \label{eq.def_fract}
D_t^\beta f(t)=\frac{1}{\Gamma(1-\beta)} \int_0^t(t-\tau)^{-\beta} f^{\prime}(\tau) \ud\tau,
\end{align}

where $\Gamma(\cdot)$ denotes the gamma function, and $f^{\prime}(\tau)$ is the first-order derivative of $f$.
The broader definition for any $\beta>0$ is deferred to <ref>.
The Caputo fractional derivative inherently integrates the entire history of the system through the integral term, emphasizing its non-local nature. For $s>\max \left\{0, s_0\right\}$, the Laplace transform of the Caputo fractional derivative is given by <cit.>[Theorem 7.1]:

\begin{align}
\mathcal{L}\left\{D_t^\beta f(t)\right\}=s^\beta F(s)-s^{\beta-1} f(0). \label{eq.fract_lap}
\end{align}

Comparing <ref>, it is evident that the Caputo derivative serves as a generalization of the first-order derivative. The alteration in the exponent of $s$ comes from the memory-dependent property in <ref>. As $\beta \to 1$, the Laplace transform of the Caputo fractional derivative converges to that of the traditional first-order derivative. When $\beta=1$, $D_t^1f=f^{\prime}$ is uniquely determined through the inverse Laplace transform <cit.>.

In summary, from the frequency domain using the Laplace transform, we observe that the Caputo fractional derivative can be seen as a natural extension of the traditional first-order derivative. For vector-valued functions, the fractional derivative is defined component-wise for each dimension.

§.§ Integer-Order Continuous GNNs

We denote an undirected graph as  $\mathcal{G}=(\mathcal{V}, \bW)$ without self-loops, where $\mathcal{V}$ is the set of $|\mathcal{V}| = N$ nodes.  The feature matrix  $\bX = \left(\left[\bx_{1}\right]\T, \cdots, \left[\bx_{N}\right]\T\right)\T \in \mathbb{R}^{N \times d}$ consists of rows $\bx_{i} \in \mathbb{R}^d$ as node feature vectors and $i$ is the node index. The $N \times N$ matrix $\bW := \left(W_{ij}\right)$ has elements $W_{ij}$ indicating the edge weight between the $i$-th and $j$-th node with $W_{ij} = W_{ji}$.
The following integer-order continuous GNNs leverage ODEs to facilitate information propagation amongst graph nodes, where features evolve as $\bX(t)$, starting from the initial condition $\bX(0) = \bX$.

GRAND: Inspired by the heat diffusion equation, GRAND <cit.> utilizes the following nonlinear autonomous dynamical system:

\begin{align}
\frac{\ud \bX(t)}{\ud t} = (\bA(\bX(t)) - \bI)\bX(t). \label{eq.grand_main}
\end{align}

where $\bA(\bX(t))\in\Real^{N\times N}$ is a learnable, time-variant attention matrix, calculated using the features $\bX(t)$, and $\bI$ denotes the identity matrix. The feature update outlined in <ref> is referred to as the GRAND-nl version (due to the nonlinearity in $\bA(\bX(t))$).
We define $d_i=\sum_{j=1}^n W_{i j}$ and let $\bD$ be a diagonal matrix with $D_{ii} = d_i$. The random walk Laplacian is then represented as $\bL = \bI - \bW\bD^{-1}$.
In a simplified context, we employ the following linear dynamical system:

\begin{align}
\frac{\ud \bX(t)}{\ud  t} = ( \bW \bD^{-1} - \bI)\bX(t) = -\bL \bX(t). \label{eq.gra_dif_L}
\end{align}

The feature update process in <ref> is the GRAND-l version. For implementations of <ref>, one may direct set $\mathbf{W D}^{-1}=\mathbf{A}(\mathbf{X}(0))$ as a column-stochastic attention matrix, rather than using a plain weight. Notably, in this time-invariant setting, the attention weight matrix, reliant on the initial node features, stays unchanged throughout the feature evolution period.

GRAND++ <cit.> adds a source term to GRAND, enhancing learning in scenarios with limited labeled nodes.
GraphCON <cit.> employs a second-order ODE, which is equivalent to two first-order ODEs, drawing inspiration from oscillator systems.
CDE <cit.> incorporates convection-diffusion equations into GNNs to address heterophilic graph challenges.
GREAD <cit.> introduces a reaction term in the GRAND model, improving its application to heterophilic graphs and formulating a diffusion-reaction equation within GNNs.
The detailed formulation for each model is presented in <ref> due to space constraints.

§ FRACTIONAL-ORDER GRAPH NEURAL DYNAMICAL NETWORK

In this section, we introduce the FROND framework, a novel approach that augments traditional integer-order continuous GNNs by incorporating fractional calculus. We elucidate the fractional counterparts of several well-established integer-order continuous GNNs, including GRAND, GRAND++, GraphCON, CDE, and GREAD, as referenced in <ref>. We provide a detailed study of the fractional extension of GRAND, and present insights into the inherent memory mechanisms in our framework through a random walk interpretation. Our theoretical findings suggest a potential mitigation of oversmoothing due to the model's slow algebraic convergence to stationarity. Subsequently, we outline the numerical FDE solvers required to implement FROND.

§.§ Framework

Consider a graph $\mathcal{G}=(\mathcal{V}, \bW)$
as defined in <ref>. Analogous to the implementation in traditional integer-order continuous GNNs, a preliminary learnable encoder function $\varphi:\calV\to\Real^d$ that maps each node to a feature vector can be applied. Stacking all the feature vectors together, we obtain $\bX\in\Real^{N\times d}$. Employing the Caputo fractional derivative outlined in <ref>, the information propagation and feature updating dynamics in FROND are characterized by the following FDE:

\begin{align}
D_t^\beta \bX(t)= \calF(\bW,\bX(t)), \quad \beta>0, \label{eq.main_object}
\end{align}

where $\beta$ denotes the fractional order of the derivative, and $\calF$ is a dynamic operator on the graph like the models presented in <ref>. The initial condition for <ref> is set as $\bX^{[\lceil\beta\rceil-1]}(0) =\ldots = \bX(0)= \bX$ consisting of the preliminary node features, with $\bX^{[i]}(t)$ denoting the $i$-th order derivative and $\lceil \beta \rceil$ is the smallest integer not less than $\beta$, akin to the initial conditions seen in integer-order ODEs.[See <ref>. We mainly consider $\beta \in (0,1]$ and the initial condition is $\bX(0)= \bX$.]
Similar to integer-order continuous GNNs, we set an integration time parameter $T$ to get $\bX(T)$. The final node embeddings for downstream tasks are then decoded using a learnable decoder $\psi(\bX(T))$.

When $\beta=1$, <ref> reverts to the class of integer-order continuous GNNs, with the infinitesimal variation of features dependent only on their present state.
Conversely, when $\beta<1$, the Caputo fractional derivative <ref> dictates that the updating process for features encompasses their entire history, not just the present state.
This paradigm facilitates memory-dependent dynamics in the framework.

For further insights into memory dependence, readers are directed to <ref>, which discusses time discretization techniques for numerically solving the system.
It illustrates how, akin to integer-order neural ODE models, time consistently acts as an analog to the layer index and how the nonlocal properties of fractional derivatives facilitate nontrivial dense or skip connections between layers.
In <ref>, when the dynamic operator $\calF$ is designated as the diffusion process in <ref>, we offer a memory-dependent non-Markovian random walk interpretation of the fractional graph neural diffusion process. Here, as $\beta\rightarrow 1$, the non-Markovian random walk increasingly detaches from the path history, becoming a Markovian walk at $\beta=1$, which is related to the normal diffusion process <cit.>. The parameter $\beta$ provides flexibility to adjust the extent of memorized dynamics embedded in the framework.
From a geometric perspective, as discussed in <ref>, the information propagation dynamics in fractal graph datasets might be more suitably described using FDEs. Choosing a non-integer $\beta$ could reveal the degree of fractality in graph datasets.

§.§.§ FROND Model Examples

When $\calF$ in <ref> is specified to the dynamics depicted in various integer-order continuous GNNs (cf. <ref>), we formulate FROND GNN variants such as F-GRAND, F-GRAND++, F-GREAD, F-CDE, and F-GraphCON, serving as fractional differential extensions of the original GNNs.

F-GRAND: Mirroring the GRAND model, the fractional-GRAND (F-GRAND) has two versions. The F-GRAND-nl version employs a time-variant FDE as follows:

\begin{align}
D_t^\beta \bX(t) = (\bA(\bX(t)) - \bI)\bX(t), \quad 0<\beta\le 1. \label{eq.frac_grand-nl}
\end{align}

It is computed using $\bX(t)$ and the attention mechanism derived from the Transformer model <cit.>. The entries of $\bA(\bX(t))=(a(\bx_i,\bx_j))$ are given by:

\begin{align}
a(\bx_i,\bx_j) = \mathrm{softmax} \left(\left\{\frac{(\mathbf{W}_K \bx_{i}\T )\T \mathbf{W}_Q \bx_{j}\T}{\bar{d}_k}\right\}\right).
\end{align}

In this formulation, $\mathbf{W}_K$ and $\mathbf{W}_Q$ are the learned matrices, and $\bar{d}_k$ signifies a hyperparameter related to the dimensionality of $\bW_K$. In parallel, the F-GRAND-l version stands as the fractional differential extension of <ref>:

\begin{align}
D_t^\beta \bX(t) = - \bL \bX(t), \quad 0<\beta\le 1. \label{eq.frac_gra_dif_L}
\end{align}

Recall that the initial condition for F-GRAND-nl and F-GRAND-l is $\bX(0)=\bX$ due to $\beta\in(0,1]$.

F-GRAND++, F-GREAD, F-CDE, and F-GraphCON: Due to space constraints, we direct the reader to <ref> for detailed formulations. Succinctly, they represent the fractional differential extensions of GRAND++, GraphCON, CDE, and GREAD. To highlight FROND's compatibility and its potential to enhance the performance of existing integer-order continuous GNNs across a variety of datasets, exhaustive experiments are provided in <ref>.

§.§ Random Walk Perspective of F-GRAND-l

The established Markov interpretation of GRAND-l <ref>, as outlined in <cit.>, aligns with F-GRAND-l <ref> when $\beta=1$. We herein broaden this interpretation to encompass a non-Markovian random walk that considers the walker's complete path history when $\beta$ is a non-integer, thereby elucidating the memory effects inherent in FROND.
In contrast to the Markovian walk, whose distribution converges exponentially to equilibrium, our strategy assures algebraic convergence, revealing F-GRAND-l's efficacy in mitigating oversmoothing as evidenced in <ref>.

To begin, we discretize the time domain into time instants as $t_n=n \sigma, \sigma>0, n=0,1,2, \ldots,$ where $\sigma$ is assumed to be small enough to ensure the validity of the approximation.
Let ${\bR(t_n)}$ be a random walk on the graph nodes $\{\bx_j\}_{j=1}^{N}$ that is not necessarily a Markov process and $\bR(t_{n+1})$ may depend on the path history $(\bR(t_0),\bR(t_1),\ldots,\bR(t_n))$ of the random walker. For convenience, we introduce the coefficients $c_k$ for $k\geq1$ and $b_n$ for $n\geq0$ from <cit.>, which are used later to define the random walk transition probability:

\begin{align}
c_k(\beta)=(-1)^{k+1}\binom{\beta}{k}=\left|\binom{\beta}{k}\right|, \quad b_n(\beta)=\sum_{k=0}^n(-1)^k\binom{\beta}{k}, \label{eq.ckbm}
\end{align}

where the generalized binomial coefficient  $\binom{\beta}{k}=\frac{\Gamma(\beta+1)}{\Gamma(k+1) \Gamma(\beta-k+1)}$ and the gamma function $\Gamma(\cdot)$ are employed in the definition of the coefficients. The sequences $c_k$ and $b_n$ consist of positive numbers, not greater than 1, decreasing strictly monotonically to zero (see supplementary material for details) and satisfy $\sum_{k=1}^{n}c_k+b_n=1$.
Using these coefficients, we define the transition probabilities of the random walk starting from $\bx_{j_0}$ as

\begin{align}
&\P(\bR(t_{n+1})=\bx_{j_{n+1}} \given \bR(t_0)=\bx_{j_0},\bR(t_1)=\bx_{j_{1}},\ldots,\bR(t_n)=\bx_{j_{n}} ) \nn
& =
\begin{cases}
c_1-\sigma^\beta & \text { if staying at current location with } j_{n+1}= j_{n}, \\
\sigma^\beta\frac{ W_{ j_{n} j_{n+1}}}{d_{j_n}}& \text { if jumping to neighboring nodes with } j_{n+1}\ne j_{n},\\
c_{n+1-k} & \text { if revisiting historical positions with } j_{n+1} = j_{k}, 1\le k\le n-1, \\
b_n & \text { if revisiting historical positions with } j_{n+1} =j_0.
\end{cases} \label{eq.random_pro}
\end{align}

This formulation integrates memory effects, considering the walker's time, position, and path history. The transition mechanism of the memory-inclusive random walk between $t_n$ and $t_{n+1}$ is elucidated as follows:
Suppose the walker is at node $j_{n}$ at time $t_n$, having a full path history $(j_{0},j_{1},\ldots,j_{n})$. We generate a  random number $\rho\in[0,1)$ uniformly, and divide the interval $[0, 1)$ into adjacent sub-intervals with lengths $c_1, c_2, \ldots, c_n, b_n$. We further subdivide the first interval (with length $c_1$) into sub-intervals of lengths $c_1 - \sigma^\beta$ and $\sigma^\beta$.

* If $\rho$ is in the first interval with length $c_1$, the walker either moves to a neighbor $j_{n+1}=k$ with probability $\sigma^\beta\frac{ W_{j_{n} k}}{d_{j_n}}$ or remains at the current position with probability $c_1-\sigma^\beta$.

* For $\rho$ in subsequent intervals, the walker jumps to a previously visited node in the history $(j_{0},j_{1},\ldots,j_{n-1})$, specifically, to $j_{n+1-k}$ if in $c_k$, or to $j_0$ if in $b_n$.

When $\beta<1$, the random walk can, with positive probability, revisit its history, restricting extensive drift.
We denote $\P(\bR(t_{n}))$ as the probability column vector, with its $j$-th element given as $\P(\bR(t_{n}) = \bx_{j})$.
Additionally, we specify $_i\P(\bR(t_{n}))$ to indicate the situation where the random walker initiates from the $i$-th node, i.e., $\bR(0)=\bx_{i}$, with probability 1. In this case, the initial probability vector $_i\P(\bR(0))$ is represented as a one-hot vector with the $i$-th entry marked as 1.
Using the technique from <cit.>, we can prove the following:

Consider the random walk defined in <ref>,  with the step size $\sigma$ and number of steps $n$. Under the conditions that $n \to \infty$ and $n\sigma = t$, the limiting probability distribution $\bP(t)\coloneqq\lim_{n \to \infty}\P(\bR(t_{n}))$ satisfies <ref>. In other words,

\begin{align}
D_t^\beta \bP(t)=-\mathbf{L} \bP(t) \label{eq.zdfva}
\end{align}

Considering that initial conditions and dimensions affect the solutions of FDEs, $\bP(t)$ and $\mathbf{X}(t)$ are not equivalent. However, due to the linearity of FDEs, the following conclusion is straightforward:

Under the conditions that $n \to \infty$ and $n\sigma = t$,  we have $\lim_{n \to \infty}\sum_i {}_i\mathbb{P}(\bR(t_{n}))\bx_i = \bX(t)$, i.e., $\sum_i {}_i\bP(t)\bx_i = \bX(t)$ with ${}_i \bP(t)\coloneqq\lim_{n \to \infty}{}_i\P(\bR(t_{n}))$ , where $\bX(t)$ is the solution to <ref> with the initial condition $\bX(0)= \bX$.

<ref> relate F-GAND-l <ref> to the non-Markovian random walk in <ref>, illustrating memory dependence in FROND.  As $\beta\rightarrow 1$, this process reverts to the Markovian random walk found in GRAND-l <cit.> in <ref>. This underscores the FROND framework's capability to apprehend more complex dynamics than integer-order continuous GNNs.

\begin{align}
&\P(\bR(t_{n+1})=\bx_{j_{n+1}} \given \bR(t_0)=\bx_{j_0},\bR(t_1)=\bx_{j_{1}},\ldots,\bR(t_n)=\bx_{j_{n}}) \vspace{-0.2cm}\label{eq.conv_random_pro}\\
&= \P (\bR(t_{n+1})=\bx_{j_{n+1}} \given \bR(t_n)=\bx_{j_{n}} ) =
\begin{cases}
1 -\sigma & \text{\scriptsize if staying at current location with } j_{n+1}= j_{n} \\
\sigma\frac{ W_{ j_{n} j_{n+1}}}{d_{j_n}} & \text{\scriptsize if jumping to neighbors with } j_{n+1}\ne j_{n}
\end{cases}\nonumber
\end{align}

since we have that all these coefficients vanishing except $c_1=1$, i.e.,

\begin{align}
c_1 =1,\quad \lim_{\beta\rightarrow 1} c_k(\beta)=0, \quad k\ge 2,\quad \lim_{\beta\rightarrow 1} b_n(\beta)=0, \quad n\ge 1. \label{eq.ckbm_limit}
\end{align}

§.§.§ Oversmoothing Mitigation of F-GRAND-l Compared to GRAND-l

The seminal research <cit.>[Corollary 3. and Remark 1] has highlighted that, when considering a GNN as a layered dynamical system, oversmoothing is a broad expression of the exponential convergence to stationary states that only retain information about graph connected components and node degrees.
Under certain conditions, the stationary distribution for the Markovian random walk <ref> is given by $\bpi=( \frac{d_1}{\sum_{j=1}^{N}d_j}, \ldots, \frac{d_N}{\sum_{j=1}^{N}d_j})$ <cit.>, with an exponentially rapid convergence rate $\left\| \P(\bR(t_{n}))- \bpi\T\right\|_{2} \sim O(e^{-r'n})$ [We use the asymptotic order notations from <cit.> in this paper.], where $r'>0$ relates to the eigenvalues of the matrix $\bL$ <cit.>, and $\|\cdot\|_2$  denotes the $\ell^2$ norm.
This behavior extends to the continuous limit, akin to a first-order linear ODE solution, exhibiting exponential convergence with some $r>0$:

\begin{align}
\left\| \bP(t)- \bpi\T\right\|_{2} \sim O(e^{-rt}).
\end{align}

In contrast, we next prove that the non-Markovian random walk <ref> converges to the stationary distribution at a slow algebraic rate, thereby helping to mitigate oversmoothing.
As $\beta\rightarrow 0$, the convergence is expected to be arbitrarily slow. In real-world scenarios where we operate within a finite horizon, this slower rate of convergence may be sufficient to alleviate oversmoothing, particularly when it is imperative for a deep model to extract distinctive features instead of achieving exponentially fast convergence to a stationary equilibrium.

Under the assumption that the graph is strongly connected and aperiodic, the stationary probability for the non-Markovian random walk <ref>, with $0<\beta<1$, is still $\bpi$, which is unique. This mirrors the stationary probability of the Markovian random walk as defined by <ref> when $\beta=1$. Notably, when $\beta<1$, the convergence of the distribution (distinct from $\boldsymbol{\pi}$) to $\boldsymbol{\pi}$  is algebraic:

\begin{align}
\left\|\bP(t) - \bpi\T \right\|_2  \sim  \Theta( t^{-\beta}). \label{eq.algebraic_conv}
\end{align}

<ref> indicate that $\bX(t)=\sum_i {}_i\bP(t)\bx_i$, as the solution to F-GRAND-l <ref>, converges to $\sum_i \bpi\T\bx_i = \bpi\T \sum_i \bx_i$ at a slow algebraic rate since $\left\|_i\bP(t) - \bpi\T \right\|_2  \sim  \Theta( t^{-\beta})$ for all $i$. Notably, $\bpi\T \sum_i \bx_i$ forms a rank 1 invariant subspace under the dynamics of <ref>, due to $\bpi$ being stationary. This underscores the difference in convergence rates, contrasting the slow algebraic rate in our case with the fast exponential rate <cit.>.

§.§ Solving FROND

every picture/.style=line width=0.75pt

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (83.13,100.77) – (82.64,154.97) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (137.09,100.77) – (136.6,154.97) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (245.03,101.47) – (244.54,155.67) ;

[color=rgb, 255:red, 208; green, 2; blue, 27   ,draw opacity=1 ][line width=2.25]    (299,101.47) – (298.51,155.67) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (83.13,100.77) .. controls (98.83,93.15) and (126.71,94.31) .. (137.09,100.77) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (83.13,100.77) .. controls (128.51,75.05) and (184.28,93.73) .. (191.06,101.01) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (137.09,100.77) .. controls (152.79,93.15) and (180.68,94.31) .. (191.06,100.77) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (191.06,101.47) .. controls (206.76,93.85) and (234.65,95.02) .. (245.03,101.47) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (245.03,101.47) .. controls (260.73,93.85) and (288.62,95.02) .. (299,101.47) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (137.09,100.77) .. controls (182.48,75.05) and (238.24,93.73) .. (245.03,101.01) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (191.06,100.77) .. controls (236.45,75.05) and (292.21,93.73) .. (299,101.01) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (83.13,100.54) .. controls (140.2,51.12) and (234.65,82.06) .. (245.03,101.47) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (137.09,100.77) .. controls (194.17,51.35) and (288.62,82.29) .. (299,101.71) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (83.13,100.77) .. controls (146.5,17.84) and (296.71,69.8) .. (299,101.47) ;

(218.45,127.82) ;
[shift=(218.45,127.82), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

(191.47,127.71) ;
[shift=(191.47,127.71), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

(164.5,127.59) ;
[shift=(164.5,127.59), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (29.16,100.77) – (28.67,154.97) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (29.16,100.77) .. controls (44.86,93.15) and (72.74,94.31) .. (83.13,100.77) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (29.16,100.19) .. controls (74.54,74.47) and (130.31,93.15) .. (137.09,100.42) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (29.16,100.77) .. controls (86.23,51.35) and (180.68,82.29) .. (191.06,101.71) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.55 ][line width=0.75]    (29.16,100.77) .. controls (92.53,17.84) and (242.74,69.8) .. (245.03,101.47) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (29.16,100.77) .. controls (122.21,-5.62) and (294.91,63.84) .. (299,101.94) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (427.13,103.44) – (426.64,157.63) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (481.09,103.44) – (480.6,157.63) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (589.03,104.14) – (588.54,158.33) ;

[color=rgb, 255:red, 208; green, 2; blue, 27   ,draw opacity=1 ][line width=2.25]    (643,104.14) – (642.51,158.33) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (427.13,103.44) .. controls (442.83,95.81) and (470.71,96.98) .. (481.09,103.44) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (481.09,103.44) .. controls (496.79,95.81) and (524.68,96.98) .. (535.06,103.44) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (535.06,104.14) .. controls (550.76,96.51) and (578.65,97.68) .. (589.03,104.14) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (589.03,104.14) .. controls (604.73,96.51) and (632.62,97.68) .. (643,104.14) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (535.06,103.44) .. controls (580.45,77.72) and (636.21,96.4) .. (643,103.67) ;

(562.45,130.49) ;
[shift=(562.45,130.49), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

(535.47,130.37) ;
[shift=(535.47,130.37), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

(508.5,130.26) ;
[shift=(508.5,130.26), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (373.16,103.44) – (372.67,157.63) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (373.16,103.44) .. controls (388.86,95.81) and (416.74,96.98) .. (427.13,103.44) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (373.16,102.86) .. controls (418.54,77.13) and (474.31,95.81) .. (481.09,103.09) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (373.16,103.44) .. controls (430.23,54.02) and (524.68,84.96) .. (535.06,104.37) ;

[color=rgb, 255:red, 74; green, 144; blue, 226   ,draw opacity=1 ] (532,111.25) .. controls (531.96,115.92) and (534.27,118.27) .. (538.94,118.31) – (577.13,118.66) .. controls (583.8,118.72) and (587.11,121.08) .. (587.06,125.75) .. controls (587.11,121.08) and (590.46,118.78) .. (597.13,118.84)(594.13,118.81) – (635.44,119.19) .. controls (640.11,119.23) and (642.46,116.92) .. (642.5,112.25) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ]   (27.67,154.97) – (314.67,156.32) ;
[shift=(316.67,156.33), rotate = 180.27] [color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (10.93,-3.29) .. controls (6.95,-1.4) and (3.31,-0.3) .. (0,0) .. controls (3.31,0.3) and (6.95,1.4) .. (10.93,3.29)   ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ]   (372.67,157.63) – (659.67,158.99) ;
[shift=(661.67,159), rotate = 180.27] [color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (10.93,-3.29) .. controls (6.95,-1.4) and (3.31,-0.3) .. (0,0) .. controls (3.31,0.3) and (6.95,1.4) .. (10.93,3.29)   ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (373.16,103.44) .. controls (466.21,-2.96) and (638.91,66.51) .. (643,104.61) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.55 ][line width=0.75]    (373.16,103.44) .. controls (436.53,20.51) and (586.74,72.46) .. (589.03,104.14) ;

(18,161.4) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( 0)$;

(70.67,160.4) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{1})$;

(126.67,161.07) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{2})$;

(234.67,161.07) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{n-1})$;

(289.67,161.73) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{n})$;

(362,163.4) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( 0)$;

(414.67,162.4) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{1})$;

(470.67,163.07) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{2})$;

(578.67,163.07) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{n-1})$;

(634.67,163.73) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{n})$;

(495,131.67) node [anchor=north west][inner sep=0.75pt]  [color=rgb, 255:red, 74; green, 144; blue, 226   ,opacity=1 ] [align=left] memory window width $\displaystyle K$ ;

Diagrams of fractional Adams–Bashforth–Moulton method with full (left) and short (right) memory.

The studies by <cit.> introduce numerical solvers specifically designed for integer-order neural ODE models.
Our research, in contrast, engages with fractional-order ODEs, entities inherently more intricate than integer-order ODEs.
To address the scenario where $\beta$ is non-integer, we introduce the fractional explicit Adams–Bashforth–Moulton solver, incorporating three variants employed in this study: the basic predictor discussed in this section, the predictor-corrector elaborated in <ref>, and the short memory principle detailed in <ref>. Additionally, we present one implicit L1 solver in <ref>. These methods exemplify how time still acts as a continuous analog to the layer index and elucidate how memory dependence manifests as nontrivial dense or skip connections between layers (see <ref>), stemming from the non-local properties of fractional derivatives.

Basic Predictor: We first employ a preliminary numerical solver  called “predictor” <cit.>  through time discretisation.  Let $h$ be a small positive discretization parameter. We have

\begin{align}
{}_{\mathrm{P}}\bX^{(k)}= \sum_{j=0}^{\lceil\beta\rceil-1} \frac{t_{k}^j}{j!} \bX^{[j]}(0)+\frac{1}{\Gamma(\beta)} \sum_{j=0}^{k-1} \mu_{j, k} \calF(\bW,\bX^{(j)}),\label{eq.pre}
\end{align}

where $\mu_{j, k}= \frac{h^\beta}{\beta}\left((k-j)^\beta-(k-1-j)^\beta\right)$, $k$ denotes the discrete time index (iteration), and $t_k=kh$ represents the discretized time steps. $\bX^{(k)}$ is the numerical approximation of  $\bX(t_k)$.
When $\beta=1$, this method simplifies to the Euler solver in <cit.> as $\mu_{j, n}\equiv h$, yielding ${}_{\mathrm{P}}\bX^{(k)} = \bX^{(k-1)}  + h\calF(\bW,\bX^{(k-1)})$.
Thus, our basic predictor can be considered as the fractional Euler method or fractional Adams–Bashforth method, which is a generalization of the Euler method used in <cit.>. However, when $\beta<1$, we need to utilize the full memory $\{\calF(\bW,\bX^{(j))}\}_{j=0}^{k-1}$.
The block diagram in <ref> shows the basic predictor and the short memory variant, highlighting the inclusion of nontrivial dense or skip connections in our framework.
A more refined visualization is conveyed in <ref>, elucidating the manner in which information propagates through layers and the graph's spatial domain.

§ EXPERIMENTS

We execute a series of experiments to illustrate that continuous GNNs formulated within the FROND framework using $D_t^\beta$ outperform their traditional counterparts based on integer-order derivatives.
Importantly, our primary aim is not to achieve state-of-the-art results, but rather to demonstrate the additional effectiveness of the FROND framework when applied to existing integer-order continuous GNNs.
In the main paper, we detail the impressive results achieved by F-GRAND, particularly emphasizing its efficacy on tree-structured data, and F-CDE, highlighting its proficiency in managing large heterophilic datasets. We also validate the slow algebraic convergence, as discussed in <ref>, by constructing deeper GNNs with non-integer $\beta<1$. To maintain consistency in the experiments presented in the main paper, the basic predictor solver is used instead of other solvers when $\beta<1$.

More Experiments In the Appendix: The <ref> section provides additional details covering various aspects such as experimental settings, described in <ref>, the computational complexity of F-GRAND in <ref>, and analysis of F-GRAND’s robustness against adversarial attacks in <ref>.
Furthermore, results related to other FROND-based continuous GNNs are extensively presented in the <ref>. In the main paper, we utilize the basic predictor, as delineated in <ref>, while the exploration of its variants is reserved for the <ref>. Additional insights into the optimal fractional-derivative order $\beta$ and fractality in graph datasets are explored in Section <ref>.

§.§ Node Classification of F-GRAND

Datasets and splitting.
We utilize datasets with varied topologies, including citation networks (Cora <cit.>, Citeseer <cit.>, Pubmed <cit.>), tree-structured datasets (Disease and Airport <cit.>), coauthor and co-purchasing graphs (CoauthorCS <cit.>, Computer and Photo <cit.>), and the ogbn-arxiv dataset <cit.>. We follow the same data splitting and pre-processing in <cit.> for Disease and Airport datasets. Consistent with experiment settings in GRAND <cit.>, we use random splits for the largest connected component of each other dataset.  We also incorporate the large-scale Ogbn-Products dataset <cit.> to demonstrate the scalability of the FROND framework, with the results displayed in <ref>.

For a comprehensive performance comparison, we select several prominent GNN models as baselines, including GCN <cit.>, and GAT <cit.>. Given the inclusion of tree-structured datasets, we also incorporate well-suited baselines: HGCN<cit.> and GIL <cit.>. To highlight the benefits of memorized dynamics in FROND, we include GRAND <cit.> as a special case of F-GRAND with $\beta=1$.
In line with <cit.>, we examine two F-GRAND variants: F-GRAND-nl <ref> and F-GRAND-l <ref>.
Graph rewiring is not explored in this study.
Where available, results from the paper <cit.> are used.

Node classification results(%) for random train-val-test splits. The best and the second-best results are highlighted in  and , respectively.

Method       Cora          Citeseer                Pubmed          CoauthorCS                 Computer                Photo              CoauthorPhy      ogbn-arxiv     Airport     Disease

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        92.8$\pm$1.0        72.2$\pm$0.3      81.6$\pm$0.6       69.8$\pm$0.5

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.0      85.7$\pm$20.3        92.5$\pm$0.90       73.7$\pm$0.1       81.6$\pm$0.4       70.4$\pm$0.5

HGCN      78.7$\pm$1.0      65.8$\pm$2.0      76.4$\pm$0.8      90.6$\pm$0.3       80.6$\pm$1.8       88.2$\pm$1.4       90.8$\pm$1.5        59.6$\pm$0.4     85.4$\pm$0.7       89.9$\pm$1.1

GIL        82.1$\pm$1.1       71.1$\pm$1.2        77.8$\pm$0.6       89.4$\pm$1.5        –       89.6$\pm$1.3     –     –        91.5$\pm$1.7     90.8$\pm$0.5

GRAND-l         83.6$\pm$1.0             73.4$\pm$0.5        78.8$\pm$1.7       92.9$\pm$0.4         83.7$\pm$1.2      92.3$\pm$0.9        93.5$\pm$0.9       71.9$\pm$0.2       80.5$\pm$9.6       74.5$\pm$3.4

GRAND-nl         82.3$\pm$1.6             70.9$\pm$1.0        77.5$\pm$1.8       92.4$\pm$0.3          82.4$\pm$2.1      92.4$\pm$0.8        91.4$\pm$1.3       71.2$\pm$0.2       90.9$\pm$1.6       81.0$\pm$6.7

F-GRAND-l        84.8$\pm$1.1             74.0$\pm$1.5        79.4$\pm$1.5       93.0$\pm$0.3        84.4$\pm$1.5      92.8$\pm$0.6        94.5$\pm$0.4      72.6$\pm$0.1          98.1$\pm$0.2      92.4$\pm$3.9

$\beta$ for F-GRAND-l    0.9      0.9     0.9    0.7    0.98     0.9     0.6     0.7     0.5    0.6

F-GRAND-nl        83.2$\pm$1.1       74.7$\pm$1.9        79.2$\pm$0.7     92.9$\pm$0.4
84.1$\pm$0.9      93.1$\pm$0.9        93.9$\pm$0.5        71.4$\pm$0.3       96.1$\pm$0.7       85.5$\pm$2.5

$\beta$ for F-GRAND-nl      0.9    0.9    0.4     0.6    0.85   0.8    0.4     0.7     0.1    0.7

oversmoothing mitigation.

tableGraph classification results.

3c|POL     3cGOS

Feature     Profile     word2vec     BERT      Profile     word2vec     BERT

GraphSage     77.60$\pm$0.68     80.36$\pm$0.68     81.22$\pm$4.81      92.10$\pm$0.08      96.58$\pm$0.22       97.07$\pm$0.23

GCN     78.28$\pm$0.52      83.89$\pm$0.53      83.44$\pm$0.38       89.53$\pm$0.49      96.28$\pm$0.08      95.96$\pm$0.75

GAT     74.03$\pm$0.53     78.69$\pm$0.78       82.71$\pm$0.19      91.18$\pm$0.23       96.57$\pm$0.34     96.61$\pm$0.45

GRAND-l     77.83$\pm$0.37     86.57$\pm$1.13      85.97$\pm$0.74       96.11$\pm$0.26     97.04$\pm$0.55      96.77$\pm$0.34

F-GRAND-l     79.49$\pm$0.43      88.69$\pm$0.37      89.29$\pm$0.93       96.40$\pm$0.19      97.40$\pm$0.03      97.53$\pm$0.14

tableNode classification accuracy of F-GRAND-l under different value of $\beta$ when time $T=8$.

$\beta$      0.1      0.3      0.5      0.7      0.9     1.0

Cora     74.80$\pm$0.42      77.0$\pm$0.98        79.60$\pm$0.91       81.56$\pm$0.30       82.68$\pm$0.64      82.37$\pm$0.59

Airport     97.09$\pm$0.87       95.80$\pm$2.03      91.66$\pm$6.34      84.36$\pm$8.04     78.73$\pm$6.33      78.88$\pm$9.67

Performance.  The results for graph node classification are summarized in <ref>, which also report the optimal $\beta$ obtained via hyperparameter tuning. Consistent with our expectations, F-GRAND surpasses GRAND across nearly all datasets, given that GRAND represents a special case of FROND with $\beta=1$. This underscores the consistent performance enhancement offered by the integration of memorized dynamics. This advantage is particularly noticeable on tree-structured datasets such as Airports and Disease, where F-GRAND markedly outperforms the baselines. For instance, F-GRAND-l outperforms both GRAND and GIL by approximately 7% on the Airport dataset. Interestingly, our experiments indicate a smaller $\beta$ (signifying greater dynamic memory) is preferable for such fractal-structured datasets, aligning with previous studies on FDEs in biological and chemical systems <cit.>. Further discussion on $\beta$ and its relation to the fractal dimension of graph datasets can be found in <ref>.

§.§ Graph Classification of F-GRAND

We employ the Fake-NewsNet datasets <cit.>, constructed from Politifact and Gossipcop fact-checking data. More details can be found in the <ref>. This dataset features three types of node features: 768-dimensional BERT features, and 300-dimensional spaCy features, both extracted using pre-trained models, and 10-dimensional profile features from Twitter accounts. The graphs in the dataset exhibit a hierarchical tree structure.
From <ref>, we observe that F-GRAND consistently outperforms GRAND with a notable edge on the POL dataset.

§.§ Oversmoothing of F-GRAND

To validate that F-GRAND mitigates the oversmoothing issue and performs well with numerous layers, we conducted an experiment using the basic predictor in the Adams Bashforth Moulton method as defined in <ref>. This allows us to generate architectures of varying depths. In this context, we utilize the fixed data splitting as described in <cit.>.
As illustrated in <ref>, optimal performance on the Cora dataset is attained with a network depth of 64 layers.
When compared to GRAND-l, F-GRAND-l maintains a consistent performance level across all datasets as the number of layers increases, with virtually no performance drop observed up to 128 layers.
This observation is consistent with our expectations, given that <ref> predicts a slow algebraic convergence. In contrast, GRAND exhibits a faster rate of performance degradation particularly on the Airport dataset. Further details on oversmoothing mitigation are in <ref>.

§.§ Ablation Study: Selection of $\beta$beta

In <ref>, we investigate the influence of $\beta$ across various graph datasets. Notably, for the Cora dataset, a larger $\beta$ is optimal, whereas, for tree-structured data, a smaller $\beta$ is preferable. This suggests that the quantity of memorized dynamics should be tailored to the dataset's topology, and a default setting of memoryless graph diffusion with $\beta=1$ may not be optimal. More comprehensive details concerning the variations in $\beta$ can be found in the appendix, specifically in <ref>.

§.§ More integer-order continuous GNNs in FROND framework

Our FROND framework can be seamlessly applied to various other integer-order continuous GNNs, as elaborated in <ref>. Specifically, here we outline the node classification results of FROND based on the CDE model in <ref>. It is evident from the results that F-CDE enhances the performance of the CDE model across almost all large heterophilic datasets. The optimal $\beta$ is determined through hyperparameter tuning. When $\beta=1$, F-CDE seamlessly reverts to CDE, and the results from the original paper are reported.
Additionally, we conduct comprehensive experiments detailed in <ref>. The results for F-GRAND++, F-GREAD, and F-GraphCON are available in <ref>, <ref>, and <ref>, respectively. Collectively, these results demonstrate that our FROND framework can significantly bolster the performance of integer-order continuous GNNs, without introducing any additional training parameters to the backbones.

Node classification accuracy(%) of large heterophilic datasets

Model     Roman-empire     Wiki-cooc     Minesweeper     Questions     Workers     Amazon-ratings

CDE       91.64$\pm$0.28      97.99$\pm$0.38        95.50$\pm$5.23       75.17$\pm$0.99      80.70$\pm$1.04      47.63$\pm$0.43

F-CDE        93.06$\pm$0.55      98.73$\pm$0.68        96.04$\pm$0.25      75.17$\pm$0.99      82.68$\pm$0.86      49.01$\pm$0.56

$\beta$ for F-CDE      0.9    0.6     0.6     1.0     0.4   0.1

§ CONCLUSION

We have introduced FROND, a novel graph learning framework that incorporates Caputo fractional derivatives to capture long-term memory in the graph feature updating dynamics. This approach has demonstrated superior performance compared to various traditional integer-order continuous GNNs. The resulting framework represents a significant advancement in graph representation learning, addressing key challenges in the field, such as oversmoothing. Our results highlight the potential of fractional calculus in enabling more effective graph learning algorithms.

§ ACKNOWLEDGMENTS AND DISCLOSURE OF FUNDING

This research is supported by the Singapore Ministry of Education Academic Research Fund Tier 2 grant MOE-T2EP20220-0002, and the National Research Foundation, Singapore and Infocomm Media Development Authority under its Future Communications Research and Development Programme. The computational work for this article was partially performed on resources of the National Supercomputing Centre, Singapore (https://www.nscc.sg).
Xuhao Li is supported by the National Natural Science Foundation of China (Grant No. 12301491) and the Anhui Provincial Natural Science Foundation (Grant No. 2208085QA02).
To improve the readability, parts of this paper have been grammatically revised using ChatGPT <cit.>.

This appendix complements the main body of our paper, providing additional details and supporting evidence for the assertions made therein. The structure of this document is as follows:

* We discuss related work in <ref>.

* We offer a concise review of fractional calculus in <ref>.

* We include more solver details and variants in <ref>.

* We present dataset statistics, experimental settings, and additional experimental results in <ref>.

* We introduce more dynamics within the FROND framework in <ref>.

* We provide proofs for all theoretical assertions made in the main paper in <ref>.

* We discuss the limitations of our work and its broader impact in the final section of this supplementary material.

§ RELATED WORK

Fractional Calculus and Its Applications

The field of fractional calculus has seen a notable surge in interest recently due to its wide-ranging applications across various domains. These include, but are not limited to, numerical analysis <cit.>, viscoelastic materials <cit.>, population growth models <cit.>, control theory <cit.>, signal processing <cit.>, financial mathematics <cit.>, and particularly in the representation of porous and fractal phenomena <cit.>.
Within these contexts, FDEs have been developed as a powerful extension to the conventional integer-ordered differential equations, offering a resilient mathematical framework for system analysis <cit.>. To illustrate, in studies related to diffusion processes, researchers have utilized fractional calculus for delineating various natural and synthetic systems, from protein diffusion in cellular membranes <cit.>, to animal migration patterns <cit.>, human mobility networks <cit.>, and even biological phenomena pertinent to respiratory tissues and neuroscience <cit.>.
Interestingly, the occurrence of subdiffusion, as modeled by FDEs, has been observed in scenarios where diffusing entities encounter intermittent obstructions due to the complex geometrical structure or interaction dynamics of the environment <cit.>.

Within the realm of deep learning, <cit.> proposes a novel approach to GNN parameter optimization using the fractional derivative. This marks a significant shift from the conventional integer-order derivative employed in optimization algorithms like SGD or Adam <cit.> with respect to the weights. The essence of their work fundamentally differs from ours, which focuses on the fractional-order evolution of node embeddings, not gradient optimization.
A detailed examination of the study by <cit.> is pivotal as it adopts fractional derivatives instead of the standard first-order derivatives during the weight updating phase of a GNN in the gradient descent.
Specifically, attention is drawn to equation (16) in <cit.>, elucidating that the fractional derivative is operational on the loss function. This stands in stark contrast to the FROND framework proposed in this work. As delineated in equation <ref> of our paper, the fractional derivative is applied to the evolving node feature, representing an implementation of a fractional-order feature updating process, thereby showcasing a clear distinction in the application of fractional derivatives.

Additionally, <cit.> incorporates insights from fractional calculus and its L1 approximation of the fractional derivative to craft a densely connected neural network. Their aim is to adeptly handle non-smooth data and counteract the vanishing gradient problem. While our research operates within a similar sphere, we have introduced fractional calculus into integer-order continuous GNNs. Our work examines the potential of fractional derivatives in node embedding evolution to address the oversmoothing issue and establishes a connection to non-Markovian dynamic processes.
Our framework paves the way for a new class of GNNs, enabling a wide spectrum of learnable feature-updating processes influenced by memory effects.

From the perspective of physics-informed machine learning, another line of research is dedicated to crafting neural networks rooted in physical laws to solve fractional PDEs. A pioneering work in this domain is the Fractional Physics Informed Neural Networks (fPINNs) <cit.>. Subsequent research, such as <cit.>, has evolved in this direction. It is worth noting that this line of research is starkly different from our problem formulation.

Integer-Order Continuous GNNs

Recent research has illuminated a fascinating intersection between differential equations and neural networks. The concept of continuous dynamical systems as a framework for deep learning has been initially explored by <cit.>.
The seminal work of <cit.> introduces neural ODEs with open-source solvers to model continuous residual layers, which has subsequently been applied to the field of GNNs. By utilizing neural ODEs, we can align the inputs and outputs of a neural network with specific physical laws, enhancing the network’s explainability <cit.>.
Additionally, separate advancements in this domain have led to improvements in neural network performance <cit.>, robustness<cit.>, and gradient stability <cit.>. In practical applications, neural ODEs are demonstrating superior performance <cit.>.
In a similar vein, <cit.> models continuous residual layers in GCN, leveraging neural ODE solvers to produce output. Further, the work of <cit.> proposes a model that considers a continuum of GNN layers, merging discrete topological structures and differential equations in a manner compatible with various static and autoregressive GNN models. The study <cit.> introduces GODE, which enables the modeling of continuous diffusion processes on graphs. It also suggests that the oversmoothing issue in GNNs may be associated with the asymptotic stability of ODEs.
Recently, GraphCON <cit.> adopts the coupled oscillator model that preserves the graph's Dirichlet energy over time and mitigates the oversmoothing problem.
In <cit.>, the authors modeled information propagation as a diffusion process of a substance from regions of higher to lower concentration. The Beltrami diffusion model is utilized in <cit.> to enhance rewiring and improve the robustness of the graph.
The study by <cit.> introduces general sheaf diffusion operators to regulate the diffusion process and maintain non-smoothness in heterophilic graphs, leading to improved node classification performance. Meanwhile, ACMP <cit.> is inspired by particle reaction-diffusion processes, taking into account repulsive and attractive force interactions between particles. Concurrently, the graph CDE model <cit.> is crafted to handle heterophilic graphs and is inspired by the convection-diffusion process. GRAND++ <cit.> leverages heat diffusion with sources to train models effectively with a limited amount of labeled training data.
Concurrently, GREAD <cit.> articulates a GNN approach, which is premised on reaction-diffusion equations, aiming to negotiate heterophilic datasets effectively. In another development, the continuous GNN as an ODE <cit.> encapsulates a graph spatial domain rewiring, leveraging the fractional order of the graph Laplacian matrix, presenting a substantial advancement in understanding graph structures.
We also recommend that interested readers refer to the recent survey <cit.> on continuous GNNs for a more thorough summarization.

Our FROND extends the above integer-order continuous GNNs by incorporating the Caputo fractional derivative. The models mentioned can be reduced from our unified mathematical framework, with variations manifesting from the choice of the dynamic operator $\calF(\bW,\bX(t))$ in <ref> and as $\beta$ equals 1 in the fractional derivative operator $D_t^\beta$.

Skip Connections in GNNs

The incorporation of skip or dense connections within network layers has been a transformative approach within deep learning literature. Initially popularized through the ResNet architecture <cit.>, this strategy introduces shortcut pathways for gradient flow during backpropagation, thereby simplifying the training of more profound networks. While this architectural design has been instrumental in improving Convolutional Neural Networks (CNNs), it has also been employed in GNNs to bolster their representational capacity and mitigate the vanishing gradient problem.
For example, the Graph U-Net <cit.> employs skip connections to enable efficient information propagation across layers. Similarly, the Jump Knowledge Network <cit.> implements a layer-aggregation mechanism that amalgamates outputs from all preceding layers, a strategy reminiscent of the dense connections found in DenseNet <cit.>.
Furthermore, the work <cit.> introduces GCNII, an extension of the standard GCN model that incorporates two simple techniques, initial residual and identity mapping, to tackle the oversmoothing problem. Expanding on the idea of depth in GNNs, <cit.> propose DeepGCNs, an innovative architecture that employs residual/dense connections along with dilated convolutions.
The work <cit.> suggests that gradient-flow message passing neural networks may be able to deal with heterophilic graphs provided that a residual connection is available.
The paper <cit.> proposes a spatial domain rewiring and focuses on long-range interactions. DRew in <cit.> does not adhere to any ODE evolutionary structure. Additionally, the skip connection in the vDRew from <cit.> specifically links an $n-k$-th layer to the $n$-th layer. This design is fundamentally different from our FDE approach.

By incorporating fractional calculus and memory effects into our framework, we not only offer a new perspective on understanding the structural design of skip connections in GNNs as a discretized fractional dynamical system, but we also establish a foundation for the development of more versatile and powerful mechanisms for graph representation learning.

§ REVIEW OF CAPUTO TIME-FRACTIONAL DERIVATIVE

We appreciate the need for a more accessible explanation of the Caputo time-fractional derivative and its derivation, as the mathematical intricacies may be challenging for some readers in the GNN community. To address this, we are providing a more comprehensive background in this section. In the main paper, we briefly touched upon fractional calculus, with a particular focus on the Caputo fractional derivative that has been employed in our work. In this appendix, we aim to provide a more detailed overview of it and explain why it is widely employed in applications. We have based our FROND framework on the assumption that the solution to the fractional differential equation exists and is unique. The appendix provides explicit conditions for this, which are automatically satisfied in most neural network designs exhibiting local Lipschitz continuity. To simplify, these conditions are akin to those for ordinary differential equations, a common assumption implicitly made in integer-order continuous GNNs such as GRAND <cit.>, GraphCON <cit.>, GRAND++ <cit.>, GREAD <cit.> and CDE <cit.>.

§.§ Caputo Fractional Derivative and Its Compatibility of Integer-order Derivative

In the main paper, our focus is predominantly on the order $\beta \in (0,1]$ for the sake of simplification. The Caputo fractional derivative of a function $f(t)$ over an interval $[0,b]$, of a general positive order $\beta \in (0,\infty)$, is defined as follows:

\begin{align}
\label{Cap_Frac}
D_t^\beta f(t)=\frac{1}{\Gamma(\lceil \beta \rceil-\beta)} \int_0^t(t-\tau)^{\lceil \beta \rceil-\beta-1} f^{[\lceil \beta \rceil]}(\tau) \mathrm{d} \tau,
\end{align}

Here, $\lceil \beta \rceil$ is the smallest integer greater than or equal to $\beta$, $\Gamma(\cdot)$ denotes the gamma function, and $f^{[\lceil \beta \rceil]}(\cdot)$ denotes the $\lceil \beta \rceil$-order derivative of $f(\cdot)$. Within this definition, it is presumed that $f^{[\lceil \beta \rceil]} \in L^1[0,b]$, i.e., $f^{[\lceil \beta \rceil]}$ is Lebesgue integrable, to ensure the well-defined nature of $D_t^\beta f(t)$ as per (<ref>) <cit.>. For a vector-valued function, the Caputo fractional derivative is defined on a component-by-component basis for each dimension, similar to the integer-order derivative. For ease of exposition, we discuss only the scalar case here, although all the following results can be generalized to vector-valued functions.
The Laplace transform for a general order $\beta \in (0,\infty)$ is presented in Theorem 7.1 <cit.> as:

\begin{align}
\label{L_Cap_Frac}
\mathcal{L} D_t^\beta f(s)=s^\beta \mathcal{L} f(s)-\sum_{k=1}^{\lceil \beta \rceil} s^{\beta -k} f^{[k-1]}(0).
\end{align}

where we assume that the Laplace transform $\mathcal{L}f$ exists on $[s_0,\infty)$ for some $s_0\in\mathbb{R}$. In contrast, for the integer-order derivative $f^{[\beta]}$ where $\beta$ is a positive integer, we also have the formulation (<ref>), with the only difference being the range of $\beta$. Therefore, as $\beta$ approaches some integer, the Laplace transform of the Caputo fractional derivative converges to the Laplace transform of the traditional integer-order derivative.  As a result, we can conclude that the Caputo fractional derivative operator generalizes the traditional integer-order derivative since their Laplace transforms coincide when $\beta$ takes an integer value. The inverse Laplace transform specifies the uniquely determined $D_t^\beta f=f^{[\beta]}$ when $\beta$ is an integer (in the sense of almost everywhere <cit.>).

Under specific reasonable conditions, we can directly present this generalization as follows. Suppose $f^{[\lceil \beta \rceil]}(t)$ (<ref>) is continuously differentiable. In this context, integration by parts can be utilized to demonstrate that

\begin{align}
% \begin{split}
D_t^\beta f(t) & =\frac{1}{\Gamma(\lceil\beta\rceil-\beta)}\Bigg(-\left[f^{[\lceil\beta\rceil]}(\tau) \frac{(t-\tau)^{\lceil\beta\rceil-\beta}}{\lceil\beta\rceil-\beta}\right]\bigg|_0^t+\int_0^t f^{[\lceil\beta\rceil+1]}(\tau) \frac{(t-\tau)^{\lceil\beta\rceil-\beta}}{\lceil\beta\rceil-\beta} \ud \tau\Bigg) \nn
& =\frac{t^{\lceil\beta\rceil-\beta} f^{[\lceil\beta\rceil]}(0)}{\Gamma(\lceil \beta \rceil-\beta+1)}+\frac{1}{\Gamma(\lceil \beta \rceil-\beta+1)} \times\int_{0}^{t}(t-\tau)^{\lceil \beta \rceil-\beta} f^{[\lceil\beta\rceil+1]}(\tau) \ud \tau.
% \end{split}
\end{align}

As $\beta\rightarrow \lceil \beta \rceil$, we have

\begin{equation}
\begin{split}
\lim_{\beta\rightarrow \lceil \beta \rceil} D_t^\beta f(t)
&= f^{[\lceil \beta \rceil]}(0)+\int_0^t f^{[\lceil\beta\rceil+1]}(\tau) \mathrm{d} \tau\\
&= f^{[\lceil \beta \rceil]} (0)+ f^{[\lceil \beta \rceil]}(t) -  f^{[\lceil \beta \rceil]}(0) \\
&= f^{[\lceil \beta \rceil]}(t).
\end{split}
\end{equation}

In parallel to the integer-order derivative, given certain conditions (<cit.>[Lemma 3.13]), the Caputo fractional derivative possesses the semigroup property as illustrated in <cit.>[Lemma 3.13]:

\begin{align}
D_t^{\varepsilon} D_t^n f= D_t^{n+\varepsilon} f.
\end{align}

Nonetheless, it is crucial to recognize that, in general, the Caputo fractional derivative does not exhibit the semigroup property, a characteristic inherent to integer-order derivatives, as detailed in <cit.>[Section 3.1].
The Caputo fractional derivative also exhibits linearity, but does not adhere to the same Leibniz and chain rules as its integer counterpart. As such properties are not utilized in our work, we refer interested readers to <cit.>[Theorem 3.17 and Remark 3.5.].

§.§ Comparison between Riemann–Liouville and Caputo Derivative

Another well-known fractional derivative is the Riemann–Liouville derivative, which, however, sees less use in practical applications (see <ref> for more insights). In this section, we offer a succinct introduction to the Riemann–Liouville derivative and compare it with Caputo's definition. The Riemann–Liouville fractional derivative is given as

\begin{align}
\widehat{D}_t^\beta f(t):=\frac{1}{\Gamma(\lceil \beta \rceil-\beta)} \frac{\mathrm{d}^{\lceil \beta \rceil}}{\mathrm{d}t^{\lceil \beta \rceil}} \int_0^t(t-\tau)^{\lceil \beta \rceil-\beta-1} f(\tau) \mathrm{d} \tau
\end{align}

Here again, we make the assumption that sufficient conditions are satisfied to ensure well-definedness (refer to <cit.>[section 2.2] for details).

We compare the Taylor expansion for the two definitions of fractional derivatives, namely the Riemann-Liouville and Caputo derivatives, with the conventional integer-order derivative. This comparison allows us to clearly highlight the distinctions among the differential equations defined under these three different approaches.

$\bullet$ Classical Integer-order Taylor Expansion: <cit.>[Theorem 2.C] Assuming that $f$ has absolutely continuous $(m-1)$-st derivative, we have that for $t\in[0, b]$,

\begin{align}
f(t)=\sum_{k=0}^{m-1} \frac{t^k}{k !} \ddfrac{^k f(0)}{t^k} + J^m\frac{\mathrm{d}^m}{\mathrm{d}t^m} f(t)  \label{eq.R7}
\end{align}

where $J^n f(t):=\frac{1}{\Gamma(n)} \int_0^t(t-\tau)^{n-1} f(\tau) \ud \tau$. Note that here, $k$ is an integer.

$\bullet$ Riemann-Liouville Fractional Taylor Expansion:  <cit.>[Theorem 2.24] Let $n>0$ and $m=\lfloor n\rfloor+1$. Assume that $f$ is such that $J^{m-n} f$ has absolutely continuous $(m-1)$-st derivative. Then,

\begin{align}
f(t)=\frac{t^{n-m}}{\Gamma(n-m+1)} J^{m-n} f(0)+\sum_{k=1}^{m-1} \frac{t^{k+n-m}}{\Gamma(k+n-m+1)} \widehat{D}_t^{k+n-m} f(0)+J^n \widehat{D}_t^n f(t) . \label{eq.R8}
\end{align}

Note that in the case $n \in \mathbb{N}$ we have $m=n+1$ and $\Gamma(n-m+1)=\Gamma(0)=\infty$, and the first term before the sum vanishes. Hence, we recover the classical result. For general $n$, the order in $\widehat{D}_t^{k+n-m}$ is not a integer.

$\bullet$ Caputo Fractional Taylor Expansion: <cit.>[Theorem 3.8.] Assume that $n \geq 0, m=\lceil n\rceil$, and $f$ has absolutely continuous $(m-1)$-st derivative. Then

\begin{align}
f(t)=\sum_{k=0}^{m-1} \frac{t^k}{k !}D_t^k f(0) + J^n  D_t^n f(t).  \label{eq.R9}
\end{align}

Note the order in $D_t^k$ is an integer. If we compare <ref>, it becomes evident that the Caputo derivative closely resembles the classical integer-order derivative in terms of Taylor expansion. This fact influences the initial conditions for the differential equations introduced in the following section.

§.§ (Caputo) Fractional Differential Equation

In this section, we first compare the initial conditions for FDEs under the Riemann-Liouville and Caputo definitions. Following this, we present the precise conditions for the existence and uniqueness of the solution to the fractional differential equation. These conditions closely align with those of ordinary differential equations, which are widely assumed by integer-order continuous GNNs <cit.>.

§.§.§ Riemann-Liouville Case
Drawing from the Riemann-Liouville fractional Taylor expansion, let us assume that $e$ is a given function with the property that there exists some function $g$ such that $g=\widehat{D}_t^\beta  e$. The solution of the Riemann-Liouville differential equation of the form

\begin{align}
\widehat{D}_t^\beta f=g   \label{eq.R10}
\end{align}

is given by

\begin{align}
f(t)=e(t)+\sum_{j=1}^{\lceil \beta \rceil} c_j t^{n-j}   \label{eq.R11}
\end{align}

where $c_j$ are arbitrary constants. In other words, to uniquely determine the solution from <ref>, we need to know the value of $\widehat{D}_t^{k+n-m} f(0)$. This is akin to a $k$ order ordinary differential equation where the initial conditions are assumed as $\frac{\mathrm{d}^k}{\mathrm{d}t^k} f(0)$, with the distinction that the order in $\widehat{D}_t^{k+n-m}$ is not an integer.

§.§.§ Caputo Case
Similarly, if $e$ is a given function with the property that $e=D_t^\beta g$ and if we intend to solve

\begin{align}
D_t^\beta f =g \label{eq.R12}
\end{align}

then we find

\begin{align}
f(t)=e(t)+\sum_{j=1}^{\lceil \beta\rceil} c_j t^{\lceil \beta \rceil-j} \label{eq.R13}
\end{align}

once more, with $c_j$ as arbitrary constants. Thus, to obtain a unique solution, it is natural to prescribe the values of integer order derivatives $f(0), D_t^1 f(0), \ldots, D_t^{\lceil \beta\rceil-1} f(0)$ in the Caputo setting, mirroring traditional ordinary differential equations.

§.§.§ Existence and Uniqueness of the (Caputo) Solution

Next, we delve into a general Caputo fractional differential equation, presented as follows:

\begin{align}
D_t^\beta y(t)=g(t, y(t)) \label{eq.R14}
\end{align}

conjoined with suitable initial conditions. As hinted in <ref>, the initial conditions take the form:

\begin{align}
D_t^k y(0)=y_0^{(k)}, \quad k=0,1, \ldots, \lceil \beta\rceil-1. \label{eq.R15}
\end{align}

The following theorem addresses the existence and uniqueness of solutions:

* Caputo existence and uniqueness theorem: <cit.>[Theorem 6.8] Let $y_0^{(0)}, \ldots, y_0^{(m-1)} \in \mathbb{R}$ and $h^*>0$. Define the set $G:=\left[0, h^*\right] \times \mathbb{R}$
and let the function $g: G \rightarrow \mathbb{R}$ be continuous and fulfill a Lipschitz condition with respect to the second variable, i.e.,

\begin{align*}
\left|g\left(x, y_1\right)-g\left(x, y_2\right)\right| \leq L\left|y_1-y_2\right|
\end{align*}

for some constant $L>0$ independent of $x, y_1$, and $y_2$. Then there uniquely exists function $y \in C\left[0, h^*\right]$ solving the initial value problem <ref>.

For a point of reference, we also provide the well-known Picard–Lindelöf uniqueness theorem for first-order ordinary differential equations.

* Picard–Lindelöf theorem <cit.>[Page 8] Let $D \subseteq \mathbb{R} \times \mathbb{R}^n$ be a closed rectangle with $\left(t_0, y_0\right) \in \operatorname{int} D$, the interior of $D$. Let $g: D \rightarrow \mathbb{R}^n$ be a function that is continuous in $t$ and Lipschitz continuous in $y$. Then, there exists some $\varepsilon>0$ such that the initial value problem

\begin{align*}
y^{\prime}(t)=g(t, y(t)), \quad y\left(t_0\right)=y_0 .
\end{align*}

has a unique solution $y(t)$ on the interval $\left[t_0, t_0+\varepsilon\right]$.

This allows us to draw parallels between the existence and uniqueness theorem of the Caputo fractional differential equation and its integer-order ordinary differential equation equivalent. We also remind readers that standard neural networks, as compositions of linear maps and pointwise non-linear activation functions with bounded derivatives (such as fully-connected and convolutional networks), satisfy global Lipschitz continuity with respect to the input. For attention neural networks, which are compositions of softmax and matrix multiplication, we observe local Lipschitz continuity.
To see this, suppose $\mathbf{v}=\operatorname{softmax}(\mathbf{u}) \in \mathbb{R}^{n \times 1}$. Then

\begin{align*}
\frac{\ud \mathbf{v}}{\partial \mathbf{u}}=\operatorname{diag}(\mathbf{v})-\mathbf{v} \mathbf{v}^{\top}=\left[\begin{array}{cccc}
v_1\left(1-v_1\right) & -v_1 v_2 & \ldots & -v_1 v_n \\
-v_2 v_1 & v_2\left(1-v_2\right) & \ldots & -v_2 v_n \\
\vdots & \vdots & \ddots & \vdots \\
-v_n v_1 & -v_n v_2 & \ldots & v_n\left(1-v_n\right)
\end{array}\right].
\end{align*}

For bounded input, we have a bounded Jacobian. All the integer-order continuous GNN works, such as recent contributions like <cit.> assume the uniqueness of the ODE solutions. This means that all the integer-order continuous GNNs can be extended by our FROND framework with fractional dynamics.

§.§ Reasons for Choosing Caputo Derivative

We now explain the reasons behind our preference for the Caputo fractional derivative:

* As previously discussed, Caputo fractional differential equations align with integer-order differential equations concerning initial conditions.

* The Caputo fractional derivative maintains a more intuitive resemblance to the integer-order derivative and satisfies the significant property of equating to zero when applied to a constant. This property is not satisfied by the Riemann-Liouville fractional derivative. Refer to <cit.>[Example 2.4. and Example 3.1.] for further clarification.

* Given its widespread application in the literature for practical use cases, numerical methods for solving Caputo fractional differential equations have been meticulously developed and exhaustively analyzed <cit.>.

§ NUMERICAL SOLVERS FOR FROND

We remind readers that numerous methods for training neural ODEs, and consequently updating the weights $\theta$ in the neural network have been proposed. These include the autodifferentiation technique in PyTorch <cit.>, the adjoint sensitivity method <cit.>, and Snode <cit.>.
In our work, we employ the most straightforward autodifferentiation technique for training FROND with fractional neural differential equations, leveraging the numerical solvers outlined in <cit.>.
While we plan to investigate more sophisticated techniques for training FROND in future work, we have open-sourced our current solver implementations in <https://github.com/zknus/torchfde>. We believe these will serve as valuable tools for the GNN community, encouraging the advancement of a unique class of GNNs that incorporate memory effects.

In traditional integer-order continuous GNNs <cit.>, the time parameter $t$ serves as a continuous analog to GNN layers, resembling the concept of neural ODEs <cit.> as continuous residual networks.
Time discretization plays a crucial role in many numerical solvers for neural ODEs. For example, the explicit Euler scheme reduces neural ODEs to residual networks with shared hidden layers <cit.>. More sophisticated discretization methods, such as adaptive step size solvers <cit.>, provide accurate solutions but require additional computational resources.

Unlike prior studies, our work involves fractional-order ODEs, which are more complex than ODEs when the derivative order $\beta$ takes non-integer values.
We present the fractional Adams–Bashforth–Moulton method with three variants utilized in this work, demonstrating how the time parameter continues to serve as a continuous analog to the layer index and how the non-local nature of fractional derivatives leads to nontrivial dense or skip connections between layers. Additionally, we also present one implicit L1 solver for solving FROND when $\beta$ is not an integer. It is worth noting that various neural ODE solvers remain applicable for FROND when $\beta$ is an integer.

We first recall the FROND framework

\begin{align*}
D_t^\beta \bX(t)= \calF(\bW,\bX(t)), \quad \beta>0,
\end{align*}

where $\beta$ denotes the fractional order of the derivative, and $\calF$ is a dynamic operator on the graph like the models presented in <ref>. The initial condition is set as $\bX^{[\lceil\beta\rceil-1]}(0) =\ldots = \bX(0)= \bX$ consisting of the preliminary node features, akin to the initial conditions seen in ODEs.

§.§ Basic predictor

Referencing <cit.>, we first employ a preliminary numerical solver  called “predictor” through time discretisation $t_j=jh$, where the discretisation parameter $h$ is a small positive value:

\begin{align}
{}_{\mathrm{P}}\bX^{(k)}= \sum_{j=0}^{\lceil\beta\rceil-1} \frac{t_{k}^j}{j!} \bX^{[j]}(0)+\frac{1}{\Gamma(\beta)} \sum_{j=0}^{k-1} \mu_{j, k} \calF(\bW,\bX^{(j)}),
\end{align}

where $\mu_{j, n}= \frac{h^\beta}{\beta}\left((n-j)^\beta-(n-1-j)^\beta\right)$, $k$ denotes the discrete time index (iteration), and $t_k=kh$ represents the discretized time steps. $\bX^{(k)}$ is the numerical approximation of  $\bX(t_k)$.
When $\beta=1$, this method simplifies to the Euler solver in <cit.> as $\mu_{j, n}\equiv h$, yielding ${}_{\mathrm{P}}\bX^{(k)} = \bX^{(k-1)}  + h\calF(\bW,\bX^{(k-1)})$.
Thus, our basic predictor can be considered as the fractional Euler method or fractional Adams–Bashforth method, which is a generalization of the Euler method used in <cit.>. However, when $\beta<1$, we need to utilize the full memory $\{\calF(\bW,\bX^{(j))}\}_{j=0}^{k-1}$.

The block diagram of this basic predictor, shown in <ref>, reveals that our framework introduces nontrivial dense or skip connections between layers.
A more refined visualization is conveyed in <ref>, elucidating the manner in which information propagates through layers and the graph's spatial domain.

§.§ Predictor-corrector

The corrector formula from <cit.>, a fractional variant of the one-step Adams-Moulton method, refines the initial approximation using the predictor ${}_{\mathrm{P}}\bX^{(k)}$ as follows:

\begin{align}
\bX^{(k)}= \sum_{j=0}^{\lceil\beta\rceil-1} \frac{t_{k}^j}{j !} \bX^{[j]}(0) +\frac{1}{\Gamma(\beta)} \sum_{j=0}^{k-1} \eta_{j, k} \calF(\bW,\bX^{(j)}) +\frac{1}{\Gamma(\beta)} \eta_{k, k} \calF(\bW,{}_{\mathrm{P}}\bX^{(k)}). \label{eq.cor}
\end{align}

Here we show the coefficients $\eta_{j, n}$ in the predictor-corrector variant <ref> from <cit.>:

\begin{align}
\eta_{j, k}(\beta)=\frac{h^\beta}{\beta(\beta+1)} \times
\begin{cases}
(k-1)^{\beta+1}-(k-1-\beta)k^\beta & \text { if } j=0, \\
(k-j+1)^{\beta+1}+(k-1-j)^{\beta+1}-2(k-j)^{\beta+1} & \text { if } 1 \le j \le k-1, \\
1 & \text { if } j=k.
\end{cases} \label{eq.solver_coeffi}
\end{align}

§.§ Short memory principle

When $T$ is large, computational time complexity becomes a challenge due to the non-local nature of fractional derivatives. To mitigate this, <cit.> suggest leveraging the short memory principle to modify the summation in <ref> to $\sum_{j=n-K}^{n-1}$. This corresponds to employing a shifting memory window with a fixed width $K$. The block diagram is depicted in <ref>.

every picture/.style=line width=0.75pt

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (83.13,100.77) – (82.64,154.97) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (137.09,100.77) – (136.6,154.97) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (245.03,101.47) – (244.54,155.67) ;

[color=rgb, 255:red, 208; green, 2; blue, 27   ,draw opacity=1 ][line width=2.25]    (299,101.47) – (298.51,155.67) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (83.13,100.77) .. controls (98.83,93.15) and (126.71,94.31) .. (137.09,100.77) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (83.13,100.77) .. controls (128.51,75.05) and (184.28,93.73) .. (191.06,101.01) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (137.09,100.77) .. controls (152.79,93.15) and (180.68,94.31) .. (191.06,100.77) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (191.06,101.47) .. controls (206.76,93.85) and (234.65,95.02) .. (245.03,101.47) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (245.03,101.47) .. controls (260.73,93.85) and (288.62,95.02) .. (299,101.47) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (137.09,100.77) .. controls (182.48,75.05) and (238.24,93.73) .. (245.03,101.01) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (191.06,100.77) .. controls (236.45,75.05) and (292.21,93.73) .. (299,101.01) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (83.13,100.54) .. controls (140.2,51.12) and (234.65,82.06) .. (245.03,101.47) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (137.09,100.77) .. controls (194.17,51.35) and (288.62,82.29) .. (299,101.71) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (83.13,100.77) .. controls (146.5,17.84) and (296.71,69.8) .. (299,101.47) ;

(218.45,127.82) ;
[shift=(218.45,127.82), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

(191.47,127.71) ;
[shift=(191.47,127.71), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

(164.5,127.59) ;
[shift=(164.5,127.59), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (29.16,100.77) – (28.67,154.97) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (29.16,100.77) .. controls (44.86,93.15) and (72.74,94.31) .. (83.13,100.77) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (29.16,100.19) .. controls (74.54,74.47) and (130.31,93.15) .. (137.09,100.42) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (29.16,100.77) .. controls (86.23,51.35) and (180.68,82.29) .. (191.06,101.71) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.55 ][line width=0.75]    (29.16,100.77) .. controls (92.53,17.84) and (242.74,69.8) .. (245.03,101.47) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (29.16,100.77) .. controls (122.21,-5.62) and (294.91,63.84) .. (299,101.94) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (427.13,103.44) – (426.64,157.63) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (481.09,103.44) – (480.6,157.63) ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (589.03,104.14) – (588.54,158.33) ;

[color=rgb, 255:red, 208; green, 2; blue, 27   ,draw opacity=1 ][line width=2.25]    (643,104.14) – (642.51,158.33) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (427.13,103.44) .. controls (442.83,95.81) and (470.71,96.98) .. (481.09,103.44) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (481.09,103.44) .. controls (496.79,95.81) and (524.68,96.98) .. (535.06,103.44) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (535.06,104.14) .. controls (550.76,96.51) and (578.65,97.68) .. (589.03,104.14) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (589.03,104.14) .. controls (604.73,96.51) and (632.62,97.68) .. (643,104.14) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (535.06,103.44) .. controls (580.45,77.72) and (636.21,96.4) .. (643,103.67) ;

(562.45,130.49) ;
[shift=(562.45,130.49), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

(535.47,130.37) ;
[shift=(535.47,130.37), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

(508.5,130.26) ;
[shift=(508.5,130.26), rotate = 0] [color=rgb, 255:red, 0; green, 0; blue, 0   ][fill=rgb, 255:red, 0; green, 0; blue, 0   ][line width=0.75]      (0, 0) circle [x radius= 1.34, y radius= 1.34]   ;

[color=rgb, 255:red, 241; green, 168; blue, 156   ,draw opacity=1 ][line width=2.25]    (373.16,103.44) – (372.67,157.63) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (373.16,103.44) .. controls (388.86,95.81) and (416.74,96.98) .. (427.13,103.44) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (373.16,102.86) .. controls (418.54,77.13) and (474.31,95.81) .. (481.09,103.09) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.56 ][line width=0.75]    (373.16,103.44) .. controls (430.23,54.02) and (524.68,84.96) .. (535.06,104.37) ;

[color=rgb, 255:red, 74; green, 144; blue, 226   ,draw opacity=1 ] (532,111.25) .. controls (531.96,115.92) and (534.27,118.27) .. (538.94,118.31) – (577.13,118.66) .. controls (583.8,118.72) and (587.11,121.08) .. (587.06,125.75) .. controls (587.11,121.08) and (590.46,118.78) .. (597.13,118.84)(594.13,118.81) – (635.44,119.19) .. controls (640.11,119.23) and (642.46,116.92) .. (642.5,112.25) ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ]   (27.67,154.97) – (314.67,156.32) ;
[shift=(316.67,156.33), rotate = 180.27] [color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (10.93,-3.29) .. controls (6.95,-1.4) and (3.31,-0.3) .. (0,0) .. controls (3.31,0.3) and (6.95,1.4) .. (10.93,3.29)   ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ]   (372.67,157.63) – (659.67,158.99) ;
[shift=(661.67,159), rotate = 180.27] [color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (10.93,-3.29) .. controls (6.95,-1.4) and (3.31,-0.3) .. (0,0) .. controls (3.31,0.3) and (6.95,1.4) .. (10.93,3.29)   ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (373.16,103.44) .. controls (466.21,-2.96) and (638.91,66.51) .. (643,104.61) ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=0.55 ][line width=0.75]    (373.16,103.44) .. controls (436.53,20.51) and (586.74,72.46) .. (589.03,104.14) ;

(18,161.4) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( 0)$;

(70.67,160.4) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{1})$;

(126.67,161.07) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{2})$;

(234.67,161.07) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{n-1})$;

(289.67,161.73) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{n})$;

(362,163.4) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( 0)$;

(414.67,162.4) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{1})$;

(470.67,163.07) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{2})$;

(578.67,163.07) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{n-1})$;

(634.67,163.73) node [anchor=north west][inner sep=0.75pt]  [font=]  $\mathbf{X}( t_{n})$;

(495,131.67) node [anchor=north west][inner sep=0.75pt]  [color=rgb, 255:red, 74; green, 144; blue, 226   ,opacity=1 ] [align=left] memory window width $\displaystyle K$ ;

Diagrams of fractional Adams–Bashforth–Moulton method with full (left) and short (right) memory.

every picture/.style=line width=0.75pt

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (190.39,67.6) – (316.6,45.3) – (315.25,98.07) – (189.03,120.37) – cycle ;

[fill=rgb, 255:red, 74; green, 144; blue, 226   ,fill opacity=1 ] (287.1,76.9) .. controls (286.37,75.4) and (287.66,73.85) .. (289.98,73.44) .. controls (292.31,73.02) and (294.79,73.89) .. (295.52,75.39) .. controls (296.25,76.88) and (294.96,78.43) .. (292.64,78.85) .. controls (290.32,79.27) and (287.84,78.39) .. (287.1,76.9) – cycle ;

(239.42,85.71) .. controls (238.69,84.22) and (239.97,82.66) .. (242.3,82.25) .. controls (244.62,81.83) and (247.1,82.7) .. (247.84,84.2) .. controls (248.57,85.69) and (247.28,87.24) .. (244.96,87.66) .. controls (242.63,88.08) and (240.15,87.2) .. (239.42,85.71) – cycle ;

(240.13,86.9) – (219.11,100.76) ;

(210.62,81.14) – (238.8,83.77) ;

[color=rgb, 255:red, 2; green, 122; blue, 250   ,draw opacity=0.4 ]   (249.8,83.83) – (262.59,81.47) – (264.6,81.09) – (287.23,76.91) ;
[shift=(247.84,84.2), rotate = 349.51] [fill=rgb, 255:red, 2; green, 122; blue, 250   ,fill opacity=0.4 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

(202.21,82.65) .. controls (201.47,81.16) and (202.76,79.61) .. (205.09,79.19) .. controls (207.41,78.77) and (209.89,79.64) .. (210.62,81.14) .. controls (211.36,82.63) and (210.07,84.18) .. (207.74,84.6) .. controls (205.42,85.02) and (202.94,84.15) .. (202.21,82.65) – cycle ;

(211.79,102.98) .. controls (211.05,101.48) and (212.34,99.93) .. (214.67,99.51) .. controls (216.99,99.1) and (219.47,99.97) .. (220.2,101.46) .. controls (220.94,102.96) and (219.65,104.51) .. (217.32,104.93) .. controls (215,105.34) and (212.52,104.47) .. (211.79,102.98) – cycle ;

[color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ]   (192,149.5) – (349.34,254.39) ;
[shift=(351,255.5), rotate = 213.69] [color=rgb, 255:red, 0; green, 0; blue, 0   ,draw opacity=1 ][line width=0.75]    (10.93,-3.29) .. controls (6.95,-1.4) and (3.31,-0.3) .. (0,0) .. controls (3.31,0.3) and (6.95,1.4) .. (10.93,3.29)   ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (233.99,97.3) – (360.2,75) – (358.85,127.77) – (232.63,150.07) – cycle ;

[fill=rgb, 255:red, 74; green, 144; blue, 226   ,fill opacity=1 ] (330.7,106.1) .. controls (329.97,104.6) and (331.26,103.05) .. (333.58,102.64) .. controls (335.91,102.22) and (338.39,103.09) .. (339.12,104.59) .. controls (339.85,106.08) and (338.56,107.63) .. (336.24,108.05) .. controls (333.92,108.47) and (331.44,107.59) .. (330.7,106.1) – cycle ;

(283.02,114.91) .. controls (282.29,113.42) and (283.57,111.86) .. (285.9,111.45) .. controls (288.22,111.03) and (290.7,111.9) .. (291.44,113.4) .. controls (292.17,114.89) and (290.88,116.44) .. (288.56,116.86) .. controls (286.23,117.28) and (283.75,116.4) .. (283.02,114.91) – cycle ;

(283.73,116.1) – (262.71,129.96) ;

(254.22,110.34) – (282.4,112.97) ;

[color=rgb, 255:red, 2; green, 122; blue, 250   ,draw opacity=0.4 ]   (293.4,113.03) – (306.19,110.67) – (308.2,110.29) – (330.83,106.11) ;
[shift=(291.44,113.4), rotate = 349.51] [fill=rgb, 255:red, 2; green, 122; blue, 250   ,fill opacity=0.4 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

(245.81,111.85) .. controls (245.07,110.36) and (246.36,108.81) .. (248.69,108.39) .. controls (251.01,107.97) and (253.49,108.84) .. (254.22,110.34) .. controls (254.96,111.83) and (253.67,113.38) .. (251.34,113.8) .. controls (249.02,114.22) and (246.54,113.35) .. (245.81,111.85) – cycle ;

(255.39,132.18) .. controls (254.65,130.68) and (255.94,129.13) .. (258.27,128.71) .. controls (260.59,128.3) and (263.07,129.17) .. (263.8,130.66) .. controls (264.54,132.16) and (263.25,133.71) .. (260.92,134.13) .. controls (258.6,134.54) and (256.12,133.67) .. (255.39,132.18) – cycle ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (279.59,126.1) – (405.8,103.8) – (404.45,156.57) – (278.23,178.87) – cycle ;

[fill=rgb, 255:red, 74; green, 144; blue, 226   ,fill opacity=1 ] (376.7,134.5) .. controls (375.97,133) and (377.26,131.45) .. (379.58,131.04) .. controls (381.91,130.62) and (384.39,131.49) .. (385.12,132.99) .. controls (385.85,134.48) and (384.56,136.03) .. (382.24,136.45) .. controls (379.92,136.87) and (377.44,135.99) .. (376.7,134.5) – cycle ;

(329.02,143.31) .. controls (328.29,141.82) and (329.57,140.26) .. (331.9,139.85) .. controls (334.22,139.43) and (336.7,140.3) .. (337.44,141.8) .. controls (338.17,143.29) and (336.88,144.84) .. (334.56,145.26) .. controls (332.23,145.68) and (329.75,144.8) .. (329.02,143.31) – cycle ;

(329.73,144.5) – (308.71,158.36) ;

(300.22,138.74) – (328.4,141.37) ;

[color=rgb, 255:red, 2; green, 122; blue, 250   ,draw opacity=0.4 ]   (339.4,141.43) – (352.19,139.07) – (354.2,138.69) – (376.83,134.51) ;
[shift=(337.44,141.8), rotate = 349.51] [fill=rgb, 255:red, 2; green, 122; blue, 250   ,fill opacity=0.4 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

(291.81,140.25) .. controls (291.07,138.76) and (292.36,137.21) .. (294.69,136.79) .. controls (297.01,136.37) and (299.49,137.24) .. (300.22,138.74) .. controls (300.96,140.23) and (299.67,141.78) .. (297.34,142.2) .. controls (295.02,142.62) and (292.54,141.75) .. (291.81,140.25) – cycle ;

(301.39,160.58) .. controls (300.65,159.08) and (301.94,157.53) .. (304.27,157.11) .. controls (306.59,156.7) and (309.07,157.57) .. (309.8,159.06) .. controls (310.54,160.56) and (309.25,162.11) .. (306.92,162.53) .. controls (304.6,162.94) and (302.12,162.07) .. (301.39,160.58) – cycle ;

[color=rgb, 255:red, 4; green, 116; blue, 248   ,draw opacity=0.4 ]   (339.12,104.59) .. controls (378.42,112.07) and (376.4,117.04) .. (382.32,132.12) ;
[shift=(383,133.8), rotate = 247.48] [fill=rgb, 255:red, 4; green, 116; blue, 248   ,fill opacity=0.4 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

[color=rgb, 255:red, 4; green, 116; blue, 248   ,draw opacity=0.4 ]   (295.52,75.39) .. controls (332.84,84.19) and (330.87,88.96) .. (336.19,101.22) ;
[shift=(337,103), rotate = 244.8] [fill=rgb, 255:red, 4; green, 116; blue, 248   ,fill opacity=0.4 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (324.39,154.6) – (450.6,132.3) – (449.25,185.07) – (323.03,207.37) – cycle ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (328.39,157.6) – (454.6,135.3) – (453.25,188.07) – (327.03,210.37) – cycle ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (333.39,160.1) – (459.6,137.8) – (458.25,190.57) – (332.03,212.87) – cycle ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (337.77,163.45) – (463.99,141.15) – (462.63,193.91) – (336.42,216.21) – cycle ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (342.27,165.95) – (468.49,143.65) – (467.13,196.41) – (340.92,218.71) – cycle ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (347.27,169.45) – (473.49,147.15) – (472.13,199.91) – (345.92,222.21) – cycle ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (353.29,175.06) – (479.5,152.76) – (478.15,205.53) – (351.93,227.83) – cycle ;

[fill=rgb, 255:red, 74; green, 144; blue, 226   ,fill opacity=1 ] (450,183.86) .. controls (449.27,182.36) and (450.56,180.81) .. (452.88,180.4) .. controls (455.21,179.98) and (457.69,180.85) .. (458.42,182.35) .. controls (459.15,183.84) and (457.86,185.39) .. (455.54,185.81) .. controls (453.22,186.23) and (450.74,185.35) .. (450,183.86) – cycle ;

(402.32,192.67) .. controls (401.59,191.17) and (402.87,189.62) .. (405.2,189.21) .. controls (407.52,188.79) and (410,189.66) .. (410.74,191.16) .. controls (411.47,192.65) and (410.18,194.2) .. (407.86,194.62) .. controls (405.53,195.04) and (403.05,194.16) .. (402.32,192.67) – cycle ;

(403.03,193.86) – (382.01,207.72) ;

(373.52,188.1) – (401.7,190.73) ;

[color=rgb, 255:red, 2; green, 122; blue, 250   ,draw opacity=0.4 ]   (412.7,190.79) – (425.49,188.43) – (427.5,188.05) – (450.13,183.86) ;
[shift=(410.74,191.16), rotate = 349.51] [fill=rgb, 255:red, 2; green, 122; blue, 250   ,fill opacity=0.4 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

(365.11,189.61) .. controls (364.37,188.12) and (365.66,186.57) .. (367.99,186.15) .. controls (370.31,185.73) and (372.79,186.6) .. (373.52,188.1) .. controls (374.26,189.59) and (372.97,191.14) .. (370.64,191.56) .. controls (368.32,191.98) and (365.84,191.11) .. (365.11,189.61) – cycle ;

(374.69,209.94) .. controls (373.95,208.44) and (375.24,206.89) .. (377.57,206.47) .. controls (379.89,206.06) and (382.37,206.93) .. (383.1,208.42) .. controls (383.84,209.92) and (382.55,211.47) .. (380.22,211.89) .. controls (377.9,212.3) and (375.42,211.43) .. (374.69,209.94) – cycle ;

[color=rgb, 255:red, 155; green, 155; blue, 155   ,draw opacity=1 ][fill=rgb, 255:red, 255; green, 255; blue, 255   ,fill opacity=1 ] (398.29,205.1) – (524.5,182.8) – (523.15,235.57) – (396.93,257.87) – cycle ;

[fill=rgb, 255:red, 74; green, 144; blue, 226   ,fill opacity=1 ] (495,213.9) .. controls (494.27,212.4) and (495.56,210.85) .. (497.88,210.44) .. controls (500.21,210.02) and (502.69,210.89) .. (503.42,212.39) .. controls (504.15,213.88) and (502.86,215.43) .. (500.54,215.85) .. controls (498.22,216.27) and (495.74,215.39) .. (495,213.9) – cycle ;

(447.32,222.71) .. controls (446.59,221.22) and (447.87,219.66) .. (450.2,219.25) .. controls (452.52,218.83) and (455,219.7) .. (455.74,221.2) .. controls (456.47,222.69) and (455.18,224.24) .. (452.86,224.66) .. controls (450.53,225.08) and (448.05,224.2) .. (447.32,222.71) – cycle ;

(448.03,223.9) – (427.01,237.76) ;

(418.52,218.14) – (446.7,220.77) ;

[color=rgb, 255:red, 2; green, 122; blue, 250   ,draw opacity=0.4 ]   (457.7,220.83) – (470.49,218.47) – (472.5,218.09) – (495.13,213.91) ;
[shift=(455.74,221.2), rotate = 349.51] [fill=rgb, 255:red, 2; green, 122; blue, 250   ,fill opacity=0.4 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

(410.11,219.65) .. controls (409.37,218.16) and (410.66,216.61) .. (412.99,216.19) .. controls (415.31,215.77) and (417.79,216.64) .. (418.52,218.14) .. controls (419.26,219.63) and (417.97,221.18) .. (415.64,221.6) .. controls (413.32,222.02) and (410.84,221.15) .. (410.11,219.65) – cycle ;

(419.69,239.98) .. controls (418.95,238.48) and (420.24,236.93) .. (422.57,236.51) .. controls (424.89,236.1) and (427.37,236.97) .. (428.1,238.46) .. controls (428.84,239.96) and (427.55,241.51) .. (425.22,241.93) .. controls (422.9,242.34) and (420.42,241.47) .. (419.69,239.98) – cycle ;

[color=rgb, 255:red, 4; green, 116; blue, 248   ,draw opacity=1 ]   (293.02,73.39) .. controls (497.47,7.08) and (513.84,135.95) .. (501.69,209.64) ;
[shift=(501.5,210.75), rotate = 279.65] [fill=rgb, 255:red, 4; green, 116; blue, 248   ,fill opacity=1 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

[color=rgb, 255:red, 4; green, 116; blue, 248   ,draw opacity=0.4 ]   (458.42,182.35) .. controls (494.12,191.28) and (495.16,197.09) .. (497.5,208.57) ;
[shift=(497.88,210.44), rotate = 257.98] [fill=rgb, 255:red, 4; green, 116; blue, 248   ,fill opacity=0.4 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

[color=rgb, 255:red, 4; green, 116; blue, 248   ,draw opacity=1 ]   (334.58,102.64) .. controls (474.62,75.43) and (504.7,118.32) .. (499.96,209.06) ;
[shift=(499.88,210.44), rotate = 273.19] [fill=rgb, 255:red, 4; green, 116; blue, 248   ,fill opacity=1 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

[color=rgb, 255:red, 4; green, 116; blue, 248   ,draw opacity=1 ]   (383.58,130.64) .. controls (492.9,114.91) and (495.95,152.49) .. (497.83,208.73) ;
[shift=(497.88,210.44), rotate = 268.11] [fill=rgb, 255:red, 4; green, 116; blue, 248   ,fill opacity=1 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

[color=rgb, 255:red, 4; green, 116; blue, 248   ,draw opacity=1 ]   (290.73,75.78) .. controls (349.79,58.6) and (386.51,101.48) .. (383.25,131.95) ;
[shift=(383,133.8), rotate = 279.29] [fill=rgb, 255:red, 4; green, 116; blue, 248   ,fill opacity=1 ][line width=0.08]  [draw opacity=0] (7.2,-1.8) – (0,0) – (7.2,1.8) – cycle    ;

(184,220) node [anchor=north west][inner sep=0.75pt]  [color=rgb, 255:red, 0; green, 0; blue, 0   ,opacity=1 ] [align=left] time discretization;

(190.77,103.94) node [anchor=north west][inner sep=0.75pt]  [font=]  $\bX(0)$;

(234.37,134.14) node [anchor=north west][inner sep=0.75pt]  [font=]  $\bX(t_1)$;

(280.37,161.94) node [anchor=north west][inner sep=0.75pt]  [font=]  $\bX(t_2)$;

(353.67,211.9) node [anchor=north west][inner sep=0.75pt]  [font=]  $\bX(t_{n-1})$;

(398.67,242.94) node [anchor=north west][inner sep=0.75pt]  [font=]  $\bX(t_{n})$;

Model discretization in FROND with the basic predictor solver. Unlike the Euler discretization in ODEs, FDEs incorporate connections to historical times, introducing memory effects. Specifically, the dark blue connections observed in FDEs are absent in ODEs. The weight of these skip connections correlates with $\mu_{j, k}(\beta)$ as detailed in <ref>.

§.§ L1 Solver

The L1 scheme is one of the most popular methods to approximate the Caputo fractional derivative in time. It utilizes a backward differencing method for effective approximation of derivatives. Referencing <cit.>, we have the L1 approximation of Caputo fractional derivative as follows:

\begin{align*}
D^{\beta}_t \bX^{(k)} \approx \mu \sum_{j=0}^{k-1} R_{k,j}^{\beta}(\bX^{(j+1)}-\bX^{(j)})
\end{align*}

where $h$ is the temporal step size,

\begin{align*}
\mu=\frac{1}{h^{\beta}\Gamma(2-\beta)},\qquad R_{k,j}^{\beta}=(k-j)^{1-\beta}-(k-j-1)^{1-\beta}, \qquad 0\leq j \leq k-1.
\end{align*}

Applying the L1 solver for our problem, we obtain

\begin{align*}
\mu \sum_{j=0}^{k-1} R_{k,j}^{\beta}(\bX^{({j+1})}-\bX^{(j)}) = (\bA(\bX^{(k)})-\bI)\bX^{(k)}.
\end{align*}

Manipulating the above equation, we have

\begin{align*}
\bX^{(k)}-\frac{1}{\mu}(\bA(\bX^{(k)})-\bI)\bX^{(k)})= \bX^{(k-1)}-\sum_{j=0}^{k-2} R_{k,j}^{\beta}(\bX^{(j+1)}-\bX^{(j)})
\end{align*}

The above formula is an implicit nonlinear scheme. To solve it without calculating the inversion of a matrix, we propose the following iteration method:

* Compute a basic approximation of $\bX(t_k)$ with the following formula:

\begin{align*}
{}_{\mathrm{P}}\bX^{(k)}-\frac{1}{\mu}(\bA(\bX^{(k-1)})-\bI)\bX^{(k-1)}= \bX^{(k-1)}-\sum_{j=0}^{k-2} R_{k,j}^{\beta}(\bX^{(j+1)}-\bX^{(j)}).
\end{align*}

* Substitute the above ${}_{\mathrm{P}}\bX^{(k)}$ into the implicit scheme to update $\bX^{(k)}$:

\begin{align}
\label{eq.l1solver}
\bX^{(k)}-\frac{1}{\mu}(\bA({}_{\mathrm{P}}\bX^{(k)})-\bI){}_{\mathrm{P}}\bX^{(k)}= \bX^{(k-1)}-\sum_{j=0}^{k-2} R_{k,j}^{\beta}(\bX^{(j+1)}-\bX^{(j)}).
\end{align}

The step (2) can be repeated multiple times to obtain an accurate approximation of $\bX(t_k)$.

§ DATASETS, SETTINGS AND MORE EXPERIMENTS FOR F-GRAND MODEL

§.§ Datasets

The dataset statistics used in <ref> are provided in <ref>. Following the experimental framework in <cit.>, we select the largest connected component from each dataset, except for the tree-like graph datasets (Airport and Disease). However, for the study of oversmoothing, we use a fixed data splitting approach over the entire datasets, as described in <cit.>.

Dataset Statistics used in <ref>

Dataset     Type     Classes     Features     Nodes     Edges

Cora     citation     7     1433     2485     5069

Citeseer     citation     6     3703     2120     3679

PubMed     citation     3     500     19717     44324

Coauthor CS     co-author     15     6805     18333     81894

Computers     co-purchase     10     767     13381     245778

Photos     co-purchase     8     745     7487     119043

CoauthorPhy     co-author      5     8415     34493     247962

OGB-Arxiv     citation     40     128     169343     1166243

Airport     tree-like     4     4     3188     3188

Disease     tree-like    2    1000    1044    1043

Dataset and graph statistics used in <ref>

Dataset     Graphs (Fake)     Total Nodes     Total Edges     Avg. Nodes per Graph

Politifact (POL)     314 (157)     41,054     40,740     131

Gossipcop (GOS)     5464 (2732)     314,262     308,798     58

§.§ Graph Classification Details

We use the Fake-NewsNet datasets from <cit.>, constructed based on fact-checking information obtained from Politifact and Gossipcop.
The dataset incorporates four distinct node feature categories, including 768-dimensional BERT features and 300-dimensional spaCy features, which are derived using pre-trained BERT and spaCy word2vec models, respectively. Additionally, a 10-dimensional profile feature is extracted from individual Twitter accounts' profiles. Each graph within the dataset is characterized by a hierarchical tree structure, with the root node representing the news item and the leaf nodes representing Twitter users who have retweeted said news.
An edge exists between a user node and the news node if the user retweeted the original news tweet, while an edge between two user nodes is established when one user retweets the news tweet from another user. This hierarchical organization facilitates the analysis of the spread and influence of both genuine and fabricated news within the Twitter ecosystem. The datasets statistics are summarized in <ref>.

§.§ Implementation Details

Our FROND framework adheres to the experimental settings of the foundational integer-order continuous GNNs, diverging only in the introduction of fractional derivatives in place of integer derivatives. In implementing FROND, we employ one fully-connected (FC) layer on the raw input features to obtain the initial node representations $\bX(0)$. Subsequently, we utilize another FC layer as the decoder function to process the FDE output, $\bX(T)$, for executing downstream tasks.
For more detailed information regarding the hyperparameter settings, we kindly direct the readers to the accompanying supplementary material, which includes the provided code for reproducibility. Our experiments were conducted using NVIDIA RTX A5000 graphics cards.

§.§ Large scale Ogbn-products dataset

In this section, we extend our evaluation to include another large-scale dataset, Ogbn-products, adhering to the experimental settings outlined in <cit.>. For effective handling of this large dataset, we employ a mini-batch training approach, which involves sampling nodes and constructing subgraphs, as proposed by GraphSAINT <cit.>.
Upon examination, we observe that F-GRAND-l demonstrates superior performance compared to both GRAND-l and the GCN model, although it falls slightly short of the performance exhibited by GraphSAGE.
This outcome could potentially be attributed to the insufficient dynamic setting in <ref>. As such, the more advanced dynamic $\calF(\bW,\bX(t))$ in <ref> may require additional refinement.

Node classification accuracy(%) on Ogbn-products dataset

Model      MLP     Node2vec     Full-batch GCN     GraphSAGE     GRAND-l     F-GRAND-l

Acc       61.06$\pm$0.08     72.49$\pm$0.10     75.64$\pm$0.21     78.29$\pm$0.16     75.56$\pm$0.67     77.25$\pm$0.62

§.§ Performance of Different Solver Variants

In this work, we introduce two types of solvers with distinct variants. We evaluate the performance of these variants in <ref>. Specifically, we run F-GRAND on the Cora and Airport datasets with $h=1$ and $T=64$. The solver variants perform comparably. For the Cora dataset, the fractional Adams–Bashforth–Moulton method with a short memory parameter of $K=10$ performs slightly worse than the other variants. However, it demonstrates comparable performance to other solver variants on the Airport dataset.

Node classification accuracy(%) under different solver when time $T=64$

Predictor<ref>     Predictor-Corrector <ref>     Short Memory     Implicit L1

Cora($\beta=0.6$)      83.44$\pm$0.91       83.45$\pm$1.09      81.51$\pm$1.07     82.85$\pm$1.08

Airport($\beta=0.1$)      97.41$\pm$0.42        96.85$\pm$0.36       97.23$\pm$0.59      96.06$\pm$1.59

Node classification accuracy based on memory $K$ on the Cora dataset when time $T=40$.

memory  $K$     1     5     10     15     20     25     30      35     40

Accuracy (%)      74.9$\pm$0.8      80.8$\pm$0.8     83.3$\pm$1.1     83.9$\pm$1.2      84.2$\pm$1.1     84.1$\pm$1.2      84.5$\pm$1.1      84.1$\pm$1.1     84.8$\pm$1.1

Inference (ms)      9.81       17.53       24.97       32.03     38.79      42.99      45.27     48.70     48.35

§.§.§ Further Clarification On Two Accuracies

This section aims to clarify potential ambiguities surrounding the term “accuracy” by distinguishing between “task accuracy” and “numerical accuracy.” Task accuracy pertains to the performance of GNNs on tasks such as node classification. In contrast, numerical accuracy relates to the precision of numerical solutions to FDEs, a critical concern in mathematics.

For example, generally, a larger $K$ value in the Short Memory solver might enhance both numerical and GNN task accuracy. However, it comes with the trade-off of demanding more computational resources. Furthermore, the two accuracies are related, but not equivalent to each other.
For added clarity, we conducted an ablation study on the Cora dataset, keeping all parameters constant except for the memory parameter $K$. The outcomes of this study are detailed in <ref>. Our observations indicate that while increasing the value of $K$ can improve numerical accuracy and potentially GNN task accuracy, the computational cost also rises. Notably, the gains in task accuracy plateau beyond a $K$ value of 15.

We also remind the readers that in the literature, to solve FDEs, there exist other more numerically accurate solvers like <cit.> that use higher convergence order. In general, these kinds of solvers can theoretically reduce computation cost and memory storage, as we can obtain the same numerical accuracy using larger step sizes compared to lower-order solvers. It does not aim to improve GNN task accuracy as we can take smaller step sizes to achieve this, but it may be helpful for other performances like computation cost and memory storage reduction. In our paper, we focus on task accuracy. Therefore, classical solvers are used in our work.
Nonetheless, more numerically accurate solvers could potentially benefit other applications of fractional dynamics, particularly when GNNs are utilized to simulate and forecast real physical systems.

§.§ Computation Time

It should be emphasized that our FROND framework does not introduce any additional training parameters to the backbone integer-order continuous GNNs. Instead, we simply modify the integration method from standard integration to fractional integration.

In this section, we report the inference time of the different solver variants in <ref>. For comparison, we consider the neural ODE solver for $\beta=1$, which includes Euler, RK4, Implicit Adams, and dopri5 methods as per in the paper <cit.>. We observe that when $T=4$, the inference time required by the FROND solver variants is similar to that of the ODE Euler solver. However, for larger $T=64$, the basic Predictor <ref> solver requires more inference time than Euler and is comparable to RK4. For more accurate approximation solver variants <ref> incorporating the corrector formula, <ref> show that these methods require more computational time as the number of iterations increases. While the advantages of these solvers might not be pronounced for GNN node classification tasks, they could provide benefits for other applications of fractional dynamics, such as when GNNs are used to simulate and forecast real physical systems.

Average time under different solvers when time $T=4$ and hidden dimension is 64 on Cora dataset

Predictor<ref>     Predictor-Corrector<ref>     Short Memory     Implicit L1     Euler     RK4     Implicit Adams     dopri5

Inference time (ms)     0.98     1.67     0.98     0.62     0.96     2.06     3.20      11.91

Average time under different solvers when time $T=64$ and hidden dimension is 64 on Cora dataset

Predictor<ref>     Predictor-Corrector<ref>     Short Memory     Implicit L1     Euler     RK4     Implicit Adams     dopri5

Inference time (ms)     44.46     160.92      30.26      221.74      12.16     42.66      103.46      66.15

Average time of <ref> with correctors, used to refine the approximation, when time $T=4$ and hidden dimension is 64 on the Cora dataset.

Predictor-Corrector <ref>       1      3      5     10

Inference time (ms)     1.67     3.31      4.74      8.34

Implicit-L1 <ref>       1      3      5     10

Inference time (ms)     0.62     1.04     1.48      2.55

Average time of <ref> with correctors, used to refine the approximation, when time $T=64$ and hidden dimension is 64 on the Cora dataset.

Predictor-Corrector <ref>       1      3

Inference time (ms)     160.92     442.88

Implicit-L1 <ref>       1      3

Inference time (ms)     221.74       441.60

§.§ Continued Study of Oversmoothing

To corroborate that FROND mitigates the issue of oversmoothing and performs well with an increasing number of layers, we conducted an experiment employing the basic predictor with up to 128 layers in the main paper. The results are presented in <ref>. For this experiment, we utilized the fixed data splitting approach for the Cora and Citeseer dataset without using the Largest Connected Component (LCC)  as described in <cit.>.

In the supplementary material, we further probe oversmoothing by conducting experiments with an increased number of layers, reaching up to 256. The results of these experiments are illustrated in <ref>. From our observations, F-GRAND-l maintains a consistent performance level even as the number of layers escalates.
This contrasts with GRAND-l, where there is a notable performance decrease with the increase in layers. For instance, on the Cora datasets, the accuracy of GRAND-l drops from 81.29% with 4 layers to 73.37% with 256 layers. In stark contrast, our F-GRAND-l model exhibits minimal performance decrease on this dataset. On the Airport dataset, F-GRAND-l registers a slight decrease to 94.91% with 256 layers from 97.0% with 4 layers. However, the performance of GRAND-l significantly drops to 53.0%. These observations align with our expectations, as <ref> predicts a slow algebraic convergence rate, while GRAND exhibits a more rapid performance degradation.

Additionally, we note that the optimal number of layers for F-GRAND is 64 on the Cora and Airport datasets, whereas on the Cirtesser dataset, the best performance is achieved with 16 layers.

oversmoothing mitigation under fixed data splitting without LCC

Dataset     Model      4     8     16     32     64      80     128    256

GCN      81.35$\pm$1.27     15.3$\pm$3.63     19.70$\pm$7.06     21.86$\pm$6.09     13.0$\pm$0.0     13.0$\pm$0.0     13.0$\pm$0.0      13.0$\pm$0.0

GAT      80.95$\pm$2.28     31.90$\pm$0.0     31.90$\pm$0.0     31.90$\pm$0.0     31.90$\pm$0.0    31.90$\pm$0.0     31.90$\pm$0.0      31.90$\pm$0.0

GRAND-l      81.29$\pm$0.43     82.95$\pm$0.52     82.48$\pm$0.46     81.72$\pm$0.35     81.33$\pm$0.22     81.07$\pm$0.44     80.09$\pm$0.43     73.37$\pm$0.59

F-GRAND-l     81.17$\pm$0.75     82.68$\pm$0.64     83.05$\pm$0.81     82.90$\pm$0.81     83.44$\pm$0.91     82.85$\pm$0.89     82.34$\pm$0.83

GCN      68.84$\pm$2.46     61.58$\pm$2.09     10.64$\pm$1.79     7.7$\pm$0.0     7.7$\pm$0.0     7.7$\pm$0.0     7.7$\pm$0.0     7.7$\pm$0.0

GAT      65.20$\pm$0.57     18.10$\pm$0.0      18.10$\pm$0.0     18.10$\pm$0.0    18.10$\pm$0.0     18.10$\pm$0.0     18.10$\pm$0.0     18.10$\pm$0.0

GRAND-l        70.68$\pm$1.23     70.39$\pm$0.68     70.18$\pm$0.56     68.90$\pm$1.50     68.01$\pm$1.47      67.44$\pm$1.25      63.45$\pm$2.86     56.98$\pm$1.26

F-GRAND-l       70.68$\pm$1.23     71.04$\pm$0.68     71.08$\pm$1.12     70.83$\pm$0.90     70.27$\pm$0.86     70.50$\pm$0.76     70.32$\pm$1.67      71.0$\pm$0.45

GCN      84.77$\pm$1.45     74.43$\pm$8.19     62.56$\pm$2.16     15.27$\pm$0.0     15.27$\pm$0.0     15.27$\pm$0.0      15.27$\pm$0.0      15.27$\pm$0.0

GAT      83.59$\pm$1.51     67.02$\pm$4.70     46.56$\pm$0.0     46.56$\pm$0.0     46.56$\pm$0.0     46.56$\pm$0.0     46.56$\pm$0.0     46.56$\pm$0.0

GRAND-l     80.53$\pm$9.59      79.88$\pm$9.67     76.24$\pm$3.80     68.67$\pm$4.02     62.28$\pm$10.83      50.38$\pm$2.98     57.96$\pm$11.63     53.0$\pm$14.85

F-GRAND-l     97.0$\pm$0.79     97.09$\pm$0.87      96.97$\pm$0.84     96.50$\pm$0.60     97.41$\pm$0.42     96.53$\pm$0.74     97.03$\pm$0.55     94.91$\pm$3.72

§.§ Ablation Study: Selection of $\beta$ Continued

In the main paper, we explore the impact of the fractional order parameter $\beta$ across a variety of graph datasets, with the results of these investigations presented in <ref>. More comprehensive details concerning the variations in $\beta$ can be found in <ref>.

Node classification accuracy(%) under different value of $\beta$ when time $T=8$.

$\beta$      0.1     0.2     0.3     0.4     0.5     0.6     0.7     0.8     0.9     1.0

Cora     74.80$\pm$0.42     76.10$\pm$0.34     77.0$\pm$0.98      77.80$\pm$0.75      79.60$\pm$0.91      80.79$\pm$0.58       81.56$\pm$0.30        82.44$\pm$0.51       82.68$\pm$0.64      82.37$\pm$0.59

Airport     97.09$\pm$0.87       96.67$\pm$0.91     95.80$\pm$2.03     94.04$\pm$3.62     91.66$\pm$6.34     89.24$\pm$7.87     84.36$\pm$8.04     79.29$\pm$6.01      78.73$\pm$6.33      78.88$\pm$9.67

§.§ Robustness Against Adversarial Attacks

Despite the significant advancements GNNs have made in inference tasks on graph-structured data, they are recognized as being susceptible to adversarial attacks <cit.>. Adversaries, aiming to deceive a trained GNN, can either introduce new nodes into the graph during the inference phase, known as an injection attack <cit.>, or manipulate the graph's topology by adding or removing edges, termed as a modification attack <cit.>.
In this section, we present preliminary experiments assessing the robustness of our model against adversarial attacks. Specifically, we carry out graph modification adversarial attacks using the Metattack method <cit.>. Our approach adheres to the attack setting described in Pro-GNN <cit.>, and we utilize the perturbed graph provided by the DeepRobust library <cit.> to ensure a fair comparison. The perturbation rate, indicating the proportion of altered edges, is incrementally adjusted in 5% steps from 0% to 25%.

The results of these experiments are presented in <ref>.
It should be noted that the impact of Meta-attacks with higher strengths detrimentally affects the performance of all models under test. However, our FROND-nl model consistently demonstrates enhanced resilience against adversarial attacks compared to the baselines, including GRAND-nl.
For instance, at a perturbation rate of 25%, F-GRAND-nl outshines the baselines by an estimated margin of 10-15% on the Cora dataset.

Comprehensive testing against various adversarial attack methods and a theoretical understanding are detailed in our recent work <cit.>.

Node classification accuracy (%) under modification, poisoning, non-targeted  attack (Metattack) in transductive learning.

Dataset     Ptb Rate(%)       GGN     GAT     GRAND-nl      F-GRAND-nl

0      83.50$\pm$0.44      83.97$\pm$0.65     83.14$\pm$1.06     83.48$\pm$1.08

5      76.55$\pm$0.79      80.44$\pm$0.74     80.54$\pm$1.17      80.25$\pm$0.90

10      70.39$\pm$1.28      75.61$\pm$0.59     76.59$\pm$1.21     77.94$\pm$0.48

15      65.10$\pm$0.71     69.78$\pm$1.28     71.62$\pm$1.39      75.14$\pm$1.16

20      59.56$\pm$2.72     59.94$\pm$0.92     57.52$\pm$1.20      69.04$\pm$1.13

25      47.53$\pm$1.96     54.78$\pm$0.74     53.70$\pm$1.91      63.40$\pm$1.44

0      71.96$\pm$0.55     73.26$\pm$0.83     71.40$\pm$1.08      70.14$\pm$0.83

5      70.88$\pm$0.62     72.89$\pm$0.83     70.99$\pm$1.12      70.0$\pm$1.72

10      67.55$\pm$0.89     70.63$\pm$0.48     68.83$\pm$1.31      68.64$\pm$1.11

15      64.52$\pm$1.11      69.02$\pm$1.09     66.78$\pm$0.92      67.90$\pm$0.41

20      62.03$\pm$3.49     61.04$\pm$1.52     58.95$\pm$1.33      65.84$\pm$0.75

25      56.94$\pm$2.09     61.85$\pm$1.12      60.52$\pm$1.29      66.50$\pm$1.16

§.§ Comparison between Riemann-Liouville (RL) derivative and Caputo Derivative

The underlying rationale for opting for the Caputo derivative over the Riemann-Liouville (RL) derivative is extensively delineated in <ref>. However, a supplementary experiment was conducted utilizing the RL derivative in lieu of the Caputo derivative, the results of which are documented in  <ref>.
It can be observed that the task accuracies for both approaches are very similar.
Further investigations on the use of different fractional derivatives and how to optimize the whole model architecture to adapt to a particular choice will be explored in future work.

Comparison between RL-GRAND-l (using Riemann-Liouville derivative) and the original F-GRAND-l (using Caputo derivative).

Method     Cora     Citeseer     Pubmed      CoauthorCS      Computer       Photo       CoauthorPhy      Airport     Disease

GRAND-l     83.6$\pm$1.0     73.4$\pm$0.5     78.8$\pm$1.7      92.9$\pm$0.4         83.7$\pm$1.2      92.3$\pm$0.9        93.5$\pm$0.9     80.5$\pm$9.6     74.5$\pm$3.4

RL-GRAND-l     84.6$\pm$1.2     74.2$\pm$1.0     80.1$\pm$1.2     92.8$\pm$0.3      87.4$\pm$1.1      93.3$\pm$0.7      94.1$\pm$0.3     96.2$\pm$0.2     90.7$\pm$1.3

F-GRAND-l     84.8$\pm$1.1      74.0$\pm$1.5      79.4$\pm$1.5     93.0$\pm$0.3        84.4$\pm$1.5      92.8$\pm$0.6        94.5$\pm$0.4     98.1$\pm$0.2      92.4$\pm$3.9

§.§ Fractal Dimension of Graph Datasets

Comparison between  the estimated fractal dimension, the best order $\beta$ and the $\delta$-hyperbolicity

Dataset     Disease     Airport     Pubmed     Citeseer     Cora

fractal dimension      2.47     2.17     2.25     0.62     1.22

best $\beta$ (F-GRAND-l)    0.6     0.5     0.9      0.9    0.9

best $\beta$ (F-GRAND-nl)     0.7    0.1     0.4      0.9    0.9

$\delta$-hyperbolicity      0.0    1.0    3.5    4.5    11.0

In <ref>, using the Compact-Box-Burning algorithm from <cit.>, we compute the fractal dimension for some datasets that have moderate sizes. As noted in <ref>, there is a clear trend between $\delta$-hyperbolicity (as referenced in <cit.> for assessing tree-like structures—with lower values suggesting more tree-like graphs) and the fractal dimension of datasets. Specifically, a lower $\delta$-hyperbolicity corresponds to a larger fractal dimension. As discussed in <ref>, we believe that our fractional derivative $D_t^\beta$ effectively captures the fractal geometry in the datasets.
Notably, we discerned a trend: a larger fractal dimension typically corresponds to a smaller optimal $\beta$.

The fractal dim of datasets. We use the Compact-Box-Burning algorithm in <cit.> to compute the log-log slope (fractal dim) of the box size and the minimum number of boxes needed to cover the graph.

§ MORE DYNAMICS IN FROND FRAMEWORK

§.§ Review of Graph ODE Models

GRAND++: The work by <cit.> introduces graph neural diffusion with a source term, aimed at graph learning in scenarios with a limited quantity of labeled nodes. This approach leverages a subset of feature vectors, those associated with labeled nodes, indexed by $\calI$, and considered “trustworthy” to act as a source term. It adheres to <ref>, incorporating an additional source term, facilitating the propagation of information from nodes in $\calI$ to node $i$.

\begin{align}
\frac{\ud \bX(t)}{\ud  t} =F(\bX(t))+s(\{\bx_{i}\}_{i\in\calI})
\end{align}

Here, $\calI$ denotes the set of source nodes, $s(\cdot)$ represents a source function, and $F(\cdot)$ embodies the function depicting the right-hand side of <ref>. The model is manifested in two variations, respectively denoted as GRAND++-nl and GRAND++-l.

GraphCON: Inspired by oscillator dynamical systems, GraphCON <cit.> is defined through the employment of second-order ODEs. It is crucial to highlight that, for computation, the second-order ODE is decomposed into two first-order ODEs:

\begin{align}
\begin{aligned}
& \frac{{\ud} \mathbf{Y}(t)}{{\ud} t} = \sigma ({\mathbf{F}_{\theta}(\mathbf{X}(t), t)) - \gamma\mathbf{X}(t) -\tilde{\alpha}\mathbf{Y}(t)}, \
& \frac{{\ud} \mathbf{X}(t)}{{\ud} t} = \mathbf{Y}(t), \label{eq.graphcon_diff}
\end{aligned}
\end{align}

where $\sigma(\cdot)$ is the activation function, $\mathbf{F}_{\theta}(\mathbf{X}(t), t)$ is the neural network function with parameters $\theta$, $\gamma$ and $\tilde{\alpha}$ are learnable coefficients, and $\mathbf{Y}(t)$ is the velocity term converting the second-order ODE to two first-order ODEs.

Analogous to the GRAND model, the GraphCON model is also available in both linear (GraphCON-l) and non-linear (GraphCON-nl) versions concerning time. The differentiation between these versions is determined by whether the function $\mathbf{F}_{\theta}$ undergoes updates based on time $t$.

CDE: With the objective of addressing heterophilic graphs, the paper <cit.> integrates the concept of convection-diffusion equations (CDE) into GNNs, leading to the proposition of the neural CDE model: This innovative model incorporates a convection term and introduces a unique velocity for each node, aiming to preserve diversity in heterophilic graphs. The corresponding formula is illustrated in <ref>.

\begin{align}
\frac{\ud \bX(t)}{\ud t} =(\mathbf{A}(\mathbf{X}(t))-\mathbf{I}) \mathbf{X}(t) + \operatorname{div}(\bV(t) \circ \bX(t))
\label{eq.cde_diff}
\end{align}

In this equation, $\bV(t)$ represents the velocity field of the graph at time $t$, $\operatorname{div}(\cdot)$ denotes the divergence operator as defined in the paper <cit.>, and $\circ$ symbolizes the element-wise (Hadamard) product.

To address the challenges posed by heterophilic graphs, the authors in <cit.> present the GREAD model. This model enhances the GRAND model by incorporating a reaction term, thereby formulating a diffusion-reaction equation within GNNs. The respective formula is depicted in <ref>, and the paper offers various alternatives for the reaction term.

\begin{align}
\frac{\ud \bX(t)}{\ud t} = -\alpha \bL(\bX(t)) + \alpha r(\bX(t))
\label{eq.gread_diff}
\end{align}

In this equation, $r(\bX(t))$ represents the reaction term, and $\alpha$ is a trainable parameter used to balance the impact of each term.

§.§ F-GRAND++

Building upon the GRAND++ model <cit.>, we define F-GRAND++ as follows:

\begin{align}