Distributed-Order Fractional Graph Operating Network
Authors: Kai Zhao
Xuhao Li
Qiyu Kang
Feng Ji
Qinxu Ding
Yanan Zhao
Wenfei Liang
Wee Peng Tay
Venue: NeurIPS 2024 spotlight (2024)
Keywords: Graph Neural Networks
Primary area: graph_neural_networks

Abstract: We introduce the Distributed-order fRActional Graph Operating Network (DRAGON), a novel continuous Graph Neural Network (GNN) framework that incorporates distributed-order fractional calculus. 
Unlike traditional continuous GNNs that utilize integer-order or single fractional-order differential equations, DRAGON uses a learnable probability distribution over a range of real numbers for the derivative orders. 
By allowing a flexible and learnable superposition of multiple derivative orders, our framework captures complex graph feature updating dynamics beyond the reach of conventional models.
We provide a comprehensive interpretation of our framework's capability to capture intricate dynamics through the lens of a non-Markovian graph random walk with node feature updating driven by an anomalous diffusion process over the graph. 
Furthermore, to highlight the versatility of the DRAGON framework, we conduct empirical evaluations across a range of graph learning tasks. The results consistently demonstrate superior performance when compared to traditional continuous GNN models. The implementation code is available at \url{https://github.com/zknus/NeurIPS-2024-DRAGON}.

Distributed-Order Fractional Graph Operating Network
Kai Zhao' Xuhao Li Qiyu Kang", Feng Ji',
Qinxu Ding', Yanan Zhao', Wenfei Liang', Wee Peng Tay! Nanyang Technological University, Annui University, Singapore University of Social Sciences
Abstract
We introduce the Distributed-order |RActional Graph Operating Network (DRAGON), novel continuous Graph Neural Network (GNN) framework that incorporates distributed-order fractional calculus. Unlike traditional continuous GNNs that utilize integer-order or single fractional-order differential equations, DRAGON uses a learnable probability distribution over range of real numbers for the derivative orders. By allowing a flexible and learnable superposition of multiple derivative orders, our framework captures complex graph feature updating dynamics beyond the reach of conventional models. We provide a comprehensive interpretation of our framework's 's capability to capture intricate dynamics through the lens of non-Markovian graph random walk with node feature updating driven by an anomalous diffusion process over the graph. Furthermore, to highlight the versatility of the DRAGON framework, we conduct empirical evaluations across range of graph learning tasks. The results consistently demonstrate superior perfor- mance when compared to traditional continuous GNN models. The implementation code is available at https: htpp:////tthiin github.......... com www.cccom.sa.aeee.. http//knussuuee /Ner (00---224-DRRAAA
Introduction
Graph Neural Networks (GNNs) have been developed to handle graph-structured data, which prevalent in domains such social networks [1], traffic networks [2], and molecular structures [3]. The fundamental principle of GNNs is to learn representations of nodes or entire graphs that encompass both the attributes of individual nodes and the topology of their connections. This objective is accomplished through a method known as message passing or information propagation, whereby each node aggregates information from its neighbors and possibly itself, over multiple iterations or layers [4]. Recent developments in the GNN landscape have increasingly embraced the principles of continuous dynamical systems for information propagation, discussed in [5]. This trend is exemplified in works such as CGNN [6], GRAND [7), GRAND++ [8]. GraphCON [9], Beltrami [10], GREAD [11], CDE [12], and HANG [13), which employ ordinary partial differential equations (ODEs/PDEs) on graphs for feature aggregation. Within these continuous GNN models,
the differential operator de=
d" is typically constrained to integer values of a, primarily 1 or
Two directions have been proposed recently based on the aforementioned continuous GNN models to enhance their capabilities. One approach TDE-GNN [14], which proposes to learn higher integer- order temporal dependencies for continuous GNN models. The other approach is FROND [15), which incorporates graph neural Fractional-ordee Differential Equations (FDEs), extending the conventional d" integer-order derivative drer to encompass a positive real number a. This adaptation not only bolsters the model's efficacy but also enhances its adversarial robustness by varying the value of a [16].
First two authors contributed equally to this work. Corresponneecce to: Qiyu Kang <kang0080@ e entuu.ddu.sg>.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

TDE-GNN, however, is limited to utilizzig integer-order ODEs and does not account for the non- local memory effects inherent in fractional-order differential operators. These operators [17] have been developed to overcome the limitations of their traditional integer-order counterparts when modeling complex real-world dynamics. The key difference between fractional and integer operators can be grasped from microscopic random walk perspective as shown in [15, 18). For instance, traditional integer-order diffusion PDEs, which model diffusive transport in homogeneous porous media, typically ignore the waiting times between particle movements. However, these models struggle when applied to solute diffusion in heterogeneous porous media, prompting the introduction of fractional-order operators to better handle these complexities [19, In fractional scenarios, particles may remain at their current position, delaying jumps to subsequent locations with fading waiting times and leading to non Non-MMrkovian process. In contrast, traditional integer- order differential equations are typically used to model Markovian movement of particles, as the derivative df(()) dt limar-o /0+A0-1(0) captures the local rate of function changes. On the other hand, although FROND utilizes a fractional-order a and demonstrates performance improvement, its capacity for feature updating dynamics remains constrained by limited temporal dependencies with single a. Moreover, the optimized performance of FROND is achieved through extensive fine-tuning of across various graph datasets. Observations from Fig. indicate that performance can fluctuate significantly as the value of the fractional order varies from 0 to 1.
The distributed-order fractional differential operator has gained recognition in fractional calculus for its capacity to model complex dynamics that traditional differential equations with integer or single fractional orders cannot sufficiently capture [21]. Inspired by this advancement, we introduce a novel continuous GNN framework named the Distributed-order fRActional Graph Operating Network (DRAGON), which extends beyond existing frameworks like TDE-GNN and FROND. Rather than designating a single, constant & with extensive fine-tuning, DRAGON employs a learnable measure over range [a, b) for a. The foundation of our framework is the distributed- fractional differential operator [22]:
Jr D' f(t)
which can be perceived as the limiting case of E, w(a)) a weighted summation over derivatives of multiple orders with weight w(-) (we employ this more common notation 1 instead of d dt henceforth). Notably, unlike TDE-GNN, which restricts o; to integer values, DRAGON allows for a continuous range of values, significantly broadening its application scope and flexibility in modeling. This operator also addresses the limitations of the single fractional-orde operator D employed in FROND, which still has a restricted capacity to model the intricacies of feature updating dynamics. From the perspective of random walk in a diffusion process, a single D dictates that the waiting time between particle jumps follows a fixed power-law distribution xx L-a-1 for 0 < < 1. In contrast, DRAGON adopts a more flexible approach, enabling broader range of waiting times across multiple temporal scales. In this paper, we demonstrate the efficacy of the DRAGON framework in modeling more intricate non-Markovian node feature updating dynamics in graph-based data. We provide evidence that DRAGON can approximate any given waiting time probability distribution pertinent to graph random walks, thus showcasing its advanced capability in capturing complex feature dynamics.
Main contributions. Our objective is to develop a general continuous GNN framework that enhances flexibility in graph feature updating dynamics. Our key contributions are summarized as follows:
We propose a generalized continuous GNN framework that incorporates distributed-order fractional derivatives, extending previous continuous GNN models into a unified approach. Specifically, our framework treats these models as special cases with (aa) taking single positive real value for [7,8, 11, 13, 15] or multiple integer values [9, 14]. Our approach facilitates flexible and learnable node feature updating dynamics stemming from the superposition of dynamics across various derivative orders.
From a theoretical standpoint, we present the mar-Mvvkovian graph random walk with flexible waiting time for DRAGON, presuming that the feature updating dynamics adhere to a diffusion principle. This exposition elucidates the rationale behind the flexible feature updating dynamics. Through empirical assessments, we test the DRAGON-enhanced versions of several prominent continuous GNN models. Our findings consistently demonstrate their outperformance. This under-

---

scores the DRAGON framework's 's potential as an augmentation to amplify the effectiveness of range of continuous GNN models.
2 Preliminaries and Related Work
This paper focuses on developing a new GNN framework centered around distributed-order fractional dynamic processes. In this section, provide concise introduction to the key concepts in fractional calculus. Throughout the paper. we adopt certain standard assumptions to ensure problem well- posedness. For instance, the well-definedness of integrations, the existence and uniqueness of the differential equation solution [23, and the allowance for interchange between summation and limit via the monotone dominated convergence theorem [25] are all assumed.
2.1 Fractional Derivative
The single fractional-orde operator D in the distributed-order fractional operator in (1) can assume various definitions. In this study, we start off with the Marchaud-Wai Weyl fractional derivative m D", recognized for its efficacy in elucidating the fading memory phenomena [26-28], which we will discuss further in Sections 2.1.1 and 3.2.
Remark However, practical engineering implementations, the Caputo fractional derivative e is more commonly utilized [15, 17]. Due to space limitations, the introduction of Caputo's 's derivative deferred the Appendix and will be subsequently employed in Section 3.3 solve DRAGON.
The Marchaud- Weyl and Caputo definitions are equivalent under certain constraints [17.2] 29].
For any a E (0, 1), the Marchaud-Wer Weyl a-order derivative of a function f. defined over the real line, specified point t defined as [29]:
a
f(t) - 1 f((t-T) dr,
MD f(t) M(1 a)
J
71+a
where r'()) is the Gamma function. For sufficiently smooth functions, according to [29], we have
f(t At) At) f(t) lim M DD df(t) lim di At-0
0-1-
Al
is evident from (2) that the Merchaud-Weyl fractional derivative is a nonlocal operator and accounts for the past values of f within the (- range, indicative of its memory effect. In terms of probability, the related non-Markovian processes for fractional systems are characterized by state evolution that depends
just the current state, but also historical states [18]. As - 1 in (3), the operator reverts to the traditional first-order derivative, representing the local change rate of the function with respect time.
2.1.1 Non-Markovian Random Walk Interpretation
We elucidate fractional-order derivatives by linking them to one-dimensional heat diffusion and memory-decaying Mon-Markovian random walks [28]. Assuming random walker moves along an axis with infinitesimal intervals of space Ar > 0 and time AT > O, the walker moves a distance of Ax from the current point x in either direction with equal probability and waits at each location for random period of time, a positive integer multiple of Ar. This introduces randomness in the waiting times between steps. We aim to compute u(x, t), the probability of the walker arriving position at time f. The waiting time distribution, tha (n), is given by power-law function n -(1+a) with da > 0 chosen to ensure Ene1 (a) = The law of total probability is expressed as: [
u(x; Σ u( u( Ar,t ndr) + u(a Az,t mdr) (n..
Here, the terms within brackets denote the probability of arriving at x from either neighboring points, Ar or x + Ar. each with probability 1/2. The sum n accounts for the possibility that the walker could have remained stationary for an extended period nAr with waiting time probability do (n). After subtracting Σ - Via (n))n(x( ndT) from both sides and rearranging, we obtain:
u(x, t) - u(x, 1 AT) (AT) 2d, (ar)a Σ
(Ax)²
8zu(2, t Bouz,,tnArrúú.lC ()))
(4)
(+a
3

---

where the aecond-order incremental quotient is defined as:
02u(1, u(x Az, t) + u(x Ar, - 2u(xxt)
(Az)
In the limit as Az, AT + and assuming that (Ax)² du(Ar)
+ ka (alal for a positive hia [28], we
obtain the time-fractional diffusion equation:
ki Myz:
M DD
(5)
2
where the summations the left-hand side of (4) converge to the integration (2). As (5) reverts to the standard heat diffusion equation:
Ōu(x, t) A
Ôt 2 Mgz-
(6)
Consequently, the aforementioned non-Markovian random walk with fading memory simplifies to the Markovian random walk, thereby eliminating the memory effects.
2.2 Integer-Order Continuous GNN Models
We denote undirected graph
(V, W). where V is the set of | N nodes and X ([x]"",. ..... [xx]")* E Nxd consists of rows X; E R as node feature vectors. The N x N adjacency matrix = (W1z) has elements Wij indicating the edge weight between i-th and j-th nodes with Wa; = W jä- In the subsequent GNNs inspired by dynamic processes, we let (() = (t)] ...., (t))") E IR Nxd be the features at time with X(0) = X serving as the initial condition. The time here acts analog to the layer index [7,30]. Typically, these dynamics can be described by:
dx(t)
F(W, F(W,((()))
dt
The function F is specifically tailored for graph dynamics as illustrated Appendix F. For instance, in the GRAND model, F is defined follows:
GRAND [7]: Drawing from the standard heat diffusion equation, GRAND formulates the following feature updating dynamics:
dx(() (A(x(t)) I).()).
(8)
df
where A A(X(() is learnable attention or fixed normalized matrix, and I is identity matrix.
2.3 Fractional-Order Continuous GNN Models
Recently, the paper [15] introduces FROND, extending traditional integer-order graph neural differen- tial equations such as (8), (40) and (42) to fractional-order equations. The framework is formalized as
D X F(W, X()),,
0,
where fE represents the graph dynamics. Further, the study [16] explores the robustness of FROND, demonstrating its ability to enhance the resilience of integer-order continuous GNNs under perturbations. This underscores the potential applications of FROND in various domains.
2.4 Motivation: Advanced Dynamics Modeling Capability
To intuitively understand the versatility and efficacy of the DRAGON framework in learning dynamics, consider three classical stress- stres--strain constitutive mod- els for viscoelastic solids: the single-order Maxwell model [31], the multi-order Zener model [32], and the distributed-orre Kelvin- Kelvin-Voigt model [33]. Using the FROND and DRAGON franeworks, we develop Neu- ral (NN) methods to predict future states based on current observations.
Table 1: Comparison of MSE for the Maxwell, Zener, and Kelvin- foint mod- els using FROND-NN and DRAGON-NN frameworks.
Model
FROND-NN
DRAGON-NN
Maxwell [31]
2.0 10
5.6 x 10-5
Zener (32)
3.6 x 10
3.5 x 10-8
Kelvin- Voigt [33]
3.3 x 10
1.4 x 10

---

The detailed descriptions and implementation specifics can be found in Appendix G.1. The results presented in Table demonstrate that the DRAGON framework excels in fitting not only the multi- order model but also in capturing the dynamics of single-order and distributed-order models. We observe that the ARAGON framework achieves a Mean Squared Error (MSE) that ten times smaller than that of the FROND method across all three models. This highlights the DRAGON framework's s exceptional ability to effectively learn and adapt to a diverse range of dynamics, surpassing the capabilities of FROND.
3 DRAGON Framework
In this section, we introduce our general DRAGON framework for GNNs, with random walk interpretation that elucidates the underlying mechanics when a specific diffusion-inspired system is utilized. Subsequently, we discuss numerical techniques for solving DRAGON. The versatility of our framework is highlighted by its capacity to encapsulate a broad spectrum of existing continuous GNN architectures, while simultaneously nurturing the development of more flexible continuous GNN designs within the research community in the future.
3.1 Framework
DRAGON generalizes the current integer-order and fractional-order continuous GNNs as it uses a learnable probability distribution over a range of real numbers for the fractional derivative orders. Consider a graph G = (V, W) composed of |V] = N nodes with W being adjacency matrix as defined in Section 2.2. Similar to the approach used integer-order continuous GNN models [5. 15] as presented in Section 2.2, we apply a preliminary learnable encoder function V - R that maps each node to a feature vector. After stacking all these feature vectors, we obtain X e R.N xd Employing the distributed-order fractional derivative outlined in (I), the feature dynamics in DRAGON are characterized by the following graph dynamic equation:
[ D*X(() = F(W, X(t)),
(10)
where [a, b) denotes the range of the order a, μ is a learnable measure of a, and F is dynamic operator on the graph as illustrated in Appendix F.
Remark 2. In practical engineering settings, the Caputo fractional derivative, represented by cD is commonly used [15, 17]. When leveraging the Caputo definition for the fractional derivative, as detailed in Section 3.3. the initial condition for (10) is given by (0) = X, where (()(()) denotes the n-th order derivative at t = O, encompassing the initial node features for all integers n E IN [0 [b]] [23]. Here, [-1 is the ceiling function, and this setup ensures a unique solution [23]. For instance, when [a, b) = [0, 1), we define the initial condition (() X.
This framework generalizes prior continuous GNN models, encompassing them special instances. Specifically, with p(a)= = 5(a= - 1), where is the Dirac delta function, (10) simplifies to a local first-order differential equation like [7,8, 10-13]. When [a, bj = [0. we may obtain distributed- order fractional wave propagation GNN model [21]. which generalizes the second-order GraphCON model (40). When μ(a) = ō(a a.) for a, IR+ (10) reduces to the FROND framework (9). Additionally, when j adopts discrete distribution over multiple integers, the model corresponds to TDE-GNN [14].
Following previous works, we set an integration time parameter T' to obtain X(T). The final node embeddings, employed for subsequent downstream tasks, can be decoded ((X(T)), where symbolizes learnable decoder function.
3.2 Non-Markovian Graph Random Walk with Flexible Memory
In this subsection, we provide a non-Markovian graph random walk interpretation for DRAGON under specific anomalous diffusion setting, where the dynamic operator J(W, X())) in (10) is set as (t)) - (t) in (8) with a fixed constant matrix A. More specifically, we obtain the following linear distributed-order FDE:
[ MD X(t) du((a) LX((),
(II)
5

---

where we set A = WD a WD -I is the random walk Laplacian. Here, D is a diagonal matrix with D d, the degree of node For clarity, without loss of generality, similar to the approach in Section 2.1.1, we interpret (() as N-dimensional probability or concentration vector P(()) over the graph nodes V at time The Marchaud-Wey m D employed (11) helps expedite the exposition of the subsequent random walk, drawing analogy from the one-dimensional random walk discussed in Section 2.1.1.
For every individual value Oa E (0, 1), we consider a random walker navigating over graph G with an infinitesimal interval of time År > We assume that there no self-loop in the graph topology. The dynamics of the random walk are characterized as follows:
1. The walker is expected to wait at the current location for random period of time. The distribution of waiting times, Va (nz). is given by power-law function da, n -1++a) with de. chosen to ensure Σn=1 Va. (n) 1.
2.
Upon deciding to make a jump, the walker can either move from the current node to a neighboring node with a probability of (Ar)"oda, ( Way if + j. Alternativel, with probability of 1 - (AT) (Ar)"oda. dad,ITa(aaa (T(-))))) it will remain at the current node i.
waiting time is wude
tovioit provioue Sudie
waiting time 2Ar
wailing time Ar
time discretization
24
Ar
jump from ueighbos
Figure Visualization of the Non-Markovian Graph Random Walk. The diagram illustrates the walker's decision-making process during the walk. After waiting for random duration nA, the walker may either remain on the current node or proceed to a neighborhood node. This reflects the flexible, memory-influenced dynamics of the walker's movement.
We denote (t; a.), the probability of the walker being at node time with a specific order oa and p(a)) = 8(a a.).. The law of total probability is expressed as:
Wij di
Py(t;ao) Σ
Σ
Pt(t nAT; a.) ranalralaolll
(12)
n=1
iCV #2
+ (o-l--AT;;are (Ar) [ (T(-a)))) taa (n).
In this equation, the summation over n accounts for the possibility that the walker may have remained stationary for period of nAT, with waiting time probability of th. (n). Fig. provides visualiza- tion of the non-Markovian graph random walk. For more explanation of the non-Markovian random
walker graphs, please refer to Appendix E. From (12), we derive Theorem 1.
Theorem 1. Given p(a)) 5((x ata) where an (0,1) and AT the probability vector whose j-th element is (1; a.), solves (11). That is to say, we have
establish that PII:; oo)).
S:' M (t; du(()) = LTP(t:
(13)

---

In Theorem 1, we present the graph random walk interpretation for the fractional anomalous diffusion equation (11) under the condition that μ =(aa a.). This condition represents single-term fractional scenario similar to FROND. At its core, this type of random walk is non- Markovian, underscoring the importance of the entire walk history.
Remark 3.
From the discussion above. for a specific Coo the waiting time is steered by the power-law distribution x n n=(+++1) Moreover, the distributed-order fractional operator can be interpreted as a flexible superposition of the dynamics behaviors embodied by individual fractional-order operators. This generalization reframes the interpretation of graph random walk and enables more nuanced dynamics that accommodate diverse waiting times. Although it is feasible to formulate a random walk interpre- tation where the waiting time is linked to the measure and converges to the solution of (11), this approach relies on the intricate stopping time technique [34][Sec 7.5] and may sacrifice flexibility in waiting time insights. Instead, we propose more modest conclusion, demonstrating that a weighted sum of Vo, (n) can approximate any waiting time, highlighting the capability of our framework in comparison to FROND.
Theorem 2. Let Co(N) be the space of functions on the natural numbers N vanishing al oo, i.e., f (NI if and only limn-+oo f(n) O. Assume the sequence (am))== is strict increasing in [0, 11, then the span of {Tam. m 2 1} is dense in Co(N) in the sense of uniform convergence.
Remark 4. Theorem demonstrates the DRAGON framework's 's ability approximate any wailing time distribution for graph random walkers, offering flexibility in modeling feature updating dynamics with varying extents of memory incorporation. This highlights the advantage of using DRAGON for deploying learnable and flexible feature updating dynamics. In contrast, FROND is confined to fixed waiting time distribution, limiting its adaptability in modeling feature updating over time.
3.3 Solving DRAGON
Previous continuous GNNs have leveraged neural ODE solvers [30] when μr 6(a - 1). For example, in the explicit Euler scheme, neural ODEs are effectively reduced residual networks with shared hidden layers [30]. Addressing the challenge of solving the distributed-orre DRAGON, the standard approach involves discretizing it into multi-term FDE. This is achieved by using a quadrature formula to approximate the integral term [21,23]. As articulated in Sections 2.1 and 1. we follow the convention in the fractional calculus literature for real-world applications and employ the Caputo definition c in this section. This choice is intuitive, as it seamlessly incorporates initial conditions into the problem as previously discussed under (10). The initial step is to approximate (10) as follows:
FDE (10) given by
Σ De F(W,
(14)
where aj E [a, | 1....... are distinct interpolation points and wy are weights associated with the measure p. Reflecting the learnable nature of μa, Wy is directly set to be a learnable parameter in our implementation.
The next step is to solve the multi-term FDE presented in (14). According to the approach outlined in [17, Theorem 8. the multi-term FDE can be transformed into system of single-order equations DDT, where 7 1/M and M is the least common multiple of the denominators of axo a...... On when these coefficients are rational numbers. The classical fractional ddlass-ashforth-MMulldd method can then be applied to solve the resulting system of single-order equations [15, (55. This method is generalization of the Euler scheme for ODEs to fractional scenarios (see Appendix C.3 for a detailed explanation).
An alternative approach involves directly approximating the fractional derivative operators as demon- strated in [36]. This discretization method can then be used to derive iterative methods for solving the multi-term FDE given in (14). Detailed procedures for this method are provided in Appendix C.4. Additionally, the approximation error analysis of the numerical solvers discussed in Appendix D.
3.4 DRAGON GNNs
In Section 2.2 and Appendix F, several continuous GNNs, such as (8), (40) and (42), which employ integer-order derivatives, are introduced. We now extend these dynamical systems to operate under
7

---

our proposed DRAGON framework, which generalizes the scenarios to involve distributed-orde fractional derivatives. More specifically, we present the following GNNs, which will be utilized in Section 4 to show the advantages of our framework over various graph benchmarks.
1. D-GRAND: By extending (8), we get
L. D dplaa (A(X(()) - ())(())
(15)
2. D-GraphCON: By extending (40), we get
L.
= ((F - (t)(
(16)
3. D-CDE: By extending (42), we get
[ D D'X(()du(a)) (A(X(1) ()x(() - div(V(()) X())),
(17)
where div(V(()) (()) is given in (43) and (44).
Depending the method used to compute the matrix A in (15), the D-GRAND model can be categorized into two versions: linear (D-GRAND-I) and non-linear (D-GRAND-n)), Similarly, based on the computation of F, in (16), the D- D-GraphCON model also has two versions: linear (D- GraphCON-I) and non-linear (D-GraphCON-nl) Detailed explanations are provided in Appendix F.1.
Experiments
Our approach aims to enhance the capabilities of continuous GNN models by flexibly combining graph dynamics across different derivative orders. To achieve this, we have integrated DRAGON into several existing continuous GNN models and assessed their performance. Specifically, we conduct experiments on our proposed D-GRAND (15), D-GraphCON (16), and D-CDE (17) in this section, as well as D-GREAD and D-GRAND++ in Appendix I.3 and Appendix 1.4.
4.1 Implementation Details
In our approach, we employ fully connected (FC) layer as the encoder, Po V - IRd to determine the initial values for DRAGON. Subsequently, another FC layer ( serves as the decoder, transforming the output of DRAGON for downstream tasks. Most existing con- tinuous GNNs are first-order or can be transformed into first-order representations of certain dynamic processes across graphs [9]. The FROND framework also restricts the fractional order to the range (0, 1). maintaining iden- tical initial conditions to those utilized in the original models. Given these considerations, we mainly restrict aj values between [0.1] in our implementation, while also balancing computational costs. The parameter a, is selected to evenly divide the entire range, aiming to comprehensively cover values between [0, 1]. Typically, set the number of a; in (14) to 10. We also explore the empirical results when oy exceeds 1, as shown in Appendix I.5. For sensitivity analysis of the number and value of aj. we refer the readers to Appendix 1.6. Details on the datasets used can be found in Appendix G.2.
Table 2: Numerical results for various meth- ods on LRGB tests. Method
Peptides-fune Test AP
Peptides-Strue Test MAE J
GCN [37)
0.5930=0.000 0.0033
0.349610.0013
GCNII (38)
0.5543:.00078
0.3471110.000
GINE [39]
0.549810.0079
0.34471000005
GatedGCN (40)
0.586410.00777
0.3420:1.00013
Transformer-LapPPE [41]
0.6326+-n.01
0.2529-40.0016
SAN+LapPE [42]
0.638410.0121
0.2683110.043
SAN+RWSE [43]
0.6039
0.2545+0000122
GCN+DRew [44]
0.699620.0076
0.2781
07781+0.0002
PathNN [45]
0.6816:0000026
0.2545
DRGNN [46]
0.6586=0.0042
0.2495++0.0015
GRAND-I
0.6962:0.0015
0.286710.0000
F-GRAND-I
0.7126=00.002
0.2677+40.0014
D-GRAND-1
0.7571=0.0014
0.2461:10.0014
4.2 Long Range Graph Benchmark
As illustrated in Remark 4, the DRAGON framework exhibits a distinctive intrinsic property: its ability to capture flexible memory effects, which is crucial for modeling long-range dependencies
8

---

in graph data [44]. To empirically validate this capability, we conduct experiments using the Long- Range Graph Benchmark (LRGB) [47]. Specifically, we focus on the Peptides molecular graphs dataset, performing graph classification on the Peptides- func dataset and graph regression based on the 3D structure of peptides in the Peptides-struct dataset. The performance metrics used are Average Precision (AP) for classification and Mean Absolute Error (MAE) for regression tasks. From Table 2, it is evident that the DRAGON framework outperforms the other methods on these two long-range graph datasets, even when compared to state-of-the-art (SOTA) techniques. Notably, DRAGON achieves an improvement of approximately 4-6% over traditional continuous GNNs like GRAND-I and F-GRAND-I. This demonstrates DRAGON's capability to effectively capture long-range dependencies in graph data.
4.3 Node Classification
4.3.1 Homophilic Graph Datasets
In our evaluation on homophilic datasets, we leverage a diverse set of datasets including citation networks (Cora [48], Citeseer [49]. Pubmed [50]). tree-structured datasets (Disease and Airport [51], as well as coauthor and co-purchasing graphs (CoauthorCS [52]. Computer and Photo [53]). For the Disease and Airport datasets, we follow the data partitioning and preprocessing procedures as described in [51]. For all other datasets, we adopt random splits for the largest connected component (LCC), line with the approach detailed in [7].
Table 3: Node classification results(%) for random train-val-test splits. The best result of each continuous GNN family is highlighted in red.
Methed
Cora
Citeseer
Pubmed
CoanthorCS Computer
Photo
CoauthorPyy
Airport
Disease
GCN [37]
1.3
1,9 7.88
91.1-0.5
82.6-14
91.24-12
92.81100
81.610066
69.8.10.5
GAT [54]
11.
1..4+19
78.7+2.3
90.5-0.6
78.0=190
85.7+20.3
92.5+099
81 6+:400
70.4+0.5
HGCN [51]
78.7+100
65.8=2.0
16.4=08
90.6=0.3
80.64188
88.2+14
908811.5
85.4+0.77
89.9+111
GIL [55]
82.1+1.1
71.1+1.2
77.800.6
89.4=1.5
89.6+1..
91.54177
90.8+00.5
GRAND-I
83.6+100
73.4-00.5
7.78+1.7
92.9-0.4
83.7+11.2
92.3-10.9
93.5++0.99
80.54-9.6
74.5+-3.4
F-GRAND-1
84811.1
74.0.1.5
84411.
92.810.6
94.5100.4
98.1100.2
92.41399
D-GRAND-1
8,31+1133 7.
79.61-2.6
93.2 0.3
87.3
93.1
94.6-10.2
98.5
93.21-2.5
GRAND ml
82.3.11.6
70.9.1.0
77.51188
92.4- 0.3
82.4.1.11
92.410.88
91.4.1.3
90.911.6
81.0.66.7
I-GRAND-al
83.211.1 74.7219
79.220.77
92.9=04
84.11009
93.11:0.9
93.91005
96.110077
85.5122.5
G.GRAND ml
83
74.8 1.6
7.6
93.1 0.3
87.1.1.0
93.4
94.3.10.6
97.7.10.4
89.3
GraphCON-1
81911.7
72.92211
78.822.
92.3-033
84.9:0055
90.81188
93.9110.4
68.61:2.1
87.54:4.1
F-GraphCON-1
84.6214
75.3=111
1..3+13
92.8=04
86244088
93.3+110
94.12:0.5
97.3400.5
92.142.8
D-GraphCON-J
84.6:11.3
74.42114
80.7:116
92.9-0.3
86.9:1100
93.7100.
94,3:10.5
98.3:10022
93.3:12.11
GmphCON-al
83.2+1144
78.241.8
1.
8..7-099
79.24-1.1
85.5+-2.3
93.14-0.3
74.14227
65.74-5.9
F-GraphCON-nl
83.9+1.2
1.
1.
90.4=0.6
83.6+2.22
94.1--0.7
93 93.0--0.6
97.3+008
86.94-4.0
D-GraphCON-al
84.2:11.2
74.0=2.1
79.5=1.1
92.0=0.2
87.1110
93.8:00.8
94.0:10.4
98.31003
91.4116
Table 4: Node classification results(%). The best and the second result for each criterion are highlighted in red and blue, respectively.
Method Audi
Roman empire
Wiki cooe
Minesweeper 0.01
Questions 0.02
Workers 0.09
Amazon ratings
-0.05
-0.03
0.14
ResNet (56)
65.71=0.44
89.3620.71
50.95::1.12
70.10=:.775
73.08++1.29
45.70=20.69
H2GCN [57]
68.09=0.29
89.24:1.32
89.95:10.38
66.66:1
81.76:10.68
41.36:20.47
CPGNN [58]
63.78=0050
8..84-0066
71.27+11.14
2..09-+2.3
7..44+1.80
4..36+00.35
GPR-GIN [59]
73.37=0.68
91. 91.9010.788
81.79:10.98
73.41::1.24
10.59::1.15
43.90=10.48
GloGNN [60]
63.85=0049
0..4
1.35
6..15:t1.92
73.90:10.95
37.28:1066
FAGCN [61]
70.53=0.99
91 0.37
89.69:10.0
1.50
81.87400.4
46.32-2250
GBK GBK-GNN [62]
75.87=0.43
97.81 0.32
83.56:10.84
12.98:11.05
78.06:10.91
43.47:10.51
ACM-GCN [63]
6..3
1.04
90.47:1.557
OOM
78.2540.78
38.51#3.38
GRAND
71.60=0.58
92.03-20446
76.67:1..8
70.67+:1.2
75.33-20.84
45.05=0065
GraphBel [10]
69.47=40.37
90.30+00.50
76.51++1.03
70.79+10.99
73.02-+0.92
03.63-+0.42
Diag-NSD [64]
0..50=00.67
92.06:1..40
89,5911061
15
79.81110.99
07.96:10.20
ACMP [65]
01.27
92.68-00.37
76.15+1112
1..3
75.03-1-.92
0..76=10.52
TDE-GNN (14)
64.29-0..58
84.95
61.15:2.24
6..94:11..9
75.13-10.81
40.33:11.37
CDE [12)
91.64=0.28
97.9910.38
95 03.50:55.23
75.17=10.99
10.70:11.04
47.63-10.43
F-CDE [15]
93.06=0.555
98.73-00668
96.04-10.25
75.17-10.99 99
82.68-10.86
49.01=0.56
D-CDE
93.87=0.41
98.58=0.12
96.47-11.89
75.53-10.98
83,02-10.86
49.43::1.26

---

4.3.2 Heterophilic Graph Datasets
For evaluating performance on heterophilic datasets, we utilize six datasets introduced in [66], with details provided in Appendix G.2. As highlighted in [66], these datasets are characterized by lower adjusted homophily Aadj. indicating a higher degree of heterophily. In our experimental setup with these heterophilic datasets, we follow the data splitting strategy described in [66], dividing the data into 50% for training. 25% for validation, and 25% for testing.
4.3.3 Performance of DRAGON framework
As shown in Table 3, for homophilic datasets such as citation networks, coauthor networks, and co- purchasing networks, our DRAGON framework enhances the performance of continuous backbones like GRAND and GraphCON. This demonstrates the ability of our DRAGON framework to seamlessly integrate with existing continuous GNNs and improve their performance. Notably, on tree-structured datasets, our DRAGON framework significantly boosts the performance of both GRAND and GraphCON. In particular, on the Airport dataset, our DRAGON framework excels, achieving a 7% performance increase compared to the GIL model specifically designed for this type of tree-like dataset. Compared to FROND, our DRAGON framework shows improvements most datasets. The results of the graph node classification on heterophilic datasets are presented Table 4. As indicated in Table 4, the proposed D-CDE model with our DRAGON framework improves the performance of the original CDE and F-CDE models on five out of the six datasets. This underscores the ability of DRAGON to capture flexible memory effects as proved in Theorem 2, highlighting its enhanced capability in modeling complex feature updating dynamics.
4.4 Model Complexity
For the
method (25), the numerical solution is computed iteratively for E = T/h time steps, where h represents the discretization size and T the integration time. This process involves repeated computation of J (W, X,) for each iteration. By storing intermediate we can express the total computational time complexity
function evaluation values {F
across the process as - + O())) where (k) indicates the computational overhead from summing and weighting the terms each step. Here, C represents the complexity of computing F. This yields total cost of O EEC + E²). If fast algorithm for the convolution computations is available, we typically require OEE log E) for the convolution [67], resulting in O EC + E log E). If the cost of weighted summing is minimal, the complexity is reduced to O (EC). For the Grünwald- Letnikov method (32), the computational complexity is the same as that of the method (25).
The term C' denotes the computational complexity of the function F. For instance, setting F to the GRAND model results in C' = [E|d, where E|| represents the edge set size and d the dimensionality of the features [7]. Alternatively, using the GREAD model results in C' = O([[) + [E2))d + |E (daax) where |Ez| accounts for the number of two-hop edges, and dmax is the maximum degree among nodes [1]]. More details of the computation cost can be found in Appendix H.
5 Conclusion
We introduce the DRAGON framework, which incorporates distributed-order fractional derivatives into continuous GNNs. DRAGON advances the field by employing a learnable distribution of fractional derivative orders, surpassing the constraints of existing continuous GNN models. This approach eliminates the need for fine-tuning the fractional order, as required in FROND, and enriches the dynamics and representational capacity of existing continuous GNN models. We also provide a flexible random walk interpretation. Through rigorous empirical testing, DRAGON has demonstrated not only its adaptability but also its consistent outperformance compared to other continuous GNN models. Consequently, DRAGON establishes itself as a powerful framework for advancing graph- related tasks.
10

---

Acknowledgments and Disclosure of Funding
This research is supported by the National Research Foundation, Singapore and Infocomm Media Development Authority under its Future Communications Research and Development Programme. Xuhao Li is supported by National Natural Science Foundation of China (Grant No. 12301491). The computational work for this article was partially performed on resources of the National Supercom- puting Centre, Singapore https://www.nscc.sgg To improve the readability, parts of this paper have been grammatically revised using ChatGPT [68].
References
[1] S. Wu, F. Sun, W. Zhang. X. Xie, and B. Cui, "Graph neural networks in recommender systems: a survey," ACM Computing Surveys, vol. 55, no. 5, pp. 1-37, 2022.
[2] W. Jiang and J. Luo, "Graph neural network for traffic forecasting: A survey," Expert Systems with Applications, vol. 207, p. 117921, 2022.
[3] Y. Wang, J. Wang, Z. Cao, and A. Barati Farimani, "Molecular contrastive learning of repre- sentations via graph neural networks," Nature Machine Intelligence, vol. 4, no. 3, pp. 279-287, 2022.
[4] Feng, Y. Chen, F. Li, A. Sarkar, and M. Zhang, "How powerful are k-hop message passing graph neural networks," Advances in Neural Information Processing Systems, vol. 35, pp. 2022.
[5] A. Han, D. Shi, L. Lin, and Gao, "From continuous dynamics to graph neural networks: Neural diffusion and beyond," arXiv preprint arXiv:2310.00121, 2023.
[6] L.-P. Xhonneux, M. Qu, and Tang. "Continuous graph neural networks," in Proc. International Conference Machine Learning, 2020, pp. 10 404-11441.
[7] B. Chamberlain, J. Rowbottom, M. Gorinova, Bronstein, Webb, and Rossi, "Grand: Graph neural diffusion," in Proc. lnt. Conf. Mach. Learn., 2021, pp. 1407-1418.
[8] M. Thorpe, H. Xia, T. Nguyen, T. Strohmer, A. Bertozzi, S. Osher, and B. Wang, "Grand++: Graph neural diffusion with source term," in Proc. International Conference Learning Repre- sentations, 2022.
[9] T. K. Rusch, B. Chamberlain, J. Rowbottom, S. Mishra, and M. Bronstein, "Graph-coupled oscillator networks," " in Proc. International Conference Machine Learning, 2022.
[10] Y. Song. Q. Kang, S. Wang, K. Zhao, and W. P. Tay, "On the robustness of graph neural diffusion to topology perturbations," in Advances Neural Information Processing Systems, 2022.
[11] J. Choi, S. Hong, N. Park, and S.-B. Cho, "Gread: Graph neural reaction-diffusitn equations."
in Proc. International Conference Machine Learning, 2023.
[12] K. Zhao, Q. Kang, Y. Song, R. She, S. Wang, and W. P. Tay, "Graph neural convaction-diffusion with heterophily," in Proc. International Joint Conference on Artificial Intelligence, 2023.
[13] K. Zhao, Q. Kang, Y. Song, R. She, S. Wang, and W. Tay, "Adversarial robustness in graph neural networks: A hamiltonian approach," in Advances Neural Information Processing Systems, 2023.
[14] M. Eliasof, E. Haber, E. Treister, and C.-B. B. Schönlieb, "On the temporal domain of differential equation inspired graph neural networks," in International Conference on Artificial Intelligence and Statistics. PMLR, 2024, pp. 1792-1800.
[15] Q. Kang, K. Zhao, Q. Ding, F. Ji, X. Li, W. Liang, Y. Song, and W. Tay, "Unleashing the potential of fractional calculus in graph neural networks with FROND," in Proc. International Conference on Learning Representations, Vienna, Austria, 2024.
[16] Q. Kang, K. Zhao. Y. Song. Y. Xie. Y. Zhao, Wang, R. She, and W. P. Tay, "Coupling graph neural networks with fractional order continuous dynamics: A robustness study," in Proc. AAAI Conference on Artificial Intelligence, Vancouver, Canada, Feb. 2024.
[17] K. Diethelm, The analysis of fractional differential equations: an orpiieation-oriented exposi- tion using differential operators of Caputo type. Lect. Notes Math., 2010, vol. 2004.
11

---

[18] R. Gorenflo and F. Mainardi, "Fractional diffusion processes: probability distributions and continuous time random walk," " in Process. Long--Ragge Correlations: Theory Appl., 2003, pp. 148-166.
[19] C. Jonescu, A. Lopes, D. Copor, J. T. Machado, and J. H. Bates, "The role of fractional calculus in modeling biological phenomena: review," Communications in Nonlinear Science and Numerical Simulation, vol. 51, pp. 141-159, 2017.
[20] D. Krapf, "Mechanisms underlying anomalous diffusion in the plasma membrane," Current Topics Membranes, vol. 75. pp. 167- 2015.
[21] W. Ding, Patnaik, S. Sidhardh, and Semperlotti, "Applications of distributed-order frac- tional operators: A review," Entropy, vol. 23, no. 1, 110, 2021.
[22] M. Caputo, "Mean fractional-orde--eervvvi University of Ferrara, vol. 41, no. 1, pp. 73-84, 1995.
differential equations and filters," Annals of the
[23] K. Diethelm and N. J. Ford, "Numerical analysis for distributed-orde differential equations." J. Comput. Appl. Math., vol. 225, no. 1, pp. 96-104, 2009.
[24] K. Diethelm and N. Ford, "Analysis of fractional differential equations." Journal of Mathemati- cal Analysis and Applications, vol. 265, no. 2. pp. 229-248, 2002.
[25] Billingsley, Convergence of probability measures. John Wiley Sons, 2013.
[26] S. G. Samko, "Fractional integrals and derivatives," Theory Appl., 1993.
[27] A. Bernardis, F. J. Martín- Reyes, P. R. Stinga, and L. Torrea, "Maximum principles, extension problem and inversion for nonlocal one-sided equations," J. Differ. Equ., vol. 260, no. 7, pp. 6333-6362, 2016.
[28] R. Stinga, "Fractional derivatives: Fourier, elephants, memory effects, viscoelastic materials and anomalous diffusions," arXiv preprint 1rXiv:2212.02279, 2022.
[29] F. Ferrari, "Weyl and marchaud derivatives: A forgotten history," Mathematics, vol. 6, no. 1, p. 6, 2018.
[30] R. T. Chen, Y. Rubanova, Bettencourt, and D. Duvenaud, "Neural ordinary differential equations." arXiv preprint arXiv: /80600366, 2018
[31] Y. Rossikhin and M. Shitikova, "A new method for solving dynamic problems of fractional derivative viscoelasticityy" Int. J. Eng. Sci., vol. 39, pp. 149-176, 2001.
[32] T. Atanackovic, Konjik, L. Oparnica, and Zorica, "Thermodynamical restrictions and wave propagation for a class of fractional order viscoelastic rods," Abstr. Appl. Anal., vol. 2011, 2011.
[33] B. Stankovic and Atanackovie, "Dynamics of a rod made of generalized kelvin-voigt visco- elastic material," J. Math. Anal. Appl., vol. 268, pp. 550-563, 2002.
[34] M. M. Meerschaert and A. Sikorskii, Stochastic models for fractional calculus. Walter de Gruyter GmbH & Co KG, 2019, vol. 43.
[35] K. Diethelm, N. J. Ford, and A. D. Freed, "Detailed error analysis for a fractional adams method," Numer. Algorithms,, vol. 36, pp. 31-52, 2004.
[36] D. Baleanu, K. Diethelm, E. Scalas, and J. Trujillo, Fractional calculus: models and numerical methods. World Scientific, 2012, vol.
[37] N. Kipf and M. Welling, "Spri-superviiee classification with graph convolutional networks." in Proc. Ini. Conf. Learn. Representutionss 2017.
[38] M. Chen, Z. Wei, Z. Huang, B. Ding, and Y. Li, "Simple and deep graph convolutional networks," Proc. Inf. Conf. Mach. Leurn, 2020, pp. 1725-1735.
[39] W. Hu, B. Liu, J. Gomes, M. Zitnik, Liang, V. Pande, and Leskovec, "Strategies for pre-training graph neural networks," 2020.
[40] X. Bresson and T. Laurent, "Residual gated graph convnets," 2018.
[41] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, F.. Kaiser, and
I. Polosukhin, "Attention all you need," Advances in neural information processing systems, vol. 30, 2017.
12

---

[42] D. Kreuzer, D. Beaini, W. Hamilton, V. Létourneau, and P. Tossou, "Rethinking graph trans- formers with spectral attention." " Advvances in Neural Information Processing Systems, vol. 34,
pp. 21 618-21 21 629, 2021.
[43] V. P. Dwivedi, A. T. Luu, T. Laurent, Y. Bengio, and X. Bresson, "Graph neural networks with learnable structural and positional representations," 2022.
[44] B. Gutteridge, X. Dong, M. M. Bronstein, and F. Di Giovanni, "Drew: Dynamically rewired message passing with delay," in International Conference on Machine Learning. PMLR, 2023,
pp. 12252-12267 267.
[45] G. Michel, G. Nikolentzos, J. F. Lutzeyer, and M. Vazirgiannis, "Path neural networks: Expres- sive and accurate graph neural networks," " in International Conference on Machine Learning. PMLR, 2023, pp. 24 747-244775
[46] J. M. Baker, Q. Wang, M. Berzins, T. Strohmer, and B. Wang, "Monotone operator theory- inspired message passing for learning long-range interaction on graphs," in International Conference on Artificial Intelligence and Statistics. PMLR, 2024, pp. 2233-2241.
[47] V. P. Dwivedi, I.. Rampášek, M. Galkin, A. Parviz, G. Wolf, A. Luu, and D. Beaini, "Long range graph benchmark," 2023.
[48] A. McCallum, K. Nigam, J. D. M. Rennie, and K. Seymore, "Automating the construction of internet portals with machine learning," Inf. Retrieval, vol. 3, pp. 127-163, 2004.
[49] P. Sen, G. Namata, M. Bilgic, L. Getoor, B. Galligher, and T. Eliassi-Rad, "Collective classifica- tion in network data," Al Magazine, vol. 29, no. 3, p. 93, Sep. 2008
[50] G. M. Namata, London, L.. Getoor, and B. Huang, "Query-driven active surveying for collective classification,"" in Workshop Min. Learn. Graphs, 2012.
[51] I. Chami, Z. Ying. C. Ré, and J. Leskovec, "Hyperbolic graph convolutional neural networks," in Advances Neural Inf. Process. Syst. 2019.
[52] O. Shchur, M. Mumme, A. Bojchevski, and S. Günnemann, "Pitfalls of graph neural network evaluation," " Relational Representation Learn. Workshop, Advances Neural Inf. Process. Syst... 2018.
[53] J. McAuley, C. Targett, Q Shi, and A. van den Hengel, "Image- limage-based recommendations on styles and substitutes," in Proc. Int. ACM SIGIR Conf. Res. Develop. Inform. Retrieval, 2015, p. 43-52.
[54] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lid, and Y. Bengio, "Graph attention networks," in Proc. Int. Conf. Learn. Representations, 2018, pp. 1-12.
[55] S. Zhu, S. Pan, C. Thou, J. Wu, Y. Cao, and B. Wang, "Graph geometry interaction learning," in Advances Neural Information Processing Systems, 2020.
[56] K. He, X. Zhang, S. Ren, and Sun, "Deep residual learning for image recognition," in Proc. Conf. Comput. Vision Pattern Recognition, 2016.
[57] J. Zhu, Y. Yan, L. Zhao, M. Heimann, L. Akoglu, and D. Koutra, "Beyond homophily in graph
neural networks: Current limitations and effective designs," Advances in Neural Inf. Process. Syst. vol. 33, pp. 7793-7804, 2020.
[58] Zhu, R. A. Rossi, A. Rao, T. Mai, Lipka, N. K. Ahmed, and D. Koutra, "Graph neural networks with heterophily." in Proc. AAAI Conf. Artif. Intell., 2021, pp. 11 168-11 176.
[59] E. Chien, J. Peng, P. Li, and O. Milenkovic, "Adaptive universal generalized pagerank graph neural network." " in Int. Conf. Learn. Representations, 2021.
[60] X. Li, R. Zhu, Y. Cheng, C. Shan, S. Luo, D. Li, and W. Qian, "Finding global homophily in graph neural networks when meeting heterophily," arXiv preprint 22Xi::2055003008 2022.
[61] D. Bo. X. Wang, C. Shi, and H. Shen, "Beyond low- -reeueeyy information in graph convolu- tional networks," in Proc. AAAI Conf. Artif. Intell., 2021, pp. 3950-3957.
[62] L. Du, X. Shi, Q. Fu, X. Ma. H. Liu, S. Han, and D. Zhang, "Gbk-gnn: Gated bi-kernel graph neural networks for modeling both homophily and heterophily." in Proc. ACM Web Conf. 2022, 2022, pp. 1550-1558.
[63] S. Luan, C. Hua, Q. Lu, J. Zhu, M. Zhao, S. Zhang, X.-W Chang, and D. Precup, "Revisiting heterophily for graph neural networks," arXiv preprint 110.0022 2022.
13

---

[64] C. Bodnar, F. D. Giovanni, B. P. Chamberlain, P. Lio, and M. M. Bronstein, "Neural sheaf diffusion: A topological perspective on heterophily and oversmoothing in GNNs," in Advances Neural Inf. Process. Syst. 2022.
[65] Y. Wang, K. Yi, X. Liu, Y. G. Wang, and S. Jin, "Acmp: Allen-cahn message passing with attrac- tive and repulsive forces for graph neural networks," in Proc. Int. Conf. Learn. Representations, 2023.
[66] O. Platonov, D. Kuznedelev, A. Babenko, and L. Prokhorenkova, "Characterizing graph datasets for node classification: Beyond Nomophily-heterophh dichotomy," arXiv preprint arXiv:2209.08 1777 2022.
[67] M. Mathieu, M. Henaff, and Y. LeCun, "Fast training of convolutional networks through ffts." arXiv preprint 1312.581,, 2013.
[68] OpenAI, "Chatgpt-4," 2022, available at: https//wwwwoppnni..oom (Accessed: 10 April 2024).
[69] A. M. Cohen, Inversion Formulae and Practical Results. Boston, MA: Springer US, 2007, pp.
23 44. [Online]. Available: https://doi.or// httpg/10.100//9777 (007/978-0-387) 3--33887-855-8.. 2_
[70] K. Diethelm, N. J. Ford, and A. D. Freed, "Detailed error analysis for a fractional adams method,"
Numer. Algorithms, vol. 36, pp. 31-52, 2004.
[71] G.-h. Gao and Z.-z. Sun, "Two alternating direction implicit difference schemes for two- dimensional distributed-oriee fractional diffusion equations." Journal of Scientific Computing. vol. 66, pp. 1281-1312, 2016.
[72] A. Quarteroni, R. Sacco, and F. Saleri, Numerical mathematics. Media, 2010, vol.
Springer Science & Business
[73] I. Podlubny, Fractional Disferential Equations. Academic Press, 1999.
[74] B. Jin, B. Li, and Zhou, "Correction of high-order bdf convolution quadrature for fractional evolution equations," SIAM Journal on Scientific Computing, vol. 39, no. 6, pp. A3129-A3152, 2017.
[75] Y. Dou, K. Shu, C. Xia, P. S. Yu, and L.. Sun, "User preference-aware fake news detection," in Proc. International ACM SIGIR Conference on Research and Development in Information Retrieval, 2021.
[76] H. Pei, B. Wei, K. C.-C. Chang, Lei, and B. Yang, "Geom-GCN: Geometric graph convolu- tional networks," " in Proc. International Conference Learning Representations, 2020.
[77] Y. Yan, M. Hashemi, K. Swersky, Y. Yang, and D. Koutra. "Two sides of the same coin: Heterophily and oversmoothing in graph convolutional neural networks," in Proc. IEEE Inter- national Conference on Data Mining. IEEE, 2022, pp. 1287-1292.
[78] D. Lim, F. Hohne, X. Li, S. L.. Huang, V. Gupta, O. Bhalerao, and S. N. Lim, 'Large scale learning on non-homophilous graphs: New benchmarks and strong simple methods," pp. 887- 20 902, 2021.
[79] M. Chen, Z. Wei, Z. Huang. B. Ding, and Y. Li, "Simple and deep graph convolutional networks," in Proc. International Conference Machine Learning, 2020.
[80] B. Chamberlain, J. Rowbottom, D. Eynard, F. Di Giovanni, X. Dong, and M. Bronstein, "Beltrami flow and neural diffusion on graphs," in Advances Neural Inf. Process. Syst. 2021, pp.
1 1594-1609.
[81] Di Giovanni, J. Rowbottom, B. P. Chamberlain, T. Markovich, and M. M. Bronstein, "Graph neural networks as gradient flows," arXiv preprint arXiv:2206. 10991, 2022.
[82] W. Hu, M. Fey, M. Zitnik, Y. Dong, H. Ren, B. Liu, M. Catasta, and J. Leskovec, "Open graph benchmark: Datasets for machine learning graphs," 2021.
[83] H. Zeng, H. Zhou, A. Srivastava, R. Kannan, and V. Prasanna, "Graphsaint: Graph sampling based inductive learning method," 2020.
[84] Almira, "Müntz type thcorems I," Surveys in Approximation Theory, vol. 2007.
[85] K. Kong, J. Chen, J. Kirchenbauer, R. Ni, C. B. Bruss, and T. Goldstein, "Goat: A global transformer on large-scale graphs," in International Conference on Machine Learning. PMLR, 2023, pp. 1770-17390.
14

---

A Introduction
This supplementary material complements the main body of our paper by providing additional details and supporting evidence for the assertions made therein. The structure of this document is organized as follows:
1. A comprehensive background on fractional calculus is detailed in Appendix B.
2. Details of the FDE solvers used in paper are outlined in Appendix C, along with the corresponding approximation error analysis for the solvers in Appendix D.
3. Additional explanations for the non-Markovian random walk interpretation are provided in Appendix E.
4. An extended introduction to traditional integer-order continuous GNNs from the literature is presented Appendix F.
5. Additional implementation details, dataset specifics, and model complexity are elaborated in Appendices G and H.
6. More experimental results are available in Appendix I.
7. Theoretical results from the main paper are rigorously proven in Appendix J.
8. Limitations and broader impacts are discussed in Appendix K.
0.95
Airport
Amszzn
0.90
0.85
0.80
0.75
0.1
0.2
0.3
Figure 2: Variation of test accuracy with fractional order a in the FROND model
B Caputo Fractional Derivative
In our study, we introduce two definitions for fractional derivatives. While the elegance and inter- pretability of the Marchaud- Weyl derivative, especially its connection to random walks, is thoroughly discussed in the main paper, the practical realm of engineering often gravitates towards the Caputo fractional derivative, denoted as c [17]. Our alignment with the fractional calculus literature leads us adopt the Caputo definition in Section 3.3. This preference stems from the inherent advantage of the Caputo derivative: it naturally integrates initial conditions, elaborated in (10). The two definitions are equivalent under certain constraints [17, 29].
Below, we explore further the details of the Caputo fractional derivative to provide readers with deeper understanding. For notational simplicity this supplementary material, except Appendix we use D interchangeably with D as we solely focus on the Caputo definition in this context.
The Caputo fractional derivative of function f(t) over interval [0, 7T, of general positive order (0, oo). is defined as follows:
D DC [[[a] a) S ( ((-r) (1-7)(0l--01 (
1
(18)
Here, [a]] is the smallest integer greater than or equal a. I(-) symbolizes the gamma function, and flla)) ([(a)(()) signifies the [a]-order derivative of f. Within this definition, it is presumed that
15

---

fllall E L' L' i.e., f...i,f.fla f[(=]] is Lebesgue integrable, to ensure the well-defined nature of D f(t as per (18) [17]. When addressing vector-valued function, the Caputo fractional derivative is defined on compponent-bycommopn ease of exposition, we explicitly handle the scalar case here, although all following results can be generalized to vector-valued functions. The Laplace transform for general order a E (0, oo) is presented in [17, Theorem 7.1] as:
basis for each dimension, similar to the integer-order derivative. For
Σ
(0).
(19)
where assume that Lf exists on [so, oc) for some so R. In contrast, for the integer-order derivative fla when is a positive integer, also have the formulation (19), with the only difference being the range a. Therefore, as a approaches some integer, the Laplace transform of the Caputo fractional derivative converges the Laplace transform of the traditional integer-order derivative. As a result, we can conclude that the Caputo fractional derivative operator generalizes the traditional integer-order derivative since their Laplace transforms coincide when takes integer value. Furthermore, the inverse Laplace transform indicates the uniquely determined D = fla) (in the sense of almost everywhere [69]).
Under specific reasonable conditions, we can directly present this generalization as follows. We suppose (t) (18) is continuously differentiable. In this context, integration by parts can be utilized to demonstrate that
(- [l (t (t--)]al-a (1- + 1 fllalt (t-rt))1)) (1-r)(al- [a]]-a[[aal]]
f(t( = [[[a] a)
[a -a dr
[a]
/ (t- r)(al-a (la))
[[[a]- - a [([al-a+1)
(20)
When a + [a], we get the following
D f()(t)= f (' (()d)+
lim a-
= =flla'(0)) + f(lall(t) - (l
(21)
= (t).
In parallel to the integer-order derivative, given certain conditions ( [17,, Lemma 3.13)), the Caputo fractional derivative possesses the semigroup property:
DDD = D S.
(22)
Note, however, that in general, the Caputo fractional derivative does not possess semigroup property [17, Lemma 12]. The Caputo fractional derivative also exhibits linearity, but does not adhere to the same Leibniz and chain rules as its integer counterpart. As such properties are not utilized in our work, we refer interested readers to [17, Theorem 3.17 and Remark 3.5.J. We believe the above explanation facilitates understanding the relation between the Caputo derivative and its generalization of the integer-order derivative.
C Numerical Solvers for FDEs
In this section, we introduce basic single-term FDEs along with techniques for solving them. We also discuss multi-term FDEs and describe methods to convert them into single-term FDEs. In our paper, we approximate the distributed-orde two techniques to solve the multi-term FDE (14): one technique directly uses the single-term FDE solver, while the other approximates each fractional differential operator. For conditions necessary for the existence and uniqueness of solutions for single- and multi-term FDEs, we direct interested readers to [17, Chapter 6 and 8] and [15].
FDE (10) using the multi-term FDE (14). We present
16

---

C.1 Single-Term Solver
single-term FDE is represented as:
D() s(1,v(())
(23)
where the initial conditions take the form:
D*y(0) yo [c]
0, [a]
(24)
[k] representing the k-order derivative at point 0. with Vo
Our approach to solving (23) is based on the fractional Adams Bashforrthhhhhh Bashforth-Moulton method described in [70]. The basic predictor yk+1 is expressed as:
[al-1
Vk++
61 9o
Σ (f Uz). j=0
(25)
r(a)
Here, denotes the current iteration time step index in the discretization process, h is the step size or time interval between successive approximations with tj = hj. and yj is the numerical
approximation of (t,).. [-] represents the ceiling function, and when 0 S 1, [a] 1. The coefficients bj, 1j, are defined as follows:
ht" byAt+1 a
(
(26)
Using this predictor, it is possible derive corrector term improve the accuracy of the solver. Nonetheless, we omit this corrector term in this work and leave its detailed exploration and implica- tions for DRAGON subsequent studies,
C.2 Convert Multi-Term to Single-Term
We reference a theorem from [17] which provides a method to transform multi-term FDEs into their single-term counterparts, specifically when dealing with rational numbers.
Theorem 3. [17, Theorem 8.1.] Consider the equation
D y(x) = f (x, v(x), D()y(a)
... D ;x())
(27)
subject to the initial conditions
y 56 (0) = yo
[5]
j --
n2 1,
O, n for all 2, 2....... and 0 n1 S 1. Assume that ny E for all 3= 1,2, ..., ., k. define M to be the least common multiple of the denominators of 71, 12: ....... n2 and set
where ny -
7 1/M and N = Mng-
Then this initial value problem is equivalent to the system of equations
Do = (x)) (zvi(x) = (²(x),
(28)
D] (2),
D7 (x))) (r. (o(x), (,b
together with the initial conditions
[5/M]]
y0
Sj/M E No.
(0
0,
else.
in the following sense:
17

---

(a) Whenever Y == (yo.... -1) with E C( [0. T for some >> 0 is the solution of the system (28), the function y == yo solves the multi-term equation initial value problem (27). Here, the notation C [0, T) denotes the space of functions that have continuous m-th derivative.
(b) Whenever y &C(ne [0, [] is a solution of the multi- term initial value problem (27), the YN-1)'
vector function Y := .Uo.)..NN11) = (y. (y.Dīy D. y...... ............y..... solves the maltidimensional initial value problem (28).
C.3 Solution Strategy I for (14)
Utilizing the theorem mentioned earlier from [17], we can address the solution (14) as presented in the main manuscript. Specifically, express (14) as
n-1 Wr D x(t) = F(W, (()) Σ Dx
(29)
==
Subsequently, the single-term solver (25) Theorem can be employed to solve this equation.
C.4 Solution Strategy II for (14)
Consider the general multi-term (or more precisely, n-termm) fractional differential equation:
w y() = s(t,v()))
(30)
==
with initial condition (0) Vou where wy are coefficients, ay E (0, 1) are fractional orders, and (t) a given function.
Divide the interval [0, 7] into E equally spaced points with step size h:
=i, 0, 1,2,...... F,
where h I. The Grünwall-Letnikov approximation for the fractional derivative Doy(() with E 10, is given by:
Dty(t,)
ha Σ 1)
[v((t-))= vol.
(31)
k=
where (i) the binomial coefficient for non-integer a:
I(a + 1) T(k+ a) -
The fractional derivative Dasy(L) for each oj can be approximated as:
y(t) R
1-1)*
ai [e((ee))) vol.
has
k=
We then combine the terms for the multi-term FDE:
(aj)
w ha:
1-1)
1)
(v(ti vol
k=
equivalently,
n
Σ
111)
03 k
[u(ti-s) - voo + Σ
Wi ha, y(t,)) Σ
W
W
y0
ha
haz
7=0
k=
j=0
3=0
=((t-1.0((;=1)))
18

---

Finally, denoting the approximation of y(fz) as at each iteration, for each from 1 to E = T/h, we update the numerical solution Yi using:
y1 144-1,3 + f=0 Wj 3
yo
21= "i "in
-
(32)
Σド= W/ KJ
This provides step-by-step approach to iteratively update the solution of the n-term FDE using the Grünwald-Letmikor approximation for fractional derivatives.
Substituting the (32) into (14), we obtain the numerical solution:
X F(W, x1) + Ey Win X Σド= Wij no
X
(33)
E5 Wj nー
where X; is numerical approximation of (t:).
D Approximation Error
As discussed in Section 3.3, solving the distributed-order FDE specified in (10) involves two primary steps:
1. Discretizing the distributed-order derivative using classical quadrature rule. For instance, assuming w(a) ((x),, the application of the composite Trapezoid rule [71, 72] yields:
Aa
DO Dt) 2
w(ao) Doo X(t) + w(a;)) wlan) X(x
Σ
((Aaa)))
(34)
where Aa = (b -b- aa/n and 0j = a a jAa. After omitting smaller terms, this approximation leads to the multi-term FDE presented in (14).
2. Solving (14) using the fractional Adams-Bal Bashnfor-Bhhhhhhh-- Moullon method as described in (25) or the trinwalld-etnnker method specified in (32).
Therefore, the approximation еттог of the true solution comprises the numerical quadrature error in Step and the numerical solver enror in Step The quadrature error directly evidenced by (34). To address the solver error, we consider the general n-term FDE as detailed in (30).
For the fractional Haans-Bashforth-Moonn transformed into a system of single-term equations. This system is then solved using the method specified in (25). The approximation error for this solver is quantified as follows [35]:
method described in (25), the multi-term FDEs
((t)) - nn = O(h (nintan(()))))
max j=
(35)
where Vj denotes the value the solution time t, as computed by the numerical method, and y(t,) represents the exact solution at time 1j. h is the step size.
For the Grünwald-Letnilke method detailed in (32), we apply the Grünwald-LLtniiko approximation [73] to each fractional derivative d y(t,, which is computed as:
1 Σ(-1)
( ((((i-)) -vo]
O(h).
hay
k=
Utilizing correction techniques detailed in [74], the approximation error is calculated as:
max 1=
lb(t,) yE = O(h).
(36)
Thus, the total error is cumulative measure of the approximation emrors from both Step and Step 2.
19

---

E Non-Markovian Graph Random Walk Interpretation
Section 3.2 details the dynamics of the random walk. For enhanced clarity, here we include the corre- sponding transition probability representation for the Non-Markovian random walker at time t, which explicitly accounts for node positions throughout the entire path history (.... g(t nAT)..... q(t Ar)). Here, q(t) represents the walker's s position on the graph nodes 21. [V] at time t. This model ensures that all historical states influence transitions, emphasizing the model's 's non-Markovian nature. We consider random walker navigating over graph 9 with an infinitesimal interval of time AT O. We assume that there no self-loop in the graph topology.
For every individual value (t E (0, 1), the transition probability of the random walk dynamics as described above Fig. is characterized follows:
P (g(t)= = jt .....((t AA)) ji-ndr...... a( - (T) - j:-
(1 K) Pla (n)
if revisiting historical positions g(t - nÂT) with j: Je-nAr, i.e., the walker's wait time når and stays the same node,
K ( dis-ar
) Dh (n)
if jumping from historical positions ji-nAr to jt, i.e., the walker's wait. time nÂr and jumps to ji-nAr 's neighbour je
(37)
where K == (AT) day [--a))] is normalization coefficient, ji-nAr is the node index visited at time nAA, and th (n) the probability that the walker's 's waiting time is nÂT. For specific oa, the waiting time tha, (n) follows power-law distribution n (a,++1,, Additionally, our distributed- order fractional operator f D* (((dd((((() acts as a flexible superposition of the dynamics driven by individual fractional-order operators D This approach allows for nuanced dynamics that adapt to diverse waiting times. Theorem 2 demonstrates its capability to approximate any waiting time distribution (n) for graph-baser random walkers, thereby providing versatility modeling feature updating dynamics with varied memory incorporation levels.
F
Integer-Order Continuous GNNs
F.
GRAND and GraphCON
For the general GRAND model, the governing equation is given by:
dX(t)
(A(X()) (I).
-I)()
(38)
dl
In the case of GRAND-I, the adjacency matrix A(((I) remains constant throughout the integration process, i.e., A(X(t)) = (X)))).
For GRAND-nl, nn. the adjacency matrix A(X(t)) is time-varying and is calculated using X(t) with the attention mechanism. The entries of A(X(t)) are given by:
xX)) W Qxy
a(xi, xj) softmax
(39)
dp
where W k and Wq are learned matrices, and dy is hyperparameter determining the dimension of Wt.
GraphCON 19]: Influenced by oscillator dynamical systems, GraphCON is given by the following second-order differential equation
d²x
- - dX((t dt2
(40)
dt
where Fo(-) represents a learnable I-neighborhood coupling function, is an activation function, and 7 and / are adjustable parameters. Equivalently, we have
aY
(FF (x))) -X()) ax(() Y(1),
(41)
20

---

In the case of GraphCON-L, similar to GRAND-I, Fo(X(t),t) A(X(t)) A() For GraphCON-nl similar GRAND-nI, (x(())) ((X(t)), where ((X(t) is still obtained from (39).
F.2 Other Continuous GNNS
Heterophilic CDE [12]: Based the convection-diffusioo equation, Heterophilic CDE includes both diffusion and convection term to address information propagation from heterophilic neighbors:
dx(() (A(X(t) (t) dtv(V(t) o ((t)))
(42)
= where (ty(t) E Rd is the velocity vector associated with each edge (i. at time t, V(t) (V (Vy(t))anee (E is the edge set containing all the pairs (i,5) s.t. Waj # 0) and
i-th row of (div(V(() x())))
E
(43)
for each node E V. The velocity V (t) is given by
V v (l) (M(x,(1) (()),,
(44)
with Mi is a learnable matrix and a denotes an activation function.
GREAD: To tackle the challenges associated with heterophilic graphs, the paper [11] introduces the GREAD model. This model extends the GRAND framework by incorporating a reaction term, thereby establishing diffusion- equation for GNNs. The governing equation for this model expressed as:
dx((t -X))) -r(((())) dt
(45)
where r (((( is a reaction term, a and B3 are trainable parameters designed to balance each term.
GRAND++: Building upon the GRAND model, the paper [8] presents the GRAND++ model. This enhancement adds source term to the original GRAND framework, aimed at addressing challenges associated with training on limited labeled data. The differential equation used in GRAND++ is:
dx(() = ((X(t)) + C(O)
(46)
dt
where L(X(t)) denotes the graph Laplacian matrix, and C(00 represents a subset of X(0), consisting only of nodes identified as "Irustworthy".
G
Implementation Specifics and Dataset Details
G.1 Example Details
The distributed-order fractional model is a natural generalization of single-term as well multi- term fractional order models. It is more powerful and practical in applications. Due to multiscale characteristics in some physics problems, single-term fractional order model fails to capture this feature. Though multi-term fractional order models can capture multiscale properties, they are unsuitable in applications where the number of terms and corresponding fractional orders are unknown. However, the tistributed-orded fractional model is capable of dealing with multiscale characteristics and does not require knowing the number of terms and corresponding fractional orders priori. Graph data has a complex nature as it is from the real world. Therefore, it is natural to use a distributed-order fractional model.
Kelvin- Voigt model [33]:
(1) Erd [ DPe(I)da.
Maxwell model [31]:
(t) E De
21

---

Zener model [32]:
Gg )o(1) = a, (e {() + (1 + do (())
(1 + bo
bo
For simplicity, take E. = Exc = 1, = 1,0, = b. = 0.1,c, = 1/4 and a = 0.3 in Maxwell model,
0.2, 3 0.6 in Zener model. The toy data is generated for common (t)) = cos(t) and c(0) = 0.5. We generate the data through an open source package FractionalD- iftEqqji (https://scifrraxx. waape//raac..gggaaattiiaaaaa FractionalDiffTq, //istalle// that totally driven by Julia and licensed with MIT License. We follow standard setups and apply built-in algorithms. Specifically, we choose PIEX algorithm which is explicit method for Maxwell and Zener models, and DOMatriXDiccrete algorithm that is a strip matrix method for Kelvin- Voigt model.
For the implementation of FROND-NN and DRAGON-NN, split the data into 80% training and 20% testing sets. We construct identical two-layer neural networks with activation functions for both FROND and DRAGON. Using the current observations, predict the next 10 points in the trajectory on the test data and calculate the MSE.
G.2 Dataset Details
In this subsection, we detail the statistics of the datasets utilized in this paper, as illustrated in Tables 5 to 7. The datasets span various domains and scales, providing a comprehensive evaluation of DRAGON's performance.
Table 5: Dataset statistics used in Table 3
Dataset
Type
Classes
Features
Nodes
Edges
Cora
citation
7
1433
2485
5069
Citeseer
citation
6
3703
2120
3679
PubMed
citation
500
19717
44324
Coauthor CS
co-author
15
6805
18333
81894
Computers
so-purchase
10
767
13381
245778
Photos
co-purchase
745
7487
119043
CoauthorPhy
co-author
5
8415
34493
247962
Airport
vike
4
4
3188
3188
Disease
tree-like
1000
1044
1043
Table 6: Dataset statistics of 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 7: Dataset and graph statistics used in Table 10
Dataset
Graphs (Fake)
Total Nodes
Total Edges
Avg. Nodes per Graph
Politifact (POL)
314 (157)
41,054 314,262
40.740 308,798
131 58
Gossipcop (GOS)
5464 (2732)
H Time Complexity
In this section, we discuss the time complexity of the model, as detailed in Table and Table It is observed that the DRAGON framework exhibits computational costs comparable to those of
22

---

traditional continuous GNN models. All experiments are conducted on NVIDIA GeForce RTX 3090 or A5000 GPUs with 24GB of memory.
Table 8: Inference time of models on the Cora dataset: integral time 7
and step size of
Model
D-GRAND-I
D-GRAND-al
D-GraphCON-J
D-GraphCON-nl
D-CDE
Inf. Time(ms)
3.78
7.21
4.18
7.80
13.68
Model
F-GRAND-I
F-GRAND-nl
F.GraphCON-I
F.GraphCON-nl
F-CDE
Inf. Time(ms)
3.29
6.62
4.15
7.37
13.18
Model
GRAND-1
GRAND-nl
GraphCON-1
GraphCON-nl
CDE
Inf. Time(ms)
2.06
5.33
3.32
6.86
12.23
Table 9: Training time per epoch on the Cora dataset: integral time T = 10 and step size of
Model
D-GRAND-I
D-GRAND-nl
D-GraphCON-I D.GraphCON-nl
D-CDE
Train. Time(ms)
30.93
78.33
40.77
82.52
160.20
Model
F-GRAND-1 F-GRAND-nl F-GraphCONA1
F-GraphCON-u -n
F-CDE
Train. Time(ms)
29.76
70.31
37.82
73.10
148.92
Model
GRAND-I
GRAND-nl
GraphCON-1
GraphCON-nl
CDE
Train. Time(ms)
22.17
74.39
41.23
88.83
166.48
I More Experiment Results
I.1 Graph Classification
Following the experiments of FROND [15]. we perform graph classification tasks on the FakeNews- Net datasets [75]. The dataset features diverse array node features, including BERT embeddings, features derived from spaCy's pre-trained models, and profile-specific features from Twitter accounts. The performance outcomes, as detailed Table 10, reveal that the DRAGON-baan forms its counterparts, showcasing the significant enhancements brought about by the DRAGON framework. This because DRAGON enables feature updating dynamics with flexible memory effects stemming from the coexistence of multiple orders of derivatives.
model outper-
Table 10: Graph classification results
POL word2vec
GOS
Feature
Profile
BERT
Profile
word2vec
BERT
GruphSage
77.60=10.68
80.36-0.68
81.22-44.81
92.10=00.08
96.58=0.22
97.07-1.233
GCN
78.2811.52 03.89:0.53
83.44:10.8
89.53:1049
96.2810.08
95.9610.775
GAT
74.0311.53
78.69:0.78
82.71:0.19
91.1811.233
96.57-0.34
96.61:10.45
GRAND-I
77.8340.37
86.57=1.13
85.97-0.74
96.1110.26
97.04-0.555
96.7710.34
F-GRAND-I
79.49-10.43
08.69=0.37
89.29=0.93
96.40=0.19
97.40=0.03
97.53+0.14
D-GRAND-I
79.58.10.37
08.94-00355
0.54
97.14:1.322
97.62-10.06
97.83-10.17
I.2 Oversmoothing Mitigation
The FROND framework has demonstrated strong performance in mitigating the oversmoothing issue in GNNs [15]. As shown in Theorem 2, DRAGON can approximate any waiting time distribution, suggesting its potential to address the oversmoothing problem as well. To verify this, we conduct node classification experiments under different integration times, which can be viewed as the number of layers when the step size is set 1. From Table 11, we observe that the DRAGON framework maintains comparable performance across various depths, demonstrating consistent mitigation of the
23

---

oversmoothing issue. Furthermore, we find that DRAGON obviously outperforms FROND on the Pubmed dataset.
Table 11: Oversmoothing mitigation under fixed data splitting without using largest connected component (LCC). '.' indicates the numerical solvers failed.
Dalaset
Model
16
32
64
128
GCN
81.35::1.27
15. 3.630-3333
19.70:7.06
21.8626.00 09
13.0=0.00
13.00:0000 0.00
GAT
80.95:2.28
31.90:10.00
31.90:10.00
31.90:0.00
31.90:10.00
31.90 31.90100000
Cora
GRAND-I
81.29:00.3
81 31.5070..7
80.58+0.63
79.00=0.56
19.10:40.22
73,80+00.8
F-GRAND-I
81.17:10.75
82.68:0.64
82.24::1117
8143:+1001
81.33:40.88
00.60:10.98
D GRGAND
81.02:10.76
82.92:1.78
82.82.10.78
82.28:1.91
81.62:10.76
81.1710.74
GCN
68.84:12.46
61.58:22.09
10.64::1.9
7.70=0.00
7.70=0.00
7,70=0.00
Citeseer
GAT
65.20:10.57
18.10:10.00
18.10:10.00
18.10:10.00
18.10:10.00
18.10:10.00
GRAND-I
10
70.39-+0.68
70.52440.4
68 90++1.50
68.01+1.477
63.45-42.86
F-GRAND-I
70.68:+:.22
70. 70.70=11550
71.14 1.22
70.85:00.57
70,50:20.84 84
70.00= 60
D-GRAND-1
71.46:40.87
71.66-00.3
71.50-+0.88
71.38+11000 1.06
71.17:41.35
70.97:40.90
GCN
76.44:t:.52
72.66:22.84 84
39.52.
40.10:22.04
38 18.40::1.. 1.34
38.42:+1.7
GAT
76.98-4-1.23
40.70-+0.00
40.70++0.00
40.70-10.00
40.70++0.00
40.70-+-0.00
Pubmed
GRAND-I
77.9440.24
78.22-0.70 70
8.8444.544
F-GRAND-I
78.96::0.64 79.08=0.61
79.62:10477
79.0420.74
78.60:10.68
74.60
D D-GRAND-I
78.42:10.13
78.72:10.30
78. 78.80:10.82
78.56:10.62
79.28:10.26
79.50:10.55
1.3 D-GREAD
Building upon the GREAD model [11], we introduce D-GREAD with the following formulation:
J." D D"X(()du(c)) -oL(X(() ar(((())
(47)
Following the experimental setting in [11], we conduct a node classification task on three heterophilic graph datasets, adhering to the data split method described in [76]. The baseline results are directly reported from [] As shown in Table 12, the DRAGON framework significantly improves upon the corresponding continuous GNNs, achieving the best performance across all three datasets. Notably, even the GRAND model, which traditionally underperforms on heterophilic graph datasets, performs exceptionally well when integrated with the DRAGON framework. This demonstrates the DRAGON framework's 's capability to learn a wide range of temporal dynamics and seamlessly integrate with continuous GNNs.
Table 12: Node classification results (%) of heterophilic graph under fixed data splits [76]
Dataset
Texas
Coreell Wisconsin
Geom-GCN [76]
66.76.12.77
60.5413.67 64.5113.66
H2GCN [57]
84.86
87.65-4.998
82.70-15.28
GGCN [77]
84.86-4.55
86.86+33.29
85.68-6663
LINKX [78]
74.60-8.31
75.49-15.72
77.8415.81
GloGNN [60]
84.32-4.115
87.06++3.58
83.51+4.26
ACM-GCN [63]
87.84=4.40
88,43+3.22
85.14=1.00
GCNI [79]
77.57+3.83
80.39+3..40
77.86
CGNN [6]
71.35=4.00
74.31+7.26
66.22+7.69
GRAND [
77.68-7.25
79.4113.66
82.16-7.09
BLEND [80]
83.24
84.12.3..56
85.9511.822
Sheaf [64]
85.05
89.41-1.774
84.86-1.771
GRAFT [8]]
84.38
87.9
83.2416.49
GREAD [11]
88.92-3.72
89.41+3..0
86.49-7.115
F-GREAD
89.46=3.74
89,57+3.36
86.89+4.11
D-GRAND
86.49-3.20
90.39-13.97
90.0
D-GREAD
9..54
90.98
89.4611.22
L.4 D-GRAND++
Expanding on the GRAND+ model [8], we introduce D-GRAND++ with the following formulation:
/
D ((((t)dp(e) L(((()) C(0)
(48)
24

---

We adhere to the experimental framework outlined in the GRAND++ study, focusing specifically on the model's efficacy limited- label scenarios. The key difference in our approach is the integration of DRAGON framework. Our results in Table 13 clearly show that D-GRAND++ not only consistently outperforms the baseline GRAND++ across various tests but also shows competitive performance with F-GRAND+++
Table 13: Node classification results (%) under limited-laber scenarios
Model
pre class
Cora
Citeseer
Pubmed
Coanthor'S
Computer
Photo
GRAND+
54.94-116.09
58,95449.59
65.94-14.887
60.30+11.50 55
67.65-00.77
83.12-40.78
F-GRAND++
87,31
19,11
25.98=2.72
67.71+1991
67.65=0.37
83,12:00,88
D-GRAND++
55.8428.79
60.0:88.22
65.80:22.88
69.59:3.81
67.84=0.21
83,00:0.64
GRAND++
66.92:410.04
64.98+8.31
69.3144.87
16.53++185
74.47=1
83.71+0.90
F-GRAND++
70.09=8.36
64.98:8.31
69.37:5.30
77.97+2.33
08.85=0.96
83.71=00.0
D.GRAND++
71.21
6..10.6..83
69.97:15.28
82.24:1 2.59
79.15 0.82
83.59:11.24
GRAND++
77,801:4446
70.03++3.63
71.99:11.1
84,834:0.84
88.3341121 0.564=0.56
F-GRAND++
1.666
70.26.12.36
53.78:15.67
86.09:1.2.08
82.64 0.56
88.56.10.67
D-GRAND++
79.18- 1.22
70.83-1-.998
73.57-4-..8
88.46-100.5
0.728
8.799-00.1
GRAND++
80, 80.86.22.99 99
72.3442242
75.13=33.88
86,94+0.44 46
82.99=0.81
90.65+1.199
F-GRAND+
82.73:10.81
73,52=:1.44
77.15:12.88
87.85
83.26=0.41
91.15::0.52
D D.GAND++
82,94 1.332 1.32
74.18:10.40
77.63:1
89.52.1 0.35
83.65 0,94
91.37.10.55
GRAND++
20
12.95:11.37
73.53:13.31
79.16:1:1.37
90,80 34
85.73=0.500
93.55
GRAND++
20
14.57::100
1.778
79.96:11.66
91.03-10.77
85.78 0.43
93.55 1..38
D-GRAND+
20
84.41-00.96
73.99-42.60
79.39+ 1.42
91.98:10.33
85.81 0.669
93.28-+0.30
1,5 DGAAD
Our framework accommodates any floating value for a. Nonetheless, for the experiments presented in our main paper, we have specified a E [0, to ensure fair comparison by maintaining identical initial conditions to those utilized the original models.
For instance, we can let a range between 0 and 2, leading to the D(oscillation--GRAAND modll model:
GRAND: dtx(t) = (A(X()) (t) de
D-GRAND: fo' (I)d()))(a) (A(X(1)) ())
D(oscillation) -ariooDD-----AA---- 1 D'X(())((a) = (A( II)(()
In contrast to GRAND and D- D-GRAND, which employ the initial condition X(O)= = X. D(oscillation)- GRAND is characterized as an oscillation differential equation and adopts the initial condition (0(0)= X(0) X. However, comparing this model to GRAND F-GRAND makes it challenging to ascertain whether performance differences arise from the varied initial conditions or the incorpo- ration of distributed fractional derivatives. To preserve the D-GRANU as diffusion-type equation with the same initial condition as its counterparts, GRAND and F-GRAND, we limit to the range a 1.
We showcase preliminary results for Dioscillation-GRRADD in Table 14. The findings reveal that Seciilaation-GGRANA does not outperform GRAND or D-GRAND on these datasets, suggesting that increasing the value of a does not contribute positively to these tasks and instead elevates the model's s complexity.
Table 14: Comparison between GRAND, D-GRAND, and Dioseillatioo)-GRAAND
Cora
Citeseer
Pubmed
Airport
Disease
GRAND
1.0
73.
78.8+
80.5+99.6
3..4
D-GRAND
1.3
1.15+
79.6
9.1
93.
DIOSADD
82.6+:1.6
7.88+
78.1422.3
93.5
8.4
25

---

1.6 Sensitivity Analysis
As demonstrated in our main paper, significant advantage of the DRAGON framework its ability to learn the optimal through the adjustment of weights wy in (14). We analyze the impact of varying the number n in (14) on the final accuracy. The findings, illustrated in Table and Table 16, reveal that test accuracy remains stable despite changes in n. underscoring DRAGON's robustness against parameter selection. This stability highlights the framework's considerable improvements, as compared with the FROND framework results depicted in Fig. 2.
Table 15: Learned w, of Airport dataset
aj
Accuracy
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
1.61
1.53
1.44
1,34
1.22
1.07
0.90
0.63
0.39
0.07
98.50-+0.13
1.59
1.43
1.22
X
0.94
0.48
X
98.38:20.15
X
2.15
X
1,31
X
0.58
X
0.17
X
0.02
97.70-H0.54
1.11
1.29
X
X
-0.01
-0.03
98.06=0.39
1.85
X
0.95
X
0.35
X
98.38+00.15
X
3.09
X
X
1.69
X
0.87
X
98. 08.42
X
3.58
2.32
0.84
18-12:4
Table 16: Learned wy of Roman-empire dataset.
Accuracy
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
1.62
1.09
0.59
0.16
-0.20
-0.46
-0.62
-0.68
-0.62
-0.38
93.87=t0.41
1.30
X
0.43
X
-0.26
-0,70
-0,77
93.52.+0.40
X
1.29
X
0.46
X
-0.21
-0.64
X
-0.59
93.50=+042
0.58
0.45
-0.07
-0.09
93. 10.30333
0.64
0.31
0.12
93.22-+0.36
X
X
0.73
X
X
X
0.23
X
0.06
93.09-+0.46
X
X
X
X
0.67
0.34
0.08
03.0
1.7 Large Scale Ogb-Products dataset
To showcase the scalability of the DRAGON framework large-scale datasets, we expand our evaluation to include the Ogb-products dataset, following the experimental protocols detailed in [82]. To manage this extensive dataset effectively, we adopted mini-batci training strategy that involves sampling nodes and constructing subgraphs, as introduced by GraphSAINT [83]. The outcomes presented in Table demonstrate that the DRAGON- model outperforms others, highlighting DRAGON's efficiency and scalability.
1.8
Hy Hyperparameters
The hyperparameters employed in Table 4 are detailed in Table 18. For the hyperparameters pertaining to all other experiments, they will be disclosed alongside the code release.
J Proofs of Results
In this section, we provide detailed proofs of the results stated in the main paper.
Table 17: Node classification accuracy(%) on rgb-products dataset
Model
MLP
Node2vec Full-baich GCN
GraphSAGE
GRAND-I
F-GRAND-I
D-GRAND-I
Acc
61.06:00088
72,49:1.110
75.64=0.21
78.2910.16
75.56+00.7
76.6110078
78.81:10.99
26

---

Table 18: Hyper-parameters used in Table 4
Dataset
Model
Ir
weight decay
indrop
dropout
hidden dim
step
Roman-empire
D-CDE
0.005
0.0001
0.4
0.2
80
0.2
Wiki-conc
D-CDE
0.001
0.0001
0.4
Minesweeper
D-CDE
0.005
0.0001
0.2
64
Questions
D-CDE
0.005
0.0001
0.4
16
Workers
D-CDE
0.005
0.0001
0
0.4
64
0.5
Amazonn-rtinng
D-CDE
0.001
0.0001
0.2
0,4
128
4
0.2
J.1 Proof Theorem
Proof. Noting that En Va. (n) = 1, we subtract En=1 Va. (n)lP nAT;a,) from both sides
of (12) to yield the following:
(t) -Py(t - nAT; ao)) 1la0 (n)
(AT) daw |T(-a)))
Σ
Σ
P,(1-nAri
Wij
P2(1 nár; (ao
dl
it]
=(Ar) a-)a
Lt - nAt; an}l, Vaa (n).
==1
Divide both sides by (AT) da, (T(-o))) have
1
Py (t; a.) Pt n AT (nAr)
[T(aa))]
Σ [LP(L TAT(a=) Va.
Let År 0 and switch the limit and the summation according dominated convergence theorem (we assume the conditions are satisfied), we have
1
(11(00)) - Pt - T;0a)
dr T 1l+0
T[( [T-aa))
= =LPP(t; anll
Since [[1 - a.) = (-oo), according to (2), we have
DDa a.) = LP(t;a.).
The proof is now complete.
J.2 Proof of Theorem
Proof. Let = 0; + It is obvious that Σen 1/ri co. Let C(O, 1] be the space of continuous function on the interval [0, 1] with the co-norm. By the Müntz-Szász theorem [84], the span of { r, E R} is dense in C'(0, 1].
Consider any f E (N).. We define function 7 E C(O, 1] associated with f as follows. We set F(O) = 0, ((/nn) = f(n), e N. We then linearly interpolate between 1/n + and 1/n for any n 2 to obtain f on the remaining points of [o. 1]. Apart from 0, the function f is piecewise linear is vanishing at
and hence continuous. It is also continuous at
According to the first paragraph, for any > 0, we can find and coefficients (w,, S i N) such that
N 7(x)
Σ
"i
for any x E [0, 1).
0
27

---

Letting O, we see that |200| < e. Therefore, for any n E N. we have
f(m)
we Paa (ni)
[(=)
Wi
72
N
f
10
+
72
i=0
22c.
where w/ is defined s.t. w(da;- The proof is now complete
K
Limitations and Broader Impacts
This paper proposes a generalized framework. DRAGON, that enhances existing continuous GNNs. However, its application is currently limited to continuous GNNs. For other types of GNNs, such as graph transformers [85], they need to be transformed into the formulation differential equations before being combined with DRAGON. A future direction address this limitation is to develop a more general DRAGON framework that does not rely on numerical solvers. Regarding broader commitment ethical standards
impacts, the future societal impact of this work depends and responsible use. crucial to ensure that advancements lead to positive outcomes without compromising individual rights contributing to inequality.
L Contribution Statement
The concept of DRAGON was initially proposed by Feng Ji and the framework is fully developed by Qiyu Kang and Kai Zhao. The manuscript was written collaboratively by Kai Zhao, Xuhao Li, and Qiyu Kang. Theoretical support for FDE was provided by Feng Ji, Qiyu Kang, Xuhao Li, and Qinxu Ding. Kai Zhao was responsible for writing the implementation code and organizing the experiments. Wenfei Liang and Yanan Zhao provided experimental support. Guidance throughout the process was provided by Wee Peng Tay.
'The proof is based on the answers to a question posted by on MathOverflow http:://mathhoverf1... net/ htt///usttooo/4411//////// wwweeettteeeeeeeeeeeeeeeeeeeeeee trass- without- 0440 446221
28

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The paper's contributions and scope detailed in the abstract and introduction Section I.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We have discussed the limitations in Appendix K.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model sell-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't It acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]

---

Justification:
The theorems proposed in Section 3.2 are supported by detailed proofs provided in Appendix J.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: have included the implementation details in Appendix and the hyperpa- rameters in Appendix 1.8.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
30

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
We have provided the source code in the supplementary material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //iiss. wwtppbcccppppuueecceeeece.ccccpp Seuiidd//ooeeuublla sssobiioissilsy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification: We have included the implementation details in Appendix and the hyperpa- rameters in Appendix 1.8.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer:[Yer Yes]
Justification:
We report mean and standard deviation values in our main experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
31

---

The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out of range (c.g. negative emror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We report the computing cost Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: / /ttp///eerrpp.eeeee//ee http///lcc//hicccciideeeee
Answer: [Yes]
Justification:
We conform, in every respect, with the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification:
We have discussed broader impacts in Appendix K.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
32

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Lieenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We cite the original paper that produced the code package dataset.
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
33

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets,
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
34


Reviews:
[{"review_id": "LOzQFt246G", "reviewer": "Reviewer_uHnh", "summary": "This paper presented a new type of continuous GNN that extends and unified many continuous GNN variants. The paper mainly generalize the distributed-order fractional derivatives to continuous GNN dynamics, and it now supports a mixing of fractional derivatives within a range of continuous orders. The author also presents many explanations for the presented generalization of continuous GNNs, with theoretical support and connections to non-markov random walk. Extensive experiments that coverage node-level and graph-level tasks are conducted, showing consistent improvement over existing baselines.", "questions": "Except the questions listed above, can the author provide more computational complexity report? Table 8 and 9 only studies cora, I would like to see the report on large datasets. (large graph & large number of graphs)", "limitations": "The author has clearly stated.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The presented framework that relies on distributed-order fractional derivatives generalizes many existing methods in the continuous GNN area, providing good theoretical support and insight for future research. \n2. While some terms and math are not familiar to me, the author present the framework and main designs greatly and easy to follow. Motivations and theorem conclusions are good to understand. \n3. Extensive experiments are conducted, and the performance improvement is notable and consistent.", "weaknesses": "1. While the author has conducted extensive experiments covering various domains, the hyperparameter selection strategy is not mentioned. Given the framework introduces many additional hyperparameters, the author should clearly state how they choose all hyperparameters, whether it is solely based on validation dataset. \n2. The author mentioned how to solve DRAGON that contains the distributed-order fractional differential equation. It seems that the final computation equation is just the equation 36 mentioned in appendix. As the equation contains approximations, can the author clearly state the approximation error to the true solution? It seems that the equation 36 is a simple iterative based equation, and I'm not sure whether solving the distributed-order differential equation is that simple. \n3. The complexity analysis is not very detailed. C and E is kind of vague.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weakness 1: Hyperparameters chosen\n \n**Response:** We appreciate the reviewer's comments.\n\nLike GRAND, CDE, and FROND, we employ a grid search on the validation dataset to optimize common hyperparameters such as hidden dimensions, learning rate, weight decay, and dropout rate.\nSome details are in Table 18 in our manuscript.\n\nHere we also clarify the parameters $\\alpha_j$ in our implementation, distinct from other continuous GNNs. \nAs detailed in Sec. 4.1, we restrict $\\alpha_j$ to $[0,1]$, with values evenly spaced. \nWe consistently apply the same 10 levels of $\\alpha_j$, from 0.1 to 1.0, across all datasets, each associated with learnable weight $w_j$.  Our approach, therefore, avoids the need to extensively fine-tune the derivative order $\\alpha$, unlike the FROND model. Furthermore, in Sec. G6, we demonstrate that this discretization is not sensitive to the number of $\\alpha_j$ levels used.  \n\n## Weakness 2: Numerical solver and the approximation error\n\n**Response:** Solving the distributed-order FDE Eq.(10) consists of two steps:\n\n- **Step 1:** Discretizing the distributed-order derivative by classical quadrature rule. \n    For example, suppose $w(\\alpha)=\\mu'(\\alpha)$, applying composite Trapezoid rule [2-3] gives \n   \\begin{align*}\n           \\displaystyle\\int_a^b D^\\alpha X(t)d\\mu(\\alpha) = \\frac{\\Delta \\alpha}{2} \\left[w(\\alpha_0) D^{\\alpha_0}X(t) + 2\\sum_{j=1}^{n-1} w(\\alpha_j)D^{\\alpha_j} X(t) + w(\\alpha_n)D^{\\alpha_n}X(t)\\right]+ O((\\Delta\\alpha)^2),   \\quad  \\quad  \\quad   \\quad \\quad \\text{(S1)} \n   \\end{align*}\n    where $\\Delta\\alpha=(b-a)/n$ and $\\alpha_j=a+j\\Delta\\alpha$. After omitting small terms, we get the multi-term fractional differential equation Eq.(14).\n\n- **Step 2:** Solving the  multi-term fractional differential equation Eq.(14) using fractional Adams–Bashforth–Moulton method Eq.(28) or Grünwald-Letnikov method Eq.(36).\n\n\n__Therefore, the approximation error to the true solution consists of numerical quadrature error in Step 1 and numerical solver error in Step 2.__ The former one is clear from (S1). For the latter, we address the main idea by considering the general $n^{th}$ order multi-term fractional differential equation: $\\sum_{j=1}^n w_j D^{\\alpha_j} y(t) = f(t) $ with initial condition $y(0) = y_0$. We have\n\n> **1).** For fractional Adams–Bashforth–Moulton method Eq.(28), the multi-term fractional differential equations are equivalently transformed into a system of single-term equations, which is then addressed using the solver for single terms as detailed in Eq.(28). The approximation error for this solution is quantified as follows [1]:\n\\begin{align*}\n    \\max_{j=0,1,\\ldots,N} \\left| y(t_j) - y_j \\right| = O(h^{1+\\min\\{\\alpha_j\\}}),  \\quad  \\quad \\quad  \\quad \\quad   \\quad \\quad \\text{(S2)} \n\\end{align*}\nwhere $y_j$ denotes the value of the solution at time $t_j$ as computed by the numerical method, and $y(t_j)$ represents the exact solution at time $t_j$, $h$ is the step size.\n\n> **2).** For Grünwald-Letnikov method Eq.(36), we apply Grünwald-Letnikov approximation [4] for each fractional derivative $D^{\\alpha_j} y(t)$ that is given by:\n\\begin{align*}\n    D^{\\alpha_j} y(t_i) = \\frac{1}{h^{\\alpha_j}} \\sum_{k=0}^{i} (-1)^k \\binom{\\alpha_j}{k} [y(t_{i-k}) - y_0] +O(h). \n\\end{align*}\nUsing techniques in [5], we get the approximation error below\n\\begin{align*}\n\\max_{j=0,1,\\ldots,N} \\left| y(t_j) - y_j \\right| = O(h). \\quad  \\quad  \\quad   \\quad \\quad \\text{(S3)} \n\\end{align*}\n\nThe total error $E_{\\text{total}}$ is actually a combination of approximation error in Step 1 and Step 2.   \n\n[1].Diethelm K, et al. Detailed error analysis for a fractional Adams method.\n\n[2].Gao G, et al. Two alternating direction implicit difference schemes for two-dimensional distributed-order fractional diffusion equations.\n\n[3].Quarteroni A, et al. Numerical mathematics. \n\n[4].Podlubny I. Fractional differential equations Elsevier, 1998.\n\n[5].Jin B, et al. Correction of high-order BDF convolution quadrature for fractional evolution equations  \n\n## Weakness 3: Complexity analysis details\n\n**Response:**  Thank you for your feedback. We clarify here the notation and detail the computational complexity for our DRAGON using the above two numerical solvers.\n\n\n> The term $E = \\frac{T}{h}$ quantifies the discretization (iteration) steps necessary for the integration process. Here, $T$ represents the integration time, $h$ is the step size, and $E$ denotes the total number of iterations required.\n\n> The term $C$ denotes the computational complexity of function $\\mathcal{F}$. For instance, setting $\\mathcal{F}$ to the GRAND model for evaluating the soft adjacency matrix results in $C =O( |\\mathcal{E}| d)$, where $|\\mathcal{E}|$ represents the edge set size and $d$ the dimensionality of the features (cf. [6]). Alternatively, using the GREAD model results in $C=O((|\\mathcal{E}| + |\\mathcal{E}_2|)d + |\\mathcal{E}| d\\_{\\text{max}})$, where $|\\mathcal{E}\\_2|$ counts the two-hop edges, and $d\\_{\\text{max}}$ is the maximum degree among nodes (cf. [7]).    \n\nFor a detailed analysis of the **model complexity**, please refer to the top **global response section.**\n\n[6] Chamberlain, Ben, et al. \"Grand: Graph neural diffusion.\" ICML, 2021.\n\n[7] Choi, Jeongwhan, et al. \"Gread: Graph neural reaction-diffusion networks.\" ICML 2023.\n\n## Question 1: More computational complexity report\n\nWe follow the suggestions and prove the computational complexity on large datasets including Ogbn-arxiv as well as on a large number of graphs within the Peptides-func and Peptides-struct datasets. The results of this analysis are detailed in **Table R1, R2, and R3** in the top **global response section**. These results demonstrate that while our framework slightly increases computational costs compared to baseline continuous GNN models, it remains feasible for large graph datasets applications.", "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": "LBXrdV81B3", "reviewer": "Reviewer_mRZs", "summary": "The paper proposes the a novel framework called DRAGON, which uses a learnable probability distribution over a range of real numbers for the fractional order in graph dynamics, generalizing previous continuous GNN models. The paper provides a non-Markovian graph random walk interpretation for the DRAGON framework, assuming the feature updating dynamics adheres to a diffusion principle. It also proves that the DRAGON framework can approximate any waiting time distribution for graph random walks. The paper integrates DRAGON into several existing continuous GNN models and conducts experiments on various graph benchmarks. The results show that the DRAGON framework outperforms other methods on long-range graph datasets and improves the performance of continuous backbones on homophilic and heterophilic datasets.", "questions": "see the above weaknesses", "limitations": "the author has adequately addressed the limitations in this paper.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The paper introduces the DRAGON framework, which incorporates distributed-order fractional calculus into continuous GNNs, a novel approach that goes beyond the traditional integer-order or single fractional-order differential equations used in previous GNN models.\nThe use of a learnable measure over a range of real numbers for the fractional order is a unique feature of DRAGON, which allows for a flexible superposition of multiple derivative orders and captures complex graph feature updating dynamics beyond the reach of conventional models.\n\n2. The paper is clearly written and well-structured. The paper provides a non-Markovian graph random walk interpretation for DRAGON, which is a new perspective in understanding the dynamics of graph neural networks. The experimental results are presented in a clear and organized way, with tables and figures that help readers easily understand the performance of the DRAGON framework compared to other methods.", "weaknesses": "1. The DRAGON framework involves distributed-order fractional calculus and a learnable measure, which may increase the complexity of the model and make it more difficult to understand and implement. This complexity could potentially limit its adoption in practical applications. \n\n2. The paper mainly focuses on a few datasets for evaluation, and the experiments could be further expanded to include more diverse and challenging datasets.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weakness 1: Implementation and complexity details   \n\n**Response:**  Thank you for your insightful comments. We are pleased to provide additional explanation regarding the implementation and complexity of our framework. \n\nSolving the distributed-order FDE Eq.(10) consists of two steps:\n\n- **Step 1:** Discretizing the distributed-order derivative by classical quadrature rule. \n    For example, suppose $w(\\alpha)=\\mu'(\\alpha)$, applying composite Trapezoid rule [2-3] gives \n   \\begin{align*}\n           \\displaystyle\\int_a^b D^\\alpha X(t)d\\mu(\\alpha) = \\frac{\\Delta \\alpha}{2} \\left[w(\\alpha_0) D^{\\alpha_0}X(t) + 2\\sum_{j=1}^{n-1} w(\\alpha_j)D^{\\alpha_j} X(t) + w(\\alpha_n)D^{\\alpha_n}X(t)\\right]+ O((\\Delta\\alpha)^2),   \\quad  \\quad  \\quad   \\quad \\quad \\text{(S1)} \n   \\end{align*}\n    where $\\Delta\\alpha=(b-a)/n$ and $\\alpha_j=a+j\\Delta\\alpha$. After omitting small terms, we get the multi-term fractional differential equation Eq.(14).\n\n- **Step 2:** Solving  the multi-term fractional differential equation Eq.(14) using fractional Adams–Bashforth–Moulton method Eq.(28) or Grünwald-Letnikov method Eq.(36).\n\nFor additional details on **error analysis**, we kindly direct the reviewer to **our Response to Weaknesses 2 for Reviewer uHnh**. For the **model complexity**, please refer to the top **global response section**.\n\nWe hope this explanation alleviates concerns regarding the complexity and illustrates the DRAGON framework's potential for broad applicability.  \n\n## Weakness 2:  More diverse and challenging datasets\n\n**Response:** Thank you very much for your thoughtful suggestion. We value the chance to further elucidate the aims and outcomes of our study. Our research demonstrates how the DRAGON framework can enhance various continuous GNN models by integrating distributed-order fractional derivatives, thus boosting performance. Our experimental setup follows the methods outlined in the continuous GNNs literature.\n\nAdditionally, in the manuscript, we have broadened our tests to include the Long Range Graph Benchmark, which features datasets such as Chemistry, to underscore the unique benefits of our DRAGON.  \n\nThe broad applicability of the DRAGON framework is notable, with potential uses in fields such as traffic forecasting [1] and learning system dynamics [2]. In the rebuttal, we have adapted integer-order differential equations in [1] to distributed-order fractional derivatives, expanding our framework's reach. Specifically, we have applied the DRAGON framework to the STG-NCDE model in preliminary time-series traffic forecasting tests, as demonstrated in Table R4. This preliminary result shows our framework's adaptability and effectiveness even without extensive parameter tuning.\n\n[1] Choi J, et al. Graph neural controlled differential equations for traffic forecasting. AAAI 2022\n\n[2] Huang Z, et al. Coupled graph ODE for learning interacting system dynamics. KDD 2021.\n\n**Table R4: Forecasting error(%) for traffic time-series data**  \n\n|         | PeMSD4  | PeMSD4 | PeMSD4 | PeMSD7 | PeMSD7 | PeMSD7 |\n|--------------|------------|-------------|-------------|------------|-------------|-------------|\n| Model        | MAE | RMSE | MAPE | MAE | RMSE | MAPE |\n| STGODE       | 20.84      | 32.82       | 13.77       | 22.59      | 37.54       | 10.14       |\n| STG-NCDE     | 19.21      | 31.09       | 12.76       | 20.53      | 33.84       | 8.80        |\n| D-STG-NCDE   | 19.05      | 31.03       | 12.49       | 20.26      | 33.29       | 8.39        |", "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": "hxgwHsqUsK", "reviewer": "Reviewer_smsh", "summary": "This paper introduces a continuous Graph Neural Network (GNN) framework DRAGON that leverages distributed-order fractional calculus. Unlike traditional continuous GNNs that utilize integer-order or single fractional-order differential equations, DRAGON employs a learnable probability distribution over a range of real numbers for the derivative orders. This approach allows for a flexible superposition of multiple derivative orders, enabling the model to capture complex graph feature updating dynamics. provide a non-Markovian graph random walk interpretation under a specific anomalous diffusion setting. The authors conduct extensive experiments across various graph learning tasks and achieve competitive results.", "questions": "1. How does the computational complexity and memory usage of DRAGON compare to other continuous GNN models in practice, especially for very large graphs? I have seen the time complexity comparison on Cora, but are there more experiments on a larger dataset since cora is small.\n2. Can the authors provide their code in order for better reviewing?\n3. in table 4, the performance of D-CDE is comparable to F-CDE. the improvements are not obvious.\n4. Can the authors provide more explanations for the intrinsic property 'its ability to capture flexible memory effects' mentioned in the paper?", "limitations": "The authors discuss limitations in the Appendix I, like limitation to continuous GNNs. However, it is better that the authors have a discussion in the main context. Besides, the scalability and computational efficiency of the proposed framework in large-scale applications may need further discussion.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The use of distributed-order fractional calculus in GNNs is novel to me,  and expands the capabilities of continuous GNN models.\nThey extend previous dynamical systems to operate under DRAGON framework and generalizes the scenarios to involve distributed-order fractional derivatives\nThe authors conduct experiments on different tasks and plenty o datasets and conduct sensitivity analysis to better evaluate their proposed framework", "weaknesses": "1. Can you provide more visualizations and explanations for non-Markovian graph random walk interpretation?\n2. lack of code for reviewing and reproducibility\n3. The paper does not thoroughly address the scalability of the proposed framework in real-world applications with very large graphs.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## Weakness 1 & Question 4:  Explanation of non-Markovian graph random walk with flexible memory  \n\n**Response:**   Thank you for your valuable comments and suggestions.  \n\n**1. Further explanation for non-Markovian graph random walk with flexible memory:**\nLines 211-223 of the manuscript detail the dynamics of the random walk. \nFor enhanced clarity, here we include the corresponding transition probability $\\mathbb{P}$ representation for the non-Markovian random walker at time $t$, which explicitly accounts for node positions throughout the entire path history $(\\ldots, \\mathbf{R}(t-n \\Delta \\tau), \\ldots, \\mathbf{R}(t-\\Delta \\tau))$, as shown in Eq.(S4). Here, $\\mathbf{R}\\left(t\\right)$ represents the walker's position on the graph nodes $\\\\{\\mathbf{x}_j\\\\}\\_{j=1}^{|\\mathcal{V}|}$ at time $t$. This model ensures that all historical states influence transitions, emphasizing the model's non-Markovian nature. We consider a random walker navigating over graph $\\mathcal{G}$ with an infinitesimal interval of time $\\Delta \\tau>0$. We assume that there is no self-loop in the graph topology. The transition probability of the random walker is characterized as follows:  \n\n\n\\begin{align}\n\\begin{aligned}\n& \\mathbb{P}\\left(\\mathbf{R}(t)=\\mathbf{x}\\_{j_t} \\mid \\ldots, \\mathbf{R}(t-n \\Delta \\tau)=\\mathbf{x}\\_{j_{t-n \\Delta \\tau}}, \\ldots, \\mathbf{R}(t-\\Delta \\tau)=\\mathbf{x}\\_{j_{t-\\Delta \\tau}}\\right) \\\\\\\\\n& = \\begin{cases}\n\\left(1- K \\right)\\psi_\\alpha(n) & \\text{if revisiting historical positions $\\mathbf{R}(t-n\\Delta \\tau)$ with $j_{t}=j_{t-n\\Delta \\tau}$, i.e., the walker's wait time is $n\\Delta \\tau$ and stays at the same node,} \\\\\\\\\n\\left(K\\frac{W_{j_{t-\\Delta \\tau} j_{t}}}{d_{j_{t-\\Delta \\tau}}}\\right)\\psi_\\alpha(n)  & \\text{if jumping from historical positions $j_{t-n\\Delta \\tau}$ to $j_{t}$, i.e., the walker's wait time is $n\\Delta \\tau$ and jumps to $j_{t-n\\Delta \\tau}$'s neighbour $j_{t}$}.\n\\end{cases} \\quad  \\quad  \\quad   \\quad \\quad \\text{(S4)} \n\\end{aligned}  \n\\end{align}\nwhere $K\\coloneqq(\\Delta \\tau)^\\alpha d_\\alpha|\\Gamma(-\\alpha)|$ is a normalization coefficient, $j_{t-n\\Delta \\tau}$ is the node index visited at time $t-n\\Delta \\tau$, and $\\psi_\\alpha(n)$ is the probability that the walker's waiting time is $n\\Delta\\tau$. \nFor a specific $\\alpha$, the waiting time $\\psi_\\alpha(n)$ follows a power-law distribution $\\propto n^{-\\left(\\alpha+1\\right)}$. Additionally, our distributed-order fractional operator $\\int D^\\alpha \\mathbf{X}(t) \\mathrm{d} \\mu(\\alpha)$ acts as a flexible superposition of the dynamics driven by individual fractional-order operators $D^\\alpha$. This approach allows for nuanced dynamics that adapt to diverse waiting times. Theorem 2 in our DRAGON framework demonstrates its capability to approximate any waiting time distribution $f(n)$ for graph-based random walkers, thereby providing versatility in modeling feature updating dynamics with varied memory incorporation levels.  \n\n**2. Visualization of the random walk:**\n\nWe kindly refer the reviewer to the attached one-page PDF for the visualization.  \n\n## Weakness 2 & Question 2:  Missing implementation code (already included in our initial submission)  \n\n**Response:** We would like to clarify that the implementation code has already been included in \"Supplementary Material\" of our initial submission.   \n\n## Weakness 3 & Question 1:  Scalability of DRAGON\n\n**Response:** Thank you for your comments. We would like to clarify that our manuscript includes experimental results on the ogbn-products dataset, which is a large-scale graph with 2,449,029 nodes and 61,859,140 edges, as presented in **Table 17** of our manuscript.\n\nRegarding the scalability of our DRAGON framework, we introduce two numerical solvers, detailed in Eq.(32) and (36). For further details on error analysis and complexity, we kindly refer the reviewer to **our Response to Weaknesses 2 for Reviewer uHnh** and the **model complexity** at the top **global response section**. \n\nAdditionally, in the rebuttal, we have conducted a computational complexity comparison for the ogbn-arxiv dataset in **Table R1**. These results demonstrate that while our framework slightly increases computational costs compared to baseline continuous GNN models, it remains feasible for large graph datasets applications.    \n\n## Question 3:  Performance of D-CDE  \n\n**Response:** \n\nRegarding classification performance, we acknowledge in Table 4 that our D-CDE marginally outperforms F-CDE. However, it is crucial to emphasize a key advantage: in the baseline FROND framework, the fractional derivative order $\\alpha$ is a hyperparameter that requires tuning for each dataset. This tuning is effort-intensive and must be repeated to identify the optimal $\\alpha$ for each new dataset. As shown by our results in Fig.1. of the manuscript, performance significantly deteriorates with non-optimal $\\alpha$ values, highlighting FROND's sensitivity to this parameter.\n\nIn contrast, our DRAGON framework utilizes a set of fractional orders $\\alpha_j$ and learns their weights $w_j$ automatically, removing the need for manual tuning of $\\alpha$. This method provides a more practical and robust solution compared to FROND, especially in scenarios where manual hyperparameter tuning is impractical or costly.", "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:
We sincerely thank all the reviewers for their insightful comments and valuable suggestions. We greatly appreciate the feedback and have thoughtfully addressed each point in our detailed responses.

In this "global" response, we provide a detailed analysis of the computational complexity of the DRAGON framework as suggested by reviewers uHnh, mRZs, and smsh. Additionally, we include further numerical results that substantiate our responses.

Moreover, in the attached PDF file, we present visualizations of the non-Markovian graph random walk, as recommended by Reviewer smsh.

**Model Complexity**  
In our manuscript, the time complexity of the DRAGON framework is discussed in Section 4.4. Here, we provide a more detailed analysis concerning the two numerical solvers employed within our study.

For the Adams–Bashforth–Moulton method (Equation 28), we compute $\mathbf{X}\_{k+1} = \mathbf{X}\_0 + \frac{1}{\Gamma(\alpha)} \sum_{j=0}^{k} b_{j, k+1} \mathcal{F}(\mathbf{W}, \mathbf{X}\_j)$. This process necessitates repeated computation of $\mathcal{F}(\mathbf{W}, \mathbf{X}\_j)$ at each iteration. Direct computation leads to a complexity of $\mathcal{O}(C E^2)$. If we save the 
intermedia function evaluation values $\\{\mathcal{F}(\mathbf{W}, \mathbf{X}\_j)\\}\_j$, the total computational complexity over the entire process can be expressed as $\sum_{k=0}^E (C + O(k))$, where $O(k)$ represents the computational overhead of summing and weighting the $k$ terms at each step. We adopt this strategy in our code implementation. If the cost of weighted summing is minimal, the complexity reduces to $O(E C)$. For the Grünwald-Letnikov method (Equation 36), the computational complexity is $O(EC)$ since it requires no repeated computation of $\mathcal{F}(\mathbf{W}, \mathbf{X}_j)$, with only one $\mathcal{F}$ computation per iteration.  

> The term $E = \frac{T}{h}$ quantifies the discretization (iteration) steps necessary for the integration process. Here, $T$ represents the integration time, $h$ is the step size, and $E$ denotes the total number of iterations required.

> The term $C$ denotes the computational complexity of function $\mathcal{F}$. For instance, setting $\mathcal{F}$ to the GRAND model for evaluating the soft adjacency matrix results in $C = |\mathcal{E}| d$, where $|\mathcal{E}|$ represents the edge set size and $d$ the dimensionality of the features (cf. [6]). Alternatively, using the GREAD model results in $C=O((|\mathcal{E}| + |\mathcal{E}_2|)d + |\mathcal{E}| d\_{\text{max}})$, where $|\mathcal{E}\_2|$ counts the two-hop edges, and $d\_{\text{max}}$ is the maximum degree among nodes (cf. [7]).

This refined complexity analysis will be included in the revised paper version.

[6] Chamberlain, Ben, et al. "Grand: Graph neural diffusion." ICML, 2021.

[7] Choi, Jeongwhan, et al. "Gread: Graph neural reaction-diffusion networks." ICML 2023.

**Table R1:** Computation time of models on the Ogbn-arxiv dataset

| Model         | D-GRAND-l | D-GRAND-nl | D-GraphCON-l | D-GraphCON-nl |
|---------------|-----------|------------|--------------|---------------|
| Inf. Time(s)  | 0.083     | 0.139      | 0.141        | 0.196         |
| Train. Time(s)| 0.33      | 0.67       | 0.57         | 0.92         |


| Model         | F-GRAND-l | F-GRAND-nl | F-GraphCON-l | F-GraphCON-nl |
|---------------|-----------|------------|--------------|---------------|
| Inf. Time(s)  | 0.047     | 0.108      | 0.062        | 0.123         |
| Train. Time(s)| 0.14      | 0.53       | 0.50         | 0.59         |

| Model         | GRAND-l   | GRAND-nl   | GraphCON-l   | GraphCON-nl   |
|---------------|-----------|------------|--------------|---------------|
| Inf. Time(s)  | 0.038     | 0.099      | 0.044        | 0.105         |
| Train. Time(s)| 0.10      | 0.51       | 0.15         | 0.55           |  

***

**Table R2:** Computation time of models on the Peptides-func dataset

| Model         | D-GRAND-l | F-GRAND-l | GRAND-l |
|---------------|-----------|-----------|---------|
| Inf. Time(s)  | 0.324     | 0.275     | 0.259   |
| Train. Time(s)| 2.853     | 2.298     | 2.034   |

***

**Table R3:** Computation time of models on the Peptides-struct dataset

| Model         | D-GRAND-l | F-GRAND-l | GRAND-l |
|---------------|-----------|-----------|---------|
| Inf. Time(s)  | 0.338     | 0.269     | 0.253   |
| Train. Time(s)| 2.923     | 2.258     | 2.008   |

Meta-review: This paper presented DRAGON, a new type of continuous GNN that extends and unifies many continuous GNN variants. 
DRAGON provides a non-Markovian graph random walk interpretation. Unlike traditional continuous GNNs that utilize integer-order or single fractional-order differential equations, DRAGON employs a learnable probability distribution over a range of real numbers for the derivative orders.

The reviews are quite enthusiastic and the responses of the authors are quite solid.
Decision: Accept (spotlight)