Distributed-Order Fractional Graph Operating Network

Authors: kai  zhao
xuhao  li
qiyu  kang
†  feng
qinxu  ding
yanan  zhao
wenfei  liang
wee peng tay
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 04:38:29

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}.

Main paper: 
 introduction: 
 Graph Neural Networks (GNNs) have been developed to handle graph-structured data, which is prevalent in domains such as 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, as 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 or partial differential equations (ODEs/PDEs) on graphs for feature aggregation. Within these continuous GNN models, the differential operator d α dt α is typically constrained to integer values of α, primarily 1 or 2. Two directions have been proposed recently based on the aforementioned continuous GNN models to enhance their capabilities. One approach is TDE-GNN [14] , which proposes to learn higher integerorder temporal dependencies for continuous GNN models. The other approach is FROND [15] , which incorporates graph neural Fractional-order Differential Equations (FDEs), extending the conventional integer-order derivative d α dt α to encompass a positive real number α. This adaptation not only bolsters the model's efficacy but also enhances its adversarial robustness by varying the value of α [16] . TDE-GNN, however, is limited to utilizing integer-order ODEs and does not account for the nonlocal 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 a 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, 20] . In fractional scenarios, particles may remain at their current position, delaying jumps to subsequent locations with fading waiting times and leading to a non-Markovian process. In contrast, traditional integer-order differential equations are typically used to model Markovian movement of particles, as the derivative df (t) dt = lim ∆t→0 f (t+∆t)-f (t) ∆t captures the local rate of function changes. On the other hand, although FROND utilizes a fractional-order α and demonstrates performance improvement, its capacity for feature updating dynamics remains constrained by limited temporal dependencies with a single α. Moreover, the optimized performance of FROND is achieved through extensive fine-tuning of α across various graph datasets. Observations from Fig. 2 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 a range [a, b] for α. The foundation of our framework is the distributed-order fractional differential operator [22] :EQUATIONwhich can be perceived as the limiting case of i w(α i )D αi f (t), a weighted summation over derivatives of multiple orders with weight w(•) (we employ this more common notation D α instead of d α / dt α henceforth). Notably, unlike TDE-GNN, which restricts α i 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-order operator D α employed in FROND, which still has a restricted capacity to model the intricacies of feature updating dynamics. From the perspective of a random walk in a diffusion process, a single D α dictates that the waiting time between particle jumps follows a fixed power-law distribution ∝ t -α-1 for 0 < α < 1. In contrast, DRAGON adopts a more flexible approach, enabling a 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 µ(α) taking a 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 non-Markovian 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 potential as an augmentation to amplify the effectiveness of a range of continuous GNN models. 
preliminaries and related work: 
 This paper focuses on developing a new GNN framework centered around distributed-order fractional dynamic processes. In this section, we provide a concise introduction to the key concepts in fractional calculus. Throughout the paper, we adopt certain standard assumptions to ensure problem wellposedness. For instance, the well-definedness of integrations, the existence and uniqueness of the differential equation solution [23, 24] , and the allowance for interchange between summation and limit via the monotone or dominated convergence theorem [25] are all assumed.The single fractional-order operator D α in the distributed-order fractional operator in (1) can assume various definitions. In this study, we start off with the Marchaud-Weyl fractional derivative M D α , recognized for its efficacy in elucidating the fading memory phenomena [26] [27] [28] , which we will discuss further in Sections 2.1.1 and 3.2. Remark 1. However, in practical engineering implementations, the Caputo fractional derivative C D α is more commonly utilized [15, 17] . Due to space limitations, the introduction of Caputo's derivative is deferred to the Appendix B and will be subsequently employed in Section 3.3 to solve DRAGON. The Marchaud-Weyl and Caputo definitions are equivalent under certain constraints [17, 29] .For any α ∈ (0, 1), the Marchaud-Weyl α-order derivative of a function f , defined over the real line, at a specified point t is defined as [29] :EQUATIONwhere Γ(•) is the Gamma function. For sufficiently smooth functions, according to [29] , we haveEQUATIONIt is evident from (2) that the Marchaud-Weyl fractional derivative is a nonlocal operator and accounts for the past values of f within the (-∞, t) 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 not just on the current state, but also on 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 to time.We elucidate fractional-order derivatives by linking them to one-dimensional heat diffusion and memory-decaying non-Markovian random walks [28] . Assuming a random walker moves along an axis with infinitesimal intervals of space ∆x > 0 and time ∆τ > 0, the walker moves a distance of ∆x from the current point x in either direction with equal probability and waits at each location for a random period of time, a positive integer multiple of ∆τ . This introduces randomness in the waiting times between steps. We aim to compute u(x, t), the probability of the walker arriving at position x at time t. The waiting time distribution, ψ α (n), is given by a power-law functiond α n -(1+α) with d α > 0 chosen to ensure ∞ n=1 ψ α (n) = 1.The law of total probability is expressed as:u(x, t) = ∞ n=1 1 2 u(x -∆x, t -n∆τ ) + 1 2 u(x + ∆x, t -n∆τ ) ψ α (n).Here, the terms within brackets denote the probability of arriving at x from either neighboring points, x -∆x or x + ∆x, each with probability 1/2. The sum over n accounts for the possibility that the walker could have remained stationary for an extended period n∆τ with a waiting time probability ψ α (n). After subtracting ∞ n=1 ψ α (n)u(x, t -n∆τ ) from both sides and rearranging, we obtain:EQUATION)where the second-order incremental quotient is defined as:δ 2 u(x, t) = u(x -∆x, t) + u(x + ∆x, t) -2u(x, t) (∆x) 2 .In the limit as ∆x, ∆τ → 0 and assuming that (∆x) 2 dα(∆τ ) α → k α |Γ(-α)| for a positive k α [28] , we obtain the time-fractional diffusion equation:EQUATIONwhere the summations on the left-hand side of (4) converge to the integration (2). As α → 1 -, (5) reverts to the standard heat diffusion equation:EQUATIONConsequently, the aforementioned non-Markovian random walk with fading memory simplifies to the Markovian random walk, thereby eliminating the memory effects.We denote an undirected graph as G = (V, W), where V is the set of |V| = N nodes andX = [x 1 ] ⊺ , • • • , [x N ] ⊺ ⊺ ∈ R N ×d consists of rows x i ∈ R 1×das node feature vectors. The N × N adjacency matrix W := (W ij ) has elements W ij indicating the edge weight between the i-th and j-th nodes with W ij = W ji . In the subsequent GNNs inspired by dynamic processes, we letX(t) = [x 1 (t)] ⊺ , . . . , [x N (t)] ⊺ ⊺∈ R N ×d be the features at time t with X(0) = X serving as the initial condition. The time t here acts as an analog to the layer index [7, 30] . Typically, these dynamics can be described by:EQUATIONThe function F is specifically tailored for graph dynamics as illustrated in Appendix F. For instance, in the GRAND model, F is defined as follows: GRAND [7] : Drawing from the standard heat diffusion equation, GRAND formulates the following feature updating dynamics:EQUATIONwhere A(X(t)) is a learnable attention or fixed normalized matrix, and I is an identity matrix.Recently, the paper [15] introduces FROND, extending traditional integer-order graph neural differential equations such as ( 8), ( 40) and ( 42) to fractional-order equations. The framework is formalized asEQUATIONwhere F represents the graph dynamics. Further, the study in [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. 3.6 × 10 -2 3.5 × 10 -3 Kelvin-Voigt [33] 3.3× 10 -3 1.4 × 10 -4To intuitively understand the versatility and efficacy of the DRAGON framework in learning dynamics, we consider three classical stress-strain constitutive models for viscoelastic solids: the single-order Maxwell model [31] , the multi-order Zener model [32] , and the distributed-order Kelvin-Voigt model [33] . Using the FROND and DRAGON frameworks, we develop Neural Network(NN) methods to predict future states based on current observations.The detailed descriptions and implementation specifics can be found in Appendix G.1. The results presented in Table 1 demonstrate that the DRAGON framework excels in fitting not only the multiorder model but also in capturing the dynamics of single-order and distributed-order models. We observe that the DRAGON framework achieves a Mean Squared Error (MSE) that is ten times smaller than that of the FROND method across all three models. This highlights the DRAGON framework's exceptional ability to effectively learn and adapt to a diverse range of dynamics, surpassing the capabilities of FROND. 
dragon framework: 
 In this section, we introduce our general DRAGON framework for GNNs, with a 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.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 in integer-order continuous GNN models [5, 15] as presented in Section 2.2, we apply a preliminary learnable encoder function φ : V → R d that maps each node to a feature vector. After stacking all these feature vectors, we obtain X ∈ R N ×d . Employing the distributed-order fractional derivative outlined in (1), the feature dynamics in DRAGON are characterized by the following graph dynamic equation:EQUATIONwhere [a, b] denotes the range of the order α, µ is a learnable measure of α, and F is a dynamic operator on the graph as illustrated in Appendix F. Remark 2. In practical engineering settings, the Caputo fractional derivative, represented by C D α , 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 X [n] (0) = X, where X [n] (0) denotes the n-th order derivative at t = 0, encompassing the initial node features for all integers n ∈ N ∩ [0, ⌈b⌉] [23] . Here, ⌈•⌉ is the ceiling function, and this setup ensures a unique solution [23] .For instance, when [a, b] = [0, 1], we define the initial condition as X(0) = X.This framework generalizes prior continuous GNN models, encompassing them as special instances. Specifically, with µ(α) = δ(α -1), where δ is the Dirac delta function, (10) simplifies to a local first-order differential equation like [7, 8, [10] [11] [12] [13] . When [a, b] = [0, 2], we may obtain a distributedorder fractional wave propagation GNN model [21] , which generalizes the second-order GraphCON model (40) . When µ(α) = δ(α -α o ) for α o ∈ R + , (10) reduces to the FROND framework (9) . Additionally, when µ adopts a 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 as ζ(X(T )), where ζ symbolizes a learnable decoder function.In this subsection, we provide a non-Markovian graph random walk interpretation for DRAGON under a specific anomalous diffusion setting, where the dynamic operator F(W, X(t)) in ( 10) is set as (A(X(t)) -I)X(t) in ( 8) with a fixed constant matrix A. More specifically, we obtain the following linear distributed-order FDE:EQUATIONwhere we set A = WD -1 and L := WD -1 -I is the random walk Laplacian. Here, D is a diagonal matrix with D ii = d i , the degree of node i. For clarity, without loss of generality, similar to the approach in Section 2.1.1, we interpret X(t) as a N -dimensional probability or concentration vector P(t) over the graph nodes V at time t. The Marchaud-Weyl M D α employed in (11) helps expedite the exposition of the subsequent random walk, drawing an analogy from the one-dimensional random walk discussed in Section 2.1.1.For every individual value α o ∈ (0, 1), we consider a random walker navigating over graph G with an infinitesimal interval of time ∆τ > 0. We assume that there is 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 a random period of time. The distribution of waiting times, ψ αo (n), is given by a power-law function d αo n -(1+αo) with d αo > 0 chosen to ensure∞ n=1 ψ αo (n) = 1. 2.Upon deciding to make a jump, the walker can either move from the current node i to a neighboring node j with a probability of (∆τ The diagram illustrates the walker's decision-making process during the walk. After waiting for a random duration n∆, 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.) αo d αo |Γ(-α o )| Wij di if i ̸ = j.We denote P j (t; α o ), the probability of the walker being at node j at time t with a specific order α o and µ(α) = δ(α -α o ). The law of total probability is expressed as:EQUATIONIn this equation, the summation over n accounts for the possibility that the walker may have remained stationary for a period of n∆τ , with a waiting time probability of ψ αo (n). Fig. 1 provides a visualization of the non-Markovian graph random walk. For more explanation of the non-Markovian random walker on graphs, please refer to Appendix E. From (12) , we can derive Theorem 1.Theorem 1. Given µ(α) = δ(α -α o ) where α o ∈ (0, 1) and ∆τ → 0, we establish that P(t; α o ), the probability vector whose j-th element is P j (t; α o ), solves (11) . That is to say, we haveEQUATIONRemark 3. In Theorem 1, we present the graph random walk interpretation for the fractional anomalous diffusion equation (11) under the condition that µ = δ(α -α o ). This condition represents a 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.From the discussion above, for a specific α o , the waiting time is steered by the power-law distribution ∝ n -(αo+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 interpretation 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 a more modest conclusion, demonstrating that a weighted sum of ψ αi (n) can approximate any waiting time, highlighting the capability of our framework in comparison to FROND. Theorem 2. Let C 0 (N) be the space of functions on the natural numbers N vanishing at ∞, i.e.,f ∈ C 0 (N) if and only lim n→∞ f (n) = 0. Assume the sequence (α m ) ∞ m=1 is strict increasing in [0, 1], then the span of {ψ αm , m ≥ 1} is dense in C 0 (N)in the sense of uniform convergence. Remark 4. Theorem 2 demonstrates the DRAGON framework's ability to approximate any waiting 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 a fixed waiting time distribution, limiting its adaptability in modeling feature updating over time.Previous continuous GNNs have leveraged neural ODE solvers [30] when µ = δ(α -1). For example, in the explicit Euler scheme, neural ODEs are effectively reduced to residual networks with shared hidden layers [30] . Addressing the challenge of solving the distributed-order FDE (10) given by DRAGON, the standard approach involves discretizing it into a 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 3.1, we follow the convention in the fractional calculus literature for real-world applications and employ the Caputo definition C D α 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:EQUATIONwhere α j ∈ [a, b], j = 0, 1, . . . , n, are distinct interpolation points and w j are weights associated with the measure µ. Reflecting the learnable nature of µ, w j 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.1], the multi-term FDE can be transformed into a system of single-order equations C D γ , where γ := 1/M and M is the least common multiple of the denominators of α 0 , α 1 , . . . , α n when these coefficients are rational numbers. The classical fractional Adams-Bashforth-Moulton method can then be applied to solve the resulting system of single-order equations [15, 35] . This method is a 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 demonstrated 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 is discussed in Appendix D.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 our proposed DRAGON framework, which generalizes the scenarios to involve distributed-order 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 getEQUATION2. D-GraphCON: By extending (40), we getEQUATION3. D-CDE: By extending (42), we getEQUATIONwhere div(V(t) • X(t)) is given in ( 43) and (44) .Depending on the method used to compute the matrix A in ( 15), the D-GRAND model can be categorized into two versions: linear (D-GRAND-l) and non-linear (D-GRAND-nl). Similarly, based on the computation of F θ in ( 16), the D-GraphCON model also has two versions: linear (D-GraphCON-l) and non-linear (D-GraphCON-nl). Detailed explanations are provided in Appendix F.1.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes]Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes]We have provided the source code in the supplementary material.Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
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 I.4. GCN [37] 0.5930±0.0023 0.3496±0.0013 GCNII [38] 0.5543±0.0078 0.3471±0.0010 GINE [39] 0.5498±0.0079 0.3447±0.0045 GatedGCN [40] 0.5864±0.0077 0.3420±0.0013Transformer+LapPE [41] 0.6326±0.0126 0.2529±0.0016 SAN+LapPE [42] 0.6384±0.0121 0.2683±0.0043 SAN+RWSE [43] 0.6439±0.0075 0.2545±0.0012 GCN+DRew [44] 0.6996±0.0076 0.2781±0.0028 PathNN [45] 0.6816±0.0026 0.2545±0.0032 DRGNN [46] 0.6586±0.0042 0.2495±0.0015 GRAND-l 0.6962±0.0015 0.2867±0.0009 F-GRAND-l 0.7126±0.0024 0.2677±0.0014 D-GRAND-l 0.7571±0.0014 0.2461±0.0014In our approach, we employ a fully connected (FC) layer as the encoder, φ : V → R d , 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 continuous 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 identical initial conditions to those utilized in the original models. Given these considerations, we mainly restrict α j values between [0,1] in our implementation, while also balancing computational costs. The parameter α j is selected to evenly divide the entire range, aiming to comprehensively cover values between [0, 1]. Typically, we set the number of α j in (14) to 10. We also explore the empirical results when α j exceeds 1, as shown in Appendix I.5. For a sensitivity analysis of the number and value of α j , we refer the readers to Appendix I.6. Details on the datasets used can be found in Appendix G.2.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 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-l and F-GRAND-l. This demonstrates DRAGON's capability to effectively capture long-range dependencies in graph data.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), in line with the approach detailed in [7] . 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 h adj , 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.As shown in Table 3 , for homophilic datasets such as citation networks, coauthor networks, and copurchasing 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 on most datasets. The results of the graph node classification on heterophilic datasets are presented in 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.For the Adams-Bashforth-Moulton 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 F (W, X j ) for each iteration. By storing intermediate function evaluation values {F (W, X j )} j , we can express the total computational time complexity across the process asE k=0 (C + O(k)), where O(k) indicates the computational overhead from summing and weighting the k terms at each step. Here, C represents the complexity of computing F. This yields a total cost of O EC + E 2 . If a fast algorithm for the convolution computations is available, we typically require O(E 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 inC = O((|E| + |E 2 |)d + |E|d max ),where |E 2 | accounts for the number of two-hop edges, and d max is the maximum degree among nodes [11] . More details of the computation cost can be found in Appendix H. 
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 graphrelated tasks. 
experimental setting/details: 
 Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: We have included the implementation details in Appendix G and the hyperparameters in Appendix I.8. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer:[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, confidence intervals, or statistical significance tests, at 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 run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions should be given (e.g., Normally distributed errors). • It 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 a 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 or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). 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.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative 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 is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are 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 to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a 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 to 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 a best faith effort.12. Licenses for existing assets Question: 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?Answer: [Yes] Justification: We cite the original paper that produced the code package or dataset.Guidelines:• The answer NA means that the paper does not use 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 a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.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, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a 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. 
 
 Appendix: 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 our 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 in Appendix F. In our study, we introduce two definitions for fractional derivatives. While the elegance and interpretability 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 D α [17] . Our alignment with the fractional calculus literature leads us to adopt the Caputo definition in Section 3.3. This preference stems from the inherent advantage of the Caputo derivative: it naturally integrates initial conditions, as 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 a deeper understanding. For notational simplicity in this supplementary material, except in Appendix J, we use D α interchangeably with C D α , as we solely focus on the Caputo definition in this context.The Caputo fractional derivative of a function f (t) over an interval [0, T ], of a general positive order α ∈ (0, ∞), is defined as follows:EQUATIONHere, ⌈α⌉ is the smallest integer greater than or equal to α, Γ(•) symbolizes the gamma function, and f [⌈α⌉] (τ ) signifies the ⌈α⌉-order derivative of f . Within this definition, it is presumed thatf [⌈α⌉] ∈ L 1 [0, T ],i.e., f [⌈α⌉] is Lebesgue integrable, to ensure the well-defined nature of D α f (t) as per ( 18) [17] . When addressing a vector-valued function, the Caputo fractional derivative is defined on a component-by-component basis for each dimension, similar to the integer-order derivative. For ease of exposition, we explicitly handle the scalar case here, although all following results can be generalized to vector-valued functions. The Laplace transform for a general order α ∈ (0, ∞) is presented in [17, Theorem 7 .1] as:EQUATION)where we assume that Lf exists on [s 0 , ∞) for some s 0 ∈ R. In contrast, for the integer-order derivative f [α] when α is a positive integer, we also have the formulation (19) , with the only difference being the range of α. Therefore, as α approaches some integer, the Laplace transform of the Caputo fractional derivative converges to the Laplace transform of the traditional integer-order derivative. As a result, we can conclude that the Caputo fractional derivative operator generalizes the traditional integer-order derivative since their Laplace transforms coincide when α takes an integer value. Furthermore, the inverse Laplace transform indicates the uniquely determined D α f = f [α] (in the sense of almost everywhere [69] ).Under specific reasonable conditions, we can directly present this generalization as follows. We suppose f [⌈α⌉] (t) (18) is continuously differentiable. In this context, integration by parts can be utilized to demonstrate thatEQUATIONWhen α → ⌈α⌉, we get the followingEQUATIONIn parallel to the integer-order derivative, given certain conditions ( [17, Lemma 3.13]), the Caputo fractional derivative possesses the semigroup property:EQUATIONNote, however, that in general, the Caputo fractional derivative does not possess semigroup property [17, Lemma 3.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.]. We believe the above explanation facilitates understanding the relation between the Caputo derivative and its generalization of the integer-order derivative.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-order FDE (10) using the multi-term FDE (14) . We present 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] .A single-term FDE is represented as:EQUATIONwhere the initial conditions take the form:EQUATIONwith y[k] 0 representing the k-order derivative at point 0. Our approach to solving (23) is based on the fractional Adams-Bashforth-Moulton method described in [70] . The basic predictor y k+1 is expressed as:EQUATIONHere, k denotes the current iteration or time step index in the discretization process, h is the step size or time interval between successive approximations with t j = hj, and y j is the numerical approximation of y(t j ). ⌈•⌉ represents the ceiling function, and when 0 < α ≤ 1, ⌈α⌉ = 1. The coefficients b j,k+1 are defined as follows:EQUATIONUsing this predictor, it is possible to derive a corrector term to improve the accuracy of the solver. Nonetheless, we omit this corrector term in this work and leave its detailed exploration and implications for DRAGON to subsequent studies.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 equationEQUATIONsubject to the initial conditionsy [j] (0) = y [j]0 , j = 0, 1, . . . , ⌈n k ⌉ -1, where n k > n k-1 > . . . > n 1 > 0, n j -n j-1 ≤ 1 for all j = 2, 3, . . . , k and 0 < n 1 ≤ 1. Assume that n j ∈ Q for all j = 1, 2, . . . , k, define M to be the least common multiple of the denominators of n 1 , n 2 , . . . , n k and set γ := 1/M and N := M n k .Then this initial value problem is equivalent to the system of equationsEQUATIONtogether with the initial conditionsy j (0) = y [j/M ] 0 , if j/M ∈ N 0 , 0, else ,in the following sense:(a) Whenever Y := (y 0 , . . , y N -1 )T with y 0 ∈ C ⌈n k ⌉ [0, T ] for some c > 0 is the solution of the system (28), the function y := y 0 solves the multi-term equation initial value problem (27) . Here, the notation C m [0, T ] denotes the space of functions that have a continuous m-th derivative. Utilizing the theorem mentioned earlier from [17] , we can address the solution of ( 14) as presented in the main manuscript. Specifically, we can express (14) asEQUATIONSubsequently, the single-term solver (25) and Theorem 3 can be employed to solve this equation.Consider the general multi-term (or more precisely, n-term) fractional differential equation:EQUATIONwith initial condition y(0) = y 0 , where w j are coefficients, α j ∈ (0, 1) are fractional orders, and f (t) is a given function.Divide the interval [0, T ] into E equally spaced points with step size h:t i = ih, i = 0, 1, 2, . . .where h = T E . The Grünwald-Letnikov approximation for the fractional derivative D α y(t) with α ∈ (0, 1) is given by:EQUATIONwhere α k is the binomial coefficient for non-integer α:α k = Γ(α + 1) Γ(k + 1)Γ(α -k + 1).The fractional derivative D αj y(t i ) for each α j can be approximated as:D αj y(t i ) ≈ 1 h αj i k=0 (-1) k α j k [y(t i-k ) -y 0 ].We then combine the terms for the multi-term FDE:n j=0 w j 1 h αj i k=0 (-1) k α j k [y(t i-k ) -y 0 ] = f (t i-1 , y(t i-1 )), or equivalently,n j=0 w j 1 h αj i k=1 (-1) k α j k [y(t i-k ) -y 0 ] + n j=0 w j 1 h αj y(t i ) - n j=0 w j 1 h αj y 0 = f (t i-1 , y(t i-1 ))Finally, denoting the approximation of y(t i ) as y i at each for each i from 1 to E := T /h, we update the numerical solution y i using:y i = f (t i-1 , y i-1 ) + n j=0 w j 1 h α j y 0 - n j=0 w j 1 h α j i k=1 (-1) k αj k [y i-k -y 0 ] n j=0 w j 1 h α j (32)This provides a step-by-step approach to iteratively update the solution of the n-term FDE using the Grünwald-Letnikov approximation for fractional derivatives.Substituting the (32) into ( 14), we obtain the numerical solution:X i = F(W, X i-1 ) + n j=0 w j 1 h α j X 0 - n j=0 w j 1 h α j i k=1 (-1) k αj k [X i-k -X 0 ] n j=0 w j 1 h α j (33)where X i is numerical approximation of X(t i ).As discussed in Section 3.3, solving the distributed-order FDE as specified in (10) involves two primary steps:1. Discretizing the distributed-order derivative using a classical quadrature rule. For instance, assuming w(α) = µ ′ (α), the application of the composite Trapezoid rule [71, 72] yields:EQUATIONwhere ∆α = (b -a)/n and α j = a + j∆α. After omitting smaller terms, this approximation leads to the multi-term FDE presented in (14). 2. Solving (14) using the fractional Adams-Bashforth-Moulton method as described in (25) or the Grünwald-Letnikov method as specified in (32) .Therefore, the approximation error of the true solution comprises the numerical quadrature error in Step 1 and the numerical solver error in Step 2. The quadrature error is directly evidenced by (34) . To address the solver error, we consider the general n-term FDE as detailed in (30) .For the fractional Adams-Bashforth-Moulton method described in (25) , the multi-term FDEs are 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] :EQUATIONwhere 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.For the Grünwald-Letnikov method detailed in (32), we apply the Grünwald-Letnikov approximation [73] to each fractional derivative D αj y(t), which is computed as:D αj y(t i ) = 1 h αj i k=0 (-1) k α j k [y(t i-k ) -y 0 ] + O(h).Utilizing correction techniques detailed in [74] , the approximation error is calculated as:EQUATIONThus, the total error is a cumulative measure of the approximation errors from both Step 1 and Step 2.Section 3.2 details the dynamics of the random walk. For enhanced clarity, here we include the sponding transition probability representation for the non-Markovian random walker at time t, which explicitly accounts for node positions throughout the entire path history (. . . , q(t -n∆τ ), . . . , q(t -∆τ )). Here, q(t) represents the walker's position on the graph nodes {1, 2, . . . , |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 G with an infinitesimal interval of time ∆τ > 0. We assume that there is no self-loop in the graph topology.For every individual value α o ∈ (0, 1), the transition probability of the random walk dynamics as described above Fig. 1 is characterized as follows:P (q(t) = j t | . . . , q(t -n∆τ ) = j t-n∆τ , . . . , q(t -∆τ ) = j t-∆τ ) =            (1 -K) ψ αo (n)if revisiting historical positions q(t -n∆τ ) with j t = j t-n∆τ , i.e., the walker's wait time is n∆τ and stays at the same node,K Wj t-∆τ j t dj t-∆τψ αo (n) if jumping from historical positions j t-n∆τ to j t , i.e., the walker's wait time is n∆τ and jumps to j t-n∆τ 's neighbour j t .(EQUATIONwhere K := (∆τ ) αo d αo |Γ(-α o )| is a normalization coefficient, j t-n∆τ is the node index visited at time t -n∆τ , and ψ αo (n) is the probability that the walker's waiting time is n∆τ . For a specific α o , the waiting time ψ αo (n) follows a power-law distribution ∝ n -(αo+1) . Additionally, our distributedorder fractional operator D α X(t)dµ(α) 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 f (n) for graph-based random walkers, thereby providing versatility in modeling feature updating dynamics with varied memory incorporation levels.For the general GRAND model, the governing equation is given by:EQUATIONIn the case of GRAND-l, the adjacency matrix A(X(t)) remains constant throughout the integration process, i.e., A(X(t)) = A(X(0)).For GRAND-nl, 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:EQUATIONwhere W K and W Q are learned matrices, and dk is a hyperparameter determining the dimension of W k .GraphCON [9] : Influenced by oscillator dynamical systems, GraphCON is given by the following second-order differential equationEQUATIONwhere F θ (•) represents a learnable 1-neighborhood coupling function, σ is an activation function, and γ and β are adjustable parameters. Equivalently, we haveEQUATIONIn the case of GraphCON-l, similar to GRAND-l, F θ (X(t), t) = A(X(t)) = A(X(0)).GraphCON-nl, similar to GRAND-nl, F θ (X(t), t) = A(X(t)), where A(X(t)) is still obtained from (39) .Heterophilic CDE [12] : Based on the convection-diffusion equation, Heterophilic CDE includes both a diffusion and convection term to address information propagation from heterophilic neighbors:EQUATIONwhere V ij (t) ∈ R d is the velocity vector associated with each edge (i, j) at time t, V(t) = {V ij (t)} (i,j)∈E , (E is the edge set containing all the pairs (i, j) s.t.EQUATIONfor each node i ∈ V. The velocity V ij (t) is given byEQUATION)with M is a learnable matrix and σ denotes an activation function.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 a diffusion-reaction equation for GNNs. The governing equation for this model is expressed as:EQUATIONwhere r(X(t)) is a reaction term, α and β are trainable parameters designed to balance each term.Building upon the GRAND model, the paper [8] presents the GRAND++ model. This enhancement adds a 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:EQUATIONwhere L(X(t)) denotes the graph Laplacian matrix, and C(0) represents a subset of X(0), consisting only of nodes identified as "trustworthy".The distributed-order fractional model is a natural generalization of single-term as well as multiterm 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 distributed-order fractional model is capable of dealing with multiscale characteristics and does not require knowing the number of terms and corresponding fractional orders a 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] :σ(t) = Eτ γ 1 0 D α ϵ(t)dα.• Maxwell model [31] :σ(t) = E ∞ τ α D α ϵ(t).• Zener model [32] :(1 + a o b o )σ(t) = a o D α ϵ(t) + c o (1 a o b o )D β ϵ(t).For simplicity, takeE = E ∞ = 1, τ = 1, a o = b o = 0.1, c o = 1/4 and α = 0.3 in Maxwell model, α = 0.2, β = 0.6 in Zener model. The toy data is generated for a common σ(t) = cos(t) and ϵ(0) = 0.5. We generate the data through an open source package FractionalD-iffEq.jl (https://scifracx.org/FractionalDiffEq.jl/stable/) that is 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 an explicit method for Maxwell and Zener models, and DOMatrixDiscrete algorithm that is a strip matrix method for Kelvin-Voigt model.For the implementation of FROND-NN and DRAGON-NN, we 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, we predict the next 10 points in the trajectory on the test data and calculate the MSE.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. In this section, we discuss the time complexity of the model, as detailed in Table 8 and Table 9 .It is observed that the DRAGON framework exhibits computational costs comparable to those of traditional continuous GNN models. All experiments are conducted on NVIDIA GeForce RTX 3090 or A5000 GPUs with of memory. Following the experiments of FROND [15] , we perform graph classification tasks on the FakeNews-Net datasets [75] . The dataset features a diverse array of 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 in Table 10 , reveal that the DRAGON-based model outperforms its counterparts, showcasing the significant enhancements brought about by the DRAGON framework. This is because DRAGON enables feature updating dynamics with flexible memory effects stemming from the coexistence of multiple orders of derivatives. 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 to 1. From Table 11 , we observe that the DRAGON framework maintains comparable performance across various depths, demonstrating consistent mitigation of the oversmoothing issue. Furthermore, we find that DRAGON obviously outperforms FROND on the Pubmed dataset. Building upon the GREAD model [11] , we introduce D-GREAD with the following formulation:EQUATIONFollowing 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 [11] . 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 capability to learn a wide range of temporal dynamics and seamlessly integrate with continuous GNNs. Expanding on the GRAND++ model [8] , we introduce D-GRAND++ with the following formulation:EQUATIONWe 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++. Our framework accommodates any floating value for α. Nonetheless, for the experiments presented in our main paper, we have specified α ∈ [0, 1] to ensure a fair comparison by maintaining identical initial conditions to those utilized in the original models.For instance, we can let α range between 0 and 2, leading to the D(oscillation)-GRAND model: GRAND: dX(t) dt = (A(X(t)) -I)X(t) D-GRAND:1 0 D α X(t)dµ(α) = (A(X(t)) -I)X(t) D(oscillation)-GRAND: 2 0 D α X(t)dµ(α) = (A(X(t)) -I)X(t)In contrast to GRAND and D-GRAND, which employ the initial condition X(0) = X, D(oscillation)-GRAND is characterized as an oscillation-type differential equation and adopts the initial condition X ′ (0) = X(0) = X. However, comparing this model to GRAND or F-GRAND makes it challenging to ascertain whether performance differences arise from the varied initial conditions or the incorporation of distributed fractional derivatives. To preserve the D-GRAND as a diffusion-type equation with the same initial condition as its counterparts, GRAND and F-GRAND, we limit α to the range 0 < α ≤ 1.We showcase preliminary results for D(oscillation)-GRAND in Table 14 . The findings reveal that D(oscillation)-GRAND does not outperform GRAND or D-GRAND on these datasets, suggesting that increasing the value of α does not contribute positively to these tasks and instead elevates the model's complexity. As demonstrated in our main paper, a significant advantage of the DRAGON framework is its ability to the optimal α through the adjustment of weights w j in (14) . We analyze the impact of varying the number n in (14) on the final accuracy. The findings, illustrated in Table 15 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 . To showcase the scalability of the DRAGON framework to 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 a mini-batch training strategy that involves sampling nodes and constructing subgraphs, as introduced by GraphSAINT [83] . The outcomes presented in Table 17 demonstrate that the DRAGON-based model outperforms others, highlighting DRAGON's efficiency and scalability.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.In this section, we provide detailed proofs of the results stated in the main paper. Letting x = 0, we see that |w 0 | < ϵ. Therefore, for any n ∈ N, havef (n) - N i=1 w ′ i • ψ αi (n) = f ( 1 n ) - N i=1 w i • 1 n ri ≤ f ( 1 n ) - N i=0 w i • 1 n ri + ϵ ≤2ϵ.where w ′ i is defined s.t.w i = w ′ i d αi .The proof is now complete 3 .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 of differential equations before being combined with DRAGON. A future direction to address this limitation is to develop a more general DRAGON framework that does not rely on numerical solvers. Regarding broader impacts, the future societal impact of this work depends on a commitment to ethical standards and responsible use. It is crucial to ensure that advancements lead to positive outcomes without compromising individual rights or contributing to inequality.The 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 in 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 well-specification, asymptotic approximations only holding locally). The authors should reflect on 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 a few datasets or with a 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, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a 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 be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a 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, or as supplemental material.The proof is based on the answers to a question posted by F. Ji on MathOverflow https://mathoverflow. net/questions/446194/stone-weierstrass-without-the-subalgebra-condition/446221# 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 
2. 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. 
3. The complexity analysis is not very detailed. C and E is kind of vague.
questions: 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: limitations: The author has clearly stated.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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. 

2. The paper mainly focuses on a few datasets for evaluation, and the experiments could be further expanded to include more diverse and challenging datasets.
questions: questions: see the above weaknesses
limitations: limitations: the author has adequately addressed the limitations in this paper.
paper_summary: paper_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.
strength_weakness: strength_weakness: 1. Can you provide more visualizations and explanations for non-Markovian graph random walk interpretation?
2. lack of code for reviewing and reproducibility
3. The paper does not thoroughly address the scalability of the proposed framework in real-world applications with very large graphs.
questions: 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.
2. Can the authors provide their code in order for better reviewing?
3. in table 4, the performance of D-CDE is comparable to F-CDE. the improvements are not obvious.
4. Can the authors provide more explanations for the intrinsic property 'its ability to capture flexible memory effects' mentioned in the paper?
limitations: 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.

Comments: comment: Dear Reviewer smsh,

May we kindly ask if you would consider revising the overall rating based on our responses?

Sincerely,

The Authors
Dear Reviewer uHnh,

Thank you for your insightful questions, which have contributed to improving our manuscript. As the rebuttal period is drawing to a close, we kindly request your feedback on whether our responses have sufficiently addressed your concerns. We value your guidance and look forward to your further comments.

Sincerely,

The Authors
Dear Reviewer uHnh,

Thank you for your insightful questions, which have contributed to the enhancement of our paper.

Could you please let us know if our responses have addressed your concerns adequately? We appreciate your guidance and look forward to your feedback.

With sincere gratitude,

The Authors
Dear Reviewer mRZs,

We sincerely appreciate the time and expertise you have dedicated to reviewing our paper. We deeply value your support.

May we kindly ask for your confirmation on whether our rebuttal has adequately addressed your concerns? Your feedback is crucial to us.

Thank you once again for your invaluable contributions and for sharing your expertise.

With sincere gratitude,

The Authors
Thank you for the new feedback! 

We would first like to clarify potential ambiguities surrounding the term "error" by distinguishing between "performance error" and "numerical error": 

> Performance error refers to the efficacy of GNNs in tasks such as node classification, where the focus is on the model's ability to correctly predict outcomes. 

> Numerical error, on the other hand, concerns the accuracy of numerical solutions to fractional differential equations compared to the "true" solution trajectory, which is a key issue in computational mathematics.  

The two errors are related, but not fully equivalent to each other.  The approximation error in (S2) and (S3) refers specifically to the "numerical error."
For further clarity, we then conducted ablation studies on the Cora dataset to observe how errors change with respect to the step size $h$.  


> (1). In the first ablation study, we fixed all other parameters and varied the step size $h$ used in each experiment. We trained the model using different step sizes $h$ during the training phase and tested it with the corresponding step sizes during the test phase. The results are presented in **Table R5**. According to (S2) and (S3), the numerical approximation error should be large when $h$ is large. However, the results from Table R5 indicate that while classification performance deteriorates with larger step size $h$, it does not degrade to an unreasonable level and still maintains adequate classification performance. This occurs because both the training and testing phases follow the same discretization procedure, and although both are far from the true FDE solution, the loss function is designed to minimize the final classification error, resulting in satisfactory performance.  

> (2). In the second ablation study, we maintained fixed parameters as in the first study but changed our approach by training the model exclusively with a small step size $h=0.1$. During the testing phase, we varied the step size $h$. The results are presented in **Table R6**. According to (S2) and (S3), the numerical approximation error should be minimal when $h$ is small, and the solution with $h=0.1$ can be presumed close to the "true" solution trajectory of the FDE.
We noted that when $h < 1$ and remains relatively small, the model still achieves good classification performance. This occurs because the "numerical error" is still minimal, keeping the numerical solution close to the "true" solution trajectory.
However, as $h$ increases beyond this range, the "numerical error" grows significantly, diverging from both the "true" solution trajectory and the approximate solution with $h=0.1$. Consequently, classification performance deteriorates substantially.

**Table R5: Step size and classification accuracy(\%) (training and test)**

| Step size | 5     | 2     | 1     | 0.5   | 0.2   | 0.1   |
|-------------|-------|-------|-------|-------|-------|-------|
| Solver Eq.(28)   | 80.24 | 82.98 | 83.11 | 83.18 | 83.15 | 83.11 |
| Solver Eq.(36)    | 80.03 | 82.33 | 82.57 | 82.91 | 83.11 | 83.19 |


**Table R6:  Step size and classification accuracy(\%) (test phase)**

| Step size  | 5     | 2     | 1     | 0.5   | 0.2   | 0.1   |
|------------|-------|-------|-------|-------|-------|-------|
| Solver Eq.(28)   | 35.23 | 65.48 | 74.42 | 80.91 | 81.02 | 83.11 |
| Solver Eq.(36)   | 39.80 | 61.62 | 73.91 | 76.14 | 80.41 | 83.19 |
Thank you for increasing the contribution score. We appreciate this positive change! However, for clarification, we would like to note that in typical NeurIPS review processes, reviewers often adjust the overall rating rather than the contribution score alone.
Thank you for your code, responses and additional experiments. I have raised the contribution score from 2 to 3.
Thank you for these detailed response and additional measurement. Your algorithm seems to be well designed: with both theoretical support and limited computational overhead. Based on (s2) and (s3), can you give some insight to the measurement of the real error in experiments? Like how the error changes and affect the performance with respect to h in real experiments.