SNODE: Spectral Discretization of Neural ODEs for System Identification
Authors: Alessio Quaglino
Marco Gallieri
Jonathan Masci
Jan Koutník
Venue:  (2020)
Keywords: Recurrent neural networks
system identification
neural ODEs

Abstract: This paper proposes the use of spectral element methods \citep{canuto_spectral_1988} for fast and accurate training of Neural Ordinary Differential Equations (ODE-Nets; \citealp{Chen2018NeuralOD}) for system identification. This is achieved by expressing their dynamics as a truncated series of Legendre polynomials. The series coefficients, as well as the network weights, are computed by minimizing the weighted sum of the loss function and the violation of the ODE-Net dynamics. The problem is solved by coordinate descent that alternately minimizes, with respect to the coefficients and the weights, two unconstrained sub-problems using standard backpropagation and gradient methods. The resulting optimization scheme is fully time-parallel and results in a low memory footprint. Experimental comparison to standard methods, such as backpropagation through explicit solvers and the adjoint technique \citep{Chen2018NeuralOD}, on training surrogate models of small and medium-scale dynamical systems shows that it is at least one order of magnitude faster at reaching a comparable value of the loss function. The corresponding testing MSE is one order of magnitude smaller as well, suggesting generalization capabilities increase.

Published as a conference paper at ICLR 2020
SNODE: SPECTRAL DISCRETIZATION OF NEURAL ODES FOR SYSTEM IDENTIFICATION
Alessio Quaglino; Marco Gallieri,
Jonathan Masci, Jan Koutník
NNAISENSE, Lugano, Switzerland
(alessio, marco, jonathan, jan @nnaisense.com
ABSTRACT
This paper proposes the use of spectral element methods (Canuto et al. 1988) for fast and accurate training of Neural Ordinary Differential Equations (ODE-Nets; Chen et al. 2018) for system identification. This is achieved by expressing their dynamics as a truncated series of Legendre polynomials. The series coefficients, as well the network weights, are computed by minimizing the weighted sum of the loss function and the violation the ODE-Net dynamics. The problem solved by coordinate descent that alternately minimizes, with respect to the coefficients and the weights, two unconstrained sub-problems using standard backpropagation and gradient methods. The resulting optimization scheme is fully time-paralles and results in a low memory footprint. Experimental comparison to standard methods, such as backpropagation through explicit solvers and the adjoint technique (Chen et al., 2018), on training surrogate models of small and medium-scale dynamical systems shows that is at least one order of magnitude faster at reaching comparable value of the loss function. The corresponding testing MSE is one order of magnitude smaller as well, suggesting generalization capabilities increase.
INTRODUCTION
Neural Ordinary Differential Equations (ODE-Nets: Chen et al. 2018) can learn latent models from observations that are sparse in time. This property has the potential to enhance the performance of neural network predictive models in applications where information is sparse in time and it is important to account for exact arrival times and delays. complex control systems and model- based reinforcement learning, planning over a long horizon is often needed, while high frequency feedback necessary for maintaining stability (Franklin et al. 2014). Discrete-time models, including RNNs (Jain & Medsker 1999), often struggle to fully meet the needs of such applications due to the fixed time resolution. ODt-Nets have been shown to provide superior performance with respect to classic RNNs time series forecasting with sparse training data. However, learning their parameters computationally intensive. In particular, ODE-Nets are memory efficient but time inefficient. this paper, we address this bottleneck and propose novel alternative strategy for system identification.
be
Summary of contributions. We propose SNODE, a compact representation of ODE-Nets for system identification with full state information that makes use of higher-order approximation of its states by means of Legendre polynomials. This is outlined in Section 4 In order to find the optimal polynomial coefficients and network parameters, we develop a novel optimization scheme, which does not require to solve an ODE at each iteration. The resulting algorithm is detailed in Section and is based on backpropagation (Linnainmaal 1970, Werbos 1981 Lecun 1988) and automatic differentiation (Paszke et al. 2017). The proposed method is fully parallel with respect to time and its approximation error reduces exponentially with the Legendre polynomial order (Canuto et al. 1988).
Summary of numerical experiments. In Section 5 our method is tested
6-state vehicle problem, where it at least one order or magnitude faster in each optimizer iteration than explicit and adjoint methods, while convergence is achieved in third of the iterations. At test time, the MSE is reduced by one order magnitude. In Section the method is used for 30-state system consisting
"Corresoonding author.

---

Published as a conference paper at ICLR 2020
of identical vehicles, coupled via known collision avoidance policy. Again, method converges in a third of the iterations required by backpropagation thourgh solver and each iteration 50x faster than the fastest explicit scheme.
NEURAL ORDINARY DIFFERENTIAL EQUATIONS
The minimization of scalar-valued loss function that depends the output of ODE-Net can be formulated as a general constrained optimization problem:
min MER
L(t, x(t)) dt,
(I)
s.t.
i(t) = f(t, xt), u(t);0),
x(to)
z0.
where xzt) X is the state, u(t) E is the input, the loss and ODE functions L and f are given, and the parameters / have to be learned. The spaces X and U are typically Sobolev (e.g. Hilbert) spaces expressing the smoothness of =(t) and u(t) (see Section 8). Equation { can be used to represent several inverse problems, for instance in machine learning, estimation, and optimal control (Stengel)
1994 Law al., 2015 Ross 2009). Problem ( be solved using gradient-basee optimization through several time-stepping schemes for solving the ODE, (Chen et al.) (018 Gholami et al., 2019) have proposed use the adjoint method when is neural network. These methods are typically relying on explicit time-stepping schemes (Butcher & Wanner 1996). Limitations of these approaches are briefly summarized:
Limitations of backpropagation through an ODE solver. The standard approach for solving this problem is to compute the gradicnts 0L/00 using backpropagation through a discrete approximation of the constraints, such Runge-Kutte Kutta methods (Runge 1895 Butcher & Wanner 1996) multi- step solvers (Raissi et al. 2018). This ensures that the solution remains feasible (within a numerical tolerance) at each iteration of gradient descent method. However, it has several drawbacks: 1) the memory cost of storing intermediate quantities during backpropagation can be significant, 2) the application of implicit methods would require solving nonlinear equation each step, 3) the numerical error can significantly affect the solution, and 4) the problem topology can be unsuitable for optimization (Petersen et 2018).
Limitations of adjoint methods. ODE-Nets (Chen et al. 2018) solve (I) using the adjoint method, which consists of simulating dynamical system defined by an appropriate augmented Hamiltonian (Ross 2009). with additional state referred to as the adjoint. In the backward pass the adjoint ODE is solved numerically to provide the gradients of the loss function. This means that intermediate states of the forward pass do not need to be stored. An additional step of the ODE solver is needed for the backward pass. This suffers from few drawbacks. 1) the dynamics of either the hidden state or the adjoint might be unstable, due to the symplectic structure of the underlying Hamiltonian system, referred to as the curse of sensitivity in (Ross) 2009); 2) the procedure requires solving a differential algebraic equation and a boundary value problem which is complex, time consuming, and might not have a solution (Ross & Karpenko 2012).
Limitations of hybrid methods. ANODE (Gholami et al., 2019) splits the problem into time batches, where the adjoint is used, storing in memory only few intermediate states from the forward pass. This allows to improve the robustness and generalization improvement could be obtained using reversible integrators. However, its computational cost is of the same order of the adjoint method and it does not offer further opportunities for parallelization.
the adjoint method. A similar
3 RELAXATION OF THE SUPERVISED LEARNING PROBLEM
Our algorithm is based on two ingredients: i) the discretization of the problem using spectral elements leading SNODE, detailed i in Section[4] and ii) the relaxation of the ODE constraint from (I) enabling efficient training through backpropagation. The latter can be applied directly at the continuous level and significantly reduces the difficulty of the optimization, as shown in examples.
2

---

Published as a conference paper at ICLR 2020
The problem in ( is split into two smaller subproblems: one finds the trajectory z(t) that minimizes an unconstrained relaxation of [T) The other trains the network weights 9 such that the trajectory becomes solution of the ODE. Both are addressed using standard gradient descent and backpropa- gation. In particular, a fixed number of ADAM or SGD steps is performed for each problem in an alternate fashion, until convergence. In the following, the details of each subproblem discussed.
Step 0: Initial trajectory.
The initial trajectory x=()) is chosen by solving the problem
f1 L(t, x()) dt,
min x(()X
s.t. x(1o) J0,
If this problem does not have a unique solution, a regularization term is added. For quadratic loss, closed-form solution is readily available. Otherwise, prescribed number of SGD iterations used.
Step 1:
Coordinate descent on residual. Once the initial trajectory x² (t) is found, 0 is computed by solving the unconstrained problem:
/" RR(1, (1),0) di, R(1,x(1), = -100))00)))0000),;;0000000)00 BER"
min
the value the residual at the optimum 8" is smaller than a prescribed tolerance, then the algorithms stops. Otherwise, steps and are iterated until convergence.
Step 2: Coordinate descent on relaxation. Once the candidate parameters 0" are found, the trajectory updated by minimizing the relaxed objective:
min (()EX 10
LL(t, + Rit, x(), dt,
(4)
s.t. x(to) = zn-
Discussion. The proposed algorithm can be seen as alternating coordinate gradient descent on the relaxed functional used in problem (4). i.e., by alternating a minimization with respect to x()) and I O, multiple minima can exist, since each choice of the parameters 0 would induce a different dynamics rtt), solution of the original constraint. For 7 * 0, the loss function in (4) trades-off off the ODE solution residual for the data fitting, providing unique solution. The choice of y implicitly
introduces satisfaction tolerance <(r), i.e., similar to regularized regression (Hastie et al. 2001
implying that (it(t) ftt, x(t;: 0)l S c(7). Concurrently, problem (3) reduces residual.
SNODE - HIGH-ORDER DISCRETIZATION OF THE RELAXED PROBLEM
In order to numerically solve the problems presented in the previous section, discretization of ar()) needed. Rather than updating the values at time points t, from the past to the future, we introduce a compact representation of the complete discrete trajectory by means of the spectral element method.
Spectral approximation. We start by representing the scalar unknown trajectory, x(t), and the known input, u(t), truncated series:
x(t) Σ m.),
u(t) Σ (t)(t).
(5)
where 24, u, R and (t) (t( sets of given basis functions that span the spaces X,, X and U, C U. In this work, we use orthogonal Legendre polynomials of order (Canuto et al. 1988) for (i)(C), where is hyperparameter, and the cosine Fourier basis for (e()), where is fixed.
3

---

Published as a conference paper at ICLR 2020
Collocation and quadrature. In order to compute the coefficients 2; of (5). we enforce the equation at a discrete set Q of collocation points tg. Here, we choose p + Gauss-Lobatto nodes, which include to This directly enforces the initial condition. Other choices are also possible (Canuto et al. [188). Introducing the vectors of series coefficients 2 = {æ: ) and of evaluations quadrature points x((4) (e(t)))yEC: the collocation problem can be solved in matrix form
I1 = M x(to), Mgi == (t,)).
(6)
We approximate the integral as sum of residual evaluations over Q Assuming that z<0) Jo- the integrand at all quadrature points io can be computed as component-wis norm
R(te, 0) = IDM 1DMM-rx(tg) (((2)) f (tg. x(te), utt)): Da: := (tt))))
Fitting the input data. For the case when problem admits unique solution, we propose a new direct training scheme, Ó-SNODE, which is summarized in Algorithm|] [ In general, a least-squares approach must be used instead. This entails computing the integral ( which can be done by evaluating the loss function L at quadrature points tg. If the input data is not available at Eg.. we approximate the integral by evaluating L at the available time points. The corresponding alternating coordinate descent scheme a-SNODE is presented in Algorithm 2 consequences of a low-data scenario on this approach.
In the next sections, we study the
We use fixed numbers N₂ and N₂ of updates for, respectively, 0 and se(t). Both are performed with standard routines, such as SGD. In our experiments, we use ADAM to optimize the parameters and an interpolation order p 14, but any other orders and solvers are possible.
Algorithm 1 6-SNODE training
Algorithm 2 a-SNODE training
Input: M, D from -6-77
Input: M, D. 0 E RT" 2ÿ E ERP,7>0 0. while ([(2]]) + (
x; < arg min, TM.))
while R(Mxj.0 0²) do
for -0... N₂ do
0* <- ADAM, ( [R(Ma).0]
zi + SGD. (xT.(x)) R(Mx,
end while
end for
Output: 8*
for i 0... N. do 0- -- ADAM, [M(M²7,0) 0)]
end for end while Output: 0
Ease of time parallelization. If R(t,) = 0 enforced explicitly at q E Q, then the resulting discrete system can be integrators can only be made parallel across the different components of x(),, the assembly of the residual can be done in parallel also across time. This massively increases the parallelization capabilities of the proposed schemes compared to standard training routines.
implicit time-stepping method of order p. However, while ODE
Memory cost. If an ODE admits a regular solution, with regularity r >> p, in the sense of Hilbert spaces, i.e., of the number of square-integrable derivatives, then the approximation error of the SNODE converges exponentially with p (Canuto et al., 1988). Hence, it produces a very compact representation of an ODE-Net. Thanks to this property, p is typically much lower than the equivalent number of time steps of explicit or implicit schemes with a fixed order. This greatly reduces the complexity and the memory requirements of the proposed method, which can be evaluated at any t via (5) by only storing few x; coefficients.
Stability and vanishing gradients. The forward Fuler method known to have small region of convergence. In other words, integrating very fast dynamics requires a very small time step, dt, in order to provide accurate results. In particular, for the solver error to be bounded, the eigenvalues of the state Jacobian of the ODE need to lie into the circle of the complex plane centered at (-1, 0) with radius 1/dt (Ciccone et al. 2018 Isermann 1989). Higher-ordes explicit methods, such as Runge- Kutta (Runge) 1895), have larger but still limited convergence regions. Our algorithms on the other hand are implicit methods, which have a larger region of convergence than recursive (explicit) methods (Hairer et al. 1993). We claim that this results in more stable and robust training. This
4

---

Published as a conference paper at ICLR 2020
claim is supported by experiments. Reducing the time step improve the Euler accuracy but can still lead to vanishing or exploding gradients (Zilly et al. [2016] Goodfellow et al. 2016). In show that methods do not suffer from this problem.
Appendix C
Experiments setup and hyperparameters. For all experiments, a common setup was employed and no optimization of hyperparameters of 100 were used. Learning rates were set to 10 for ADAM (for all methods) and 10-3 for SGD (for a-SNODE). For the a-SNODE method, ? ADAM algorithms at each epoch, outlined in Algorithm 2 The initial trajectory was perturbed as Ji = za + c. ~ ( 0.1). This perturbation prevents the exact convergence of Algorithm during initialization, allowing to perform the alternating coordinate descent algorithm.
performed. Time horizon 7 = 10s and batch size
and 10 iterations were used for the SGD and
5 MODELING OF PLANAR VEHICLE DYNAMICS
Let us consider the system
ij = J(nu,, Mú + d(o) C(u)u= u,
(8)
where E I are the states, u (Fr O, Tay) is the control, C(v) is the Coriolis matrix, d(v) is the (linear) damping force, and (() encodes the coordinate transformation from the body to the world frame (Fossen 2011). A gray-box model is built using neural network for each matrix
- (m: C(v) fo(u; d(v) fa(u;,
Each network consists of two layers, the first with a tanh activation. Bias is excluded for fc and fa. For fs, so)()) and (os are used as input features, where Φ the vehicle orientation. When inserted in (8) these discrete networks produce ODE-Net that is a surrogate model of the physical system. The trajectories the system and the learning curves are shown in Appendix
6-SNODE. MSE = 0.109
(b) S-SNDE. MSE 0.019
(c) DoPr5. MSE = 0.307
(d) Euler. MSE 0.476
Figure 1: Vehicle dynamics testing in the high-data regime. Each plot compares the true trajectory (dotted red) with the prediction from the surrogate (solid blue). The time horizon 5x than in training. The shown methods are: (a)-(b) the proposed 8-SNODE and a-SNODE schemes, (c) fifth- order Hrrmand-Prince and (d) Explicit Euler. Models trained with our methods have better test MSE and forecast more precisely.
Comparison of methods in the high-data regime. In the case of 8-SNCDE and a-SNODE, only + 1 points are needed for the accurate integration of the loss function, if such points coincide with the Gauss-Lobatto quadrature points. We found that 100 equally-spaced points produce comparable result. Therefore, the training performance of the novel and traditional training methods were compared by sampling the trajectories at 100 equally-spaced time points. Table Talshows that 5-SNODE outperforms EKPR-DoPrS. r5 by factor of 50, while producing significantly improved generalization. The speedup reduces 20 for a-SNODE, which however yields further reduction of the testing MSE by factor of 10, as can be seen in Figure ]
Comparison the low-data regime. The performance of the methods was compared using fewer time points, randomly sampled from uniform distribution. For the baselines, evenly-spaced points

---

Published as a conference paper at ICLR 2020
Table I: Simulation results for the vehicle dynamics example. The standard backpropagation
schemes are denoted by BKPR, while the adjoint is denoted by ADJ. The Formand-Pringe method of order with adaptive time step denoted by DoPr5, 5, while Euler denotes explicit Euler with time step equal to the data sampling time. For 3KPR-DoPr" DoPr5 we set rtol=10-² atol=10 Experiments were performed on a i9 Apple laptop with 32GB RAM. Reported times are per iteration and were averaged over the total number of iterations performed.
(a) High-data regime
(b) Low-data regime
Training
Test
Training
Test
Time [ms]
Lter
Final loss
MSE
% data
# Iter
Final loss
MSE
6-SNODE
480
0.011
0.109
D-SNODE
50
140
0.010
0.029
a-SNODB
134
100
0.011
0.019
a-SNODE
25
190
0.009
0.052
SKFR-Euler
97
1200
0.047
0.476
BKFR-Buler
50
1200
0.163
1,056
RKPR-D6FY5
RKFR-Ful
0.643
3,771
ADJ-suler
208
1200
0.047
0.492
BKFR-DOFr5
50
1160
0.011
0.35
AnJ-DDoPrS
3516
1140
0.011
0.35
RKFR-DoFr5
25
1180
0.011
0.291
were used. Table shows that a-SNODE preserves a good testing MSE, at the price of an increased number of iterations. With only 25% of data, a-SNODE is 10x faster than EKPR-DoPr'S Moreover, its test MSE is 1/7 than BKPR R DoPr5 and up to 1/70 than time step of DoPr3 improves significantly the baseline but is unable match the accuracty of the proposed methods. The adjoint method produced the same results as the backprop (+2%).
Buler, showing that the adaptive
6 LEARNING MULTI-AGENT SIMULATION
Consider the multi-agent system consisting of N₂ kinematic vehicles:
No
1/s = ([n:) tanh
w + Kc(T)) +
Σ Ko(nem No. jif
where 17 E R3N. are the states (planar position and orientation), (1), is the coordinate transform from the body to world frame, common to all agents. The agents velocities are determined by known arbitrary control and collision avoidance policies, respectively, K, and Ko plus some additional high frequency measurable signal = w(t). shared by all vehicles. The control laws are non-linear and are described in detail in Appendix B We wish learn their kinematics matrix by means of a neural network as in Section
3
The task is simpler here, but the resulting ODE has 3NN states, coupled by
Ko. We simulate N₂ 10 agents in series.
(a) 5-SNODE MSE- = 1.7 e-4
(b) A-SNODE MSE = 9.1 e-5
(c) BKPR-Euler MSE = 1.3 e-3
(d) ADJ-Rule er MSE = 5.4
Figure 2: Multi-agent testing high-data regime. Each plot compares the true trajectory (dotted
red) with the prediction from the surrogate (solid blue). The time horizon is 4x longer than in training. The shown methods are: (a)-(b) the proposed 8-SNODE and a-SNODE schemes, and (c)-(d) Explicit Euler. Models trained with our methods have better test MSE and forecast more precisely.
Comparison of methods with full and sparse data. The learning curves for high-data regime are in Figure 3 For method a-SNODE, training terminated when the loss in (4 is less than L R. with L = 0. 11 and R 0.01. For the case of 20% data, we set = 0.01. Table summarizes results. 0-SNODE the fastest method, followed by a-SNODE, which is the best performing. Iteration time of

---

Published as a conference paper at ICLR 2020
BKPR-Euler is 50x slower, with 14x worse test MSE. in between BKPR-Euler and our methods. Random sampligg of the data by 50% and 20% (evenly-spaced for the baselines) makes ADJ-Eulerr fall back the most. BKFR-DoPr" failed to find a time step meeting the tolerances, therefore they were increased to rtol= 10-5, atol= 10-7
is the slowest but its test MSE
Since the loss continued to increase, training terminated at 200 epochs. ADJ-DoPrS failed to compute gradients. Test trajectories are in Figure 2 Additional details are in Appendix B
Table 2: Simulation results for the multi-agent example. Experiments were performed on a 2.9 GHz Intel Core i7 Apple laptop with 16GB of RAM. For the o-SNODE method, each iteration consists of 10 steps of SGD and 10 of ADAM. In high-data regime, a-SNODE converges in 5x less, 2.6x faster, iterations than Euler with 14x smaller test MSE. Method å-SNODE is even 12x faster, with similar performance. Gains are comparable in low-data regime. Method BKPR - failed.
(a) High-data regime
(b) Low-data regime
Training
Training
Time [ms]
leer
Loss
MSB
data
lter
Loss
MSE
106
500
0.12
1.7
0-SNOD2
160
0.22
1.5
o- S-SNODE
1999
180
0.12
9.1
a SNODD
20
175
0,23
e.1
EXPA Buler
5330
1200
0.177
1-3
EXPR Bules
50
990
0.43
28
EXPA-DOPE5
9580
Fail 200
21707
EaPK-Buler
20
990
0.645
1.6.3
ADJ-Buler
8610
990
0.104
5.4
ADJ-Bul or
50
2990
0.53
94c33
ADJ-DeFr5
ADJ-Buler
0.63
Robustness of the methods. The use of a high order variable-step method (DoPr5), providing an accurate ODE solution, does not however lead to good training results. In particular, the loss function continued to increase over the iterations. On the other hand, despite being nearly 50 times slower than methods, the fixed-step forward Euler solver was successfully used for learning the dynamics of a 30-state system in the training configuration described in Appendix B One should however note that, in this configuration, the gains for the collision avoidance policy K. (which couples the ODE) were set to small values. This makes the system simpler and more stable than having larger gain. As result, if one attempts to train with the test configuration from Appendix BB where the gains are increased and the system is more unstable, then backpropagating trough Euler simply fails. Comparing Figures and it can seen that the learning curves our methods are unaffected by the change in the gains, while EKPR-Euler and ADJ-Bule: fail to decrease the loss.
7 RELATED WORK
RNN training pathologies. One of the first RNNs to be trained successfully were LSTMs (Greff et al. 2017), due to their particular architecture. Training an arbitrary RNN effectively is generally difficult as standard RNN dynamics can become unstable or chaotic during training and this can cause the gradients to explode and SGD to fail (Pascanu al. 2012). When RNNs consist of discretised ODEs, then stability of SGD intrinsically related to the size of the convergence region of the solver (Ciccone al. 2018). Since higher-order and implicit solvers have larger convergence region (Hairer
10
101
A-SMODE
di S-ODDD EKPR-Euler
a a-SNODE
Log)Loss)
102.
LogiLoss) Loss)
102
EXPR-Eullor
ADJ-Eular
AD.--Biler
10".
10a.
200
400
(200
200
1000
2000
400J
6000
8000
Time [s]
Figure 3: Multi-agent learning using different methods. Training loss vs. iterations (left) and execution time (right). Our methods converge one order of magnitude faster and lower loss. Note that changing the simulation parameters can make Euler unstable (see Figure 4p.
7

---

Published as a conference paper at ICLR 2020
10".
104.
Log(Loss)
10"-
6-80008
Log(Loss)
10'
SMCDE
EKPR Fer
6-SMODE
ADJ-Fuler
0-SKODR
10/
102
RMPR-Ful
ADJ-Buule
0:
200
100
(20
800
1000
2000
4000
6000
8000
# Iterations
Time [6]
Figure 4: Multi-agent learning with increased controller gains using different methods. Training loss vs. iterations (left) and execution time (right). Our methods remain stable and their learning curves are the same as Figure 3
and ADJ-Buler are not able to train in this configuration. The proposed algorithms are inherently more stable and robust than the baselines.
et al. 1993), following (Pascanu et al.
2012) it can be argued that our method has the potential to mitigate instabilities and hence make the learning efficient. This is supported by our results.
Unrolled architectures. In GGrves 2016), an RNN has been used with a stopping criterion, for iterative estimation with adaptive computation time. Highway (Srivastava al. 2015) and residual networks (He et af., (215)) have been studied in (Greff et al. 2016) as unrolled estimators. this context, (Haber & Ruthotto 2017) treated residual networks
autonomous discrete-ODEs and investigated their stability. Finally, in (Ciccone et al.) 2018) discrete-time non-autonomous ODE based residual networks has been made explicitly stable and convergent to input-dependant equilibrium, then used for adaptive computation.
Training stable ODEs. In (Haber & Ruthotto 2017 Ciccone et al. 2018), ODE stability conditions where used to train unrolled recurrent residual networks. Similarly, when using our method on be enforced by projecting the state weight matrices, A, into the Hurwitz stable space: i.e. A < O. At test time, overall stability will also depend on the solver (Durran 2010; Isermann 1989). Therefore, high order variable step method (e.g. DoPr5) should be used at test time in order to minimize the approximation error.
(Ciccone et al. 2018) ODE stability
Dynamics and machine learning. A physics prior neural network was used by (Jia al. 2018) in the form of consistency loss with data from a simulation. In (De Avila Belbute-Peree et al. 2018), differentiable physics framework was introduced for point mass planar models with contact dynamics. (Ruthotto & Haber, 2018) looked at Partial Differential Equations (PDEs) to analyze neural networks. while (Raissi Karniadakis 2018 Raissi al. (017)) used Gaussian Processes (GP) to model PDEs. The solution of linear ODE was used in (Soleimani et al. 2017) in conjunction with a structured multi-output GP to model patients outcome of continuous treatment observed at random times. (Pathak et al. 2017) predicted the divergence rate of a chaotic system with RNNs.
8 SCOPE AND LIMITATIONS
Test time and cross-validation At test time, since the future outputs are unknown, an explicit integrator needed. For cross-validation, the loss needs instead to be evaluated on a different dataset. In order to do available during cross-validation, a corresponding polynomial representation the form can be found and the relaxed loss 44 can be evaluated efficiently.
needs to solve the ODE forward in time. However, since the output data is
Nonsmooth dynamics. We have assumed that the ODE-Net dynamics has a regularity r > in order to take advantage of the exponential convergence of spectral methods, i.e., that their approximation error reduces as Ohh),, where is h is the size of the window used to discretize the interval. However, this might not be true in general. In these cases, the optimal choice would be to a hp-spectral approach (Canuto et al. 1988), where is reduced locally only near the discontinuities. This is very closely related to adaptive time-stepping for ODE solvers.
8

---

Published as a conference paper at ICLR 2020
Topological properties, convergence, and better generalization. There are few theoretical open questions stemming from this work. We argue that one reason for the performance improvement shown by algorithms is the fact that the set of functions generated by fixed neural network topology does not posses favorable topological properties for optimization, as discussed (Petersen et al. 2018). Therefore, the constraint relaxation proposed in this work may improve the properties of the optimization space. This is similar to interior point methods and can help with accelerating the convergence but also with preventing local minima. One further explanation is the fact that the proposed method does not suffer from vanishing nor exploding gradients, as shown in Appendinx Moreover, our approach very closely resembles the MAC scheme, for which theoretical convergence results are available (Carreira-Perpinan & Wang 2014
Synchronous vs Asynchronous. The proposed method can be used for an arbi- trary cascade of dynamical systems they can be expressed as single ODE, When only the final state of one ODE (or its trajectory) is fed into the next block, e.g. as in (Ciccone et al. 2018), the method could be extended by means of 2.M smaller optimizations, where M is the number of ODEs.
Multiple ODEs:
Hidden states. Latent states do not appear in the loss, so training and particularly initializing the polynomial coefficients is more difficult. A hybrid approach is to warm-start the optimizer using few iterations of backpropagation. We plan to investigate a full spectral approach in the future.
REFERENCES
J.C. Butcher and G. Wanner. Runge-Kutte methods: some historical notes. Applied Nu- merical Mathematics, 22(1-3):113-51, November 1996. ISSN 01689274. doi: 10.1016/ 80168-9274496)00048-77 URL p11/801689274960000
www.auu...llleeeeeeeeeeeee...eee com/ con/retrieve/
C.G. Canuto, M.Y. Hussaini, A. Quarteroni, and T.A. Zang. Spectral Methods in Fluid Dynamics. Springer- Verlag, 1988.
Miguel Carreira-Perpinan and Weiran Wang. Distributed optimization of deeply nested systems, In Artificial Intelligence and Statistics, pp. 10-19, 2014.
Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K. Duvenaud. Neural Ordinary Differential Equations. NeurIPS, 2018
Marco Ciccone, Marco Gallieri, Jonathan Masci, Christian Osendorfer, and Faustino Gomez. Nets-net: Stable deep networks from non -utoomooooouu differential equations. In NeurIPS, 2018
Filipe De Avila Belbute- Peres, Kevin A. Smith, Kelsey Allen, Joshua B. Tenenbaum, and Zico Kolter. End-to-end differentiable physics for learning and control. In NeurIPS, 2018
Dale R. Durran. Numerical Methods for Fluid Dynamics. Springer New York, 2010. doi: 10.1007/ 978- 1-4419-6412-0. URL O 1110777978811-441-641--------
Thor Fossen. Handbook of marine craft hydrodynamics and motion control. John Wiley & Sons, 2011.
Gene F. Franklin, J. David Powell, and Abbas Emami-Naeini. Feedback Control of Dynamic Systems (7th Edition). Pearson, 2014. ISBN 0133496597.
Amir Gholami, Kurt Keutzer, and George Biros. ANODE: Unconditionally Accurate Memory- Efficient Gradients for Neural ODEs. Preprint at arXiv:1902. 10298, 2019.
Ian Goodfellow, Yoshua Bengio, and Aaron Courville, Deep Learning. MIT Press, 2016. http:
77 7/www. www.deeepeear..nooooooooo
Alex Graves. Adaptive Computation Time for Recurrent Neural Networks. arXiv:1603.00 1603.08983 [cs], March 2016. (tRL:ttt:://aatt http:/awwiii.g....aa.oo//. haag//ab///ee /0t...60.... 08983
arXiv: 1603,0898,
Klaus Greff, Rupesh Kumar Srivastava, and Jürgen Schmidhuber. Highway and residual networks learn unrolled iterative estimation. CoRR, abs/1612.07771, 2016. URL http: // /rrxiv.ogg/
abs/ 1 61 2. 07771

---

Published as a conference paper at ICLR 2020
Klaus Greff, Rupesh Kumar Srivastava, Jan Koutník, Bas R. Steunebrink, and Jürgen Schmidhuber. LSTM: A Search Space Odyssey. IEEE Transactions on Neural Networks and Learning Systems, 28(10):2222 222:22222,,22, October 2017. ISSN 2162-237XX, 2162-2388. doi: 10. 10.1109/TNNLS.001 2582924. URL
0400 04069
arXiv: 1503.04069.
Eldad Haber and Lars Ruthotto. Stable Architectures for Deep Neural Networks. arXiv preprint arXiv:/705.082: 170..0334, 2017. URL.httts::/arxiv..oom abs/ 705. 03341
E. Hairer, S. P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations I (2Nd Revised. F: Nonstiff Problems. Berinner-Verlagee Berlin, Heidelberg. 1993. ISBN 0-387-56670-8.
Trevor Hastie, Robert Tibshirami, and Jerome Friedman. The Elements of Statistical Learning. Springer Series Statistics. Springer New York Inc., New York, NY, USA, 2001.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep Residual Learning for Image Recognition, 15 033385 [cs], December 2015. URL
03385 arXiv: 1512.03385.
Roger A. Hom and Charles R. Johnson. Matrix Analysis. Cambridge University Press, New York, NY, USA, 2nd edition, 2012. ISBN 0521548233, 9780521548236.
Sergey loffe and Christian Szegedy. Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift. arXiv:1502.03166 [cs], February 2015. URL http:
7/arxiv. org/ 111111022..11111.1
arXiv: 1502.03167.
Rolf Isermann. Digital Control Systems. Springer Berlin Heidelberg, 1989. doi: 10.1007/
118-3-642-86417-2
htttps::////:: ottp//...../011.0////111111//110
0042-86417-9
L. Jain and L. R. Medsker. Recurrent Neural Networks: Design and Applications (International Series on Computational Intelligence). CRC Press, 1999. ISBN 0849371813.
Xiaowei Jia, Anuj Karpatne, Jared Willard, Michael Steinbach, Jordan S. Read, Paul C. Hanson, Hilary A. Dugan, and Vipin Kumar. Physics guided recurrent neural networks for modeling dynamical systems: Application to monitoring water temperature and quality lakes. CoRR, abs/1810.02280,, 2018.
org/ hbtp/ab//aaa 1810. 02880
Kody Law, Andrew Stuart, and Kostas Zygalakis. Data assimilation. Cham, Switzerland: Springer, 2015.
Yann Lecun. A theoretical framework for back-propagation. In D. Touretzky, G. Hinton, and T. Se- jnowski (eds.), Proceedings of the 1988 Connectionist Models Summer School, CMU, Pittsburg. PA. pp. 21-28. Morgan Kaufmann, 1988.
Linnainmaa. The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors. Master's thesis, Univ. Helsinki, 1970.
Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. CoRR, abs/1211.5063, 2012.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang. Zachary De Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. In NIPS 2017 Workshop Autodiff, 2017.
Jaideep Pathak, Zhixin Lu, Brian R. Hunt, Michelle Girvan, and Edward Ott. Using Machine Learning to Replicate Chaotic Attractors and Calculate Lyapunov Exponents from Data. Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(12) 102, December 2017. ISSN 1054- 1500, 1089-7682. doi: 10.1063/1.50100 /tt.itppp///aaiii.oo/////////OO0 arXiv: 1710.07313.
07313
Philipp Petersen, Mones Raslan, and Felix Voigtlaender. Topological properties propertisllpoppettiess the set of functions generated by neural networks of fixed size. arXiv preprint (ote::/006005555 2018.
10

---

Published as a conference paper at ICLR 2020
Maziar Raissi and George Em Karniadakis. Hidden Physics Models: Machine Learning of Nonlinear Partial Differential Equations. Journal of Computational Physics, 357:125-141, March 2018. ISSN 00219991. doi: 10. 0..10166.jpp2201111..3999 http://////atttttx///////// www.wwwww...g......aaoereeeebeea www.aa////.......a 0057880888888888 00588 arXiv: 1708.00588.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Numerical Gaussian Processes for Time-dependent and Non-linear Partial Differential Equations. 1..3v.7703..1000 [cs, math, stat], March 2017. URL htt.: /ttp:://axxivv .rg..
org/ aaa//b////// 10230
arXiv: 1703.10230.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Multistep Neural Networks for Data- driven Discovery of Nonlinear Dynamical Systems, arxiv:1801.0111 [nlin, physicx:hyssics, stat], January 2018. URL
11w...11111100010 01236
arXiv: 1801.01 1.01.00236.
I. Michael Ross. A primer Pontryagin's principle in optimal control, volume 2. Collegiate publishers San Francisco, 2009.
I. Michael Ross and Mark Karpenko. A review of pseudospectral optimal control: From theory to flight. Annual Reviews in Control, 36(2):182-19) 2012.
C. Runge. Ueber die numerische Aufli von Differentialglecchungen Mathematische Annalen, 46(2):167-177, June 1895. ISSN 0025-5831, 1432-1807. doi: (0.1007/BF01444800. URL
Lars Ruthotto and Eldad Haber. Deep Neural Networks Motivated by Partial Differential Equa- tions. arXiv: arXiv:/804.00 1804.04222 [es, math, staf], April 2018. URL http:://axxv.....t
1804
04272 arXiv: 1804.04272.
Bruno Siciliano, Lorenzo Sciavicco, Luigi Villani. and Giuseppe Oriolo. Robotics: Modelling, Planning and Control (Advanced Textbooks in Control and Signal Processing). Springer, 2008
Hossein Soleimani, Adarsh Subbaswamy, and Suchi Saria. Treatment- Models for Counter- factual Reasoning with Continuous-time Continuous valuue Interventions. arXiv:1704.0200 [cs. stat], April 2017.
org/ abs 1704. 02038
arXiv: 1704.02038.
Rupesh Kumar Srivastava, Klaus Greff, and Jürgen Schmidhuber.
Highway Networks.
arXiv:1505.0030 00387 [cs], May 2015. URL http: http:://axxi..or//aapp... 0055//5../1510.0000 00387 1505.00387.
arXiv:
Robert F. Stengel. Optimal control and estimation. Dover books on advanced mathematics. Dover Publications, New York, 1994. ISBN 978-0-486-68200-66
P. J. Werbos. Applications of advances in nonlinear sensitivity analysis. Proceedings of the 10th IFIP Conference, 31.8 4.9, NYC, pp. 762-770, 1981.
Julian Georg Zilly, Rupesh Kumar Srivastava, Jan Koutník, and J urgen Schmidhuber. Recurrent Highway Networks. arXiv: [cs], July 2016. URL
1607. 03474 arXiv: 1607.03474.

---

Published as a conference paper at ICLR 2020
APPENDIX
A ADDITIONAL MATERIAL FOR THE VEHICLE DYNAMICS EXAMPLE
The model is formulated in a concentrated parameter form Siciiliano et al. 2008 We follow the notation of (Fossen 2011). Recall the system definition:
J(n). Mi + d(v) + CCU) u,
where U FSS are the states, namely, the I, and coordinates in fixed (world) frame, the vehicle orientation with respect this this frame, and the body-frame velocities, Ug. Vy. and angular rate, w. The input is set of torques in the body-frame, (Fz, 0, Tzy). The Kinematic matrix is
(os sin(e)
(in cos ())
0
J(n) =
1
the mass matrix is
m
0
M
0
0
0
0
where the vehicle mass and / represents the rotational intertia. The Coriolis matrix
mi
C(v)
miw
0
0
0
and the damping force d(v) = kgv. We set m = 1 and Kς 1. The input, u, comes from
Fourier series with fundamental amplitude 1.
Loss after seconds
101
102
100
10-1
10
500
1000
(a)
(b) 1 100
(c) 300
(d) Loss function
Loss after seconds
101
10'
10:
107
10
10
500
1000
(0) 20
(f) 100
(g) 300
(h) Loss function
Figure 5: Vehicle dynamics learning in high-data regime. Top row: 8-SNODE method. Bottom row: fifth-order Dormand-Prine method (Butcher & Wanner, 1996) with backpropagation. First three columns: comparison of true trajectories (red) with the prediction from the surrogate (black) at different iterations of the optimization. Last column: loss function at each iteration. 8-SNODE has faster convergence than Euler.
12

---

Published as a conference paper at ICLR 2020
Loss after 21.971 seconds
10'
109
10
(a) = 20
ib) 60
(c) i= 100
(d)
100
Loss after 116 6.6944 seconds
10r
10a
107
10"
(c) = 20
(1) i 100
(g) i 300
(h) Loss function
Figure 6: Vehicle dynamics learning in high-data regime. Top row: a-SNODE method. Bottom
Explicit Euler method with back-propagation First three columns: comparison of true trajecto- ries (red) with the prediction from the surrogate (black) at different iterations of the optimization. Last column: value of the loss function at each iteration.
row:
Figure 7: Vehicle dynamics learning with perturbed input data for the a-SNODE method. Left: comparison of true trajectories (red) with the input data (black). Right: comparison in the time domain for the states of the first sample of the batch.
B ADDITIONAL MATERIAL FOR THE MULTI-AGENT EXAMPLE
The multi- agent simulation consists of N₂ kinematic vehicles:
N. Σ Ko(m.)))) Jt
ni. = (,,), U, tanh w + K.(T)) +
N.
where 7) E R3N. are the states for each vehicle, namely, the æi, ys positions and the orientation Φ; of vehicle the world frame, while e; R2N. are the controls signals, the form of linear and angular velocities, P4, wi. The kinematics matrix is
cos (i
0 sin
((n:
0
1
13

---

Published as a conference paper at ICLR 2020
The agents velocities are determined by known arbitrary control and collision avoidance policies, respectively, Kg and Kg. In particular:
ky
]. Ko(n, 11) = - -kvoe - ] kooddrs
Ka(m:)
kooØ;
where 60, atan2 ((-3i): (-xi)) oi and d0.5 atan2 ((ys Vj); (x; - xzy) Oi.
Training configuration. We set
ky 0.05 Ke 0.1, Kyo 0.001,
Kgo 0.01, 4, 0.01.
The signal v w(t) is generated by Fourier series with fundamental amplitude 0.1.
Test configuration.
We change Ko(ni, 75)
kine to kipndors
Kis (1/ mj)
with Koo 0.05 and k2o 0.1. We also set to
0.
(a) a-SNODE 50% data
(b) Euler 50% data
(c) D-SNODE 20% data
(d) Buler 20% data
Figure 8: Multi-agent testing in low-data regime.
C GRADIENT ANALYSIS
Consider the classic discrete-tine RNN:
x(t 1) f(x(t),
(10)
Then, given a loss -= Eitto ((x(t),, the following gradients are used during training
ac al(x(t)) dr(()
(11)
00
dx(()
00
where, for any the chain rule gives:
(x(++1) 0f(x(t), u(t);
Or(()
+ (((())
(12)
00
00
00
(xx(t) f(x(t), 0(
dx(1)
Iteration of (12) is the main principle of backpropagation through time (Goodfellww al. 2016). A
known fact is that iterating (12) is similar to geometric series. Therefore, depending on the spectral
radius of the Jacobian, (x)))) it can result in vanishing (p 1) or exploding (p > gradients (Zilly et al. 2016).
We can now demonstrate that, by avoiding function compositions, our approach immune to the exploding gradient problem. In particular, our gradients are fully time-independent and their accuracy not affected by the sequence length. Recall the ODE:
*(1) (())))
(13)
The following result obtained:
We consider single batch element: f and in R==x where ny is the dimensionality.
14

---

Published as a conference paper at ICLR 2020
Theorem 1. Assume the full Jacobian of (x((t), u(t); has finite spectral radius for all Then, the norms of the gradients used in Algorithm anni Algorithm 2 also finite for all t.
We will prove the theorem for Algorithm since it includes the relevant parts of Algorithm
Proof. Given the Legendre basis functions, (t). i 1, ........p.. define the ODE residual as:
r(z(t), 0) Σ 2,1i(t) - (),ds(t).
for which the residual loss is (() ||(xx(t), 6)1? Then, Algorithm consists of the concurrent minimization of the relaxed loss function:
((((t))) + R(x((),
with respect to the coefficients z, the Legendre polynomial given 8", and of the residual loss:
R(x(t),0)dl
with respect to 0 given the set of coefficients x i = I....,p. For the residual loss gradients are:
aR(x*(t), - 2r Of (Ein of (t)) u(t):
(14)
00
00
where there is no recursion over the previous values x(t 7), since the basis functions (t) are given and the points x are treated as data. By assumption, the Jacobian of fx((t), 0) has finite singular values. Hence, by standard matrix norm identities (see Chapter ofrn & Johnson (2012) the gradient of R(x(t), 0) with respect to 0 has finite norm for all t.
Similarly, the absence of time recursions the gradient for the v; update using the relaxed loss also follows trivially from the fact that the coefficients z; are independent variables, that assume given 0*, that the basis functions s/s (4) are fixed. Then, the claim follows again from the assumption on the Jacobian of fI, u; 0).
Note that the result of Theorem I cannot easily be achieved when training ODEs using backpropaga- tion through the solver the adjoint method unless some gradient conditioning, such clipping or batch normalization (loffe & Szegedy 2013), is performed after applying f. On the other hand, our result relies only on the gradient of f being finite. This is trivial for shallow f and, for very deep f, the methods in (Toffe & Szegedy 2015 Srivastava et al. 2015] Ciccone et al. [018 He et al.) 2015) can applied just inside the architecture J. necessary. This fundamental difference with standard RNN training allows for f to have unrestricted Jacobian magnitudes which are needed to effectively model instabilities and long-short term dynamics, similar LSTMs (Greff et al. 2017).
15


Reviews:
[{"review_id": "SkelTPIi5B", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "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 work extends prior work on Neural ODEs.  From what I understand, the Neural ODE approach builds off the idea of representing the sequence of transformations of a hidden state (in residual nets, RNNs, etc.) as an ODE parameterized by a neural network.   In the original paper, the network is optimized via gradients calculated by the adjoint sensitivity method.  This paper puts forth the following contributions: a compact representation of the state transition function as a combination of Legendre polynomials, and an optimization scheme whose error is tied to the polynomial order and whose structure lends itself easily to parallelization.  The authors also demonstrate their model on an experiment on planar vehicle dynamics, in which their model is shown to have improved predictive quality and efficiency.\n\nI am inclined to accept this paper, due to its various contributions on speed and performance, with the caveat for some clarifications on how the experiments were conducted and compared.  Given these clarifications in an author response, I would be willing to increase the score.\n\nPros of the paper:\n\t1) Trajectory predictions on the two planar vehicle dynamics experiments was impressive.\n\t2) The proposed representation of the state transition dynamics is indeed more memory-efficient, and its approximation error is modulatable by the hyperparameter order \\p.\n\t3) Speedup due to parallelization is substantial.\nCons of the paper:\n\t1) The experiments did not display a proper comparison against the hybrid method mentioned in Section 1.  The experiments also did not compare against adjoint methods in the multi-agent example, or in the low-data regime for the single-agent example.  Instead, the experiments mostly highlighted the problems with direct backpropagation through the ODE solver, which is already well-known to have issues in robustness and stability.  While it is nice to have empirical results that showcase this, a more comprehensive comparison against current adjoint methods would be more interesting, especially in the multi-agent example.\n\t2) It slightly detracts from the cleanliness of the story that we must first create an initial trajectory, before performing our coordinate descent.  \n\nQuestions and Points of Confusion:\n\t1) What intuits the choices of the Legendre polynomial as your set of basis functions, and the Gauss-Lobatto scheme to select collocation points, instead of alternative candidates?\n\t2) In Section 5, it was mentioned that \"100 equally-spaced points produce a comparable result\" to the Gauss-Lobatto quadrature points.  Furthermore, it was mentioned that these evenly-spaced collocation points were used for experiments - was this the case for all experiments?  If so, then what purpose does Gauss-Lobatto play in the paper? \n\t3) In Figure 1, were the DoPr5 and Euler plots shown for the backpropagation or adjoint method?  If it was the backpropagation method, the plots for the adjoint method would be highly interesting to show.\n\t4) From what I understand from Section 4, when we perform step 0 and step 2 to update the trajectory, we simply optimize the coefficients \\x_i and \\u_i directly to minimize the current objective as both \\x(t) and \\u(t) are represented as a combination of Legendre polynomials.  However, in Equation 8, for the planar vehicle dynamics, \\u is deterministically generated from the current states and network weights.  It makes sense to add structure to \\u, as we need a way to make sure the inputs can indeed generate the trajectory of \\x.  Does this mean the spectral method is not performed for \\u as was detailed in Equation 5?\n\t5) I am confused about how the gray-box models are built in Section 5.  It states that \"For \\f_J, sin(\\phi) and cos(\\phi) are used as input features, where \\phi is the vehicle orientation.\"  Does that mean that only \\phi from \\eta is passed in as an input, both sin(\\phi) and cos(\\phi) are passed in as inputs, or the entire current \\eta is passed in as input?  And is the output a new \\eta, which is then structured into the 3x3 matrix \\J(\\eta) in Appendix A?  Or does it simply output the matrix directly for the given value of \\phi.  I am confused because it is written that \\J(\\eta) is equal to \\f_J(\\eta;\\theta_J), but then in the appendix it is written that \\J(\\eta) is equal to the matrix - so where is the network?  This confusion extends to the other gray-box models \\C(v) and \\d(v).\n\t6) It is written in Equation 3 that the residual also takes in the input \\u(t).  However, in the experiments, namely Equation 8, it does not appear like \\u is used to calculate the residual at all.\n\t7) What motivated the use of the planar vehicle dynamics experiment to showcase your model?  Were there other baselines or benchmarks you considered or attempted?\n", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly."}, "extra_text": {}}, {"review_id": "H1lLFOChKS", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "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 work proposes a new approach for the evolution of Neural ODEs for particle systems. The authors suggest to replace the traditional backpropagation through ODEs or the recent adjoint method for backpropagation and instead solve the problem as an alternating optimization scheme. In particular it is suggested that using spectral methods (Legendre's polynomials) first compute a minimizer of the trajectory (optimizing a trajectory x(t)) given the Loss/Langrangian (that is based on the data times of training). After an initial trajectory is computed, follow an alternating minimization, where in the first step, minimize the discrepancy between the network (that describes the time change of the ODE) and the time derivative of the current trajectory. In the second step, re-compute the trajectory with the new updated network and modified lagrangian to update the network parameters again. This two step optimization is applied back and forth until a residual condition is reached, i.e. the loss is small. The network's parameters in the two stages are optimized via SGD and ADAM respectively. Further, to perform the required numerical integration in each step the authors apply Gaussian quadrature and turn the overall optimization scheme into the repeating application of a finite number of gradient updates in an alternating fashion (as has become popular in recent Deep Learning approaches e.g. GANs). The authors test the new method on different particle systems' trajectories and observe a numerical speed-up and improved accuracy as compared with previous approaches.\n\nAdmittingly, I am not an expert in Neural ODEs and am not too familiar with the literature investigating neural networks for modelling differential equations and control. The paper, however is well written and the presentation is very nice in my opinion. It is hard for me to judge how original some of the ideas presented but from my perspective they seem quite solid.\n\nOverall, I am currently voting for weak accept, for solid presentation and content, but with with the following problems:\n\nIt seems that Neural ODEs are most beneficial, over the alternatives, when used with irregular data times and or sparse number of time points. This is a point that is mostly missing in the discussion and the experimental section, from my understanding the experimental section uses equally spaced intervals. If this is the case, I do not find them sufficient and I would hope to see how does this method perform when trained with sparser and or irregular time points. Currently the method presented is illustrated as an effective algorithm for noisy ODE solvers. For this reason, I am also wondering whether this is also the right venue to present this interesting work.\n\nOther small things:\nIn the second sentence\n\"ODE-Nets have been shown to provide superior performance with respect to classic\nRNNs on time series forecasting with sparse training data.\"\nIt could be nice to provide a reference illustrating the improved performance of Neural ODEs over RNNS on a time series forecasting task.\n", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}, {"review_id": "BJg92wGiYS", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "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": "TITLE\nSNODE: Spectral Discretization of Neural ODEs for System Identification\n\nREVIEW SUMMARY\nExceptionally clear and well written paper that demonstrates a strong improvement on an important problem. Novel, timely, and of broad interest.\n\nPAPER SUMMARY\nThe paper presents a new method for estimating parameters in a neural ODE based on a polynomial representation of trajectories and alternating updates of trajectories and neural net parameters.\n\nCLARITY\nThe presentation is exceptionally clear.\n\nORIGINALITY\nTo my knowledge the proposed method is novel.\n\nSIGNIFICANCE\nThe paper demonstrates a strong and practically significant improvement in learning, and I expect the results will be of interest to everybody working in this area.\n\nFURTHER COMMENTS\n\nIn eq. 2, \"X\" (blackboard typeface) is not defined. At this point, it is not clear how x(t) is represented.\n\n\"trades-off\" -> trades off\n\n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This work proposes using spectral element methods to speed up training of ODE Networks for system identification. The authors utilize truncated series of Legendre polynomials to analyze the dynamics and then conduct experiments that shows their proposed scheme achieves an order of magnitude improvement in training speed compared to baseline methods. Reviewers raised some concerns (e.g. empirical comparison against adjoint methods in the multi-agent example) or asked for clarifications (e.g. details of time sampling of the data). The authors adequately addressed most of these concerns via rebuttal response as well as revising the initial submission. At the end, all reviewers recommended for accept based on contributions of this work on improving training speed of ODE Networks. R4 hopes that some of the additional concerns that are not yet reflected in the current revision, be addressed in the camera ready version. 
Decision: Accept (Poster)