Distributed-Order Fractional Graph Operating Network

Weakness 1 & Question 4: Explanation of non-Markovian graph random walk with flexible memory

Response: Thank you for your valuable comments and suggestions.

1. Further explanation for non-Markovian graph random walk with flexible memory: Lines 211-223 of the manuscript detail the dynamics of the random walk. For 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:

$

\begin{aligned} & \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) \\ & = \begin{cases} \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,} \\ \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}$}. \end{cases} \quad \quad \quad \quad \quad \text{(S4)} \end{aligned}

$

where $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$. For 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.

2. Visualization of the random walk:

We kindly refer the reviewer to the attached one-page PDF for the visualization.

Weakness 2 & Question 2: Missing implementation code (already included in our initial submission)

Response: We would like to clarify that the implementation code has already been included in "Supplementary Material" of our initial submission.

Weakness 3 & Question 1: Scalability of DRAGON

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.

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

Additionally, 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.

Question 3: Performance of D-CDE

Response:

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

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

Weakness 1: Implementation and complexity details

Response: Thank you for your insightful comments. We are pleased to provide additional explanation regarding the implementation and complexity of our framework.

Solving the distributed-order FDE Eq.(10) consists of two steps:

• Step 1: Discretizing the distributed-order derivative by classical quadrature rule. For example, suppose $w(\alpha)=\mu'(\alpha)$, applying composite Trapezoid rule [2-3] gives

$

       \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)} 

$

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

• 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).

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

We hope this explanation alleviates concerns regarding the complexity and illustrates the DRAGON framework's potential for broad applicability.

Weakness 2: More diverse and challenging datasets

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.

Additionally, 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.

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

[1] Choi J, et al. Graph neural controlled differential equations for traffic forecasting. AAAI 2022

[2] Huang Z, et al. Coupled graph ODE for learning interacting system dynamics. KDD 2021.

Table R4: Forecasting error(%) for traffic time-series data

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

Weakness 1: Hyperparameters chosen

Response: We appreciate the reviewer's comments.

Like 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. Some details are in Table 18 in our manuscript.

Here we also clarify the parameters $\alpha_j$ in our implementation, distinct from other continuous GNNs. As detailed in Sec. 4.1, we restrict $\alpha_j$ to $[0,1]$, with values evenly spaced. We 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.

Weakness 2: Numerical solver and the approximation error

Response: Solving the distributed-order FDE Eq.(10) consists of two steps:

• Step 1: Discretizing the distributed-order derivative by classical quadrature rule. For example, suppose $w(\alpha)=\mu'(\alpha)$, applying composite Trapezoid rule [2-3] gives

$

       \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)} 

$

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

• 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).

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

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]:

$

\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)} 

$

where $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.

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:

$

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

$

Using techniques in [5], we get the approximation error below

$

\max_{j=0,1,\ldots,N} \left| y(t_j) - y_j \right| = O(h). \quad \quad \quad \quad \quad \text{(S3)}

$

The total error $E_{\text{total}}$ is actually a combination of approximation error in Step 1 and Step 2.

[1].Diethelm K, et al. Detailed error analysis for a fractional Adams method.

[2].Gao G, et al. Two alternating direction implicit difference schemes for two-dimensional distributed-order fractional diffusion equations.

[3].Quarteroni A, et al. Numerical mathematics.

[4].Podlubny I. Fractional differential equations Elsevier, 1998.

[5].Jin B, et al. Correction of high-order BDF convolution quadrature for fractional evolution equations

Weakness 3: Complexity analysis details

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.

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 =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]).

For a detailed analysis of the model complexity, please refer to the top global response section.

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

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

Question 1: More computational complexity report

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

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)

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