Unleashing the Potential of Fractional Calculus in Graph Neural Networks with FROND
Authors: Qiyu Kang
Kai Zhao
Qinxu Ding
Feng Ji
Xuhao Li
Wenfei Liang
Yang Song
Wee Peng Tay
Venue: ICLR 2024 spotlight (2024)
Keywords: graph neural network
Primary area: learning on graphs and other geometries & topologies

Abstract: 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 \url{https://github.com/zknus/ICLR2024-FROND}.

Published as a conference paper at ICLR 2024
UNLEASHING THE POTENTIAL OF FRACTIONAL CAL- CULUS IN GRAPH NEURAL NETWORKS WITH FROND
Qiyu Kang!*; Kai Zhao'; Qinxu Ding", Feng Ji', Xuhao Li Wenfei Liang', Yang Song".
Wee Peng Tay
Nanyang Technological University 'Singapore University of Social Sciences
3 AAnhui University #C3 AI, Singapore
ABSTRACT
We introduce the FRactional-Orde 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 em- ploys the Caputo fractional derivative to leverage the non-local properties of fractional calculus. This approach enables the capture of long-term dependen- cies 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 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 consis- tently improved performance and underscoring the framework's s potential effective extension to enhance traditional continuous GNNs. The code available at https: http//gtttutuuu com/ zknuz DAtEnn//IIIR2222----1 FROND.
INTRODUCTION
Graph Neural Networks (GNNs) have excelled in diverse domains, e.g., chemistry (Yue et al., 2019). finance (Ashoor et al., 2020), and social media (Kipf & Welling, 2017; Zhang et al., 2022; Wu et al., 2021). The message passing scheme (Feng et al., 2022), where features are aggregated along edges and iteratively propagated through layers, crucial for the success of GNNs. Over the past few years, numerous types of GNNs have been proposed, including Graph Convolutional Networks (GCN) (Kipf & Welling, 2017). Graph Attention Networks (GAT) (Veličković et al., 2018), and GraphSAGE (Hamilton et al., 2017). Recent works, such (Chamberlain et al., 2021c; Thorpe et al., 2022; Rusch et al., 2022; Song et al., 2022; Choi et al., 2023; Zhao et al., 2023a; Kang et al., 2023), 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 (Han et al., 2023).
Within these integer-orde continuous GNNs, the differential operator d3 / dt has been constrained to integer values of 3, primarily or 2. However, over recent decades, the wider scientific community has explored fractional-ordes differential operators, where /3 can be any real number. These expansions have proven pivotal in various applications characterized by non-local and memory-depencee behaviors, with examples including viscoelastic materials (Bagley & Torvik, 1983), anomalous transport mechanisms (Gómez- 1982). Unlike conventional integer-ordes derivatives that measure the function's 's instantaneous rate of change and focus on the local vicinity, fractional-ordes derivatives (Tarasov. 2011 1) consider the entire historical trajectory of the function.
et al., 2016), and fractal media (Mandelbrot & Mandelbrot,
First two authors contributed equally. Cooresponndence to: Qiyu Kang wttp:/000000......uuu........

---

Published as a conference paper at ICLR 2024
We introduce the FRactional- Order graph Neural Dynamical network (FROND) framework, new approach that broadens the capabilities of traditional integer-order continuous GNNs by incorporating fractional calculus. It naturally generalizes the integer-order derivative d dt in these GNNs to accommodate any positive real number /3. This modification gives FROND the ability to incorporate memory--eepeneen dynamics for information propagation and feature updating, enabling refined graph representations and improved performance potentially. Importandy, this technique assures at least equivalent performance to integer-order models, as setting & to integer values reverts the models to their traditional integer-order forms.
Several works like (Maskey et al., 2023) 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 fractfractional derivatives for updating graph node features, modeled as memory-inclusive dynamical process. Other works like (Liu et al., 2022) 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 Appendix A.
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 pervasive simiasitllrrtty across different scales, hinting at underlying fractal behavior (Song et al., 2005; Kim et al., 2007; Masters, 2004). It has been well-established that dynamical processes with self-similarity dispersion of heat or mass over these structures is best modeled using fractional diffusion equations (Diaz-Diaz & Estrada, 2022). Further investigations have revealed a direct connection between the fractal dimension of these structures and the order 3 in fractional derivatives d dt (Nigmatullin, 1992; Tarasov, 2011). This revelation births a compelling insight: the optimal 33 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.
such fractal media are more accurately described using FDEs. For instance, the
Main contributions. Our objective in this paper is to formulate generalized tractional-orde continuous GNN framework. Our key contributions are summarized as follows:
We propose a novel, generalized continuous GNN framework that incorporates non-local fractional derivatives d / dtP_ This framework generalizes the prior class of integer-order continuous GNNs, subsuming them as special instances with 3 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 fcature-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 exponentially swift, we establish that in FROND, convergence follows 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 (Chamberlain et al., 2021c; Thorpe et al., 2022; Rusch et al., 2022; Song et al., 2022: Choi et al., 2023; Zhao et al., 2023a), substantiate this claim. Through detailed ablation studies, we provide insights into the choice of numerical schemes and parameters.
2 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 Appendix B.
2.1 CAPUTO FRACTIONAL DERIVATIVE
The literature offers various fractional derivative definitions, notably by Riemann, Liouville, Chapman, and Caputo (Tarasov, 2011). Our study leverages the Caputo fractional derivative, due to the reasons listed in Appendix B.4. The traditional first-order derivative of scalar function ft) represents the
2

---

Published as a conference paper at ICLR 2024
local rate of change of the function at point, defined as: df(t)
dt lima f(+ Let F(s)
denote the Laplace transform of ft), assumed to exist on [Su. o0) for some so E R. Under certain
conditions (Korn & Korn, 2000), the Laplace transform of df(1 is given by:
dt
{{ df((t) }
c{
dt J - =F(s) f(0)
(1)
The Caputo fractional derivative of order 3 E (0, 1] for a function f(t) is defined as follows:
D r(1 3) (f-
where r(-) denotes the gamma function, and (r( is the first-order derivative of f. The broader definition for any 3 > deferred to Appendix B. The Caputo fractional derivative inherently integrates the entire history of the system through the integral term, emphasizing its non-local nature. For 8 > max {0, s.), the Laplace transform of the Caputo fractional derivative is given by (Diethelm, [010)[Theorem 7.1]:
c = Comparing (1) and (3), it is evident that the Caputo derivative serves as generalization of the first-order derivative. The alteration in the exponent of comes from the memory memory-dependent property in (2). As 3 - 1, the Laplace transform of the Caputo fractional derivative converges to that of the traditional first-order derivative. When 3 = 1. Dif Laplace transform (Cohen, 2007).
is uniquely determined through the inverse
In summary, from the frequency domain using the Laplace transform, we observe that the Caputo fractional derivative can be vector-valued functions, the fractional derivative is defined component-wis for each dimension.
natural extension of the traditional first-order derivative. For
2.2 STTEGER-ORDEE CONTINUOUS GNNS
We denote undirected graph G (V, W) without self-loops, where V is the set of |/ N nodes. The feature matrix X = ([x1]" TT ... [xxw])" E R Nxd consists of rows x; E Rd as node feature vectors and is the node index. The N N matrix W == (W,,) has elements W., indicating the edge weight between the i-th and j-th node with Wa; ij Wze. The following integer-order continuous GNNs leverage ODEs to facilitate information propagation amongst graph nodes, where features evolve as (t), starting from the initial condition X(0)= = X.
GRAND: Inspired by the heat diffusion equation, GRAND (Chamberlain et al., 2021c) utilizes the following nonlinear autonomous dynamical system:
dx(() = (A(X(I))I)))) - di
where A(X(t) E NXN is a learnable, time- variant attention matrix, calculated using the features X(t,, and I denotes the identity matrix. The feature update outlined in (4) is referred to as the GRAND-nl version (due to the nonlinearity in A(X(t))). We define d; = ΣØ= Way and let D be a diagonal matrix with Di d,. The random walk Laplacian is then represented as L = I WD-1
In a simplified context, we employ the following linear dynamical system:
dX(() (WD I)x(() LX)
dt
The feature update process in (5) is the GRAND-I version. For implementations of (5), one may direct set WD-1 (X)( as column-stochasse 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++ (Thorpe et al., 2022) adds source term to GRAND, enhancing learning in scenarios with limited labeled nodes. GraphCON (Rusch et al., 2022) employs a second-order ODE, which is equivalent to two first-order ODEs, drawing inspiration from oscillator systems. CDE (Zhao et al., 2023a) incorporates convection-diffusion equations into GNNs to address heterophilic graph challenges. GREAD (Choi et al., 2023) introduces reaction term in the GRAND model, improving its application to heterophilic graphs and formulating a diffusion-reactinn The detailed formulation for each model is presented in Appendix E. due to space constraints.
equation within GNNs.
3

---

Published as a conference paper at ICLR 2024
3 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 coun- terparts of several well-cstablished integer-ordes continuous GNNs, including GRAND, GRAND++, GraphCON, CDE, and GREAD, referenced in Section 2.2. We provide detailed study of the fractional extension of GRAND, and present insights into the inherent memory mechanisms in our framework through random walk interpretation. Our theoretical findings suggest potential mitiga- tion of oversmoothing due to the model's slow algebraic convergence to stationarity. Subsequently, we outline the numerical FDE solvers required implement FROND.
3.1 FRAMEWORK
Consider a graph 9 = (V. W) as defined in Section 2.2. Analogous to the implementation in traditional integer-order continuous GNNs, a preliminary learnable encoder function y V + Rd that maps each node to feature vector can be applied. Stacking all the feature vectors together, we obtain X E RNNx xd Employing the Caputo fractional derivative outlined in Section 2.1, the information propagation and feature updating dynamics in FROND are characterized by the following FDE:
Dix(1)= = F(W, X()), 8>0 0,
(6)
where 3 denotes the fractional order of the derivative, and F is a dynamic operator on the graph like the
models presented in Section 2.2. The initial condition for (6) is set as xl( (0) = =... x()) X
consisting of the preliminary node features, with (t) denoting the i-th order derivative and [3] is
the smallest integer not less than 3, akin to the initial conditions to integer-order continuous GNNs, we set an integration time parameter T to get X The final node embeddings for downstream tasks are then decoded using learnable decoder o6 (X()).
in integer-order ODEs. Similar
When B = 1, (6) reverts to the class of integer-order continuous GNNs, with the infinitesimal variation of features dependent only on their present state. Conversely, when < 1. the Caputo fractional derivative (2) dictates that the updating process for features encompasses their entire history, not just the present state. This paradigm facilitates memory deeeeddd dynamics in the framework.
For further insights into memory dependence, readers are directed to Section 3.3, 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 Section 3.2, when the dynamic operator F is designated the diffusion process in (5), we offer a memory-dependeet non-Markovian random walk interpretation of the fractional graph neural diffusion process. Here, 3 + 1, the Mon-Markovian random walk increasingly detaches from the path history, becoming a Markovian walk at 3 = 1, which is related to the normal diffusion process (Thorpe et al., 2022). The parameter B provides flexibility to adjust the extent of memorized dynamics embedded in the framework. From a geometric perspective, as discussed in Section 1, the information propagation dynamics in fractal graph datasets might be more suitably described using FDES. Choosing a non-integer 3 could reveal the degree of fractality in graph datasets.
3.1.1 FROND MODEL EXAMPLES
When F in (6) is specified to the dynamics depicted in various integer-order continuous GNNs (cf. Section 2.2), we formulate FROND GNN variants such 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-GRAD (F-GRAND) has two versions.
The F-GRAND-nl version employs a time-variant FDE as follows:
D - (A(X()) I)) 0 0<351.
See Appendix B.3.2, We mainly consider / E (0. 1] and the initial condition is X(0) X.
4

---

Published as a conference paper at ICLR 2024
It is computed using (t) and the attention mechanism derived from the Transformer model (Vaswani et al., 2017). The entries of (t)(t) = (a)x, are given by:
a(x,,3, x,) softmax
(8)
In this formulation, W and We are the learned matrices, and de signifies a hyperparameter related to the dimensionality of W K- In parallel, the F-GRAND-I version stands as the fractional differential extension of (5):
x( = (L)((), 0 0<BS1.
Recall that the initial condition for F-GRAND-nl and F-GRAND-I is X - X due to 3 E (0.
F-GRAND++, F-GREAD, F-CDE, and F-GraphCON: Due to space constraints, direct the reader to Appendix E for detailed formulations. Suceinctly, 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-ordet continuous GNNs across a variety of datasets, exhaustive experiments are provided in Section 4 and Appendix E.
3.2 RANDOM WALK PERSPECTIVE OF F-GRAND-L
The established Markov interpretation of GRAND-I (5), as outlined in (Thorpe et al., 2022), aligns with F-GRAND-I (9) when 3 1. We herein broaden this interpretation to encompass a non- Markovian random walk that considers the walker's complete path history when 3 is non-integer, thereby elucidating the memory effects inherent FROND. In contrast to the Markovian walk, whose distribution converges exponentially to equilibrium, our strategy assures algebraic convergence, revealing F-GRAND-I's efficacy in mitigating oversmoothing as evidenced in Section 4.3.
To begin, we discretize the time domain into time instants as ln = no, a > 0, n. 0. 2... where is assumed to be small enough to ensure the validity of the approximation. Let R((,) be a random walk on the graph nodes {x, ح that is not necessarily a Markov process and R(tn+1) may depend on the path history (R(ta) R(11), (fff)........ R(,,)) ) of the random walker. For convenience, we introduce the coefficients for k 2 and b, for n 2 0 from (Gorenflo et al., 2002), which are used later to define the random walk transition probability:
<a (3) (-1)
1 K
(
7 Σ (-)*
B 1:
(10)
k=g
where the generalized binomial coefficient employed in the definition of the coefficients. The sequences cj and b, consist of positive numbers, not greater than 1, decreasing strictly monotonically zero (see supplementary material for details) and satisfy ==== Cs + bn = 1. Using these coefficients, we define the transition probabilities of the random walk starting from Xjo
r( (T++))(0=++)) and the gamma function r(-) are
T(3+1)
P(R(tn+1) = j
R(ta) Xjo: R(tt) = X),. R(t,)) = =x,,)
C1 0
if staying at current location with fn+1 = jnn
Wininti の dyn
if jumping to neighboring nodes with jn+1 jn: revisiting historical positions with ja+1 = jk, 1 S < 1
(1I)
Cn+1-b bn
if revisiting historical positions with fa+1 = ju.
This formulation integrates memory effects, considering the walker's time, position, and path history. The transition mechanism of the memory- -ncluuive random walk between and t,+1 is elucidated as follows: Suppose the walker at node ja at time Ln having a full path history fi, ..)......
We generate random number p E [0, 1) uniformly, and divide the interval [0, 1) into adjacent sub-intervals with lengths C2........ Ca, big. We further subdivide the first interval (with length c1)
into sub-intervals of lengths e " and "
1. If p is in the first interval with length c1, the walker either moves to a neighbor jut1 k with probability 3 Wiah
remains the current position with probability e₁ g².
dyn
2. For p in subsequent intervals, the walker jumps to a previously visited node the history (jo: fi.... - In specifically, ja+1 if in Cg. jo if in by.
5

---

Published as a conference paper at ICLR 2024
When 3 < 1, the random walk can, with positive probability, revisit its history, restricting exten- sive drift. We denote ((R(t,)) as the probability column vector, with its j-th element given as P(R(t,,) = x,). Additionally, we specify ,P ,P(R(t,)) to indicate the situation where the random walker initiates from the i-th node, i.e., R()) = X;, with probability 1. In this case, the initial one-hot vector with the i-th entry marked 1. Using the technique from (Gorenflo et al., 2002), we can prove the following:
probability vector; ,PP(R (0)) is represented
Theorem 1. Consider the random walk defined in (11), with the step size and number of steps n. Under the conditions that n + oo and no = t. the limiting probability distribution P(t) = oom P(R(t,)) satisfies (9). In other words.
D. -PP())-L(())
(12)
Considering that initial conditions and dimensions affect the solutions of FDEs, P(t) and ((t are not equivalent. However, due to the linearity of FDEs, the following conclusion is straightforward: Corollary 1. Under the conditions that n + 0o and no = t, we have lim, Σ. (,,,(((,,)))... X(t), ie.. Ex P())x; = X(t with (() == lim ,P ((((tt)) where X(t) is the solution to (9) with the initial condition (() X.
Remark 1. Theorem I and Corollary I relate F-GAND-I (9) to the non-Markovian random walk in (11), illustrating memory dependence in FROND. As B + 1, this process reverts to the Marko- vian random walk found in GRAND-I (Thorpe et al., 2022) in (13). This underscores the FROND
framework's 's capability to apprehend more complex dynamics than integer-order continuous GNNs.
P 1((((1+1) = Xjn+1
(to) = Xjo: R(t)) = Xi......
R(,,) x,.)
(13)
1
if staying at current location with In+11 In
= =I(((( = Xjn+i R(tn)= = xzn.)
Winint1 dym
juupping to neighbors with Jn+1 + In
since we have that all these coefficients vanishing except 01 = 1, i.e. (B) 0,
=1 lim
#22 2,
lim (g) =0, 0, 12
(14)
8-1
8-1
3.2.1 OVERSMOOTHING MITIGATION OF F-GRAND-I. COMPARED TO GRAND-L
The seminal research (Oono & Suzuki, 2020)(orrlllary 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 (13) is given by π ( (g EX 1 بو "3 (Thorpe et al., 2022), with an exponentially
dx
O(e-r'n) where r' >> O relates to the eigenvalues of the matrix L (Chung, 1997), and || ll2 denotes the E² norm. This behavior extends to the continuous limit, akin first-order linear ODE solution, exhibiting exponential convergence with some > 0:
rapid convergence rate (I.(R(())) -
[P(t)= - π 122 ~ O(e"").
(15)
In contrast, we next prove that the non- Mar-Markoiinn random walk (11) converges to the stationary
distribution at a slow algebraic rate, thereby helping to mitigate oversmoothing. As 3 + 0, the convergence is expected to be arbitrarily slow. In real-world scenarios where we operate within finite horizon, this slower rate of convergence may be sufficient to alleviate oversmoothing, particularly when it is imperative for deep model to extract distinctive features instead of achieving exponentially fast convergence to a stationary equilibrium.
Theorem 2. Under the assumption that the graph is strongly connected and aperiodic, the stationary probability for the non-Markovian random walk (11), with 0 < 3 < 1, is still πT, which is unique, This mirrors the stationary probability of the Markovian random walk as defined by (13) when 3 == 1. Notably, when B < 1, the convergence of the distribution (distinct from π)) to π is algebraic:
[PP(I) c 112 ((--)).
(16)
= - E, as the solution to F- GRAND-I (9), converges E, π TXx; = -- Σx x; at a slow algebraic rate since [l:P(t) T=Il ll2 ~
Remark 2. Corollary I and Theorem 2 indicate that X(t)=
@(L-9) for all i. Notably, π " forms a rank invariant subspace under the dynamics of (9), due to π being stationary, This underscores the difference in convergence rates, contrasting the slow algebraic rate in our case with the fast exponential rate (Oono & Suzuki, 2020; Zhao et al., 2023b)
We use the asymptotic order notations from (Notations, 2023) in this paper.

---

Published as a conference paper at ICLR 2024
3.3 SOLVING FROND
memory window width K
[(0)
Xi6.]
X(ss
Xif..1) Xit.)
X(0)
X(ti)
X(()
X(h=)) X(te)
Figure I: Diagrams of fractional Adams Bashforth Moulton method with full (left) and short (right) memory.
The studies by (Chen et al.. 2018b; Quaglino et al., 2019; Yan et al., 2018) introduce numerical solvers specifically designed for integer-order neural ODE models. Our research, in contrast, engages with fractional-orde ODEs, entities inherently more intricate than integer-order ODEs. To address the scenario where 3 is non-integer, we introduce the fractional explicit Adams-Badd Bashns--saffohh--M solver, incorporating three variants employed in this study: the basic predictor discussed in this section, the predictor-corrector elaborated in Appendix C.2, and the short memory principle detailed in Appendix C.3. Additionally, methods exemplify how time still acts as a continuous analog to the layer index and elucidate how memory dependence manifests as nontrivial dense skip connections between layers (see Figs. 1 and 4), stemming from the non-local properties of fractional derivatives.
present one implicit L1 solver in Appendix C.4. These
Basic Predictor: We first employ preliminary numerical solver called "predictor" (Diethelm et al., 2004) through time discretisation. Let h be a small positive discretization parameter. We have
[Ø] -1
PX
Σ
k x (0) +
1 ((3) Fo,
(17)
where 13,8 h2 -) )))). le denotes the discrete time index (iteration), and
t kh represents the discretized time steps. (k) is the numerical approximation of X(tx). When 3 = 1, this method simplifies to the Fuler solver in (Chen et al., 201 Chamberlain et al., 2021c) as Mj.n. =/ h. yielding p X (k) X (A-
hF((W X (k)))). Thus, our basic predictor can be considered as the fractional Euler method or fractional Adams Bashforthh method, which is generalization of the Euler method used in (Chen et al., 2018b; Chamberlain et al., 2021c). However, when 3 < 1, we need to utilize the full memory (F(W, (F(WX(() } The block diagram in Fig. 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 Fig. 4, elucidating the manner in which information propagates through layers and the graph's 's spatial domain.
4 EXPERIMENTS
We execute series of experiments to illustrate that continuous GNNs formulated within the FROND framework using DJ outperform their traditional counterparts based integer-order derivatives. Importantly, our primary aim is not to achieve state-of-the-an results, but rather to demonstrate the additional effectiveness of the FROND framework when applied existing integer-order continuous GNNs. In the main paper, we detail the impressive results achieved by F-GRAND, particularly empha- sizing 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 Theorem 2, by constructing deeper GNNs with non- non-integer presented in the main paper, the basic predictor solver is used instead of other solvers when 3 < 1.
1. To maintain consistency the experiments
The Appendix D section provides additional details covering various aspects such as experimental settings, described in Appendices D.1 to D.3, the computational complexity F-GRAND in Appendix D.6, and analysis of F-GRAND's s robustness against adversarial attacks in Appendix D.9. Furthermore, results related to other FROND-based continuous GNNs are extensively presented in the Appendix E. In the main paper, we utilize the basic predictor, as delineated in (17), while the exploration of its variants is reserved for the Appendix D.5. Additional insights into the optimal fractional- Section Appendix D.11.
More Experiments In the Appendix:
order 3 and fractality in graph datasets are explored in
7

---

Published as a conference paper at ICLR 2024
4.1 NODE CLASSIFICATION OF F-GRAND
Datasets and splitting. We utilize datasets with varied topologies, including citation networks (Cora (McCallum et al., 2004), Citeseer (Sen et al., 2008), Pubmed (Namata et al., 2012)), tree-structured datasets (Disease and Airport (Chami al., 2019)), coauthor and co-purchasing graphs (CoauthorCS (Shchur et al., 2018), Computer and Photo (McAuley et al., 2015)), and the ogbn-arxiv dataset (Hu et al., 2020). We follow the same data splitting and pre-processing in (Chami et al., 2019) for Disease and Airport datasets. Consistent with experiment settings in GRAND (Chamberlain et al., 2021c), we use random splits for the largest connected component of each other dataset. We also incorporate the large-scale Ogbn-Products dataset (Hu et al., 2021) to demonstrate the scalability of the FROND framework, with the results displayed in Table 7. Methods. For a comprehensive performance comparison, we select several prominent GNN models as baselines, including GCN (Kipf& & Welling, 2017), and GAT (Veličković et al., 2018). Given the inclusion of tree-structured datasets, also incorporate well-suited baselines: HGCN(Chami et al., 2019) and GIL (Zhu et al., 2020b). To highlight the benefits of memorized dynamics in FROND, include GRAND (Chamberlain et al., 2021c) line with (Chamberlain et al., 2021c), we examine two F-GRAND variants: F-GRAND-nl (7) and F-GRAND-I (9). Graph rewiring not explored in this study. Where available, results from the paper (Chamberlain et al., 2021c) are used.
special case of F-GRAND with / 1. In
Table I: Node classification results(%) for random train-val val-test splits. The best and the second-best results are highlighted in red and blue, respectively.
Method
Com
Citesee
Pubmed
CoanthorCS
Computer
Photo
CoauthrPPyy
oghn-arxiv
Airport
Disense
GCN 81.5113
7199199
77.812.9
91.1.10.5
826112
91.2.1122
92.81.00
02.2-033
$161066
69.8..0.5
GAT 81841.
7142219
78.742.3
90.5:0.6
780=1990
85.7#20.3
92.22090
737=00.1
62..4
20.4=0.5
HGCN
78.711.0
65.81220
$06:18
88.211.4
90.811.5
59.620.4
85.420.7
89.111.1
GIL
82,111.1
71111:2
77.110.6
89411.
89.611.
91.11177
90.8.00.5
GRAND=1 GRAND-al
83
7344005
788++17
92.9+014
837=1.2
92,3+0.9
93.5+099
719-02
80.5+966
74.5+13.4
8.
70.91:0
77.5.1..8
92440.3
824422.1
92.410.88
91,411.3
71.2:002
90.9166
81.0.16.7
F-GRANDAI
848811
74.0=115
79.4=1.5
9300:003
$4,4=155
92.8=066
94,5+0.4
72.600.1
98.1=0022
92.4=3.9
for I
0.9
0.9
0.9
0.7
0.98
0.9
0.6
0.5
0.6
F. F-GRAND
83.211.1
74.7.19
79.2.1077 0.4
92.91044
$4.1.0.9
93.1.00.9
93,910.55
F1.4103
96.110.7 0.1
85 85.512.55
3 for F-GRAND-al
0.9
0.9
0.6
0.85
0.8
0.4
0.7
0.7
Table 2: Graph classification results.
POL
Feature
word2xec
BERT
Prodie
Accuraty Pancentage
80-
GraphSagge
67.6010068
10:66
St. 21.22.4881
02.0000:008
96,58.00 2...222
7828:00.2
95.96100755
74.03.055
38.990.78
82.7110119
91.850025
96.57.0.54
96.011.0455
Datacets:
GRANTS
7783-0337
16.7+1.13
N5.07-0.74
96.11-0.26
97,04-0.55
96:77+0.34
Coro
G
79.492043
38.6910.37
89.29..0.3
96.40:0.19
97.40:2..0
97.53.0.14
40-
Aiport
Table 3: Node classification accuracy of F- GRAND-I under different value of B when time T 8.
Hodals
GEN
0.1
0.5
0.7
09
Cura
74.8010422 77,01098
79.6010.91
81.56110300
82.68
82.37.0.39
Layers
Airoot 97.09+0077
|
95 05.80-2.05 91.6666.34
84.36+804
78.73-6.33
78.88+9.67
Figure 2: oversmoothing mitigation.
The results for graph node classification are summarized in Table 1, which also report the optimal 3 obtained via hyperparameter tuning. Consistent with our expectations, F-GRAND surpasses GRAND across nearly all datasets, given that GRAND represents special case of FROND 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- I outperforms both GRAND and GIL by approximately 7% on the Airport dataset. Interestingly, our experiments indicate a smaller 3 (signifying greater dynamic memory) is preferable for such innctal-trruttred datasets, aligning with previous studies on FDEs in biological and chemical systems (Nigmatullin, 1986; Ionescu et al., 2017). Further discussion on 3 and its relation to the fractal dimension of graph datasets can be found in Section 4.4 and Appendix D.11.
Performance.
with 3 =
8

---

Published as a conference paper at ICLR 2024
4.2 GRAPH CLASSIFICATION OF F-GRAND
We employ the Neke-NewsNet datasets (Dou et al., 2021), constructed from Politifact and Gossipcop fact-checking data. More details be found in the Appendix D.2. This dataset features three types of node features: 768-dimensional BERT features, and 300-dimensional spaCy features, both extracted using pre-trained models, and (0-dimensional profile features from Twitter accounts. The graphs in the dataset exhibit a hierarchical tree structure. From Table 2, we observe that F-GRAND consistently outperforms GRAND with notable edge on the POL. dataset.
4.3 OVERSMOOTHING 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 (17). This allows us to generate architectures of varying depths. In this context, utilize the fixed data splitting as described in (Chami et al., 2019). As illustrated in Fig. 2, optimal performance the Cora dataset is attained with network depth of 64 layers. When compared to GRAND-I, F-GRAND-I 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 Theorem predicts a slow algebraic convergence. In contrast, GRAND exhibits faster rate of performance degradation particularly on the Airport dataset. Further details on oversmoothing mitigation are in Appendix D.7.
4.4 ABLATION STUDY: SELECTION OF
In Table we investigate the influence of 3 across various graph datasets. Notably, for the Cora dataset, larger 3 is optimal, whereas, for tree-structured data, smaller 3 is preferable. This suggests that the quantity of memorized dynamics should tailored to the dataset's topology, and a default setting of memoryless graph diffusion with 8 1 may not be optimal. More comprehensive details concerning the variations in 33 can be found the appendix, specifically in Table 15.
4.5
MORE INTEGER-ORDEE CONTINUOUS GNNS IN FROND FRAMEWORK
Our FROND framework can be seamlessly applied to various other integer-order continuous GNNs, as elaborated in Appendix F. Specifically, here we outline the node classification results of FROND) based the CDE model in Table 4. It is evident from the results that F-CDE enhances the perfor- mance of the CDE model across almost all large heterophilic datasets. The optimal / is determined through hyperparameter tuning. When 3 1, F-CDE seamlessly reverts CDE, and the results from the original paper are reported. Additionally, we conduct comprehensive experiments detailed in Appendix E. The results for F-GRAND++, F-GREAD, and F-GraphCON are available in Ta- ble 19, Table 23, and Table 25, 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.
Table 4: Node classification accuracy(%) of large heterophilic datasets
Model
Roman-empire
Wiki-cooe
Minesweeper
Questions
Workers
Amazon-ratings
CDE
91. 00.648 93.06:10.55
97.99-10.38
95.50-t5.23
17:1999999
80.70 1.04
47.63=0.43
F-CDE
98.73-10.68 01.6
9..04
75.17
82.6
49.01=0.56
8 for F-DDE
0.9
0.6
1.0
0.4
0.1
5 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 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.

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS AND DISCLOSURE OF FUNDING
This research is supported by the Singapore Ministry of Education Academic Research Fund Tier 2 grant MOE-T2HP20220-00
and the National Research Foundation, Singapore and Infocomm Media Development Authority under its Future Communications Research and Development Pro- gramme. The computational work for this article was partially performed resources of the National Supercomputing Centre, Singapore httpss//www.nscc..gg 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 OpenAl (2022).
REFERENCES
Ricardo Almeida, Nuno RO Bastos, and M Teresa T Monteiro. Modeling some real phenomena by fractional differential equations. Mathematical Methods in the Applied Sciences, 39(16): 4846-4855, 2016.
Uri Alon and Eran Yahav. On the bottleneck of graph neural networks and its practical implications. In Proc. Int. Conf. Learn. Representations, 2021.
Harbir Antil, Ratna Khatri, Rainald Löhner, and Deepanshu Verma. Fractional deep neural network via constrained optimization. Mach. Learn.: Sci. Technol., 2(1):015003, 2020.
Haitham Ashoor, Xiaowen Chen, Wojciech Rosikiewicz, Jiahui Wang, Albert Cheng, Ping Wang, Yijun Ruan, and Sheng Li. Graph embedding and unsupervised learning predict genomic sub- compartments from hic chromatin interaction data. Nat. Commun., 11, 2020.
Kendall Atkinson, Weimin Han, and David E Stewart. Numerical solution of ordinary differential equations. John Wiley & Sons, 2011.
Pedro HC Avelar. Anderson Tavares, Mareo Gori, and Luis C Lamb. Discrete and continuous deep residual learning over graphs. arXiv preprint arXiv: 1919599 2019.
Ronald L Bagley and PJ Torvik. A theoretical basis for the application of fractional calculus to viscoelasticity. Rheology, (7(3):001-210 1983.
Deyu Bo, Xiao Wang, Chuan Shi, and Huawei Shen. Beyond low-frequency information in graph con- volutional networks. In Proceedings of the AAAI Conference Artificial Intelligence, volume 35, pp. 301333950995 2021.
Cristian Bodnar, Francesco Di Giovanni, Benjamin Paul Chamberlain, Pietro Liò, and Michael M. Bronstein. Neural sheaf diffusion: A topological perspective on heterophily and oversmoothing in GNNs. In Advances Neural Inf. Process. Syst., 2022.
Dirk Brockmann, Lars Hufnagel, and Theo Geisel. The scaling laws of human travel. Nature, 439 (7075):462-465, 2006.
Ben Chamberlain, James Rowbottom, Maria I Gorinova, Michael Bronstein, Stefan Webb, and Emanuele Rossi. Grand: Graph neural diffusion. In Proc. Int. Conf. Mach. Learn., pp. 1407-1418,
2021a.
Benjamin Chamberlain, James Rowbottom, Davide Eynard, Francesco Di Giovanni, Xiaowen Dong, and Michael Bronstein. Beltrami flow and neural diffusion on graphs. In Advances Neural Inf. Process. Syst., pp. 1594-1609, 2021b.
Benjamin Paul Chamberlain, James Rowbottom, Maria Goronova, Stefan Webb, Emanuele Rossi, and Michael M Bronstein. Grand: Graph neural diffusion. In Proc. Int. Conf. Mach. Learn., 2021c
Ines Chami, Zhitao Ying, Christopher Ré, and Jure Leskovec. Hyperbolic graph convolutional neural networks. In Advances Neural Inf. Process. Syst. 2019.
Jinyin Chen, Yangyang Wu, Xuanheng Xu, Yixian Chen, Haibin Zheng, and Qi Xuan. Fast gradient attack network embedding. ArXiv, 2018a.
10

---

Published as a conference paper at ICLR 2024
Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li. Simple and deep graph convolutional networks. In Proc. Int. Conf. Mach. Learn., pp. 1725-1735, 2020.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In Advances Neural Inf. Process. Syst., 2018b.
Eli Chien, Jianhao Peng, Pan Li, and Olgica Milenkovic. Adaptive universal generalized pagerank graph neural network. arXiv preprint arXiv:2006. 07988, 2020.
Jeongwhan Choi, Scoyoung Hong, Noseong Park, and Sung-Bat Bac Cho. Gread: Graph neural reaction- diffusion networks. In Proc. Int. Conf. Mach. Learn., 2023.
Fan RK Chung. Spectral graph theory, volume 92, American Mathematical Soc., 1997.
Alan M. Cohen. Inversion Formulae and Practical Results, pp. 23-44. Springer US, Boston, MA, 2007.
Bernard D Coleman and Walter Noll. Foundations of linear viscoelasticity. Rev. Modern Phys., 33 (2):239, 1961.
Weihua Deng. Short memory principle and a predictor-correcto approach for fractional differential equations. J. Comput. Appl. Math., 206(1):174-188, 2007.
Francesco Di Giovanni, James Rowbottom, Benjamin P Chamberlain, Thomas Markovich, and Michael M Bronstein. Graph neural networks as gradient flows. arXiv preprint arXiv:2206.1099) 2022.
Francesco Di Giovanni, James Rowbottom, Benjamin Paul Chamberlain, Thomas Markovich, and Michael M Bronstein. Understanding convolution on graphs via energies. Tran. Mach. Learn. Res., 2023.
Fernando Diaz-Diaz and Ernesto Estrada. Time and space generalized diffusion equation on graph/networks. Chaos, Solitons and Fractals, 156:1 1,,,,,
Kai Diethelm. The analysis of fractional differential equations: an application-oriented exposition using differential operators of Caputo type, volume 2004. Springer, 2010.
Kai Diethelm and Neville J Ford. Analysis of fractional differential equations. J. Math. Anal. Appl., 265(2):229-244 2002.
Kai Diethelm, Neville Ford, and Alan Freed. Detailed error analysis for fractional adams method. Numer. Algorithms, 36:31-52, 2004.
Yingtong Dou, Kai Shu, Congying Xia, Philip S. Yu, and Lichao Sun. User preference-aware fake news detection. In Proc. Int. ACM SIGIR Conf. Res. Develop. Inform. Retrieval, 2021.
Jian Du, Shanghang Zhang, Guanhang Wu, José M. F. Moura, and Soummya Kar. Topology adaptive graph convolutional networks. ArXiv, hts//770..10700 2017.
Lun Du, Xiaozhou Shi, Qiang Fu, Xiaojun Ma, Hengyu Liu, Shi Han, and Dongmei Zhang. Gbk- gnn: Gated bi-kernel graph neural networks for modeling both homophily and heterophily. In Proceedings of the ACM Web Conference 2022, pp. 1550-1558, 2022.
Emilien Dupont, Arnaud Doucet, and Yee Whye Teh. Augmented neural odes. In Advances Neural Inf. Process. Syst., pp. 1-11, 2019.
Jiarui Feng, Yixin Chen, Fuhai Li, Anindya Sarkar, and Muhan Zhang. How powerful are k-hop message passing graph neural networks. Advances Neural Information Processing Systems, 35: 4776-4790, 2022.
Guang-hua Gao and Zhi-zhong Sun. A compact finite difference scheme for the fractional sub- diffusion equations. Journal of Computational Physics, 230(3):586-595 2011.
Hongyang Gao and Shuiwang Ji. Graph u-nets. Proc. Int Conf. Mach. Learn., pp. 2083-2092, 2019.
11

---

Published as a conference paper at ICLR 2024
José Francisco Gómez-Aguil Margarita Miranda-Hemnnndzz, MG López-Lópezz Victor Manuel and Dumitru Baleanu. Modeling and simulation the fractional space-time diffusion equation. Commun. Nonlinear Sci. Numer. Simul., 30(1-3):115-127 2016.
Rudolf Gorenñlo and Francesco Mainardi. Fractional diffusion processes: probability distributions and continuous time random walk. In Process. Long-Range Correlations: Theory Appl., pp. 148-166. Springer, 2003.
Rudolf Gorenflo, Francesco Mainardi, Daniele Moretti, and Paolo Paradisi. Time fractional diffusion:
discrete random walk approach. Nonlinear Dynamics, 19:129-143, 2002.
Alessio Gravina, Davide Bacciu, and Claudio Gallicchio. Anti-symmetric dgn: A stable architecture for deep graph networks. In Proc. Int. Conf. Learn. Representations, 2022.
Ling Guo, Hao Wu, Xiaochen Yu, and Tao Zhou. Monte carlo fpinns: Deep learning method for forward and inverse problems involving high dimensional fractional partial differential equations. Comput. Methods Appl. Mechanics Eng., 400: 15523, 2022.
Kyle B Gustafson, Basil S Bayati, and Philip A Eckhoff. Fractional diffusion emulates human mobility network during a simulated disease outbreak. Frontiers Ecology Evol., 5:35, 2017.
Benjamin Gutteridge, Xiaowen Dong, Michael M Bronstein, and Francesco Di Giovanni. Drew: Dynamically rewired message passing with delay. Proc. Int. Conf. Much. Learn., pp. 12252- 12267, 2023.
Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse Problems, 34 (1):1-23, December 2017.
William L. Hamilton, Rex Ying, and Jure Leskovec. Inductive representation learning large graphs. In Advances Neural Inf. Process. Syst., 2017.
Andi Han, Dai Shi, Lequan Lin, and Junbin Gao. From continuous dynamics to graph neural networks: Neural diffusion and beyond. arXiv preprint arXiv:2310.1012/, 2023.
Philip Hartman. Ordinary differential equations. SIAM, 2002.
Kaiming He. Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proc. Conf. Comput. Vision Pattern Recognition, 2016.
Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, New York, 2012.
Weihua Hu, Matthias Fey, Marinka Zitnik, Yuxiao Dong. Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec. Open graph benchmark: Datasets for machine learning on graphs. arXiv:2005.0066 2020.
Weihua Hu, Matthias Fey, Marinka Zitnik, Yuxiao Dong, Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec. Open graph benchmark: Datasets for machine learning graphs, 2021.
Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Weinberger. Densely connected convolutional networks. In Proc. Conf. Comput. Vision Puttern Recognition, 2017..
Hussain Hussain, Meng Cao, Sandipan Sikdar, Denis Helic, Elisabeth Lex, Markus Strohmaier, and Roman Kern. Adversarial inter-group link injection degrades the faimess of graph neural networks. arXiv preprint arXiv:2209.05%) 2022.
Tonescu, A Lopes, Dana Copot, JA Tenreiro Machado, and Jason HT Bates. The role of fractional calculus in modeling biological phenomena: A review. Commun. Nonlinear Sci. Numer. Simul., 51:141-159, 2017.
Rana Javadi, Hamid Mesgarani, Omid Nikan, and Zakieh Avazzadeh. Solving fractional order differential equations by using fractional radial basis function neural network. Symmetry, 15(6): 1275, 2023.
12

---

Published as a conference paper at ICLR 2024
Bangti Jin, Buyang Li, and Zhi Zhou. Correction of high-order bdf convolution quadrature for fractional evolution equations. SIAM J. Sci. Comput., 39())33319933333111 2017.
Wei Jin, Yao Ma, Xiaorui Liu, Xianfeng Tang, Suhang Wang, and Jiliang Tang. Graph structure learning for robust graph neural networks. In Proc. Int. Conf. Knowl. Discovery Data Mining, pp. 66-74, 2020.
Qiyu Kang, Yang Song, Qinxu Ding, and Wee Peng Tay. Stable neural ODE with Lyapunov-stable equilibrium points for defending against adversarial attacks. In Advances Neural Inf. Process. Syst., 2021.
Qiyu Kang, Kai Zhao, Yang Song, Sijie Wang, and Wee Peng Tay. Node embedding from neural Hamiltonian orbits in graph neural networks. Proc. International Conference on Machine Learning, pp. 15786-15808, 2023.
Qiyu Kang, Kai Zhao, Yang Song, Yihang Xie, Yanan Zhao, Sijie Wang, Rui She, and Wee Peng Tay. Coupling graph neural networks with fractional order continuous dynamics: A robustness study. In Proc. AAAI Conference on Artificial Intelligence, Vancouver, Canada, 2024.
IS Kim, K-I Goh, G Salvi, E Oh, B Kahng, and D Kim. Fractality in complex networks: Critical and supercritical skeletons. Physical Review E, 75(1):016110, 2007.
Diederik P Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint 1:XXv::4422.99000 2014.
Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In Proc. Int. Conf. Learn. Representations, 2017.
Granino Arthur Korn and Theresa M Korn. Mathematical handbook for scientists and engineers: definitions, theorems, and formulas for reference and review. Courier Corporation, 2000.
Diego Krapf. Mechanisms underlying anomalous diffusion in the plasma membrane. Current Topics Membranes, 75:167-207 2015.
Guohao Li, Matthias Muller, Ali Thabet, and Bernard Ghanem. Deepgens: Can gens go deep as cnns? Proc. Int. Conf. Learn. Representations, pp. 9267-9276, 2019.
Guohao Li, Chenxin Xiong, Ali Thabet, and Bernard Ghanem. Deepergen: All you need to train deeper gens, arXiv preprint wwXiv:2006.07739, 2020a.
Xiang Li, Renyu Zhu, Yao Cheng, Caihua Shan, Sigiang Luo, Dongsheng Li, and Weining Qian. Finding global homophily in graph neural networks when meeting heterophily. In International Conference on Machine Learning, pp. 13242-13256. PMLR, 2022.
Yaxin Li, Wei Jin, Han Xu, and Jiliang Tang. Deeprobust: A pytorch library for adversarial attacks and defenses. arXiv preprint 200:::20000014499 2020h.
Derek Lim, Felix Hohne, Xiuyu Li, Sijia Linda Huang. Vaishnavi Gupta, Omkar Bhalerao, and Ser Nam Lim. Large scale learning on non non-homophilous graphs: New benchmarks and strong simple methods. Advances in Neural Information Processing Systems, 34:20887-20902. 2021.
Zijian Liu, Yaning Wang, Yang Luo, and Chunbo Luo. A regularized graph neural network based on approximate fractional order gradients. Mathematics, 10(8):1320, 2022.
Sitao Luan, Chenqing Hua, Qincheng Lu, Jiaqi Zhu, Mingde Zhao, Shuyuan Zhang, Xiao-Wen Chang, and Doina Precup. Revisiting heterophily for graph neural networks. Advances in neural information processing systems, 35:1 13:1332-1375. 2022.
Chunwan Lv and Chuanju Xu. Error analysis of a high order method for time-fractional diffusion equations. SIAM J. Sci. Comput., A87744 2016.
J Tenreiro Machado, Virginia Kiryakova, and Francesco Mainardi. Recent history of fractional calculus. Communications in nonlinaaar science and numerical simulation, (160):1140-1155 2011.
13

---

Published as a conference paper at ICLR 2024
Benoit B Mandelbrot and Benoit B Mandelbrot. The fractal geometry of nature, volume 1. WH freeman New York, 1982.
Sohir Maskey, Raffaele Paolino, Aras Bacho, and Gitta Kutyniok. A fractional graph laplacian approach oversmoothing. arXiv preprint arXiv:2305. 2023.
Barry R Masters. Fractal analysis of the vascular tree in the human retina. Annu. Rev. Biomed. Eng., 6:427-452, 2004.
Julian McAuley, Christopher Targett, Qinfeng Shi, and Anton van den Hengel. Image-based baage-baedd rec- ommendations styles and substitutes. In Proc. Int. ACM SIGIR Conf. Res. Develop. Inform. Retrieval, pp. 43-52, 2015.
Andrew McCallum, Kamal Nigam, Jason D. M. Rennie, and Kristie Seymore. Automating the construction of internet portals with machine learning. Inf. Retrieval, 3:127-163, 2004.
Federico Monti, Davide Boscaini, Jonathan Masci, Emanuele Rodola, Jan Svoboda, and Michael M Bronstein. Geometrie deep learning graphs and manifolds using mixture model cnns. In Proc. Conf. Comput. Vision Pattern Recognition, pp. 5115-5124, 2017.
Galileo Mark Namata, Ben London, Lise Getoor, and Bert Huang. Query query-driven active surveying for collective classification. In Workshop Mining Learn. Graphs, 2012.
Ravil' Rashidovich Nigmatullin. Fractional integral and its physical interpretation. Theoretical and Mathematical Physics, 90(3):242-255) 1992.
RR Nigmatullin. The realization of the generalized transfer equation in medium with fractal geometry. Physica status solidi (b), 133(1):425 430, 1986.
Bachmann- Order Notations. Big notation, 2023. URL https://enn wttp://di...ii.iiia........ oxt/wwiki/Biq___nottat
Accessed: Sep 1, 2023.
Kenta Oono and Taiji Suzuki. Graph neural networks exponentially lose expressive power for node classification. In Proc. Int. Conf. Learn. Representations, 2020.
OpenAL. Chatgpt-4, 2022. Available at: https: http://wwwpp.....ce openai. .com com (Accessed: 10 April 2024).
Guofei Pang, Lu Lu, and George Em Kariaadkis. fpinns: Fractional physics-informed neural
networks. SIAM Sci. Comput., 4144)A2603-A2266, 2019.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. In Advances Neural Inf. Process. 2017,2217
Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang. Geom-gen: Geometric graph convolutional networks. arXiv preprint arXiw:2002.05287, 2020.
Igor Podlubny. Fractional-order systems and fractional-order controllers. Institute of Experimental Physics, Slovak Academy of Sciences, Kosice, 12(3):1-18, 1994.
Igor Podlubny. Fractional Differential Equations. Academic Press, 1999.
Michael Poli, Stefano Massaroli, Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Graph neural ordinary differential equations. arXiv preprint 171.52191..07532, 2019.
Alessio Quaglino, Marco Gallieri. Jonathan Masci, and Jan Koutník. Snode: Spectral discretization of neural odes for system identification. In Proc. Int. Conf. Learn. Representations, 2019.
Ahmed Gomaa Radwan, Ahmed Elwakil, and Ahmed M Soliman. Fractional-order sinusoidal oscillators: design procedure and practical examples. IEEE Tran. Circuits and Syst., 55(7): 2008.
Konstantin Rusch, Ben Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael Bronstein. Graph-coupled oscillator networks. In Proc. Int. Conf. Mach. Learn., 2022.
14

---

Published as a conference paper at ICLR 2024
Enrico Scalas, Rudolf Gorenflo, and Francesco Mainardi. Fractional calculus and continuous-time finance. Physica A: Statistical Mechanics and its Applications, 284(1-4):376-33 2000.
Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Eliassi-Rad. Collective classification in network data. Al Magazine, 29(3):93, Sep. 2008.
Oleksandr Shchur, Maximilian Mumme, Alcksandar Bojchevski, and Stephan Günnemann. Pitfalls of graph neural network evaluation. Relational Representation Learning Workshop, Advances Neural Inf. Process. Syst., 2018.
Rui She, Qiyu Kang, Sijie Wang. Wee Peng Tay, Yong Liang Guan, Diego Navarro Navarro, and Andreas Hartmannsgruber, Image aatch-matching with graph-based learning in street scenes. IEEE Trans. Image Process., 32:3465-3480, 2023a.
Rui She, Qiyu Kang, Sijie Wang, Yuán-Rui Yáng, Kai Zhao, Yang Song, and Wee Peng Tay. Robustmat: Neural diffusion for street landmark patch matching under challenging environments. IEEE Trans. Image Process., 2023b.
Rui She, Qiyu Kang, Sijie Wang, Wee Peng Tay, Kai Zhao, Yang Song, Tianyu Geng, Yi Xu, Diego Navarro Navarro, and Andreas Hartmannsgruber. PointDifformer: Robust point cloud registration with neural diffusion and transformer. IEEE Transactions on Geoscience and Remote Sensing, 62:1 - 15, 2024a. doi: 10.
Rui She, Sijie Wang, Qiyu Kang. Kai Zhao, Yang Song, Wee Peng Tay, Tianyu Geng, and Xingchao Jian. PosDifINet: Positional neural diffusion for point cloud registration in large field of view with perturbations. In Proc. AAAI Conference on Artificial Intelligence, Vancouver, Canada, 2024b.
Chaoming Song, Shlomo Havlin, and Hernan A Makse. Self-similarity of complex networks. Nature, 433(7024):392-395, 2005.
Chaoming Song, Lazaros K Gallos, Shlomo Havlin, and Hernán A Makse. How calculate the fractal dimension of a complex network: the box covering algorithm. J. Stat. Mech. Theory Exp., 2007(03):P03006, 2007.
Yang Song, Qiyu Kang, Sijie Wang. Kai Zhao, and Wee Peng Tay. On the robustness of graph neural diffusion to topology perturbations. In Advances Neural Inf. Process. Syst., 2022.
Didier Sornette. Critical phenomena in natural sciences: chaos, fractals, selforganization and disorder: concepts and tools. Springer Science & Business Media, 2006.
Zhi-zhong Sun and Xiaonan Wu. fully discrete difference scheme for diffusion-wave system. Applied Numerical Mathematics, 56(2):193-209, 2006.
Vasily E Tarasov. Fractional dynamics: applications of fractional calculus to dynamics of particles, fields and media. Springer Science & Business Media, 2011.
Matthew Thorpe, Hedi Xia, Tan Nguyen, Thomas Strohmer, Andrea Bertozzi, Stanley Osher, and Bao Wang. Grand++: Graph neural diffusion with a source term. In Proc. Int. Conf. Learn. Representations, 2022.
Wen Yi Tian, Han Zhou, and Weihua Deng. A class of second order difference approximations for solving space fractional diffusion equations. Math. Comput., 84(294): 1 1703-1277 2015.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention all you need. Advances in neural information processing systems, 30, 2017
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks. In Proc. Int. Conf. Learn. Representations, pp. 1-12, 2018.
Jihong Wang, Minnan Luo, Fnu Suya, Jundong Li, Zijiang Yang, and Qinghua Zheng. Scalable attack on graph data by injecting vicious nodes. Data Mining Knowl. Discovery, pp. - 27, 2020.
Shupeng Wang, Hui Zhang, and Xiaoyun Jiang. Fractional iiysics-ifformed neural networks for time-fractional phase field models. Nonlinear Dyn., 1103))2715-2739, 2022a
15

---

Published as a conference paper at ICLR 2024
Sijie Wang, Qiyu Kang. Rui She, Wee Peng Tay, Andreas Hartmannsgruber, and Diego Navarro Navarro. RobustLoc: Robust camera pose regression in challenging driving environments. In Proc. AAAI Conference Artificial Intelligence, Feb. 2023.
Yuelin Wang, Kai Yi, Xinliang Liu, Yu Guang Wang, and Shi Jin. Acmp: Allen-cahn message passing with attractive and repulsive forces for graph neural networks. In Proc. Int. Conf. Learn. Representations, 2022b.
Marein Waniek, Tomasz Michalak, Michael J. Wooldridge, and Talal Rahwan. Hiding individuals and communities in social network. Nature Human Behaviour, 2(1):139-144: 2018.
Ee Weinan. A proposal on machine learning via dynamical systems. Commun. Math. Statist., 1(5): 1-11, 2017.
Wikipedia. Hardy-littlewood tauberian theorem, 2023. URL hwtph:tpp////eee wi www.iiiiiiee pedia a wedi org/
wiki/Ha $22880 wew.ittllleee ottpe///.....tttoooooooooooaaaae
Accessed: Sep 1, 2023.
Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S. Yu. A comprehensive survey on graph neural networks. IEEE Trans. Neural Netw. Learn. Syst.. 32(1): 4-24. 2021.
Louis-Pasce Xhonneux, Meng Qu, and Jian Tang. Continuous graph neural networks. Proc. Iut. Conf. Mach. Learn., pp. 10432-10441, 2020.
Keyulu Xu, Chengtao Li, Yonglong Tian, Tomohiro Sonobe, Ken-ichi Kawarabayashi, and Stefanic Jegelka. Representation learning on graphs with jumping knowledge networks. In Proc. Int. Conf. 2018.
Mach. Learn.. pp.
Hanshu Yan, Jiawei Du, Vincent YF Tan, and Jiashi Feng. On robustness of neural ordinary differential equations. In Advances Neural Inf. Process. Syst., pp. 1-13, 2018.
Yujum Yan, Milad Hashemi, Kevin Swersky, Yaoqing Yang, and Danai Koutra. Two sides of the same coin: Heterophily and oversmoothing graph convolutional neural networks. In 2022 IEEE International Conference on Data Mining (ICDM), pp. 1287-1292. IEEE, 2022.
Xiang Yue, Zhen Wang. Jingong Huang. Srinivasan Parthasarathy, Soheil Moosavinasab, Yungui Huang, Simon Lin, Wen Zhang, Ping Zhang, and Huan Sun. Graph embedding on biomedical networks: methods, applications and evaluations. Bioinformatics, (2444):121-1255 2019.
Santos B Yuste and Luis Acedo. An explicit finite difference method and new von neumann-type stability analysis for fractional diffusion equations. SIAM Journal on Numerical Analysis, 42(5): 1862-1874, 2005.
Hanqing Zeng, Hongkuan Zhou, Ajitesh Srivastava, Rajgopal Kannan, and Viktor Prasanna. Graph- saint: Graph sampling based inductive learning method, 2020.
Ziwei Zhang, Peng Cui, and Wenwu Zhu. Deep learning graphs: A survey, IEEE Trans. Knowl. Data Eng., 34(1):249-270, Jan 2022.
Kai Zhao, Qiyu Kang, Yang Song, Rui She, Sijie Wang, and Wee Peng Tay. Graph neural convection- diffusion with heterophily. In Proc. Inter. Joint Conf. Artificial Intell., Macao, China, 2023a.
Kai Zhao, Qiyu Kang, Yang Song, Rui She, Sijie Wang, and Wee Peng Tay. Adversarial robustness in graph neural networks: A Hamiltonian energy conservation approach. In Advances in Neural Information Processing Systems, New Orleans, USA, 2023b.
Qinkai Zheng, Yixiao Fei, Yanhao Li, Qingmin Liu, Minhao Hu, and Qibo Sun. Kdd cup 2020 ml track 2 adversarial attacks and defense academic graph 1st place solution, 2022. URL https:/ / thub. wtt.ccoam///////...e Sttnnttnnnisnnn wam.isslss///.mm aaas//DDDD /DD_//DD_CCCC CUP 2000p22000_00010002210022 21F220MMTracc2___P_22_ SPETT. Accessed:
May 1, 2022.
Jiong Zhu, Yujun Yan, Lingxiao Zhao, Mark Heimann, Leman Akoglu, and Danai Koutra. Beyond homophily in graph neural networks: Current limitations and effective designs. Advances in neural information processing systems, 33:7793-7804, 2020a.
16

---

Published as a conference paper at ICLR 2024
Jiong Zhu, Ryan A Rossi, Anup Rao, Tung Mai, Nedim Lipka, Nesreen K Ahmed, and Danai Koutra. Graph neural networks with heterophily. In Proceedings of the AAAI conference on artificial intelligence, volume 35, pp. 11168-1117,,,,11111111111111-111 2021.
Shichao Zhu, Shirui Pan, Chuan Zhou, Jia Wu, Yanan Cao, and Bin Wang. Graph geometry interaction learning. In Advances Neural Inf. Process. Syst., 2020b.
Juntang Zhuang, Nicha Dvorekk Xiaoxiao Li, and James S Duncan. Ordinary differential equations on graph networks. 2019.
Xu Zou, Qinkai Zheng, Yuxiao Dong, Xinyu Guan, Evgeny Kharlamov, Jialiang Lu, and Jie Tang. Tdgia: Effective injection attacks graph neural networks. In Proc. Int. Conf. Knowl. Discovery Data Mining, pp. 2461-2471,,0222
Daniel Zügner and Stephan Günnemann. Adversarial attacks on graph neural networks via meta learning. Proc. Int. Conf. Learn. Representations, 2019.
Daniel Zügner, Amir Akbarnejad, and Stephan Günnemann. Adversarial attacks on neural networks for graph data, In Proc. Int. Conf. Knowl. Discovery Data Mining, 2018.
17

---

Published as a conference paper at ICLR 2024
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:
1. We discuss related work in Appendix A.
2. We offer concise review of fractional calculus in Appendix B.
3. We include more solver details and variants in Appendix C.
4. We present dataset statistics, experimental settings, and additional experimental results in Appendix D.
5. We introduce more dynamics within the FROND framework in Appendix E.
6. We provide proofs for all theoretical assertions made in the main paper in Appendix F.
7. We discuss the limitations of our work and its broader impact in the final section of this supplementary material.
A RELATED WORK
Fractional Calculus and Its Applications
The field of fractional calculus has seen notable surge interest recently due to its wide-ranging applications across various domains. These include, but are not limited to, numerical analysis (Yuste & Acedo, 2005), viscoelastic materials (Coleman & Noll, 1961), population growth models (Almeida et al., 2016), control theory (Podlubny, 1994), signal processing (Machado et al., 2011), financial mathematics (Scalas et al., 2000), and particularly in the representation of porous and fractal phenomena (Nigmatullin. 1986; Mandelbrot & Mandelbrot, 1982; lonescu et al., 2017). Within these contexts, FDEs have been developed as a powerful extension to the conventional integer- differential equations, offering resilient mathematical framework for system analysis (Diethelm & Ford, 2002). 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 (Krapf, 2015), to animal migration patterns (Brockmann et al., 2006), human mobility networks (Gustafson et al., 2017), and even biological phenomena pertinent to respiratory tissues and neuroscience (Ionescu et al., 2017). 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 (Diaz-Diaz & Estrada, 2022; Sornette, 2006).
Within the realm of deep learning, (Liu et al., 2022) 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 (Kingma & Ba, 2014) with respect to the weights. The essence of their work fundamentally differs from ours, which focuses the fractional-order evolution of node embeddings. not gradient optimization. A detailed examination of the study by (Liu et al., 2022) is pivotal as it adopts fractional derivatives instead of the standard first-order derivatives during the weight updating phase of GNN in the gradient descent. Specifically, attention is drawn to equation (16) in (Liu et al.. 2022), 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 (6) of our paper, the fractional derivative applied to the evolving node feature, representing an implementation of a fractional- feature updating process, thereby showcasing a clear distinction in the application of fractional derivatives.
Additionally. (Antil al., 2020) incorporates insights from fractional calculus and its L1 approx- imation of the fractional derivative to craft densely connected neural network. Their aim is to adeptly handle non-smooth data and counteract the vanishing gradient problem. While our research operates within 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 connection to non-Markoviar dynamic processes. Our framework paves the way for a new class GNNs, enabling a wide spectrum of learnable feature- updating processes influenced by memory effects.
From the perspective of physics- to crafting neural networks rooted in physical laws solve fractional PDEs. A pioneering work
physics-informed informed machine learning, another line of research is dedicated
18

---

Published as a conference paper at ICLR 2024
in this domain is the Fractional Physics Informed Neural Networks (fPINNs) (Pang al., 2019). Subsequent research, such as (Guo et al., 2022; Javadi et al., 2023; Wang et al., 2022a), 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 framework for deep learning has been initially explored by (Weinan, 2017). The seminal work of (Chen et al., 2018b) 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 's explainability (Weinan, 2017; Chamberlain et al., 2021c). Additionally, separate advancements in this domain have led to improvements in neural network performance (Dupont et al., 2019), robustness(Yan) (an et al., 2018; Kang et al., 2021), and gradient stability (Haber & Ruthotto, 2017: Gravina et al., 2022). In practical applications, neural ODEs are demonstrating superior performance (She et al., 2024a;b; 2023b; Wang et al., 2023; She et al., 2023a). In a similar vein, (Avelar et al., 2019) models continuous residual layers in GCN, leveraging neural ODE solvers to produce output. Further, the work of (Poli et al., 2019) proposes 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 (Zhuang et al., 2019) 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 (Rusch al., 2022) adopts the coupled oscillator model that preserves the graph's ' Dirichlet energy over time and mitigates the oversmoothing problem. In (Chamberlain et al., 2021a), the authors modeled information propagation as a diffusion process of a substance from regions of higher to lower concentration. The Beltrami diffusion model utilized in (Chamberlain et al., 202 bb; Song et al., 2022) to enhance rewiring and improve the robustness of the graph. The study by (Bodnar et al., 2022) 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 (Wang et al., 2022b) is inspired by particle reaction-diffusee interactions between particles. Concurrently, the graph CDE model (Zhao et al., 2023a) is crafted to handle heterophilic graphs and is inspired by the convection -iifttiionn process. GRAND++ (Thorpe al., 2022) leverages heat diffusion with sources to train models effectively with limited amount of labeled training data. Concurrently, GREAD (Choi et al., 2023) articulates a GNN approach, which is premised on craftiondliffusion equations. aiming negotiate heterophilic datasets effectively. In another development, the continuous GNN as an ODE (Maskey et al., 2023) 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 (Han et al., 2023) on continuous GNNs for a more thorough summarization.
processes, taking into account repulsive and attractive force
Our FROND extends the above integer-order continuous GNNs by incorporating the Caputo fracoffrctionn fractional derivative. The models mentioned can be reduced from our unifted mathematical frumework, with variations manifesting from the choice of the dynamic operator F(W, (t)) in (6) and as B equals 1 in the fractional derivative operator D---
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 (He et al., 2016), 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 (Gao & Ji, 2019) employs skip connections to enable efficient information propagation across layers. Similarly, the Jump Knowledge Network (Xu et al., 2018) implements aayer-aggregattin mechanism that amalgamates outputs from all preceding layers, strategy reminiscent of the dense connections found in DenseNet (Huang et al., 2017). Furthermore, the work (Chen et al., 2020) introduces GCNII, an extension of the standard GCN model that incorporates
19

---

Published as a conference paper at ICLR 2024
two simple techniques, initial residual and identity mapping, tackle the oversmoothing problem. Expanding on the idea of depth in GNNs, (Li et al., 2019; 2020a) propose DeepGCNs, an innovative architecture that employs residual/dense connections along with dilated convolutions. The work (Di Giovanni et al., 2023) suggests that gradient- message passing neural networks may be able to deal with heterophilic graphs provided that a residual connection is available. The paper (Gutteridge et al., 2023) proposes a spatial domain rewiring and focuses on long-range interactions. DRew in (Gutteridge al., 2023) does not adhere to any ODE evolutionary structure. Additionally, the skip connection in the vDRew from (Gutteridge et al., 2023) specifically links an k-th layer to the n-th layer. This design fundamentally different from our FDE approach.
By incorporating fractional calculus and memory effects into our framework, we not only offer 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.
B 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 more detailed overview of it and explain why it widely employed in applications. We have based our FROND framework 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 GRAND (Chamberlain et al., 2021c), GraphCON (Rusch et al., 2022), GRAND++ (Thorpe et al., 2022), GREAD (Choi et al., 2023) and CDE (Zhao al., 2023a).
the assumption that the solution to the fractional differential equation
B.1 CAPUTO FRACTIONAL DERIVATIVE AND ITS COMPATIBILITY OF WNTEGER-ORDER DERIVATIVE
In the main paper, our focus is predominantly on the order 3 E (0, 1] for the sake of simplification. The Caputo fractional derivative of a function f(t) over an interval [O, b]. of a general positive order B E (0, oo), is defined as follows:
1 D2 -9) (-00-))00-0-11100000
(18)
Here, [.8] is the smallest integer greater than or equal to 3, I'(-) denotes the gamma function, and (l)) denotes the [3] -order derivative of f(-). Within this definition, presumed that f[TA1l E L [0, 6). i.e., f[T@1]; is Lebesgue integrable, to ensure the well-defined nature of D ft) as per (18) (Diethelm, 2010). For vector-valued function, the Caputo fractional derivative is defined on a component- component-y-ccomoneen bomoopnetttt 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 general order B E (0, oo) is presented in Theorem 7.1 (Diethelm, 2010) as:
[ Σ
( &
8-k llk (0))
(19)
where we assume that the Laplace transform Lf exists on [so, xc) for some s0 E R. In contrast, for the integer-order derivative f[A] where is a positive integer, we also have the formulation (19), with the only difference being the range of B. Therefore, 3 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
20

---

Published as a conference paper at ICLR 2024
generalizes the traditional integer-order derivative since their Laplace transforms coincide when takes integer value. The inverse Laplace transform specifies the uniquely determined D f( when |3 is an integer (in the sense of almost everywhere (Cohen, 2007)).
Under specific reasonable conditions, (t) (18) is continuously differentiable. In this context, integration by parts can be utilized to demonstrate that
directly present this generalization follows. Suppose
(- ((T (t [8) (t-T)
[8] - + (' (1101+2 (t (6-77 -T
D. [[[B] - 3)
dr
[3] -
T([ - [([5) - B+1] ( (1- - 7)181 -)1811-1(11++11 (T)
(20)
As 3 + [3], we have
f(t = fI) + 1''
( ((lo)+)(()drr
lim 8 1-1
-0) (l) (l
(21)
(t)
In parallel to the integer-order derivative, given certain conditions ((Diehhel, 2010)[Leemma 3.13)), the Caputo fractional derivative possesses the semigroup property as illustrated in (Diethelm, 2010)[Lemma 3.13]:
D
(22)
Nonetheless, it is crucial to recognize that, in general, the Caputo fractional derivative does not exhibit the semigroup property, a characteristic inherent integer-order derivatives, as detailed in (Diethelm, 2010)[Section 3. 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, refer interested readers to (Diethelm, Theoree 3.17 and Remark 3.5.].
B.2 COMPARISON BETWEEN RIEMANN-LIOUVILLE AND CAPUTO DERIVATIVE
Another well-known fractional derivative is the Riemann-Loouville derivative, which, however, sees less use in practical applications (see Appendix B.4 for more insights). In this section, we offer succinet introduction to the Riemann-Liouville derivative and compare it with Caputo's 's definition. The
fractional derivative given as
dl0
S (1
D f(t) = r[[ - B) de m
(23)
Here again, we make the assumption that sufficient conditions are satisfied to ensure well-definedness (refer (Diethelm, 2010][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.
Classical Integer-order Taylor Expansion: (Diethelm, 2010)[Th 2.C] Assuming that f has absooltely continuous (m- 1)-st derivative, we have that for E [0, b).
m
th d f(0)
+ J d
f(t) =
dim ftt
(24)
h!
dek
where J" f²)( = Pon (t- f(T) dr. Note that here, is an integer.
Riemann-Liouve Fractional Taylor Expansion: (Diethelm, 2010)[hheorem 2.24] Let n > and m [n] 1. Assume that is such that j f has absolutely continuous (m-1)-st derivative.
21

---

Published as a conference paper at ICLR 2024
Then,
t r(n m
f + m
f(t)
f(0)
r(k - kk f(0) + J Dz (25)
Note that in the case n E N we have m and F(n m 1) [()) oc, and the first term
before the sum vanishes. Hence, we recover the classical result. For general n. the order in D
not a integer.
Caputo Fractional Taylor Expansion: (Diethelm, 2 3.8.] Assume that n 2 0, m
[n], and f has absolutely continuous (m - 1)-st derivative. Then
f(t)
Jr k!
(26)
k=
Note the order in D is integer. If compare (24) to (26), 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.
B.3 (CAPUTO) FRACTIONAL DIFFERENTIAL EQUATION
In this section, we first compare the initial conditions for FDEs under the Riemann Liouvvile 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 (Chamberlain al., 2021c; Rusch et al., 2022; Thorpe et al., 2022; Choi et al., 2023; Zhao et al., 2023a).
B.3.1 RIEMANN-LIOUVILLE CASE
Drawing from the Riemann- Riemil--iouvvlle fractional Taylor expansion, let assume that e is given function with the property that there exists some function g such that g = D.e. The solution of the Riemann-Liouville differential equation of the form
D
(27)
is given by
[ cyt
f(t) = c(f) +
(28)
where ej are arbitrary constants. In other words, to uniquely determine the solution from (25), we need to know the value of D f(0). This is akin to a k order ordinary differential equation a where the initial conditions are assumed as a f(),, with the distinction that the order in D.+nnm is not integer.
B.3.2 CAPUTO CASE
Similarly, if e is given function with the property that Dig g and if we intend to solve
DeS g
(29)
then we find
10 cyt [M] -j
f(t)= e(t) +
(30)
j=1
more, with ej as arbitrary constants. Thus, obtain unique solution, is natural prescribe the (
values of integer order derivatives f0), Dj (O).....
D 18111 -1 f(0) in the Caputo setting, mirroring
traditional ordinary differential equations.
22

---

Published as a conference paper at ICLR 2024
B.3.3 EXISTENCE AND UNIQUENESS OF THE (CAPUTO) SOLUTION
Next, we delve into a general Caputo fractional differential equation, presented follows:
D()(() g(l, v())
(31)
conjoined with suitable initial conditions. As hinted in (29) and (30), the initial conditions take the form:
(0) = (k) k = 0, ...... [ -
(32)
The following theorem addresses the existence and uniqueness of solutions:
Caputo existence and uniqueness theorem: (Diethelm, 2010][Theorem 6.8] Let [m-1) E I and h* > 0. Define the set G := [0, h*] x R and let the function g g:G- be continuous and fulfiil a Lipschitz condition with respect the second variable. i.c.
No ..... Vu
(z, g (z, v2)1 S L /11 -22
for some constant L > independent of z; 1/1, and 1/2- Then there uniquely exists function y E C [0, solving the initial value problem (31) and (32).
For point of reference, also provide the well-known (ieard-Linded uniqueness theorem for first-order ordinary differential equations.
Picard-Lindelff theorem (Hartman, [Page 8] Let D C R R be closed rectangle with (fo. yo) E int D, the interior of D. Let g D - R" be a function that continuous in t and Lipschitz continuous in y. Then, there exists some 0 such that the initial value problem
(t) ((t))(t)), y (tv) yu-
has a unique solution y(t) on the interval [Lo, to +
This allows 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, compositions of lincar 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 this, suppose softmax(u) Rax1 Then
01 (1 (1-
-1 U2
U1 Un
dv du
U201
U2
U2Unn
diag(v) TT
TV
Uy, V1
Un 02
Un (1
Un))
For bounded input, we have a bounded Jacobian. All the integer-order continuous GNN works, such recent contributions like (Chamberlain et al., 2021c; Rusch al., 2022; Thorpe al., 2022; Choi et al., 2023; Zhao et al., 2023a) assume the uniqueness of the ODE solutions. This means that all the integer-order continuous GNNs be extended by FROND framework with fractional dynamics.
B.4 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.
2. The Caputo fractional derivative maintains a more intuitive resemblance to the integer- order derivative and satisfies the sigmificant property of equating to zero when applied to constant. This property is not satisfied by the Riemann-Liouwy (Diethelm, 2010)[Example 2.4. and Example 3.1.] for further clarification.
fractional derivative. Refer
3. 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 (Diethelm, 2010: Diethelm et al., 2004; Deng, 2007).
23

---

Published as a conference paper at ICLR 2024
C NUMERICAL SOLVERS FOR FROND
We remind readers that numerous methods for training neural ODEs, and consequently updating the weights I in the neural network have been proposed. These include the autodifferentiation technique in PyTorch (Yan et al., 2018; Paszke et al., 2017), the adjoint sensitivity method (Chen et al., 2018b), and Snode (Quaglino al., 2019). In our work, we employ the most straightforward autodifferentiation technique for training FROND with fractional neural differential equations, leveraging the numerical solvers outlined in (Diethelm, 2010; Diethelm al., 2004; Deng. 2007). While we plan to investigate more sophisticated techniques for training FROND in future work, we have open-sourced our current solver implementations in ht tps: / httptuub.b...///. com/ ceenus//orchhe serve as valuable tools for the GNN community, encouraging the advancement of unique class of GNNs that incorporate memory effects.
We believe these will
In traditional integer-order continuous GNNs (Chamberlain et al., 2021c; Thorpe et al., 2022; Rusch al., 2022: Song et al., 2022; Choi al., 2023; Zhao al., 2023a), the time parameter serves as a continuous analog to GNN layers, resembling the concept of neural ODEs (Chen et al., 2018b) as continuous residual networks. Time discretization plays a crucial role in many numerical solvers for neural ODEs. For example, the explicit Buler scheme reduces neural ODEs to residual networks with shared hidden layers (Chen et al., 2018b). More sophisticated discretization methods, such as adaptive step size solvers (Atkinson et al., 2011), 1), 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 3 takes non-integer values. We present the fractional Bastuss-BBasffrtt-Muuluoo method with three variants utilized in this work, demonstrating how the time parameter continues to serve as continuous analog to the layer index and how the non-local nature of fractional derivatives leads to nontrivial dense or skip connections between layers. Addi- tionally, we also present one implicit LI solver for solving FROND when 3 is not an integer. It is worth noting that various neural ODE solvers remain applicable for FROND when 3 is an integer.
We first recall the FROND framework
x( F(W, X()), 8 > O,
where 3 denotes the fractional order of the derivative, and F is a dynamic operator on the graph like the models presented in Section 2.2. The initial condition set as X (0) =... = =()))= = X
consisting of the preliminary node features, akin to the initial conditions seen in ODEs.
C.1 BASIC PREDICTOR
Referencing (Diethelm et al., 2004), we first employ a preliminary numerical solver called "predictor" through time discretisation t, jh, where the discretisation parameter is a small positive value:
--
k-1 Σ Fo, (W)
(1) (b)
Σ
x (1)
(33)
j!
T(3)
==
where Mj,n- ((n-j)) -1 k denotes the discrete time index (iteration), and
j) t = kh represents the discretized time steps. (1) is the numerical approximation of X(t2). When 1, this method simplifies to the Euler solver in (Chen al., 2018b: Chamberlain al., 2021c) as Mj.n h. yielding (1) X (X +hF(W (k)))). Thus, our basic predictor can be considered as the fractional Fuler method or fractional Adams-Bashfone method, which is a generalization of the Fuler method used in (Chen et al., 2018b; Chamberlain al., 2021c). However, when 3 1,
need utilize the full memory (F(W. x k-1
1==0
The block diagram of this basic predictor, shown in Fig. 3, reveals that our framework introduces nontrivial dense or skip connections between layers. A more refined visualization is conveyed in Fig. 4, elucidating the manner in which information propagates through layers and the graph's 's spatial domain.
24

---

Published as a conference paper at ICLR 2024
C.2 PREDICTOR- PREDICTOR-CORRECTON CORRECTOR
The corrector formula from (Diethelm et al., 2004), a fractional variant of the one-step Adams- Moulton method, refines the initial approximation using the predictor p. p as follows:
[81-1 Σ j=0
(
(1)
1 (0) + ((3) Σ DrFF(WW xú) +
((3) Th F(W,
(34)
Here we show the coefficients 7j.n in the predictor-correctorr variant (34) from (Diethelm et al., 2004):
(k 1)0+1 (k 33
if 0,
h
(B) B(3 +1)
(k ( + (A: j]] -2(k - j)
if j k
if = k.
(35)
C.3 SHORT MEMORY PRINCIPLE
When T large, computational time complexity becomes a challenge due to the non-local nature of fractional derivatives. To mitigate this, (Deng, 2007; Podlubny, 1999) suggest leveraging the short memory principle modify the summation in (17) and (34) to Ej=n-k- This corresponds to employing shifting memory window with a fixed width K. The block diagram is depicted in Fig. 3.
memocy window width K
X(0)
X(6]
X(ta
Xito-i) Xit.)
X(ta)
X(t..)) X(C.)
Figure 3: Diagrams fractional Adams Bashforth-Moulonn method with full (left) short (right) memory.
X10)
X/ti
XIS)
X(t
time discretization
Figure 4: 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 FDFs are absent in ODEs. The weight of these skip connections correlates with H).x (B) as detailed (17).
C.4 SOLVER
The L1 scheme is of the most popular methods to approximate the Caputo fractional derivative in time. It utilizes a backward differencing method for effective approximation of derivatives. Refer- encing (Gao & Sun, 2011; Sun & Wu, 2006), we have the I.1 approximation of Caputo fractional
25

---

Published as a conference paper at ICLR 2024
derivative as follows:
k-1
x
R
x(1)
~
ki
1=0
where h is the temporal step size,
1
R kij (k- - 1) 1
0Sisk-1.
h2 3)
Applying the LI solver for our problem, we obtain
2
(x(tt-XXX) (A(X(1)) ()x
Manipulating the above equation, we have
k-2
(k) - =
(*-11)
Σ
R -x
1=0
The above formula is an implicit nonlinear scheme. To solve it without calculating the inversion of matrix, we propose the following iteration method:
(1) Compute a basic approximation of X(tx) with the following formula:
k-2
PX (R)
1 (A(( ()x = (1- Σ
n
R (x(+1)
1
(2) Substitute the above PX (k) into the implicit scheme to update ((k:;
k-2
x (6)
(A( 1)px (k-1) Σ
R
- -
(36)
The step (2) can be repeated multiple times to obtain accurate approximation of X(tx).
D DATASETS, SETTINGS AND MORE EXPERIMENTS FOR F-GRAND MODEL.
D.1 DATASETS
The dataset statistics used Table are provided Table Following experimental framework in (Chamberlain et al., 2021c), 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 fixed data splitting approach over the entire datasets, described in (Chami et al., 2019).
D.2 GRAPH CLASSIFICATION DETAILS
We use the Fake-NewsNet datasets from (Dou et al., 2021), constructed based on fact-checking information obtained from Politifact and Gossipcop. The dataset incorporates four distinct node feature categories, including 768-dimensiona. BERT features and 500-dimensional spaCy features, which are derived using pre-trained BERT and spaCy word2vee models, respectively. Additionally, 10-dimensional profile feature is extracted from individual Twitter accounts' profiles. Each graph within the dataset characterized by 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 the user retweeted the original news tweet, while an edge between two nodes is established when This hierarchical organization facilitates the analysis of the spread and influence both genuine and fabricated news within the Twitter ecosystem. The datasets statistics summarized in Table 6.
retweets the news tweet from another user.
26

---

Published as a conference paper at ICLR 2024
Table 5: Dataset Statistics used in Table
Dataset
Type
Classes
Features
Nodes
Edges
Cora
citation
7
1433
2485
5069
Citeseer
citation
3703
2120
3679
PubMed
citation
3
500
19717
44324
Coauthor CS
co-author
15
6805
18333
81894
Computers
co- to-purchase
10
767
13381
245778
Photos
co-pune purchase
8
745
7487
119043
CoauthorPhy
co-author
5
8415
34493
247962
OGB-Arxiv
citation
40
128
169343
1 1166243
Airport
tree-like
4
4
3188
3188
Disease
tree-like
2
1000
1044
1043
Table 6: Dataset and graph statisties used in Table 2
Dataset Politifact (POL)
Graphs (Fake)
Total Nodes
Total Edges
Avg. Nodes per Graph
314 (157)
41,054 314,262
40,740
131
Gossipcop (GOS)
5464 (2732)
308,798
58
D.3 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. implementing FROND, we employ one fully-connected (FC) layer on the raw input features to obtain the initial node representations X (0). Subsequently, as the decoder function to process the FDE output, X()), for executing downstream tasks. For more detailed information regarding the hyperparameter settings, kindly direct the readers the accompanying supplementary material, which includes the provided code for reproducibility. Our experiments were conducted using NVIDIA RTX A5000 graphics cards.
utilize another FC layer
D.4 LARGE SCALE OGBN-PRODUCT, DATASET
In this section, we extend our evaluation to include another large-scale dataset, Ogbn-products. adhering the experimental settings outlined in (Hu et al., 2021). For effective handling of this large dataset, we employ mini-batch training approach, which involves sampling nodes and constructing subgraphs. proposed by GraphSAINT (Zeng et al., 2020). Upon examination, F-GRAND-I demonstrates superior performance compared to both GRAND-I 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 (9). As such, the more advanced dynamic F(W, X (t)) in (6) may require additional refinement.
observe that
Table 7: Node classification accuracy(%) on Ogbn-products dataset
Model
MLP
Node2vec
Full-batch GCN GraphSAGE
GRAND-1 F-GRAND-I
Acc
61.06
72.49+0.10
75.64+00.21
78.29++0.16
75.56++0.67 77.25+0.62
D.5 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 Table 8. Specifically, we run F-GRAND on the Cora and Airport datasets with h and = 64. The solver variants perform comparably. For the Cora dataset, the fractional Adams Bashforth-Moulte worse than the other variants. However, it demonstrates comparable performance to other solver variants the Airport dataset.
method with short memory parameter of K 10 performs slightly
27

---

Published as a conference paper at ICLR 2024
Table 8: Node classification accuracy(%) under different solver when time 7 64
Predictor(17)
Predictor-Correcton (34)
Short Memory
Implicit L1
Cora(B3 0.6) Airport(3 = 0.1)
83.44+00.91
83,45 1.09
81.51+11.07 97.23+00.59
62.85+1.008 96.06+ 1.59
97.41+0.42
96.85 0.36
Table 9: Node classification accuracy based memory K on the Cora dataset when time T
40.
memory K
15
20
25
40
Accuracy (%)
74.9
80.8+:0.88
83.3:t1.11
83.9
84.2+11.1
84.1=1.2
84.5=11
84.1+11.1
8..8+1.1
Inference (ms)
9.81
17.53
24.97
32,03
38.79
42.99
45.27
48.70
48.35
D.5.1 FURTHER CLARIFICATION ON TWo ACCURACIES
This section aims 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 node classification. In contrast, numerical accuracy relates the precision of numerical solutions to FDEs, 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 Table 9. 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 the literature, to solve FDEs, there exist other more numerically accurate solvers like (Jin et al., 2017; Tian al., 2015; Lv & Xu, 2016) that higher convergence order. In general, these kinds of solvers can theoretically reduce computation cost and memory storage, obtain the same numerical accuracy using larger step sizes compared to lower-order solvers. 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.
D.6 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 fractional integration.
In this section, we report the inference time of the different solver variants in Tables 10 to 13. For comparison, consider the neural ODE solver for 3 1. which includes Euler, RK4, Implicit Adams, and dopri5 methods as per in the paper (Chen et al., 201 2018b). We observe that when T = 4, the inference time required by the FROND solver variants is similar that of the ODE Fuler solver. However, for larger 7' = 64, the basic Predictor (17) solver requires more inference time than Euler and is comparable RK4. For more accurate approximation solver variants (34) and (36) incorporat- ing the corrector formula, Tables 12 and 13 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.
D.7 CONTINUED STUDY OF OVERSMOOTHING
To corroborate that FROND mitigates the issue of oversmoothing and performs well with an increasing number of layers, conducted an experiment employing the basic predictor with up to 128 layers in
28

---

Published as a conference paper at ICLR 2024
Table 10: Average time under different solvers when time 7 4 and hidden dimension 64 on Cora dataset
Predictor(17)
34)
Short Memory
Implicit L.I
Fuler
RK4
Implicit Adams
dopri5
Inference time (ms)
0.98
1.67
0.98
0.62
0.96
2.06
3.20
11.91
Table 11: Average time under different solvers when time T 64 and hidden dimension is 64 on Cora dataset
Predictor(L')
Predictar-Coreettoft Correcton33)
Short Memory
L1
Euler
KK4
Implicit Adamss
dogriS
Inference
46 160.92
30.26 26
221.74
42.66
103.45
66.15
Table 12: Average time of (34) and (36) with correctors, used to refine the approximation, when time 4 and hidden dimension is 64 on the Cora dataset.
Predictor-Corrector (34)
5
10
Inference time (ms)
1.67
3.31
4.74
8.34
Implicit-L.1 (36)
3
5
Inference time (ms)
0.62
1.04
1.48
2.55
Table 13: Average time of (34) and (36) with correctors, used refine the approximation, when time
T' = 64 and hidden dimension is 64 on the Cora dataset.
Predictor-Correctorr (34)
3
Inference time (ms)
160.92
442.88
Implicit-LI (36)
3
Inference time (ms)
221.74
441.60
the main paper. The results are presented in Fig. 2. 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 (Chami et al., 2019).
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 Table 14. From our observations, F-GRAND-I maintains a consistent performance level even as the number of layers escalates. This contrasts with GRAND-I, where there is notable performance decrease with the increase in layers. For instance, on the Cora datasets, the accuracy of GRAND-I drops from 81.29% with 4 layers to 73.37% with 256 layers. In stark contrast, our F-GRAND-I model exhibits minimal performance decrease on this dataset. On the Airport dataset, F-GRAND-I registers slight decrease 94.91% with 256 layers from 97.0% with 4 layers. However, the performance of GRAND-1 significantly drops to 53.0%. These observations align with our expectations, as Theorem 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.
D.8 ABLATION STUDY: SELECTION OF 3 CONTINURD
In the main paper, explore the impact of the fractional order parameter 3 across variety of graph datasets, with the results of these investigations presented in Table 3. More comprehensive details concerning the variations in 3 can be found in Table 15.
29

---

Published as a conference paper at ICLR 2024
Table 14: oversmoothing mitigation under fixed data splitting without LCC
Model
GCN
$1.35=1.27
15.3
19. 19.70:1:7000 06
21 89
13,0:10.
13010000
130100
13,0:100
Corn
GRAND
GAT
01.29-0.43
20.95=2.28
82.95:00.52
11:00=000
82.48:10.16
31 10:00:000
61.72:1:0.55
31.50:100
$1.33:10.22
31.90=00.0
81,07:00.44
31,90:0000
8009:1043
31.00200
73.37:10.59
31,902000
F-GRAND-1
01.17=0.75
82,68:0.64
83.05::0.81
K2 82.90:1881
83.4620.91
82.85::0.89
62.341:0.33
81.74:20.53
GCN
68.84=2.46
61,5822.09
10.64::179
7.7400
7.74000
7.740..0
1:7000
7.7.000
Cheseer
GRANDA
GAT
70.08=122
0.57
70,39=0.68
18:10=000
20.18-00556
18.10+300
18:00000
68.1-:1447
18.10:0.0
67.44+
18.10400
63.45-42266
18.00+00
56.9-+
18:10=0000
F-GRAND--1
70.68 1.23
71 71.0420668
71.00 1.12
70
70.27:20.86
70.50
71.02:445
GCN
14.77=1455
74.43=8.199
62.5622.16 2616
15.27:0.0
15.27=00.0
15.274:0.0
15:27:000
15.27=0.00
Airport
GRAND
GAT
33.59=151 1.51
$0.53=9.58
67.02-4.700
19.88-9967
1:.24:133500
46.56+00000
40.56-00000
67 1:402
62.29 61.28:110.8 10,83
0.656=0000
50.38::..55
46.26--000
57.96-:11.63
#6.6-+00UU
14.0
46.36=0000
F-GRADD
97.0200799
97.0928877
96977-0.44
9650000.00
97.41-0422
96.53-20.74
97.03400.55
34.71
Table 15: Node classification accuracy(%) (%) under different value of / when time T
8.
8
0.1
0.2
0.3
0.4
0.5
066
0.7
0.5
1.0
Coru
7480+0042
76.10=0.34
77.0+-0.88
07.80=0.75
79.90=0..11
50.79:00.58
81.56400.00
82.4++051
92.63
82.37.00.99
Aippst
97.09+0.67
9666-0091
05.80++22.3
94.0
91.56-+6..4
79
84.36+-8004
79.29-6.011
78.73+6.33
78.88+9667
D.9 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 (Zügner et al., 2018). Adversaries, aiming deceive trained GNN, can either introduce phase, known as an injection attack (Wang et al., 2020; Zheng et al., 2022; Zou et al., 2021; Hussain et al., 2022), manipulate the graph's 's topology by adding or removing edges, termed as a modification attack (Chen et al., 2018a; Waniek et al., 2018; Du et al., 2017). 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 (Zügner & Günnemann, 2019). Our approach adheres to the attack setting described in Pro-GNN (Jin et al., 2020), and we utilize the perturbed graph provided by the DeepRobust library (Li et al., 2020b) to ensure fair comparison. The perturbation rate, indicating the proportion of altered edges, is incrementally adjusted in 5% steps from 0% to 25%
nodes into the graph during the inference
The results of these experiments are presented in Table 16. 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, 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 theoretical understanding are detailed in our recent work (Kang et al., 2024).
Table 16: Node classification accuracy (%) under modification, poisoning, non-targeted attack (Metattack) in transductive learning.
Dataset
Ptb Rate(*)
GGN
GAT
GRAND-ul
F-GRAND-nl
0
83.50=0.44
83.97=0.65
83.14=1.06
83.48+11.00
76.55=0..79
80.44=0.74
80.54=1.17
80.25++0.00
10
7..3
75.61
7..59
77.94.10.48
Cora
15
65.10:10.71
6..78
71.62+1.39
14.14
20
59.56++2.72
59.94-0.92
1..5
69.04411.13
25
17.53
54.78=0.74
53.70=1.91
63.40
0
71.5610.55
73.26=0.83
1..40=1.08
70.1440..8
Citeseer
70.88+00.62
72.89-0.83
70.99=1.12
70.0+1.72
10
67.55+0.88
70.63+0.48
68.83-1.31
68.64+1.111
15
64.52+1.111
1.002+1.09
66.78-0.92
67.90+0.41
20
63.03
61.04
18.95
65,84.10.75
25
56.9411..09
61.85
60.52.1.29
66.50.11.16
30

---

Published as a conference paper at ICLR 2024
D.10 COMPARISON BETWEEN DERIVATIVE
(RL) DERIVATIVE AND CAPUTO
The underlying rationale for opting for the Caputo derivative over the fiemann-Liouviinn (RL) derivative is extensively delineated in Appendix B.4. However, supplementary experiment was conducted utilizing the RI. derivative in lieu of the Caputo derivative, the results of which are documented in Table 17. 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.
Table 17: Comparison between RL-GRAND-I (using Riemann- Liouville derivative) and the original F-GRAND-I (using Caputo derivative).
Method
Cora
Citeseer
Pubmed
CoauthorCSS
Computer
Photo
ComuthorPhy
Airport
Discuse
GRAND4
83.6110
73.42:..5
78.8+1.7
92.910.4
83.7=1.2
92.3:0.9
93.5=0.9
80.549.6
74.5:3..
RL- GRAND-
84.61122
74,2110
80.1112
92.810.3
87.4=1.1
93.310.7
94.110.3
90.71113 96.2.0.22
F-GRANDA
84.8+111
74.0411.
79.4
93.041.3
84.4 1.5
92.8-10.6
94.5-10.4
98.1-10.2
92444399
D.11 FRACTAL DIMENSION OF GRAPH DATASETS
Table 18: Comparison between the estimated fractal dimension, the best order B and the 6- hyperbolicity
Dataset
Disease
Airport
Pubmed
Citeseer
Cora
fractal dimension
2.47
2.17
2.25
0.62
1.22
best 3 (F-GRAND-I)
0.6
0.5
0.9
0.9
0.9
best 3 (F-GRAND-nl)
0.7
0.1
0.4
0.9
0.9
9-hyperbolicity
0.0
1.0
3.5
4.5
11.0
In Fig. 5, using the Compact-Bow Beeett-ox-Bumiing algorithm from (Song et al., 2007), we compute the fractal dimension for some datasets that have moderate sizes. As noted in Table 1, there is clear trend between d-hyperbolicity (as referenced (Chami et al., 2019) for assessing tree-like structures-wer lower values suggesting more tike graphs) and the fractal dimension of datasets. Specifically, lower 9-hyperbolicity corresponds to larger fractal dimension. As discussed in Sections and 4, we believe that our fractional derivative D effectively captures the fractal geometry in the datasets. Notably, we discerned trend: larger fractal dimension typically corresponds to a smaller optimal B.
E MORE DYNAMICS IN FROND FRAMEWORK
E. REVIEW GRAPH ODE MODELS
The work by (Thorpe et al., 2022) introduces graph neural diffusion with source term, aimed at graph learning in scenarios with a limited quantity of labeled nodes. This approach leverages subset of feature vectors, those associated with labeled nodes, indexed by and considered "trustworthy" to act as source term. It adheres to (4) and (5), incorporating an additional source term, facilitating the propagation of information from nodes I to node i.
GRAND++:
dx(t) = F(X(t)) + s s(xx)er) dt
(37)
Here, I denotes the set of source nodes, s(.) represents source function, and F(-) embodies the function depicting the right-hand side of (4) and (5). The model is manifested in two variations, respectively denoted as GRAND++nl and GRAND++11
GraphCON: Inspired by oscillator dynamical systems, GraphCON (Rusch et al., 2022) is defined through the employment of second-order ODEs. It is crucial to highlight that, for computation, the
31

---

Published as a conference paper at ICLR 2024
Airport Slope: 2.17
Disease Slope: 2.47
Pubmed Slope: 2.25
Cora Slope:
Claaseer Slope: 0.62
Number of Boxes
Size
Figure 5: The fractal dim of datasets. We use the Compact-Boo Box- Box-Burning algorithm in (Song et al., 2007) to compute the log-log slope (fractal dim) of the box size and the minimum number of boxes needed to cover the graph.
second-order ODE is decomposed into two first-order ODEs:
dY dt - o - - -
dX(t) Y(t),
(38)
dl
where a--) is the activation function, Fg (x()), is the neural network function with parameters ? and â are learnable coefficients, and (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-I) and non-linear (GraphCON-nl) versions concerning time. The differentiation between these versions determined by whether the function Fş undergoes updates based on time t.
With the objective of addressing heterophilic graphs, the paper (Zhao et al., 2023a) 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 preserve diversity in heterophilic graphs. The corresponding formula is illustrated in (39).
CDE:
dx(t) (A(X()) ()x(t) + div (V(() X(t))
(39)
df
In this equation, (it) represents the velocity field of the graph time t, div(-) denotes the divergence operator as defined in the paper (Chamberlain et al., 2021c; Song et al., 2022), and symbolizes the element-wise (Hadamard) product.
GREAD: To address the challenges posed by heterophilic graphs, the authors in (Choi et al., 2023)
present the GREAD model. This model enhances the GRAND model by incorporating reaction term, thereby formulating diffusion-reaction equation within GNNs. The respective formula is depicted in (40), and the paper offers various alternatives for the reaction term.
dX (t) = aL(((t)) (x)))
(40)
dt
In this equation, r(X(t)) represents the reaction term, and a is a trainable parameter used to balance the impact of each term.
32

---

Published as a conference paper at ICLR 2024
E.2 F-GRAND++
Building upon the GRAND++ model (Thorpe et al., 2022), we define F-GRAND++ as follows:
((x(()==(((())) s(xx)ez
(41)
We follow the same experimental settings as delineated in the GRAND++ paper. Given that the primary focus of GRAND++ is the model's performance under limited-label scenarios, our experi- ments also align with this setting. The sole distinction lies in the incorporation of fractional dynamics. Within this framework,
substitute the ordinary differential equation dx(() used in GRAND++
with our FROND fractional derivative (t). The optimal is determined through hyperparameter tuning. When
1, F-GRAND++ seamlessly reverts GRAND++, and the results from the original paper are reported. Our observations distinctly indicate that the Fractional-RRANDA con- sistently surpasses the performance of the original GRAND++ in nearly all scenarios. We also present the complete comparison results in Table 20, where it is evident that F-GRAND++ demonstrates greater effectiveness in learning with low labeling rates compared GRAND++, GRAND, and other baseline methods.
Table 19: Node classification results (%) under limited-label scenarios
Model
pre class
Cora
Citeseer
Pubmed
CoauthorCS
Computer
Photo
GRAND++
1
$4.94.1. 16.09
19.59
65.94
60.30.1
67.65 0.337
83.12.10.78
F-GRAND++
1
57.31 -8.89
69.11 0.95
2.72
67.71 0.7
67.65
83.12.10.78
B
0.95
0.85
1.0
1.0
GRAND++
2
66.92.1 10.04
64.98 1.331
69.31-1..88
1..5
74.47:11.44
83.71
F-GRAND++
2
70.09=8.36
64.98+8.31
69.37-5.36
77.97=2.35
78.85=0.96
83.71=0.90
0.9
1.0
0.95
0.5
0.8
1.0
GRAND++ F-GRAND++
77. 4.44
70.03433.6
71.99=1.11
84.83=00.84
82.64-0..56
18.33+1.21
78.79111.66
70.26
7.678
86.09.12.09
82.64
88.5610.66
5
0.9
0.8
0.9
0.8
1.0
0.75
GRAND++
80.86:1.2.9
72.34:1:2.2
33.18
86.94:10.46
0..89
90.65
F-GRAND++
10
82.7310.88
73.52.1.444
77.15
87.85.1..44
83.26.10.41
91.1510.55
B
0.95
0.9
0.95
0.6
0.7
0.95
GRAND+
1..7
33.53:13.31
1.337
0..8
93.55 0..7
F-GRAND++
20
84.57=1007
1..81=1.78
79.96=1.66 1.68
91.03=0072
85.78=0.43
93.55=0038
B
0.9
0.85
0.95
0.9
0.9
1.0
E.3 F-CDE
Drawing inspiration from the graph neural CDE model (Zhao et al., 2023a), we further define the F-CDE model as follows:
()x(() + div(V(()) x())
(42)
In this expression, V(t) represents the velocity field of the graph time t. The divergence operator, div(.), is defined as per the formulation given in (Song et al., 2022), and symbolizes the element- wise (Hadamard) product.
We follow the same experimental setting in the CDE paper(Zhao et al., 2023a). Given that the primary focus of CDE is on evaluating model performance on large heterophilic datasets, our ex- periments are also conducted under similar conditions. The statistics for the dataset are available in Table 21. The sole distinction our approach lies in incorporating fractional dynamics; we achieve this by replacing the ODE used in CDE with FROND fractional derivative. The complete com- parison results in Table 22 conspicuously reveal that Fractional CDE exhibits superior performance compared to the conventional CDE and other baselines across various datasets.
33

---

Published as a conference paper at ICLR 2024
Table 20:
Full table: Classification accuracy of different GNNs trained with different number of labeled data per class (#per class) on six benchmark graph node classification tasks. The highest accuracy is highlighted in bold for each number of labeled data per class. These results show that F-GRAND++ is more effective in learning with low-labeling rates than GRAND++ and GRAND. Where available, baseline results are cited from (Thorpe al., 2022).
F-RRAN++
78.79
70.26
2,36
73.38
2.00
0.56
1:
0.81
77.15
$4.57 1.07
74.81
1.78
79,96
1.48
91.03
0.72
85.78
10.43
93.55
0.38
54,91 1 15.99
55.95
1 9.50
65.94 4.87
00.30
1.50
67.65
0.37
83.11 0.78
66.92
10.04
64.98
8.31
09.31
4.87
76.53 = 1.85
76.47 1: 1.48
83.71
1: 0.90
GRANDH
77,90
4.46
70,0 4: 3.63
71.99 4: 1.91
54.83 0.54
82.64 4: 0.50
89.33
10
$0.86 + 2.99
72.54
2.42
75.13 + 5.88
86-94
0.46
82.09
0.81
90.65
20
$2.95
1.37
73.53 + 3.31
79.16
1.37
90.80
0.34
85.33 + 0.50
93.55
0.38
59.15
1.66
11.16
5.58
1.85
2.10
1.09
78,76
00.85
14.01
59.06
9.76
60.45
16.20
83.61 1.49
76.90 1.49
$3.61 1.
GCN
73 86 1.7
67.24
4.19
7.93
86.65
0:43
82.47
0.97
$8.86 1.56
10
78.82 5.38
72.18
3.47
72.59 + 3.19
88.00 0.50
82.53 + 0.74
90.41 # 0.35
20
52 07 2,03
74.21 + 2:00
76.89
3.07
91 09
0.55
82.94 + 1:54
91.05 0.11
47.85 1: 15.38
50.31
14.27
55,04
1231
51.13 3.24
27.14 1: 7.81
13.55 5.15
59.33 13.55
$5.55
9.19
60.24
14.44
63.12 6.00
65.07 8.86
76,89
GAT
71 (4 5.74
67.57
500
68.54
5.75
71.65 4.53
71.43 7.34
$3.01
10
76. 31 4.87
71.35
1: 4.92
72.44
3.50
74.71
3.35
76.04 1: 0.35
42 2.38
2.28
2.90
1
79.95
2.88
80.05
43.00 14.11
11.45
12.71
M1.35 1.35
27
GruphSage
91.33
79.98
47.72 15335
39.13
56,47
$5.99 5:17
23.78
34,72
60.85 +: 14.01
48.52 2: 9-52
61.03 6,93
76.57 4.06
38.19
3.72
43.013 +
MoNet (Monti at.., 2017)
73,86 7.07
66
67.92:45 2.50
K7 02 1.67
59
4.75
71.
10
78,82 5.38
68.08 1: 6.29
71.24 1: 1.54
58.76 0.49
68.66 1 3.30
78.66 3.17
20
$2.00
2.03
71.52 4.11
76,89 1.75
90.31 0:41
73,66 1 2.87
$8.61 1. 1.18
Table 21: Dataset statistics used in Table 4
Dataset
Nodes
Edges
Classes
Node Features
Roman-empire
22662
32927
18
300
Wiki-cooc
10000
2243042
5
100
Minesweeper
10000
39402
2
7
Questions
48921
153540
2
301
Workers
11758
519000
2
10
Amaon-ratings
24492
93050
5
300
Table 22: Full table: Node classification accuracy(%6) %) of large heterophilic datasets.
Model
Roman-empire
Wiki come
Mineswesper
Questions
Woekers
Amazon ratings
ResNet
65. 14.71++0.44
89. 36-60.71
50.95+1.122
70.10+0775
73.08+ 1.28
45 70-00.699
H2GCN/Zha 2020a)
66.66+
CPGNNiZha et 2021)
63.78:10.50
84.84:30.66
71 1.14
67.09::2.63
72.44:0.80
44. 34.36:00.35
GPR UPR--NNNiChien et 2020)
73 37++0.68
91 90
81 79+00.9
73.41=:1.24
70.59:11155
43 11.0-00.44
GloGNN(L 2022)
63 05.54.0..9
98.49-00.45
62 03.53+1133
67.15+:192
73:9==0.5
17.28-00.6
FAGCN(Bo 202))
70.53:10.99
88:30.37
89
77,04:1.56
81.87:20.94
32:12.5550
GIR--GNN(Di 2022)
75 87-40044
97. 81-33333
X3 56+0084
72,99+-1.00 1.05
78.06=90991
41 41.47:0.51
ACM- LCM-GCNiLumn 2022)
68. 1..35511..95
87.48:11.00 106
90.1710.57
0OM
78.25:10.78
38,51:13,38
GRAND AADhChamberlain
92.03-10.46
76.67:40.8
20.67+:1288
75.33+0084
15.05-29066
SomphiBel(Sogg et al., 2022)
69 07.44771
90. 30:40.50
76.51
70.79:00999
73.02:00:2
43 13:63:00.42
Diag-NSD(RRdiar 2012)
77 50++0.67
92.06-+-044
89 59+00.61
69.25+11.15
79.81-101999
37 96---02220
ACMP al.,
02.68-11.377
CDE
91. 01.64-4-0.28
97. 99.99-0..38
95 50:45.23
75.17::0.99
80. 20.:0
47. 47.63-10.43
FCDD
93.06- 1-.55
V8.73-H166
96.04--0255
75.17:0.99
82.6N=-0.66
49.01
F-CDE
0.9
0.6
34

---

Published as a conference paper at ICLR 2024
E.4 F-GREAD
Our FROND framework is also extendable to the GREAD model (Choi et al., 2023), as defined in (43).
x -a)((()) +
aar
(43)
where r r(X(()) represents reaction term, and is trainable parameter used emphasize each term.
We adhere to the same experimental setting outlined in the GREAD paper (Choi et al., 2023), concentrating exclusively on heterophilic datasets. We choose the Blurring-Sharpnning (BS) as the reaction term formulate both GREAD-BS and F-GREAD-BS, as GREAD-BS exhibits strong performance according to Table 4 in the GREAD paper (Choi et al., 2023). The results presented in Table 23 (refer to Table 24 for comprehensive comparisons with other basclines) demonstrate that our FROND framework enhances the performance of GREAD across all examined datasets.
Table 23: Node classification accuracy(%) of heterophilic datasets
Model
Chameleon
Squirrel
Film
Texas
Wisconsin
GREAD-BS
11.3
59.22+ 1.00
37. 17.9
89.41 88.92+:3.72 92.97+4.399
F-GREAD-BS
7..45+1.98
38.78
9..89++3.80
B
0.9
0.9
0.8
0.9
0.9
Table 24: Full table: Node classification accuracy(%) of heterophilic datasets
Model
Chameleon
Squirrel
Film
Texas
Wisconsin
Geom Cooei al., 2020)
60.00=2.81
38.15=00.92
1.15
66.76:2.72
64.51
H2GCN(Zhu et al., 2020a)
60.11=2.15
36.48=1.86
1.07
84.86=7.23
87.65:1.998
GGCN(Yan et al., 2022)
71.14=1.80
15.17
11.56
84.86-44.55
86.86++3.29
LINKX(Iim et al., 2021)
68.42= 1.38
61.81+1.80
16.55
74.60+8.37
75.49+45.72
GloGNN(I.i et al., 2022)
69.78-2.42
17.54
37.3
84.32-44.15
87.06+3..53
ACM-GCN(Luan et al., 2022)
66.93
14.40.11.88
3..28-11009
87.8411.40
88.43113.22
GCNlI(Chen et al., 2020)
63.86-3.04
38.47
37.4411.30
77.5743.83
80.39-13.40
CGNN(Xhonneux al., 2020)
46.89=1.66
29.24=11009
35.95--0.86
71.35+4.05
7..31
GR GRAND(C (Cammberlain et al., 2021a)
54.67=2.54
10.05=1.550
35.62+
75.6
79.41
BLENDDChambeele et al., 2021b)
60.11-2.00
13.06
35.63
83.24
84.12.13.56
(oonar et al., 2022)
68.04 1.58
56.34.11.2
37.81
85.05
89.41114.74
GRAFF(Di Giovanni et al., 2022)
71,08 1.75
14.52
36.09
88.38
87.45
GREAD-BS
71.38=1.31
59.22=11.4
17.90=
88.92-3..72
89 39.41
F-GREAD-BS 8
71.45=1.98
60.86
38.28
92.97-4.39 0.9
3..59-3.80
0.9
0.9
0.8
0.9
E.5 F-GRAPHCON
We also incorporate the following fractional-order oscillators dynamics, inspired by (Radwan et al., 2008; Rusch et al., 2022):
D = a (Fo(X,t) -X aY x
(44)
which represent the fractional dynamics version of GraphCON (Rusch et al., 2022). We denote this as F-GraphCON, with two variants, F-GraphCON-GCN and F-GraphCON-GAT Here, Fo is set as GCN and GAT, in the setting described in (Rusch et al., 2022). We refer readers (Rusch et al., 2022) for further details. Notably, when 3 = 1, F-GraphCON simplifies to GraphCON, devoid of memory functionality.
35

---

Published as a conference paper at ICLR 2024
Table 25: Node classification accuracy(%) based on GraphCON model
Cora
Citeseer
Pubmed
Airport
Disease
GraphCON-GON
81.9+11.7
72.9+2.1
78.8+
68.6+2..1
87.5+t4.1
GraphCON-GAT
83.2
73.24:188
1..4
74.122.7
65.715.9
B-GGrahhCON-GCN
1.4
7.13++1.1
8..3+1.3
0.5
92.1 2.8
0.9
0.8
0.9
0.1
0.1 86.914.0
F-GraphCON-GAT
83.911.2
13.4
1..4
97.3 0.1
0.7
0.9
1.0
0.1
Table 26: Full table: Node classification accuracy (%) based on GraphCON model.
Cora
Citeseer
Pubmed
Airport
Disease
GCN
1..
11.9
77.8
81.6=0..6
69.8+0.55
GAT
81.8++1.3
1..
78.7+2..3
81.6++0.4
70.4+-0.5
HGCN
7..7+110
65.8+2.0
76.4+-0.8
85.4++0.7
89.9+1111
GIL
82.1
71.1
77.810.6
91.5+11.
90.8
GRAND-I
83.611.0
73.4
78.811.77
80.519.6
74.5+13.4
GRAND-nl
82.311.6
70.9
77.5
90.9
81.016.7
GraphCON-GON
81.911.7
72.942.1
78.8+2.6
68.6+2.1
87.511.11
GraphCON-GAT
83.2#1.4
73.2
79.4
74.1#2.7
65.7+5.9
F-GraphCON-GNN
84.61:14
75.311111
80.3+1.3
97.310.5
92.1=2.8
0.9
0.8
0.9
0.1
0.1
F-GraphCON-GAT
83.9+1122
1.34+ 13.4++1.5 0.9
97.3++0.8 0.1
86.9+-400
0.7
1.0
0.1
Table 27: Node classification accuracy(%) of undirected graphs based on F-FLODE model
Film
Squirrel
Chameleon
FLODE
37.16
64.23
13.60
F-FLODE
37.95-11.27
65.53 1.83
74.17
0.8
0.9
0.9
Table 28: Node classification accuracy(%) of directed graphs based on F-FLODE model
Film
Squirrel 74.03+1.55
Chameleon
FLODE
1..4
77.98
F-FLODE
37.9711.11 0.9
7..03 0.9
78.51
0.9
E.6 F-FLODE
In the work of (Maskey et al., 2023), the authors introduce the FLODE model, which incorporates fractional graph shift operators within integer-order continuous GNNs. Specifically, instead of utilizing Laplacian matrix L, they employ the fractional power of L, denoted as L (see (45)). Our research diverges from this approach, focusing on the incorporation of frectfrationnl derivative DJ for updating graph node features in a memory-inclusive dynamical process. It is pivotal to differentiate the term "fractional" as used in our work from that in (Maskey et al., 2023). they signify fundamentally distinct concepts in the literature. Fundamentally, FLODE differs from our work in key aspects:
36

---

Published as a conference paper at ICLR 2024
FLODE employs the fractional (real-valued) power of L, namely L The feature evolution model used by FLODE, specifically in its first heat diffusion-type variant, is given by:
dx = -L dt
(LOOE)
This is graph spatial domain rewiring technique, as L introduces dense connections compared to L. As a result, FLODE introduces space-based long-range interactions during the feature updating process.
In contrast, our FROND model incorporates the time-fractional derivative Dj to update graph node features memory- dynamical process. In this context, time acts as continuous counterpart to the layer index, leading to significant dense skip connec- tions between layers due to memory dependence. Thus, FROND induces time/layer-baset long-range interactions in the feature update process. Note that FLODE does not utilize time-fractional derivatives. Our method is not only compatible with various integer-order continuous GNNs, including FLODE (see (F-FLODE)), but also extends them graph FDE models.
We next introduce the F-FLODE model, which utilizes time-fractional derivatives for updating graph node features FLODE:
xt) -L
(F-FLODE)
where denotes the symmetrically normalized adjacency matrix. The a-fractional power of the graph Laplacian, L given by:
=
(45)
In this formulation, U, E, and are obtained from the SVD decomposition of L UEVH and E RR represents the order. The channel mixing matrix 4, a symmetric matrix, follows the setting in (Maskey al., 2023).
Following the experimental setup outlined in (Maskey et al., 2023), we present our results in Tables 27 and 28, demonstrating that our FROND framework enhances the performance of FLODE across all evaluated datasets. Note the difference in the equations in (FLODE) and (F-FLODE), where the two are equivalent only when the FLODE model as special case when 3 = 1. Our experimental results indicate that F-FLODE outperforms FLODE with the optimal 3 o 1 in general.
1. This example illustrates that the FROND framework encompasses
F PROOFS OF RESULTS
In this section, we provide detailed proofs of the results stated in the main paper.
F.1 PROOF OF THEOREM
Proof. We observe that for 8 < they possess the propertics, the coefficients ck. bm defined in (10) satisfying the following properties (Gorenflo et al., 2002).
Ch = 1, =
0
1=1
60 = 1, bm =
Ck. 1 = bn > b
b2 bg
+ 0.
k=1
k=
From the definition of the transition probability (11), we have
37

---

Published as a conference paper at ICLR 2024
P(R(In+1) = xx)
(,(P(R((t) = xni + en P((())) = xxn +... + 022P(R((nn xn)+
+ + - -o))P((((())) = xn) + ' 00
WiAh
P(R(In) =xx)
dj
j=1
-,,((((10) = -xh) + + (.,(((())) = Xn. +...
+ 02
Pe ((R(((n(1) =xn)+
N
Wigh P(R((()) = =x;). (46)
+ (R)(( = X1) - (P((((n) = xn) +
Σ
dj
==1
By rearranging, have
P(R(fx+1)= = xxi)
Σ
C&P
(((((tn1)) =x)) b,,P(R((to) = -xn)
k=1
3 ())
(11)k+] (3)
={-1)
(()((n-1)) =x,) (11))&t1
P( (1(( xt)
k.
k=1
n
1-I)*
P(R =xn)
k
k=
B
Σ 1) &
P(R P(R = =x,)= Σ
(-1)*
P(R xni)
k
33)
1- *
[P(((tn+1-R) xh) P(R xx,]]
k)
= - ((R(tn) = x1) + Σ P()(((())==) = =x,).
d2
Dividing both sides of the final equality by 0,, it follows that
P(R(tn+11- x) xx) P(R Xt.)
k=
(((tt,) = X + j P(((t) xx))
(47)
From the Seiinwald-Letnkeee fractional derivatives formulation (Podlubny, 1999) (2.54)), the limit of LHS of (47) is
F(R(tn+1-s)= xn) P(R xn) = D.P(((() = x = PD))
0
k
To K=D
(48)
where Ptt) lim ((R(tn)) and [D(P(t)) denotes the h-th element of the vector. On the other hand, the RHS of (47) is
W P(R(tt)) xx)) = [L)(((((())))
-P(((t,) xx) +
(49)
where ((((t,)) is the probability (column) vector with j-th element being P(R(t,) = xj), and (-LI((((t))]] denotes the h-th element of the vector (Rt.))
Putting them together, we have
P()) = -LP(t
(50)
since assume l,, in the limit. The proof of Theorem I is now complete.
38

---

Published as a conference paper at ICLR 2024
F.2 PROOF OF COROLLARY
It directly follows from the the linearity of FDEs and ( X = Σ, P(0)x, where recall that
the initial probability vector .P(R(0)) .P(0) is represented as one-hot vector with the i-th entry marked as 1.
F.3 PROOF OF THEOREM
Before presenting the formal proof, we aim provide additional insights and intuition regarding the algebraic convergence from two perspectives.
Fractional Random Walk Perspective: In a standard random walk, a walker moves to new position each time step without delay. However, in fractional random walk, which is more reflective of our model's behavior, the walker has probability of revisiting past positions. This revisitation is not arbitrary; power-law distribution with a long tail. This characteristic fundamentally changes the walk's dynamics, introducing memory component and leading to slower, algebraic rate of convergence. This behavior is intrinsically different from normal random walks, where the absence of waiting times facilitates quicker, exponential, convergence.
governed by waiting time that follows
Analytic Perspective: From an analytic perspective, the essential slow algebraic rate pri- marily stems from the slow convergence of the Mittag-Lefiih function towards zero. To elucidate this, let us consider the scalar scenario. Recall that the Mittag-Lefh Leffler function E3 is defined as:
Ea(=) = Σ
T(J +
1=
for values of z where the series converges. Specifically, when 3 = 1,
80
Ei(2) j=0 r 1
=
exp(2)
j!
corresponds the well-known exponential function. According to [AI, Theorem 4.3. the eigenfunctions of the Caputo derivative are expressed through the Mittag-Leffte function. In more precise terms, if define y(t) as
y(() = E. (-At"), 12 0,
it follows that
(tt( = -
Notably, when 8 this reduces dexp(-At) --xexp(-XX We examine the behavior of Eg (-X From (Diethelm, 2010)[Theorem 7.3.], when 0 < 8 < 1, it is noted that: (a) The function y(t) is completely monotonic on (0, 20). (b) As x > oo,
-ー
v(() = XT(1 B) (1+0(1))
Thus, the function Eg St converges to zero at a rate θ (t- Our paper extends this to the general high dimenssonl case by replacing the scalar À with the Laplacian matrix L, wherein the eigenvalues of L play a critical role analogous to in the scalar case.
For diagonalizable Laplacian matrix L, the proof essentially reverts to the scalar case as outlined above (refer to (56) in our paper). However, in scenarios where L non- diagonalizable and has general Jordan normal form, it becomes necessary to employ the Laplace transform technique to demonstrate that the algebraic rate remains valid (refer to the context between (56) and (58)).
39

---

Published as a conference paper at ICLR 2024
Proof. We first prove the stationary probability π
by induction. Assume
Σ d
Σ 21=1 2z
that for i I,..., n, the probability distribution P(R(t,)) always equals π.. For i n + 1, from (46), it follows that
(P(((1+)))] = P(R((n+1)) = =xh)
=b,P(R((to) =x,) + Σ C+
(((((((+)) xx.)
N
Wihh ((((t))) =x,)
P((((t)))) = xx)) +
0
dj
/=1
N
Wih
= =TAbR* +
Thcck
TAO
T
dj
j=1
dj
0 Wih
(bn
(k) T
d2
dj
Wih
Thh T
1=1 2==1" d2
dp
8
8
T. π k
j=1 dj
This proves the existence of stationary probability. The uniqueness follows from this observation:
if ((((t)))) = n' * π., we do not have (P(()))) = (t)))) since otherwise it indicates that the
Markov chain defined by
P(R((n+1) = jn41
R(fo) =
.(, (..)) ................................ R(t..) -xj.)
(51)
= =(([(t++1) = X, [txx)) = xx))
(52)
=P(R((22 =x; [R,)) x:)
(53)
[ee c 8 + b1 if staying at current location with j i
@ Wa d,
if jumping to neighboring nodes with j *
(54)
has stationary distribution other than æ, which contradicts the assumption of a strongly connected and aperiodic graph.
We next establish the algebraic convergence as 0 < 1.
is evident that for the matrix WD given that it is column stochastic and the graph is strongly connected and aperiodic, the Ferron-Frobenius theorem Horn & Johnson (2012)[Lemm) 8.4.3., Theorem 8.4.4] confirms that the value is the unique eigenvalue of this matrix that equals its spectral radius, which is also 1. Consequently, it follows that the matrix L = I WD- has eigenvalue of 0, with all other eigenvalues possessing positive real parts. Considering the Jordan canonical form of L, denoted as L SJS-1 it is observed that J contains a block that consists solely of single while the other blocks are characterized by eigenvalues λ possessing positive real parts.
WLOG, we assume that the dimension of X R in (9), as this is consistent with handling the probability vector P(t) described in (12). We rewrite it as
D -JY
(55)
where S X(()) = (t)) E N representing a transformation of the feature space, and the trans- formed initial condition is defined as S ( = (().
the matrix is diagonalizable, then its Jordan canonical form becomes a diagonal matrix, with the diagonal elements representing the eigenvalues of L. In this scenario, the differential equation can be decoupled into a set of independent equations, each described by
D.Y - Yz(t).
(56)
40

---

Published as a conference paper at ICLR 2024
Here, Y signifies the k-th component of the vector Y. According to Diethelm (2010][h 4.3.], the solution to each differential equation the given context is represented as:
(4) (-))
(57)
where Eg(-) is the Mittag-Leffler function define as Eạ(2) = Eÿ= r(f++1) and r(-) is the gamma function. This formulation leads to two important observations:
1. For the index such that the eigenvalue X O, the solution simplifies (t) Y,(0). This corresponds to a stationary vector in the original space when transformed back to (()..
2. According to Podlubny (1999) 1.4.], for indices + j. since 2ę has a positive real part. the convergence to zero is characterized by the following order:
(a(() ((
Asymptotically, this indicates that all components Y (t), except Y, (t),,, will converge to zero at an algebraic rate. In terms of (t), this translates into convergence towards stationary vector in the eigenspace corresponding to the eigenvalue 0, while components associated with other eigenspaces diminish at an algebraic rate.
If the matrix J is not diagonal, the entries of (t) corresponding to distinct Jordan blocks in J remain uncoupled. Therefore, it suffices to consider single Jordan block corresponding non- zero eigenvalue dg- In this case, employing the Laplace transform technique becomes useful for demonstrating that the algebraic rate of convergence remains valid. We assume the Jordan block J (As), associated with ÅR, is of size It follows that for this Jordan block we have
Y. = Au Yi (2))),))),
:
2 Y ( - (t), D. Y m- (t) -
D Y (t) = (t,,,,
which can be solved from the bottom up. Beginning with the last equation, we obtain:
Ym (4) = Ym (0)()))))~>>a')) O(1-)).
Further, the differential equation for Y (()) is given by:
ml) = - ~ Y (1) Y E0)
Applying the Laplace transform and referring to (3), we obtain:
C (DDYY ( = ym) 3 m- (0)
where (s) is the Laplace transform of Y y mI)) For the right-hand side the differential equation, we have L (Ax Y mm- = A (s). A Additionally, the Laplace transform of the Mittag-Leffler function Eg (->xte) known to be 5x Podlubny (1999)(eq 1.80J. Consequently, the equation in the Laplace domain represented as:
88-1
3 1 3 y- (0) = Ym Ym Sl
+ A
Rearranging this equation to isolate (s) yields:
Yn (8) 88-18 m- (0) Y (m)
+
8 + A
As s > 0, follows that (s) (ss 8-1). Applying the same process recursively, we find that Y((s) --1 for all 1, .,... m. Invoking the Hardy-Littlewood Tauberian theorem Wikipedia (2023), we can conclude that for all indices = 1...., m, the following relationship holds:
(58) Consequently, we can deduce that, akin to the scenarios involving diagonalizable matrices, the feature components associated with other eigenspaces in condiigonaalizale cases also diminish at an algebraic rate.
Y,(t) = (
The proof now is complete.
41

---

Published as a conference paper at ICLR 2024
LIMITATIONS
Our research proposes an advanced graph diffusion framework that integrates time-fractional deriva- tives, effectively encompassing many GNNs. Nonetheless, it presents certain limitations. crucial element we have overlooked is the application of the fractional derivative in the spatial domain. In fractional diffusion equations, this implies substituting the standard second-order spatial derivative with a Riesz-Felles derivative (Gorenflo & Mainardi, 2003), thus modeling random walk with space-based long-range jumps. Incorporating such ppacefrractionll diffusion equation within GNNs could potentially alleviate issues like the bottleneck and over-squashing highlighted in (Alon & Yahav, 2021). This represents current limitation of our work and suggests compelling future research trajectory that merges both time and space fractional derivatives in GNNs.
BROADER IMPACT
The introduction of FROND holds significant potential for applications such sensor networks, transportation, and manufacturing. FROND's s ability to encapsulate long-term memory neural dy- namical processes enhance the representation of complex interconnections, improving predictive modeling and efficiency. This could lead to more responsive sensor networks, optimized routing in transportation, and improved visibility into manufacturing process networks. However, the advent of FROND and similar models may also have mixed labor implications. While these technologies might render certain repetitive tasks obsolete, potentially displacing jobs, they may also generate new opportunities focused on developing and maintaining such advanced systems. Moreover, the shift from mundane tasks could enable workers focus more on strategic and creative roles, enhancing job satisfaction and productivity. It's paramount that the deployment of FROND is done ethically, with ample support for reskilling those whose roles may be affected. This helps ensure that the broader impact of this technology is beneficial to society as whole.
42


Reviews:
[{"review_id": "vBxXO0Zs4z", "reviewer": "Reviewer_2jGF", "summary": "The paper proposed a fractional variation for graph diffusion methods and its derivative methods. The author provides numerical solutions, theoretical support, and experimental data to support their claim that fractional variation performs better than vanilla graph diffusion methods.", "questions": "N/A", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The method is novel and is a good direction for exploring graph neural diffusion methods.\n- The paper is detailed and easy to read\n- The paper has extensive comparisons between methods\n- The paper answers the question about its computation cost with detailed experiments in appendix.\n\nOverall, this paper is an updated version of a paper I've reviewed before. The authors have answered all my questions in this version. I think this paper starts from a nice idea and contains all the details required, so I would recommend acceptance.", "weaknesses": "N/A", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "9xaQSqiSYf", "reviewer": "Reviewer_Qi4C", "summary": "FROND is a method that uses concepts from fractional calculus applied to GNNs. \nThe method is based on defining the Caputo derivative and a solver that integrates the ODE. \nThe authors provide ample theory and several experiments to show the benefit of adding the fractional derivative component to GNNs.", "questions": "I left questions in the review and I am looking forward to seeing the author's response.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. The paper is well written. It was easy to follow and understand.\n\n2. The authors show how the proposed method can encapsulate existing models such as GRAND or GraphCON.\n\n3. The experiments show that adding a fractional derivative is useful.", "weaknesses": "1. Missing neural ODE literature: 'Stable Architectures for Deep Neural Networks' and 'A Proposal on Machine Learning via Dynamical Systems'.\n\n2. Missing graph ODE literature: recent papers like 'Anti-Symmetric DGN: a stable architecture for Deep Graph Networks' and 'Ordinary differential equations on graph networks'.\n\n3. In one of the main contributions, it is said that \"We provide an interpretation from the perspective of a non-Markovian graph random walk when the model feature-updating dynamics is inspired by the fractional heat diffusion process. Contrasting with the traditional Markovian random walk implicit in traditional graph neural diffusion models whose convergence to the stationary equilibrium is exponentially swift, we establish that in FROND, convergence follows an algebraic rate.\". Why is it true? if $\\beta=2$ then the process is not diffusive at all. Rather, it is oscillatory, as shown in GraphCON (Rusch et al.)\n\n4. In section 2.3, the authors should also discuss FLODE ('A Fractional Graph Laplacian Approach to Oversmoothing') which is very similar to this work and also uses fractional calculus.\n\n5. In section 3.1 the authors discuss the initial conditions of the ODE. It is not clear to me how do you initialize $\\beta$ time steps. From the text I can infer that it is the same condition as the input features. Is that was was actually done? If so, does it make sense from an ODE perspective? Have the authors tried other initialization procedures?\n\n6. The authors mention that here only $\\beta$ is only considered between 0 and 1. I wonder why. How would your model behave theoretically and practically if it larger than 1?\n\n7. I am not sure it is correct that the model can have global or 'full path' properties if $\\beta$ is smaller than 1. For example, I think it is fair to say that if $beta$ is indeed smaller than 1, then a second order process as in GraphCON cannot be represented by the model. \n\n8. The experiments indeed show that the proposed method improves compared to baselines produced by the authors, but they are quite narrow and show a partial picture of the current state of the art and existing methods. I would expect that the authors compare their work (experimentally) with other methods like FLODE, CDE, GRAND++, as well as other recent methods like ACMII-GCN++ ('Is Heterophily A Real Nightmare For Graph Neural Networks To Do Node Classification?') or DRew ('DRew: Dynamically Rewired Message Passing with Delay').\n\n9. The authors state that the proposed method can be applied to any ordinary differential equation GNN, so can the authors please also show the results when applied to other baseline methods as discussed in the paper?\n\n10. $\\beta$ is a hyperparameter. What would happen if you learn it? how will it influence the results and the stability of your solver? Is there a principles way to choosing the hyperparameter?\n\n11. I am not certain that the method is novel, as it was also shown in 'Fractional Graph Convolutional Networks (FGCN) for Semi-Supervised Learning'.\n\n12. Missing literature about skip connections in GNNs: 'Understanding convolution on graphs via energies'.\n\n13. A general point - the focus of the paper is the mitigation of oversmoothing. But, also as the authors state, there are many methods that already do it. Then my question is what is the added value of using this mechanism ? Also, another important issue is oversquashing in graphs. Can the authors discuss how and if would the proposed method can help with that issue?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Dtn17G90n4", "reviewer": "Reviewer_1hAW", "summary": "The paper extends the graph neural ODE framework by allowing fractional (Caputo) derivatives in the time variable. By leveraging on global information in the fractional derivatives, the authors prove a slow mixing theorem that prevents oversmoothing of node features. Experimentally, the authors demonstrate this fractional calculus framework over different graph ODE models and show they achieve good performances.", "questions": "Is it possible to get a clearer distinction between the approach of this work, which is modeling node feature evolution through the layers as a FDE, versus Maskey et al. 2023, which proposes using fractional graph Laplacian in the usual ODE framework. It appears that Maskey et al. 2023 approach also tackles oversmoothing via long-range dependency of the dynamics, which is the main theoretical justification of the current work as well. More specifically, are there simple examples in which one framework strictly encapsulates another? Most importantly, what is the advantage of using this framework over Maskey et al. 2023 framework? \n\nI understand fast mixing of graph random walk results of Chung 1997 and its dependence on various factors, such as eigenvalues of the adjacency/Laplacian. However, it is not immediately clear to me that the same fast rate carries over to graph neural ODE (which has some kind of skip connection across depths). Can this be explained more thoroughly?\n\nIs it possible to give a proof sketch/intuition of Theorem 2, in particular, why should we expect the slow algebraic rate? It is also interesting that the rate is tight. Following the proof in the appendix seems to suggest that this come from a deeper result by Mainardi but there is not a lot of intuition there.\n\nMinor issues:\nM1: F(s) in equation (1) is not defined (I assume it is the Laplace transform of f). The variable s is also not defined (I assume it is the variable of the transformed function). \nM2: There should also be conditions under which the Laplace transform exists (and the Laplace transform of the derivative)\n\n---\nAfter the rebuttal phase, the authors have addressed my concerns and I am raising my score to an 8.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Novelty: To the best of my knowledge, this is the first approach to directly generalize graph neural ODE to fractional derivatives and demonstrate its applicability in real-world datasets.\n\nFlexibility: The framework is general enough to be incorporated to a wide range of existing graph neural ODE in the literature, such as GRAND, GRAND++, GREAD, etc.\n\nExperiments: Empirical study conducted is extensive and results are explained comprehensively. Showing competitiveness of the new framework over existing ones in many different dimensions.", "weaknesses": "I have not studied the appendix closely so it is possible that some of these questions are addressed there. \n\nIt is not immediately clear how this current approach quantitatively/qualitatively compares to existing approaches that exploit long-range memory in the modeling process, for instance Maskey et al. 2023 (see question 1).\n\nDue to the long-range memory information, I believe it is expected that this approach is more computationally heavy than traditional neural graph ODE. It would be more complete to have a discussion of this increased cost, if there are any, as well as techniques used to overcome it. This is crucial in scaling the approach to larger datasets. \n\nI am very willing to raise my score if these issues are addressed sufficiently.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "aDYeC4S5nQ", "reviewer": "Reviewer_7EpT", "summary": "The paper introduces the Fractional-Order graph Neural Dynamical network (FROND), a novel learning framework that enhances traditional graph neural ordinary differential equation (ODE) models by integrating the time-fractional Caputo derivative. This incorporation allows FROND to capture long-term memories in feature updating due to the non-local nature of fractional calculus, addressing the limitation of Markovian updates in existing graph neural ODE models and promising improved graph representation learning.", "questions": "What is the computational complexity of FROND? What is the T chosen for each experiment? and how is the short memory principle applied?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper is well-written, providing a clear and straightforward presentation of the content, which enhances the overall readability.\n\n2. The innovative integration of time-fractional derivatives into traditional graph ODEs is a novel approach that effectively addresses key issues like non-local interactions and over-smoothing.\n\n3. The proposal is supported by theoretical motivations.\n\n4. An extensive evaluation of the framework is presented, demonstrating its effectiveness and versatility across various settings and providing substantial empirical evidence of its performance.", "weaknesses": "The correlation between beta and fractal dimemsion is not clear.  For instance, despite Pubmed having a higher fractal dimension of 2.25 compared to Airport, the optimal beta for it is set at 0.9. This observation raises curiosity about the specific conditions or types of datasets under which FROND demonstrates significant performance improvements. Clarification on this matter would greatly enhance the reader’s understanding and application of FROND in various contexts.\n\nTo further highlight the strengths of FROND and to provide clearer guidance on its optimal application scenarios, I would recommend conducting additional evaluations on datasets that necessitate long-range interactions[1]. \n\nThe content in section 3.3 offers valuable insights, and I believe it could be enriched with additional technical details and formulations related to the graph layer. This enhancement would aid readers in developing a more comprehensive and profound understanding of the model.\n\n\nDrawing a more explicit connection between fractal characteristics and FROND’s efficacy, particularly in handling tree-like data, would contribute to a more coherent narrative and justification for the framework. I kindly suggest expanding on this aspect.\n\n\n[1]Dwivedi, Vijay Prakash, et al. \"Long range graph benchmark.\" Advances in Neural Information Processing Systems 35 (2022): 22326-22340.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (spotlight)