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.
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.
I 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.
Pros of the paper:
	1) Trajectory predictions on the two planar vehicle dynamics experiments was impressive.
	2) The proposed representation of the state transition dynamics is indeed more memory-efficient, and its approximation error is modulatable by the hyperparameter order \p.
	3) Speedup due to parallelization is substantial.
Cons of the paper:
	1) 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.
	2) It slightly detracts from the cleanliness of the story that we must first create an initial trajectory, before performing our coordinate descent.  
Questions and Points of Confusion:
	1) 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?
	2) 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? 
	3) 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.
	4) 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?
	5) 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).
	6) 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.
	7) What motivated the use of the planar vehicle dynamics experiment to showcase your model?  Were there other baselines or benchmarks you considered or attempted?