Data-Driven Higher Order Differential Equations Inspired Graph Neural Networks

Authors: moshe  eliasof
eldad  haber
eran  treister
carola-bibiane  schönlieb
Venue: ICLR_cc_2024_Workshop_AI4DiffEqtnsInSci
Published: 2024-02-11 07:53:35

Abstract: A recent innovation in Graph Neural Networks (GNNs) is the family of Differential Equation-Inspired Graph Neural Networks (DE-GNNs), which leverage principles from continuous dynamical systems to model information flow on graphs with built-in properties. However, existing DE-GNNs rely on first or second-order temporal orders. In this paper, we propose a neural extension to those pre-defined temporal dependencies. We show that our model, called TDE-GNN, can capture a wide range of temporal dynamics that go beyond typical first or second-order methods, and provide use cases where existing temporal models are challenged. We demonstrate the benefit of learning the temporal dependencies using our method rather than using pre-defined temporal dynamics on several graph benchmarks.

Main paper: 
 introduction: 
 In recent years, it has been shown that GNNs can be viewed as the discretization of Ordinary Differential Equations (ODE), providing interpretable behavior, such as smoothing (Poli et al., 2019; Chamberlain et al., 2021) , energy conservation (Eliasof et al., 2021; Rusch et al., 2022) , anti-symmetry (Gravina et al., 2023; 2024) , pattern formation (Wang et al., 2022; Choi et al., 2023b) , and more. We refer to this family of architectures as DE-GNNs. Many works in the field of DE-GNNs consider stationary data, that is, data that is not time dependent. As such, most works focus on the spatial interactions between nodes, while employing first-order temporal dynamics. As we show in Example 1, this can be limiting. Therefore, we study the importance of the higher-order DE-GNNs, and propose a novel mechanism to model the of the underlying ODEs of GNN layers in a data-driven fashion.The goal of this work is to develop and study a novel mechanism to model the temporal domain of DE-GNNs. Our approach, called TDE-GNN , is based on learning (i) the temporal order, and, (ii) the temporal dependency in a data-driven fashion. The temporal order defines the order of the underlying ODE as favored by the data, while the temporal dependency specifies the relationship between intermediate DE steps. To the best of our knowledge, this is the first work to study the temporal domain of DE-GNNs, in the sense that it learns higher-order DEs in a general manner. All DE-GNNs known to us, utilize either first or second-order time dependencies. In other words, existing DE-GNNs assume fixed temporal behavior, which is constant in time and is pre-defined. This shortcoming, as we show later, can be rather limiting when complex phenomena are to be modeled. Our TDE-GNN can also be viewed as an extension for Residual Networks (He et al., 2016) , which is aimed to incorporate node features from previous time steps (i.e., layers) in a learnable manner. Pendulum y 1 (t) coordinate vs. time (c) The prediction performance of naive, 1st, and 2nd order models. Higher-order models offer improved predictions. 
mathematical background and motivation: 
 Notations. We consider a graph G = (V, E), where V is a set of n nodes, and E ⊆ V × V is a set of m edges. The i-th node is associated with a possibly time-dependent hidden feature vector f i (t) ∈ R k . Let F (t) = [f 0 (t), . . . , f n-1 (t)] ⊤ be a matrix of the node state (features) at time t.Differential Equations Inspired GNNs (DE-GNNs). The basic idea of the DE-GNN family of architectures is to view GNNs as the discretization of the following ODE:EQUATIONwith initial condition F (t = 0) = F (0) , and s(F (t); G) is a spatial operator that depends on the graph G and the node features F (t). Specifically, it is common to employ graph diffusion, combined with a channel mixing operator implemented by a multilayer perceptron (MLP) . Some examples of such methods were proposed in Chamberlain et al. (2021) ; Eliasof et al. (2021) , and others. Because we focus on the temporal component of the ODE in this paper, we employ a similar spatial term that combines diffusion and channel mixing, as discussed later. Then, the graph ODE in Equation ( 1) is discretized in time, until time T , typically with the forward Euler method. The chosen discretization times are considered as GNN layers, with a total of L time steps with step size h such that T = hL.Order Matters. The tasks considered in this work predict node values. The common theme is that we view them as the prediction of the time and space evolution of the node features, given past and current node features. A popular approach is to treat the problem by combining a GNN with time series mechanisms such as LSTM (Hochreiter & Schmidhuber, 1997) or GRU (Cho et al., 2014) .While such techniques have shown promising results, we provide Example 1, where a standard GNN-LSTM is challenged, in the sense that it does not perform better than a naive solution. We attribute this shortcoming to the basic assumption of models like LSTM, that future predictions can be based on the previous state, effectively assuming first-order dynamics, which may not be sufficient to model higher-order phenomena, as shown below.Example 1. (Nonlinear Pendulum) We consider the problem of a nonlinear pendulum as a graph with two nodes, v 0 fixed at (0, 0), and, v 1 , located at (x 1 (t), y 1 (t)), illustrated in Figure 1a .Evaluating the coordinates of node v 1 at time t can be done by solving the Newtonian mechanics that define the pendulum's motion:EQUATIONwhere q(F ) is a gradient of the energy that characterizes the behavior of the pendulum. We discretize Equation (2) using the leapfrog method (Ascher, 2008) , to generate a time series data of the pendulum vertices locations. Recall that node v 0 is static, and remains in (0,0), while v 1 moves according to Equation (2). We plot the location of v 1 in Figure 1b . We define a task whose inputs are observed locations of the pendulum nodes, and the goal is to predict the future locations of nodes. We consider four possible prediction models to address this task: (i) A naive prediction model, oftentimes called a consistency model, that simply outputs the latest available state. Formally, it is given by F (l+1) = F (l) . (ii) A GNN-LSTM, similar to Seo et al. (2018) , (iii) a second order GNN-LSTM, and, (iv) our TDE-GNN limited to second order for a fair comparison, to be defined later in Section 3.We report the obtained prediction mean squared error (MSE) compared to the ground-truth data in Figure 1c . We observe that the first-order GNN-LSTM model is as limited as the naive model of consistency. To understand the limitations of GNN-LSTM, it is key to recall Equation (2), and see that a pendulum's motion involves a second-order system. However, the first-order GNN-LSTM mechanism considers only the latest state (node features) F (l) . Indeed, when considering a secondorder GNN-LSTM, one can obtain improved performance. Finally, we see that a second-order TDE-GNN offers further prediction performance improvement, as plotted in Figure 1c .This example demonstrates that the order (length) of the history used for prediction is important. If too short of a history is used, it may be impossible to accurately predict the behavior of systems with order higher than the available history length. Since for many problems, the order is unknown and can vary in time, we allow TDE-GNN to learn the order from the data. 
learning higher-order de-gnns: 
 As discussed, previous works have so far mostly considered first order time dependent GNNs as in Equation ( 1), whose forward Euler discretization reads:EQUATIONwhere F (l) ∈ R n×k are the node features at the l-th layer, h is a positive step size, and s is the spatial term specified s in Appendix F.2. Here, we study the temporal dynamics, and introduce a TDE-GNN layer that stems from the following ODE, with a maximal order of the hyperparameter o ≥ 1:EQUATIONaccompanied by the initial conditionsEQUATIONThe forward Euler discretization of Equation (4) yields our TDE-GNN layer:EQUATIONHere, we define c(H(l) o ) = [c 1 (H (l) o ), . . . , c o (H (l)o )] ∈ R o , which are learned weights based on previous node features of up to order o, formally denoted by:EQUATIONwhere ∥ denotes the stacking operation. It is important to note that Equation (6) models the dynamics of the l-th layer with o -1 previous layers, and therefore yields a discretization of an ODE of order o.For stability of computation, and for interoperability, we demand that o p=1 c p = 1. This constraint ensures that the coefficients approximate derivatives of the data up to o-th order (Evans, 1998) . In Appendix C we describe the implementation of learning the temporal coefficients c(H(l) o ).Understanding the learned coefficients c(H(l) o ).It is important to note that Equation (6) extends the idea of residual networks. Typical residual networks (as in ResNet (He et al., 2016) ) can be obtained by Equation ( 6) with o = 1 and c 1 = 1, which yields the forward Euler method that is often used to discretize diffusive GNNs (Chamberlain et al., 2021; Eliasof et al., 2021) . Also, second order oscillatory GNNs as proposed in Eliasof et al. (2021) ; Rusch et al. (2022) , can be implemented by Equation ( 6) with o = 2 and c = [c 1 , c 2 ] = [2, -1]. Overall, our TDE-GNN can implement, as well as extend, both of these types of architectures by learning higher-order dynamics with adaptive coefficients c(H (l) o ). Those extensions allow our TDE-GNN to model a diverse family of dynamics that cannot be obtained with the aforementioned methods. We provide Example 2 that shows how a third-order DE-GNN is implemented by our method. Thus, treating the temporal domain in DE-GNNs using our learnable framework allows us to reveal and understand the order of the underlying time-dependent process in a data-driven fashion. Furthermore, as shown in Example 1, such a treatment can be crucial to accurately model data that stems from higher-order phenomena.To demonstrate the efficacy of TDE-GNN, we experiment with two tasks: (i) node classification reported in Appendix D, and, (ii) spatio-temporal node forecasting, on several benchmarks. We provide benchmark details and statistics in Appendix G.2. The hyperparameters are determined using a grid search, as discussed in Appendix G.3. Because we propose two possible implementations of the temporal learning mechanism in Section Appendix C, we denote the direct parameterization variant by TDE-GNN D and the attention based parameterization variant by TDE-GNN A . A detailed description of the TDE-GNN architectures is given in Appendix F.3. We also experiment with a baseline model that we call DE-GNN and is implemented according to Equation ( 6) with o = 1 and c 1 = 1, that is, it considers only first-order dynamics, similarly to existing GNNs inspired by DEs. Our code is available at https://github.com/MosheEliasof/TDE-GNN.We experiment with spatio-temporal datasets, where the goal is to forecast future node values, given time-series data. We use the Chickenpox-Hungary, PedalMe-London, and Wikipedia-Math datasets from Rozemberczki et al. (2021b) . We use incremental training, mean-squared-error (MSE) loss, and testing procedure from Rozemberczki et al. (2021b) . We report the prediction performance, in terms of MSE, in Table 1 , considering recent methods like DCRNN (Li et al., 2018) , GConv (Seo et al., 2018) , GC-LSTM (Chen et al., 2018a) , DyGrAE (Taheri et al., 2019; Taheri & Berger-Wolf, 2019) , EGCN Pareja et al. (2020) , A3T-GCN (Zhu et al., 2020a) , T-GCN (Zhao et al., 2019) , MPNN LSTM (Panagopoulos et al., 2021) , and AGCRN (Bai et al., 2020) . Our results are reported in Table 1 , showing improvement over existing temporal GNN models, as well as the baseline of DE-GNN. This result further highlights the importance of learning higher-order dynamics offered by our TDE-GNN. The influence of the order o. As shown in Example 1, having sufficiently high order can be crucial to modeling complex data. In that example, we have demonstrated this significance via a synthetic task where the order is known. We now supplement this study by reporting the obtained performance as a function of the order o on real-world datasets, where the exact order of the underlying process that generated the data is unknown. To provide a comprehensive study of the impact of o, we report the results on both the node classification and spatio-temporal node forecasting tasks considered in this work, on several datasets. The results are reported in Figure 2 . Our results indicate, in congruence with Example 1, that higher-order models can improve performance compared to using a first-order model only. We also note that interestingly, for the Chickenpox-Hungary dataset, we see that a second-order model performs almost as well as a third, fourth, or fifth-order model. While we do not know the exact order of such a real-world dataset, the empirical results may hint at the actual underlying process of the spread of the Chickenpox disease in this dataset being second order.Inspecting c(H(l) o ).In Example 1, the underlying process of the data is known to be second-order, we now inspect and analyze the obtained coefficients c for a varying order o ∈ {2, 3, 4, 5}. As we show in Appendix G.5, it is possible to verify that the learned coefficients yield a valid discretization of the second-derivative operator ∂ 2 F ∂t 2 , revealing the true order of the pendulum's motion.Graph Neural Networks Inspired by Differential Equations. Adopting the interpretation of convolutional neural networks (CNNs) as discretizations of ODEs and PDEs (Ruthotto & Haber, 2018; Chen et al., 2018b; Zhang et al., 2019) to GNNs, works like GCDE (Poli et al., 2019) , GODE (Zhuang et al., 2020) , GRAND (Chamberlain et al., 2021) , PDE-GCN D (Eliasof et al., 2021) , GRAND++ (Thorpe et al., 2022) and others, propose to view GNN layers as time steps in the integration of the non-linear heat equation. This perspective allows to control the diffusion (smoothing) in the network, to understand oversmoothing (Nt & Maehara, 2019; Oono & Suzuki, 2020; Cai & Wang, 2020) 2023a) that advection and diffusion can improve weather forecasting performance. However, all the considered works discussed here utilize first-order temporal dynamics, while focusing on the spatial term of the ordinary differential equation. In this paper, we explore and study the temporal domain in the context of DE-GNNs, and show its importance to model complex systems and improve performance.Example 2. (3rd Order TDE-GNN) We now draw a link between a third-order TDE-GNN (with an order hyperparameter o = 3), and a third-order ODE. Note that every set of coefficients {c 1 , c 2 , c 3 } that sum to 1, with a step size h = 1 can be spanned by the basis:EQUATION)with the constraint 3 i=1 α i = 1. Note that the vectors that multiply α 1 , α 2 and α 3 correspond to a first-order, second-order, and third-order finite difference, respectively (i.e., ∂F (t) ∂t , ∂ 2 F (t) ∂t 2 , and∂ 3 F (t)∂t 3 ). Since the basis in Equation (8) is complete, for any c 1 , c 2 , c 3 that sum to 1, there exists a 3rd order differential equation whose discretization yields the same coefficients. o ) can be any real-valued number. These requirements offer both training stability (due to the normalization in requirement (i)), and the approximation of a finite order derivative (Evans, 1998) . We now discuss two implementations that we consider in this paper, and later, in our experiments in Section 4, we compare their performance.Direct parameterization. Perhaps the most intuitive implementation is obtained by direct parameterization, where we directly learn a vector c ∈ R o , and divide it by the sum of its entries (to satisfy requirement (i)), leading to the temporal coefficients vector:EQUATIONNote that in this case, the coefficients vector c is not directly influenced by the historyH (l)o , but it is still optimized according to the history via backpropagation.Attention-based parameterization. In addition to the direct parameterization, we propose a novel mechanism that leverages an attention mechanism as in Vaswani et al. (2017) . A key feature of the attention mechanism is that it outputs a pairwise score of its input. The novelty here is to apply the attention mechanism on the temporal dimension. To this end, by collecting and appropriately shaping the node features of the previous o layers, denoted by H (l) o ∈ R o×n×k as defined in Equation ( 6), and feeding it to an attention layer (Vaswani et al., 2017) , one obtains a pairwise score map S ∈ [0, 1] o×o . The last row in S represents the temporal scores of the l-th layer with the previous layers o -1. Clearly, requirement (i) is met by the SoftMax function used in the attention mechanism in Vaswani et al. (2017) . However, a SoftMax function yields non-negative pairwise values, which do not satisfy requirement (ii). Such a limitation will prevent, for instance, the ability to implement the oscillatory equation discussed in 3 using our TDE-GNN. Therefore, we follow the same implementation as in Vaswani et al. (2017) up to the SoftMax step. Instead, we only apply a normalization of the obtained pairwise interaction map by dividing it by its sum. This procedure satisfies both requirements (i) and (ii). We provide further details about the implementation in Appendix F.1. Initial conditions. When considering high-order ODEs, the aspect of the initial conditions of the model is important (Ascher & Petzold, 1998) . We consider two use cases that are treated differently. First, when solving a stationary problem such as node classification, where only a single initial temporal condition is available, we use o MLPs to embed this single state into o states, and then use the network in Equation ( 6). This initialization, as well as the application of a TDE-GNN at the first layer, is illustrated in Figure 3 .For stationary problems (e.g., node classification), the input consists of a single time step at time T 0 . Given input features I (0) ∈ R n×kin , we embed them using an MLP to obtain the initial conditions of the ODE, denoted by F (0) ∈ R n×k . For spatio-temporal tasks (e.g., forecasting node quantities), the node features [I (0) , . . . , I (r) ] are provided at sampled times [T 0 , . . . , T r ], and embedded in latent space. The node features at the final GNN layer F (L) are then fed to a classifier to output the desired shaped prediction to be compared with the labeled data, depending on the task. Note that for time series graph problems where we have a time series as input, we use at least o historical data in order to initialize the states. In this case, the frequency of the observed input can be different than the frequency of the hidden space F (l) that discretizes the ODE, in the sense that we can use more hidden layers than observed inputs. Upon receiving at least o input observations, we embed them using an MLP to obtain o hidden initial conditions. In Figure 4 , we illustrate the described process as well as the application of a TDE-GNN layer to the inputs. Past the initialization step, in both stationary and non-stationary cases, the features update relies on previously computed hidden node features, as described in Equation ( 6).Complexity. Compared to existing DE-GNN methods, our TDE-GNN involves additional o -1 additions of previous node features, to allow modeling differential equations of order o, and achieve improved performance, as we show in Section 4. If the coefficients c are obtained using the direct parameterization, then o -1 scalar multiplications are required. If the attention based mechanism is utilized to learn and evaluate c, adding O(n • k • o 2 ) multiplications. We note that o, the order hyperparameter of TDE-GNN, is typically significantly smaller than the number of channels k and nodes n, because it is bounded by the number of layers L. In Appendix G.4 we report the training and inference runtimes of the proposed implementations.Properties of TDE-GNN. Our TDE-GNN draws inspiration from a stable discrete process of ODE integration, that generates future time values, which is regarded as the node features evolution throughout the layers. Therefore, a natural question that arises is whether the obtained network is stable. Indeed, if the proposed architecture is unstable, then it may be difficult to fit the data, or, the network may not generalize well (see Haber & Ruthotto (2017) for stability definition and a thorough discussion). To this end, we prove the following theorem in Appendix E.Theorem 1 (Stability of TDE-GNN). For the discretization of Equation (6), there exists a vector c = [c 1 , . . . , c o ] such that the discrete solution is stable.In our study in Section 4, we verify Theorem 1, and show that the learnable weights c can be interpreted as finite difference derivatives.We experiment with homophilic and non-homophilic datasets. The homophilic datasets are Cora (McCallum et al., 2000) , Citeseer (Sen et al., 2008) , and Pubmed (Namata et al., 2012) . The nonhomophilic datasets are Chameleon, Squirrel, and Film from Rozemberczki et al. (2021a) . In all cases, we use the 10 splits from Pei et al. (2020) , and report their average accuracy and standard deviation in Table 2 . We consider three types of baselines: (i) 'general' GNN architectures, such as GCN (Kipf & Welling, 2017) , GAT (Veličković et al., 2018) , GCNII (Chen et al., 2020) , Geom-GCN (Pei et al., 2020) , NSD (Bodnar et al., 2022) , GGCN (Yan et al., 2021) , H2GCN (Zhu et al., 2020b) , FAGCN (Bo et al., 2021) , GPRGNN (Chien et al., 2021) , GRAFF (Giovanni et al., 2023) , LINKX (Lim et al., 2021) , and ACMII (Luan et al., 2022) . (ii) GNNs inspired by differential equations (DEs), including: GRAND (Chamberlain et al., 2021) , PDE-GCN (Eliasof et al., 2021) , GRAND++ (Thorpe et al., 2022) , GREAD (Choi et al., 2023b) , CDE (Zhao et al., 2023) , and FLODE (Maskey et al., 2023) . The common theme of those methods is that all consider first-order temporal behavior with o = 1, c 1 = 1, except for PDE-GCN, which considers a second-order model, whereo = 2, c = [c 1 , c 2 ] = [2, -1].In contrast, our TDE-GNN can learn the vector of coefficients c with maximal order o, and unless otherwise specified, we set the order to be the number of layers in the network, i.e., o = L. The third baseline we consider is (iii) a vanilla version of our TDE-GNN, where o = 1, c 1 = 1, and we call this variant DE-GNN. Our results in Table 2 suggest that for homophilic graphs which are known to benefit from diffusion (Gasteiger et al., 2019) , TDE-GNN performs similarly to other first-order differential equations inspired GNNs, including our baseline DE-GNN, since a diffusion process can be described using a first-order ODE. We find that the significance of learning higher-order dynamics is more pronounced for non-homophilic graphs which may stem from more complex phenomena than homophilic graphs. For instance, we find that our TDE-GNN A achieves an accuracy of 78.48%, compared to the baseline vanilla, first order DE-GNN with 70.99% -a considerable improvement. We now prove Theorem 1 from the main paper.Proof: Let us write the discrete DE in Equation ( 6) explicitly asEQUATION)Similar to the proofs for multistep ODE methods (Ascher et al., 1995) , to prove the stability of the method, assuming the Jacobians of s have non-positive real part1 , it is sufficient to consider only the temporal term:F (l+1) = c o F (l-o+1) + . . . + c 1 F (l)This is a linear, constant-coefficient differential equation, and it must be stable for Equation ( 6) to be stable. Also, as common in proofs of multistep ODE methods (Ascher et al., 1995) , we start from a solution of the formF (l) = ξ l(meaning ξ to the power of l). Substituting we obtainξ l+1 = c o ξ l-o+1 + . . . c 1 ξ lDividing by ξ l-o+1 we obtain the polynomial equationξ o+1 -c 1 ξ o -. . . -c o = 0 This is a polynomial of degree o + 1 with coefficients [1, -c 1 , . . . , -c o ] where c i = 1. Let ρ(c)be the roots of the polynomial. It is straightforward to see (by substitution) that 1 is a root of the polynomial. For the solution to be stable, we need to have that:F (n) = ξ n |F (n) | ≤ |F (n-1) | |ξ n | ≤ |ξ n-1 | → |ξ| ≤ 1In multi-step methods for ODEs, this condition is referred to as the root condition. Furthermore, as shown in Ascher et al. (1995) , since the coefficients c are to be determined (learned, in the case of TDE-GNN), there always exists a set of coefficients such that the root condition is satisfied. □Remark 1 While it is difficult to verify the root condition analytically, it is possible to compute it numerically, thus revealing the order of the process we learn, as we also show in G.5.We now describe the implementation of TDE-GNN A , i.e., the TDE-GNN implemented by an attention mechanism. Namely, to learn the dynamics between the node features at a current layer l and the previous o -1 layers, we utilize a multi-head self-attention mechanism (Vaswani et al., 2017) that assigns scores between the considered layer l and the o -1 layers. The difference in our implementation compared to a standard attention module as in Vaswani et al. (2017) is that we remove the SoftMax normalization step, as discussed in C, and it is required to allow both positive and negative numbers. We denote the attention mechanism by MHA. The MHA computes a score for each pair of layers (l i , l j ) ∈ o × o. As an input to the MHA, we use the history feature tensor as described in Equation ( 7), which is comprised of the stacking of the current and previous (layer-wise) o -1 node features. Then, the output of the attention module is given by:EQUATIONThe (l i , l j )-th entry in So represents the connection between the l i -th and l j -th layers. Specifically, the last row of S (l) o represents the connection between the current layer l and the previous o -1 layers. Therefore we define the unnormalized layer coefficients vector as the last row of S(l) o , that can be extracted using the following Python notations:EQUATIONIn order to satisfy condition (i) that demands the weights to sum to 1, described in C, we also add a normalization step, such that the coefficients are defined as:EQUATIONThe temporal mechanism developed in this paper is generic and can possibly be applied to various DE-GNNs, and is the novelty of our work. However, a GNN inspired by differential equations, and therefore also our TDE-GNN, is not complete without the spatial term that propagates node features across the graph. As discussed in Equation (1), our spatial aggregation function s is based on the combination of a channel mixing operation realized by an MLP and the symmetric normalized graph Laplacian L = D -1 2 (D -A)D -1 2 where D is the degree matrix, and A is the adjacency matrix of the graph G. Formally, the spatial aggregation is given by:EQUATION)where σ is a non-linear activation function, ReLU in our implementation, and h is a positive step size, and W (l) ∈ R k×k are the learnable weights of the MLP.Therefore, substituting the prescribed temporal in Equation ( 14) into Equation ( 6) leads to our TDE-GNN layer, as follows:EQUATIONF.3 ARCHITECTURE DETAILS Node classification architecture. We now elaborate on the TDE-GNN architecture for stationary problems, as used in our node classification experiments. The overall architecture flow is similar to standard GNN architectures for node classification, such as GCN (Kipf & Welling, 2017) and GCNII (Chen et al., 2020) . It is composed of initial embedding layers e 1 , . . . , e o , followed by L TDE-GNN layers, and a classifier implemented by a linear layer denoted by e out . The complete flow of this architecture is described in Algorithm 1. We train this architecture on node classification datasets by minimizing the cross-entropy loss between the ground-truth node labels Y and the predicted node labels Ỹ.Algorithm 1 TDE-GNN for stationary problems with order o.Input: Node features I (0) ∈ R n×kin Output: Predicted node labels Ỹ ∈ R n×kout 1: procedure TDE-GNN 2:I (0) ← Dropout(I (0) , p) 3:F (-o+1) = e 1 (I (0) ); F (-o+2) = e 2 (I (0) ); . . . ; F (0) = e o (I (0) ) 4:Initialize history tensor H Return Ỹ 14: end procedure Spatio-Temporal Node Forecasting. The typical task in spatio-temporal datasets is to predict future quantities (e.g., driving speed) given several previous time steps (also called frames). Formally, one is given an input tensor I in temporal = [I (0) , . . . I (r) ] ∈ R n×rkin , where r is the number of input (observed) time frames, and the goal is to predict a time frames ahead, i.e., the ground-truth is given by I gt temporal = [I (r+1) , . . . , I (r+a) ] ∈ R n×akin . This is in contrast to stationary datasets such as Cora (McCallum et al., 2000) , where input node features I (0) ∈ R n×kin are given, and the goal is to fit to some ground-truth Y ∈ R n×kout which can also be of different dimensionality in its output space. In this context, a stationary dataset can be thought of as setting r = a = 1 for the non-stationary settings. We show the overall flow of our TDE-GNN architecture for non-stationary problems in Algorithm 22 .In this architecture, we update the hidden state feature matrix F Similarly to Attention models (Vaswani et al., 2017) , we incorporate time embedding based on the concatenation of sine and cosine function evaluations with varying frequencies multiplied by the time of the input frames, as input to our TDE-GNN, denoted by T emb ∈ R n×rk temb , where we choose the number of frequencies to be 10, and by the concatenation of both sine and cosine lead to k temb = 20. We note that the time embedding is computed in a pre-processing fashion. To initialize the hidden feature matricesF (0) state , F (0)hist , we embed the input data I in temporal , concatenated with T emb , using two fully connected layers denoted by e state and e hist .During training, we minimize the mean squared error (MSE) between the ground truth future node quantities and the predicted quantities by TDE-GNN, similar to the training procedure of the rest of the considered methods in Table 1 . Specifically, following Rozemberczki et al. (2021b) , the goal is to predict the node quantities of the next time frame given 4 previous time frames.Algorithm 2 TDE-GNN for non-stationary problems with order o. Initialize history tensor Ho according to Equation (7).7:for l = 0 . . . L -1 do 8: The pendulum's motion in Example 1 is modeled by a time-varying frequency pendulum, such that q(F ; t) = sin(ω(t)F ), where ω(t) = 1 -0.04 sin(t).EQUATIONNode classification datasets. We report the statistics of the datasets used in our node classification experiments in Table 3 . All datasets are publicly available, and appropriate references to the data sources are provided in the main paper. 5 . Also, unless otherwise specified, in all experiments, we use L = 8 layers, and for node classification datasets, we use o = L. For the spatio-temporal datasets, we use o = r, i.e., the order is set to be equal to the number of historical data given by the task. In addition to the complexity discussion in the main paper, we provide the measured runtimes in Table 6 . Learning the temporal order and dynamics requires additional computations compared to the vanilla baseline of DE-GNN, however, it also offers improved performance as we show in our experiments in Section 4. We report the measured training and inferences runtimes, and the number of parameters on the Cora dataset in Table 6 . We measure the runtimes using an Nvidia-RTX3090 with 24GB of memory, which is the same GPU used to conduct our experiments. We now conduct an analysis of the learned coefficients c for the pendulum problem in Example 1.For convenience, we present the learned coefficients again, in Table 7 . Our analysis consists of two parts: (i) stability, and, (ii) consistency. Stability analysis. Following the derivations in Theorem 1 proof presented in Appendix E, we examine the root conditions of the learned coefficients. In Table 8 we report the absolute values of the characteristic polynomial with the coefficients from Table 7 . We note that for orders 2, 3, and 4, stability is obtained. However, for o = 5 the coefficients have one unstable mode. This result suggests that one should not use o = 5 for the pendulum problem in Example 1. However, orders lower than 5 yield stability. We believe that a stable fifth-order model is possible to be learned from the data, however, the incorporation of the root condition to the learning process requires adding constraints to the learning process and therefore is beyond the scope of this work. Table 8 : The absolute value of the roots r 1 , . . . , r 5 of the characteristic polynomial with coefficients from Table 7 .Consistency analysis. There are two conditions that verify the consistency of the learned coefficients.The first condition requires that the sum of the coefficients c equals to 1, which implies that if the spatial term in Equation ( 15), the future state (node features) is equal to a weighted average of the previous o -1 states. This implies that a constant solution can always be achieved. The second condition requires that the application of the coefficients c to a known discrete function yields a consistent approximation to its derivatives of some order.Note that condition (i) is satisfied by our construction of c. The second condition can be verified numerically, as we show now. To this end, we discretize the function y = sin(2πt) in the interval [0, 1]. We then apply the stencils based on the learned coefficients c as shown in Table 7 . As can be depicted in Figure 5 , the application of the stencils to the discrete function y(t) yields a scaled version of the second derivative of y(t) (which also equals to y, that is, ∂ 2 y(t) ∂t 2 = βy(t)). This result shows that our TDE-GNN is able to reveal the true order that describes the pendulum's motion, which is a second-order process.CONCLUSIONSIn this paper, we showed that incorporating higher-order interactions can be crucial to model data that arises from complex phenomena. This understanding motivated us to develop a novel architecture called TDE-GNN, that offers the modeling and understanding of higher-order ODE behaviors in GNNs, as well as improved downstream task performance.If the Jacobian has a positive real part, then the underlying ODE is unstable, and therefore its discretization is also unstable.In Algorithm 2, ⊕ denotes channel-wise concatenation. 
 
 Appendix:  


Decision: Excellent work that proposes TDE-GNNs for learning high order temporal dynamics using GNNs. TDE-GNNs model the temporal dynamics via neural ODEs, with maximal order >= 1 for modeling longer and more complex dependencies. The paper contains extensive details on the methods and is very well written.

Reviews: review: title: title: This paper proposes TDE-GNNs for learning high order temporal dynamics using GNNs. The results are promising, which supports the motivation of TDE-GNN and demonstrates its effectiveness.
main_review: main_review: This paper proposes TDE-GNNs for learning high order temporal dynamics using GNNs. TDE-GNNs model the temporal dynamics via neural ODEs, with maximal order >= 1 for modeling longer and more complex dependencies.

Pros
1. This paper is well written and easy to follow.
2. Extensive empirical studies on various datasets and comparisons to baselines demonstrate the effectivenss of proposed TDE-GNN.
3. The ablation on order $o$ shows that the performance consistently improve as $o$ increase, which well supports the motivation of TDE-GNN.
4. Implementation details are provided.

Cons
1. It seems that with such discretization in Eq(6), it is hard for TDE-GNN to incorporate advantages from neural ODEs, i.e., getting a continuous solution trajectory can be difficult.
title: title: Very nice, seminal work.
main_review: main_review: Just be mindful about some typos. 
For instance, a the end of the first paragraph in the introduction, "to model the of".