Title: ORDINARY DIFFERENTIAL EQUATIONS ON GRAPH NETWORKS

Abstract: Recently various neural networks have been proposed for irregularly structured data such as graphs and manifolds. To our knowledge, all existing graph networks have discrete depth. Inspired by neural ordinary differential equation (NODE) for data in the Euclidean domain, we extend the idea of continuous-depth models to graphs, and propose graph ordinary differential equation (GODE). The derivative of hidden node states are parameterized with a graph neural network, and the output states are the solution to this ordinary differential equation. Noticing that NODEs are typically trained with the adjoint method, with the advantages of adaptive evaluation and a free-form continuous invertible model; however, their performance on benchmark image classification tasks is significantly inferior to discrete-layer models. We show the reason is adjoint method generates inaccurate gradient due to numerical error in reverse-mode integration. We then propose a memory-efficient framework for free-form ODEs with accurate gradient estimation, which is fundamental to deep-learning models. Furthermore, when applied to invertible blocks, our method achieves constant memory cost with GODE (NODE). We show our method for free-form ODEs generalizes to various model structures and achieves high accuracy for both NODE and GODE in benchmark tasks.

1. INTRODUCTION:

Convolutional neural networks (CNN) have achieved great success in various tasks, such as image classification (He et al., 2016) and segmentation (Long et al., 2015) , video processing (Deng et al., 2014) and machine translation (Sutskever et al., 2014) . However, CNNs are limited to data that can be represented by a grid in the Euclidean domain, such as images (2D grid) and text (1D grid), which hinders their application in irregularly structured datasets.

A graph data structure represents objects as nodes and relations between objects as edges. Graphs are widely used to model irregularly structured data, such as social networks (Kipf & Welling, 2016) , protein interaction networks (Fout et al., 2017) , citation and knowledge graphs (Hamaguchi et al., 2017) . Early works use traditional methods such as random walk (Lovász et al., 1993) , independent component analysis (ICA) (Hyvärinen & Oja, 2000) and graph embedding (Yan et al., 2006) to model graphs, however their performance is inferior due to the low expressive capacity.

Recently a new class of models called graph neural networks (GNN) (Scarselli et al., 2008) were proposed. Inspired by the success of CNNs, researchers generalize convolution operations to graphs to capture the local information. There are mainly two types of methods to perform convolution on a graph: spectral methods and non-spectral methods. Spectral methods typically first compute the graph Laplacian, then perform filtering in the spectral domain (Bruna et al., 2013) . Other methods aim to approximate the filters without computing the graph Laplacian for faster speed (Defferrard et al., 2016) . For non-spectral methods, the convolution operation is directly performed in the graph domain, aggregating information only from the neighbors of a node (Duvenaud et al., 2015; Atwood & Towsley, 2016) . The recently proposed GraphSAGE (Hamilton et al., 2017 ) learns a convolution kernel in an inductive manner.

To our knowledge, all existing GNN models mentioned above have a structure of discrete layers. The discrete structure makes it hard for the GNN to model continuous diffusion processes (Freidlin & Wentzell, 1993; Kondor & Lafferty, 2002) in graphs. The recently proposed neural ordinary differential equation (NODE) (Chen et al., 2018) views a neural network as an ordinary differential equation (ODE), whose derivative is parameterized by the network, and the output is the solution to this ODE. We extend NODE from the Euclidean domain to graphs and propose graph ordinary differential equations (GODE) , where the message propagation on a graph is modeled as an ODE.

NODEs are typically trained with adjoint method. NODEs have the advantages of adaptive evaluation, accuracy-speed control by changing error tolerance, and are free-form continuous invertible models (Chen et al., 2018; Grathwohl et al., 2018) . However, to our knowledge, in benchmark image classification tasks, NODEs are significantly inferior to state-of-the-art discrete-layer models (error rate: 19% for NODE vs 7% for ResNet18 on CIFAR10) (Dupont et al., 2019; Gholami et al., 2019) . In this work, we show this is caused by error in gradient estimation during training of NODE, and propose a memory-efficient framework for accurate gradient estimation. We demonstrate our framework for free-form ODEs generalizes to various model structures, and achieves high accuracy for both NODE and GODE in benchmark tasks. Our contribution can be summarized as follows:

1. We propose a framework for free-form NODEs to accurately estimate the gradient, which is fundamental to deep-learning models. Our method significantly improves the performance on benchmark classification (reduces test error from 19% to 5% on CIFAR10). 2. Our framework is memory-efficient for free-form ODEs. When applied to restricted-form invertible blocks, the model achieves constant memory usage. 3. We generalize ODE to graph data and propose GODE models. 4. We demonstrate improved performance on different graph models and various datasets.

2. RELATED WORKS:



2.1. NEURAL NETWORKS AND DIFFERENTIAL EQUATIONS:

There have been efforts to view neural networks as differential equations. Lu (2017) viewed a residual network as a discretization of a differential equation and proposed several new architectures based on numerical methods in ODE solver. Haber & Ruthotto (2017) proposed a stable architecture based on analysis of the ODE. Chen et al. (2018) proposed neural ordinary differential equation (NODE), which treats the neural network as a continuous ODE. NODE was later used in a continuous normalizing flow for generative models (Grathwohl et al., 2018) .

There have been many studies on the training of NODE. The adjoint method has long been widely used in optimal control (Stapor et al., 2018) and geophysical problems (Plessix, 2006) , and recently applied to ODE (Chen et al., 2018) . Dupont et al. (2019) proposed augmented neural ODEs to improve the expressive capacity of NODEs. However, to our knowledge, none of the methods above discusses the inaccurate gradient estimation issue; empirical performances of NODE in benchmark classification tasks are significantly inferior to state-of-the-art discrete-layer models.

2.2. GRAPH NEURAL NETWORKS:

GNNs can be divided into two categories: spectral methods and non-spectral methods. Spectral GNNs perform filtering in the Fourier domain of a graph, thus need information of the whole graph to determine the graph Laplacian. In contrast, non-spectral GNNs only consider message aggregation around neighbor nodes, therefore are localized and generally require less computation (Zhou et al., 2018) .

We first briefly introduce several spectral methods. Bruna et al. (2013) first introduced graph convolution in the Fourier domain based on the graph Laplacian, however the computation burden is heavy because of non-localized filters. Henaff et al. (2015) incorporated a graph estimation procedure in spectral networks and parameterized spectral filters into a localized version with smooth coefficients. Defferrard et al. (2016) used Chebyshev expansion to approximate the filters without the need to compute the graph Laplacian and its eigenvectors, therefore significantly accelerated the running speed. Kipf & Welling (2016) proposed to use a localized first-order approximation of graph convolution on graph data and achieved superior performance in semi-supervised tasks for node classification. Defferrard et al. (2016) proposed fast localized spectral filtering on graphs. graphs. GraphSAGE (Hamilton et al., 2017) samples a fixed size of neighbors for each node for fast localized inference. Graph attention networks (Veličković et al., 2017) learn different weights for different neighbors of a node. The graph isomorphism network (GIN) (Xu et al., 2018) has a structure as expressive as the Weisfeiler-Lehman graph isomorphism test.

2.3. INVERTIBLE BLOCKS:

Invertible blocks are a family of neural network blocks whose forward function is a bijective mapping. Therefore, the input to a bijective block can be accurately reconstructed from its outputs. Invertible blocks have been used in normalizing flow (Rezende & Mohamed, 2015; Dinh, 2016; Kingma & Dhariwal, 2018; Dinh et al., 2014; Kingma et al., 2016) , where the model is required to be invertible in order to calculate the log-density of data distribution. Later on, Jacobsen et al. (2018) used bijective blocks to build invertible networks. Gomez et al. (2017) proposed to use invertible blocks to perform back propagation without storing activation, which achieves a memory-efficient network structure. They were able to discard activation of middle layers, because each layer's activation can be reconstructed from the next layer with invertible blocks.

3. ACCURATE GRADIENT ESTIMATION FOR TRAINING OF NODE:



3.1. FROM DISCRETE MODELS TO CONTINUOUS MODELS:

We first consider discrete-layer models with residual connection (He et al., 2016) , which can be represented as:

EQUATION

) where x k is the states in the kth layer; f k (•) is any differentiable function whose output has the same shape as its input.

When we add more layers with shared weights, and let the stepsize in Eq. 1 go to infinitesimal, the difference equation turns into a neural ordinary differential equation (NODE) (Chen et al., 2018) :

EQUATION

We use z(t) in the continuous case and x k in the discrete case to represent hidden states. f (•) is the derivative parameterized by a network. Note that a key difference between Eq. 1 and 2 is the form of f : in the discrete case, different layers (different k values) have their own function f k ; while in the continuous case, f is shared across all time t.

The forward pass of model with discrete layers can be written as:

EQUATION

) where K is the total number of layers. Then an output layer (e.g. fully-connected layer for classification) is applied on x K .

The forward pass of a NODE is:

EQUATION

where z(0) = input and T is the integration time, corresponding to number of layers K in the discrete case. The transformation of states z is modeled as the solution to the NODE. Then an output layer is applied on z(T ). Integration in the forward pass can be performed with any ODE solver, such as the Euler Method, Runge-Kutta Method, VODE solver and Dopris Solver (Milne & Milne, 1953; Brown et al., 1989; Ascher et al., 1997) .

3.2. BACK-PROP WITH ADJOINT METHOD IS SENSITIVE TO NUMERICAL ERROR:

The adjoint method is widely used in optimal process control and functional analysis (Stapor et al., 2018; Pontryagin, 2018) . We follow the method by (Chen et al., 2018) . Denote model parameters as θ, which is independent of time. Define the adjoint as: where L is the loss function. Then we have

EQUATION

EQUATION

with detailed proof from optimization perspective in appendix F. Then we can perform gradient descent to optimize θ to minimize L. Eq. 6 is a reverse-time integration, which can be solved with any ODE solver (Chen et al., 2018) . To evaluate ∂f (z(t),t,θ) ∂θ , we need to determine z(t) by solving Eq. 2 reverse-time (Directly storing z(t) during forward pass requires a large memory consumption, because the continuous model is equivalent to an infinite-layer model). To summarize, in the forward pass we solve Eq. 2 forward in time; in the backward pass, we solve Eq. 2 and 6 reverse in time, with initial condition determined from Eq. 5 at time T .

We give an intuition why the reverse-time ODE solver causes inaccurate gradient in adjoint methods. The backward pass (Eq. 6) requires determining f (z(t), t, θ) and ∂f (z(t),t,θ) ∂θ , which requires determining z(t) by solving Eq. 2 reverse-time. As shown in Fig. 1 (a, b ), the hidden state solved forward-time (z(t i )) and the hidden state solved reverse-time (h(t i )) may not be equal; this could be caused by the instability of reverse-time ODE, and is represented by the mismatch between z(t) (dashed curve) and h(t) (solid curve). Error h(t) -z(t) will cause error in gradient dL dθ . Proposition 1 For an ODE in the form dz (t) dt = f (z(t), t), denote the Jacobian of f as J f . If this ODE is stable both in forward-time and reverse-time, then Re(λ i (J f )) = 0 ∀i, where λ i (J f ) is the ith eigenvalue of J f , and Re(λ) is the real part of λ.

Detailed proof is in appendix C. Proposition 1 indicates that if the Jacobian of the original system Eq. 2 has eigenvalues whose real-part are not 0, then either the reverse-time or forward-time ODE is unstable. When |Re(λ)| is large, either forward-time or reverse-time ODE is sensitive to numerical errors. This phenomenon is also addressed in Chang et al. (2018) . This instability affects the accuracy of solution to Eq. 2 and 6, thus affects the accuracy of the computed gradient.

3.3. MEMORY-EFFICIENT DIRECT BACK-PROPAGATION THROUGH ODE SOLVER:

The adjoint method might be sensitive to numerical errors when solving the ODE in reverse-time.

To resolve this, we propose to directly back-propagate through the ODE solver.

As in Fig. 1 (a), the ODE solver uses discretization for numerical integration, evaluated at time points {t 0 , t 1 , ...t N }. Fig. 1(c ) demonstrates the direct back-propagation with accurate hidden states h(t i ), which can be achieved with two methods: (1) the activation z(t i ) can be saved in cache for back-prop, but requires huge memory; or (2) we can accurately reconstruct z(t i ) by re-building the computation graph directly at evaluated time points {t i }. Since the model is evaluated at the same time points t i in forward-time, it's guaranteed that z(t i ) = h(t i ). Therefore direct back-prop is accurate, regardless of the stability of Eq. 2.

Similar to the continuous case, we can define the adjoint with discrete time. Then we have:

a i = ∂L ∂z(t i ) , a i+1 = a i ∂z(t i+1 ) ∂z(t i ) , dL dθ = N i=1 a i ∂z(t i ) ∂θ (7)

where a i is the adjoint for the ith step in discrete forward-time ODE solution. Eq. 7 can be viewed as a numerical discretization of Eq. 6. We show Eq. 6 can be derived from an optimization perspective. Detailed derivations of Eq. 6-7 are in appendix E and F.

Algorithm 1: Algorithm for accurate gradient estimation in ODE solver for free-form functions Define model dz(t) dt = f (z(t), t), where f is a free-form function. Denote integration time as T . Forward (f, T, z 0 , tolerance): 

t = 0, z = z 0 state 0 = f.state dict(),

z(t i+1 ) = step(f, z(t i ), t i , step = t i+1 -t i ) For t i in {t N , t N -1 , ..., t 1 , t 0 } : Determine a(t i ) = ∂L ∂z(ti) and ∂z(ti) ∂θ dL dθ = N i=1 a i ∂z(ti) ∂θ

return dL dθ Details of our method are summarized in Algorithm 1. We discuss its properties below:

Summary of the algorithm (1) During forward pass, the solver performs a numerical integration, with the stepsize adaptively varying with error estimation. (2) During forward pass, the solver outputs the integrated value, and the evaluation time points {t i }. All middle activations are deleted to save memory. (3) During backward pass, the solver re-builds the computation graph, by directly evaluating at saved time points, without adaptive searching. (4) During backward pass, the solver performs a numerical version (Eq. 7) of reverse-time integration (Eq. 6).

Support for free-form continuous dynamics:

There's no constraint on the form of f . Therefore, our algorithm is a generic method.

Memory consumption analysis (1) Suppose f has N f layers, the number of forward evaluation step is N t on average, and the evaluations to adaptively search for an optimal stepsize is K. A naive solver will take O(N f × N t × K), while our method consumes O(N f × N t ) because all middle activations are deleted during forward pass, and we don't need to search for optimal stepsize in backward pass.

(2) If we perform step-wise checkpoint method, where we only store z(t i ) for all t i , and compute the gradient ∂z(ti+1) ∂z(ti) for one t i at a time, then the memory consumption can be reduced to O(N f +N t ). (3) Since the solver can handle free-form functions, it can also handle restricted form invertible block (see below). In this case, we don't need to store z(t i ), and the memory consumption can reduce to O(N f ).

More memory-efficient solver with invertible blocks Restricting the form of f to invertible blocks (Gomez et al., 2017) allows for O(N f ) memory consumption. For invertible blocks, input x is split into two parts (x 1 , x 2 ) of the same size (e.g. x has shape N × C, where N is batch size, C is channel number; we can split x into x 1 and x 2 with shape N × C

2 ). The forward and inverse of a bijective block can be denoted as:

   y 2 = ψ x 2 , F (x 1 ) y 1 = ψ x 1 , G(y 2 )    x 1 = ψ -1 y 1 , G(y 2 ) x 2 = ψ -1 y 2 , F (x 1 ) (8)

where the output of a bijective block is denoted as (y 1 , y 2 ) with the same size as (x 1 , x 2 ). F and G are any differentiable neural networks, whose output has the same shape as the input. ψ(α, β) is a differentiable bijective function w.r.t α when β is given; ψ -1 (α, β) is the inverse function of ψ.

Theorem 1 If ψ(α, β) is a bijective function w.r.t α when β is given, then the block defined by Eq. 8 is a bijective mapping.

Proof of Theorem 1 is given in appendix D. Based on this, we can apply different ψ functions for different tasks. Since x can be accurately reconstructed from y, there's no need to store activations, hence is memory-efficient. Details for back-prop without storing activation are in appendix B.

4. NEURAL ORDINARY DIFFERENTIAL EQUATION ON GRAPH NETWORKS:

We first introduce graph neural networks with discrete layers, then extend to the continuous case and introduce graph ordinary differential equations (GODE).

4.1. MESSAGE PASSING IN GNN:

As shown in Fig. 2 , a graph is represented with nodes (marked with circles) and edges (solid lines). We assign a unique color to each node for ease of visualization. Current GNNs can generally be represented in a message passing scheme (Fey & Lenssen, 2019) :

EQUATION

)

EQUATION

EQUATION

where x u k represents states of the uth node in the graph at kth layer and e u,v represents the edge between nodes u and v. N (u) represents the set of neighbor nodes for node u. ζ represents a differentiable, permutation invariant operation such as mean, max or sum. γ (k) and φ (k) are differentiable functions parameterized by neural networks.

For a specific node u, a GNN can be viewed as a 3-stage model, corresponding to Eq. 9-11:

(1) Message passing, where neighbor nodes v ∈ N (u) send information to node u, denoted by message (v,u) . The message is generated from function φ(•), parameterized by a neural network.

(2) Message aggregation, where a node u aggregates all messages from its neighbors N (u), denoted as aggregation u . The aggregation function ζ is typically permutation invariant operations such as mean and sum, because graphs are invariant to permutation. (3) Update, where the states of a node are updated according to its original states x u k-1 and aggregation of messages aggregation u , denoted as γ(•).

4.2. CONTINUOUS-TIME MODELS ON GRAPHS:

We can convert a discrete-time GNN to continuous-time GNN by replacing f in Eq. 2 with the message passing process defined in Eq. 9 to 11, which we call graph ordinary differential equation (GODE). A diagram of GODE is shown in Fig. 2 . Because GODE is an ODE in nature, it can capture highly non-linear functions, thus has the potential to outperform its discrete-layer counterparts.

We demonstrate that the asymptotic stability of GODE could be related to the over-smoothing phenomena (Li et al., 2018 ). It's demonstrated that graph convolution is a special case of Laplacian smoothing (Li et al., 2018) , which can be written as Y = (I -γ D-1/2 L D-1/2 )X where X and Y are the input and output of a graph-conv layer respectively, Ã = A + I where A is the adjacency matrix, and D is the corresponding degree matrix of Ã, and γ is a positive scaling constant.

When modified from a discrete model to a continuous model, the continuous smoothing process is:

dX dt = -γ D-1/2 L D-1/2 X (12)

Since all eigenvalues of the symmetrically normalized Laplacian are real and non-negative, then all eigenvalues of the above ODE are real and non-positive. Suppose all eigenvalues of the normalized Laplacian are non-zero. In this case, the ODE has only negative eigenvalues, hence the ODE above is asymptotically stable (Lyapunov, 1992). Hence as time t grows sufficiently large, all trajectories are close enough. In the experiments, this suggests if integration time T is large enough, all nodes (from different classes) will have very similar features, thus the classification accuracy will drop. 

5. EXPERIMENTS:



5.1. DATASETS:

To evaluate our method on general NODE, we conducted experiments with a CNN-NODE on two benchmark image classification tasks (CIFAR10 and CIFAR100) (Krizhevsky et al., 2009) .

We also evaluated our method on benchmark graph datasets, including 2 bioinformatic graph classification datasets (MUTAG and PROTEINS), 2 social network graph classification datasets (IMDB-BINRAY, REDDIT-BINARY) (Yanardag & Vishwanathan, 2015) , and 3 citation networks (Cora, CiteSeer and PubMed). For graph classification tasks, different from the experiment settings in Xu et al. (2018) , we input the raw dataset into our models without pre-processing. For node classification tasks, we performed transductive inference and strictly followed the train-validation-test split by Kipf & Welling (2016) , where less than 6% nodes are used as training examples. Details of datasets are summarized in appendix A.

5.2. MODEL STRUCTURES:

For image classification tasks, we directly modify a ResNet18 into its corresponding NODE model. For each block, the function is dz(t) dt = f (z(t)) where f is a sequence of conv -bn -reluconv -bn -relu layers. f is the same as residual branch in ResNet, and it can be replaced with any free-form functions. For tasks on graph datasets, GODE can be applied to any graph neural network by simply replacing f in Eq. 2 with corresponding structures (free-form functions), or replacing F, G in Eq. 8 with other structures (invertible blocks). To demonstrate that GODE is easily generalized to existing structures, we used several different GNN architectures, including the graph convolutional network (GCN) (Kipf & Welling, 2016), graph attention network (GAT) (Veličković et al., 2017) , graph network approximated with Chebyshev expansion (ChebNet) (Defferrard et al., 2016) , and graph isomorphism network (GIN) (Xu et al., 2018) . For a fair comparison, we trained GNNs with different depths of layers (1-3 middle layers, besides an initial layer to transform data into specified channels, and a final layer to generate prediction), and reported the best results among all depths for each model structure.

On the same task, different models use the same hyper-parameters on model structures, such as channel number. For graph classification tasks, we set the channel number of hidden layers as 32 for all models; for ChebNet, we set the number of hops as 16. For node classification tasks, we set the channel number as 16 for GCN and ChebNet, and set number of hops as 3 for ChebNet; for GAT, we used 8 heads, and set each head as 8 channels.

For every GNN structure, we experimented with different number of hidden layers (1,2,3), and calculated the mean and variance of accuracy of 10 runs.

5.3. COMPARISON OF BACK-PROPAGATION METHODS:

We compared the adjoint method and direct back-propagation on the same network, and demonstrated direct back-prop generates higher accuracy. For CNN-NODE on classification tasks, we directly modify a ResNet18 into NODE18, and report resuls in Table . 1; for graph networks, we train a GODE model with a GCN to parameterize the derivative, and report results in Table . 2.

Empirical performance Direct back-propagation consistently outperformed the adjoint method for both tasks. This result validates our analysis on the instability of the adjoint method, which is intuitively caused by the instability of the reverse-time ODE. On image classification tasks, compared to adjoint method, our training method reduces error rate of NODE18 from 19% (37%) to 5%(23%) on CIFAR10 (CIFAR100). Furthermore, NODE18 has the same number of parameters as ResNet18, but outperforms deeper networks such as ResNet101 on both datasets. Our method also consistently outperforms the adjoint method on several benchmark graph datasets, as shown in 

5.4. GENERAL BIJECTIVE BLOCKS:

We demonstrate that bijective blocks defined as Eq. 8 can be easily generalized: F and G are general neural networks, which can be adapted to different tasks; ψ(α, β) can be any differentiable bijective mapping w.r.t. α when β is given. We demonstrate two examples of ψ:

(1) additive, forward is η = ψ(α, β) = α + β, inverse is α = ψ -1 (η, β) = η -β; (2) linear sigmoid, forward is η = ψ(α, β) = α × sigmoid(β), inverse is α = ψ -1 (η, β) = η/sigmoid(β).

Results for different ψ are reported in Table 3 . Note that we experimented with different depths and reported the best accuracy for each model, and performed a paired t-test on results from GODE and their discrete-layer counterparts. Most GODE models outperformed their corresponding discretelayer models significantly, validating the effectiveness of GODE; different ψ functions behaved similarly on our node classification tasks, indicating the continuous-time model is more important than coupling function ψ. We also validate the lower memory cost, with details in appendix B.

5.5. RESULTS ON GRAPH CLASSIFICATION TASK:

Results for different models on graph classification tasks are summarized in Table 4 . We experimented with different structures, including GCN, ChebNet and GIN; for corresponding GODE models (marked with ODE), we tested both free-form (marked with "free") and invertible block (marked with "INV"). We performed paired t-test comparing GODE and its discrete-layer counterparts. For most experiments, GODE models performed significantly better. This indicates the continuous process model might be important for graph models.

5.6. IMPACT OF INTEGRATION TIME:

For a NODE and GODE model, during inference, we test the influence of integration time. Results are summarized in Table . 5. When integration time is short, the network does not gather sufficient information from neighbors; when integration time is too long, the model is sensitive to over-smooth issue, as discussed in Sec. 4.2. We observe accuracy drop in both cases.

6. CONCLUSIONS:

We propose GODE, which enables us to model continuous diffusion process on graphs. We propose a memory-efficient direct back-propagation method to accurately determine the gradient for general free-form NODEs, and validate its superior performance on both image classification tasks and graph data. Furthermore, we related the over-smoothing of GNN to asymptotic stability of ODE.

Our paper tackles the fundamental problem of gradient estimation for NODE; to our knowledge, it's the first paper to improve accuracy on benchmark tasks to comparable with state-of-the-art discrete layer models. It's an important step to apply NODE from theory to practice.

A DATASETS:

We perform experiments on various datasets, including citation networks (Cora, CiteSeer, PubMed), social networks (COLLAB, IMDB-BINARY, REDDIT-BINARY), and bioinformatics datasets (MUTAG, PROTEINS). Details of each dataset are summarized in Table 1 . We explain the structure and conduct experiments for the invertible block here.

Structure of invertible blocks Structure of invertible blocks are shown in Fig. 1 . We follow the work of Gomez et al. (2017) with two important modifications: (1) We generalize to a family of bijective blocks with different ψ in Eq. 8 in the main paper, while Gomez et al. (2017) restrict the form of ψ to be sum.

(2) We propose a parameter state checkpoint method, which enables bijective blocks to be called more than once, while still generating accurate inversion.

The algorithm is summarized in Algo. 2. We write the pseudo code for forward and backward function as in PyTorch. Note that we use "inversion" to represent reconstructing input from the output, and use "backward" to denote calculation of the gradient. To reduce memory consumption, in the forward function, we only keep the outputs y 1 , y 2 and delete all other variables and computation graphs. In the backward function, we first "inverse" the block to calculate x 1 , x 2 from y 1 , y 2 , then perform a local forward and calculate the gradient ∂[y1,y2] ∂[x1,x2] .

Experiments In this section we demonstrate that our bijective block is memory efficient. We trained a GODE model with bijective blocks, and compared the memory consumption using our memory-efficient function as in Algo. 2 and a memory-inefficient method as in conventional backpropagation. Results were measured with a batchsize of 100 on MUTAG dataset. 

∂L ∂[x 1 ,x 2 ] = ∂L ∂[y 1 ,y 2 ] ∂[Y 1 ,Y 2 ] ∂[X 1 ,X 2 ] ∂L ∂[θ F ,θ G ] = ∂L ∂[y 1 ,y 2 ] ∂[Y 1 ,Y 2 ] ∂[θ F ,θ G ] delete Y1, Y2, X1, X2 return ∂L/∂[x1, x2], ∂L/∂[θF , θG] C PROOF FOR PROPOSITION 1

Proposition 1 For an ODE in the form dz(t) dt = f (z(t), t), denote the Jacobian of f as J f . If this ODE is stable both in forward-time and reverse-time, then Re(λ i (J f )) = 0 ∀i, where λ i (J f ) is the ith eigenvalue of J f , and Re(λ) is the real part of λ.

Proof Denote s = T -t, where T is the end time. Notice that the reverse-time in t is equivalent to forward-time in s. Therefore, we have forward-time ODE:

EQUATION

and reverse-time ODE:

EQUATION

Therefore, we have λ(J f ) = -λ(J g ). For both forward-time and reverse-time ODE to be stable, the eigenvalues of J need to have non-positive real part. Therefore

EQUATION

The only solution is Reλ

EQUATION

D PROOF FOR THEOREM 1:

Theorem 1 For bijective block whose forward and reverse mappings are defined as

F orward(x 1 , x 2 ) =    y 2 = ψ x 2 , F (x 1 ) y 1 = ψ x 1 , G(y 2 ) Reverse(y 1 , y 2 ) =    x 1 = ψ -1 y 1 , G(y 2 ) x 2 = ψ -1 y 2 , F (x 1 )

If ψ(α, β) is a bijective function w.r.t α when β is given, then the block is a bijective mapping.

Proof To prove the forward mapping is bijective, it is equivalent to prove the mapping is both injective and surjective.

Injective We need to prove, if F orward(x 1 , x 2 ) = F orward(x 3 , x 4 ), then x 1 = x 3 , x 2 = x 4 .

The assumption above is equivalent to

EQUATION

Since ψ(α, β) is bijective w.r.t α when β is given, from Eq.( 6), we have x 1 = x 3 . Similarly, condition on x 1 = x 3 and Eq.( 5), using bijective property of ψ, we have x 2 = x 4 . Therefore, the mapping is injective.

Surjective We need to prove ∀ [y 1 , y 2 ], ∃ [x 1 , x 2 ] s.t. F orward(x 1 , x 2 ) = [y 1 , y 2 ].

Given y 1 , y 2 , we construct

EQUATION

Then for the forward function, given bijective property of ψ, apply F orward and Reverse defined in the proposition statement,

z 2 = ψ(x 2 , F (x 1 )) = ψ ψ -1 y 2 , F (x 1 ) , F (x 1 ) = y 2 (8) z 1 = ψ(x 1 , G(y 2 )) = ψ ψ -1 y 1 , G(y 2 ) , G(y 2 ) = y 1 (9) Therefore we construct x 1 , x 2 s.t. F orward(x 1 , x 2 ) = [y 1 , y 2 ].

Therefore the mapping is surjective. Therefore is bijective.

E DERIVATION OF GRADIENT IN DISCRETE CASE:

We use a figure to demonstrate the computation graph, and derive the gradient from the computation graph.

The loss is L, forward pass is denoted with black arrows, gradient back-propagation is shown with red arrows. We use p to denote each path from θ to L, corresponding to all paths in red that goes from L to θ. In this section we derive the gradient of parameters in an neural-ODE model from an optimization perspective. Then we extend from continuous cases to discrete cases.

Notations With the same notations as in the main paper, we use z(t) to denote hidden states z at time t. Denote parameters as θ, and input as x, target as y, and predicted output as ŷ. Denote the loss as J(ŷ, y). Denote the integration time as 0 to T .

Problem setup The continuous model is defined to follow an ODE: dz(t) dt = f (z(t), t, θ), s.t. z(0) = x (12)

We assume f is differentiable, since f is represented by a neural network in our case. The forward pass is defined as:

EQUATION

The loss function is defined as: J(ŷ, y) = J(z(T ), y) (14) We formulate the training process as an optimization problem:

EQUATION

For simplicity, Eq. 15 only considers one ODE block. In the case of multiple blocks, z(T ) is the input to the next ODE block. As long as we can derive dLoss dθ and dLoss dz(0) when dLoss dz(T ) is given, the same analysis here can be applied to the case with a chain of ODE blocks.

Lagrangian Multiplier Method We use the Lagrangian Multiplier Method to solve the problem defined in Eq. 15. For simplicity, only consider one example (can be easily extended to multiple examples cases), the Lagrangian is

EQUATION

Karush-Kuhn-Tucker (KKT) conditions are necessary conditions for an solution to be optimal. In the following sections we start from the KKT condition and derive our results.

Derivative w.r.t. λ At optimal point, we have δL δλ = 0. Note that λ is a function of t, we derive the derivative from calculus of variation.

Consider a cotninuous and differentiable perturbation λ(t) on λ(t), and a scalar , L now becomes a function of , L( ) = J(z(0) + T 0 f (z(t), t, θ), y) + At optimal λ(t), dL d | =0 = 0 for all continuous differentiable λ(t). Therefore, dz(t) dt -f (z(t), t, θ) = 0, ∀t ∈ (0, T )

From continuous to discrete case To derive corresponding results in discrete cases, we need to replace all integration with finite sum.

In discrete cases, the ODE condition turns into:

EQUATION

from Eq. 31, we can get:

EQUATION

Re-arranging terms we have:

[(- ∂L ∂z i+1

) -(-

∂L ∂z i )]/[t i+1 -t i ] + (- ∂L ∂z i+1 ) ∂f (z i , t i , θ) ∂z i = 0 (33)

which is the discrete version of Eq. 29. Which also corresponds to our analysis in Eq. 10 and 11.

:



Figures and Tables:

Figure 1: Comparison of two methods for back-propagation on NODE. As in figure (a), the ODE solver is discretized at points {t0, t1, ..., tN } during forward pass. Black dashed curve shows hidden state solved in forward-time, denoted as z(t).Figure (b)shows the adjoint method, red solid line shows the hidden state solved in reverse-time, denoted as h(t). Ideally z(t) = h(t) and dashed curve overlaps with solid curve; however, the reverse-time solution could be numerically unstable, and causes z(t) = h(t), thus causes error in gradient. Figure(c) shows the direct back-propagation through ODE solver. In direct back-propagation, we save evaluation time points {t0, t1, ...tN } during forward pass; during backward pass, we re-build the computation graph by directly evaluating at the same time points. In this way, z(ti) = h(ti). Since the hidden state can be accurately reconstructed, the gradient can be accurately evaluated.

Figure (b)  shows the adjoint method, red solid line shows the hidden state solved in reverse-time, denoted as h(t). Ideally z(t) = h(t) and dashed curve overlaps with solid curve; however, the reverse-time solution could be numerically unstable, and causes z(t) = h(t), thus causes error in gradient. Figure(c) shows the direct back-propagation through ODE solver. In direct back-propagation, we save evaluation time points {t0, t1, ...tN } during forward pass; during backward pass, we re-build the computation graph by directly evaluating at the same time points. In this way, z(ti) = h(ti). Since the hidden state can be accurately reconstructed, the gradient can be accurately evaluated.

Figure 2: Diecrete-time and continuous-time models on a graph. Nodes are represented with circles, and each node is represented with a unique color. Edges are represented with solid lines. For discrete-time models in (a), the hidden states of nodes are updated with discrete steps. For continuous-time models in (b), hidden states of each node evolves continuously with time. The dynamics of nodes are represented with dashed lines, with the same color as corresponding nodes.

Figure 1: Structure of bijective blocks. F and G can be any differentiable neural network whose output has the same shape as its input. Blue dot (Orange diamond) represents the forward (inverse) of a bijective function, corresponding to ψ (ψ -1 ) in Eq. 8 of the main paper. Left (right) figure represents the forward (inverse) as in Eq. 8.

easy to check the conditions for Leibniz integral rule, and we can switch integral and differentiationdt -f (z(t), t, θ)]dt (18)

cache.save(state 0 ) Select initial step size h = h 0 (adaptively with adaptive step-size solver).

Table content: <table/>

Error rate on test set. "Ours" represent NODE models directly modified from ResNet18, trained with Heun-Euler solver, but tested with different solvers. "Adjoint" is the result when trained and tested with adjoint method, reported by(Gholami et al., 2019). We also report results from standard ResNet.

Table content: <table><tr><td/><td/><td/><td/><td>Ours</td><td/><td/><td/><td/><td colspan="2">Literature</td></tr><tr><td/><td/><td colspan="6">Adaptive Stepsize Solvers Heun-Euler RK23 RK45 Euler RK2 Fixed stepsize solvers RK4</td><td colspan="4">adjoint-ODE ResNet18 ResNet50 ResNet101</td></tr><tr><td colspan="2">CIFAR10</td><td>4.85</td><td>4.92</td><td>5.29</td><td>5.52</td><td>5.27</td><td>5.24</td><td>19.2</td><td>6.98</td><td>6.38</td><td>6.25</td></tr><tr><td colspan="2">CIFAR100</td><td>22.66</td><td colspan="5">24.13 23.56 24.44 24.44 24.43</td><td>37.6</td><td>27.08</td><td>25.73</td><td>24.84</td></tr><tr><td/><td colspan="2">MUTAG</td><td colspan="2">PROTEINS IMDB</td><td/><td>REDDIT</td><td/><td/><td/><td/></tr><tr><td colspan="4">adjoint 68.1± 4.6 67.0±3.7</td><td colspan="3">72.1±0.4 69.5±5.9</td><td/><td/><td/><td/></tr><tr><td>ours</td><td colspan="3">80.8±8.3 73.9±3.1</td><td colspan="3">74.6±5.1 92.4±2.1</td><td/><td/><td/><td/></tr></table>

Accuracy of adjoint method and direct back-prop, for a GODE model with GCN as the derivative function.

Table content: <table/>

Table. 2.

Table content: <table><tr><td>Model</td><td/><td>ψ</td><td>Cora</td><td>CiteSeer</td><td>PubMed</td><td>Model</td><td/><td>MUTAG</td><td>PROTEINS IMDB</td><td>REDDIT</td></tr><tr><td>GCN</td><td>DISC ODE</td><td colspan="4">81.6±0.5 81.7±0.7 l sig 81.8±0.3** 72.4±0.8** 80.1±0.3** 71.6±0.3 79.2±0.1 add 72.4±0.6** 80.0±0.2**</td><td>GCN</td><td>DISC ODE</td><td colspan="2">73.3±5.2 INV 78.1±6.2** 74.7±4.3** 75.3±5.3* 72.4±3.1 74.2±3.6 free 75.1±5.3** 76.6±3.9** 73.9±4.6</td><td>85.9±1.8 89.2±3.2** 88.5±3.0**</td></tr><tr><td/><td>DISC</td><td/><td>82.9±0.3</td><td>71.7±0.8</td><td>78.9±0.3</td><td/><td>DISC</td><td>84.0±6.4</td><td>70.6±3.9</td><td>71.9±3.8</td><td>91.0±1.5</td></tr><tr><td>GAT</td><td>ODE</td><td colspan="4">add l sig 83.1±0.4* 83.3±0.3** 72.1±0.6** 79.1±0.5* 72.1±0.3** 79.0±0.5</td><td>Cheb</td><td>ODE</td><td>INV 85.0±8.3 free 86.1±6.3*</td><td>72.7±3.7** 72.0±5.0* 72.5±4.7** 73.6±4.0** 92.4±1.6** 91.2±1.5</td></tr><tr><td>Cheb</td><td>DISC ODE</td><td colspan="2">82.1±0.5 82.4±0.5* l sig 82.2±0.4* add</td><td colspan="2">70.8±0.5 71.1±0.5** 77.8±1.2** 76.6±0.8 70.8±0.6 77.0±1.1*</td><td>GIN</td><td>DISC ODE</td><td colspan="2">85.0±6.4 INV 87.7±5.2* free 87.8±5.1** 73.1±3.7 73.0±3.1 74.5±3.5*</td><td>73.3±5.1 75.4±3.9** 90.5±1.5** 89.2±2.5 74.8±4.1** 91.0±2.4**</td></tr><tr><td colspan="6">Table 3: Results on node classification tasks. We compared various discrete-layer structures (marked with DISC) and their corresponding GODE models (marked with ODE). We tested GODE model with dif-ferent ψ functions ("l sig" represents linear sigmoid). For each model, we use * (**) to mark GODE models that outperform corresponding discrete-layer baselines at a 5% (1%) significance level under paired t-test.</td><td colspan="4">Table 4: Results on graph-classification tasks. For each base model structure, discerete-layer model is marked with DISC; for corresponding GODE, we tested both free-form functions ("free"), and their in-vertible block form ("INV"). For each model, we use * (**) to mark GODE models that outperform corre-sponding discrete-layer baselines at a 5% (1%) signif-icance level under paired t-test.</td></tr><tr><td colspan="7">Integration time 0.5 Cora 80.5 81.6 80.1 80.1 79.3 77.6 68.9 35.1 1.0 1.5 2.0 5.0 10.0 20.0 100.0 CIFAR10 91.3 95.2 94.2 88.4 10.0 10.0 --</td><td colspan="3">Table 5: Accuracy of a free-form GCN-ODE on Cora and NODE18 on CIFAR10, varying with integration time.</td></tr><tr><td colspan="10">Robustness to ODE solvers We implemented adaptive ODE solvers of different orders, as shown</td></tr><tr><td colspan="10">in Table 1. HeunEuler, RK23, RK45 are of order 1, 2, 4 respectively, i.e., for each step forward in</td></tr><tr><td colspan="10">time f is evaluated 1, 2, 4 times respectively. During inference, using different solvers is equivalent</td></tr></table>

Statistics of datasets

Table content: <table><tr><td>Dataset</td><td colspan="3">Graphs Nodes Edges</td><td colspan="3">Features Classes Label rate</td></tr><tr><td>Cora</td><td>1</td><td>2,708</td><td>5,278</td><td>1,433</td><td>7</td><td>0.052</td></tr><tr><td>CiteSeer</td><td>1</td><td>3,327</td><td>4,552</td><td>3,703</td><td>6</td><td>0.036</td></tr><tr><td>PubMEd</td><td>1</td><td colspan="3">19,717 44,324 500</td><td>3</td><td>0.003</td></tr><tr><td>MUTAG</td><td>188</td><td>17.93</td><td>19.79</td><td>7</td><td>2</td><td>0.8</td></tr><tr><td>PROTEINS</td><td>1,113</td><td>39.06</td><td>72.82</td><td>3</td><td>2</td><td>0.8</td></tr><tr><td>IMDB-BINARY</td><td>1,000</td><td>19.77</td><td>96.53</td><td>-</td><td>2</td><td>0.8</td></tr><tr><td colspan="2">REDDIT-BINARY 200</td><td colspan="3">429.63 497.76 -</td><td>2</td><td>0.8</td></tr><tr><td colspan="4">B DETAILS ABOUT INVERTIBLE BLOCKS</td><td/><td/><td/></tr></table>

Memory consumption of bijective blocks.Results are summarized in Table.2. We measured the memory consumption with different depths, which is the number of ODE blocks. When depth increases from 10 to 20, the memory by conventional methods increases from 5.3G to 10.5G, while our memory-efficient version only increases from 2.2G to 2.6G. In theory, our bijective block takes O(1) memory, because we only need to store the outputs in cache, while deleting activations of middle layers. For memory-efficient network, the slightly increased memory consumption is because states of F, G need to be cached; but this step takes up minimal memory compared to input data.

Table content: <table><tr><td colspan="3">Depth Memory-efficient Conventional</td><td/></tr><tr><td>10</td><td>2.2G</td><td>5.3G</td><td>"Conventional" represents storing activation of all</td></tr><tr><td>20</td><td>2.6G</td><td>10.5G</td><td>layers in cache, "memory-efficient" represents our method in Algo. 2.</td></tr></table>

Peer review (iclr2020):
Summary The authors discussed that most graph NNs to date considered discrete layers and hence are difficult to model diffusion processes on graphs. This paper proposed the Neural ODE on a graph, termed with a graph ODE, to tackle this problem. The authors gave a sufficient condition under which the adjoint method gets instable and pointed out potential issues of training Neural ODEs using it. To overcome the instability issue, the author proposed to backpropagate errors directly at discretized points. Since the naive implementation of the direct method is memory-consuming, the authors used invertible blocks as building blocks of graph ODEs, which do not store the intermediate activations for backward propagation. Finally, the authors conducted empirical studies to see the effectiveness of the proposed method. Decision I recommend rejecting the paper weakly because I think the extension of Neural ODEs to graphs is straightforward and that the empirical study is not strong enough to compensate for the weakness of the novelty. Theoretical justification of numerical instability of Neural ODEs (Proposition 1) and its empirical verification (Section 5.4) give new insights for understanding Neural ODEs. However, if I understand correctly, the formulation of graph ODEs do not use the internal structures of graph NNs, even the fact that the underlying neural network is a graph NN. Therefore, I think the extension from Neural ODEs to graphs is a bit too straightforward. The authors proposed a method which improves the stability and memory-efficiency of training. We can apply this method to general Neural ODEs, too. In addition, we can attribute the idea of using invertible blocks to existing works (Gomez et al., 2017). Finally, regarding the empirical performance of graph ODEs, the performance improvement from existing GNNs is within the standard deviations. Therefore, I think the empirical result is not sufficiently strong to justify the novelty of applying Neural ODEs to graphs. Taking these things into account, although the authors gave a new result on Neural ODEs, I think the contribution is limited from the viewpoint of the study of graph NNs. Suggestion - As I wrote in the Decision section, the theoretical results are not restricted to graph ODEs but valid for general neural ODEs. Therefore, I think the authors do not have to restrict the application areas to graph ODEs. The possible direction of the paper is to further analysis of training neural ODEs (e.g., instability). On the other hand, if the authors are interested in the extension of neural ODEs to graphs, I expect a more detailed relationship between the neural ODEs framework and underlying GNNs. For example, I am curious how the topological information of graphs affects the graph ODEs via spectral-type graph convolution operations and what is the relationship to the oversmoothing phenomena (Li et al., 2018). - Since Theorem 1 is applicable not only graph ODEs but also Neural ODEs, it implies that ordinal Neural ODEs are also vulnerable to the instability. The experiments in Dupont (2019), which this paper referenced, pointed out the instability of Neural ODEs. I am wondering the proposed method can enhance the training of neural ODEs, too. [Li et al., 2018] Li, Qimai, Zhichao Han, and Xiao-Ming Wu. *Deeper insights into graph convolutional networks for semi-supervised learning.* Thirty-Second AAAI Conference on Artificial Intelligence. 2018.