Characteristic Neural Ordinary Differential Equation
Authors: Xingzi Xu
Ali Hasan
Khalil Elkhalil
Jie Ding
Vahid Tarokh
Venue: ICLR 2023 poster (2023)
Keywords: Neural ODE
Differential Equation
Method of characteristics

Abstract: We propose Characteristic-Neural Ordinary Differential Equations (C-NODEs), a framework for extending Neural Ordinary Differential Equations (NODEs) beyond ODEs. While NODE models the evolution of latent variables as the solution to an ODE, C-NODE models the evolution of the latent variables as the solution of a family of first-order partial differential equations (PDEs) along curves on which the PDEs reduce to ODEs, referred to as characteristic curves. This reduction along characteristic curves allows for analyzing PDEs through standard techniques used for ODEs, in particular the adjoint sensitivity method. We also derive C-NODE-based continuous normalizing flows, which describe the density evolution of latent variables along multiple dimensions. Empirical results demonstrate the improvements provided by the proposed method for irregularly sampled time series prediction on MuJoCo, PhysioNet, and Human Activity datasets and classification and density estimation on CIFAR-10, SVHN, and MNIST datasets given a similar computational budget as the existing NODE methods.
The results also provide empirical evidence that the learned curves improve the system efficiency using a lower number of parameters and function evaluations compared with those of the baselines. 

Published as a conference paper at ICLR 2023
CHARACTERISTIC NEURAL ORDINARY DIFFERENTIAL EQUATIONS
Xingzi Xu 1, 1, Ali Hasan 1. Khalil Elkhalil', Jie Ding', Vahid Tarokh'
Department of Electrical and Computer Engineering, Duke University
2 School of Statistics, University of Minnesota
{xingzi. xu, ali hasan, khalil. elkhalil, vahid. tarokh, @duke.edu di dingieumn.ee
ABSTRACT
We propose Characteristic-Neurr Ordinary Differential Equations (C-NODEs), a framework for extending Neural Ordinary Differential Equations (NODEs) beyond ODEs. While NODE models the evolution of latent variables as the solution to an ODE, C-NODE models the evolution of the latent variables as the solution of a family of first-order partial differential equations (PDEs) along curves on which the PDEs reduce to ODEs, referred to as characteristic curves. This reduction along characteristic curves allows for analyzing PDEs through standard techniques used for ODEs, in particular the adjoint sensitivity method. We also derive C- NODE-based continuous normalizing flows, which describe the density evolution of latent variables along multiple dimensions. Empirical results demonstrate the improvements provided by the proposed method for irregularly sampled time series prediction on MuJoCo, PhysioNet, and Human Activity datasets and classification and density estimation CIFAR-10, 10, SVHN, and MNIST datasets given similar computational budget as the existing NODE methods. The results also provide empirical evidence that the learned curves improve the system efficiency using a lower number of parameters and function evaluations compared with those of the baselines.
INTRODUCTION
Deep learning and differential equations share many connections, and techniques in the intersection have led to insights in both fields. One predominant connection is based on certain neural network architectures resembling numerical integration schemes, leading the development of Neural Ordinary Differential Equations (NODEs) (Chen et al., 2019b). NODEs use neural network parameterization of an ODE to learn a mapping from observed variables to latent variable that is the solution to the learned ODE. A central benefit of NODES is the constant memory cost, when backward passes are computed using the adjoint sensitivity method rather than backpropagating through individual forward solver steps. Backpropagating through adaptive differential equation solvers to train NODEs will often result in extensive memory use, as mentioned in [Chen et al.) (20195). Moreover, NODEs provide flexible probability density representation often referred to continuous normalizing flows (CNFs). However, since NODES can only represent solutions to ODEs, the class of functions is somewhat limited and may not apply more general problems that do have smooth and one-to-one mappings. To address this limitation, a series of analyses based on methods from differential equations have been employed to enhance the representation capabilities of NODEs, such as the theory of controlled differential equations (Kidger et al.) (2020), learning higher-ordee ODEs (Massaroli et al. 2021), augmenting dynamics (Dupont et al., 2019), and considering dynamies with delay terms (Zhu et al., 1021). Additionally, certain works consider generalizing the ODE case to partial differential equations (PDEs), such as in Ruthotto & Haber (2020) and [un et al.) (2019). These PDE-based methods do not use the adjoint method, removing the primary advantage of constant memory cost. This leads us to the central question motivating the work: can we combine the benefits
"Equal contribution.
1 Corresponning Author.

---

Published as a conference paper at ICLR 2023
of the rich function class of PDEs with the efficiency of the adjoint method? To do so, propose method of continuous-depth neural networks that solves a PDE over parametric curves that reduce the PDE to ODE. Such curves are known characteristics, and they define the solution of the PDE in terms of an ODE (Griffiths et al. 2013). The proposed Characteristic Neural Ordinary Differential Equations (C-NODE) learn both the characteristics and the ODE along the characteristics to solve the PDE over the data space. This allows for richer class of models while still incorporating the same memory efficiency of the adjoint method. The proposed C-NODE is also an extension of existing methods, as it improves the empirical accuracy of these methods in classification tasks, time series prediction tasks, and image quality in generation tasks.
RELATED WORK
NODE is often motivated as continuous form of
Residual Network (ResNet) (He et al., 2015), since the forward Buler integration
Single Characteristic (ODE) Multiple Characteristics (PDE)
ResNet can be interpreted scheme on the latent state (Weinan 2017). Specifically, a
ResNet is composed of multiple blocks with each block can be represented as: M++1 Me + f(u,, 0), where the evolving hidden state at time t and f f(ue, 0) is in- terpreted as the gradient at time 6, namely du/ (u)))(ue). Generalizing the model to a step size given by At re- sults in Ut+ ut + f(u,, 0) At. To adapt this model to a continuous setting, we let At - and obtain: lim (U++Ar 4./44 At du(t)/dt. The model can then At-30
Figure 1: Comparison of traditional
NODE (left) and proposed C-NODE (right). The solution to NODE is the solution
be evaluated through existing numerical integration tech- niques, as proposed by Chen et al. (2019b):
single ODE, whereas C- NODE represents series of ODEs that form the solution to a PDE. Each color C-NOE represents the solution to an ODE with a different initial condition. NODE represents a single ODE, and can only represent u(x, t) along one dimen- sion, for example, u(x 0,
ft du(t) dt alto) + f( f(u(t),t, 0) u(t1)= =((t) + f"'
dt
Numerical integration can then be treated as black box, using numerical schemes beyond the forward Euler to achieve higher numerical precision. However, since black box integrators can take an arbitrary number of interme- diate steps, backpropagating through individual steps would require too much memory since the individual steps must be saved. [Cen et al. (2019b) addressed this problem by using adjoint back- propagation, which has constant memory usage. For given loss function on the terminal state t of the hidden state )). the adjoint a(t) is governed by another ODE:
da(t) aat) dt
Of (f(u(t),t,))
0C
a(t1) au((t))
Du
that dictates the gradient with respect to the parameters. The loss L(u(t))) can then be calculated by solving another ODE (the adjoint) rather than backpropagating through the calculations involved in the numerical integration.
However, assuming that the hidden state is governed by ODE imposes limitation the expressiveness of the mapping. For example. Dupont et al. (2019) describes a notable limitation of NODES in the inability represent dynamical systems with intersecting trajectories. In response to such limitations, many works have tried to increase the expressiveness of the mapping. Dupont al. (2019) proposed solve the intersection trajectories problem by augmenting the vector space, lifting the points into additional dimensions; Zhu et al. (202T) included time delay in the equation to represent dynamical systems of greater complexity; Massaroli et al. (202T) proposed to condition the vector field on the inputs, allowing the integration limits to be conditioned on the input; Massaroli] ct al. (202I) and Norcliffe et al. (2020) additionally proposed and proved a second-order ODE system can efficiently solve the intersecting trajectories problem. We note however that the interpretation of NODE as a continuous form of ResNet is also problematic, owing to the fact that the empirical behavior of the ResNet does not match the theoretical properties (Krishnapriyan et al.) 2022 Ott jet al., 2021). As such, alternative interpretations of the process represented by ODE have been considered. In Zhang (2019), the authors considered augmentation where the augmented
2

---

Published as a conference paper at ICLR 2023
state corresponds to the parameters of the network governing the latent state. Queiruga et al. (2021 describes the latent evolution through a series of basis functions thereby allowing important concepts such BatchNorm to be effectively translated in the continuous setting, achieving state-of-the- the-ttt-rttt art performance on a variety of image classification tasks. Further improvements to performance have
been made by considering different numerical integrators (Matsubara et al. 2021 Zhuang al. 2020 2021).
On a related front, multiple works have attempted expand NODE other types of differential equation beyond ODEs. [Sun et al. (2019) employed dictionary method and expanded NODEs to a PDE case, achieving high accuracies both in approximating PDEs and in classifying real-world image datasets. However, [Sun et al. (2019) suggested that the method is unstable when training with the adjoint method and therefore unable to make use of the benefits that come with training with adjoint. PDE, but did
Zhang et al. (2018) proposed density transform approach based on the
not consider using adjoint-based training. Multiple works have expanded the stochastic differential
equations setting and developed efficient optimization methods for them including (Güler et al. 2011
Jia & Benson 2019 Kidger et al. 202 la Li et al. 2020 Liu et al. 2019 Xu et al. 2022). Kidger et al. (2020); Morrill al. (2021a) b) used ideas from rough path theory and controlled differential equations to propose a NODE architecture as a continuous recurrent neural network framework. Salvi
al. (2022) considered stochastic PDEs for spatio-temporal dynamics prediction. Additionally, Chen et al. (2020) models spatio-temporal point-processes using NODEs, and Rubanova et al. (2019):
De Brouwer et al. (2019) makes predictions on time series data using NODEs.
3 METHOD
We describe the proposed C-NODE method in this section by first providing a brief introduction to the method of characteristics (MoC) for solving PDEs with illustrative example. We then discuss the types of PDEs we can describe using this method. We finally discuss how we apply the MoC to our C-NODE framework.
3.1 METHOD OF CHARACTERISTICS
The MoC provides a procedure for transforming certain PDEs into ODEs along paths known as characteristics. In the most general sense, the method applies to general hyperbolic differential equations. We will introduce MoC using a canonical example involving the inviscid Burgers equation, and defer to Griffiths al. (2015] Chapter 9) for a more complete introduction to the topic. Let u(x t) R x R+ + R satisfy the following inviscid Burgers equation
Ou at
+ u Ou
0,
Or
where we drop the dependence on and for ease of notation. We are interested in the solution of over some bounded domain Ω C R x R... Consider parametric forms for the spatial component a(s): (0. 27 + R and temporal components t(s): [0, T + R. over the fictitious variable E [0, 7]. Intuitively, this allows us to solve an equation on curves x, t that are parameterized by a variable s which we denote (s()),((s) as the characteristic. Expanding and writing as the total derivative, we get
d ds m<<(s))((s))
Ou dr du dt at ds'
Or ds
Recall the original PDE in equation and substituting the proper terms into equation 2 for /z/ds
u, di/ds = 1, du/ds = O, we then recover equation Note that we now have a system of 3 ODEs which we can solve to obtain the characteristics as =(s) = us + Tig and ((s) = s + lo.
which are functions of initial conditions zo 4o- Finally, by solving over a grid of initial conditions
{a6 ) E A2, we can obtain the solution of the PDE over Ω. Putting all together, we have a new ODE that is written as
d u(x(ss)
+ u Ou
ds
=
Ot
dr
= 0,
3

---

Published as a conference paper at ICLR 2023
where integrate over through
(()) = [" ds T u(us + zo s)ds,
{f u(((s), (())dss
ds
using the adjoint method with boundary conditions fn- This contrasts the usual direct integration over the variable t that is done in NODE; we now jointly couple the integration through the character- istics. An example of solving this equation over multiple initial conditions is given in Figure | with the contrast to standard NODE integration.
Hyperbolic PDEs in Machine Learning To motivate using MoC for machine learning problems such as classification or density estimation, we again note that MoC. most generally applies to hyperbolic PDEs. These PDEs roughly descrihe the propagation of physical quantities through time. Such equations may be appropriate for deep learning tasks due to their ability to transport data into different regions of the state space. For instance, of transporting high-dimensional data points that are not linearly separable to spaces where they are linearly separable. Similarly, in generative modeling, we transport a base distribution data distribution.
classification task, we consider the problem
3.2 NEURAL REPRESENTATION OF CHARACTERISTICS
In the proposed method, we learn the components involved the MoC, namely the characteristics and the function coefficients. We now generalize the example given in 3.1 which involved two variables, a k-dimensional system. Specifically, consider the following nonhomogeneous boundary value problem (BVP):
du = c(z1, .... Tk. u), tu + E (11 Tk: dz; 0t (()(())0)) = un
x, t R (0) on X E R
Here, u1 & x R dependent values of u and r's. This problem is well-defined and has solution as long as du is continuous (Evans 2010).
Ih" a multivariate map, ag Rk+n - RR and c: Re+n - R are functions
Σx= a. 8x;
MoC is generally used in scalar context, but the correspondence the vector case is relatively straightforward. A proof of this can be found in Appendix C.I To begin, we decompose the PDE equation into the following system of ODEs
dr:
a (x1, o. Tk, u),
ds
dt
11
(5)
ds
du
Su dx; c(x1, Fk. u). Σ
k
(6)
ds Oz, ds
We represent this ODE system by parameterizing du./ ds and 0u/Ox; with neural networks. Conse- quently, du/ds ds is evolving according equation
Following this expansion, we arrive at
u(x(T).T) = o) du (x, u) ds
(7)
ds
dx
m(x(0),0) + [Jau] (x. (x,
(x, O2)
ds
where remove u's dependency x(s) and x's dependency for simplicity of notation. In equation]] 7 the functions J u and dx/ ds are represented as neural networks with inputs x, u and parameters O2.
4

---

Published as a conference paper at ICLR 2023
Conditioning data Previous works primarily modeled the task of classifying set of data points with fixed differential equation, neglecting possible structural variations lying in the data. Here, we condition C-NODE each data point, resulting in solving PDE with different initial condition and the hyperholic PDE interpretation of the latent variables. Consider the term given by the integrand in equation[7 The neural network representing the characteristic dx fx/dd ds is conditioned the input data E _. Define mapping g()) R + IR" and we have dæ; ds
(I............ Th. u.g(z))
By introducing g(z) in equation 8 the equation describing the characteristies changes depending on the current data point. This leads the classification task being modeled with family rather than one single differential equation and allows the C-NODE system to model dynamical systems with intersecting trajectories. This property becomes helpful in Proposition 4.1 in proving that C-NODE can represent intersecting dynamics.
3.3 TRAINING C-NODES
Having introduced the main components of C-NODEs, we can now integrate them into unified algorithm. To motivate this section, and to be consistent with part of the empirical evaluation, we will consider classification tasks with data {(z) (z)
Zj R" Yy E Zt. For instance, Zj may be an image, and yj is its class label. In the approach pursue here, the image zj is first passed through a feature extractor function O1) R + II" with parameters 1 The output of g is the feature (i) g(zz: 1) that provides the boundary condition for the PDE on uli), We integrate along
u
different characteristic curves indexed by E [0, with boundary condition (x(0) 0) uo and compute the terminal values given by equation 7
where we mentioned in Section 3.2.
1 (x(T). = ( + ( (xx u
dx
(xx, u (t) u1
(), O22 ds
(9)
ds Finally, (T),, is passed through another neural network, (u (x(7)); O3) with input u (x(T). and parameters O3 whose output are the probabilities of each class labels for image z.. The entire learning process now reduced to finding optimal weights (θ,, 02, Q3) which can be achieved by minimizing the loss
L Σ (x()
where L (.) is the corresponding loss function (e.g. cross entropy in classification). In Algorithm 2 illustrate the implementation procedure with the forward Euler method for simplicity for the framework but note any ODE solver can be used.
Combining MoC with Existing NODE Modifications As mentioned in Section the proposed C-NODES method can be used as an extension to existing NODE frameworks. In all NODE modifications, the underlying expression of f" f(t, O)dt remains the same. Modifying this expression to f Jxu(x, u; 9)dx//s(xx u, up; O)ds results in the proposed C-NODE architecture, with the size of x being hyperparameter.
PROPERTIES OF C-NODES
C-NODE has number of theoretical properties that contribute to its expressiveness. We provide some theoretical results on these properties in the proceeding sections. We also define continuous normalizing flows (CNFs) with C-NODEs, extending the CNFs originally defined with NODEs.
Intersecting trajectories As mentioned in Dupont et al. 2211 one limitation NODE is that the mappings cannot represent intersecting dynamics. We prove by construction that conditioning initial conditions allows C-NODEs to represent some dynamical systems with intersecting trajectories in the following proposition:
5

---

Published as a conference paper at ICLR 2023
Proposition 4.1. The C-NODE represent a dynamical system u(s), du/ds g(s, u) R x R + R, where when u(0) then u(1) u(0) + ((s, u)ds 0; and when u(0) then u(1) (0) So (s, u)ds 1.
Proof. See Appendix
Density estimation with C-NODEs C-NODES can also be used to define a continuous density flow that models the density of variable over space subject to the variable satisfying PDE, extending the continuous normalizing flows defined with NODEs. For NODEs, if u(t) E IR" follows the ODE
du(t)/dt = (u(()), where [(u(()) E RR" then log likelihood from Chen al. 2019b. Appendix A) is given by:
i Ologp(u(t)
df du(())
(10)
at
Similar to the change of log probability of NODEs, as in equation proposition for C-NODEs:
10
we provide the following
Proposition 4.2. Let u(s) be a finite continuous random variable with probability density function du and de, are uniformly Lipschitz
p(u(s)) and let u(s) satisfy da
Ou dæ, Assuming 8x; de öz; ds continuous u and continuous in s, then the evolution of the log probability of u follows:
du(s) Eizn
Ologp(u(s)) -tr (
0
k
Ou dæ,
Ös
Ou
Ora ds
1=1
Proof. See Appendix C.3
CNFs are continuous and invertible one-to-one mappings onto themselves, i... homeomorphisms. Zhang al. (2020) proved that vanilla NODEs are not universal estimators of homeomorphisms, but augmented neural ODEs (ANODEs) are universal estimators of homeomorphisms. We demonstrate that C-NODEs are pointwise estimators of homeomorphisms, which we formalize in the following proposition:
Y, T C IRP, initial condition up, and time T > 0, there exists a flow u(s, un) E IR" following du Ou da + du dt such that u(T, un) h(wo).
Proposition 4.3, Given any homeomorphism h Y +
Proof. See Appendix C.4
5 EXPERIMENTS
We present experiments image classification tasks, time series prediction tasks, image generation tasks on benchmark datasets, and a synthetic PDE regression task.
5.1 CLASSIFICATION EXPERIMENTS WITH IMAGE DATASETS
We first conduct experiments for classification tasks on high-dimensional image datasets, including MNIST, CIFAR-10 10, and SVHN. We provide results for C-NODE and also combine the frame- work with existing methods, including ANODEs (Dupont et al. 2019). Input Layer NODEs (IL- NODEs) (Massaroli et al.) 2021), and 2nd-Order NODES (Massaroli et al. (021). For all classification experiments, we set the encoder of input images for conditioning to be identity, i.e., g(z) = z, making the input into C-NODE the original image. This way we focus exclusively on the performance of C-NODE.
The results for the experiments using the adjoint method are reported in Table I We investigate the performances of the models on classification accuracy and the number of function evaluations (NFE) taken in the adaptive numerical integration. NFE is indicator of the model's computational complexity and can be interpreted as the network depth for the continuous NODE (Chen et al. 2019b). Using similar number of parameters, combining C-NODEs with different models consistently results

---

Published as a conference paper at ICLR 2023
in higher accuracy and mostly uses smaller numbers of NFEs, indicating better parameter efficiency. An ablation study on C-NODEs' and NODEs' parameters can be found in Appendix E.2 While the average number of function evaluations tends to be lower for C-NODE, we additionally note that, compared to ANODE, training C-NODE with the adjoint method can sometimes have decreased stability. We define instability having a NFE 1000. To get rough idea of the differences in stability, when training NODE, C-NODE, and ANODE for image classification on the SVHN dataset for forty instances, NODE appeared unstable six times, C-NODE was unstable three times, and ANODE was never unstable. We note that this was only apparent in the SVHN experiment and when considering C-NODE by itself; the average NFE decreases when adding C-NODE A-NODE and this was never experienced in the ANODE+C-NODE experiments.
Dataset
Method
Accuracy
NFE 4
Param.[K]
NODE
75.28 0.836%
131
115.444
C-NODE
82.19 + 0.478%
124
113.851
ANODE
89.8 + 0.952%
167
112.234
ANODE+C-NODE
92.23 1 0.176%
146
112.276
SVHN
2nd-Ord
88.22 1.11%
161
112.801
Ord-OOD+C
92.37 1 0.118%
135
112.843
IL-NODE
89.69 0.369%
195
113.368
IL-NODE+C-NODE
93.31 0.088%
113.752
NODE
56.30 0.742%
115.444
C-NODE
64.28 0.243%
113.851
ANODE
70.99 0.483%
112.234
ANODE+C-NODE
71.36 0.220%
112.276
CIFAR-10
2nd-Ord
70.84 0.360%
189
112.801
2nd-OOd+C-NODE
73.68 0. 0.15%%
131
112.843
IL-NODE
72.55 0. 0.238%
134
113,368
IL-NODE+C-NODE
73.78 0. 0.155%
113.752
NODE
96.90 + 0.154%
72
85.468
C-NODE
97.56 + 0.431%
72
83.041
ANODE
99.12 + 0.021%
68
89.408
ANODE+C-NODE
99.20 + 0.002%
60
88.321
MNIST
2nd-Ord
99.35 + 0.002%
52
89.552
2nd Ord-------N
99.38 + 0 0.037%
61
88.465
IL-NODE
99.33 1 0.039%
53
89.597
IL-NODE+C-NODE
99.33 + 0.001%
60
88.51
Table 1: Mean test results over runs of different NODE models over SVHN 10, and MNIST. Accuracy and NFE convergence reported. Applying C-NODE always increases models' accuracy and usually reduces models' NEE as well as the standard error.
5.2 TIME SERIES PREDICTION WITH C-NODES
We test C-NODEs, NODEs, and augmented versions of C-NODEs and NODES on the time series prediction problem using 100 sequences of each of the following datasets: the PhysioNet dataset, containing measurements from the first 48 hours of patients' admissions to ICU (Goldberger et al.) 2000 (June 13): the Hopper dataset, containing physical simulation results generated with the Hopper model from the Deepmind Control Suite (Rubanova et al.) 2019: Tassa et al.) 2018);; the Human Activity dataset in the UCI dataset, containing time series data on twelve features from five individuals walking, sitting, and lying (Dua & Grafi 2017); and the synthetic ODE dataset, generated with first-order ODE whose initial condition follows a Gaussian distribution, and Gaussian noises are added to the observations (Rubanova al. 2019). The PhysioNet and Human Activity datasets both contain stochastic processes while the Hopper dataset and the ODE dataset are deterministic. Since NODE, ANODE, and C-NODE model deterministic processes, the Hopper dataset and the ODE
7

---

Published as a conference paper at ICLR 2023
Dataset
Method
MSE 4
Param.[K]
NODE
- 10"2
76.663
C-NODE
(1.07 + x 10
76.205
PHYSIONET
ANODE
(1.15 10-
86.
C-NODE ANODE
(1.04 10
85.705
NODE
(1.998 10-1
51.042
C-NODE
(1.797 10
50.964
HUMAN ACTIVITY
ANODE
(1.140 0.048) 10-1
69.948
C-NODE ANODE
09.22 10
69.397
NODE
(5.87 + 0.23) x 10-2
44.776
C-NODE
(5.68 0.11) ) 10-2
44.218
HOPPER
ANODE
(5.68 + 0.10) x 10-
51.976
C-NODE ANODE
(4.96 + 0.16) 10-2
51.506
NODE
(4.73 1 0.17)
0.17) x 10-2
43.231
C-NODE
(4.17 + 1.34)
x 10 10-2
42.935
ODE
ANODE
(1.60 1 0.53) x 10-2
43.767
C-NODE + ANODE
(4.01 + 1.23) 10 10-2
43.556
Table Mean test results over 4 runs different NODE models over PhysioNet, Human Activity, Hopper, and ODE. Test mean square errors are reported. Applying C-NODE always reduces the test error, and mostly reduces' the standard error. Training dynamics are shown in Figures Appendix A.2
in
dataset readily fit the modeling frameworks. On the other hand, the PhysioNet and Human Activity datasets require augmenting the dynamics of the NODE models with stochasticity to model the arrival of events.
We follow the experimental setup for interpolation tasks Rubanova (2019). where we define an autoregressive model with the encoder being an ODE-RNN model and the decoder being latent differential equatior The main purpose of this experiment to compare the ODE-RNN when using C-NODE versus NODE. ODE-RNN is a standard method for including ODE modeling in time series tasks as described Chen et al. (2019b) and Rubanova et al. (2019). We consider interpolation tasks by first encoding the time series { ( of length N and computing the approximate posterior g(z0 (æ., 4) === = (Hzo ozo) done in Rubanova et al. (2019). Then, (220), (20) are computed as ODDERNN,(122) (o)), where is function that encodes the model whose states obey an ODE between observations and are updated according to new observations as described in Rubanova al.] (2019). To predict the state at an observation time t,, we sample initial states 2o from the posterior which are then decoded using another neural network. We finally average generated observations at each observation time to compute the test data errors. The results are presented in Table and suggest that C-NODE based models use slightly fewer parameters while achieving lower eтror rates than NODE models. C-NODE models the latent dynamics as first order PDE, which is natural extension of the ODE model that NODE uses.
terminal hidden states into mean and variance of the latent variable zg. ODE-RNN(-)
5.3 CONTINUOUS NORMALIZING FLOW WITH C-NODES
We compare the performance of CNFs defined with NODEs to flows defined with C-NODEs on MNIST and CIFAR-10, 10. We use the Hutchinson trace estimator to calculate the trace and use multi- scale convolutional architectures to model the density transformation as done in (Dinh et al. 2017] Grathwohl et al. ) Differential equations are solved using the Runge-Kun Kutta method of order of the Dormand- SwwtaanddPriin---e Shane---haminee solver and trained with the adjoint method. Although the Euler
2This is hased the code of Rubanova et al.
2019 provided
https::/github.
Yul wawaaaaaaaaaaaaaaaaaaaaaaa aatiaaaoooalaattt__ooo
This is based on the code that the authors of (rrttwohl et al. 2019 provided in https://github
com/ hamm/rqiian chrrgifhheFFFoor
8

---

Published as a conference paper at ICLR 2023
MNIST
CIFAR-10
Model
B/D 4
Param. +
NFE
B/D 4
Param.
NFE 4
Real NVP (Dinh al.
2017
1.05
N/A
3.49
N/A
Glow (Kingma & Dhariwal
2018
1.06
N/A
3.35
44.0M
RQ-NSF (Durkan al. 2019)
3.38
25.2M
Res. Flow (Chen al.
2019a
0.97
16.6M
3.28
25.2M
CP-Flow Huang al.
2021
1.02
2.9M
3.40
1.9M
NODE
1.00
335.1K
1350
3.49
410.1K
1847
C-NODE
0.95
338.0K
1323
3.44
406.0K
1538
Table 3: Experimental results generation tasks, with NODE, C-NODE, and other models. B/D indicates Bits/dim. Using similar number parameters, C-NODE outperforms NODE on all three datasets, and has significantly lower NFE when training for CIFAR-10.
forward method faster, experimental results show that its fixed step size often leads to negative Bits/Dim, indicating the importance of adaptive solvers. As shown in table figure using similar number of parameters, experimental results show that CNFs defined with C-NODEs perform better than CNFs defined with NODEs in terms of Bits/Dim, well having lower variance, and using a lower NFE on both MNIST and CIFAR-10.
5.4 PDE MODELING WITH C-NODES
We consider regression example for a hyperbolic PDE with a known analytical solution. Since NODEs assume that the latent state is only dependent on a scalar (namely time), they cannot model dependencies that vary over multiple spatial variables required by most PDEs. We modify the assumptions used in the classification and density estimation experiments where the boundary conditions were constant as in equation
We approximate the following BVP:
u
Dz
u, = 2t,
(11)
u(22
1 x 2
which has an analytical solution given by u(x, t) 2x exp(t)
2exp(!)+I' We generate a training dataset by
randomly sampling 200 points (x. f), E t j0, 1). as well as values u(x, t) at those points. We test C-NODE and NODE on 200 points randomly sampled as (x, t) E [1 [o. 1]. For this experiment, C-NODE uses 809 parameters while NODE uses 1185 parameters. We quantify the differences in the representation capabilities by examining how well each method can represent the PDE. C-NODE deviates 8.05% from the test set, while NODE deviates 30.52%. Further experimental details can found in Appendix(A.4. A.4.I
DISCUSSION
We describe an approach for extending NODEs to the case of PDEs by solving series of ODEs along the characteristics of a PDE. The approach applies to any black-box ODE solver and can combine with existing HODE-based frameworks. We empirically showcase its efficacy on classification tasks while demonstrating its success in improving convergence using Euler forward method without the adjoint method. C-NODEs also consistently achieve lower testing MSEs over different time series prediction datasets, while having lower standard errors. Additionally, C-NODEs empirically achieve better performances on density estimation tasks, while being more efficient with the number of parameters and using lower NFEs. C-NODE's efficiency over physical modeling is also highlighted with additional experiments. Discussion on limitations can be found in Appendix B
ACKNOWLEDGMENTS
This work was supported in part by the Office of Naval Research (ONR) under grant number N00114-1-1-2590. AH was supported by NSF-GRFP.

---

Published as a conference paper at ICLR 2023
REFERENCES
Ricky T. Q. Chen, Jens Behrmann, David Kristjanson Duvenaud, and JömnHenrik Jacobsen. Residual flows for invertible generative modeling. arXiv preprint arXiv:1906.c 02735, 2019a.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary
differential equations. arXiv preprint arXiv:1806.07366 2019b.
Ricky TQ Chen, Brandon Amos, and Maximilian Nickel. Neural spatio-temport point processes. arXiv preprint arXiv:2011.0458) 2020.
Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. Gru-ode- Gruu-odebbyes: Continuous modeling of sporadically time series. Advances in neural information processing systems, 32, 2019.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real nvp. arXiv preprint
2017.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http://archiv/ lcs. wcw.udu/ml
Emilien Dupont, Amaud Doucet, and Yee Whye Teh. Augmented neural ODEs. arXiv preprint 190.000904.0681., 2019.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Advances in neural information processing systems, 32, 2019.
Lawrence Evans. Partial Differential Equations. American Mathematical Society, 2010.
A. L. Goldberger, L. A. N. Amaral, L. Glass, J. M. Hausdorff, P. Ch. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, and H. E. Stanley. PhysioBank, PhysioToolkit, and PhysioNet: Compo- nents of a new research resource for complex physiologie signals. Circulation, 101(23)e215-<222,, 2000 (June 13). Circulation Electronic Pages: http://eirc. http://chh.....oaooaaaaoooeen/// PMID:1085218: 8; doi: 1116://1...R.111.23...
Will Grathwohl, Ricky T. Q Chen, Jesse Bettencourt, Ilya Sutskever, and David Kristjanson Duve- naud. Ffjord: Free-form continuous dynamics for scalable reversible generative models. arXiv preprint ar arXiv:1810.01367 2019.
David F Griffiths, John W Dold, and David J Silvester. Essential partial differential equations. Springer, 2015.
Batuhan Güler, Alexis Laignelet, and Panos Parpas. Towards robust and stable deep learning algo- rithms for forward backward stochastic differential equations. arXiv preprint arXiv: 1910. 16223, 2019.
Kaiming He. Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. arXiv preprint arXiv:1512.03385 2015.
Ralph Howard. The gronwall inequality, 1998. URL bttp://peoole..aahhscc..o
howardd/kotes/groowall...
Chin-Wei Huang, Ricky T. Q. Chen, Christos Tsirigotis, and Aaron C. Courville. Convex potential flows: Universal probability distributions with optimal transport and convex optimization. arXiv preprint 20X0::20/2.05942,, 2021.
Junteng Jia and Austin R Benson. Neural jump stochastic differential equations. Advances in Neural Information Processing Systems, 32, 2019.
Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equations for irregular time series. arXiv preprint arXiv:2005.0891 2020.
Patrick Kidger, James Foster, Xuechen Li, and Terry Lyons. Efficient and accurate gradients for neural sdes. arXiv preprint arXiv:2105. 13493, 2021a.
10

---

Published as a conference paper at ICLR 2023
Patrick Kidger, James Foster, Xuechen Li, and Tenry J Lyons. Neural sdes dimentionll gans. In International Conference on Machine Learning, pp. 5453-5463. PMLR, 2021b.
Durk P Kingma and Prafulla Dhariwal. Glow: Generative flow with invertible 1x1 convolutions. Advances in neural information processing systems, 31, 2018.
Aditi S. Krishnapriyan, A. F. Queiruga, N. Benjamin Erichson, and Michael W. Mahoney. Learning continuous models for continuous physics. ArXiv, abs/2202.08494, 2022.
Xuechen Li, Ting-Kam Leonard Wong, Ricky T. Q. Chen, and David Duvenaud. Scalable gradients for stochastic differential equations. arXiv preprint vrXiv:2001.013 328, 2020.
Xuanqing Liu, Tesi Xiao, Si Si, Qin Cao, Sanjiv Kumar, and Cho-Jui Hsieh. Neural sde: Stabilizing 2019.
neural ode networks with stochastic noise. arXiv preprint
Stefano Massaroli, Michael Poli, Jinkyoo Park, Atsushi Yamashita, and Hajime Asama. Dissecting neural odes. arXiv preprint arXiv:2002.0807/) 2021.
Takashi Matsubara, Yuto Miyatake, and Takaharu Yaguchi. Symplectic adjoint method for exact gradient of neural ode with minimal memory. Advances in Neural Information Processing Systems, 34:20772-20784, 2021.
James Morrill, Patrick Kidger, Lingyi Yang, and Terry Lyons. Neural controlled differential equations for online prediction tasks, 2021a.
James Morrill, Cristopher Salvi, Patrick Kidger, James Foster, and Terry Lyons. Neural rough differential equations for long time series. arXiv preprint arXiv:2009.0829) 2021b.
Alexander Norcliffe, Cristian Bodnar, Ben Day. Nikola Simidjievski, and Pietro Liö. On second order behaviour in augmented neural odes. arXiv preprint arXiv:2006.0722 2020.
Katharina Ott, Prateck Katiyar, Philipp Hennig, and Michael Tiemann. Resnet after all: Neural [ode]s and their numerical solution. In International Conference on Learning Representations. 2021. URL BRR.itttps///peereere/
net. fww.m/orrmmmmmo d==.11d55222222222552 wtt.:///eeooouuooooo.....eexooo.
Alejandro Francisco Queiruga, N. Benjamin Erichson, Liam Hodgkinson, and Michael W. Mahoney. Stateful ODE-nets using basis function expansions. In A. Beygelzimer, Y. Dauphin, P. Liang, and Wortman Vaughan (eds.), Advances in Neural Information Processing Systems, 2021. URL htttps://openreve ew. net htt./orrtt//rer (errummid=-yyyp (iddkvitttt rmA
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows. In ICML, 2015.
Yulia Rubanova, Ricky TQ Chen, and David K Duvenaud. Latent ordinary differential equations for ireeularly-samplee time series. Advances in neural information processing systems, 32, 2019.
I Lars Ruthotto and Eldad Haber. Deep neural networks motivated by partial differential equations. Journal of Mathematical Imaging and Vision, 62(3):352-364, 2020.
Cristopher Salvi, Maud Lemercier, and Andris Gerasimovics. Neural stochastic partial differential equations: Resolution-invariant learning of continuous spatiotemporal dynamics. arXiv preprint arXiv:2110.10249, 2022.
Robert Strichartz. The way of analysis. Jones and bartlett mathematics, 2000.
Yifan Sun, Linan Zhang, and Hayden Schaeffer. Neupde: Neural network based ordinary and partial differential equations for modeling time-dependent data. arXiv preprint 1XXi:19008.03999 2019.
Yuval Tassa, Yotam Doron, Alistair Muldal, Tom Erez, Yazhe Li. Diego de Las Casas, David Budden, Abbas Abdolmaleki, Josh Merel, Andrew Lefrancq. Timothy P. Lillicrap, and Martin A. Riedmiller. Deepmind control suite. ArXiv, abs/180 100 100600, 2018.

---

Published as a conference paper at ICLR 2023
E. Weinan. A proposal machine learning via dynamical systems. Communications in Mathematics and Statistics, 5, 2017.
Winnie Xu, Ricky TQ Chen, Xuechen Li, and David Duvenaud. Infinitely deep bayesian neural net- works with stochastic differential equations. In International Conference on Artificial Intelligence and Statistics, pp. 721-738. PMLR, 2022.
Han Zhang, Xi Gao, Jacob Unterman, and Tom Arodz. Approximation capabilities of neural odes and invertible residual networks. In ICML, 2020.
Linfeng Zhang, Weinan E, and Lei Wang- Monge-ampère flow for generative modeling. arXiv preprint arXiv: 10188, 2018.
Tianjun Zhang. Zhewei Yao, Amir Gholami, Kurt Keutzer, Joseph Gonzalez, George Biros, and Michael Mahoney. Anodev2: A coupled neural ode evolution framework. arXiv preprint 1906.099000 2019.
Qunxi Zhu, Yao Guo, and Wei Lin. Neural delay differential equations. In International Confer- ence on Learning Representations, 2021. URL
Q1 jmmmz 2M2
Juntang Zhuang, Nicha Dvomek, Xiaoxiao Li, Sekhar Tatikonda, Xenophon Papademetris, and James Duncan. Adaptive checkpoint adjoint method for gradient estimation in neural ode. In International Conference Machine Learning, pp. 11639-11649, PMLR, 2020.
Juntang Zhuang, Nicha C Dvornek, sekhar tatikonda, and James s Duncan. (MALI): memory efficient and reverse accurate integrator for neural (ode) Is. In International Conference Learning Representations, 2021.
A EXPERIMENTAL DETAILS
Implementation details of this paper can be found at NeuralPDE. SibralPDE.git
A.1 EXPERIMENTAL DETAILS OF CLASSIFICATION TASKS
We report the average performance over five independent training processes, and the models are trained for 100 epochs for all three datasets. We also report the training dynamics of C-NODE and NODE using the adjoint sensitivity method and the euler backpropagation, as shown in Figure 2
As shown in Figure using the Euler solver, it appears that C-NODEs converge faster than the vanilla NODEs (usually in one epoch) while generally having a more stable training process with smaller variance. Additionally, on experiments with MNIST, C-NODEs converge in only one epoch, while NODEs converge in roughly 15 epochs. This provides additional empirical evidence on the benefits of training using the characteristics.
The input for 2nd-Ord. NODE, and C-NODE are the original images. In the IL-NODE, we transform the input to a latent space before the integration by the integral; that is, we raise the exhx dimensional input image into the k(c+p)xkxw dimensional latent feature spaceF)
For SVHN and CIFAR-10, we assume x E R3 i.e., the Jacobian Jxu (Qu/Ox,, du/ 00.00000//000 0u/Ora). We model each partial derivative 0u/Or; with a separate convo- lutional network. The network architecture for the network modeling the partial derivatives is as shown in Table The network architecture for the network modeling dx ds is as shown Table 6 The architecture in Tables 6lae used for both CIFAR-10 10 and SVHN. Note that the network archi- tecture for MNIST differs slightly due the lower dimensionality of MNIST. The hyperparameters used are as shown in T
TThis is based on the code of Massaroli et al. 2021 provide in torchöýn
com/Di fFFqMI./
12

---

Published as a conference paper at ICLR 2023
130
130
20
20
20
20
20
20
42
2.5-
2.0
2.0
20
20
2.0
2.20
2.0
1.5
42
20
20
0.75
C11A4
(a) Adjoint training:
(b) Backprop through Euler training:
Figure 2: Red: NODE, Blue: C-NODE. Training dynamics of different datasets with adjoint in Fig. and with Euler in Fig. 2b averaged over five runs. The first column is the training process of SVHN. the second column of CIFAR-IO, and the third column is of MNIST. By incorporating the C-NODE method, we achieve a more stable training process in both CIFAR-10 10 and SVHN, while achieving higher accuracy. Full-sized figure in supplementary materials.
We decode the result after performing the continuous transformations along characteristics curves, back to the Rexhxw dimensional object space. Combining this with the C-NODE can seen as solving PDE on the latest features of the images rather than the images directly.
Unlike ODEs, we take derivatives with respect different variables in PDEs. For a PDE with variables, this results in the constraint of the balance equations
i, E (1, A),
dz,I;"
This can be satisfied by defining the k-th derivative with neural network, and integrate k - times to get the first order derivatives. Another way of satisfying the balance equation to drop the dependency on the variables, i.e., Vi E (1,2,..., k),
du dr;
When drop the dependency, all higher order derivatives zero, and the balance equations satisfied.
All experiments were performed NVIDIA RTX 3090 GPUs cloud cluster.
Data
MNIST
SVHN
CIFAR10
Model
CNN
CNN
CNN
# of PDE Dimensions
2
3
3
Optimizer
Adam
Adam
Adam
Learning Rate
1.00E-3
1.00E-3
1.00E-3
Weight Decay
5.00E-04
5.00E-04
5.00E-04
Table 4: Training hyperparameters for image classification.
A.2 EXPERIMENTAL RESULTS AND DETAILS TIME SERIES PREDICTIONS
A.2.1 EXPERIMENTAL DETAILS OF TIME SERIES PREDICTIONS ON REAI-WORLDD DATASETS
We provide more detailed explanation of the time-series experiments which are based the ODE-RNN framework described in Rubanova et al. (2019). Our experiments test the effectiveness of NODE, C-NODE, ANODE, and their combinations under the ODE-RNN framework by computing
13

---

Published as a conference paper at ICLR 2023
Operation Layer
Kernel Size
Stride
Padding
Dilation
Size of Output
Input
N/A
N/A
N/A
N/A
x 32 x 32
Convolution Layer
3
1
0
1
60 X 32 x 32
ReLU
N/A
N/A
N/A
N/A
60 x x
Convolution Layer
3
1
0
1
60 x 32 X 32
ReLU
N/A
N/A
N/A
N/A
60 x 32 x 32
Convolution Layer
3
1
1
x 32 x 32
Table Network architecture of the network modeling terms the Jacobian.
Operation Layer
Kernel Size
Stride
Padding
Dilation
Size of Output
Input
N/A
N/A
N/A
N/A
4 x 32 32
Convolution Layer
1
I
32 32
ReLU
N/A
N/A
N/A
N/A
8 x 32 x 32
Convolution Layer
1
O
32 32
ReLU
N/A
N/A
N/A
N/A
x 32 x 32
Convolution Layer
1
1
0
1
3 x 32 32
Flatten
N/A
N/A
N/A
N/A
3072 X
Linear
N/A
N/A
N/A
N/A
X 1
ReLU
N/A
N/A
N/A
N/A
x 1
Table Network architecture of the network modeling the characteristic curve dx/ds.
evaluation metrics under the different ODE models. The ODE-RNN framework involves combining the strengths of both a neural ODE and RNN to initially embed the time series history into a latent distribution parameterized by an RNN and then decode the predicted latent distribution into the original data space. The main components of the method are illustrated in Algorithm [ Forecasts computed by integrating the latent space in time according the neural ODE model with initial condition zo distributed according to the parameterization by the RNN. We define the approximate N (M oz.). where P ds (ODEE-NNN, & ((2,, o)) To bypass the RNN requirement of a fixed observation rate, the ODE-RNN model uses states that obey an ODE in between observations and are updated new observations. Given set of time series data {x; fa ).M we embed the time series using the ODE-RNN. We then pass the output of the ODE-RNN model through function e(-) to get the mean N and the variance a of the posterior of 20- To predict the value at timestamp T, sample K initial conditions 2n from the posterior (=,, 4 o)), integrate the latent ODE model until timestamp and use the average the K integrations the result.
posterior of 2g as (20o [I, t4 1.V
For all experiments, we follow the experimental setup as described in inthtss:: http:/ggitthba com/ wuliaRubanova/aatee
We experiment with the ODE-RNN framework using NODE, C-NODE, ANODE, and their combinations as the ODE backbone of both the ODE-RNN model and the latent ODE. Training NODE follows the original setup, with the dimension of the ODE model in the ODE-RNN being 20, the number of units per layer in each of GRU update networks being 100, the number of units per layer in ODE function being around 100, the number of layers in ODE function in generative and the recognition ODE both being 1.
We use a C-NODE with dimensionality of 8. The network architecture details are given in Tables 7 8 The training hyperparameters are given in Table The number of units per layer in the network describing dx ds 12. For the network describing Ju/Ox, the dimension the ODE model in the ODE-RNN is 20, the number of units per layer in each of GRU update networks is 100, the number of units per layer in the ODE function is tuned to match the number of parameters in the NODE models, the number of layers in ODE function in generative and the recognition ODE is I.
14

---

Published as a conference paper at ICLR 2023
Algorithm 1 Algorithm for ODE-RNN model(Rubanova et al. 2019)
Input: Data points and their timestamps { !=1,...
ho 0
for i in 1, .... do
M ODDSSovee[ee h (4)))44)
h; RNNCell(h) x44)
end for O2 OutputNN(h:.) for all 1, N Return: {oi N: hN
0.022
C-NODE
C-NODE
0.020
NODE
0.020
NODE
0.018
Test MSe
Test MSE
0.016
0.014
0.014
0.012
0.010
0.010
60 Epoch
60
80
100
Epoch
(a) Without Augmentation;
(b) With Augmentation;
Figure 3: Red: NODE, Blue: C-NODE. Training dynamics of CNODE, NODE, and their augmented versions on Physionet dataset. C-NODE achieves lower testing MSE, and has a more stable training dynamics
Operation Layer
Input Features
Output Features
Linear Layer
20
70
Tanh
N/A
N/A
Linear Layer
70
70
Tanh
N/A
N/A
Linear Layer
70
160
Table 7: Network Structure of 0u/Ox, when using ODE dataset
0.7
C.NDDE NODE
0.6
C-NODE
NODE
0.6
0.5
0.5
M5R
0.4
0.3
0.2
0.1
20
40
60 Epoch
100
40
60 Epoch
80
100
(a) Without Augmentation;
(b) With Augmentation;
Figure 4: Red: NODE. Blue: C-NODE, Training dynamics of CNODE, NODE, and their augmented versions on Human Activity dataset. C-NODE achieves lower testing MSE, and has a more stable training dynamics
15

---

Published as a conference paper at ICLR 2023
0.09
C-VODE NODE
C-NODE
NODE
0.08
0.08
Test MSE
0.07
0.07
0.06
0.06
0.05
20
40
60
80
100
20
40
60
80
100
Epoch
Epoch
Without Augmentation;
(b) With Augmentation;
Figure 5: Red: NODE, Blue: C-NODE. Training dynamics of CNODE, NODE, and their augmented versions Hopper dataset. C-NODE achieves lower testing MSE, and has stable training dynamics
0.525
C-NODE
C-NODE
NODE
0.500
0.50
Test MSE
0.475
Test MSE
0.450
0.45
0.425
0.40
0.400
0.375
Epoch
Epoch
(a) Without Augmentation;
(b) With Augmentation
Figure Red: NODE, Blue: C-NODE. Training dynamics CNODE, NODE, and their augmented
versions toy ODE dataset. C-NODE achieves significantly lower testing MSE
Operation Layer
Input Features
Output Features
Linear Layer
2
Tanh
N/A
N/A
Linear Layer
12
Tanh
N/A
N/A
Linear Layer
12
Table 8: Network Structure of dx/ds when using ODE dataset
Data
Physionet
Activity
Hopper
Periodic
# of PDE Dimensions
8
8
8
7
Optimizer
Adamax
Adamax
Adamax
Adamax
Learning Rate
1.0E-2
1.0E-2
1.0E-2
1.0E-2
Weight Decay
0.0
0.0
0.0
0.0
Table 9; Training hyperparameters for time series analysis
16

---

Published as a conference paper at ICLR 2023
A.2.2 EXPERIMENT RESULTS OF TIME SERIES PREDICTIONS ON SYNTHETIC DATASET
We C-NODES, ANODES, and NODES synthetic time series prediction problem. We define a function by u(x, t) =
2exp(f)+1 and we sample = u(x, t) + 0,1e,, where E ~ >V 1) over
2z exp(t)
E [1.2], t E [0, 1] to generate the training dataset. We test the performance on & E [n, n + 1] with
n {0, 1, 5}.
Time
[0,1]
[1,2]
[2.3]
[3,4]
[4,5]
[5,6]
NODE
0.0322
0.1764
0.4681
0.8093
1.1911
1.6202
ANODE
0.0428
0.0629
0.1248
0.2778
0.5360
0.9252
C-NODE
0.0270
0.0365
0.0582
0.1474
0.3300
0.6054
Table 10: Time series prediction results for NODE, ANODE, and C-NODE different time intervals. Errors are testing mean squared errors. Across all time intervals, C-NODE outperforms NODE and ANODE.
We also test C-NODEs, NODEs, and ANODEs on time series prediction with different levels of noise. Specifically, using the function above, form training and testing dataset with ~ (0, m E {0, 1, .j). We test the performance on the time period E [0, 1].
Noise Level
0
2
3
4
5
NODE
0.0326
0.1784
0.7886
1.9685
3.7530
6.5533
ANODE
0.04
0.1984
0.6035
1.0574
1.4850
2.0593
C-NODE
0.0267
0.1011
0.3294
0.7148
1.2856
2.0834
Table 11:
Time series prediction results for NODE, ANODE, and C-NODE different noise levels. Errors are testing mean squared errors.
A.2.3 EXPERIMENTAL DETAILS OF TIME SERIES PREDICTIONS ON SYNTHETIC DATASET
We consider the task of predicting u(x, t) the initial condition of (1,0).
at different times i, over I E [1, We specify 2.e4TT
2zz.e
We use a dimensional C-NODE network. The result is calculated with
=(z,t)=u(1,, = u(1,0) + / ==1
Ou dz, ds.
Dz, ds
NODE is calculated with
((s,t)==(1)) + f
Su dt.
ōt
The experiment results are given in table II
In our experiments, C-NODEs use 1221 parameters, ANODES use 1270 parameters, NODEs use 1290 parameters.
All experiments were performed on NVIDIA RTX 3080 GPUs on local machine.
A.3 EXPERIMENTAL DETAILS OF CONTINUOUS NORMALIZING FLOWS
We report the average performance over four independent training processes. As shown in Figure 8 compared to NODE, using a C-NODE structure improves the stability of training, as well as having better performance. Specifically, the standard errors for C-NODEs MNIST, SVHN, and CIFAR-10 are 0.37%, 0.51 %, and 0.24% respectively, and for NODEs the standard errors on MNIST, SVHN and CIFAR-10 1.07%, 07%, 0.32%, and 0.22% respectively. The training dynamics of CNFs MNIST
17

---

Published as a conference paper at ICLR 2023
2.2
NODE
2.0
1.8
1.6
1.2
1.0
40
Epoch
Figure 7: Red: NODE. Blue: C-NODE. Training dynamics of CNFs on MNIST dataset with adjoint
method. We present Bits/dim of the first 50 training epochs.
Bits/dim
2
C-NODE NODE
1
Epoch
Bits/din /dim
3.0
2.5
Epoch
6
Bits/dim
4
0
20
40
60
80
Epoch
Figure 8: The training process averaged over 4 runs of C-NODE and NODE. The first row the results on MNIST, the second row are the results on SVHN, the third row are the results on CIFAR-10.
are reported in Figure 7 while having more stable training dynamics.
C-NODE consistently achieves lower Bits/dim across the training process,
The training hyperparameters are provided in Table 14 The experiments are developed using code adapted from the code that the authors of Grathwohl et al. (2019) provided in https:: github....... com/ rt qicnen/////eee ford
The network structures are provided in Tables
All experiments were performed on NVIDIA RTX 3090 GPUs cloud cluster.
18

---

Published as a conference paper at ICLR 2023
Operation Layer
Input Channels
Output Channels
Kernel Size
Stride
Padding
Convolutional Layer
8
(3,3)
(1,1)
(1,1)
Tanh
N/A
N/A
N/A
N/A
N/A
Convulutional Layer
70
(3,3)
(1,1)
(1,1)
Tanh
N/A
N/A
N/A
N/A
N/A
Convolutional Layer
71
70
(3,3)
(1,1)
(1,1)
Tanh
N/A
N/A
N/A
N/A
N/A
Convolutional Layer
71
8
(3,3)
(1,1)
(1,1)
Tanh
N/A
N/A
N/A
N/A
N/A
Convolutional Layer
3
(3,3)
(1,1)
(1,1)
Table 12: Network structure 0u/Ox; when using CIFAR- 10 dataset
Oper. Layer
Input Chan.
Output Chan.
Kern.
Stride
Pad.
Input Fea.
Output Fea.
Cony. Layer
4
12
(3,3)
(1,1)
(1,1)
N/A
N/A
SiLU
N/A
N/A
N/A
N/A
N/A
N/A
N/A
Conv. Layer
12
12
(3,3)
(1,1)
(1.1)
N/A
N/A
SiLU
N/A
N/A
N/A
N/A
N/A
N/A
N/A
Conv. Layer
12
3
(3,3)
(1,1)
(1,1)
N/A
N/A
Flatten
N/A
N/A
N/A
N/A
N/A
N/A
N/A
Linear Layer
71
8
(3,3)
(1,1)
(1,1)
3072
Table 13: Network structure of dx ds when using CIFAR-10 dataset
Data
MNIST
SVHN
CIFAR-10
# of PDE Dimensions
2
Optimizer
Adam
Adam
Adam
Learning Rate
1.00E-3
1.00E-3
1.00E-3
Weight Decay
0.0
0.0
0.0
Table 14: Training hyperparameters for continuous normalizing flow
19

---

of The result is calculated as Published as a conference paper at ICLR 2023
A.4 EXPERIMENTAL DETAILS OF PDE MODELING
A.4.1 2-DIMENSIONAL BURGER'S S EQUATION
We want to solve the initial value problem
Hhe
Bu
u
08
u,
0)
2x.
15
where the exact solution is u(x,t) t) =
(2e"+1] Our dataset's input are 200 randomly sampled points
2ne
(ze, (), [1,]] [0, 1), and the dataset's outputare the exact solutions at those points.
For the C-NODE architecture, we define four networks: (N, t) for as (t) for the characteristic path (z(s),t(s)) (N)(x) for the initial condition. The result calculated in four steps:
Ou N22 6) for Ru 2
1. Integrate Au So du(z(s)))(s) ds Se
Bu dt
du dr ds NN₂ NN₂ NN
de
BE
N Na[l]ds before.
2. Given z, t, solve equation L + N N * = x for iteratively, with 4++1
I N ((N)))) tt. to initialized be
3. Calculate initial value u(x(0), t(0)) (N)
4. u(z, t) = Au + t(())),(()))
For the NODE architecture, we define one network: NN1 (z, t) for &u u(x, I) fo कe di So dtt.
All experiments were performed on NVIDIA RTX 3080-TI GPUs on a local machine.
A.4.2 100-DIMENSIONAL CONVECTION EQUATION
We additionally experiment on solving a 100-dimensional convection equation given by:
du div (t) (t) t)) 9
x(x, 0) cxp( 2 x ((0)
0.01t
m(t)= sin
where x E R 100 u R 100 R+ R. / R4 100
This equation has an analytical solution u(x, t) = exp( -xp(-32lx -
We generate training dataset with 1000 points, and testing dataset with 1000 points. We uniformly sample x, with each æ; E [0, 0.5) and / j0, 10). We calculate the output using the analytical solution.
We define three networks for C-NODE. The first network models dx/ds R01 + R 101, the second network models Juu: R001 + 101 the third network I' R00 - models the initial condition u(x, 0). We first integrate (x, () using the first two networks. The output of these networks is then input to the third network to obtain the value of the solution. The total number of parameters for the networks is 11611.
We define two networks for NODE. The first network models du/ds : 101 - 10101 the second network T 101 R models the initial condition u(x, 0). We first integrate (x, using the first network. The output 101 We put the output into the second network, and arrive at the output. The number of parameters used here is 14214.
We test C-NODE and NODE on the dataset with Gaussian noises at different magnitudes. As shown in Table 15 C-NODE performs 14.2% better than NODE when there is noise, 44.4% better
20

---

Published as a conference paper at ICLR 2023
0.50
NODE C-NODE
Computation Time
0.040
NODE C-NODE
0.45
Computation
0.035
0.40
0.35
0.030
0.30
0.025
10
20
30
40
40
(a) Computation Time With Adjoint Method;
(b) Computation Time With Euler Forward Method;
Figure 9: Red: NODE. Blue: C-NODE. Computation time each epochs using the adjoint method and Euler forward integration. C-NODE uses slightly more computation time for Euler forward integration, and the same amount of time when using the adjoint method.
when there's 's 10% noise, 39.18% better when there's 's 20% noise, and 40.1% better when there's 's 30% noise. The architecture of C-NODE and NODE are as given Tables I6 optimized with the AdamW optimizer, with learning rate of x T0 and weight decay of 5 x 10-4
18] The models
We also report computation time of NODE and C-NODE using the adjoint method and the Euler forward integration. As shown in Figure 9 C-NODE uses slightly more computation time for Euler forward integration, and the amount of time when using the adjoint method. Although each call of C-NODE uses slightly more time, NODE uses bigger number of functional evaluations. This results in roughly a similar computation time for NODE and C-NODE when integrated using the adjoint method.
note on the difference between plysics-informed neural networks (PINNs) PINNs (Raissi et al. (2019) is method that solves PDE using regularization approach by minimizing regular- ization term enforcing the PDE. Comparing this to C-NODE, there are a few apparent differences. First, PINNs is very general framework that accommodate almost any type of PDE. This contrasts with C-NODE which only computes solutions to hyperbolic PDEs. On the other hand, PINNs generally scale poorly with dimension due to the difficulty in computing high dimensional derivatives. Besides, PINNs require the exact form of the PDEs being solved, whereas C-NODE does require this information.
Noise Level
10%
20%
30%
NODE
0.0500
0.0644
0.0675
0.0688
C-NODE
0.0438
0.0446
0.0485
0.0491
Table 15: C-NODE performs better than NODE on all noise levels. C-NODE uses 1161 parameters, NODE uses 14215 parameters.
Operation Layer
Input Channels
Output Channels
Linear Layer
101
16
Softplus
N/A
N/A
Linear Layer
16
16
Softplus
N/A
N/A
Linear Layer
16
101
Table 16: Network structure of J
21

---

Published as a conference paper at ICLR 2023
Operation Layer
Input Channels
Output Channels
Linear Layer
102
16
Softplus
N/A
N/A
Linear Layer
16
16
Softplus
N/A
N/A
Linear Layer
16
101
Table 17: Network structure of dx/ds ds
Operation Layer
Input Channels
Output Channels
Linear Layer
101
32
ReLU
N/A
N/A
Linear Layer
32
32
ReLU
N/A
N/A
Linear Layer
32
Table 18: Network structure of Regressor
Operation Layer
Input Channels
Output Channels
Linear Layer
101
40
Softplus
N/A
N/A
Linear Layer
40
40
Softplus
N/A
N/A
Linear Layer
40
101
Table 19: Network structure of NODE
22

---

Published as a conference paper at ICLR 2023
B LIMITATIONS
There are several limitations to the proposed method. The MoC only applies to hyperbolic PDEs, and only consider first-order semi-linear PDEs in this paper. This may be a limitation since this is specific class of PDEs that does not model all data. We also did not enforce any particular structure to prevent characteristics from intersecting, which may result in shock waves and rarefactions. However, we believe that this is unlikely to happen due to the high dimensionality of the ambient space. As noted in the experiments section, C-NODE have decreased stability compared ANODE when defined by exceeding an extreme threshold of number of function evaluations.
C APPROXIMATION CAPABILITIES OF C-NODE
Proposition C.1 (Method of Characteristics for Vector Valued PDEs). Let u(x), rx) R T be the solution of a first order semilinear PDE a bounded domain Q C RK of the form
du
Σ Eai(1,, Tki u) dr.
c(21, c(21,-.,, Tk, u) on (x1,...,, vk) XEΩ.
(12)
Additionally, let a (ai,...,,.) functions. Define a system of ODEs as
)²
Bk+n - R Te - Ra be Lipschitz continuous
dx
(s) a(x(s), U(s)
de
du (s) c(x(s), U(s)
x(0) = Xo, an
u(xo) = uo
U(0) = up
where Xo and Un define the initial condition, ON is the boundary of the domain Q. Given initial
conditions Xg, Ug. the solution of this system of ODEs U(s) [a, + Rd is equal to the solution of the PDE in Equation equation [quation[22llon
the characteristic curve defined by x(s), i.e. u(x(s) s)) = Us).
The union of solutions U(s) for all Xg E 00 is equal to the solution of the original PDE in Equation equation 72 for all x S2
Lemma C.2 (Gronwall's s Lemma (Howard 1998)). Let U C R be an open set. Let f:r::: x (0, T
Rn be a continuous function and lei h1, h T] - U satisfy the initial value problems:
dh,(t) [(h, ((,,(1).). h2(0) x1:
dt
dh2(t) (h2(1).)) 1), h2(0) X2-
there exists non-negative constant C such that for all [0, [[f(hz(t), ((r f) < C'I h (t)l,,,
where || is the Euclidean norm. Then, for all t E [0, T)
[l2(()) h2)ll c x1l1
C.1 PROOF OF PROPOSITION
This proof is largely based on the proof for the univarate case provided a We extend for the vector valued case.
Proof. For PDE on with & input, and an n--imennsionll output, we have a, RR+= + R, and c pk+n + R", In proposition C.I
Ou day
E R
we look at PDEs in the following form
Ou
a (1)
c(z),...,, Lk. u).
(13)
Ox:
for_
quasilines
23

---

Published as a conference paper at ICLR 2023
Defining and substituting (z .
(a1 (.1,...,a)) az)", and Jacobian J(u(x))
0u
3u
E Ru*k into Equation equation 12 result in
[u((x))a(xx u) c(x, u).
(14)
From proposition C.I
the characteristic curves are given by
dx. a.. u)
ds
and the ODE system is given by
ds () = a(x(s), ((s)),
(15)
dU
ds (s) c(x(s). U()).
(16)
Define the difference between the solution to equation 16 and the PDE in equation 12
(s) = |lu(x(s)) Us)] (u(x(s)) U(s))* (u(x(s)) U(s)),
Differentiating ((s) with respect to s and plugging in equation we get
(' = dA(s)
ds 2(u(x(s)) - U(s) (J(ux' (s) U'(s)
2[u(x(()) - U(s] [(u)a(((s)) (s), U(s)) c(x(s), U(s))].
(17)
- c(2 Fk. u) 0. Plugging this equality into equation 14 gives us E 0z(22 (211.... u au equation 17 and rearrange terms, we have
0z:
A'ss) 2fu(x(s)) - U(s) ([())a(((s)) U()) c(((s) U(s)] [[())((((s)))))) c(x(s),
Combining terms, have
A 2(u U) ([Ju))(UU) (())a() ccu))
2(u (J(u) [a(U) a(u)] + (c(U) c(u)))
Applying triangle inequality, we have
||'l|| < 2 (
a(u)l + ||c(()) - c(u)ll).
By the assumption in proposition C.I and c are Lipschitz continuous. By Lipschitz continuity. we have [a(U) a(u))ll < Al|u U| and |c(U) - c(u)) < B B-u--Ujj U|]. for some constants co, since both and Jacobian J are continuous
and B in R_.. Also, for compact set [0, su]. mapping. J(u) also compact. Since subspace of IR" is compact and only it is closed and bounded, Juu is bounded (Strichartz 2000). Thus, [[(i)]] M for some constant M in R+.
Define C = 2(AM + B), we have
[A(ss)] < 2(AM]|u U|| + Bl U - Ull
1-l--o -
-s)((s)]]-
From proposition C.I we have u(x(0)) U(0). As proved above, we have
du(x(s))
dU(s) = (A)))) < C+)A())I
ds
ds
where C' < oc. 80. Thus, by lemma
we have
|A(s) < ( (CI]((0)] e [[(()))) U(O) 0.
This further implies that U(s) u(x())) so long and are Lipschitz continuous.
24

---

Published as a conference paper at ICLR 2023
1.0
0.8
0.0
0.6
0.4
U
0.2
0.8
1.0
Figure 10: Comparison of C-NODEs and NODEs. C-NODEs (solid blue) learn family of integration paths conditioned on the input value, avoiding intersecting dynamics. NODEs (dashed red) integrate along ID line that is not conditioned on the input value and can not represent functions requiring intersecting dynamics.
C.2 PROOF OF PROPOSITION
Proof. Suppose have C-NODE given by du ds dr ds Write out specific functions for these terms to match the desired properties of the function. Define initial condition u00 0) = ug. By setting dr (s, uo. 0) L, ds
Ou dr
du dt Øt ds
dt (s, uo 0) uu.
Ou
Ou (u(2r, (), 22
(w(x; 1),0)
Øt
have the ODE and solution,
du
1 2uo
ds
u(s; un) (1
2uo)
0
1 - 0
-1
To be specific, we can represent this system with the following family PDEs:
Ou
Ou
+
2uo-
dr
Ot
We can solve this system to obtain a function that has intersecting trajectories. The solution is visualized in Figure [0] which shows that C-NODE can be used to learn and represent this function G. should be noted that this not the only possible solution to function G, as when 0t/Os O, we fall back to NODE system with the dynamical system conditioned the input data. In this conditioned setting, we can then represent G by stopping the dynamics at different times t as in Massaroli et al. 2021
C.3 PROOF OF PROPOSITION
The proof uses the change of variables formula for particle that depends vector rather than scalar and it follows directly from the proof given (Chen et al. 2019b] Appendix A). We provide the full proof for completeness.
25

---

Published as a conference paper at ICLR 2023
Proof. We initially assume that Ex= Ou dag is Lipschitz continuous in and continuous in so every initial value problem has unique solution (Evans 2010). We additionally assume u(s) is bounded.
dz;
We show that the probability flow satisfies
Op(u(s))
du dx; dzg ds ==1
tr
0s
Ou
Define T. u(s + c). The discrete change of variables states that u1 f(o) log log p(uo) log det 0f (Rezende Mohamed 2015).
Bug
Taking the limit of the time difference between ug and u1, and using the definition of derivatives,
Olog p(u())
log logp(u(s -lgg (og(((()))
lim
Ot
logp(u(s)) log det 2 T.(((t)) log p( u(s))
lim
c+0
log det T.T((u(s))
Bu
lim c->
log det a Te (u()))
lim
Du
e-+0
dc
a log det
à
lim
log
det
T (u((s))
- lim
(u())))
de
du
c- de
du
lim
det
(u())))
<-00 det &
(u(ss) Du
dE
Du
lim,-ol
de det du
8 Tu((s)))
lim det 0 Tu(ss) du
lim
det T((((s))
c-0
dc Su
The Jacobi's formula states that if A is a differentiable map from the real numbers to x matrices, then de det A(t) = (adj(A(t)) dA(t) ). where adj denotes the adjugate matrix. Applying this relation, obtain
i Ologp(u(t))
lim s00
adj
o I..(((())) du
a a
I,(u(s))
de du
(2)r.(((1))
a 2 de Du T.(((s)))
tr
lim adj
lim
6
Du
c-+
tr
adj
Ou u(()
lim
de tu T.(((s))
+
tr
lim
T. (.((ss))
De
Ou
26

---

Published as a conference paper at ICLR 2023
Substituting T. with its Taylor series expansion and taking the limit, we have the desired result
Ölog p( (u(t)
Ø
0
du
tr
lim -I
u
+
+((2)) + O(cc) +((2)+....
0t
de du
ds
0 à
Du dr;
+0(2) + +(())) +..))
tr
lim
+ c
1 de du
Dz, ds
im1
lim
i
+
du dx;
+0())+0(0)+
de
Du
Oz; ds
0
d
Du dæ
+0(e) 0 + +(()++
tr
lim
c-0++
vu
vx; ds
k
Tu dx;
tr
Ou
Ozy ds
1=1
C.4 PROOF OF PROPOSITION
4.3
Proof. To prove proposition 4.3 we need to show that for any homcomorphism h(.), there exists a u(s, uo) E R following C-NODE system such that u(s = T, uo) h(ua)-
Without loss of generality, we suppose that = 1.
Defining C-NODE system as:
du Su de +
Su dt at de
da du
dz (s, ug) 1,
du (u(x, t)) h (uo))
8x dt
ds (s, up) = 140
Du
37 (u(I, = -1
Then, du = h(wo) ug. At 1, have
du d
u(s = 1, un) = uu(s = 0, uo)
ds
uo + L Du dx
du dt
ds
Fx ds
0t ds
/ h(wo))
uo +
1 + (-1) unds
uo + h((w)) - uwo
>h(a))-
The inverse map will be defined by integration backwards. Specifically, we have
u(s 0. un) u(s 1, 1,u0) { du ds ds
Ou dr
Ou dl Øt ds
((uo))
ds
Or ds
>(uu))
h(uo) 11 (-1) unds
h((u))
h(uo)) + uo
ug.
27

---

Published as a conference paper at ICLR 2023
Thus, for any homeomorphism h(-), there exists C-NODE system, such that forward integration for time 8 = is equivalent as applying ht)), and backward integration for time s = 1 is equivalent to applying (-).
D INCLUDING INITIAL STATE IN NODE'S INPUT
Conditioning on the initial condition allows NODE to model intersecting trajectories. We compare NODE conditioned on initial values and C-NODE's performance on image classification tasks on MNIST, SVHN, and CIFAR-10 10 datasets, using the forward Euler solver and adjoint solver.
As shown in Figure 11 by conditioning on initial values, NODE's performances are improved over all datasets when using adjoint integration, and improved SVHN dataset when using Euler forward integration. When using Euler forward integration, C-NODE performs better on all three datasets than NODE and conditioned-NODEE. C-NODE higher accuracy within fewer epochs and has lower variance throughout the training process. All methods use a comparable numbers of parameters with C-NODE using the fewest as reported in Table 20 When using adjoint method, C-NODE performs better on all datasets than NODE, and performs better on SVHN and CIFAR-10 datasets than Mnditioned-NODEE while achieving similar performance to conditioned-NOO NODE.
It is worth noting that the characteristic curves of the C-NODE also depend on the initial values, which allows C-NODE to model a dynamical system with intersecting trajectories. C-NODE and conditioned NODE are two different methods of including the initial values. Empirically, C-NODE achieves better results and is more efficient with parameters, as suggested Figure and and Table 20
Dataset
Method
Accuracy (Euler) f
Accuracy (Adaptive) f
Param. [K] 4
NODE
82.42 = 0.043
64.510 = 1.637
115.444
SVHN
COND
83.30 = 0.131
80.39 + 1.088
116.399
C-NODE
86.28 = 0.086
82.35 + 1.458
113.851
NODE
59.492 + 0.366
52.13 + 1.464
115.444
CIFAR-I0
COND
56.77 = 0.137
57.71 = 0.226
118.447
C-NODE
61.97 = 0.139
63.57 + 0.421
113.851
NODE
97.09 0.024
96.79 + 0.374
85.468
MNIST
COND
93.38 = 1.798
98.05 + 0.112
87.287
C-NODE
96.93 = 0.103
97.26 + 0.276
83.041
Table 20: Parameter counts and classification accuracy for different models and integration schemes in Conditional-NODE experiment.
All experiments were performed on NVIDIA RTX 3090 GPUs on a cloud cluster.
E ABLATION STUDY
E. ABLATION STUDY ON DIMENSION OF C-NODE
We perform an ablation study on the impact of the number of dimensions of the C-NODE we implement. This study allows us to evaluate the relationship between the model performance and the model's limit of mathematical approximating power. Empirical results show that as we increase the number of dimensions used in the C-NODE model, the C-NODE's performance first improves and then declines, due to overfitting. We have found out that information criteria like AIC and BIC can be successfully applied for dimension selection in this scenario.
In previous experiments, we represent du/Ox; with separate and independent neural networks (u, 6). Here, we represent all & functions as a vector-valued function approximate this vector-valued function with neural network c(u, the Euler solver to have better training stability when the neural network has a large number of parameters. Experiment details for the ablation study is as shown in Figures 12
"h/02z i". We 8). The model is trained using
13 14
28

---

Published as a conference paper at ICLR 2023
60
A00074
20
2.2
1.00
200
1.79
2.0
15
2.25
2.0
15
C196
MAIS1
54-9
MhI5T
(a) Backprop through Fuler training;
(b) Adjoint training;
Figure 11: Red: NODE. Blue: C-NODE. Orange: Conditional NODE Training dynamics of different datasets with adjoint in Fig- and with Fuler in Fig. 1Tb averaged over four runs. The first column the training process of SVHN, the second column is of CIFAR-10, and the third column is of MNIST. By conditioning on the initial values, NODE's performances improved on all datasets when using adjoint integration, and improved on SVHN dataset when using the Euler backpropagation.
E.2 ABLATION STUDY ON NUMBER PARAMETERS
We show C-NODE's parameter efficiency over NODE with an ablation study on the image classi- fication task the CIFAR-10 dataset. Specifically, under similar training setup, we experiment with C-NODE with 95071, 1, 55855, and 17379 parameters and experiment with NODE with 96044, 56828, and 17444 parameters. As shown in Figure [3] although C-NODE has more variance in its performance, it outperforms NODE along the whole training process in all three cases.
29

---

Published as a conference paper at ICLR 2023
90
Accuracy
80
70
16D
64D
32D
128D
60
1D
4D
16D
64D
2D
8D
32D
-22DD
Test Loss
1.8
1.6
2.0
1D
4D
16D
64D
Train Loss
2D
8D
32D
128D
1.8
1.6
20
40
60
80
100
Figure 12: The training process averaged over 4 runs of C-NODE with 1, 2, 4, 8, 16, 32, 64, 128. 256, 512, and 1024 dimensions the MNIST dataset. The first row is the accuracy of prediction, the second row is the testing error, and the third row the training error.
30

---

Published as a conference paper at ICLR 2023
85
Accuracy
80
75
1D
4D
16D
64D
2D
8D
32D
128D
1.75
1D
4D
16D
64D
Test Loss
1.70
2D
8D
32D
128D
1.65
1.60
1.7
1D
4D
16D
64D
Train Loss
2D
8D
32D
128D
1.6
1.5
20
40
60
80
100
Figure 13: The training process averaged over runs of C-NODE with 1, 2, 4, 8, 16, 32, 64, and 128 dimensions on the SVHN dataset. The first row is the accuracy of prediction, the second row the testing error, and the third row is the training error.
31

---

Published as a conference paper at ICLR 2023
60
Accuracy
50
40
1D
16D
64D
2D
32D
128D
30
1D
4D
16D
64D
2.1
2D
8D
32D
-128D
Test Loss
2.0
1.9
2.2
1D
4D
16D
64D
Train Loss
2D
8D
32D
128D
2.0
1.8
1.6
20
40
60
80
100
Figure 14: The training process averaged over runs of C-NODE with 1, 2. 4, 8, 16, 32, 64, and 128 dimensions the CIFAR- 10 dataset. The first row is the accuracy of prediction, the second row is the testing error, and the third row is the training error.
32

---

Published as a conference paper at ICLR 2023
60
60
Accuracy
40
40
40
20
20
20
0
20
40
0
20
40
20
40
Test Loss
2.2
2.2
2.2
2.0
2.0
2.0
0
20
40
20
40
20
40
Train Loss
2.25
2.25
2.25
2.00
2.00
2.00
1.75
1.75
1.75
20
40
20
40
40
95k
56k
17k
Figure 15: The training process averaged over four runs of C-NODE with 95071, 55855, and 17379 parameters on the CIFAR-10 10 dataset, and NODE with 96044, 55855, and 17379 parameters. The first row is the prediction accuracy, the second row is the testing error. and the third row is the training error. Blue lines are the results for C-NODE, and red lines the results for NODE.
33

---

Published as a conference paper at ICLR 2023
Algorithm 2 C-NODE algorithm using the forward Buler method
for each input data 2; do
extract image feature u(s - 0) - g(2z; 1) with a feature extractor neural network. procedure Integration along = 0 +
for each time step 9m do
calculate dx (x, u g Zy; 11); 22) and Jxu(x, O2). Jxu
da
calculate calculate (1m++1) = ((m)) + du (8m+1 - 8m
ds
ds
calculate x(5m+1) = x(sm) + (8m+1 -8m)
dx
end for end procedure
classify u(s 1) with neural network ((u(x((s = 1)), 02).
end for
F ALGORITHM FOR IMAGE CLASSIFICAITON WITH C-NODE
We provide algorithm for training C-NODE for the purpose of image classification.
G ALGORITHM FOR CONTINUOUS NORMALIZING FLOWS DEFINED WITH C-NODE
We additionally provide algorithms for training and sampling CNFs defined with C-NODEs.
Algorithm 3 Algorithm for training CNFs defined with C-NODE
given probability density function of p(u(s 0)) Po(-)
for each input data 2j do
u11)
3j
Given
logp(z,) - log p(u())
u(0)
procedure Integrate from + to get
log (oz)( - log p p(u()))
for each time step 8m. do
calculate dx da (x, u: g(y; 11): 22) and Jx u: O2). du de = J de tr au J de with Hutchinson trace estimator (Grathwohl al. 2019
calculate calculate
dx
u1)
u(sm) log p(u(mm))
du
calculate
(ogg(u(sm+1))
0 log pl p(u()))
(sm-1 sm).
Ge
end for
evaluate Pu(ul (u(0))
calculate log p(z,) = (log p(z,) log p(u(0)))
+ log (u())))
optimize p(z,)
with an optimization algorithm (stochastic gradient descent etc.)
end for
34

---

Published as a conference paper at ICLR 2023
Algorithm 4 Algorithm for sampling CNFs defined with C-NODE
procedure sample u(s = 0) from base distribution po(+)
procedure Integrate from - 1 to get u(s = 1)
for each time step 8m. do
calculate dx
ds (x, u; gggz((z Zy; O ); O2) and Jxu(x, u: O2).
calculate du Jxu
dx
de calculate u(sm+1) = u(mm) + du end for end procedure
ds
(sm+1 - sm).
ds
u(s 1) our sample from the CNF
35


Reviews:
[{"review_id": "mi9yIrD0Rd", "reviewer": "Reviewer_212J", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "Partial differential equations (PDEs) are ubiquitous in science and engineering, and modern deep learning frameworks for modeling PDEs are of interest for scientific machine learning. This paper leverages the method of characteristics to extend the framework of Neural Ordinary Differential Equations (Neural ODEs) to modeling (hyperbolic) PDEs. The method of characteristics allows to transform hyperbolic PDEs into ODEs. Then, tools for learning Neural ODEs can be used to model the transformed PDE. Several different problems are used to demonstrate the performance of the proposed Characteristic Neural ODE framework.", "summary_of_the_review": "In summary, the presented ideas are solid, and I feel that the innovations are novel and significant. The proposed framework connects Neural ODEs to an interesting class of PDEs that are relevant in many real-world applications. Hence, I have no doubt that this work has the potential to inspire future research in this area. Overall, this paper is of comparable quality compared to several related works that have been recently been accepted at ICLR and other ML conference.  Given that the reviewers address my comments, I recommend to accept this paper.", "strength_and_weaknesses": "This paper introduces a novel and relevant extension for the existing Neural ODE framework, and the presented experimental results are significant. The proposed framework allows to model a broad class of interesting PDEs, i.e., hyperbolic PDEs. I don’t feel that it is a restrictive limitation of the framework that it can handle only hyperbolic PDEs, since this class is of general interest in many problems.  \n\nWhile the key innovation of this paper is a new deep learning method for PDE approximation, my concern is that this framework could just be viewed as a version of augmented NODEs. However, the experiments seem to suggest that the proposed method gives good results. To me, the strength here is in the interpretation, but I'm guessing from the PDE experiment that you cannot get the NNs to represent what their analytical solution counterparts should be. Still, I think that this paper is a good step in the right direction.  Here are things that I would like to see addressed, and that I feel will help to improve the paper: \n\n\n* The presentation of the mathematical exposition can be streamlined in several places to make the paper better readable. I feel that sometimes some unnecessary terms/notations are introduced without additional explanation. For instance, the middle term in the first equation of the related work is more confusing than helpful. I don’t get why it is important to explicitly state that the right-hand side is a function of (u,t). If this notation is important to understand later concepts, then it should be mentioned / discussed with a sentence within the text. I would try to keep things simple so that the paper is more accessible to a broader audience.  \n\n\n* The related work is discussing that a residual unit can be modeled as an ODE in the limit. While in theory that is fine, this statement is problematic in practice. Typically, you can’t choose $\\Delta t$ arbitrarily small, if you want to train a model in finite time. Hence, while there is some algebraic resemblance, the connection between resnets and the forward Euler discretization of an ODE is problematic. I would suggest citing works that point out some of these issues: [1,2]. \n\n\n* Since the presented approach is a method for PDE approximation, I would suggest devoting more room for experimental results in the main body of the manuscript to illustrate the framework, i.e., move some of the experiments from the appendix to the main text and expand the discussion. Additional experiments for PDEs in the appendix would be welcome. I am saying this, because I don’t feel that the proposed method is of interest for image classification. Also, the image classification experiments ignore mention of state-of-the-art Neural ODE frameworks [3,4]. For instance [4] achieves about 90% test accuracy on CIFAR-10 with a similar parameter count, and 99.6% test accuracy on MNIST with only 18k parameters. These works should at least be cited in the related works.  \n\n\n* The quality of the plots can be improved. Instead of having solid confidence bands, I suggest using transparent confidence bands so that it can be seen how large the overlap is. \n\n\n* Lastly, I would add a short discussion about the differences / advantages of the proposed approach compared to physics-informed neural networks. I think that the proposed approach is more general for modeling PDEs. So, it would be nice to point this out.  \n\n \n\n[1] Ott, Katharina, et al. \"Resnet after all: Neural odes and their numerical solution.\" International Conference on Learning Representations. 2020. \n\n[2] Krishnapriyan, Aditi S., et al. \"Learning continuous models for continuous physics.\" arXiv preprint arXiv:2202.08494 (2022). \n\n[3] Zhang, Tianjun, et al. \"ANODEV2: A coupled neural ODE framework.\" Advances in Neural Information Processing Systems 32 (2019). \n\n[4] Queiruga, Alejandro, et al. \"Stateful ODE-Nets using basis function expansions.\" Advances in Neural Information Processing Systems 34 (2021): 21770-21781. ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "iK3bM9LsN2o", "reviewer": "Reviewer_TVDi", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The manuscript introduces a neural ODE (NODE) type model inspired by the method of characteristics from PDE theory, which is therefore referred to as characteristic NODE or C-NODE. Compared to vanilla NODEs the resulting parametric dynamical system can exhibit intersecting dynamics, which can be desirable for certain applications.\nThis result is strengthened by showing that the proposed model is more expressive than vanilla ODEs in the sense that it can compute arbitrary homeomorphisms as input-output maps. \nFurther, a formula for the evoluation of the log density under C-NODEs is presented, which makes the model applicable to density estimation problems. \nThe proposed model is tested on image classification, time series prediction, density estimation as well as approximation of the solution of a quasilinear hyperbolig PDE.\n", "summary_of_the_review": "The manuscript proposes an interesting new NODE type model, which is inspired by the method of characteristics for hyperbolic PDEs. \nOverall, I believe that the idea is intriguing and promising. In its current form, I am mostly concerned about the clarity of the paper (the introduction of the model is sloppy or at least not very clear) as well as the reproducability. I believe that both of these points should be addressed prior to a publication of the manuscript.", "strength_and_weaknesses": "**Strengths:**\n* The proposed model is based on an important theoretical tool in the theory of PDEs. \n* The expressivity of the proposed model is studied theoretically\n* The experiments showing performance gains of the model cover a variety of different problems and include comparisons against a variety of different methods.\n\n**Weaknesses:**\n\n* On the side of the experiments, I believe the computation time should be added for a reference, since the number of parameters does not completely describe the complexity of training. Further, in the limitations it is stated that the proposed model exhibits instabilities during training. A discussion regarding this should already present in the experiments section. \n\n* The experiment on the representational power when it comes to approximating solutions of PDEs seems not well desgined. As far as I see it, C-NODEs are shown to exhibit lower generalization error compared to NODEs when used to regress the solution of a hyperbolic PDE. This however does not imply that C-NODEs are more effective in approximating the regressed function. \n\n* One weakness of the paper is the motivation to use a model inspired by hyperbolic PDEs. The main argument given is its ability to „transport data into different regions in state space“, which is also highlighted by the theoretical expressivity results. Where I agree that this is indeed a favorable property, this is already the case for simple augmented NODEs. To me it is not clear whether the improvements in the experiments are actually due to an underlying hyperbolic structure in the data of the experiments. That being said, I think this is fine, since the manuscript only speculates about this and does not offer this as an explanation.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "CxP6TASxJ9", "reviewer": "Reviewer_eEjp", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes to model the evolution of latent variables in continuous time models as the solution of first-order partial differential equations. The authors builds upon the method of characteristics for solving PDEs and show that their approach provides both theoretical and practical advantages. In particular, they show that their method can model interesecting trajectories and can still benefit from the adjoint method for backpropagation. Experimentally, the method outperforms classical NODE and ANODE in a variety of tasks (classificiation, time series classification as well as NF). ", "summary_of_the_review": "This work proposes and interesting idea from PDE solving to improve NODE architecture. The experiments are broad in scope. I still have some concern regarding a mismatch between the experiments and the motivations of the paper : regarding the adjoint and the usage of the first input in the dynamics. Some work on the motivation and comparison with more relevant baseline (the one I suggested above) would be appreciated. ", "strength_and_weaknesses": "## Strengths\n\n- The method is nicely presented and the overall paper is clear. The method is also well motivated, using the method of characteristics for solving PDEs. \n\n- The experimental section is solid, with various setttings (spanning from classification, to time series modeling and normalizing flows). The results indicate improved performance of C-NODE.\n\n## Weaknesses\n\n\nFrom the paper, I read two motivations : the ability to keep the adjoint method for training and the ability to model more complex dynamics than NODE (in particular intersecting trajectories). Yet, I think these are currently under-developped in the paper.\n\n*About the adjoint*\n\n- For the adjoint method, from experience, the improved memory complexity comes at the cost of significantly higher computational time, making it therefore unpractical. Therefore, this motivation is quite weak in my opinion. What is more, the cases used in the experimients don't seem to be requiring the adjoint. method for training either. I would therefore recommend to make this adjoint argument more compelling by showing a case study where the adjoint is indeed necessary.\n\n*About the ability to model intersecting trajectories*\n\n- The ability to model intersecting trajectories seem to arise from the fact that you include the first time point $g(\\mathbf{z})$ in the input of the dynamical system function. In my sense, that alone suffices to allow for intersecting trajectories. If you agree, this comes, in my sense, with multiple implications. First, that statement should be made more clear in the text so the reader can better see where this ability comes from. Second, this would represent the most interesting baseline to compare against. Namely, you take a neural ODE and you feed the first $\\mathbf{z}$ as input, similiarly as you do for the dynamics of the characteristics. This would be, in my opinion, a more compelling argument for using the method of characteristics framework. \n- What is more, the non-intersecting property clashes with the invertibility statement of section 4.2. Indeed, I do not see how intersecting can be invertible as they won't be injective. I think this would deserve some clarification. I believe this confusion is due to what you consider is the input (if you condition on $z$, it is invertible$, but not otherwise).\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "5vtNtT1Bo4", "reviewer": "Reviewer_74Jk", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors proposed Characteristic Neural ODE (C-NODE), which extends the idea of Neural ODE to combine with PDEs, and enables the model to learn a family of models (rather than a single Neural ODE model).\n\nFurthermore, the authors demonstrated that C-NODE has a nice property, that it can model intersection trajectories (similar to augmented Neural ODE), in contrast Neural ODE losing such capability. The authors also derived the formulation for continuous normalizing flow with C-NODE (in the PDE settings).\n\nFinally, the authors validated the performance of C-NODE with experiments.", "summary_of_the_review": "The paper is novel, well presented and solid in theory. Though a little bit weak in experiment, but it could come from that the model size is small. Overall is a very good paper.", "strength_and_weaknesses": "Strength:\n1. The idea to extend Neural ODE to PDE via Characteristic equations is novel. It great expands the application of Neural ODEs.\n2. The theoretical derivation is solid. (Though I only loosely checked the proofs the theories, they look pretty convincing to me).\n3. The paper is well-presented and easy to follow, especially in the method section, the authors started from the general background of characteristic equations for PDE, then extend them to the Neural ODE case.\n\nWeakness:\n1. Could the authors expand a bit more on the difference between augmented NODE and C-NODE? The formulation of C-NODE essentially models the original NODE in a higher dimension, which looks close to augmented NODE. So I wonder what is the key difference.\n2. Experiments\n    Despite the nice theory and presentation, the paper lacks comparison with some stronger baselines. For example, [1,2,3] achieved slightly better results on CNF on Cifar10, though they are focused on numerical solvers rather than PDE, it's worth mentioning in the paper. \n\nPS: I doubt if the slightly worse result of C-NODE is because the model is smaller than in the literature. I would be happy to increase score if the author show better results with a model of comparable size.\n\n[1] Zhuang, Juntang, et al. \"Adaptive checkpoint adjoint method for gradient estimation in neural ode.\" International Conference on Machine Learning. PMLR, 2020. \n[2] Matsubara, Takashi, Yuto Miyatake, and Takaharu Yaguchi. \"Symplectic adjoint method for exact gradient of neural ode with minimal memory.\" Advances in Neural Information Processing Systems 34 (2021): 20772-20784.\n[3] Zhuang, Juntang, et al. \"MALI: A memory efficient and reverse accurate integrator for Neural ODEs.\" arXiv preprint arXiv:2102.04668 (2021).", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster