GRAND++: GRAPH NEURAL DIFFUSION WITH A SOURCE TERM
Year: 2022

Abstract:
We propose GRAph Neural Diffusion with a source term (GRAND++) for graph deep learning with a limited number of labeled nodes, i.e., low-labeling rate. GRAND++ is a class of continuous-depth graph deep learning architectures whose theoretical underpinning is the diffusion process on graphs with a source term. The source term guarantees two interesting theoretical properties of GRAND++: (i) the representation of graph nodes, under the dynamics of GRAND++, will not converge to a constant vector over all nodes even as the time goes to infinity, which mitigates the over-smoothing issue of graph neural networks and enables graph learning in very deep architectures. (ii) GRAND++ can provide accurate classification even when the model is trained with a very limited number of labeled training data. We experimentally verify the above two advantages on various graph deep learning benchmark tasks, showing a significant improvement over many existing graph neural networks.

heading: 1 INTRODUCTION
text: Graph neural networks (GNNs) are the backbone for deep learning on graphs. Recent GNN architectures include graph convolutional networks (GCNs) [30], ChebyNet [16], GraphSAGE [29], neural graph fingerprints [20], message passing neural networks [28], and graph attention networks (GATs) [54]. These graph deep networks have achieved success in many applications, including computational physics and computational chemistry [20, 28, 3], recommender systems [41, 62], and social networks [63, 47]. Hyperbolic GNNs have also been proposed to enable certain kinds of data embedding with much smaller distortion [11, 37]. See [6] for some recent advances of GNN algorithm development and applications.
A well-known problem of GNNs is that increasing the depth of GNNs often results in a significant drop in performance on various graph learning tasks. This performance degradation has been widely interpreted as the over-smoothing issue of GNNs [35, 44, 12]. Intuitively, GNN layers update the node representation by taking a weighted average of its neighbors’ features, making representations for neighboring nodes to be similar. As the GNN architecture gets deeper, all nodes’ representation will become indistinguishable resulting in over-smoothing. In Sec. 2, we briefly show that certain GNNs have a diffusive nature which makes over-smoothing inevitable. Another interesting interpretation of the GNN performance degradation is via a bottleneck [1], since a GNN tends to represent exponentially growing information from neighbors with fixed-size vectors. Several algorithms have been proposed to mitigate the over-smoothing of GNNs, including skip connection and dilated convolution [33], Jumping Knowledge [60], DropEdge [49], PairNorm [64], graph neural diffusion (GRAND) [10], and wave equation motivated GNNs [21]. Nevertheless, developing deep GNN architectures is still in its infancy compared to the development of other deep networks.
Besides suffering from over-smoothing, we notice that the accuracy of existing GNNs drops severely when they are trained with a limited labeled data. As illustrated in Fig. 1, the test accuracy of several
∗Correspond to wangbaonj@gmail.com or matthew.thorpe-2@manchester.ac.uk
celebrated GNN architectures, including GCN, GAT, and GraphSage, drops rapidly when they are trained with fewer labeled data. Moreover, the variance of classification accuracy grows significantly as number of labeled nodes drops. Indeed, semi-supervised graph learning with very low-labeling rates has been studied in the Laplace learning and graph deep learning settings, see, e.g., [36, 9, 23]; one question is can we develop new GNN architectures to improve the performance of graph deep learning in low-labeling rate regimes?

heading: 1.1 OUR CONTRIBUTION
text: With the above GNN problems in mind, we focus on developing new continuous-depth GNNs to overcome over-smoothing and boost the accuracy of GNNs with a limited number of labeled data. We first present a random walk interpretation of the GRAND model [10], revealing a potentially inevitable over-smoothing phenomenon when GRAND is implicitly very deep. Based on the random walk viewpoint of GRAND, we then propose graph neural diffusion with a source term (GRAND++) that corrects the bias arising from the diffusion process, see Sec. 5 for details. GRAND++ theoretically guarantees that: (i) under GRAND++ dynamics, the graph node features do not converge to a constant vector over all nodes even as the time goes to infinity, and (ii) GRAND++ can provide accurate prediction even when it is trained with a limited number of labeled nodes. Moreover, these theoretical results resonate with the practical advantages of GRAND++. We summarize the major practical advantages of GRAND++ below. • GRAND++ can effectively overcome the over-smoothing issue; it is remarkably more accurate
than existing GNNs when the architecture is very deep. • GRAND++ is suitable for graph deep learning when only a few nodes are labeled as training data.
Moreover, in the low-labeling rates, GRAND++ can be more accurate when the network is deeper. • GRAND++ inherits the continuous-depth merit from GRAND, which defines the network depth
implicitly and enables memory-efficient training by using the adjoint method.

heading: 1.2 RELATED WORK
text: Diffusion on graphs and continuous-depth graph neural networks. Diffusion has been defined on graphs, see, e.g., [25, 24], and used in various applications, including data clustering and dimension reduction [15, 4], image processing [27, 22, 17, 38], and semi-supervised graph nodes classification [67, 65]. From the numerical viewpoint, fast algorithms have been proposed for using diffusion on graphs to solve penalized graph cut problems [26]. The connection between GNNs and diffusion on graphs has been studied substantially. For instance, GNN has been interpreted as a diffusion process on graphs, which performs low-pass filtering on the input features [43]. Moreover, insights from the diffusion process on graphs have been used to improve the performance of GNNs, see, e.g., [2, 36, 31, 57].
Leveraging neural ordinary differential equations (ODEs) [13], continuous-depth GNNs have been proposed, see, e.g., [45, 58, 68]. One recent work is GRAND [10], which parameterizes the diffusion equation on graphs with a neural network. See Sec. 3 for a brief review of GRAND. Neural ODEs. Neural ODEs [13] are a class of continuous-depth neural networks whose depth is defined implicitly. Training neural ODEs using the adjoint method [46] is more memory efficient than training other neural networks using backpropagation. We provide a brief review of neural ODEs and the adjoint method in Appendix C. GRANDs [10] are a class of neural partial differential equations (PDEs) on graphs that can also be considered as a coupled system of neural ODEs. Furthermore, GRANDs are also trained by using the adjoint method.
Laplace learning and Poisson learning. Laplace learning has been used for semi-supervised data classification [67, 65, 56], image processing [7, 27], etc. Direct application of Laplace learning with Gaussian weights [5] or locally linear embedding weights [50] for the above tasks may cause inference inconsistency when only a limited number of graph nodes are labeled, resulting in poor performance. Several algorithms address the inference inconsistency at low labeling rate. They include up-weighting the weights of the labeled data [52] and the p-Laplacian [8, 48, 66]. In [9], the authors have proposed Poisson learning for improving Laplace learning at extremely low-labeling rate regimes. Poisson learning augments Laplace learning with a Green’s function at each labeled data, enabling accurate node classification when only a few labeled data are available. Compared to Laplace learning, Poisson learning adds Green’s function to the label of each labeled node and then performs label propagation to predict the label for unlabeled graph nodes. GRAND and GRAND++ both learn graph node representations and perform prediction by activating the node representations, which are fundamentally different from Laplace and Poisson learning.

heading: 1.3 NOTATION
text: We denote scalars by lower- or upper-case letters and vectors and matrices by lower- and upper-case boldface letters, respectively. For a matrix A, we denote its transpose as A> and its Hadamard product with another matrix B as A B, i.e., the entrywise multiplication of A and B. We write the set {1, 2, · · · , n} as [n]. We denote the probability and expectation of a given random variable x as P(x) and E[x], respectively. The meaning of other notations can be inferred from the context.

heading: 1.4 ORGANIZATION
text: The paper is organized as follows: In Sec. 2, we review diffusion equation on graphs and its connection to GNNs. In Secs. 3 and 4, we briefly review GRAND and present a random walk interpretation of GRAND, respectively. Leveraging the random walk viewpoint of GRAND, we propose GRAND++ for deep graph learning with theoretical guarantees in Sec. 5. We verify the efficacy of GRAND++ in Sec. 6. Technical proofs and more validations are provided in the appendix.

heading: 2 BACKGROUND
text: Diffusion equation on graphs. Let G = (X,W ) represent an undirected graph with n nodes, where X = ( [x(1)]>, · · · , [x(n)]>
)> ∈ Rn×d with each row x(i) ∈ Rd a feature vector and W := ( Wij ) a n × n matrix with Wij representing the similarity (edge weight) between the ith and jth feature vectors, and we assume Wij = Wji. Consider the following diffusion process that evolves the feature matrix X on the graph (see Appendix A for a brief review of calculus on graphs):
∂X(t)
∂t = div
( G(X(t), t) ∇X(t) ) , (1)
where X(t) = ( [x(1)(t)]>, · · · , [x(n)(t)]> )> ∈ Rn×d with x(i)(0) = x(i), ∇ and div are the gradient and divergence operators, respectively. The matrix G(X(t), t) is chosen such that W G is right-stochastic, i.e., each row of W G summing to 1. In the machine learning setting, we can parameterize G with learnable parameters θ which we denote by G(X(t), t, θ). The initial features are evolved under the diffusion dynamics (1) from t = 0 to T to learn the final representation X(T ) for further machine learning tasks.
In the simplest case when G(X(t), t) is only dependent on the initial node features X , i.e., G is time-independent, right-stochasticity implies ∑ jWijGij = 1 for all i, and so we focus on the
particular case when Gij = 1/di with di = ∑n j=1Wij . In this case the right-hand side of (1) reduces to the negative of the random-walk Laplacian applied to X(t) and (1) becomes
∂X(t)
∂t = div
( G(X(t), t) ∇X(t) ) = −LX(t), (2)
where L = I −D−1W := I −A (A := A(X)) is the random walk Laplacian and D is diagonal with Dii = di. See [14, 25, 24] for more about random walk Laplacian and diffusion on graphs.
Graph neural networks. Applying forward Euler discretization, with step size δt < 1, of (2) gives
X(kδt) = X((k − 1)δt)− δtLX((k − 1)δt) := L̃X((k − 1)δt), for k = 1, 2, · · · ,K, (3)
T = Kδt, and X(0) = X . Note that the matrix L̃ is the discretization of the diffusion operator, which is a special low-pass filter. Equation (3) is a prototype for motivating GNNs: by introducing weights W (k) ∈ Rd×d and a nonlinearity σ, e.g., ReLU, into (3), we have
X((k + 1)δt) = σ ( L̃X(kδt)W (k) ) . (4)
The model in (4) is similar to the well-established GCN architecture proposed in [30]. The diffusive nature of the GNN architecture in (4) further explains the over-smoothing issue of training deep GNNs; the deeper the network architecture is, the more the node features diffuse. Eventually, all nodes share similar features and become indistinguishable. See Sec. 5 for a detailed analysis.

heading: 3 A BRIEF REVIEW OF GRAND
text: GRAND is a new continuous-depth GNN proposed in [10]. It integrates a learnable encoder function φ and a learnable decoder function ψ with the neural network parameterized graph diffusion process, resulting in the prediction Y = ψ(X(T )), where X(T ) is computed as
X(T ) = X(0) + ∫ T 0 ∂X(t) ∂t dt, with X(0) = φ(X), (5)
where ∂X(t)/∂t is given by the graph diffusion equation (2). From the neural ODE perspective, we can perform forward propagation of GRAND, i.e., we solve (5), using numerical ODE solvers.
In the simplest case, when G is only dependent on the initial node features, we can rewrite (1) as
∂X(t) ∂t = ( A(X)− I ) X(t), (6)
GRAND models the diffusivity A(X) in (6) by the multi-head self-attention mechanism; potential choices of the attention function include the ones proposed in [53, 54]. More precisely, in GRAND A(X) = 1h ∑h l=1 A
l(X) with h being the number of heads and the attention matrix Al(X) = (al(xi,xj)), for l = 1, · · · , h, is computed as follows:
al(xi,xj) = exp
( LeakyReLU(al>[W lxi‖W lxj ]) )∑ k∈Ni exp ( LeakyReLU(al>[W lxi‖W lxk])
) , (7) where W l and al are learned, ‖ is the concatenation operator, and Ni is the index set of the nodes that are connected to the ith node in the graph. GRAND with the attention in (7) is called GRAND-l, that is, GRAND-l is a special case of GRAND when the diffusivity is dependent only on the initial graph node features. Time-dependent attention and graph rewiring can be integrated into GRAND, resulting in GRAND-nl and GRAND-nl-rw, respectively [10]. From the ODE viewpoint, GRAND and its variants are a class of coupled neural ODEs defined on an unweighted graph. Their merits include continuous-depth and memory-efficient training using the adjoint method [46, 13].

heading: 4 RANDOM WALK VIEWPOINT OF GRAND
text: In this section, we present a random walk interpretation of GRAND. The connection between graph random walks and the diffusion equation has been extensively studied, but we recap the key idea here to motivate the new GRAND with a source term architecture. Let {B(i)(k)}k∈N be the random walk on {x(j)(0)}nj=1 defined by, for δt ∈ [0, 1],
B(i)(0) = x(i)(0)
P ( B(i)(k + 1) = x(`)(0)|B(i)(k) = x(j)(0) ) = { 1− δt if ` = j δtWj` dj if ` 6= j (8)
where dj = ∑n `=1Wj` (assume W`` = 0 for all `). Proposition 1 below is well-known, see [67]. We provide the proof of Proposition 1 and all the subsequent theoretical results in Appendix B.
Proposition 1 Let X solve (3) and B(i) be the random walk determined by (8) where δt ∈ [0, 1]. Then x(i)(δtk) = E [ B(i)(k) ] .
Proposition 2 below gives the stationary distribution of the random walk {B(i)(k)}k∈N.
Proposition 2 Assume the graph G = (X,W ) is connected. Then, the stationary distribution of {B(i)(k)}k∈N is
π = ( d1∑n j=1 dj , . . . , dn∑n j=1 dj ) , (9)
which is independent of the starting position x(i).
Furthermore, we have the following theoretical result on the asymptotic behavior of graph node features under the GRAND dynamics given by (3).
Proposition 3 Assume the graph G = (X,W ) is connected. Then for all i = 1, · · · , n, we have
x(i)(kδt)→ x̃ := n∑ j=1 x(j)(0)πj , as k →∞.
Hence, for the case of (3), i.e., GRAND-l, we expect the output to be approximately independent of the input, due to over-smoothing. Of course, once we reintroduce the X(t) dependence back into G in (1) and (2) or into the operator A in (6) then the above arguments no longer hold. Nevertheless, the GRAND architectures are built on a principle that is ill-suited to deep networks. In the next section we introduce a source term and perform a similar random walk analysis that illustrates how the new architecture can be better suited for deep GNN architectures.

heading: 5 GRAND++: GRAPH NEURAL DIFFUSION WITH A SOURCE TERM

heading: 5.1 ALGORITHM AND FORMULATION
text: At the core of GRAND++ is the introduction of a source term into GRAND, leveraging the random walk viewpoint of the diffusion process. We take a small subset of feature vectors, indexed by I ⊆ [n], believed to be “trustworthy” for use as a source term. In particular, we use the features of labeled data. The GRAND++ dynamics are defined by a diffusion equation with a source term (we use the variable z for GRAND++-related dynamics and x for GRAND dynamics)
∂z(i)(t)
∂t = div [G(Z(t), t) ∇Z(t)](i) + ∑ j∈I δijCj (10)
where Cj is the source at feature vector of node j. Below we motivate a particular choice of Cj .
The key idea is to first characterise the bias that arises from the diffusion and use that to propose a correction via the choice of source terms Cj . Following the simplifications in (2), our diffusion equation (without the source term) follows the approximate dynamics when t 1
∂x(i)(t)
∂t = − [LX(t)](i) = −x(i)(t)︸ ︷︷ ︸
≈x̃
+ 1
di n∑ j=1 Wij x (j)(t)︸ ︷︷ ︸ ≈x̃ ≈ 0.
For i ∈ I, it transpires that choosing Ci = x(i) − x̂ (where x̂ is defined below) gives rise to a random walk interpretation that allows us to prove that the oversmoothing seen in the GRAND model is avoided.
One can in fact choose x̃ with a certain degree of freedom. If we initialise X(0) = X then we obtain x̃ = ∑n j=1 x
(j)πj (as is usual in the GRAND model). However, as the similarities are encoded in the graph weights, and the diffusion dynamics will drive it towards a non-trivial state, we can choose a different initialization than X(0) = X . Through connections with random walks we, in the next subsection, motivate an alternative initialisation
n∑ i=1 z(i)(0) = ∑ i∈I x(i) − x̂ di , where x̂ = 1 |I| ∑ j∈I x(j) (11)
with the dynamics ∂z(i)(t)
∂t = div [G(Z(t), t) ∇Z(t)](i) + ∑ j∈I δij ( x(i) − x̂ ) . (12)
For example, we could choose
z(i)(0) = { 1 di ( x(i) − x̂ ) if i ∈ I
0 otherwise, or z(i)(0) = x(i) − c,
where c = 1n (∑n i=1 x (i)− ∑ j∈I x(j)−x̂ dj ) is chosen such that (11) holds. We do not believe that the constant c (that shifts by a constant) is particularly important but it is included to provide a random walk interpretation which helps to understand the deep architecture (when T is big) behaviour of the model GRAND++. The justification for this choice will be made in Sec. 5.2. To summarize, the GRAND++ model in (12), with initial condition satisfying (11), simply adds a source term to the original GRAND model and uses a different initial condition. Therefore, the nonlinear diffusivity and graph rewiring tricks used by GRAND can be easily integrated into GRAND++. In terms of implementation, since GRAND++ merely changes the right-hand side of GRAND, which again can be regarded as a system of coupled first-order neural ODEs; we can leverage neural ODE training, testing, and inference for GRAND++ similar to GRAND.
In the next subsection we explore the random walk connection of the above model, suggesting that building a graph neural network based on the diffusion with source model does not suffer from the same degeneracy as we observed in Sec. 4 and is therefore better suited to build deep GNNs. In particular, we can write the diffusion with source model as the short time expected behaviour of a random walk and therefore we do not have the issue of reaching the stationary state (in other words passing the mixing time). Our experiments in Sec. 6 suggest the formal motivation holds and we are able to design deep GNNs.

heading: 5.2 THE RANDOM WALK PERSPECTIVE OF GRAND++
text: Let us continue to consider the simplified model in the previous subsection, i.e., assume the dynamics are governed by
∂z(i)(t)
∂t = − [LZ(t)](i) + ∑ j∈I δij ( x(i) − x̂ ) (13)
where the initial condition satisfies (11). Using the forward Euler discretisation of the above dynamics we have
z(i)(δtk) = z (i)(δt(k − 1))− δt [LZ(δt(k − 1))](i) + δt ∑ j∈I δij ( x(i) − x̂ ) , (14)
for k = 1, 2, . . . ,K where again T = Kδt.
We use the same random walk as that introduced in Sec. 4, i.e. the random walk defined by (8), but we will now only consider random walks that are initialised on the nodes indexed by I.
Proposition 4 Let Z solve (14) with the initial condition satisfying (11), and let B(i) be the random walk determined by (8). Then,∣∣∣∣z(i)(kδt)− E[ k∑
s=0
1
di ∑ j∈I ( x(j) − x̂ ) 1B(j)(s)=x(i) ]∣∣∣∣→ 0 as k →∞. Remark 1 In the limit k →∞ the term
E [ k∑ s=0 1 di ∑ j∈I x(j)1B(j)(s)=x(i) ] is formally a function of the random walk at all times. Whilst if k is very large (i.e. in comparison to the mixing time) we still have that
E [ 1
di ∑ j∈I x(j)1B(j)(k)=x(i) ] = 1 di ∑ j∈I x(j) P ( B(j)(k) = x(j) ) ︸ ︷︷ ︸
≈πi
≈ πi di ∑ j∈I x(j) (15)
and on the other hand E [ 1
di ∑ j∈I x̂1B(j)(k)=x(i) ] = x̂ di ∑ j∈I P ( B(j)(k) = x(j) ) ︸ ︷︷ ︸
≈πi
≈ πi|I|x̂ di . (16)
From the definition of x̂ we see that (15) and (16) are approximately equal. And therefore we can understand E[ 1di ∑ j∈I x̂1B(j)(k)=x(i) ] as the long time behaviour of E[ 1 di ∑ j∈I x
(j)1B(j)(k)=x(i) ] Very formally we can see that subtracting the long-time behaviour from the all-time behaviour leaves us with the short time behaviour. This provides one explanation as to why we do not expect the deep layers to be determined by the stationary state of the random walk (at which point there is little dependence on the initial layers, causing the deep layers to be approximately constant).
Remark 2 The random walk interpretation
E [ k∑ s=0 1 di ∑ j∈I ( x(j) − x̂ ) 1B(j)=x(i) ] (17)
can be considered to be dual to the random walk interpretation in Sec. 4: in Sec. 4 we released the random walker from the node of interest, whilst now we release the random walkers from nodes indexed by I and see how many of them hit the node of interest. We note also that we do not require a lower bound on the size of the set I. Indeed, if |I| is fixed whilst one takes the number of feature vectors n→∞ we still expect many properties of GRAND++, in particular Proposition 5 below, to hold. This is due to the asymptotic well-posedness of the dual random walk in low labeling rates [9].
Proposition 3 reveals that in the simple setting of (2), GRAND converges to a constant when its depth goes to infinity. However, this is not true for GRAND++ since the graph node features will not converge to a constant vector driven by the GRAND++, as shown in Proposition 5 below.
Proposition 5 Assume the graph G = (X,W ) is connected. Then z(i)(kδt) that was defined in (14) does not converge to a constant vector as a function of i as k →∞. That is, the node features will not become the same across graph nodes under the GRAND++ dynamics.
Remark 3 Proposition 5 guarantees GRAND++ is less likely to suffer from over-smoothing than GRAND, and in particular it shows that we have a non-constant deep layer limit, i.e., as t → ∞. Analysing the limit is beyond the scope of the paper but we have seen one characterisation in Proposition 4. By construction we have ∂z(i)(t)/∂t ≈ 0 for i ∈ I so one should expect that the deep layer limit is (close to) a smooth interpolation of the feature vectors labeled by I.
The continuous time model (10) is, in the special case of (13), the mean field limit of the probabilistic formulation (17). Our proposed algorithm is formulated from the mean-field limit.

heading: 6 EXPERIMENTS
text: In this section, we compare the performance of GRAND++ with GRAND and several other popular GNNs on various graph node classification tasks. We aim to show the practical advantages of GRAND++ in learning with limited labeled data and using deep architectures. Without mentioning clearly, we use the same hyperparameters that that used for GRAND in [10] for GRAND++. We provide detailed descriptions of experimental settings and datasets that are omitted in the main text in Appendix D.1. For all experiments, we run 100 splits for each dataset with 20 random seeds for each split, which are conducted on a server with four NVIDIA RTX 3090 graphics cards.
We compare the performance of GRAND++ and its nonlinear and graph rewiring variants with several popular GNNs on various graph node classification benchmarks. Except for the integration time, which measures the implicit depth of GRAND and GRAND++, we adopt the experimental settings of GRAND in [10] for GRAND++ include numerical differential equation solvers. Following [10], we study seven graph node classification datasets, namely CORA, CiteSeer, PubMed, CoauthorCS, Computer, Photo, and ogbn-arxiv; we describe these datasets in Appendix D.1.

heading: 6.1 GRAND++ IS MORE RESILIENT TO DEEP ARCHITECTURES
text: We first show that our introduced source term in (12) can improve the accuracy of GRAND-l when the architecture is deep, i.e., the integration time T in (5) is big. We denote GRAND-l with the source term as GRAND++-l. For each node classification task, we train all models using the same number of labeled nodes as in [10]. Figure 2 contrasts the performance of GRAND-l and GRAND++-l with different depths, or T , on CORA, CiteSeer, Computer, and Photo datasets. We provide the detailed results on PubMed and CoauthorCS, together with more comparisons of GRAND++-l with GRAND-l and several other celebrated GNNs include GCN, GAT, and GraphSage in Table 5 in Appendix D.2. The results in Fig. 2 and Appendix D.2 confirm that GRAND-l suffers less from over-smoothing compared to GCN, GAT, and GraphSage. Moreover, GRAND++-l performs on par with GRAND-l when the depth (T ) of the network is small, but GRAND++-l significantly outperforms GRAND-l when T is large. As T increases, the margin becomes wider, indicating that GRAND++-l can overcome over-smoothing much more effectively than GRAND-l. Note that we did not use uniform depth for GRAND-l and GRAND++-l on all datasets because the adaptive step-size ODE solver fails when T is large for some tasks.
Next, we compare GRAND-l and GRAND++-l on the ogbn-arxiv node classification task, which is a large-scale benchmark. We train two models using labeling rates of 3.0% and 5.0%, respectively; the corresponding test accuracy for GRAND-l/GRAND++-l are 65.26%/66.64% and 67.42%/67.77%, respectively. GRAND++-l outperforms GRAND-l in both labeling rates.We further compare GRAND and GRAND++ with different depth on the ogbn-arxiv task in Appendix D.6.

heading: 6.2 GRAND++ IS MORE ACCURATE WITH LIMITED LABELED TRAINING DATA
text: Besides helping to overcome over-smoothing, our theory shows that the source term can boost the accuracy of GRAND-l with lowlabeling rates. Table 1 compares the accuracy of GRAND++-l with GRAND-l, GCN, GAT, GraphSage, and MoNet, trained with different numbers of labeled data. Here, we slightly tune T for GRAND++ based on the optimal value for GRAND, see Table 4 in the Appendix for their values. We see that with few labeled data, in most tasks GRAND++-l is significantly more accurate than the other GNNs include GRANDl, confirming our theoretical insight. For Coau-
thorCS task, both GRAND-l and GRAND++-l are worse than GCN and GraphSage. Moreover, increasing the depth of GRAND++-l can improve the classification accuracy with limited training data, but this is not the case for GRAND-l, see Fig. 3. We perform the t-test in Appendix D.5 to confirm the statistical significance of the accuracy gain of GRAND++ over GRAND in Table 1.

heading: 6.3 TIME-DEPENDENT ATTENTION AND GRAPH REWIRING
text: The previous experimental results show that GRAND++-l enhances the accuracy of GRAND-l in the cases when the labeled training data is limited and when the network is deep. Here, we explore the same strategy for GRAND-nl and GRAND-nl-rw; we name the corresponding models with the new source term GRAND++-nl and GRAND++-nl-rw, respectively. Table 2 compares GRAND-nl and
GRAND-nl-rw with the corresponding model with a source term. We see that overall GRAND++-nl (GRAND++-nl-rw) outperforms GRAND-nl (GRAND-nl-rw) when the network is deep, i.e., T is big. We further study the low-labeling rate regimes in Appendix D.3.

heading: 7 CONCLUDING REMARKS
text: We propose GRAND++, which augments graph neural diffusion with a source term. We present some theory that connects the model to a random walk formulation on graphs. GRAND++ outperforms many existing GNNs for graph deep learning with very deep architectures and when the number of labeled data is limited. GRAND++ can be regarded as coupled ODE system in which each ODE has an external force term. As such, it is natural to consider if advanced techniques in accelerating training, test, and inference of neural ODEs can be leveraged to improve the efficiency and accuracy of GRAND++, in particular high-order neural ODEs [19, 61, 42, 59] and noise injection [55]. It is interesting to note that the second-order neural ODE can be connected to the wave equation in the graph setting, which can automatically bypass over-smoothing. We leave studying the second-order neural ODE on graphs as future work.

heading: 8 ACKNOWLEDGEMENT
text: This material is based on research sponsored by NSF grants DMS-1924935, DMS-1952339, DMS2027248 and NSF CCF-1934568, DOE grant DE-SC0021142, and ONR grant N00014-18-1-2527 and the MURI grant N00014-20-1-2787. MT would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme Mathematics of Deep Learning when work on this paper was undertaken (EPSRC grant number EP/R014604/1) and acknowledge support from the European Union Horizon 2020 research and innovation programmes under the Marie Skłodowska-Curie grant agreement No. 777826 (NoMADS). MT also holds a Turing Fellowship at the Alan Turing Institute.

heading: A BACKGROUND ON GRAPH DIFFERENTIAL OPERATORS
text: Let (X,W ) represent a graph where X = ([x(1)]>, . . . , [x(n)]>)> ∈ Rn×d is the matrix where each row x(i) ∈ Rd is a feature vector and W = (Wij)ni,j=1 is a n×nmatrix withWij representing the similarity (edge weight) between the ith and jth feature vector. We assume that we are dealing with an undirected graph, i.e., Wij = Wji. A Rk1 -valued function on the nodes of the graph can be represented as a matrix U ∈ Rn×k1 by U = ([u(1)]>, . . . , [u(n)]>)> and we define the inner product
〈U ,V 〉 = n∑ i=1 u(i) · v(i).
Similarly, a Rk2 -valued function on the edges can be represented as a third-order tensor U ∈ Rn×n×k2 which we write as
U =  U (1,1) · · · U (1,n) ... . . . ...
U (n,1) · · · U (n,n)  and U (i,j) ∈ Rk2 . On edge functions we use the inner product
〈U ,V〉 = 1 2 n∑ i,j=1 WijU (i,j) · V(i,j).
Multiplication between a matrix and an edge function is usually defined pointwise and we use the notation [A U ](i,j) = AijU (i,j) ∈ Rk2 , for a matrix A ∈ Rn×n and an edge function U ∈ Rn×n×k2 , to make this clear. Similarly, pointwise multiplication between two matrices A,B ∈ Rn×n is defined by [A B]ij = AijBij ∈ R. When a matrix is acting as a linear operator on a node function we use the usual matrix-vector notation and write [AU ](i) = ∑n j=1Aiju
(j) ∈ Rk1 for a matrix A ∈ Rn×n and node function U = ([u(1)]>, . . . , [u(n)]>)> ∈ Rn×k1 . In the sequel we will have k1 = k2 = d.
The gradient of a node-function U = ([u(1)]>, . . . , [u(n)]>)> ∈ Rn×d is defined as the edgefunction ∇U ∈ Rn×n×d with [∇U ](i,j) = u(j) − u(i) ∈ Rd. The divergence divV = ([[divV](1)]>, . . . , [[divV](n)]>)> ∈ Rn×d of an edge-function V ∈ Rn×n×d is defined as
[divV](i) = n∑ j=1 WijV(i,j)
for all i = 1, . . . , n. For anti-symmetric edge functions, i.e. V(i,j) = −V(j,i) for all i, j, we have that the divergence is the negative adjoint to the gradient, i.e.
〈divV,U〉 = −〈V,∇U〉.

heading: B TECHNICAL PROOFS
text: Proof: [Proof of Proposition 1] For notational convenience let us assume that x(i)(0) = x(i). Clearly
E [ B(i)(0) ] = x(i) = x(i)(0)
for all i = 1, . . . , n. Assume that
E [ B(i)(k) ] = x(i)(δtk)
for all i = 1, . . . , n. Then, E [ B(i)(k + 1) ] =
n∑ j=1 x(j)P ( B(i)(k + 1) = x(j) ) =
n∑ j=1 n∑ `=1 x(j)P ( B(`)(k) = x(j)|B(i)(1) = x(`) ) P ( B(i)(1) = x(`) ) =
n∑ j=1 n∑ `=1 x(j) ( (1− δt)1i=` + δtWi` di ) P ( B(i)(1) = x(`) ) = (1− δt)
n∑ j=1 x(j)P ( B(i)(1) = x(`) ) + δt di n∑ `=1 Wi` n∑ j=1 x(j)P ( B(`)(k) = x(j) ) = (1− δt)E [ B(i)(k) ] + δt di n∑ `=1 Wi`E [ B(`)(k)
] = (1− δt)x(i)(δtk) +
δt di n∑ `=1 Wi`x (`)(δtk)
= x(i)(δtk) + δt di n∑ `=1 Wi`Wi` ( x(`)(δtk)− x(i)(δtk) ) = x(i)(δtk)− δt [LX(δtk)](i)
= x(i)(δt(k + 1)),
as required. Proof: [Proof of Proposition 2] Let P = ( Pij ) ∈ Rn×n be the probability transition kernel, so
Pij = { 1− δt if i = j, δtWij di
if i 6= j. We have
n∑ i=1 πiPij = n∑ i=1 di∑n k=1 dk ( (1− δt1i=j + δtWij di ) = dj(1− δt)∑n k=1 dk + δt ∑n i=1Wij∑n k=1 dk
= dj∑n k=1 dk
= πj ,
as required.
Proof: [Proof of Proposition 3] The proof follows from a simple application of Propositions 1 and 2. Namely, for any i ∈ {1, . . . , n}
x(i)(kδt) = E [ B(i)(k) ] = n∑ j=1 x(j)(0)P ( B(i)(k) = x(j) ) → n∑ j=1 x(j)(0)πj = x̃
as k →∞.
Proof: [Proof of Proposition 4] Let
y(i)(k) = E  k∑ s=0 1 di ∑ j∈I ( x(j) − x̂ ) 1B(j)(s)=x(i)  .
Notice that
E [ k∑ s=0 1B(j)(s)=x(i) ]
= k∑ s=0 P ( B(j)(s) = x(i) ) = P ( B(j)(0) = x(i)
) ︸ ︷︷ ︸
δij
+ k∑ s=1 P ( B(j)(s) = x(i) )
= δij + k∑ s=1 n∑ `=1 P ( B(j)(s) = x(i)|B(j)(s− 1) = x(`) ) P ( B(j)(s− 1) = x(`) ) = δij +
k∑ s=1 n∑ `=1 ( (1− δt)δ`i + δtW`i d` ) P ( B(j)(s− 1) = x(`) ) = δij + (1− δt)
k∑ s=1 P ( B(j)(s− 1) = x(i) ) + δt n∑ `=1 W`i d` k∑ s=1 P ( B(j)(s− 1) = x(`) ) = δij + (1− δt)
k−1∑ s=0 P ( B(j)(s) = x(i) ) + δt n∑ `=1 W`i d` k−1∑ s=0 P ( B(j)(s) = x(`) ) = δij + (1− δt)E
[ k−1∑ s=0 1B(j)(s)=x(i) ] + δt n∑ `=1 W`i d` E [ k−1∑ s=0 1B(j)(s)=x(`) ] .
From the definition of Y and the above recursive relationship we have
y(i)(k) = 1
di ∑ j∈I ( x(j) − x̂ ) E [ k∑ s=0 1B(j)(s)=x(i) ]
= 1
di ∑ j∈I ( x(j) − x̂ ) δij + (1− δt) 1 di ∑ j∈I ( x(j) − x̂ ) E [ k−1∑ s=0 1B(j)(s)=x(`) ]
+ δt di n∑ `=1 W`i d` ∑ j∈I ( x(j) − x̂ ) E [ k−1∑ s=0 1B(j)(s)=x(`) ]
= 1
di ∑ j∈I ( x(j) − x̂ ) δij + (1− δt)y(i)(k − 1) + δt di n∑ `=1 W`iy (`)(k − 1)
= y(i)(k − 1) + 1 di ∑ j∈I ( x(j) − x̂ ) δij − δt [LY (k − 1)](i) .
Now we let
w(i)(k) = di ( z(i)(kδt)− y(i)(k) ) so that W satisfies
w(i)(k) = w(i)(k − 1)− δt [LW (k − 1)](i) = [PW (k − 1)](i) .
Hence, W (k) = P kW (0). Since the stationary distribution of the random walk with transition kernel is π, we have limk→∞P k = 1π>. Hence, as k → ∞, w(i)(k) → ∑n j=1 πjw (j)(0) = 0 since w(j)(0) = 0 by the choice in the initialisation of Z.
Proof: [Proof of Proposition 5] We prove the result by contradiction. Assume that there exists z such that z(i)(kδt)→ z for all i = 1, . . . , n as k →∞. Then LZ(kδt)→ 0. Since we can write
z(i)(kδt)− z(j)(kδt) = z(i)((k − 1)δt)− z(j)((k − 1)δt)
− δt ( [LZ((k − 1)δt)](i) − [LZ((k − 1)δt)](j) ) + δt
∑ `∈I ( δi`(x (i) − x̂)− δj`(x(j) − x̂) ) ,
then taking the limit k →∞ implies∑ `∈I δi`(x (i) − x̂) = ∑ `∈I δj`(x (j) − x̂)
for all i, j which is clearly not true.

heading: C NEURAL ODES AND TRAINING NEURAL ODES WITH ADJOINT METHOD
text: Neural ODEs [13] are a class of continuous-depth (-time) neural networks that are particularly suitable for learning complex dynamics from irregularly sampled sequential data, see, e.g., [13, 51, 19, 39, 42]. Mathematically, a neural ODE is the first-order ODE:
dh(t)
dt = f(h(t), t, θ), (18)
where f(h(t), t, θ) ∈ Rd is specified by a neural network parameterised by θ, e.g., a two-layer feedforward neural network. Starting from the input h(0), neural ODEs learn the representation and perform prediction by solving (18) from t = 0 to T using a numerical integrator with a given error tolerance, often with an adaptive step size solver (or adaptive solver for short) [18]. Solving (18) from t = 0 to T in a single pass with an adaptive solver requires evaluating f(h(t), t, θ) at various timesteps, with computational complexity counted by the number of forward function evaluations (forward NFEs) [13].
The adjoint sensitivity method (or adjoint method) [46], is a memory-efficient method for training neural ODEs. We regard the output h(T ) as the prediction and denote the loss between h(T ) and the ground truth as L. Let a(t) := ∂L/∂h(t) be the adjoint state, then we have (see [13, 46] for details)
dL dθ = ∫ T 0 a(t)> ∂f(h(t), t, θ) ∂θ dt, (19)
with a(t) satisfying the following adjoint ODE
da(t) dt = −a(t)> ∂ ∂h f(h(t), t, θ), (20)
which is solved numerically from t = T to 0 and also requires the evaluation of the right-hand side of (20) at various timestamps, and the backward NFEs measure the computational complexity.

heading: D EXPERIMENTAL DETAILS AND MORE EXPERIMENTAL RESULTS
text: D.1 DATASETS AND EXPERIMENTAL SETTINGS
Graph node classification dataset. Following [10], we consider the largest connected component of seven graph node classification datasets, including CORA, CiteSeer, PubMed, coauthor graph CoauthorCS, and Amazon co-purchasing graphs Computer and Photo, and a large scale ogbn-arxiv dataset. For completeness, we list the number of classes, the number of features, and the number of nodes and edges of each dataset in Table 3. More detailed information can be found in [10].
Depth of GRAND and GRAND++ for the results in Table 1. Table 4 lists the fine-tuned T for the results in Table 1. Due to the limited time, we only search around the value of optimal T for GRAND with grid spacing 0.1.
D.2 CLASSIFICATION ACCURACY OF GNNS WITH DIFFERENT DEPTHS
In this subsection, we provide detailed numbers that correspond to Fig. 1. We further compare GRAND++-l with several other GNN architectures, including GCN, GAT, and GraphSage, with different depths on a few benchmark datasets. Table 5 lists the classification accuracy of GRAND++, GRAND, and three benchmark GNNs with different depths on six graph node classification tasks. Again, we see that GRAND++ is better than the other GNN models when the networks are deep.
D.3 MORE RESULTS ON TIME-DEPENDENT ATTENTION AND GRAPH REWIRING
We further explore the effects of the source term for GRAND-nl and GRAND-nl-rw in the lowlabeling rate regimes. Table 6 compares GRAND-nl and GRAND-nl-rw with the corresponding model with a source term. We see that GRAND-nl and GRAND-nl-rw are almost always worse than the vanilla GRAND-l, consistent with the results reported in [10]. GRAND++-nl and GRAND++nl-rw cannot help learning at low labeling rates anymore. However, when the labeling rates are not low, GRAND++-nl or GRAND++-nl-rw can outperform GRAND-nl and GRAND-nl-rw, even outperform GRAND++.

heading: 1 50.55± 15.68 50.63± 17.71 48.89± 11.51 47.94± 11.06

heading: 2 65.06± 9.35 61.24± 16.19 59.96± 7.90 58.25± 11.97
text: D.4 CLASSIFICATION ACCURACY OF GNNS WITH FEWER NUMBER OF TRAINING DATA
Besides the results shown in Fig. 1, we further test the classification accuracy of more benchmark GNN architectures trained with fewer numbers of labeled data per class. Tables 7-9 list the classification accuracy, on the test set, of different benchmark GNN models when they are trained with different numbers of labeled nodes per class.
D.5 T-TEST OF THE ACCURACY IMPROVEMENT OF GRAND++ OVER GRAND
To confirm the statistical significance of the accuracy improvement of GRAND++ over GRAND in Table 1, in this subsection, we conduct t-test experiments at 0.95 confidence to compare GRAND and GRAND++ on six different benchmark graph node classification tasks. We first perform unpaired t-tests to show the improvement of GRAND++ over GRAND on low labeled datasets using the following t-score
t-score = µGRAND++ − µGRAND√
σ2GRAND++ n + σ2GRAND++ n
, (21)
where µ and σ2 are the mean and variance of the performances of each model, and n is the number of runs for each model. The t-test score are shown in Table 10.
For some entries in Table 10 that are not significant enough, we further conduct paired t-test between GRAND++ and GRAND on these specific datasets as shown in Table 11. Since a large portion of variance comes from splitting of the datasets, we pair up tests of GRAND and GRAND++ with the same splitting in this experiment. In this case, a sample of difference of size n is computed, and t-test score can be computed using the equation
D.6 TASKS FOR FURTHER EVALUATING DEEP GRAND AND GRAND++
Open graph benchmark with paper citation network (ogbn-arxiv). Ogbn-arxiv consists of 169, 343 nodes and 1, 166, 243 directed edges. Each node is an arxiv paper represented by a 128- dimensional features and each directed edge indicates the citation direction. This dataset is used for node property prediction and has been a popular benchmark to test the advantage of deep graph neural networks over shallow graph neural networks [34, 32]. Compared to the GRAND model used in [10], we reduce the hidden dimension from 162 to 81 to fit the model into the GPU in our lab.

Review summaries:
This paper studies neural-based diffusion models for graph data. It considers deep learning on graphs as a continuous diffusion process and treats GNNs as discretizations of an underlying PDE. The authors build upon an existing work (GRAND) by adding a source term to the objective function. The source term is designed to alleviate the over-smoothing problem of GNN, a well-known phenomenon that deeper GNN learns similar representation for each node. They theoretically show the effectiveness of the source term in mitigation over-smoothing problem. Furthermore, they empirically assess their method's effectiveness in learning with limited labeled data and using deep architectures on several benchmark datasets.

The paper introduces GRAND++, GRAph Neural Diffusion with a source term. As its predecessor GRAND, this method is focused on the development of a new continuous-depth GNN to tackle the over-smoothing and the bottleneck issues which are typical of deep GNNs. The authors introduce a random-walk interpretation of GRAND, which shows that when the network is very deep it is still prone to over-smoothing. To address this problem they add a source term to the problem formulation, which is shown to provide better results in very deep models and in classification tasks which are characterised by a low amount of labeled data.

In this paper, the authors propose a graph deep learning method called GRAND++. This framework based on GRAND and can work with a limited number of labeled nodes. Experiments are also conducted to demonstrate the effectiveness of the method.