Continuous Graph Neural Networks

Introduction

Graph neural networks (GNNs) have been attracting growing interest due to their simplicity and effectiveness in a variety of applications such as node classification , link prediction , chemical properties prediction , and natural language understanding . The essential idea of GNNs is to design multiple graph propagation layers to iteratively update each node representation by aggregating the node representations from their neighbours and the representation of the node itself. In practice, a few layers (two or three) are usually sufficient for most tasks , and more layers may lead to inferior performance .

A key avenue to improving GNNs is being able to build deeper networks to learn more complex relationships between the data and the output labels. The GCN propagation layer smooths the node representations, i.e. nodes near each other in the graph become more similar . This can lead to over-smoothing as we stack more and more layers, meaning that nodes representations converge to the same value leading to worse performance . Thus, it is important to alleviate the node over-smoothing effect, whereby node representations converge to the same value.

Furthermore, it is crucial to improve our theoretical understanding of GNNs to enable us to characterise what signals from the graph structure we can learn. Recent work on understanding GCN has considered GCN as a discrete dynamical system defined by the discrete layers. In addition, Chen et al. demonstrated that the usage of discrete layers is not the only perspective to build neural networks. They pointed out that discrete layers with residual connections can be viewed as a discretisation of a continuous Ordinary Differential Equation (ODE). They showed that this approach is more memory-efficient and is able to model the dynamics of hidden layers more smoothly.

We use the continuous perspective inspired by diffusion based methods to propose a new propagation scheme, which we analyse using tools from ordinary differential equations (i.e. continuous dynamical systems). Indeed we are able to explain what representations our model learns as well as why it does not suffer from the common over-smoothing problem observed in GNNs. Allowing us to build ‘deeper’ networks in the sense that our model works well with large values of time. The key factor for the resilience to over-smoothing is the use of a restart distribution as originally proposed in Pagerank in a continuous setting. The intuition is that the restart distribution helps to not forget the information from low powers of the adjacency matrix, hence enabling the model to converge towards a meaningful stationary distribution.

The main contributions of this paper are:

We propose two continuous ODEs of increasing model capacity inspired by PageRank and diffusion based methods;

We theoretically analyse the representation learned by our layer and show that as $t\to\infty$ our method approaches a stable fixed point, which captures the graph structure together with the original node features. Because we are stable as $t\to\infty$ our network can have an infinite number of ‘layers’ and is able to learn long-range dependencies;

We demonstrate that our model is memory efficient and is robust to the choice of $t$. Further, we demonstrate on the node classification against competitive baselines that our model is able to outperform many existing state-of-the-art methods.

Preliminaries

Let a graph $G=(V,E)$Throughout the paper we will only consider simple graphs be defined by vertices $V$ and edges $E\subseteq V\times V$ between vertices in $V$. It is common in graph machine learning to use an adjacency matrix $\boldsymbol{Adj}$ as an alternative characterisation. Given a node ordering $\pi$ and a graph $G$, the elements of the adjacency matrix $\boldsymbol{Adj}^{|V|\times|V|}$ are defined by the edge set $E$:

As the degree of nodes can be very different, we typically normalize the adjacency matrix as $\boldsymbol{D}^{-\frac{1}{2}}\boldsymbol{Adj}\boldsymbol{D}^{-\frac{1}{2}}$, where $\boldsymbol{D}$ is the degree matrix of $\boldsymbol{Adj}$. Such a normalized adjacency matrix always has an eigenvalue decomposition, and the eigenvalues are in the interval $$ . The negative eigenvalues can make graph learning algorithms unstable in practice, and hence we follow Kipf and Welling and leverage the following regularized matrix for characterizing graph structures:

where $\alpha\in(0,1)$ is a hyperparameter, and the eigenvalues of $\boldsymbol{A}$ are in the interval $[0,\alpha]$.

A continuous ODE throughout this paper will refer to an equation of the following form:

where $x$ may be scaler, vector valued, or matrix valued and $f$ is function that we will parametrise to define the hidden dynamic. Chen et al. showed how we can backpropagate through such an ODE equation and hence use it as a building block for a neural network.

Related Work

Neural ODE Neural ODE is an approach for modelling a continuous dynamics on hidden representation, where the dynamic is characterised through an ODE parameterised by a neural network. However, these methods can only deal with unstructured data, where different inputs are independent. Our approach extends this in a novel way to graph structured data.

GNNs. Graph neural networks are an effective approach for learning node representations in graphs. Typically, GNNs model discrete dynamics of node representations with multiple propagation layers, where in each layer the representation of each node is updated according to messages from neighbouring nodes. The majority of GNNs learn the relevant information from the graph structure $\boldsymbol{A}$ by learning finite polynomial filters $g$ to apply to the eigenvalues $\boldsymbol{\Lambda}$ of the graph Laplacian $\boldsymbol{L}=\boldsymbol{P}\boldsymbol{\Lambda}\boldsymbol{P}^{-1}$ in each propagation layer . GCN , for instance, uses a first-order Chebyshev polynomial. However, existing GNNs (e.g. GCN) have been shown to suffer from over-smoothing, i.e. node representations start to converge to the same value. Compared to these studies, we follow continuous dynamical systems to model the continuous dynamic on node representations. Moreover, our approach does not have the over-smoothing problem, because our model converges to a meaningful representation (fixed point) as $t\to\infty$.

A proposed solution in recent work is to either add a residual connection to the preceding layer or to use a concatenation of each layer . The latter approach does not scale to very deep networks due to the growing size of the representation. Chen et al. ’s approach of turning residual connections from the preceding layer into a continuous ODE has gathered much attention . In contrast to many of these works we are able to provide a theoretical justification for our ODE and our representation does not grow with depth.

A significant amount of work has focused on understanding the theoretical properties of GCN and related architectures better . Dehmamy et al. argue that to learn the topology of a graph the network must learn the graph moments—an ensemble average of a polynomial of $\boldsymbol{Adj}$. They show that GCN can only learn graph moments of a specific power of $\boldsymbol{Adj}$. To remedy this they concatenate $[\boldsymbol{h}^{(1)},\boldsymbol{h}^{(2)},\ldots,\boldsymbol{h}^{(n)}]$ the feature maps (output) of each layer. Oono and Suzuki demonstrate that deeper GCNs exponentially lose expressive power by showing that in the limit as the number of layers goes to infinity the node representation is projected onto the nullspace of the Laplacian, which implies that two nodes with identical degree in the same connected component are will have the same representation. Instead we propose a continuous dynamical system which does not suffer from this problem and demonstrate so as $t\to\infty$. NT and Maehara focus on explaining which assumptions hold such that GCN and SGC work on the common citation networks. Our work fits into the existing literature by using the Neural ODE framework to provide a novel intuition to what is needed to address loss of expressive power with depth. We do this by proposing and analysing our specific solution.

Concurrent work. Several concurrent works have developped similar ideas; proposes an ODE based on treating the GCN-layer as a continuous vector field and combines discrete and continuous layers. Other concurrent works use the Neural ODE framework and parametrise the derivative function using a 2- or 3-layer GNN directly, instead we develop a continuous message-passing layer and do not use a discrete deep neural network to parametrise the derivative. In contrast, we motivate our ODE from diffusion based methods and theoretically justify how our approach helps to address over-smoothing as well as long-range connections between nodes.

Model

The key step of the framework is to design an effective ODE for defining the continuous dynamics on node representations, and thereby modelling the dependency of nodes. We design two such ODEs of increasing model capacity based on the intuition from existing diffusion-based methods on graphs. In the first ODE, each feature channel (i.e. dimension) of node representations evolves independently (see Sec. 4.1), whereas in the second ODE we also model the interaction of different feature channels (see Sec. 4.2).

Since different nodes in a graph are interconnected, a desirable ODE should take the graph structure into consideration, and allow information to propagate among different nodes. Motivated by existing diffusion-based methods on graphs (e.g. PageRank and label propagation ), an effective way for characterizing the propagation process is to use the following step-wise propagation equations:

where we see that the representation $\boldsymbol{H}_{n}$ at step $n$ incorporates all the information propagated up to $n$ steps with the initial representation $\boldsymbol{H}_{0}$.

The discrete dynamic in Eq. (5), where $\boldsymbol{A}$ has an eigenvalue decomposition, is a discretisation of the following ODE:

with the initial value $\boldsymbol{H}(0)=(\ln\boldsymbol{A})^{-1}(\boldsymbol{A}-\boldsymbol{I})\boldsymbol{E}$, where $\boldsymbol{E}=\mathcal{E}(\boldsymbol{X})$ is the output of the encoder $\mathcal{E}$.

We provide the proof in the Supplementary material. In practice, $\ln\boldsymbol{A}$ in Eq. (6) is intractable to compute, hence we approximate it using the first order of the Taylor expansion, i.e. $\ln\boldsymbol{A}\approx(\boldsymbol{A}-\boldsymbol{I})$, which gives us:

with the initial value being $\boldsymbol{H}(0)=\boldsymbol{E}$. This is the ODE we use in our model CGNN. The intuition behind the ODE defined in Eq. 7 can be understood from an epidemic modelling perspective. The epidemic model aims at studying the dynamics of infection in a population. Typically, the model assumes that the infection of people is affected by three factors, i.e. the infection from neighbours, the natural recovery, and the natural physique of people. Suppose that we treat the latent vectors $\boldsymbol{H}(t)$ as the infection conditions of a group of people at time $t$, then the three factors can be naturally captured by three terms: $\boldsymbol{A}\boldsymbol{H}(t)$ for the infection from neighbours, $-\boldsymbol{H}(t)$ for natural recovery, and $\boldsymbol{E}$ for the natural physique. Therefore, the infection dynamics in a population can be intuitively modelled by our first-order ODE, in Eq. (6), indicating that the intuition of our ODE agrees with the epidemic model.

The ODE we use can be understood theoretically. Specifically, the node representation matrix $\boldsymbol{H}(t)$ at time $t$ has an analytical form, which is formally stated in the following proposition.

The analytical solution of the ODE defined in Eq. (7) is given by:

We prove the proposition in the Supplementary material. From the proposition, since the eigenvalues of $\boldsymbol{A}-\boldsymbol{I}$ are in the interval $[-1,0)$, as we increase $t$ to $\infty$, the exponential term $e^{(\boldsymbol{A}-\boldsymbol{I})t}$ will approach $\boldsymbol{0}$, i.e. $\lim_{t\to\infty}e^{(\boldsymbol{A}-\boldsymbol{I})t}=\boldsymbol{0}$. Therefore, for large enough $t$ we can approximate $\boldsymbol{H}(t)$ as:

Thus, $\boldsymbol{H}(t)$ can be seen as the summation of all different orders of propagated information (i.e. $\{\boldsymbol{A}^{i}\boldsymbol{E}\}_{i=1}^{\infty}$). In this way, our approach essentially has an infinite number of discrete propagation layers, allowing us to model global dependencies of nodes more effectively than existing GNNs.

Implementation. In the node classification task, our decoder $\mathcal{D}$ to compute the node-label matrix $\boldsymbol{Y}=\mathcal{D}(\boldsymbol{H}(t_{1}))$ is a softmax classifier with the ReLU activation function .

Note the parameter $\alpha$ in Eq. (2) decides the eigenvalues of $\boldsymbol{A}$, which thereby determines how quickly the higher order powers of $\boldsymbol{A}$ go to , this also means that by specifying $\alpha$ per node we can control how much of the neighbourhood each node gets to see as smaller values $\alpha$ imply that the powers of $\boldsymbol{A}$ vanish faster. In our final model we learn these parameters $\alpha$.

2 Case 2: Modelling the Interaction of Feature Channels

The ODE so far models different feature channels (i.e. dimensions of hidden representations) independently, where different channels are not able to interact with each other, and thus the ODE may not capture the correct dynamics of the graph. To allow the interaction between different feature channels, we are inspired by the success of a linear variant of GCN (i.e. Simple GCN ) and consider a more powerful discrete dynamic:

Similar to the previous section, we extend the discrete propagation process in Eq. (10) to continuous cases by viewing each $\boldsymbol{H}_{n}$ as a Riemann sum of an integral from time to time $t=n$, which yields the following proposition:

Suppose that the eigenvalue decompositions of $\boldsymbol{A},\boldsymbol{W}$ are $\boldsymbol{A}=\boldsymbol{P}\boldsymbol{\Lambda}\boldsymbol{P}^{-1}$ and $\boldsymbol{W}=\boldsymbol{Q}\boldsymbol{\Phi}\boldsymbol{Q}^{-1}$, respectively, then the discrete dynamic in Eq. (10) is a discretisation of the following ODE:

where $\boldsymbol{E}=\mathcal{E}(\boldsymbol{X})$ is the output of the encoder $\mathcal{E}$ and with the initial value $\boldsymbol{H}(0)=\boldsymbol{P}\boldsymbol{F}\boldsymbol{Q}^{-1}$, where

where $\widetilde{\boldsymbol{E}}=\boldsymbol{P}^{-1}\boldsymbol{E}\boldsymbol{Q}$.

The proof is provided in the Supplementary material. By making a first-order Taylor approximation to get rid of the matrix logarithm, we further obtain:

with the initial value being $\boldsymbol{H}(0)=\boldsymbol{E}$. This is the ODE we use in our model CGNN with weights. The ODEs of the form as in Eq. (13) have been studied in some detail in the control theory literature, where they are known as the Sylvester differential equation . Intuitively, $\boldsymbol{E}$ here would be the input into the system with the goal to get the system $\boldsymbol{H}$ into a desired state $\boldsymbol{H}(t)$. The matrices $\boldsymbol{A}-\boldsymbol{I}$ and $\boldsymbol{W}-\boldsymbol{I}$ describe the natural evolution of the system.

The ODE in Eq. (13) also has appealing theoretical properties. Specifically, $\boldsymbol{H}(t)$ has an analytical form as shown in the following proposition.

Suppose that the eigenvalue decompositions of $\boldsymbol{A}-\boldsymbol{I},\boldsymbol{W}-\boldsymbol{I}$ are $\boldsymbol{A}-\boldsymbol{I}=\boldsymbol{P}\boldsymbol{\Lambda}^{\prime}\boldsymbol{P}^{-1}$ and $\boldsymbol{W}-\boldsymbol{I}=\boldsymbol{Q}\boldsymbol{\Phi}^{\prime}\boldsymbol{Q}^{-1}$, respectively, then the analytical soluton of the ODE in Eq. (13) is given by:

where $\widetilde{\boldsymbol{E}}=\boldsymbol{P}^{-1}\boldsymbol{E}\boldsymbol{Q}$.

We prove the proposition in the Supplementary material. According to the definition of $\boldsymbol{A}$ and also our assumption about $\boldsymbol{W}$, the eigenvalues of $\boldsymbol{A}-\boldsymbol{I}$ and $\boldsymbol{W}-\boldsymbol{I}$ are in $(-1,0)$, and therefore $\Lambda_{i,i}^{\prime}<0$ for every $i$ and $\Phi_{j,j}^{\prime}<0$ for every $j$. Hence, as we increase $t$ to $\infty$, the exponential terms will approach 0, and hence for large enough $t$ we can approximate $\boldsymbol{H}(t)$ as:

Based on the above results, if $\boldsymbol{W}=\boldsymbol{I}$, then $\boldsymbol{H}(t)$ will converge to the same result as in Eq. (9), and hence the ODE defined in Eq. (6) is a special case of the ODE in Eq. (11).

Implementation. We use the same decoder as for the case when $W=I$.

where $\beta$ is a hyperparameter, and the above secondary update enables $\boldsymbol{U}$ to be close to the manifold of orthogonal matrices after each training step.

Finally, to help stabilise training we use the idea from and add auxiliary dimensions to hidden representation only during the continuous propagation process. Specifically, we double the latent representation initialising the second half of the initial representation with 0 and throwing the result away after solving the continuous ODE. This very slightly improves results, but importantly stabilises training significantly (see Supplementary material).

Discussion

Our continuous model for information propagation has several advantages over previous discrete GNNs such as GCN:

Learning global dependencies in the graph;

$\alpha$ represents the “diffusion” constant, which is learned;

Weights entangle channels continuously over time;

Insight into the role of the restart distribution $\boldsymbol{H}_{0}$.

1. Robustness with time to over-smoothing: Despite the effectiveness of the discrete propagation process, it has been shown in that the usage of discrete GCN layers can be tricky as the number $n$ of layers (time in the continuous case) is a critical choice. Theoretical work in further showed that on dense graphs as the number of GCN layers grow there is exponential information loss in the node representations. In contrast, our method is experimentally not very sensitive to the integration time chosen and theoretically does not suffer from information loss as time goes to infinity.

2. Global dependencies: Recent work has shown that to improve on GNN it is necessary to build deeper networks to be able to learn long-range dependencies between nodes. Our work, thanks to the stability with time is able to learn global dependencies between nodes in the graph. Eq. (9) demonstrates that we propagate the information from all powers of the adjacency matrix, thus we are able to learn global dependencies.

3. Diffusion constant: The parameter $\alpha$ scales the matrix $\boldsymbol{A}$ (see Eq. (2)), i.e. it controls the rate of diffusion. Hence, $\alpha$ controls the rate at which higher-order powers of $\boldsymbol{A}$ vanish. Since each node has its own parameter $\alpha$ that is learned, our model is able to control the diffusion, i.e. the weight of higher-order powers, for each node independently.

4. Entangling channels during graph propagation: The ODE with weights (Eq. (11)) allows the model to entangle the information from different channels over time. In addition, we are able to explain how the eigenvalues of the weight matrix affect the learned representation (Eq. (14)).

5. Insight into the role of the restart distribution $\boldsymbol{H}_{0}$: In both of our ODEs, Eq. (7) and (13), the derivative depends on $\boldsymbol{E}$, which equals to the initial value $\boldsymbol{H}(0)$. To intuitively understand the effect of the initial value in our ODEs, consider an ODE without $\boldsymbol{E}$, i.e. $\boldsymbol{H}^{\prime}(t)=(\boldsymbol{A}-I)\boldsymbol{H}(t)$. The analytical solution to the ODE $\boldsymbol{H}^{\prime}(t)=(\boldsymbol{A}-I)\boldsymbol{H}(t)$ is given by $\boldsymbol{H}(t)=\exp[(\boldsymbol{A}-\boldsymbol{I})t]\boldsymbol{H}(0)$. Remembering that $\boldsymbol{A}-\boldsymbol{I}$ is simply a first order approximation of $\ln\boldsymbol{A}$, we can see that the analytical solution we are trying to approximate is $\boldsymbol{H}(t)=\boldsymbol{A}^{t}\boldsymbol{H}(0)$. Thus, the end time of the ODE now determines, which power of the $\boldsymbol{Adj}$ we learn. Indeed in our experiments we show that removing the term $\boldsymbol{H}(0)$ causes us to become very sensitive to the end time chosen rather than just needing a sufficiently large value (see Fig. 2).

Experiment

In this section, we evaluate the performance of our proposed approach on the semi-supervised node classification task.

In our experiment, we use four benchmark datasets for evaluation, including Cora, Citeseer, Pubmed, and NELL. Following existing studies , we use the standard data splits from for Cora, Citeseer and Pubmed, where 20 nodes of each class are used for training and another 500 labeled nodes are used for validation. For the NELL dataset, as the data split used in is not available, we create a new split for experiment. The results are in Table 2. We further run experiments with random splits on the same datasets in Table 3. The statistics of the datasets are summarized in Table 1. Accuracy is used as the evaluation metric.

2 Compared Algorithms

Discrete GNNs: For standard graph neural networks which model the discrete dynamic of node representations, we mainly compare with the Graph Convolutional Network (GCN) and the Graph Attention Network (GAT) , which are the most representative methods.

Continuous GNNs: There is also a recent concurrent work which learns node representations through modelling the continuous dynamics of node representations, where the ODE is parameterised by a graph neural network. We also compare with this method (GODE).

CGNN: For our proposed Continuous Graph Neural Network (CGNN), we consider a few variants. Specifically, CGNN leverages the ODE in Eq. (7) to define the continuous dynamic of node representations, where different feature channels are independent. CGNN with weight uses the ODE in Eq. (13), which allows different feature channels to interact with each other. We also compare with CGNN discrete, which uses the discrete propagation process defined in Eq. (4) for node representation learning (with $n=50$).

3 Parameter Settings

We do a random hyperparameter search using the Orion framework with 50 retries. The mean accuracy over 10 runs is reported for each dataset in Table 2.

4 Results

1. Comparison with existing methods. The main results of different compared algorithms are summarized in Table 2. Compared with standard discrete graph neural networks, such as GCN and GAT, our approach achieves significantly better results in most cases. The reason is that our approach can better capture the long-term dependency of different nodes. Besides, our approach also outperforms the concurrent work GODE on Cora and Pubmed. This is because our ODEs are designed based on our prior knowledge about information propagation in graphs, whereas GODE parameterises the ODE by straightforwardly using an existing graph neural network (e.g. GCN or GAT), which may not effectively learn to propagate information in graphs. Overall, our approach achieves comparable results to state-of-the-art graph neural networks on node classification.

2. Comparison of CGNN and its variants. The ODEs in CGNN are inspired by the discrete propagation process in Eq. (4), which can be directly used for modelling the dynamic on node representations. Compared with this variant (CGNN discrete), CGNN achieves much better results on all the datasets, showing that modelling the dynamic on nodes continuously is more effective for node representation learning. Furthermore, comparing the ODEs with or without modelling the interactions of feature channels (CGNN with weight and CGNN respectively), we see that their results are close. A possible reason is that the datasets used in experiments are quite easy, and thus modelling the interactions of feature channels (CGNN with weight) does not bring much gain in performance. We anticipate CGNN with weight could be more effective on more challenging graphs, and we leave it as future work to verify this.

3. Performance with respect to time steps. One major advantage of CGNN over existing methods is that it is robust to the over-smoothing problem. Next, we systematically justify this point by presenting the performance of different methods under different numbers of layers (e.g. GCN and GAT) or the ending time $t$ (e.g. CGNN and its variants).

The results on Cora and Pubmed are presented in Fig. 2. For GCN and GAT, the optimal results are achieved when the number of layers is 2 or 3. If we stack more layers, the results drop significantly due to the over-smoothing problem. Therefore, GCN and GAT are only able to leverage information within 3 steps for each node to learn the representation. In contrast to them, the performance of CGNN is more stable and the optimal results are achieved when $t>10$, which shows that CGNN is robust to over-smoothing and can effectively model long-term dependencies of nodes. To demonstrate this we use the ODE $\boldsymbol{H}^{\prime}(t)=(\boldsymbol{A}-\boldsymbol{I})\boldsymbol{H}(0)$ with $\boldsymbol{H}(0)=\boldsymbol{E}$, which gets much worse results (CGNN w/o $\boldsymbol{H}(0)$), showing the importance of the initial value for modelling the continuous dynamic. Finally, CGNN also outperforms the variant which directly models the discrete dynamic in Eq. (4) (CGNN discrete), which demonstrates the advantage of our continuous approach.

4. Memory Efficiency. Finally, we compare the memory efficiency of different methods on Cora and Pubmed in Fig. 3. For all the methods modelling the discrete dynamic, i.e. GCN, GAT, and CGNN discrete, the memory cost is linear to the number of discrete propagation layers. In contrast, through using the adjoint method for optimization, CGNN has a constant memory cost and the cost is quite small, which is hence able to model long-term node dependency on large graphs.

Conclusion

In this paper, we build the connection between recent graph neural networks and traditional dynamic systems. Based on the connection, we further propose continuous graph neural networks (CGNNs), which generalise existing discrete graph neural networks to continuous cases through defining the evolution of node representations with ODEs. Our ODEs are motivated by existing diffusion-based methods on graphs, where two different ways are considered, including different feature channels change independently or interact with each other. Extensive theoretical and empirical analysis prove the effectiveness of CGNN over many existing methods. Our current approach assumes that connected nodes are similar (‘homophily’ assumption), we leave it for future work to be able to learn more complex non-linear relationships such as can be found in molecules or knowledge graphs .

This project is supported by the Natural Sciences and Engineering Research Council (NSERC) Discovery Grant, the Canada CIFAR AI Chair Program, collaboration grants between Microsoft Research and Mila, Amazon Faculty Research Award, Tencent AI Lab Rhino-Bird Gift Fund and a NRC Collaborative R&D Project (AI4D-CORE-06).

We also would like to thank Joey Bose, Alexander Tong, Emma Rocheteau, and Andreea Deac for useful comments on the manuscript and Sékou-Oumar Kaba for pointing out a mistake in one equation.

References

Appendix A Proof of Proposition 1 and 2

Proof of Proposition 1: The starting point is to see

where $\Delta t=\frac{t+1-0}{n+1}$ with $t=n$ and $\boldsymbol{E}=\boldsymbol{H}_{0}$ as before. So now letting $n\to\infty$ we get the following integral

However, $\boldsymbol{A}^{t+1}$ is intractable in practice to compute for non-integer $t$, hence we solve the ODE by considering the second order ODE and then integrating again.

and solving for the constant using the fact that

we will make use of a Ansatz and use the integrating factor $\exp(-(\boldsymbol{A}-\boldsymbol{I})t)$.

Appendix B Proof of Proposition 3 and 4

To prove Proposition 3 and 4, we first prove the following Lemmata:

The analytical solution of the following ODE,

where $\boldsymbol{B}$ and $\boldsymbol{C}$ have eigenvalue decompositions $\boldsymbol{P}\boldsymbol{\Lambda}\boldsymbol{P}^{-1}$ and $\boldsymbol{Q}\boldsymbol{\Phi}\boldsymbol{Q}^{-1}$ respectively, with initial value $\boldsymbol{H}(0)$ is:

with $\widetilde{\boldsymbol{D}}=\boldsymbol{P}^{-1}\boldsymbol{D}\boldsymbol{Q}$.

Proof: First note that the ODE in Eq. (32) is known as the Sylvester ODE with analytical solution:

Hence, to prove Lemma 1 it remains to solve the integral using the help our assumptions.

Assuming that $\boldsymbol{A}$ and $\boldsymbol{W}$ have eigenvalue decompositions $\boldsymbol{P}\boldsymbol{\Lambda}\boldsymbol{P}^{-1}$ and $\boldsymbol{Q}\boldsymbol{\Phi}\boldsymbol{Q}^{-1}$ respectively,

with $\widetilde{\boldsymbol{E}}=\boldsymbol{P}^{-1}\boldsymbol{E}\boldsymbol{Q}$.

Proof: We start by using the eigenvalue decompositions of $\boldsymbol{A}$ and $\boldsymbol{W}$

where $\widetilde{\boldsymbol{E}}=\boldsymbol{P}^{-1}\boldsymbol{E}\boldsymbol{Q}$. Now we can consider the integral element-wise to get the required result:

Proof of Proposition 3: For the discrete dynamic defined

$\boldsymbol{H}_{n}$ can be rewritten as follows:

Recall that as in the case where $\boldsymbol{W}=\boldsymbol{I}$ we move to a continuous setting by interpreting the equation as a Riemann integral:

To derive a corresponding ODE, we consider the derivative of $\boldsymbol{H}(t)$ with respect to $t$, yielding the ODE below:

To get the ODE in a nicer form, we consider the second-derivative of $\boldsymbol{H}$:

By integrating over $t$ in both sides of the above equation, we can obtain:

Note that the initial value of the ODE is $\boldsymbol{H}(0)$ by Lemma 2, we know that

where $\widetilde{\boldsymbol{E}}=\boldsymbol{P}^{-1}\boldsymbol{E}\boldsymbol{Q}$.

Proof of Proposition 4: Proposition 4 follows trivially from Lemma 1. QED $\Box$

Appendix C Hyperparameters & training details

Across Cora, Citeseer, and Pubmed we use a hidden dimension of 16, input dropout of 0.5 in the encoder, and weight decay of $5\times 10^{-4}$. For NELL we use a hidden dimension of 64, dropout of 0.1 (in both encoder and decoder), and weight decay of $1\times 10^{-5}$.$\gamma$ determines the weight of self-loops in the graph.

In practice, we used the augmentation proposed in Dupont et al. to stabilise training, however, it had little no effect on the final performance. Writing down the augmentation mathematically yields the following matrix differential equation for the ODE

where $\boldsymbol{O}(0)=0$. All the augmentation does is add some latent dimensions to allow the trajectory of the ODE (which cannot cross itself) a potentially simpler path.

Appendix D Running time

In Table 8 we compare the running times between our algorithms and GCN as measured on the Cora dataset using the hyperparameters in Table 4. For GCN we used the hyperparameters quoted in . For all models we used 400 epochs and measured the total running time. All results were collected on a single machine with the following specs: Intel Core i7-6700HQ CPU at 2.60GHz with 16GB of RAM with NVIDIA GeForce GTX 960M with 2GB of dedicated GPU memory.

Appendix E Comparison to GCN with skip connections

We compare with GCN with residual links and with 2, 4, 8, and 16 layers respectively. The results on Cora, Citeseer, and PubMed in both the fixed data split setting and random split setting are given in Table 9.