Graph Neural Ordinary Differential Equations-based method for Collaborative Filtering

Graph convolution network(GCN) is variant of Graph Neural Networks(GNNs). The basic idea of GCNs is to learn embeddings for nodes by smoothing features over the graph[1, 2, 3, 4]. For each layer, it aggregates the features of neighbors as the new embeddings of the target node. The propagation rule can be written as:

Where $AGG$ is aggregation function, $H_{n}$ is node embedding for the $n$th layer, $H_{n+1}$ is node embedding for $n+1$th layer, $W$ is the transformation matrix, and $A$ is adjacency matrix of the graph. However, there are studies[1] have shown that the GCN can only learn graph moments of a specific power of the Graph. To overcome this limitation, many GCN models[1, 11, 15] concatenate the node embeddings of each layer, which allows the model to capture information from different powers of $A$ and improves its ability to model complex graph structures.

Neural Ordinary Differential Equations(NODEs)[17] is a method that aims to model continuous dynamics on hidden state of neural networks. This is achieved by characterizing dynamics through an ordinary differential equation (ODE), which is parameterized by a neural network. The goal of this approach is to learn implicit differential equations from data. One of the key advantages of Neural ODEs is that they provide a flexible framework for modeling complex systems. By using ODEs, it is possible to model systems that exhibit continuous behavior. Additionally, by using neural networks to parameterize the ODEs, it is possible to capture complex patterns in the data that may be difficult to model using discrete methods. Neural ordinary differential equations calculate $h(t_{1})=h(t_{0})+\int_{t_{0}}^{t_{1}}f(h(t),\theta)dt$, where $f$ is a neural network parameterized by $\theta$ that approximates $\frac{dh(t)}{dt}$, to derive $h(t_{1})$ from $h(t_{0})$. $\theta$ is trained from data. The variable $t_{0}$ and $t_{1}$ are staring time and ending time, commonly set $t_{0}=0$ and then $t_{1}$ can be viewed as the number of layers in a neural networks, whereas it can be any non-negative real number in NODEs. In order to solve ODE problem, there are different ODE solvers that can be used to solve the ODE problem, such as Euler, RK2, RK4, and they can also be viewed as various neural network architectures. For example, the general form of the residual connection can be written as $h(t+1)=h(t)+f(h(t);\theta)$, which is same as the explicit Euler method for solving ODE problems. In general, the choice of ODE solver depends on the specific problem being solved. For instance, the explicit Euler method is simple and easy to implement, but it can be unstable. On the other hand, the fourth-order Runge-Kutta(RK4) method is more stable and accurate, but it can be computationally expensive[24, 19].

Deep learning has traditionally been dominated by discrete models. However, there are some studies[26, 17] have proposed treating neural networks as equipped with a continuum of layers, rather than a discrete set of layers. This approach enables a reformulation of the forward and backward pass as the solution to the initial value problem of an ODE. Such approaches allow direct modeling of ODEs and can guide discovery of novel general purpose deep learning models. This concept also build some connections between GCNs and traditional dynamical systems. Continuous Graph Neural Networks (CGNNs)[23] are an architecture that models continuous dynamics on node embeddings. First, it employ a neural encoder to project each node into a latent space based on its features. Afterwards, the encoded nodes are treated as the initial value $H(0)$, and an ODE is designed to define the continuous dynamics on node embeddings. This allows CGCNs to capture complex patterns of interactions between nodes over time. Finally, the node embeddings obtained at the ending time $t_{1}$ can be used for downstream applications such as node classification. Motivated by diffusion-based methods on graphs, an effective way for characterizing the propagation process is to use the following propagation equations:

where $H_{n}$ is the embeddings for all the nodes at stage $n$, and $H_{0}$ is the initial embeddings. Intuitively, each node at stage $n+1$ learns the node information from its neighbors through $AH_{n}$ and remembers its original node features through $H_{0}$. This way is to learn the graph structure without forgetting the original node features. And then, the formula can be derived as follows:

where the embeddings $H_{n}$ at stage $n$ incorporates all the information propagated up to $n$ steps with the initial embeddings $H_{0}$. With the idea of NODEs, we can replace the discrete step $n$ with a continuous variable $t$, and use an ODE to characterize such a continuous propagation dynamic. This allows us to naturally move from the discrete propagation process to the continuous case. Specifically, the embeddings $H_{t}$ at time $t$ incorporates all the information propagated up to time $t$ with the initial embeddings $H_{0}$. We can view this as a Riemann sum of an integral from time $0$ to time $t$, which leads us to the following proposition. To derive the equation using Taylor expansion, we have:

And,

where $H(0)$ represents the initial embeddings of all nodes, $H(t_{1})$ represents the node embeddings at time $t_{1}$ and $A$ represents the adjacency matrix of the graph. Typically, $t_{0}$ is set to 0 and $t_{1}$ can be treated as the number of layers. Therefore, $H(t_{1})$ can be interpreted as the node embeddings after $t_{1}$ convolutions.

Unlike Continuous Graph Neural Networks (CGNNs), which develop a continuous message-passing layer, Graph Neural Ordinary Differential Equations (GODEs)[24] parametrize the derivative function using 2 or 3 GCN layers directly. Based on graph spectral theory, the residual version of GCN layers[1, 27] are in the following form:

where $A$ is the graph Laplacian and $\sigma$ is a nonlinear activation function, $W$ is the transformation matrix, and $H_{n+1}$ is the node embedding at the $(n+1)$th layer. We denote the graph convolution operator as $G_{C}$. Where $G_{C}(H(n+1))=\sigma(AGG(W,A,H_{n}))$. The $n$-layer convolution can be written as follows:

Where $G_{C}$ is the graph convolution layer, $n$ is number of layers. And the GODEs to calculate the final node embeddings is:

where $H(t)$ represents the node embeddings at time $t$, $H(0)$ represents the initial embeddings, $GCN(n,H(t))$ represents the $n$ graph convolution layer with input node embeddings $H(t)$. $H(t_{1})$ can be treat as the final node embeddings that estimate by the ODE, with the input of $H(t)$.

There are various ODE solvers can solve the integral problem, and it can generalize various neural network architectures[17]. For instance, the explicit Euler method can be written as following step:

where $H(t)$ is node embeddings at state $t$. To estimate the solution of the differential equation, we can integrate the states up to $H(t_{1})$ using this formula, starting from a given initial value of $H(t_{0})$. The resulting $H(t)$ provide an approximation of the solution to the differential equation. Note that this equation is identical to a residual connection when $h=1$. Other ODE solvers like RK4 it is also can be viewed as dense convolution networks[17]. Other than Euler and RK4, there are many other solvers exist, such as dopri5, dopri8, fixed adams, midpoint and so on. For GODEs, higher-order ODE solvers are generally more effective, as long as the graph is dense enough to benefit from the additional computation [24]. In this study, we will only test two solvers: explicit Euler and RK4.

Collaborative filtering(CF)[7, 8, 9, 10, 11, 12, 13, 14] is a technique that can filter out items that a user might like on the basis of reactions by similar users. A typical CF paradigm involves learning user and item node embeddings, from historical interaction data, and give a top-$k$ recommendation based on the pairwise similarity of the user and item.

Let $E_{0}^{u}\in R^{N\times D}$ and $E_{0}^{i}\in R^{M\times D}$ be the initial user and item embeddings, respectively. Since user-item interactions can be represented as bipartite graphs, it is popular to adopt GCNs for CF. One of the early used GCN-based CF methods is NGCF[9], which utilizes non-linear activations and transformation matrices to map the user embedding space to the item embedding space. At each layer, user and item embeddings are combined by concatenation. However, non-linear activations and transformation matrices are not necessary for CF, as the user-item graph is sparse and non-linear GCNs are prone to overfitting. Moreover, empirical evidence suggests that transforming between user and item embedding spaces does not improve CF accuracy significantly. A simpler and more effective GCN-based CF method, called LightGCN[15], has been proposed. LightGCN removes the activation function and transformation matrix, and achieves state-of-the-art accuracy on many datasets. This suggests that linear GCNs with layer combination perform best among various design choices. The definition of the linear graph convolutional layer is as follows:

Here, $A$ represents the normalized adjacency matrix, and $E_{k}$ represents the user/item embedding matrix at the $k$-th layer. The initial embeddings, which are denoted as $E^{u}_{0}$ and $E^{p}_{0}$. The final embeddings $E^{u}_{f}$ and $E^{p}_{f}$ are calculated using layer combination, expressed as follows:

where $K$ is the number of layers, $w_{l}$ is weight for $l$th layer, and $E^{u}_{f}$ and $E^{i}_{f}$ are the final node embeddings. All CF methods learn the initial embeddings of users and items. The model prediction is defined as the inner product of user and item final embeddings:

where $E_{f}^{u}$ is the final users embeddings and $E_{f}^{i}$ is the final items embeddings. Typically, the following Bayesian personalized ranking (BPR)[10] loss is used to train the initial embedding vectors in the field of collaborative filtering.

The aggregation of messages in LightGCN, a linear GCN model, can be viewed as the integration of an ordinary differential equation (ODE). Building on this concept, LT-OCF[19] was developed as an ODE-based method. The complete model can be viewed as a continuous version of LightGCN with a learnable number of layers. Figure 1 offers a comparison of the two methods. The message passing of LT-OCF for single ODE can be written as follows:

where, $u(t)\in R^{N\times D}$ refers to user embeddings, and $i(t)\in R^{M\times D}$ refers to item embeddings at state $t$. $K$ can be treated as the number of layers. $u(K)$ and $i(K)$ are the user embeddings and item embeddings at time $K$, respectively. For the whole model, the time interval $[0,K]$ is split into $T$ parts. For instance, if $T=2$, the model will have $2$ ODE blocks. The first ODE block at time $[0,t_{1}]$ and the second ODE block at time $[t_{1},K]$. The final embedding will be a combination of the output of these two blocks, the initial embeddings, and the embeddings at time $K$. The final embeddings are calculated as follows:

where $u(0)$ and $i(0)$ are initial user embeddings and initial item embeddings respectively, $w_{i}$ is weight for each ODE block, and $K$ is the ending time or number of layers. $E^{u}_{f}$ and $E^{i}_{f}$ are final user embeddings and final item embeddings, respectively.

As we previously mentioned, the linear GCN model can be interpreted as an ODE problem [19, 17]. This means that we can view the GCN process as a continuous flow of information across the graph, where the state of each node changes continuously over time, governed by a differential equation. This concept has led to the development of Graph Neural Ordinary Differential Equations (GODEs), which extend the GCN framework to a continuous-time setting. With GODEs, we can estimate the final node embedding with one or two GCN layers. Specifically, by solving the ODEs governing the GCN process, we can obtain a continuous-time dynamical system that describes the evolution of the node embeddings over time. This can be viewed as a generalization of the static GCN model, where we can capture more complex temporal dynamics in the graph. In addition, we can also use the GODE framework to incorporate other types of information, such as node attributes or edge features. Inspired by the success of GODE in modeling continuous-time dynamics, we believe that this framework is also suitable for the collaborative filtering (CF) task. In CF, the goal is to estimate user and item embeddings. By leveraging GODEs, we can estimate these embeddings more accurately.

We propose a new method based on the idea of Graph Neural ODE, as shown in Fig 1. Instead of creating a continuous message passing layer, we parametrize the derivative function using multiple layers of GCN directly. In other words, we use the information captured by one or two GCN layers to estimate the final state of the embedding by solving an ODE problem without creating all of the layers. Unlike the concept of CGNN[18], we treat $t_{1}$ as hyperparameters, and our model doesn’t have any layer combinations, since integration can be viewed as the summation of all layers from time $0$ to $t_{1}$. The initial embeddings are directly fed into the ODEs, and the output of the ODEs becomes the final embedding. The overall formula can be written as:

In this model, $A$ is the adjacency matrix, $E^{u}_{0}$ and $E^{i}_{0}$ are the initial user and item embeddings, respectively. $W$ is a trainable weight and $n$ is the number of layers. The user and item embeddings at time $t$ are denoted as $E^{u}(t)$ and $E^{i}(t)$, respectively. The final user embeddings and item embeddings are denoted as $E^{u}_{f}$ and $E^{i}_{f}$, respectively. Unlike other GCN-based models [11, 15], where the final embeddings are a combination of all the layers. Our method aims to estimate the final embeddings using the information that captured by several GCN layers through an ODE function. Same as other methods, we use BPR[10] loss to train the embeddings and use ODE solvers, such as explicit Euler and RK4, to solve the ODE problem.

When it comes to solving ODEs, there are many options available, each with their own strengths and weaknesses. In fact, different ODE solvers may produce vastly different results even when applied to the same dataset. As a result, we believe that it is crucial to carefully consider the choice of ODE solver based on the characteristics of the dataset being used. By doing so, we can ensure that we are using an appropriate solver that will provide accurate and reliable results.

And the model prediction is defined as the inner product of the final embeddings of the user and item:

where $y_{u,i}$ is the rating matrix for all users.

Suppose the following setting in our method:i) the number of layers $n$ is $1$, ii) We use euler’s method to solve the ODE problem, iii) we set $W$ to $1$. Then our model can be written as follows:

which is same as the Continuous Graph Neural Networks(CGCN), or a LT-OCF with single ODE block ($K=1$).

In this section, we introduce our experimental environments and results. All experiments were conducted in the following hardware environments: 12th Gen Intel(R) Core(TM) i9-12900F, with one GeForce RTX 3080 GPU. We utilize the torchdiffeq PyTorch package to solve and backpropagate through the ODE solver.

Datasets. We use the public Amazon Reviews dataset[28] with four benchmark categories, including: Beauty, Health, Cell Phones, Office Product. The details of the datasets are summarized in Table I. We follow the 5-core setting as existing works on users and the same transformation[11, 15, 29, 30] of treating the existence of reviews as positives. We sort each user’s interactions chronologically and adopt the leave-one-out setting, with the last interacted item for testing and the second last interaction for validation.

Eveluation. We use the standard ranking evaluation metrics, including Recall@N and NDCG@N, to evaluate the averaged ranking performance over all users. These metrics are widely used in existing works, such as [11, 15, 29, 19]. Recall@N measures the correctness of ranking the ground-truth items in the top-N list, regardless of the ranking position. It is a widely used metric for evaluating the effectiveness of top-N recommendation algorithms. NDCG@N extends the idea of Recall@N by giving higher weights to higher ranking positions. It measures the quality of a ranking algorithm based on the graded relevance of the items in the ranked list. We only report the N=20 due to limited space.

Baselines. In total, we compare GODE-CF with various types of the state-of-the-art models:

• •

  NGCF[11] is a GCN-based CF method. It uses feature transformation and non-linear activations

• •

  LightGCN[15] are linear GCN-based CF methods.

• •

  UltraGCN[25] is an ultra-simplified formulation of GCN which skips explicit message passing and directly approximate the limit of infinite message passing layers.

• •

  layerGCN[31] is a layer-refined GCN model, that refines layer representations during information propagation and node updating of GCN.

• •

  GTN[29] a graph trend filtering networks framework to capture the adaptive reliability of the interactions between users and items. Extensive experiments on four real-world datasets demonstrate the.

• •

  LT-OCF[19] is a NODE-based method that learn the optimal architecture of the model for CF.

In this section we introduce the grid search of the hyperparameters. The learning rate is $0.001$, and the embedding size is $128$ for all models and the Normal distribution of $N(0,1)$ is used to set initial embeddings.

• •

  For LT-OCF, the regularization coefficient $\lambda$ in all method is in $\{1e-4,1e-3,1e-2\}$. The number of elements $K$ is in $\{2,3,4\}$. The number of learnable intermediate time points T is in $\{1,2,3\}$, we consider the following solvers: the explicit Euler, and RK4.

• •

  For NGCF, we search the node and message dropouts from $\{0.1,0.3,0.5\}$ and number of layers from $\{1,2,3,4,5\}$.

• •

  For GTN, we search parameters $\beta$ and $\gamma$ from $\{0.001,0.005,0.01,0.05,0.1\}$.

• •

  For UltraGCN, we search their weights from $\{1,1e-3,1e-5,1e-7\}$ and negative weights from $\{10,50,100\}$.

• •

  For layerGCN, we search number of layers from $\{1,2,3,4,5\}$

• •

  For LightGCN, we search the node and message dropouts from $\{0.1,0.3,0.5\}$ and number of layers from $\{1,2,3,4,5\}$

• •

  For GODE-CF, we search number of layers from $\{1,2,3\}$, and $t$ from $\{0.7,0.75,0.8,0.85,0.9,0.95,1\}$

In Table II, we summarize the overall accuracy in terms of Recall and NDCG.

• •

  Among GCN-based methods, GODE-CF consistently achieves the best performance across all four datasets. In particular, on Health, GODE-CF significantly outperforms the strongest GCN-based baseline, GTN, by 11% on Recall@20. For other baselines, LightGCN exhibits the best performance on Beauty, while GTN demonstrates the best performance on Health compared to other GCN-based baselines.

• •

  Among all baselines, the ODE-based model LT-OCF shows state-of-the-art accuracy for Beauty. However, our method, GODE-CF, achieves the best accuracy in all cases. The Euler solver yields the best results for all four datasets compared to RK4. We believe that the choice of solver may affect performance depending on the dataset.

We compare the training curve of LightGCN and GODE-CF for Beauty and Health. In general, our method provides a faster training speed in terms of the number of epochs than LightGCN.

In this section, we provide empirical evidence demonstrating the superiority of GODE-CF’s training efficiency compared to other baselines. We believe that it is essential to prove this convincingly, so we compare the total training time and epochs required to achieve the best performance among the models. Our experiment includes validation time in the training time to provide a comprehensive overview of the training process. To ensure a fair comparison, we use a uniform code framework that we implemented ourselves for all models. We train all models with a fixed number of 1000 epochs to eliminate the effect of varying epoch numbers. Additionally, we conduct our experiments on all four datasets using the same machine:12th Gen Intel(R) Core(TM) i9-12900F, with one GeForce RTX 3080 GPU for all compared models.

Table III shows the details of the results. The explicit Euler providing the best performance for our model for all four datasets. Since this solver is simple and fast, which is a contributing factor to GODE-CF’s short training time. By comparing the total training time and epochs required to achieve the best performance, we demonstrate that GODE-CF consistently outperforms other baselines.

In this section, we compare two different solvers for ordinary differential equations (ODEs): Euler and RK4. The Fig  3 and Fig 4 are the training curves of these two solvers. Other solvers, such as DOPRI, are almost the same as RK4, but RK4 has smaller computation complexity. Therefore, we only use RK4 and Euler’s method for this experiment. Contrary to what we thought, Euler’s method shows better accuracy for almost all cases, whereas RK4 should have better accuracy in terms of solving general ODE problems. For instance, the Euler method achieves a Recall/NDCG of 0.050791/0.019095 vs. 0.048531/0.018325 with RK4 on Cell Phones. However, for other datasets, we observe very similar patterns where Euler’s method consistently outperforms RK4. And Euler’s method provides a faster training speed in terms of the number of epochs than that of RK4 for all of the cases. This is because using the RK4 solver is more likely to result in overfitting the training data in this model.

We also conducted experiments that explored the impact of varying the number of layers inside ODE on the model’s performance. The detailed results are in Table IV.

Based on the results, we have found that a 2-layer design yields optimal results for all four cases. This is because the ODE needs to estimate the final state of the embeddings by using enough information captured by the GCN layer. In other words, the GCN layer need to provides the necessary information for the ODE to make optimal estimations.

On the other hand, a 1-layer design may not provide enough information for the ODE to make a good estimation. This is because one GCN layer may not capture all the necessary information required to estimate the final state of the embeddings.

Conversely, a 3-layer design may capture excessive information, leading to overfitting of the training data. To address this issue, we suggest incorporating the dropout technique when using a 3-layer design. Dropout can help prevent overfitting and ensure that the final state estimated by the ODE is more reliable.

By varying $t$, we also investigate how different $t$ affect the model accuracy. The detailed results are shown in Figure  5.

According to the plot, we can see that the model with only a small number of layers requires a larger value of $t$ to reach the optimal performance. Conversely, models with larger numbers of layers require a smaller value of $t$ to achieve the optimal performance. This is because models with fewer layers provide less information to the ODE, which results in the requirement of a longer time step to reach the final state. Conversely, models with a higher number of layers capture more information, which subsequently results in them being closer to the final state, thereby allowing the ODE to take smaller time steps $t$. It is important to note that the number of layers plays a crucial role in the performance of the model. An interesting finding is that, for both 2-layer and 3-layer settings, there is a time point at which the model experiences a significant performance drop. However, after this time point, the performance begins to recover. The cause of this issue remains unknown, so we will leave it for future work.

We are model has a learneable weight for each layer inside ODE. We are also interested in to investigate whether the weight has a significant impact on the model’s performance. The detailed results are shown in Table V. Based on the results, we found that the impact of weight is highly dependent on the dataset. For most cases, the model without weight has better performance than the model with weight on the Beauty, but the difference is not significant. However, for the Office, the model has a larger improvement when using weight. Whether or not to use weight depends highly on the specific case. It’s recommended to evaluate the impact of weight on the model for each specific dataset.