SNODE: SPECTRAL DISCRETIZATION OF NEURAL ODES FOR SYSTEM IDENTIFICATION
Year: 2020

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 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 (Chen et al., 2018), 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.

heading: 1 INTRODUCTION
text: 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. In complex control systems and model-based reinforcement learning, planning over a long horizon is often needed, while high frequency feedback is 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. ODE-Nets have been shown to provide superior performance with respect to classic RNNs on time series forecasting with sparse training data. However, learning their parameters can be computationally intensive. In particular, ODE-Nets are memory efficient but time inefficient. In this paper, we address this bottleneck and propose a novel alternative strategy for system identification.
Summary of contributions. We propose SNODE, a compact representation of ODE-Nets for system identification with full state information that makes use of a 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 3 and is based on backpropagation (Linnainmaa, 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 on a 6-state vehicle problem, where it is at least one order or magnitude faster in each optimizer iteration than explicit and adjoint methods, while convergence is achieved in a third of the iterations. At test time, the MSE is reduced by one order of magnitude. In Section 6, the method is used for a 30-state system consisting
∗Corresponding author.
of identical vehicles, coupled via a known collision avoidance policy. Again, our method converges in a third of the iterations required by backpropagation thourgh a solver and each iteration is 50x faster than the fastest explicit scheme.

heading: 2 NEURAL ORDINARY DIFFERENTIAL EQUATIONS
text: The minimization of a scalar-valued loss function that depends on the output of an ODE-Net can be formulated as a general constrained optimization problem:
min θ∈Rm ∫ t1 t0 L(t, x(t)) dt, (1)
s.t. ẋ(t) = f(t, x(t), u(t); θ),
x(t0) = x0,
where x(t) ∈ X is the state, u(t) ∈ U 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 x(t) and u(t) (see Section 8). Equation (1) can be used to represent several inverse problems, for instance in machine learning, estimation, and optimal control (Stengel, 1994; Law et al., 2015; Ross, 2009). Problem (1) can be solved using gradient-based optimization through several time-stepping schemes for solving the ODE. (Chen et al., 2018; Gholami et al., 2019) have proposed to use the adjoint method when f is a 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 gradients ∂L/∂θ using backpropagation through a discrete approximation of the constraints, such as Runge-Kutta methods (Runge, 1895; Butcher & Wanner, 1996) or multistep solvers (Raissi et al., 2018). This ensures that the solution remains feasible (within a numerical tolerance) at each iteration of a 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 a nonlinear equation at each step, 3) the numerical error can significantly affect the solution, and 4) the problem topology can be unsuitable for optimization (Petersen et al., 2018).
Limitations of adjoint methods. ODE-Nets (Chen et al., 2018) solve (1) using the adjoint method, which consists of simulating a dynamical system defined by an appropriate augmented Hamiltonian (Ross, 2009), with an 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 a 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 of the adjoint method. A similar 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.

heading: 3 RELAXATION OF THE SUPERVISED LEARNING PROBLEM
text: Our algorithm is based on two ingredients: i) the discretization of the problem using spectral elements leading to SNODE, detailed in Section 4, and ii) the relaxation of the ODE constraint from (1), 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 our examples.
The problem in (1) is split into two smaller subproblems: one finds the trajectory x(t) that minimizes an unconstrained relaxation of (1). The other trains the network weights θ such that the trajectory becomes a solution of the ODE. Both are addressed using standard gradient descent and backpropagation. 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 are discussed.
Step 0: Initial trajectory. The initial trajectory x(t) is chosen by solving the problem
min x(t)∈X ∫ t1 t0 L(t, x(t)) dt, (2)
s.t. x(t0) = x0,
If this problem does not have a unique solution, a regularization term is added. For a quadratic loss, a closed-form solution is readily available. Otherwise, a prescribed number of SGD iterations is used.
Step 1: Coordinate descent on residual. Once the initial trajectory x?(t) is found, θ is computed by solving the unconstrained problem:
min θ∈Rm ∫ t1 t0 R(t, x?(t), θ) dt, R(t, x(t), θ) := ‖ẋ(t)− f(t, x(t), u(t); θ)‖2 (3)
If the value of the residual at the optimum θ∗ is smaller than a prescribed tolerance, then the algorithms stops. Otherwise, steps 1 and 2 are iterated until convergence.
Step 2: Coordinate descent on relaxation. Once the candidate parameters θ? are found, the trajectory is updated by minimizing the relaxed objective:
min x(t)∈X ∫ t1 t0 γL(t, x(t)) +R(t, x(t), θ?) dt, (4)
s.t. x(t0) = x0.
Discussion. The proposed algorithm can be seen as an alternating coordinate gradient descent on the relaxed functional used in problem (4), i.e., by alternating a minimization with respect to x(t) and θ. If γ = 0, multiple minima can exist, since each choice of the parameters θ would induce a different dynamics x(t), solution of the original constraint. For γ 6= 0, the loss function in (4) trades-off the ODE solution residual for the data fitting, providing a unique solution. The choice of γ implicitly introduces a satisfaction tolerance (γ), i.e., similar to regularized regression (Hastie et al., 2001), implying that ‖ẋ(t)− f(t, x(t); θ)‖ ≤ (γ). Concurrently, problem (3) reduces the residual.
4 SNODE – HIGH-ORDER DISCRETIZATION OF THE RELAXED PROBLEM
In order to numerically solve the problems presented in the previous section, a discretization of x(t) is needed. Rather than updating the values at time points ti 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), as truncated series:
x(t) = p∑ i=0 xiψi(t), u(t) = z∑ i=0 uiζi(t), (5)
where xi, ui ∈ R and ψi(t), ζi(t) are sets of given basis functions that span the spaces Xh ⊂ X and Uh ⊂ U. In this work, we use orthogonal Legendre polynomials of order p (Canuto et al., 1988) for ψi(t), where p is a hyperparameter, and the cosine Fourier basis for ζi(t), where z is fixed.
Collocation and quadrature. In order to compute the coefficients xi of (5), we enforce the equation at a discrete set Q of collocation points tq. Here, we choose p+ 1 Gauss-Lobatto nodes, which include t = t0. This directly enforces the initial condition. Other choices are also possible (Canuto et al., 1988). Introducing the vectors of series coefficients xI = {xi}pi=0 and of evaluations at quadrature points x(tQ) = {x(tq)}q∈Q, the collocation problem can be solved in matrix form as
xI = M −1x(tQ), Mqi := ψi(tq). (6)
We approximate the integral (3) as a sum of residual evaluations over Q. Assuming that x(0) = x0, the integrand at all quadrature points tQ can be computed as a component-wise norm
R(tQ, x(tQ), θ) = ‖DM−1x(tQ)− f (tQ, x(tQ), u(tQ); θ) ‖2, Dqi := ψ̇i(tq). (7)
Fitting the input data. For the case when problem (2) admits a unique a solution, we propose a new direct training scheme, δ-SNODE, which is summarized in Algorithm 1. In general, a least-squares approach must be used instead. This entails computing the integral in (2), which can be done by evaluating the loss function L at quadrature points tq. If the input data is not available at tq, we approximate the integral by evaluating L at the available time points. The corresponding alternating coordinate descent scheme α-SNODE is presented in Algorithm 2. In the next sections, we study the consequences of a low-data scenario on this approach.
We use fixed numbers Nt and Nx of updates for, respectively, θ and x(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 δ-SNODE training Input: M,D from (6)-(7). x∗I ← argminx L(Mx) while R(Mx∗I , θ∗) > δ do θ∗ ← ADAMθ [R(Mx∗I , θ)]
end while Output: θ∗
Algorithm 2 α-SNODE training Input: M,D, θ∗ ∈ Rm, x∗I ∈ Rp, γ > 0. while γL(Mx?I) +R(Mx∗I , θ∗) > δ do
for i = 0 . . . Nx do x∗I ← SGDx [γL(Mx) +R(Mx, θ∗)] end for for i = 0 . . . Nt do θ∗ ← ADAMθ [R(Mx∗I , θ)]
end for end while Output: θ∗
Ease of time parallelization. If R(tq) = 0 is enforced explicitly at q ∈ Q, then the resulting discrete system can be seen as an implicit time-stepping method of order p. However, while ODE integrators can only be made parallel across the different components of x(t), 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.
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 xi coefficients.
Stability and vanishing gradients. The forward Euler method is known to have a 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-order 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 a more stable and robust training. This
claim is supported by our experiments. Reducing the time step can improve the Euler accuracy but it can still lead to vanishing or exploding gradients (Zilly et al., 2016; Goodfellow et al., 2016). In Appendix C, we show that our methods do not suffer from this problem.
Experiments setup and hyperparameters. For all experiments, a common setup was employed and no optimization of hyperparameters was performed. Time horizon T = 10s and batch size of 100 were used. Learning rates were set to 10−2 for ADAM (for all methods) and 10−3 for SGD (for α-SNODE). For the α-SNODE method, γ = 3 and 10 iterations were used for the SGD and ADAM algorithms at each epoch, as outlined in Algorithm 2. The initial trajectory was perturbed as x0 = x0 + ξ, ξ ∼ U(−0.1, 0.1). This perturbation prevents the exact convergence of Algorithm 1 during initialization, allowing to perform the alternating coordinate descent algorithm.

heading: 5 MODELING OF A PLANAR VEHICLE DYNAMICS
text: Let us consider the system
η̇ = J(η)v, Mv̇ + d(v) + C(v)v = u, (8)
where η, v ∈ R3 are the states, u = (Fx, 0, τxy) is the control, C(v) is the Coriolis matrix, d(v) is the (linear) damping force, and J(η) encodes the coordinate transformation from the body to the world frame (Fossen, 2011). A gray-box model is built using a neural network for each matrix
Ĵ(η) = fJ(η; θJ), Ĉ(v) = fC(v; θC), d̂(v) = fd(v; , θd).
Each network consists of two layers, the first with a tanh activation. Bias is excluded for fC and fd. For fJ , sin(φ) and cos(φ) are used as input features, where φ is the vehicle orientation. When inserted in (8), these discrete networks produce an ODE-Net that is a surrogate model of the physical system. The trajectories of the system and the learning curves are shown in Appendix A.
Comparison of methods in the high-data regime. In the case of δ-SNODE and α-SNODE, only p + 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 a 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 1a shows that δ-SNODE outperforms BKPR-DoPr5 by a factor of 50, while producing a significantly improved generalization. The speedup reduces to 20 for α-SNODE, which however yields a further reduction of the testing MSE by a factor of 10, as can be seen in Figure 1.
Comparison in the low-data regime. The performance of the methods was compared using fewer time points, randomly sampled from a uniform distribution. For the baselines, evenly-spaced points
were used. Table 1b shows that α-SNODE preserves a good testing MSE, at the price of an increased number of iterations. With only 25% of data, α-SNODE is 10x faster than BKPR-DoPr5. Moreover, its test MSE is 1/7 than BKPR-DoPr5 and up to 1/70 than BKPR-Euler, showing that the adaptive time step of DoPr5 improves significantly the baseline but it is unable to match the accuracty of the proposed methods. The adjoint method produced the same results as the backprop (±2%).

heading: 6 LEARNING A MULTI-AGENT SIMULATION
text: Consider the multi-agent system consisting of Na kinematic vehicles:
η̇i = J(ηi) tanh w +Kc(ηi) + 1 Na Na∑ j 6=i Ko(ηi, ηj)  , (9) where η ∈ R3Na are the states (planar position and orientation), J(ηi) 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, Kc and Ko plus some additional high frequency measurable signal w = w(t), shared by all vehicles. The control laws are non-linear and are described in detail in Appendix B. We wish to learn their kinematics matrix by means of a neural network as in Section 5. The task is simpler here, but the resulting ODE has 3Na states, coupled by K0. We simulate Na = 10 agents in series.
Comparison of methods with full and sparse data. The learning curves for high-data regime are in Figure 3. For method α-SNODE, training was 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 L̄ = 0.01. Table 2 summarizes results. δ-SNODE is the fastest method, followed by α-SNODE which is the best performing. Iteration time of
BKPR-Euler is 50x slower, with 14x worse test MSE. ADJ-Euler is the slowest but its test MSE is in between BKPR-Euler and our methods. Random down-sampling of the data by 50% and 20% (evenly-spaced for the baselines) makes ADJ-Euler fall back the most. BKPR-DoPr5 failed to find a time step meeting the tolerances, therefore they were increased to rtol= 10−5, atol= 10−7. Since the loss continued to increase, training was terminated at 200 epochs. ADJ-DoPr5 failed to compute gradients. Test trajectories are in Figure 2. Additional details are in Appendix B.
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 our 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 Ko (which couples the ODE) were set to small values. This makes the system simpler and more stable than having a larger gain. As a result, if one attempts to train with the test configuration from Appendix B, where the gains are increased and the system is more unstable, then backpropagating trough Euler simply fails. Comparing Figures 3 and 4, it can be seen that the learning curves of our methods are unaffected by the change in the gains, while BKPR-Euler and ADJ-Euler fail to decrease the loss.

heading: 7 RELATED WORK
text: 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 et al., 2012). When RNNs consist of discretised ODEs, then stability of SGD is intrinsically related to the size of the convergence region of the solver (Ciccone et al., 2018). Since higher-order and implicit solvers have larger convergence region (Hairer
et al., 1993), following (Pascanu et al., 2012) it can be argued that our method has the potential to mitigate instabilities and hence to make the learning more efficient. This is supported by our results.
Unrolled architectures. In (Graves, 2016), an RNN has been used with a stopping criterion, for iterative estimation with adaptive computation time. Highway (Srivastava et al., 2015) and residual networks (He et al., 2015) have been studied in (Greff et al., 2016) as unrolled estimators. In this context, (Haber & Ruthotto, 2017) treated residual networks as autonomous discrete-ODEs and investigated their stability. Finally, in (Ciccone et al., 2018) a discrete-time non-autonomous ODE based on residual networks has been made explicitly stable and convergent to an 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 (Ciccone et al., 2018) ODE stability can be enforced by projecting the state weight matrices, A, into the Hurwitz stable space: i.e. A ≺ 0. At test time, overall stability will also depend on the solver (Durran, 2010; Isermann, 1989). Therefore, a high order variable step method (e.g. DoPr5) should be used at test time in order to minimize the approximation error.
Dynamics and machine learning. A physics prior on a neural network was used by (Jia et al., 2018) in the form of a consistency loss with data from a simulation. In (De Avila Belbute-Peres et al., 2018), a 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 et al., 2017) used Gaussian Processes (GP) to model PDEs. The solution of a 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.

heading: 8 SCOPE AND LIMITATIONS
text: Test time and cross-validation At test time, since the future outputs are unknown, an explicit integrator is needed. For cross-validation, the loss needs instead to be evaluated on a different dataset. In order to do so, one needs to solve the ODE forward in time. However, since the output data is available during cross-validation, a corresponding polynomial representation of the form (5) can be found and the relaxed loss (4) can be evaluated efficiently.
Nonsmooth dynamics. We have assumed that the ODE-Net dynamics has a regularity r > p in order to take advantage of the exponential convergence of spectral methods, i.e., that their approximation error reduces as O(hp), 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 use a hp-spectral approach (Canuto et al., 1988), where h is reduced locally only near the discontinuities. This is very closely related to adaptive time-stepping for ODE solvers.
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 our algorithms is the fact that the set of functions generated by a fixed neural network topology does not posses favorable topological properties for optimization, as discussed in (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 C. Moreover, our approach very closely resembles the MAC scheme, for which theoretical convergence results are available (Carreira-Perpinan & Wang, 2014).
Multiple ODEs: Synchronous vs Asynchronous. The proposed method can be used for an arbitrary cascade of dynamical systems as they can be expressed as a 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 2M smaller optimizations, where M is the number of 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.

heading: A ADDITIONAL MATERIAL FOR THE VEHICLE DYNAMICS EXAMPLE
text: The model is formulated in a concentrated parameter form (Siciliano et al., 2008). We follow the notation of (Fossen, 2011). Recall the system definition:
η̇ = J(η)v, Mv̇ + d(v) + C(v)v = u,
where η, v ∈ R3 are the states, namely, the x, and y coordinates in a fixed (world) frame, the vehicle orientation with respect this this frame, φ, and the body-frame velocities, vx, vy , and angular rate, ω. The input is a set of torques in the body-frame, u = (Fx, 0, τxy). The Kinematic matrix is
J(η) = [ cos(φ) − sin(φ) 0 sin(φ) cos(φ) 0
0 0 1
] ,
the mass matrix is
M = [ m 0 0 0 m 0 0 0 I ] ,
where m is the vehicle mass and I represents the rotational intertia. The Coriolis matrix is
C(v) = [ 0 −mω 0 mω 0 0 0 0 0 ] ,
and the damping force is d(v) = kdv. We set m = 1 and Kd = 1. The input, u, comes from a Fourier series with fundamental amplitude 1.

heading: B ADDITIONAL MATERIAL FOR THE MULTI-AGENT EXAMPLE
text: The multi-agent simulation consists of Na kinematic vehicles:
η̇i = J(ηi)vi, vi = tanh w +Kc(ηi) + 1 Na Na∑ j 6=i Ko(ηi, ηj)  , where η ∈ R3Na are the states for each vehicle, namely, the xi, yi positions and the orientation φi of vehicle i in the world frame, while vi ∈ R2Na are the controls signals, in the form of linear and angular velocities, νi, ωi. The kinematics matrix is
J(ηi) = [ cos(φi) 0 sin(φi) 0
0 1
] .
The agents velocities are determined by known arbitrary control and collision avoidance policies, respectively, Kc and Ko. In particular:
Kc(ηi) =
[ kv
kφδφi
] , Ko(ηi, ηj) = [ −kvoe− d ls e−|δφij+π/2|
kφoδφij
] ,
where δφi = atan2 ((−yi), (−xi))− φi, and δφij = atan2 ((yi − yj), (xi − xj))− φi.
Training configuration. We set kv = 0.05 kφ = 0.1, kvo = 0.001, kφo = 0.01, ls = 0.01.
The signal w = w(t) is generated by a Fourier series with fundamental amplitude 0.1.
Test configuration. We change Ko(ηi, ηj) to:
Ktest(ηi, ηj) =
[ −kvoe− d ls
kφoδφij
] ,
with kvo = 0.05 and kφo = 0.1. We also set w = 0.

heading: C GRADIENT ANALYSIS
text: Consider the classic discrete-time RNN: x(t+ 1) = f(x(t), u(t); θ). (10)
Then, given a loss L = ∑tf t=t0 L(x(t)), the following gradients are used during training1:
∂L ∂θ = tf∑ t=t0 ∂L(x(t)) ∂x(t) ∂x(t) ∂θ , (11)
where, for any t, the chain rule gives: ∂x(t+ 1)
∂θ = ∂f(x(t), u(t); θ) ∂θ + J(x(t)) ∂x(t) ∂θ , (12)
J(x(t)) = f(x(t), u(t); θ)
∂x(t)
Iteration of (12) is the main principle of backpropagation through time (Goodfellow et al., 2016). A known fact is that iterating (12) is similar to a geometric series. Therefore, depending on the spectral radius of the Jacobian, ρ(J(x(t))), it can result in vanishing (ρ < 1) or exploding (ρ > 1) gradients (Zilly et al., 2016).
We can now demonstrate that, by avoiding function compositions, our approach is immune to the exploding gradient problem. In particular, our gradients are fully time-independent and their accuracy is not affected by the sequence length. Recall the ODE:
ẋ(t) = f(x(t), u(t); θ). (13) The following result is obtained:
1We consider a single batch element: f and x are in Rnx×1, where nx is the state dimensionality.
Theorem 1. Assume the full Jacobian of f(x(t), u(t); θ) has a finite spectral radius for all t. Then, the norms of the gradients used in Algorithm 1 and Algorithm 2 are also finite for all t.
We will prove the theorem for Algorithm 2 since it includes the relevant parts of Algorithm 1.
Proof. Given the Legendre basis functions, ψi(t), i = 1, . . . , p, define the ODE residual as:
r(x(t), θ) = p∑ i=0 xiψ̇i(t)− f ( p∑ i=0 xiψi(t), u(t); θ ) ,
for which the residual loss is R(x(t), θ) = ‖r(x(t), θ)‖2. Then, Algorithm 2 consists of the concurrent minimization of the relaxed loss function:∫ t1
t0
γL(x(t))) +R(x(t), θ?)dt,
with respect to the coefficients xi of the Legendre polynomial given θ?, and of the residual loss:∫ t1 t0 R(x(t), θ)dt,
with respect to θ given the set of coefficients x?i , i = 1, . . . , p. For the residual loss gradients are:
∂R(x?(t), θ)
∂θ = −2 r(x?(t), θ)T
∂f ( ∑p i=0 x ? iψi(t), u(t); θ)
∂θ , (14)
where there is no recursion over the previous values x(t − τ), since the basis functions ψi(t) are given and the points x?i are treated as data. By assumption, the Jacobian of f(x(t), u; θ) has finite singular values. Hence, by standard matrix norm identities (see Chapter 5 of Horn & Johnson (2012)) the gradient of R(x(t), θ) with respect to θ has a finite norm for all t.
Similarly, the absence of time recursions in the gradient for the xi update using the relaxed loss also follows trivially from the fact that the coefficients xi are independent variables, that we assume a given θ?, that the basis functions ψi(t) are fixed. Then, the claim follows again from the assumption on the Jacobian of f(x, u; θ).
Note that the result of Theorem 1 cannot easily be achieved when training ODEs using backpropagation through the solver or the adjoint method unless some gradient conditioning, such as clipping or batch normalization (Ioffe & Szegedy, 2015), is performed after applying f . On the other hand, our result relies only on the gradient of f being finite. This is trivial for a shallow f and, for a very deep f , the methods in (Ioffe & Szegedy, 2015; Srivastava et al., 2015; Ciccone et al., 2018; He et al., 2015) can be applied just inside the architecture of f , if 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 to LSTMs (Greff et al., 2017).

Review summaries:
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.