Neural Delay Differential Equations
Authors: Qunxi Zhu
Yao Guo
Wei Lin
Venue: ICLR 2021 Poster (2021)
Keywords: Delay differential equations
neural networks

Abstract:     Neural Ordinary Differential Equations (NODEs), a framework of continuous-depth neural networks, have been widely applied, showing exceptional efficacy in coping with some representative datasets.  Recently, an augmented framework has been successfully developed for conquering some limitations emergent in application of the original framework.  Here we propose a new class of continuous-depth neural networks with delay, named as Neural Delay Differential Equations (NDDEs),  and, for computing the corresponding gradients, we use the adjoint sensitivity method to obtain the delayed dynamics of the adjoint. Since the differential equations with delays are usually seen as dynamical systems of infinite dimension possessing more fruitful dynamics, the NDDEs, compared to the NODEs, own a stronger capacity of nonlinear representations.  Indeed, we analytically validate that the NDDEs are of universal approximators, and further articulate an extension of the NDDEs, where the initial function of the NDDEs is supposed to satisfy ODEs.   More importantly, we use several illustrative examples to demonstrate the outstanding capacities of the NDDEs and the NDDEs with ODEs' initial value.  More precisely, (1) we successfully model the delayed dynamics where the trajectories in the lower-dimensional phase space could be mutually intersected, while the traditional NODEs without any argumentation are not directly applicable for such modeling, and (2) we achieve lower loss and higher accuracy not only for the  data produced synthetically by complex models but also for the real-world image datasets, i.e., CIFAR10, MNIST and SVHN.   Our results on the NDDEs reveal that appropriately articulating the elements of dynamical systems into the network design is truly beneficial to promoting the network performance.

Published as a conference paper at ICLR 2021
NEURAL DELAY DIFFERENTIAL EQUATIONS
Qunxi Zhu, Yao Guo & Wei Lin
School of Mathematical Science, Research Institute of Intelligent Complex Systems and ISTBI State Key Laboratory of Medical Neurobiology and MOE Frontiers Center for Brain Science Fudan University
Shanghai 200433, China (gxzhul 6, yguo, wl injn fudan. aadan.edds edu. cn
ABSTRACT
Neural Ordinary Differential Equations (NODEs), framework of continuous- depth neural networks, have been widely applied, showing exceptional efficacy coping with some representative datasets. Recently, augmented framework has been successfully developed for conquering some limitations emergent applica- tion of the original framework. Here we propose new class of continuous-depth neural networks with delay, named as Neural Delay Differential Equations (ND- DEs), and, for computing the corresponding gradients, tivity method to obtain the delayed dynamics of the adjoint. Since the differential equations with delays usually seen as dynamical systems of infinite dimension possessing more fruitful dynamics, the NDDEs, compared to the NODEs, own a stronger capacity of nonlinear representations. Indeed, we analytically validate that the NDDEs are of universal approximators, and further articulate an exten- sion of the NDDEs, where the initial function of the NDDEs is supposed satisfy ODEs. More importantly, we use several illustrative examples to demonstrate the outstanding capacities of the NDDEs and the NDDEs with ODEs' initial value. Specifically, (1) we successfully model the delayed dynamics where the trajec- tories in the lower-dimensional phase space could be mutually intersected, while the traditional NODEs without any argumentation are not directly applicable for such modeling, and (2) we achieve lower loss and higher accuracy not only for the data produced synthetically by complex models but also for the real-world image datasets, i.e., CIFARIO, MNIST, and SVHN, Our results on the NDDEs reveal that appropriately articulating the elements of dynamical systems into the network design is truly beneficial to promoting the network performance.
the adjoint sensi-
INTRODUCTION
A series of recent works have revealed a close connection between neural networks and dynamical
systems E 2017
Li et al. 2017 Haber & Ruthotto [017]] Chang al. 2017 Li & Hao 2018 Lu
et al. 2018 et al. 2019 Chang et 2019; Ruthotto & Haber 2019 Zhang et al. 2019a Pathak
al. 2018 Fang al.) 2018 Zhu al.) 2019: Tang al. 2020). On one hand, the deep neural networks can be used to solve the ordinary/partia differential equations that cannot be easily com- puted using the traditional algorithms. On the other hand, the elements of the dynamical systems can be useful for establishing novel and efficient frameworks of neural networks. Typical examples include the Neural Ordinary Differential Equations (NODEs), where the infinitesimal time of ordi- nary differential equations is regarded as the "depth" of a considered neural network (Chen et al. 2018
Though the advantages of the NODEs were demonstrated through modeling continuous-time datasets and continuous normalizing flows with constant memory cost (Chen et al., 2018)) the lim- ited capability of representation for some functions were also studied (Dupont et al. 2019). Indeed, the NODEs cannot be directly used to describe the dynamical systems where the trajectories in the phase space mutually intersected. Also, the NODEs cannot model only few variables from some physical or/and physiological systems where the effect of time delay is
"htp://homepage. www.effudducccccco edw.cn/weilin, To whom correspondence should be addressed: Q.Z., Y.G., and W.I..

---

Published as a conference paper at ICLR 2021
inevitably present. From view point of dynamical systems theory, all these attributed to the characteristic of finite-diment dimensiimnnionn for the NODEs.
In this article, propose novel framework of continuous-depth neural networks with delay, named as Neural Delay Differential Equations (NDDEs), We apply the adjoint sensitivity method to compute the corresponding gradients, where the obtained adjoint systems are also in form of delay differential equations. The main virtues of the NDDEs include:
feasible and computable algorithms for computing the gradients of the loss function based the adjoint systems,
representation capability of the vector fields which allow the intersection of the trajectories in the nower-dimensione. phase space, and
accurate reconstruction of the complex dynamical systems with effects of time delays based on the observed time-series data.
2
RELATED WORKS
NODEs Inspired by the residual neural networks (He et al. 2016) and the other analogous frame- works, the NODES were established, which can be represented by multiple residual blocks as
ho+1 he f(hr, wx),
where h; is the hidden state of the t-th layer. fh,; w,) is differential function preserving the dimension of h, and w, the parameter pending for learning. The evolution of h, can be viewed the special case of the following equation
h441 he At f f(h,, wi)
with At = 1. As suggested in (Chen et al. achieving parameter efficiency of the NODES. This unified operation was also employed in the other neural networks, such as the recurrent neural networks (RNNs) (Rumelhart et al., 1986; Elman 1990) and the ALBERT (Lan et al. 2019). Letting At + and using the unified parameter w instead of w,, we obtain the continuous evolution of the hidden state h, as
2018), all the parameters are unified into w for
h at h dh;
lim A1-0
th dt
At
which is in the form of ordinary differential equations. Actually, the NODEs can act as feature extraction, mapping an input to a point in the feature space by computing the forward path of NODE as:
h(T) h() + [ f(e,, 1; 1, h()) input,
where h20) input the original data point or its transformation, and T the integration time (assuming that the system starts at f 0).
Under predefined loss function L()())), (Chen et al., 2018) employed the adjoint sensitivity method to compute the memory-efficient gradients of the parameters along with the ODE solvers.
More precisely, they defined the adjoint variable, A(t) aL(M(T)) whose dynamics another ODE, i.e.,
oh())
da(1) -X()) Ōs -s(h,, (:),,;:w) dt
(1)
Oh,
The gradients are computed by an integral as:
dw - ff- -X(t) Of t: w)
dL
dt.
(2)
dw
Chen et al., 2018) calculated the gradients by calling an ODE solver with extended ODEs (i.e., concatenating the original state, the adjoint, and the other partial derivatives for the parameters at each time point into a single vector). Notably, for the regression task of the time series, the loss function probably depend the state at multiple observational times, such the form of
2

---

Published as a conference paper at ICLR 2021
[(h(to), h(,x)). Under such case, must update the adjoint state instantly by adding the partial derivative of the loss at each observational time point.
As emphasized in (Dupont et al. 2019), the flow of the NODEs cannot represent some functions omnipresently emergent in applications. Typical examples include the following two-valued func- tion with argument: (1-(1) -1 1 and (Td( = 1. Our framework desires to conquer the representation limitation observed applying the NODEs.
Optimal control As mentioned above, closed connection between deep neural networks and dy- namical systems have been emphasized in the literature and, correspondingly, theories, methods and tools of dynamical systems have been employed, e.g. the theory of optimal control (E 2017 Li
et al. 2017
Haber & Ruthotto 2017 Chang et al. 2017 Li & Hao 2018 E et al. 2019 Chang
al.
2019 Ruthotto Haber 2019 Zhang et al. 2019a). Generally,
model typical task using
deep neural network and then train the network parameters such that the given loss function can be reduced by some learning algorithm. In fact, training network can seen as solving an optimal control problem on difference or differential equations (E et al. 2019). The parameters act some objective function. typical problem of optimal control on
controller with a goal of finding an optimal control to minimize/maxine
Clearly, the framework of the NODES can be formulated ODEs. Additionally. the framework of NODEs has been generalized to the other dynamical sys- tems, such the Partial Differential Equations (PDEs) (Han et al. 2018 Long et al.) 2018 2019 Ruthotto & Haber] 2019 [un et al. 2020) and the Stochastic Differential Equations (SDEs) (Lu
20188 (Sun et al. 2018 Liu et 2019), where the theory of optimal control has been pletely established. It is worthwhile to mention that the optimal control theory is tightly connected with and benefits from the method of the classical calculus of variations (Liberzon) 2011). We also
will transform our framework into an optimal control problem, and finally solve it using the method of the calculus of variations.
3 THE FRAMEWORK OF NDDES
XP
3.1 FORMULATION OF NDDES
In this section, we establish framework of continuous-ded depth neural networks. To this end, first introduce the concept of delay defer- ential equations (DDEs). The DDEs are always written in form where the derivative of given variable at time is affected not only by the cur- rent state of this variable but also the states at some previous time instants or time durations (Erneux 2009). Such kind of delayed dynamics play an important role description of the com- plex phenomena emergent in many real-world systems, such as physical, chemical, ecological, and
-1-6
-1
Figure 2: (Right) Two continuous trajectories gen- crated by the DDEs are intersected, mapping -1 (resp., 1) 1 (resp.. -1), while (Left) the ODEs cannot represent such mapping.
Rd
Rd
Rd [-t,
Rd
L(M(T))
L(MT))
delay effect
10())
((t)
ODE f(hg; w)
loss function
DDE fhu by w)
loss function
t
E [-z,0]
t=
Figure 1: Sketchy diagrams of the NODES and the NDDES, respectively, with the initial value h(O) and the initial function Φ(t). The NODEs and the NDDEs act as the feature extractors, and the following layer processes the features with a predefined loss function.

---

Published as a conference paper at ICLR 2021
physiological systems. In this article, we consider system of DDE with single time delay:
{
f(,, (; -> 0, t <<= 0,
#t h(() ott),
where the positive constant 7 the time delay and (t)) ((t) is the initial function before the time ? = 0. Clearly, in the initial problem of ODEs, only need to initialize the state of the variable at i = while we initialize the DDEs using a continuous function. Here, to highlight the difference
between ODEs and DDEs, provide simple example:
{
dr -2t-rr df z(t)
1>= 0, t <<= 0.
To.
where r 1 (with time delay) or T (without time delay). As shown in Figure 2 the DDE flow can map 1 to and to -1; nevertheless, this cannot be made for the ODE whose trajectories are intersected with each other in the t-r space in Figure
3.2 ADJOINT METHOD FOR NDDES
Assume that the forward pass of DDEs is complete. Then, we need compute the gradients in reverse-mode differentiation by using the adjoint sensitivity method (Chen et al., [018 Pontryagin
al. 1962). We consider an augmented variable, named as adjoint and defined as
(L(h(T))
X(1) (t)
where L(-) is the loss function pending for optimization. Notably, the resulted system for the adjoint in form of DDE well.
Theorem (Adjoint method for NDDEs). Consider the loss function I(.). Then, the dynamics of adjoint can be written as
dA(t) -At) 0h t,
T) @f(+h h1, w) X(o,T--(tt,, <= T
dt
Oh:
Ōh.
A(T)= = (T)
Oh(T)
(4)
where ()) is a rypical characteristic funccttto
We provide ways prove Theorem 1 which are, respectively, shown in Appendix. Using h(t) and A(t), we compute the gradients with respect to the parameters w as:
- -xt) 0h,,, (hh,, w) dw
dL. =
dt.
(5)
dw
Clearly, when the delay r approaches zero, the adjoint dynamies degenerate as the conventional case of the NODES (Chen et al. 2018).
We solve the forward pass of h and backward pass for h. À and dL by a piece-wise ODE solver, dep which is shown in Algorithm [ For simplicity, we denote by f(t) and g(t) the vector filed of h and X, respectively. Moreover, in this paper, we only consider the initial function d(t) as a constant function, i.e., (()) = ho. Assume that T = n and denote ((x(() f(k + 1), 9u (1) g(k +1)
and >xx (1) A(k +1).
In the traditional framework of the NODES, we can calculate the gradients of the loss function and recompute the hidden states by solving another augmented ODEs in a reversal time duration. How- ever, to achieve the reverse-mode of the NDDEs in Algorithm [ we need to store the checkpoints of the forward hidden states h(i 7) for i = 0, 1,.,, n, which, together with the adjoint A(t), can help to recompute h(t) backwards in every time periods. The main idea of the Algorithm I is to convert the DDEs as a piece-wise ODEs such that one can naturally employ the framework of the NODEs solve it. The complexity of Algorithm is analyzed in Appendix.
4

---

Published as a conference paper at ICLR 2021
Algorithm Piece-wiss reverse-mode derivative of an DDE initial function problem
Input:
dynamics parameters w, time delay т, start time O, stop time = n T, final state (T"),
loss gradient (L/
0L
for in range(n - 1, -1, -1;
[h(T), (T - 7), 1((T, aL
Al 8h( (h((++)) 1)r) 8w
OL
def aug dynaaiiss(hhnnn--- (1).... An (), Aa (:..)...(11))) (1), w)
return - (t) ..... fa(t), In (t)....,9(()) - xt
8ur
aL AL ODESolle(s0- aug dynamics, T, 7,0, w)
8w
return
B1
2h(0)
8w
ILLUSTRATIVE EXPERIMENTS
4.1 EXPERIMENTS ON SYNTHETIC DATASETS
Here, we use some synthetic datasets produced by typical examples to compare the performance of the NODES and the NDDEs. In (Dupont et al. (2019), it is proved that the NODEs cannot represent the function g _ R, defined by
g(x)) = { 1,
| S r1 if S x
if
-1,
if
r3:
where 0 ri <12 < r3 and || is the Euclidean norm. The following proposition show that the NDDEs have stronger capability of representation.
Proposition 1 The NDDEs can represent the function g(æ)) specified above.
To validate this proposition, we construct a special form of the NDDEs by
(t)(t) dt h -r,
t >= 0 and 1,
(! 0,
t >> 0 and i = 2,...,. d,
h(t) x.
t <<=0 0.
where r = (r1 + /22 is constant and the final time point T is supposed to be equal to the time delay 7 with some sufficient large value. Under such configurations, we can linearly separate the two clusters by some hyper plane.
In Fig. 3 we, respectively, show the original data of g(x) for the dimension d 2 and the transformed data by the DDEs (6) with the pa- rameters Fı
12 2, r3 3, T 10, and T = 10. Clearly, the transformed data by the DDEs are linearly separable. We also train the NDDEs to represent the g(x) for d = 2, whose evolutions during the training procedure are, respectively, shown in Fig. g This figure also includes the evolution of the NODEs. No- tably, while the NODEs is struggled to break apart the annulus, the NDDEs easily separate them. The training losses and the flows of the NODEs and the NDDEs are depicted, respectively, Fig. Particularly, the NDDEs achieve lower losses with the faster speed and directly separate the two clusters in the original 2-D space; how- ever, the NODEs achieve it only by increasing the dimension of the data and separating them in a higher-dimensien space.
Figure 3: (Left) The data at the time t 0 and (Right) the transformed data at sufficient large final time of the DDEs G). Here, the trans- formed data are linearly separable.
In general, we have the following theoretical result for the NDDEs, whose proof is provided in Appendix.
5

---

Published as a conference paper at ICLR 2021
OOM.
O
Figure 4: Evolutions of the NDDES (top) and the NODEs (bottom) in the feature space during the training procedure. Here, the evolution of the NODEs is directly produced by the code provided in (Dupont et al. 2019).
Neural ODE Neura DOE
Training loss
0.6
0.4
0.2
0.0
10
20
Epochs (a)
(b)
(c)
(d)
Figure 5: Presented are the training losses (a) of the NODEs and the NDDEs on fitting the function g(z) for d 2. Also presented the flows, from the data at the initial time point (b), of the NODEs (d) and the NDDES (c) after training. The flow of the NODEs is generated by the code provided in (Dupont et al. 2019)
Theorem (Universal upproximating capability of the NDDEs). For any given continuous funciitionnoiio F R. - RR,, if one can construci a neural nelwork for upproximating the map G(æ)) = T [F((x) - x), then there exists an NDDE of n r--imenssoon that can model the map x + F(x), that is, h(T) so F(x) with the initial function (bt)) for < O.
Additionally, NDDEs are suitable for fitting the time series with the delay effect in the original systems, which cannot be easily achieved by using the NODEs. To illustrate this, we use model of DDEs, written as _ A tanh(zz) (x(t) x(t - 7)) with x xt)) = Zo for O. Given the time series generated by the system, we use the NDDEs and the NODES to fit it, respectively. Figure shows that the NNDEs approach much lower loss, compared to the NODES. More interestingly, the NODEs prefer to fit the dimension of 2. and its loss always sustains at larger value, e.g. 0.25 in Figure B The main reason is that two different trajectories generated by the autonomous ODEs cannot intersect with each other in the phase space, due to the uniqueness of the solution of ODEs.
We perform experiments on another two classical DDEs, i.e., the population dynamics and the Mackey-Glass system (Erneux 2009). Specifically, the equation of the dimensionless population dynamics is å = (1 (t )) where a=(t) is the ratio the population to the carrying the growth rate and is the time delay. The Mackey-Glass system
capacity of the population, written as à B - - (tt(),, where =(t) is the number of the blood cells, G, n, T,? are the parameters of biological significance. The NODEs and the NDDEs are tested on these two dy- namies. As shown in Figure 7 very low training loss is achieved for the NDDEs while the the

---

Published as a conference paper at ICLR 2021
Training loss
1.0
Parameters
1
X1.X2
x2
0.5
0.0
0
0.5
0.0
0.5
5000
5000
1.0
Parameters
5
Training loss
X1
0
0.5
-1
-1
0.0
0
-0.5 0.0 0.5
5000
5000
X1
Iterations
Iterations
Figure 6: Comparison of the NDDEs (top) versus the NODEs (bottom) in the fitting 2-D time series with the delay effect the original system. From the left to the right, the true and the fitted time series, the true and the fitted trajectories in the phase spaces, the dynamics of the parameters in the neural networks, and the dynamics of the losses during the training processes.
loss of the NODEs does not go down afterwards, always sustaining at some larger value. As for the predicting durations, the NDDEs thereby achieve better performance than the NODEs. The details of the training could be found in Appendix.
0.3
0.3
0.3
Neural DDE ANDDE
Neural DDE ANODE
Neural DDE ANDDE
Neural DDE ANODE
0.2
Test loss
0.2
Test loss
0.2
Test loss
0.2
Training
0.1
0.1
0.1
0.0
0.0
0.0
0.0
1500
3000
0
1500
3000
1500
3000
1500
3000
0.3
0.3
0.3
0.3
Neural DDE ODE
Neural
Neural
0.2
Test loss loss
0.2
Test loss
0.2
loss
0.2
Training
ANDDE
ANODE
ANODE
ANDDE
0.1
0.1
0.1
0.0
0.0
0.0
0.0
1500
3000
0
1500 Epochs
3000
1500 Epochs
3000
0
1500
3000
Epochs
Epochs
Figure 7: The training losses and the test losses of the population dynamics (top) and the Mackey- Glass system (bottom) by using the NDDES, the NODES and the ANODES (where the augmented dimension equals to 1). The first column in the figures is the training losses. The figures from the second column to the fourth column present the test losses over the time intervals of the length т, 2r, 5T. The delays for the two DDEs are both designed as 1, respectively. The other parameters are set as r = 1.8, B = 4.0, n = 9.65, and 7 2.
4.2 EXPERIMENTS ON IMAGE DATASETS
For the image datasets, we not only use the NODES and the NDDEs but also an extension of the NDDEs, called the NODE+NDDE, which treats the initial function as ODE. In our experiments, such a model exhibits the best performance, revealing strong capabilities in modelling and feature representations. Moreover, inspired by the idea of the augmented NODEs (Dupont et al. 2019), extend the NDDEs and the NODE+NDDE to the A+NDDE and the A+NODE+NDDE, respectively. Precisely, for the augmented models, augment the original image space to higher-dimensional

---

Published as a conference paper at ICLR 2021
space, i..,, Rex -> ((c+))xxx"" where h, w, and p are, respectively, the number of channels, height, width of the image, and the augmented dimension. With such configurations of the same augmented dimension and approximately the same number of the model parameters, comparison studies on the image datasets using different models are reasonable. For the NODES, we model the vector filed f [(((t) as the convolutional architectures together with slightly different hyperparam- eter setups in (Dupont et al. 2019). The initial hidden state is set as fe(0) E Rexhxw with respect to an image. For the NDDES, we design the vector filed as f(concat (f(t), h - T))), mapping the space from IR 2xxxxx hxw to Rexhxx where concat(:,.) executes the concatenation of two tensors on the channel dimension. Its initial function is designed as a constant function, i.e., h(t) is identical to the input/image for 1 < O. For the NODE+NDDE, we model the initial function as an ODE, which follows the same model structure of the NODEs. For the augmented models, the augmented dimen- sion is chosen from the set {1, 2, 4]. Moreover, the training details could be found in Appendix, including the training setting and the number of function evaluations for each model on the image datasets.
The training processes on MNIST, CIFARIO, and SVHN are shown in Fig. g Overall, the NDDEs and its extensions have their training losses decreasing faster than the NODES/ANODEs which achieve lower training and test loss. Also the test accuracies are much higher than the that of the NODES/ANODES (refer Tab. D). Naturally, the better performance of the NDDEs is attributed to the integration of the information not only on the hidden states at the current time / but at the previous time t- F. This kind of framework is akin to the key idea proposed in (Huang al. 2017 1 where the information is processed on many hidden states. Here, we run all the experiments for 5 times independently.
Training loss Training loss Training loss
2.0
1.75
NODE+NDO
0.6
Test loss
Accuracy
1.5
Neural DDE
Neural DDE
1.50
Neural ODE
0.5
NODE NDDE
1.0
Neural DDE
0.5
1.25
0.4
Neural ODE
0
15
30
0
15
30
0
15
30
1.0
1.0
1.0
NODE+NDDE
Test loss
NDDE
Accuracy
Neural DDE
Neural DDE
0.9
0.5
0.5
NODE+NNDE
0.8
Neural DDE
Neural ODE
0.0
0.0
0.7
0
15
30
0
15
30
0
15
30
1.5
1.5
0.9
1.0
Neural DDE
NODE+NDDE
Test loss
Neural DDE
NODE+NDDE
Accuracy
0.8
Naural OnE
1.0
ODE
OODC+NDDE
0.5
0.7
Neural DDE
Neural ODE
0.0
0.5
0.6
0
15 Epochs
30
0
15 Epochs
30
0
15 Epochs
30
Figure& 8: The training loss (left column), the test loss (middle column), and the accuracy (right column) over 5 realizations for the three image sets, i.e., CIFARIO (top row), MNIST (middle row), and SVHN (bottom row).
DISCUSSION
In this section, we present the limitations of the NDDEs and further suggest several potential direc- tions for future works. We add the delay effect to the NDDEs, which renders the model absolutely irreversible. Algorithm [ thus requires storage of the checkpoints of the hidden state h(t) at every instant of the multiple т. Actually, solving the DDEs is transformed to solving the ODEs of an increasingly high dimension with respect to the ratio of the final time T and the time delay т. This definitely indicates a high computational cost. To further apply and improve the framework of the NDDEs, a few potential directions for future works are suggested, including:
8

---

Published as a conference paper at ICLR 2021
CIFARIO
MNIST
SVHN
NODE
53.92% + 0.67
96.21% + 0.66
80.66%
0.56
NDDE
55.69% + 0.39
96.22% 0.55
81.49%
0.09
NODE+NDDE
55.89% +: 0.71
97.26% - 0.22
82.60% +
AI+NIDE
56.14% + 0.48
97.89%
0.14
81.17%
1.229
AL+NDDE
56.83% + 0.60
97.83%
+.07
82.46%
0.28
AI+NODE+NDDE
57.31% + 0.61
0.07 98.16%
83.02% + 0.37
A2+NODE
57.27% + 0.46
98.25%
0.08
81.73%
1..3%-0..0 0.92
A2+NDDE
58.13% + 0.32
98.22%
0.04
82.43%
0.26
A2+NODE+NDDE
58.40% 0.31
98.26% 0.06
83.73%+0.72
A4+NODE
58.93% + 0.33
98.33%
0.12
82.72% 0.60
A4+NDDE
59.35% + 0.48
98.31% 0.03
82.87% 0.55
MANNDE+NNDDE
59.94% + 0.66
98.52% + 0.11
83.62% - 0.51 0.51
The test accuracies with their standard deviations over 5 realizations the three image datasets. In the first column, p (=1, 2, or 4) in Ap means the number of the channels of zeros into the input image during the augmentation of the image space Rexx - (g(ep)) xhx (Dupont al., 2019). For each model, the initial (resp. final) time is set as 0 (resp. 1), and the delays of the NDDES and extensions are all set 1, simply equal the final time.
Table 1:
Applications to more real-world datasets. In the real-world systems such as physical, chemical, biological, and ecological systems, the delay effects are inevitably omnipresent, truly affecting the dynamics of the produced time-series (Bocharov & Rihan 2000 Kajiwara et 20.22. The NDDEs are undoubtedly suitable for realizing model-free and accurate prediction (Quaglino et al. 2019). Additionally, since the NODEs have been generalized to the areas of Computer Vision and Natural Language Processing He et al. (2019); Yang et al.] (2019); Liu et al. (2020). the framework of the NDDEs probably be applied to the analogous areas, where the delay effects should be ubiquitous in those streaming data.
Extension of the NDDEs. A single constant time-delay in the NDDEs be further generalized to the case of multiple or/and distributed time delays (Shampine Thompson (00T). As such, the model is likely have much stronger capability to extract the feature, because the model lever- ages the information at different time points to make the decision in time. All these extensions could be potentially suitable for some complex tasks. However, such complex model may require tremendously huge computational cost.
Time-dependent controllers. From viewpoint of control theory, the parameters in the NODEs/NDDEs could be regarded as time-independent controllers, viz. constant controllers. A natural generalization way is to model the parameters as time-dependent controllers. controllers were proposed in (Zhang et al., 2019b), where the parameters w(t) were treated as an- other learnable ODE i = q(w(t), p). q +) is a different neural network, and the parameters p and the initial state w0) are pending for optimization. Also, the idea of using a neural network to generate the other one was initially conceived some earlier works including the study of the hypernetworks (Ha et al. 2016).
fact, such
CONCLUSION
In this article, we establish the framework of NDDEs, whose vector fields are determined mainly by the states at the previous time. We employ the adjoint sensitivity method to compute the gradients of the loss function. The obtained adjoint dynamics backward follow another DDEs coupled with the forward hidden states. We show that the NDDEs can represent some typical functions that cannot be represented by the original framework of NODEs. Moreover, we have validated analytically that the NDDEs possess universal approximating capability. We also demonstrate the exceptional efficacy of the proposed framework by using the synthetic data real-world datasets. All these reveal that integrating the elements of dynamical systems to the architecture of neural networks could be potentially beneficial the promotion of network performance.

---

Published as a conference paper at ICLR 2021
AUTHOR CONTRIBUTIONS
WL conceived the idea, QZ. YG WL performed the research, QZ & YG performed the mathe- matical arguments, QZ analyzed the data, and QZ & WL wrote the paper.
ACKNOWLEDGMENTS
WL was supported by the National Key RD Program of China (Grant no. 2018YFC0116600) the National Natural Science Foundation of China (Grant no. 11925103), and by the STCSM (Grant No. 18D21201000, 1951 1101444,, 19511 132000).
REFERENCES
Gennadii A Bocharov and Fathalla A Rihan. Numerical modelling in biosciences using delay differ- ential equations. Journal of Computational and Applied Mathematics, 125(1-2):133-19, 2000.
Bo Chang, Lili Meng, Eldad Haber, Frederick Tung, and David Begert. Multi-level residual net- works from dynamical systems view. arXiv preprint arXiv:1710.10335: 2017.
Bo Chang, Minmin Chen, Eldad Haber, and Ed H Chi. Antisymmetricrnn: A dynamical system view recurrent neural networks. arXiv preprint arXiv:1902.09653 2019.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In Advances in neural information processing systems, pp. 6571-6588 2018.
Emilien Dupont, Arnaud Doucet, and Yee Whye Teh. Augmented neural odes. In Advances in Neural Information Processing Systems, pp. 3140-3150, 2019.
Weinan E. A proposal on machine learning via dynamical systems. Communications in Mathematics and Statistics, 5(1):1-11, 2017.
Weinan E. Jiequn Han, and Qianxiao Li. A mean-field optimal control formulation of deep learning. Research in the Mathematical Sciences, 6(1):10, 2019.
Jeffrey L Elman. Finding structure in time. Cognitive science, 14(2):199-211, 1990.
Thomas Emeux Applied delay differential equations, volume 3. Springer Science & Business Media, 2009.
Leheng Fang. Wei Lin, and Qiang Luo. Brain-inspired constructive learning algorithms with evo- lutionally additive nonlinear neurons. International Journal of Bifurcation and Chaos, 28(05): 1850068-, 2018.
David Ha, Andrew Dai, and Quoe V Le. Hypernetworks. arXiv preprint arXiv: 06096699000000 2016.
Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse Problems, 34(1):014000 2017.
Jiequn Han, Arnulf Jentzen, and E Weinan. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115(34):8505-8511 8510, 2018.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp- 770-778, 2016.
Xiangyu He, Zitao Mo. Peisong Wang, Yang Liu, Mingyuan Yang, and Jian Cheng. Ode-inspired network design for single image super-resolution. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 1732-1741, 2019.
Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger. Densely connected convolutional networks. Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 4700-4708, 2017.
10

---

Published as a conference paper at ICLR 2021
Tsuyoshi Kajiwara, Toru Sasaki, and Yasuhiro Takeuchi. Construction of Iyapunov functionals for delay differential equations in virology and epidemiology. Nonlinear Analysis: Real World Applications, (000)1802-18200 2012.
Zhenzhong Lan, Mingda Chen, Sebastian Goodman, Kevin Gimpel, Piyush Sharma, and Radu Sori- cut. Albert: A lite hert for self-supervised learning of language representations. arXiv preprint arXiv:1909. 19942, 2019.
Qianxiao Li and Shuji Hao. An optimal control approach deep learning and applications to discrete-weight neural networks. arXiv preprint arXiv:1803.012: 2018.
Qianxiao Li, Long Chen. Cheng Tai, and E Weinan. Maximum principle based algorithms for deep learning. The Journal of Machine Learning Research, 18(1):5998-66) -6266, 2017.
Daniel Liberzon. Calculus of variations and optimal control theory: concise introduction. Prince- ton University Press, 2011.
Xuanqing Liu, Tesi Xiao, Si Si, Qin Cao, Sanjiv Kumar, and Cho-Jui Hsich. Neural sde: Stabilizing neural ode networks with stochastic noise. arXiv preprint arXiv:1906002355 2019. 2019.
Xuanqing Liu, Hsiang-Fu Fu Yu, Inderjit Dhillon, and Cho-Jui Hsich. Learning to encode position for transformer with continuous dynamical model. arXiv preprint arXiv:2003.09229, 2020.
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. Pde-net: Learning pdes from data. In International Conference Machine Learning, pp. 3208-32011 2018.
Zichao Long, Yiping Lu, and Bin Dong. Pde-net 2.0: Learning pdes from data with a numeric- symbolic hybrid deep network. Journal of Computational Physies, 399: 399:108925, 2019.
Yiping Lu, Aoxiao Zhong, Quanzheng Li, and Bin Dong. Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations. International Conference on Machine Learning. pp. 3276-3285. PMLR, 2018.
Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott. Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach. Physical Review Letters, 120(2), 2018
LS Pontryagin, VG Boltyanskij, RV Gamkrelidze, and EF Mishchenko. The mathematical theory of optimal processes. 1962.
Alessio Quaglino, Marco Gallieri, Jonathan Masci, and Jan Koutník. Snode: Spectral discretization of neural odes for system identification. arXiv preprint arXiv:1906.07033. 2019.
David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams. Learning representations by back- propagating errors. nature, 323(6888) 1986.
Lars Ruthotto and Eldad Haber. Deep neural networks motivated by partial differential equations. Journal of Mathematical Imaging and Vision, pp. 1-13, 2019.
Lawrence F Shampine and Skip Thompson. Solving ddes in matlab, Applied Numerical Mathemat- ics, 37(4):441-45 4588, 2001.
Qi Sun, Yunzhe Tao, and Qiang Du. Stochastic training of residual networks: a differential equation viewpoint. arXiv preprini
2018
Yifan Sun, Linan Zhang, and Hayden Schaeffer. Neupde: Neural network based ordinary and partial differential equations for modeling time-dependent data. In Mathematical and Scientific Machine Learning, pp. 352-372. PMLR, 2020.
Yang Tang, Jürgen Kurths, Wei Lin, Edward Ott, and Ljupco Kocarev. Introduction to focus is- sue: When machine learning meets complex systems: Networks, chaos, and nonlinear dynamics. Chaos: An Interdisciplinary Journal of Nonlinear Science, 30315, 2020.

---

Published as a conference paper at ICLR 2021
Guandao Yang, Xun Huang, Zekun Hao, Ming- Yu Liu, Serge Belongie, and Bharath Hariharan. Pointflow: 3d point cloud generation with continuous normalizing flows. In Proceedings of the IEEE International Conference Computer Vision, pp. 4541-45551 2019.
Dinghuai Zhang, Tianyuan Zhang, Yiping Lu, Zhanxing Zhu, and Bin Dong. You only propagate once: Painless adversarial training using maximal principle. arXiv preprint arxh::1955.0087. (3), 2019a.
Tianjun Zhang, Zhewei Yao, Amir Gholami, Joseph E Gonzalez, Kurt Keutzer, Michael W Ma- honey, and George Biros. Anodev2: A coupled neural ode framework. In Advances in Neural Information Processing Systems, pp. 5151-5161, 2019b.
Qunxi Zhu, Huanfei Ma, and Wei Lin. Detecting unstable periodic orbits based only on time series: When adaptive delayed feedback control meets reservoir computing. Chaos, 29(9):93125-93122 2019.
A APPENDIX
A.1 THE FIRST PROOF THEOREM !
Here, we present direct proof of the adjoint method for the NDDEs. For the neat of the proof, the following notations are slightly different from those in the main text.
Let xt)) obey the DDE written
(t) (()(t),((t),)))) ((t), u(t) = x(t -)), r), E 10, xx())= zo, E [-r,0].
whose adjoint state is defined
aL
x(1) == Ox(()'
(8)
where L is the loss function, that is, : ((X())), After discretizing the above DDE, we have
xxt + At) =x(t) + At ((),
=x()) + At f(=(t), o))
x + At) = =x(t + T) + At ((((((((++)))) ytt +r), 0(t +)),
=x(l + -) + Al =())) -))-
According to the definition of A() and applying the chain rule, have
OL
Oxi( + At)
aL
Ex(t + + At)
A(t) ōr(t At)
(1)
(x(()
Oz(t + At)
θx(t)
dx + At)
+ + + (t(
X(t + At)
dx(()
dr(()
(10)
-(t+ At) + At fc(=(t))((t)
+ +X(++ At) (a((t + T), 0), +))) X(o.T=((())
= A(E At) + At [A(e + At) fa(t) + X(t + At) [t
(t(++ +T) x10))=-((())
which implies,
x(t) lim X(t + At) - A(1)
At -t(( fa(t) - X(t +T) ft +)) 7) x(o,-=-((t)
(11)
where (-) derived as
characteristic function. For the parameter (t)), the result can analogously
aL
aL dx(t + At)
00(t)
dzit + At) 00(t)
(12)
At X(t + At) (t)
12

---

Published as a conference paper at ICLR 2021
In this article, (tt) is considered to be constant variable, i.e., 8(t) 8, which yields:
al
lim E At A(t At) ft(
00
(() fo(t)dt
(13)
-X(t) fo(t)dt.
To summarize, we get the gradients with respect to z(0) and 0 in form of augmented DDEs. These DDEs are backward in time and written in an integral form of
=(0) x(T)+ { ((z.y.0)dt.
0L
OL 0
+ fr -X(t) f2(t) +(+ -f(t (t))
(14)
(x(
(x(T) 0
0L
A(t) fo(t)dt.
00
T
A.2 THE SECOND PROOF OF THEOREM I
Here, we mainly employ the Lagrangian multiplier method and the calculus of variation to derive the adjoint method for the NDDEs. First, we define the Lagrangian by
" A(0(==((x(()))))) d))
C = L((T() +
(15)
We thereby need to find the so-called Karush-Kuln -uhh--ukkerrrrr Tucker (KKT) conditions for the Lagrangian, which is necessary for finding the optimal solution of 9. To obtain the KKT conditions, the calculus of variation is applied by taking variations with respect to A(t) and (t)).
For A()), let Â(t) be a continuous and differentiable function with a scalar e. We add the perturbation (()) to the Lagrangian L, which results in a new Lagrangian, denoted by
((e) = L(x(T)) + (. (A(1) + A(()) - [(x(0, 0).
(16)
In order to obey the optimal condition for A(t), we require the following equation
d.(cc)
de f
x(t)( - f(x()) d)
(17)
be zero. Due the arbitrariness of A(t), we obtain
)-f(c(t). 0, Vt e [0,T],
(18)
which is exactly the DDE forward in time.
Analogously, we can take variation with respect to xxf) and let (tt) be a continuous and differen- tiable function with scalar e. Here, .itt) for t E [-7, 0]. We also denote by ((t) Ø(t 7)
fort E [0, T]. The new Lagrangian under the perturbation of &(t) becomes
C(c) := L(x(T) Cf((T) + l'' A(t) ddx(t) + cb(t)
f(x(t) + cz((). v(t) + co(t), 0) dt.
dt
(19)
13

---

Published as a conference paper at ICLR 2021
We then compute the dCfe) de which gives
dc(c)
8(T) + L' A(t) ØL
dā(() - (
(z(())), o)()))
-=0 = Ox(T)
de
dt
OL
T(T) ((1)@(1)) !'" -5(1) dd(t) (())(((())00)) ()²(1) X(1)((((11) v()), () - A()) ((T) + {-¹ -8(1) dA(t) X(O_((())
integration by parts
dx
dl
+
ØL
(T(TT)
di
L (())s v()), r)dt
(x(T) + +X(T) a(T)+ 1-' -a(1) dd(t) (0))(((()))000)00()0)
OL
dl
A(t' TT)fy(a(t + +T),u(t' + r), (t)))((t))dt
ØL
X(T) B(T)+ 1'' -a(t) dx(l) (0))
(x(T)
dt
L' (at' +T), ()) ()))))
variation on the interval [-r,0] 0]
OL
+ +X()) B(T)
(T((T)
11 t(t)
(-aMt) dx(t)
(-d di --,-::/:0000000-----------------
de
(20)
Notice that dC(c)
c= 0 is satisfied for all continuous differrttaleeeeeeee i(f) at the optimal x0t). Thus,
we have
dt(t)
aL
dt - + +))
A(T) = 0x
(21)
Therefore, the adjoint state follows a DDE as well.
A.3 THE PROOF OF THEOREM 2
The validation for Theorem is straightforward. Consider the NDDEs in the following form:
{ dhe fh: w)), (((t) x,
t>= t t<=0 0
where т equals to the final time T. Due to h(t) = x for t S 0. then the vector field of the NDDES in the interval (0, T is constant. This implies h(T)
T f(x; w). Assume that the neural network f(x; w) is able to approximate the map G(x:) = [F(x) - xr). Then, we have
h(T) x + [F(()) - x = F(x). The proof is completed.
A.4 COMPLEXITY ANALYSIS OF THE ALGORITHM I
As shown in (Chen et al. 2018), the memory and the time costs for solving the NODEs are O(1) and OLL) where L is the number of the function evaluations in the time interval [0. 77. More precisely, solving the NODEs is memory efficient without storing any intermediate states of the evaluated time points in solving the NODEs. Here, intend to analyze the complexity of Algorithm T It should be noting that the fftee-ofthhe-ee for the DDE software not as advanced as that for the ODE software. Hence, solving DDE is much difficult compared with solving the ODE. There exists several DDE
14

---

Published as a conference paper at ICLR 2021
solvers, such the popular solver, the dde23 provided by MATLAB. However, these DDE solvers usually need to store the history states to help the solvers access the past time state h(t - -))- Hence, the memory cost of DDE solvers is (H)), where H is the number of the history states. This is the major difference between the DDEs and the ODEs, as solving the ODEs memory efficient, i.e., O(1). In Algorithm 1 propose piece-wise ODE solver to solve the DDEs. The underlying idea is to transform the DDEs into the piece-wise ODEs for every T time periods such that one can naturally employ the framework of the NODES. More precisely, we compute the state at time kt by solving an augmented ODE with the augmented initial state, i.e., concatenating the states at time 1)T into single vector. Such a method has several strengths and has weaknesses as well. The strengths include:
T, 0, T, ..... (k
One can easily implement the algorithm using the framework of the NODEs, and
becomes quite memory efficient (n). where only need store small number of the forward states, h(0).... (o))....ii h..r), to help the algorithm compute the adjoint and the gradients in reverse mode. Here, the final T is assumed to be not very large compared with the time delay т, for example, = nr with a small integer n.
Algorithm [
Notably, we chosen n = (i.e., = r = 1.0) of the NDDEs in the experiments on the image datasets, resulting quite memory efficient computation. Additionally, the weaknesses involve:
Algorithm[mmay suffer from high memory cost, if the final T is extremely larger than the T (i.e., many forward states, (0),, h(nr), are required to be stored), and
the increasing augmented dimension of the ODEs may lead to the heavy time cost for each evaluation in the ODE solver.
B EXPERIMENTAL DETAILS
B. SYNTHETIC DATASETS
For the classification of the dataset of concentric spheres, we almost follow the experimental con- figurations used in Dupont et al. 2019). The dataset contain 2000 data points in the outer annulus and 1000 data points in the inner sphere. We choose 11 = 0.5, r2 1.0, and r3 = 1.5. We solve the classification problem using the NODE and the NDDE, whose structures are described, respectively, as:
1. Structure of the NODEs: the vector field is modeled as
f(z) Wout Re...u.WWWW.................
where Wi,, 32x2 W I 32 and W Wout 2x32;
2. Structure of the NDDEs: the vector field is modeled as
f(x(1 - T)) = W ReLL x(( -))),
where W.,. 33222 W E 32x233 and Wout 2R323
For each model, we choose the optimizer Adam with le-3 as the learning rate, 64 as the batch size, and as the number the training epochs.
We also utilize the synthetic datasets to address the regression problem of the time series. The optimizer for these datasets is chosen as the Adam with 0.01 as the learning rate and the mean absolute error (MAE) the loss function. In this article, we choose the true time delays of the underlying systems for the NDDEs. We describe the structures of the models for the synthetic datasets follows.
For the DDE = A tanh(2(t)) + z(t - T)), we generate trajectory by choosing its parameters 0.5, A [[-1, 1]. -1],, xxt)
= [0,1] for t < 0. the final time T 2.5, and the sampling time period equal to 0.1. We model using the NODE and the NDDE, whose structures illustrated, respectively, as:
15

---

Published as a conference paper at ICLR 2021
1. Structure of the NODE: the vector field is modeled as
f() = Wout tanh(Wint
where in E R and Wout E R2x10
2. Structure of the NDDE: the vector field modeled as f(x(t), - 7)) Wont (n + z(t (t 7))), where Wi 10x2 and Wout R2
We use the whole trajectory to train the models with the total iterations equal 5000.
For the population dynamics á rx(()(1 - ar( - r1) with sxt) Io for t S 0 and the Mackey- Glass systems
3 11+r"(t- Tg) >x(() with xxt) 2g for S O, we choose the parameters as T1 72 = 1, = 1.8, 3 = 4. 7 = 9.65, and y = 2. We then generate 100 different trajectories within the time interval [0, 8] with 0.05 as the sampling time period for both the population dynamics and the Mackey-Glan systems. We split the trajectories into two parts. The first part within the time interval [0, 3] is used the training data. The other part is used for testing. We model both systems by applying the NODE, the NDDE, and the ANODE, whose structures are introduced, respectively, as:
x(t-T)
1. Structure of the NODE: the vector field is modeled as f(x) out tanh(W tanh(W..:))),
where Win R0x1, W E R 10x 10 and Wout IR 10
2. Structure of the NDDE: the vector field is modeled as
f(x(t), x(t 7)) = Wont tanh(W W Sanb(W,..onnca(((((..dll -T))))
where Win E 10 W E TR 10x10 and Wout E 1=
3. Structure of the ANODE: the vector field is modeled as [ (1)) Wont lanh(W tanh(W (nnn())))
where Win in IR 10x2 W R 10 Wout 22 and the augmented dimension equal to 1. Notably, we need to align the augmented trajectories with the target trajectories to regressed. To do it, we choose the data in the first component and exclude the data in the augmented component, i.e., simply projecting the augmented data into the space of the target data.
We train each model for the total 3000 iterations.
B.2 IMAGE DATASETS
The image experiments shown in this article are mainly depend on the code provided (Dupont et al. 2019). Also, the structures of each model almost follow from the structures proposed in (Dupont et al., 2201).. More precisely, we apply the convolutional block with the following strutures and dimensions
conv, filters, paddings,
ReLU activation function,
3 X 3 conv, filters, 1 paddings,
ReLU activation function,
1 X conv, filters, paddings.
where different for each model and e is the number of the channels of the input images. In Tab. 2
we specify the information for each model. As can be seen, we fix approximately the same number of the parameters for each model. We select the hyperparamters for the optimizer Adam through 1e-3 as the learning rate, 256 the batch size, and 30 as the number of the training epochs.
C ADDITIONAL RESULTS FOR THE EXPERIMENTS
16

---

Published as a conference paper at ICLR 2021
CIFARIO
MNIST
SVHN
NODE
92/107645
92/84395
92/107645
NDDE
92/107921
92/84487
92/107921
NODEANDDE
64/105680
64/82156
64/105680
AT+NODE
86/108398
87/84335
86/108398
AI+NDDE
85/107189
86/ 86/82944
85/107189
AI+NODE+NDDE
60/107218
61/83526
60/107218
A2+NODE
79/ 19/108332
81/83280
79/ 108332
A2+NDDE
78/107297
81/83473
78/10729
AZ+NODEANDDE
55/ 55/107265
57/83101
55/107265
A4+NODE
63/108426
70/84155
63/108426
A4+NDDE
62/107719
69/83237
62/107719
A4NNODE-NDDE 43/106663
49/ 49/83859
43/ 43/106663
Table 2: The number of the filters and the whole parameters in each model used for CIFARIO 10, MNIST, and SVHN. In the first column, Ap with = 2, 4 indicates the augmentation of the image space gexhxx + g(ctp)xhxw
NODE
ANODE
NDDE
16-
14
16-
Training
14-
12-
12-
10
10
10
56
06
0.95
130
1000
330
100
125
1.5
2.0
1.75
175
3.75
1.90
150
1.50
Test
122
0.50
0.25
025
0.25
0.5
10
1.5
20
03
1.0
155
20
0.
10
15
2.0
2(t 7)
x
Figure 9: The phase portraits of the population dynamics in the training and the test stages. We only show the 10 phase portraits from the total 100 phase portraits for the training and the testing time-series. Here, the solid lines correspond to the true dynamics, while the dash lines correspond to the trajectories generated by the trained models in the training duration and in the test duration, respectively.
17

---

Published as a conference paper at ICLR 2021
NODE
ANODE
NDDE
14.
14-
14-
12.
12.
12.
Training
1.0
10.
13-
0.8
3.5
03.
8..
06.
06-
0.4
44
44
04
06
10
0.4
08
2.8
10
14
04
16
12
24
12-
12-
12-
10
Test
2.8
86.
36
05
54
06
08
1.0
2.6 12
xt T)
z1t
T)
The phase portraits of the Makey-Glass systems exhibiting chaotic dynamics in the training and the stages. We only present the 10 phase portraits from the total 100 phase portraits for the training and the testing time-series. Here, the solid lines correspond to the true dynamics. while the dash lines correspond the trajectories generated by the trained models in the training duration and in the test duration, respectively.
Figure 10:
18

---

Published as a conference paper at ICLR 2021
CIFAR10
MNIST
SVHN
200
200
200
Function Evaluations
NODE+NDDE
NODE+NDDE
NODE-NDDE
NDDE
NODE
NODE
150
150
150
NODE
100
100
100
50
50
50
0
10
20
30
10
20
30
200
200
200
Function Evaluations
NL=NODE=NDDE
AL-NODE-NODE
NDDDD
NDDE
A1 NDDE
Al NODE
150
-NODE
150
150
100
100
100
50
50
0
10
20
30
10
20
30
200
200
200
# Function Evaluations
NDDE
-NDDE
NODE NDDE
A2-NDDE
A2+NDDEE
A2 -DDDO
150
150
NODE
150
100
100
100
50
50
50
0
10
20
30
10
20
30
Function Evaluations
200
200
200
-4-D NDDE
NDDE
-ODEE
ADDNDOE
A4 +DDDE
A4 -DDE
150
A1-NODE
150
A4 MODN
150
A1-NODE
100
100
100
50
50
10
20
Epachs
Epochs
Epochs
Figure 11: Evolution of the forward number of the function evaluations (NFEs) during the training process for each model on the image datesets, CIFARIO, MNIST, and SVHN.
19

---

Published as a conference paper at ICLR 2021
CIFAR10
MNIST
SVHN
# of Function Eva luaaiionns
NODE+NDDE
NODE
NDDE
200
200
150
150
100
100
50
50
NOCE
NODE
0
0
20
200
200
of Function Evalustions
AL-NODE-NDDE AI-NDDE
200
150
A1-NODE
150
150
NLDNODE+NODEE
AI-NODE-NDDE
100
AL-NDDE
100
A1+NDDE
50
50
0
10
20
30
10
20
10
20
30
200
200
# of Function Evaluations
A2-NODE IDDE
AD-NDDE
200
150
A2+NODE
150
150
100
A2+NDDE
100
A2+NDDE
100
NODE
A2-NOCE
50
50
50
0
10
20
30
0
1.0
20
30
10
20
30
# of Function Evaluations
NNDEE IDDE
AR-NDDE
200
AR+NODE+NDOE NDDE
A4 MONODE-NDOEE
100
100
A4+NOOE
N4+NODE
50
Epochs
Epochs
Epochs
Figure 12: Evolution of the backward number of the function evaluations (NFEs) during the training process for each model on the image datesets, CIFAR1O, MNIST, and SVHN.
20


Reviews:
[{"review_id": "zwmuAaW4VET", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "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": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper builds on the Neural ODE framework by using delay differential equations (DDEs) instead of ordinary differential equations (ODEs). This can help in the modeling of systems with a time delay effect, and overcome many limitations of the Neural ODE framework.\n\nUsing DDEs is a novel technique in machine learning and can complement and build on the framework of Neural ODEs and help model systems with time delay dependencies and overcome some limitations of ODEs.\n\nThe paper was well written with a clear description of the model and theory to support it as well as the algorithm to train it. The experiments are well described though they lacked the model parameters and more description on their significance.\n\nMajor points:\n\n1) Model parameters for the experiments were not listed. This can be especially important for systems with a time delay effect as in Figure 6 and 7. Was the same delay used as in the equation of each system? What happens when the delay is different? \n\n2) While some complex systems were used to demonstrate the potential of the proposed model, as in Figure 7, the abstract claims that the model can successfully model chaotic dynamics yet that was not shown in the experiments. The dynamic regime for these systems with their chosen parameters were not clear together with the effect of time delay on these dynamics. Exploring this more closely can demonstrate the full potential of NDDEs.\n\n3) The authors claim that unlike NODEs, NDDEs can overcome the problem of mutually intersected trajectories in phase space. They point to the experiment in Figure 6 as evidence of that. I feel there needs to be more discussion on this point both theoretically, proving that DDEs can indeed model intersections in phase space, and with a clearer demonstration experimentally.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "w7K9Ih8PDvT", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper contributes to a recent line of researches linked to Neural Ordinary Differential Equations (Chen et al. 2018) and variants. \nThis new family of deep neural networks (NODE) generalize ideas from Residual Networks and consider continuous dynamics of hidden units using an Ordinary Differential Equation (ODE) specified by a neural network. Computing in such networks consists in taking x=h(0) as input and define the output layer h(T) to be the solution of an ODE initial value problem at time T. The main ingredient which makes learning possible is the backpropagation algorithm through the last layer ODE solver, that relies on the adjoint sensitivity method (Pontryagin et al., 1962).\n\nThis work considers Delay Differential Equations instead of ODE, which allows to implement more complex dynamics and thus achieve estimation of more complex functions. The contributions of the paper are the following:\n(1)\ta novel deep continuous -time deep neural network defined by the following equations\nFor a given delay \\tau (never discussed), Neural Delay Differential equations define dynamics of the form:\ndh_t\\dt = f(h_t,h_{t-\\tau}, t ;w),  t \\geq 0,\nh(t) = \\phi(t), t \\leq 0\n        In this way the network can take into account a former hidden layer.\n(2)\tthe derivation of the adjoint sensitivity method for delayed equations, which seems relatively straightforward and which is backed by two proofs in the appendix.\n(3)\tA novel learning algorithm that implements the forward for h and the backward pass for h,lambda (the augmented variable) and dL/dw (L is the loss) by a piece-wise ODE solver, dealing with the different delayed states.\n(4)\tExperiments on 2D toydatasets and on classic differential equations such as Mackey-Glass are shown to exhibit the ability of DDE to cope with those dynamics, in constrast to ODE.\n(5)\tExperiments on image datasets\n(6)\t\nThe strong point of this paper is of course the proposal of the new variant of NODE, which comes with a novel algorithm and overcomes some limitations of NODE. I was interested by the examples of functions not covered by NODE and covered by NDDEand easily convinced by that. \n\nThere are some weak points in this work. Some of them could be easily improved I think, others call for further  work.\nRelatively Minor points\n-\tThe paper is not self-content and does not a very good job in explaining the context of Neural ODE. I suggest that you more clearly describe NODE, as a hypothesis class and then as a learning algorithm.\n-\tThe two uses of NDDE concern modeling time-series  or implementing a classification/regression function. Can you each time precise inputs/outputs, samples you use in training. Just for clarity, the reader guesses of course but its rather slows down the reading this absence of notation and formalization ?\nMajor points and questions to the authors:\n-\tMore importantly, as a novel algorithm is introduced, I expect to see a complexity in time analysis. As for NODE I understand that the complexity in space is favorable. \n-\tAn associated question is the role of tau, the delay. How is it chosen ? I imagine that if tau converges towards 0 we find the behaviour of NODE ? what was its value in the two real-world experiments? Was it selected by cross-validation ?  Is it too computationally heavy to consider multiple delays ? \n-\tEventually, I have doubts and questions about the real opportunities for using NDDE in the two real word datasets: yes the divergence of NODE poses problems : in these experiments no doubt that NDDE has a more stable empirical behavior, with a small variance. However, in fine, the average performance is nearly of the same order, a significant difference only on MNIST but with a very bad score considering this is an easy problem. I won’t qualify these results as exceptional and I kindly engage the authors to remain modest.\n-\tNow can we cope without delays by introducing new states in the way augmented  Neural ODE are built (Dupont et al.) ?\n-\tFor what kinds of classification problems NDDE could be more interesting than NODE – The question is certainly difficult.\n\nIn conclusion, the paper tackles a very interesting topic and I had pleasure to discover it (with the help of literature). Although the idea to consider DDE instead of ODE is incremental, its relevant and novel, and certainly promising and worth to be explored because it addresses some of the issues of NODE.\nHowever I stay on my hunger on different points, there are pending questions that require to be answered before acceptance.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "E96gRb_Z8A0", "reviewer": "AnonReviewer2", "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": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper presents a modeling class of parameterized first-order\ndifferential equations that are conditional on some delayed\npast states. This is a promising direction for the community to go\nto push past current ODE modeling limitations, but I recommend for\nrejection in the current form due to improper characterization\nand evaluation with respect to prior work (more details below).\n\n# Strengths\nTo the best of my knowledge this is a novel modeling class\nfor neural ODEs. Delay ODEs are well-studied in contexts\noutside of the machine learning community that would be\nuseful to bring into it.\n\nOn the representation capacity, delay ODEs are able to\novercome one instance of an intersecting path issue that\nthe usual first-order neural ODEs are unable to represent.\n\nThe application to population dynamics and Mackey-Glass systems\nthat are naturally represented an delay differential equations\nseems reasonable and well-motivated that baseline methods may\nnot be able to capture this modeling class.\n\n# Weaknesses\nThe strongest weaknesses of this paper are the characterization\nand evaluation in comparison to [Dupont 2019], which points\nout some of the same modeling issues and proposes an\nalternative way to fix them.\nEven though this paper cites [Dupont2019], it omits that they consider\nidentical experimental settings and does not include their results:\n\n  1. Figure 4 and 5 are almost identical to figures in [Dupont2019]\n\n  2. Table 1 presents results in MNIST and CIFAR10 and is identical\n     to the setting in [Dupont2019]. The neural ODE baselines\n     are almost consistent with the ones reported in [Dupont2019], but\n     [Dupont2019] outperforms the method presented here, and are\n     omitted from the table. [Dupont2019] additionally\n     considers the SVHN dataset.\n\n# Other comments and questions\nIs the adjoint method for the NDDEs very different than the adjoint\nmethod for NODEs? It could be insightful to intuitively explain\nthis.\n\nThe extension to more general delay differential equations seems\npromising and working out the details in this general case in\naddition to the case of having a single delay would add to the paper's\ntheory and methodological contributions.", "extra_scores": {}, "extra_text": {}}, {"review_id": "0SpVZ9DAW2v", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "-- The motivation for resorting to DDEs over ODEs is laid out reasonably well, although I think the intuition behind why NODEs cannot learn certain functions is lacking (specifically, that trajectories cannot cross due to the lack of \"memory\" provided by the delay).\n\n-- Figure 2 could be more informative. As it stands, it only indicates that the initial state of a NDDE is specified by a function over an interval, rather than a single vector, which is all that is being used.\n\n-- Algorithm 1 is challenging to understand without some knowledge of how DDEs are solved.\n\n-- In the experiments, NDDEs are only compared against NODEs, but it seems like ANODEs are the real competition. Notably in the image experiments performed in the ANODE paper they report substantially higher accuracies than NDDEs.\n\n-- This is mentioned in the discussion section, but I think the paper would benefit from an experiment involving a real-world dataset with some sort of delay component, ideally in which NDDEs outperform NODEs/ANODEs. The existing synthetic experiments are designed to show the strengths of NDDEs, while the image classification tasks don't substantially differentiate them from NODEs.\n\n-- The number of function evaluations is not reported anywhere, which seems like a significant statistic when NDDEs are substantially more complex than NODEs.\n\n-- Figure 7 could probably be improved by reducing the number of different parameter values reported (e.g. 3 columns instead of 6). The description also lists 5 parameter values but there are 6 columns.\n\n-- The discussion paragraph \"Extensions of the NDDEs\" suggests a generalization to more complex delays, but all of the models in this paper ensure that the delayed value is always the constant initial value. Also, modeling the initial function as an ODE is mentioned, and I wonder whether the implication is that this would involve learning the initial function as a NODE itself. Either way, it's not obvious how either of these would yield better representations.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper is a variant of the large growing class of Neural ODEs, and adds dependency on a time delay to the baseline, which allows to model a larger class of physical systems, in particular adding the possibility of crossing paths in phase space.

After initial evaluation, the paper was on the fence, with 2 reviewers providing favorable reviews, and 2 reviewers recommending rejection. A particular important issue raised was positioning with respect to prior art, [Dupont 2019], with some substantial overlap between the papers; requests of theoretical discussions of the class of studied systems and its properties.

Most of these remarks have been addressed by the authors, in particular positioning and experimental comparisons.

The AC judged that the paper had been sufficiently improved and recommends acceptance.
Decision: Accept (Poster)