Accelerating Neural ODEs: A Variational Formulation-based Approach
Authors: Hongjue Zhao
Yuchen Wang
Hairong Qi
Zijie Huang
Han Zhao
Lui Sha
Huajie Shao
Venue: ICLR 2025 Poster (2025)
Keywords: Neural ordinary differential equations
irregularly-sampled dynamical systems
variational formulation
acceleration
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We propose a novel approach based on variational formulation (VF) to accelerate the training of Neural ODEs for dynamical systems.

Abstract: Neural Ordinary Differential Equations (Neural ODEs or NODEs) excel at modeling continuous dynamical systems from observational data, especially when the data is irregularly sampled. However, existing training methods predominantly rely on numerical ODE solvers, which are time-consuming and prone to accumulating numerical errors over time due to autoregression. In this work, we propose VF-NODE, a novel approach based on the variational formulation (VF) to accelerate the training of NODEs. Unlike existing training methods, the proposed VF-NODEs implement a series of global integrals, thus evaluating Deep Neural Network (DNN)--based vector fields only at specific observed data points. This strategy drastically reduces the number of function evaluations (NFEs). Moreover, our method eliminates the use of autoregression, thereby reducing error accumulations for modeling dynamical systems. Nevertheless, the VF loss introduces oscillatory terms into the integrals when using the Fourier basis. We incorporate Filon's method to address this issue. To further enhance the performance for noisy and incomplete data, we employ the natural cubic spline regression to estimate a closed-form approximation. We provide a fundamental analysis of how our approach minimizes computational costs. Extensive experiments demonstrate that our approach accelerates NODE training by 10 to 1000 times compared to existing NODE-based methods, while achieving higher or comparable accuracy in dynamical systems. The code is available at https://github.com/ZhaoHongjue/VF-NODE-ICLR2025.

Published as a conference paper at ICLR 2025
ACCELERATING NEURAL ODES:
A VARIATIONAL
FORMULATION-BASED APPROACH
Hongjue Zhao', Yuchen Wang?, Hairong QF,, Zijie Huang', Han Zhao', Lui Sha', Huajie Shao? University of Illinois Urbana- ChaamaCham,,iin, "iilliam & Mary, "University of Tennessee, Knoxville, Universiity California Los Angeles
{hongjue2, hanzhao, Irs)eillinois edu, {ywang142, hshao wm. edu hqi @utk. edu, zi zijiehuang@cs. wwla.ede
ABSTRACT
Neural Ordinary Differential Equations (Neural ODEs or NODEs) excel at mod- eling continuous dynamical systems from observational data, especially when the data is irregularly sampled. However, existing training methods predomi- nantly rely on numerical ODE solvers, which are time-consuming and prone to accumulating numerical errors over time due to autoregression. In this work, we propose VE-NODE, novel approach based on the variational formulation (VF) to accelerate the training of NODEs. Unlike existing training methods, the proposed VF-NODEs implement series of global integrals, thus evaluat- ing Deep Neural Network (DNN)-based vector fields only at specific observed data points. This strategy drastically reduces the number of function evaluations (NFEs). Moreover, our method eliminates the use of autoregression, thereby re- ducing erтоr accumulations for modeling dynamical systems. Nevertheless, the VF loss introduces oscillatory terms into the integrals when using the Fourier basis. We incorporate Filon's method to address this issue. To further enhance the performance for noisy and incomplete data, we employ the natural cubic spline regression estimate a closed-form approximation. We provide fun- damental analysis of how our approach minimizes computational costs. Exten- sive experiments demonstrate that our approach accelerates NODE training by 10 to 1000 times compared to existing NODE-based methods, while achieving higher or comparable accuracy in dynamical systems. The code available at https: / thub. com hatmm//aaoooogggo jue -
INTRODUCTION
Neural ordinary differential equations (NODEs) (Chen et al., 2018) represent family of continuous- depth machine learning models. Drawing inspiration from ResNets (He et al., 2016), NODEs aim to parameterize vector fields of ODEs using deep neural networks (DNNs),
ft,
where fo x0, Rd is DNN and 0 denotes the model parameters. The continuous nature and specific inductive bias of NODES render them particularly well-suited for modeling dynamical systems from irregularly sampled time series data (Rubanova et al., 2019; Kidger al., 2020). Thus, NODEs have been widely applied to various dynamical system applications, such as multi- agent trajectory forecasting (Wen al., 2022), model-based reinforcement learning (Alvarez al., 2020), optimal control (Chi, 2024) and chemical reaction process modeling (Yin et al., 2023).
In existing training frameworks of NODEs, numerical ODE solvers play crucial role. The forward pass outcomes are directly calculated using numerical ODE solvers. For the backward pass, there are two methods commonly employed to backpropagate through ODE solvers (Kidger, 2022; Onken & Ruthotto, 2020): (1) discretize tibsrreiizethhennptiiniee, which involves directly backpropagating through operations of ODE solvers, and (2) method, utilized in (Chen et al., 2018), which introduces additional adjoint ODEs. In this method, gradients of the scalar loss function with respect to parameters NODEs are computed by solving these adjoint ODEs. More details be found in Appendix B.I.
sewiiizeethe--iiseette also known as the adjoint sensitivity

---

Published as a conference paper at ICLR 2025
(a) DDE-solver-baaed Training Methods
Solution of NODE 3
MSE Loss
Numerical ODE Salver
(R(N)) x(0))
2(0) ODESoweeJee...... a(5a).
Input
(b) VF-NODE
||
J (())
+
I.
folt, a (t))
VF Loss
Sempled Trajectory at
Filor's Method
Filon's Method
Spline Regression
Closed-fore Approximation of Estimated Trajectory
Closed form Approximation of DNN-based Vector Fields
Extra time points required by ODE solver
Sempled time points
2(1)
an 1x)
folt,
it 4)
Spline Interpoiation
Estimated Trajectory
DNN-based Vector Fields fa(t,
DNN fo
Figure 1: Comparison of VF-NODEs and ODE-solver-based methods. (a) ODE-solver-bbsed train- ing methods: ODE solvers compute solutions of the NODEs, necessitating evaluations of the DNN- based vector fields at additional data points. (b) VF-NODEs: implement global integrals numeri- cally using Filon's 's method and spline regression in the VF loss, evaluating the DNN-based vector fields only at specific data points.
Nonetheless, these ODE-solver-based training methods face two significant limitations. First, they are inherently consuminn.. The internal mechanisms within numerical ODE solvers can incur significant computational costs in solving NODEs (Lipman et al., 2022). This is attributed to the numerous evaluations of DNN-based vector fields beyond given sampled data points, as shown in Fig. 1. For the thiiiizetteenddiseretizee approach, this issue may be exacerbated by the introduction of additional adjoint ODEs. Second, existing approaches may suffer from low accuracy. On the one hand, the autoregressive nature of most numerical ODE solvers can lead to error accumulation, as discussed in Appendix B.2. On the other hand, the ceemmiie--ten--diseeee additional numerical discretization error, resulting inaccurate gradients, and potentially causing the training process to fail entirely (Gholami et al., 2019). In short, these limitations stem from the use of numerical ODE solvers.
approach incurs
To date, various approaches have been proposed to address the aforementioned limitations of NODEs. For example, to alleviate the computational bottleneck of NODEs, some works attempt to constrain the complexity of learned dynamics (Finlay et al., 2020; Kelly et al., 2020; Pal et al., 2021), while others seek to directly modify the ODE-solver-basee training process (Kidger et al., 2021: Djeumou et al., 2022: Norcliffe & Deisenroth, 2023; Matei et al., 2023). In the domain of irregularly sampled time series tasks, several models built upon neural differential equations have been proposed (Rubanova et al., 2019; Kidger et al., 2020). Nevertheless, these approaches remain heavily reliant on ODE solvers, preventing them from effectively addressing the computational bot- tleneck.
To address the challenges associated with the reliance on numerical ODE solvers during the training of NODEs, we propose the VF-NODE, novel approach that employs the variational formulation (VF) (Brunel et al., 2014; Hackbusch, 2017) to accelerate training. The proposed approach integrates VF-based loss function within the standard NODE architecture. Our motivation that the VF loss can implement a series of global integrals, thus evaluating the DNN-based vector field only at spe- cific data points in the observations. In contrast, existing ODE-solver- based training methods have to evaluate the DNN-based vector fields at extra data points due to step size settings. Consequently, our method significantly reduces the number of function evaluations (NFEs). On the other hand, these global integrals mitigate autoregression compared with ODE-solver-based training methods, eliminating error accumulation and improving prediction accuracy. However, the utilization of the VF loss poses two challenges for NODEs. (i) It could result in additional oscillatory integrals due to the introduction of Fourier basis functions like sine and cosine. (ii) It does not perform well when data is missing and noisy. To address the first challenge, incorporate Filon's method (Deaño
2

---

Published as a conference paper at ICLR 2025
et al., 2017) for oscillatory integrals. For the second challenge, we adopt the natural cubic spline regression (De Boor, 1978) to obtain a closed-form approximation for Filon's 's method. We also provide a fundamental analysis of the acceleration benefits of our approach in Section 4.3. Exten- sive experiments demonstrate that the VF-NODE achieves a 10 to 1000 times speed increase over existing baselines while maintaining higher comparable accuracy.
Our contributions are three-fold: (1) We introduce a novel approach that employs the variational formulation (VF) to significantly accelerate the training of NODEs. This method significantly re- duces the number of function evaluations (NFEs) and improves the prediction accuracy. (2) We integrate Filon's s method with natural cubic spline regression to effectively compute oscillatory in- tegrals from noisy and partially observed data within the VF loss. (3) Evaluation results on multiple dynamical systems including one real- wollw application demonstrate that our approach achieves 10 to 1000 times acceleration in training speed compared to the baselines while achieving higher or competitive accuracy.
RELATED WORK
NODEs for Dynamical Systems. NODEs have been introduced model continuous dynamical systems. Latent ODEs based on VAEs were first proposed in (Chen et al., 2018), with RNNs as encoders and NODEs as decoders. Rubanova et al. (2019) proposed ODE-RNNS, which use NODEs to simulate continuous dynamics of hidden states in RNNs, and incorporated them into Latent ODEs as encoders. Kidger et al. (2020) introduced Neural Controlled Differential Equations (NCDEs), which can be viewed as the deep limit of RNNs. However, complex mechanisms involved in these methods exacerbate the computational burden during training. In contrast, our approach built on common NODEs, requiring only simple multi-layer perceptrons (MLPs) and the VF loss. It can significantly speed up NODE training and improve the robustness against noisy data. Biloš et al. (2021) introduced Neural Flows alternative to NODEs, which directly use DNNs to model the solution of ODEs. However, they have difficulty modeling autonomous systems.
Acceleration Techniques of NODEs. Existing works accelerating the NODEs training can be di- vided into three categories. (i) Regularization -bssed methods. Some approaches try to constrain the complexity of learned dynamics using regularization techniques, including high-order derivatives regularization (Kelly et al., 2020), kinetic regularization (Finlay et al., 2020), temporal regulariza- tion (Ghosh et al., 2020) and others. Unfortunately. these methods may only accelcrate inference, and the training time may not be reduced. (ii) Design new architectures. Some works directly restrict the dynamics by designing special architectures of NODEs, for example, Heavy Ball NODEs (Xia et al., 2021) and Nesterov NODEs (Nguyen et al., 2022). In addition, model order reduction tech- niques have been utilized to directly compress DNNs in NODEs (Lehtimäki et al., 2022). However, these methods limit the expressivity of learned dynamics, making them unsuitable for dynamical sys- tems. (iii) Modify training process. Other approaches try to modify the training process of NODEs to implement acceleration, such as the IRDM (Daulbaev et al., 2020), seminorn approach (Kidger et al., 202 l), the Taylor theorem (Djeumou et al., 2022), the Gauß Gegenneg quadrature (Norcliffe & Deisenroth, 2023), and the sensitivity- existing modifications still rely on ODE solvers. Thus, the existing computational bottleneck is not tackled effectively, leading to poor acceleration performance. In contrast, our method only performs global numerical integrals, which is much more efficient than existing approaches.
gradient descent (Matei et al., 2023). Nevertheless, all
Variation Formulation of ODES. The VF method, as employed in our work, was initially intro- duced for parameter estimation of known ODEs (Brunel et al., 2014). Subsequent research, such as D-CODE (Qian et al., 2022), extended the use of VF for symbolic regression to facilitate the dis- covery of equations. This approach was primarily adopted to circumvent the need for numerically estimating derivatives from noisy data. However, it was not utilized to expedite the training of mod- els, unlike in our work. Furthermore, the evaluation of D-CODE was limited to dynamical systems sampled at regular intervals, employing the basic compound trapezoidal rule (Press, 2007) for com- puting numerical integrals in the loss function. Notably, D-CODE did not address the challenges posed by oscillatory integrals, which arise with the Fourier basis function (Deaño et al., 2017). In contrast, our application of VF to NODEs serves distinct purposes. Firstly, we leverage VF to accel- erate NODEs training. Secondly, we utilize it to avoid autoregression, thereby enhancing prediction accuracy. From application perspective, our foeus is modeling dynamical systems based on
3

---

Published as a conference paper at ICLR 2025
ineeullaly sampled data. On the technical front, integrate Filon's method to effectively handle oscillatory integrals, a significant advancement over the D-CODE approach.
3 PRELIMINARIES
Problem Setup. Consider time-series dataset in which each trajectory is represented as (((k))),
{(to, stto)), (t1, x(t)))...... (f(,)))) where 0 = to < t < .... < tx T, and x(t+) = (on (1x), (24),),... ..., x,((k)] (IR U {+))² are noisy observations, with denoting possible miss- ing values. Our goal is to speed up the training of NODEs for dynamical systems.
3.1 VARIATIONAL FORMULATION OF ODES
In this subsection, we formally introduce the variational formulation (VF) of ODEs utilized in Theorem 1. Through this formulation, we can establish a direct connection between the trajectory and the vector field f through numerical integral.
Theorem (Variational Formulation of ODEs (Brunel al., 2014; Hackbusch, 2017)). Consider E N+, T E R+, continuous functions : [0, 71 + R f [o, 7]] x IR + R² and ø E Co 71, where C² is the set of continuously differentiable functions. Here we define the functionals
c(æ, f, 0) = [/ a(()((1)d dt + [." f((,x(())(() dt.
(2)
Then æ is the solution to the ODEs f(t, x)) if and only if
c(æ,, f. Φ) = 0. vo E [0 s.t. (00) = ((T) 0
(3)
The proof of Theorem can be found in (Qian et al., 2022). This theorem establishes the necessary and sufficient conditions under which æ the solution the ODEs ă = f(t æ). Moreover, in contrast to existing training frameworks for NODEs, this connection does not rely solving the ODE numerically.
Intuitively, for NODE expressed Eq. (1),
attempt to find the optimal parameters 0 by minimizing ||c(æ, f. 1) Σ (2x f. σ) zero satisfy Eq. (3). To operationalize this concept NODE training, we introduce Theorem 2 (Qian al., 2022), which is outlined below.
Theorem 2.
Let [0. T] Rd be continuously differentiable function which satisfies i 2l, Then for the Lipschitz continuous neural network fo [0, T x Rd + Rd, where are parameters, the following limit holds.
lim Σ So. = [[[o,, - 1z) o (l, x) L-000
(4)
Fml
= where {0), d2...... , )r) are series of Hilbert orthonormal basis for L2 7] such that (x(0) (e(T) = 0, l I,...,L and [[[fe fz) (t, x) = fo" (t, x(t)) x(t))) dt.
The proof of Theorem 2 can be found in Appendix C. In general, Theorem 2 enables the appli- cation of the VF to the training NODEs. By leveraging this theorem, we convert an infinite number of constraints into series of constraints using set of orthogonal basis functions. Addi- tionally, we establish a connection between the VF and distances between functions, implying that as || fo, 1 converges to zero, the parameters in NODEs also converge to their optimal values.
3.2 FILON'S METHOD
Filon's method is designed for oscillatory integrals. Consider the integral 1 sin((t))))) de. We aim to compute this integral numerically using the available data points (h(tz) I with 0 = to < lK T.
4

---

Published as a conference paper at ICLR 2025
For general numerical integration techniques, such the Newton-Cotes formula (Press, 2007), they utilize polynomials (t)(t to approximate the entire integrand h(t) sin(wst within each time inter- val [. Consequently, these methods struggle to handle high-frequency oscillations increases, because polynomials cannot effectively match the high-order derivatives of the integrand. Further analysis on this issue is provided in Appendix D. To address this limitation, introduce Filon's method (Deaño et al., 2017), which employs polynomials solely to approximate the non- oscillatory component h(t) in the original integrand. This approach enables the estimation of the original integral as follows:
J h(t) sin(wt) dt
K-1 Σ 1m
h(t) sin(wi) dt R
Σ
qu (in(ww) dt,
(5)
where (8) is the approximate ne-th-order polynomial of h(t) over the time interval [tx. fR+1]. is important to note that the integral fu++ (t) (wt) dt can always be computed analytically. This approach ensures that the accuracy of numerical integrals is not affected by w, effectively addressing the challenge of oscillatory integrals.
4 PROPOSED METHOD
In this section, we begin by introducing the VF loss for NODEs based on Theorem 2 and addressing the computational challenges in the VF loss using Filon's 's method and natural cubic spline regression. Then, we elaborate on the steps of VF-NODEs. Finally, we provide fundamental analysis of the acceleration benefits of our approach.
4.1 COMPUTING THE VF LOSS
In Section 3.1, we systematically introduced the VF utilized in this work, and demonstrated the connection between this formulation and the distance between the vector field fe of the NODE and the ground truth vector field f. Based on Theorem and Theorem 2, the optimal parameters 0 of NODE, as expressed in Eq. (1), can be obtained by solving the following optimization problem:
N
ΣΣ cccoll SS, So. 1
or arg min
(6)
c(x11) 0z) h' Ll (t) S( (t))
where (x11)] ]==1 are trajectories in the dataset, and (x = 22/T sin(Tet/T)) which can be consid- ered as the sine Fourier basis. However, computing the VF loss poses two computational challenges. (i) The use of the sine Fourier basis introduces oscillatory terms into the integrals. (ii) It is challeng- ing to compute integrals numerically from noisy and incomplete data.
To address the first challenge, we introduce the Filon's method, which is designed for oscillatory integrals. Integrals in Eq. (6) are one-dimensional oscillatory integrals, and these integrals can be concluded as 1 sin(t) dl, without loss of generality. As discussed in Section 3.2, polynomi- als struggle to capture the nature of the oscillatory integrand f((t) sin(wt) as w increases. To tackle this issue, we leverage Filon's method, which uses polynomials to approximate the non-oscillatory part h(t)) in each interval. By focusing on h(t) alone, we eliminate the influence w on the precision of the closed-form approximation.
To deal with the second challenge, we introduce the natural cubic spline regression (De Boor, 1978), which is used build up a precise closed-form polynomial approximation for h(()) in Filon's method from noisy data. For ) with 0 = to < 4 < ... < lx = T, we aim to construct cubic polynomial (t) in each interval [fx, 1x in Eq. (5):
3
(t( a,.0 + Σ an.m (t t² = a + a (x) + ax.2(t - tx)² + ax,((t - tx) (7)
m==1
By meeting the requirements for continuity, differentiabilityy natural boundary conditions, and smoothness, we calculate the coefficients 4g.m (k O,..., K 1,m 2,33 in Eq. (7). In
5

---

Published as a conference paper at ICLR 2025
general, spline regression allows create sufficiently precise closed-form approximations while maintaining smoothness. We provide a detailed process for this calculation in Appendix E.
In summary, to compute the integral Jo ((t) sin dt, we first construct the natural cubic spline approximation of h((), denoted as T& (t), within each time interval using sampled data {(t,, {(tx)))k= k== By substituting Eq. (7) into Eq. (5), we can estimate the integral accurately by computing Dk=0 -1 ptata (1(() sin(wt) dl.
4.2 VE-NODES
The proposed VF-NODE uses the same neural architecture as the vanilla NODE, but it is trained based on the VF loss. Given trajectory {(t,, x(t))) where 0 to t tk T, present the detailed steps of our method.
Step 1: Perform natural cubic spline regression on æ((2) to get the spline coefficients ak,m.
Step Estimate the trajectory based on the spline to remove noise: x((x) ai50 +
En at.m(ts - -x ax.o (k = O,..,, 1) and @(fx) = Σm ak ax- - 1x-1))
Step 3: Evaluate the vector fields (l (((4)), k O,...,.. K.
Step 4: Perform natural cubic spline interpolation on fo(!,, x(4)) to get the spline coefficients
by.m
Step 5: Compute the VF loss ||c(æ, fo. 2 based on ak.m and by.m=
/'/ a a(t) dt S/ folt. fo(() dt
Σ cl for 12 Σ
(8)
where
S' a(t)de(t) dt
2 al Σ
3
tx+s Clett (t (t-tz) -xxx cos
v
akmm
dt
T
T
T
km0 m=0
(9)
Ld folt, ((1,)())))(()) R0
2 K
3
totl
Tlt dt.
ba.m
(1 (1)
sin
T
k=0 m=0
These steps summarized as Algorithm in Appendix F. In this process, the natural cubic spline serves two key roles in our method. (i) It builds up precise closed-form approximations the trajec- tory at), denoted as (t), and vector fields ft, att)) from noisy observations, which is essential for estimating the oscillatory integral accurately. (ii) It fills missing values in the trajectory att), which is necessary for evaluating the vector fields (t))
Remark It is important emphasize that our primary use of natural cubic spline is for construct- ing the closed-fome form approximation, which is essential for computing numerical integrals, rather than for transforming the original irregularly sampled data into regularly spaced data. While we do uti- lize this technique to fill in missing values, it is natural byproduct of constructing the closed-form approximation.
Remark 2. Because we have removed noise in observations by performing spline regression on x<(), only need perform spline interpolation on fo (+)) for obtaining closed-form approxima- tion.
4.3 FUNDAMENTAL ANALYSIS FOR THE ACCELERATION OF VF-NODES
Now we demonstrate how our method accelerates the training of NODES in two ways: (1) reduce the number of function evaluations (NFEs) and (2) improve parallelizability Consider a trajectory with M observed data points. As shown in Fig. (a), ODE-solver-based training methods require evalu- ating vector fields fo many more than M times. For example, common staptive-step-szee explicit Runge-Kutte Kutta methods (Butcher, 2016) an autoregressive formula ...1 En h E/=1 2j9; to estimate the solution of NODE expressed in Eq. (1). Here, h is the step size, J is the order

---

Published as a conference paper at ICLR 2025
of the ODE solver, 2j coefficients and 9, are evaluations of the DNN-based vector field fø- In this case, to make seeesste-forwaal prediction, the DNN-based vector field must be be evaluated for J times. Moreover, to ensure accuracy, these ODE solvers often require evaluating the DNN- based vector field at additional data points beyond those provided in the observations. As a result, the vector field fo will be evaluated > M J times using ODE-solver-basee training methods. Additionally, due to the autoregressive nature of ODE solvers, these vector fields must be evaluated step by step.
However, VF-NODES only evaluate the DNN-based vector fields for exactly M times, as shown in Fig. 1 (b). As discussed in Section 4.2, to compute global integrals in the VF loss numerically, we evaluate the vector fields fe only at the specific sampled data points in Step 3, in order to construct a closed-form spline approximation in Step 4. In addition, these vector fields can be evaluated simultaneously in Step 3. Although I. integrals need to be computed in c(æ: fo, Φb)) for i 1,...: L, these integrals share the same vector fields, allowing for efficient computation.
Summary of the Proposed VF-NODES. Our method enables the learning of parameters in NODEs without relying on numerical ODE solvers. Regarding computational efficiency, our approach re- quires only a series of global integrals for each trajectory. In these numerical integrals, the vector fields of NODEs are only evaluated at observed data points in parallel, significantly reducing the number of function evaluations and achieving better parallelizability. Additionally, our approach effectively mitigates error accumulation from autoregression, thus improving prediction accuracy.
5 EXPERIMENTS
In this section, we conduct a series of experiments to evaluate the performance of VF-NODEs. We begin by applying the method to four dynamical systems of various dimensions. Subsequently, we assess its performance on the real-world COVID-19 dataset. Lastly, we perform ablation studies to examine the model's 's key components.
Baselines. To evaluate acceleration performance, we compare our method against the following training approaches for NODE-based models: (1) the discretine-thennoptiiize approach (Dis-Opt), (2) the pptimize-then-discetiin approach (Opt-Dis) (Chen et al., 2018), and (3) the seminorm ap- proach (Kidger et al., 2021). Additionally, we compare the accuracy of our proposed method against two categories of SOTA methods. The first category includes NODE-based models: (I) Vanilla NODE (Chen et al., 2018), (2) TayNODE (Kelly et al., 2020), (3) Latent ODE with RNN encoder (Chen et al., 2018), 8), (4) ODE-RN BON- (Rubanova al., 2019), (5) Latent ODE with ODE- RNN encoder (Rubanova et al., 2019), and (6) NCDE (Kidger et al., 2020). The second category consists of Neural Flows (Biloš al., 2021), where neural networks go used directly model the solution of ODEs as x = (1, x.). This category includes models such as ResNet Flows and GRU Flows. Detailed hyperparameter settings of these models are provided in Appendix I.1. We also discuss other training methods for NODEs based on spline mtchods, as shown in Appendix H.
5.1 SIMULATION OF LOW-DIMENSIONA DYNAMICAL SYSTEMS
Dynamical Systems. We select four dynamical systems from fields such as biology, biochem- istry, genetics and epidemiology: the glycolytic oscillator (Sel'Kov, 1968), the genetic toggle switch (Gardner et al., 2000), the repressilator (Elowitz & Leibler, 2000), and the age-structured SIR model (Ram & Schaposnik, 2021). These systems vary in dimensionality; the repressilator is six dixddmmesioonal, and the age-structured SIR model is -imeennnii, while the remaining systems are two-dimensional. Further dataset details are provided in Appendix 1.2.
Datasets. For each system, except the age-structured SIR model, we generate 125 trajectories with randomly sampled initial conditions, of which 100 are used for training and 25 for validation. Each trajectory is uniformly sampled with 100 data points from UO,, with values randomly dropped at rate of r via masking. Following (Qian et al., 2022), Gaussian noise e = OR std(w(t)) is added to each observation. For baselines, each trajectory is segmented into 5 short segments to enhance the performance of baselines, while for VF-NODEs, the whole trajectory is used to improve the precision of spline regression. Additionally, we generate 25 test trajectories, each containing 200 randomly sampled points. The first 100 points are sampled from [0, T] for interpolation tasks, and the remaining 100 from UTT,20 2T) for extrapolation tasks, In experiments, T is set to 10, an to
7

---

Published as a conference paper at ICLR 2025
0.01, and the retention ratio model except for the number of trajectories: 500 are generated for training and validation (400 for training and 100 for validation), and additional 100 trajectories generated for testing.
0.8. All settings remain the same for the age-structured SIR
Evaluation on Acceleration. We compare the acceleration performance of our method against NODE-based baselines in terms of average training time per epoch using the glycolytic oscillator. As shown Fig. 2, our method significantly accelerates the training of NODE-based baselines by factor of 10 1000. This improvement is attributed to our approach, which relies solely on global numerical integration over the time domain, thereby reducing NFEs during training. Average training time per epoch on other dynamical systems are provided in Appendix G.I.
Avg Training Time Epoch (s)
10
VF-NODEs
Dis-Opt
Opt-Dis
Seminorm
10
10-1
VF-NODEs
TayNODE
NCDE
Figure 2: Average training time per epoch (second) for each method on the glycolytie model. Our method can achieve 10 to 1000 times faster than the baselines. Due to the high stability of the training speed for some methods, the uncertainty is negligible and not clearly visible in the figure.
Evaluation Prediction Error. We evaluate the prediction error of different methods in- terpolation and extrapolation tasks using mean squared error (MSE), as shown in Tables 1 and respectively. The results demonstrate that method consistently achieves higher or comparable performance compared wth the baselines. This superior performance can be attributed to two key factors: (I) the standard NODE architecture in VF-NODES aligns well with dynamical systems represented by ODEs, and (2) the use of global integrals in the VF loss mitigates autoregression, ef- fectively reducing error accumulation. Although Vanilla NODEs share the same neural architecture, they rely on ODE solvers for training, which leads to error accumulation and significantly slower training speeds compared to VF-NODEs. the case of TayNODEs, the regularization term restricts the complexity of the learned dynamics, limiting the expressivity of the model and leading to even worse performance than Vanilla NODEs. Neural Flows also perform worse than VF-NODEs, as they must explicitly model the time dependencies of solutions. In contrast, VF-NODEs can leverage ODE solvers to handle time evolution during inference. Additionally, other baselines are designed for more general irregular time series data. They do not align with the form of these dynamical systems. Detailed training settings can found in Appendix I.3.
Additionally, test method the Gompertz model (Gompertz, 1825) and the Lotka- Volteraa equations (Kingsland, 1995). Related experimental results can be found in Appendix G.2.
Table I: Testing MSE (meantstandard deviation) for interpolation tasks on dynamical systems with 80% observed data (r 0.8). Lower values indicate better performance. Here ein refers to x x104m Latent ODE refers Latent ODE with an RNN encoder. Latent ODE 2 refers to Latent ODE with an ODE-RNN encoder. The best results are highlighted in bold black, and the second- best results are highlighted in bold purple.
Glycolytic (1.51c-03)s(1.44eee) 06.66-4)))
Toggle
Repressilator
AgeSIR
Vanilla NODE
(2.25-002:+(5046-003
(7.54-00))=(6.556-44)
TayNODE
(3.20=-03)4(1.32e-033 (1.37e-02)4
(I.4R=-0))((..22-022
(3. (3.18e-0))((7.36--22
Latent ODE
(2.21e-01)(..33--00) (6.57c-01)+(1448--01
(2.336)0))99226000))
(4.18++01)(6655==000))
ODE-RNN
8.800-495)+(2330-055 (1.63e-03)H(5555eee))
(1.41e-1))(7.99e-033)
(6.7 (6.75e+00
Latent ODE 2
(1000=01)=(1.877000) (6.31c-01)+(4.66o-00)
(I336c+09 (1:52=00))(1222++000)
(7.47c-0))(8.776-002)
NCDE
(3.49c-02)4(1336ee) (3.96e-02)((333--00000
(1.22e+01 (43775e=00)
ResNet Flow
(2.84e-01)4
(7.09e01))42.21e-0111
(1.03=00))H(814eel))))
(.966)))))))
GRU Flow
(3.80=-0))((.34ee00
(2 (.45==00)((1888--0)))
(7445e000)144686-00)
(4. 10119e011)(1.22200
VT-NODE (Ours)
(6.35--05)4(2666-066
(L..9-00)):660006-05)
(1.92-02))(2.22-044)
(7.39-00))
8

---

Published as a conference paper at ICLR 2025
Table 2: Testing MSE (meandsstandard deviation) for extrapolation task 4 dynamical systems with 80% observed data (r = 0.8). Lower values indicate better performance. Here ein refers to =104.. Latent ODE refers Latent ODE with an RNN encoder. Latent ODE 2 refers to Latent ODE with an ODE-RNN encoder. The best results are highlighted in bold black, and the second- best results are highlighted in bold purple.
Glycolytic
Toggle
Repressilator
AgeSIR
Vanilla NODE
(8.79-04))((76444440 (8.14--07)+(6722--07)
(1.25x-01)=(3311e-22)
(1.99-02)=<(1.66ee
TayNODE
(4.71e-03))(3.600-00 (5.09e-02))(5009--022
(8.73=-01)=(1335-001
(4.52--0))((331e-0))
Latent ODE
(2.39e-01)a(7229e-0)) (1.48=+00)t(11288)))))))
(9. 9.11==00)==(119<+000)
(2.0e=00)x((.85ee+0
ODE-RNN
(4.68c-05)+(11115-555)
(204e-04)=(1142e2-))
(7 (7.48c+00)=(). (2.02--0))(7.51e-003)
Latent ODE
(4.96=+00)+(5700:0)) (1.6e-01)=(1.6ee-00)
(1.36c+09)s(11933))) (3099=000)((2990000)
NCDE
(.000=00))(44448-000
(1 (1.00e+0)=((738---0000
(6.73=+00
(2. (tbe=00))((8227-000000
ResNet Flow
(3:47c=000)((.822-0000
(6.32=000)111177-000
(6.56=+01)#(2223+00)
(L.95=+01)=(3329c))
GRU Flow
(7.39--0))4(2223-000)
(5002e+00)+(5222-001
(1 (1.84e+01)4(5744-00)
(6.02e=01)a( (a110090))(1899-011
VF NODE (Ours)
(1.33e-0))a(3355--05)
(4.79--0))+(5.24--88)
(1.23e-01)x(1444-00)
(2.377-02)=
5.2 REAL-WORLD APPLICATION: COVID-19 DATASET
We also evaluate the performance of VF-NODEs on the real-world COVID-19 dataset. Us- ing covsirphy (Takaya & Team, 2020- (024), we leverage data from the COVID-19 Data Hub (Guidotti & Ardia, 2020). This data can be modeled using four- dinensinnii phase-dependent SIR-F framework (Takaya & Team, 2020-2024), where the ODE parameters vary over time. VF- NODEs, however, are not designed to handle these time-varying dynamics due to the limitations of the standard NODE architecture. To address this, we apply S-R analysis (Balkew, 2010) to segment the data, allowing the ODE parameters to remain constant within each segment. In this study, we evaluate our method and the baselines using data from four countries: Japan, Italy, Norway, and India. For each country, extract 100 data points from the longest segment in single trajectory, allocating 80 for training, 10 for validation, and 10 for testing. Each data point consists of four variables: Susceptible, Infected, Recovered, and Fatal. We standardize the data for consistency and employ spline regression to smooth it. The experimental results, presented in Table 3, demonstrate that our method achieves significantly better performance compared the baselines. This superior performance can be attributed to the fact that the COVID-19 dataset conforms to the SIR-F ODE system, which aligns with the inductive bias of VF-NODEs. Vanilla NODEs outperform most other baselines because they share the same neural architecture as VF-NODEs. However, Vanilla NODEs have lower accuracy than our method due to the error accumulation caused by autoregression in ODE-solver-base training methods. Detailed training settings are provided in Appendix I.3.
Table 3: Testing MSE (meantstandard deviation) on the real-world COVID-19 dataset. Lower values indicate better performance. Here etn refers to x 110#1 Latent ODE refers to Latent ODE with an RNN encoder. Latent ODE refers to Latent ODE with an ODE-RNN encoder. The best results are highlighted in bold black.
Japan
Italy
Norway
India
Vanilla NODE
(1442==00::(744-000 (L.336-02))(13336000)
(1003-033)+66330-05))
(8.866-01)=(33766-00
TayNODE
(2002=+00=(88000-00) (3886-0)))457700033))
(1.18-02))(11088-4)))
Latent ODE 1
(1042-00)) (6.56e-01)(228876001
(1883-01)=((1100000)
(5.69--0))+1.559-000)
ODE-RNN
(84436000) (1.10e-01)+<1.122) 12e-02)
(5.58e003y<(1.44e-00)
(.09e-01)s(11744e4)))
Latent ODE 2
(1133+0))*(1.31c+00)) (3.36e-01 (227000))
(d.122-01)+(2990-000)
(4.07c-01)a( (51.e6000)
NCDE
(8.72c001)=(2446:44 (.0444))))((.100=000)
(1.27=00))((.200-00))
ResNet Flow
(9.330-0))=(16666-00) (3.696-02)+(3.200000))
(2.51c-02)+(1655--02)
(2.04o02)
GRU Flow
(1.39e+000)(2296-0002
(3.07e-03)=(3076-07) (8.63-003)((111e-044))
(2.52e-03)=(1884ee00))
VF NODE (Ours)
(1.376-01)+(4822-002) (16.4:-43))(.111-004)
(3.43--44)=11118-044)
(5.686-04))(2288-00))
Furthermore, we also apply the proposed VF-NODEs to model the temporal effects of chemotherapy on tumor volume. The experimental results are presented in Appendix G.3.

---

Published as a conference paper at ICLR 2025
5.3 ABLATION STUDIES
Lastly we conduct ablation studies to investigate the impact of various components on the perfor- mance of VF-NODES. Specifically, we examine: (I) the influence of Filon's method and spline regression and (2) the type of basis functions. The experimental settings are consistent with those described in Section 5.1.
Influence of Filon's method and spline regression. There are two key components enabling the accurate computation of the VF loss in VF-NODES: (1) Filon's 's method and (2) natural cubic spline regression. To evaluate the influence of these components, we assess the performance of VF-NODEs under the following settings: (1) using only Filon's method (Filon's 's method + natural cubic spline interpolation), (2) using only natural cubic spline regression (trapezoidal rule (Press, 2007) + natural cubic spline regression), (3) using neither (trapezoidal rule + natural cubic spline interpolation). (4) using Filon's method with Hermite cubic spline. The experimental results, shown in Table 4, demonstrate that while both components enhance performance, natural cubic spline regression has more significant impact. The inferior performance of Hermite cubic spline can be attributed to its reliance on numerical differentiation, which is ill- ill-posed on noisy, sparse data.
Table Testing MSE (meantstandard deviation) for the ablation study on Filon's method and spline regression. Lower values indicate better performance. Here een refers to x x10
Glycolytic
Toggle
Filon
Spline
Regression
Interpolation
Extrapolation
Interpolation
Extrapolation
Natural
1e++03)e(11 29c+666
(2. 3.55++00)
(9. 160
(8. (6e+37)=(l
Natural
(5. (5.76e=02)#((56660
6774002
(1.40c+00)=(1400c)
(4.35e+20)=44.366)
X Natural
(5.59c=01)#654440)) (8++++011)
(..60=011++66555000
(1 (0000e111((1.0eee
/ Hermite
(1. (1550c000)422.7c=000) (3.90=-01)+(555e+000
(3.74k-01)=(1.14cc
(9 60.66-001
VE.NODE (Ours)
(6.35--05)4(2668-666)
(1635-00))
(1.99-04)):66009-053
(4.79e-77)=(5.24e-08)
Type of basis functions. To illustrate the rationale behind using Fourier basis and Filon's 'ss method,
we test two alternative types of basis functions: (1) (be(t) (t/T)(t/T 1)(t/T - e/ L). and (2) Hermite cubic spline interpolation
three points: (0, 0), (T/2, and (T, 0). with derivatives at these points set to 0. Filon's 's method was not applied to these basis functions. The experimental results, shown in Table 5, indicate that VF-NODES with Fourier basis significantly outperform those with other basis functions. This can be attributed to the requirement that de should form a Hilbert orthonormal basis in the function space. as per Theorem 2. However, performing orthonormalization on these basis functions is computationally expensive. Other Hilbert orthonor- mal basis functions, such as orthonormal polynomials Chihara (2011), also exist but fail to meet the boundary conditions (0) (T) O. While Qian et al. (2022) used B-splines, their recursive generation (De Boor, 1978) can significantly slow down the training VF-NODEs.
Table 5: Testing MSE (meantstandard deviation) for the ablation study on the types for basis func- tions dr. Lower values indicate better performance. Here etn refers to x 10
Glycolytic
Toggle
Interpolation
Extrapolation
Interpolation
Extrapolation
1)U/ -// (/L)
(3.89e-02)448B37cc)) (L555-02)
(3.44e-01)a(3399ee))) 79c-001
(6.23--01
huudi basis
(7.06e-01)4(5.3966 (01)))) (5.296-01)4(4.99-00)
(4. (116==00)444555-000
(1 (1,62=00))(2266+000
VF-NODE (Ours)
(6.35--05)1/226880000) (1.636-04))(33355-05)
(1.69-0))))6
(4.79:-07::(5246668))
We also conduct some additional ablation studies, including the sensitivity of VF-NODEs to hyper-
parameters, the effect of sampling, and the effect of noise, which can be found in Appendix G.4.
CONCLUSION
This work introduced a novel training method based on variational formulation to accelerate the training of NODEs for dynamical systems. Our approach required only a series of global integrals in the loss computation, eliminating the need for traditional numerical ODE solvers. To overcome the challenges posed by oscillatory integrals in the VF loss, we incorporated Filon's 's method to en- hance model performance. Additionally, we developed a natural cubic spline regression to better handle noisy and incomplete data. Extensive experiments demonstrated that our method signifi- cantly accelerates NODÉs training while maintaining high accuracy. Limitations are discussed in Appendix A.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
Research reported in this paper was sponsored in part by NSF CPS 2311086, NSF CIRC 716152, and Faculty Research Grant at William Mary 141446.
ETHICS STATEMENT
Our research centered on methodological advancements in NODEs, specifically aiming to enhance the efficient modeling of dynamical systems in fields such biology, epidemiology. and genetics. This work does not involve human subjects, personal data, or any form of invasive technology. All datasets used in our experiments either generated from known ODEs sourced from publicly available datasets, ensuring compliance with their respective licenses and usage policies.
REPRODUCIBILITY STATEMENT
We are committed to promoting reproducibility in scientific research. To facilitate this, we have pro- vided comprehensive details of our experimental setups in Section and Appendix I. This includes thorough descriptions of the datasets used, model architectures, hyperparameter settings, training procedures, and evaluation metrics. We will also release our source code for reproducing all experi- ments public repository upon publication.
REFERENCES
Victor M Martinez Alvarez, Rareş Rosca, and Cristian G Fälcujescu. Dynode: Neural or- dinary differential equations for dynamics modeling in continuous control. arXiv preprint arXiv:2009.04278 2020.
Teshome Mogessic Balkew. The SIR model when S (t) is A multi-exponential function. functioon East Ten- nessee State University, 2010.
Marin Biloš, Johanna Sommer, Syama Sundar Rangapuram, Tim Januschowski, and Stephan Günnemann. Neural flows: Efficient alternative to neural odes. Advances in neural informa- tion processing systems, 34:21325-21337, 2021.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL http:
//github. com/ jax 1aaamxll///1 mll /jax
Nicolas J.-B. Brunel, Quentin Clairon, and Florence d' Alche-Buc. Parametric Estimation of Ordi- nary Differential Equations with Orthogonality Conditions. Journal of the American Statistical Association, 109(505):133-185, 2014.
John Charles Butcher. Numerical methods for ordinary differential equations. John Wiley & Sons, 2016.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Cheng Chi. Nodec: Neural ode for optimal control of unknown dynamical systems. arXiv preprint arXiv:2401.01836, 2024.
Theodore Chihara. An introduction orthogonal polynomials. Courier Corporation, 2011.
Talgat Daulbaev, Alexandr Katrutsa, Larisa Markeeva, Julia Gusak, Andrzej Cichocki, and Ivan Oseledets. Interpolation technique to speed up gradients propagation in neural odes. Advances in Neural Information Processing Systems, 33:16689-16700, 2020.
Carl De Boor. practical guide to splines, volume 27. springer verlag New York, 1978.

---

Published as a conference paper at ICLR 2025
Alfredo Deaño, Daan Huybrechs, and Arieh Iserles. Computing highty oscillatory integrals. SIAM, 2017.
DeepMind, Igor Babuschkin, Kate Baumli, Alison Bell, Surya Bhupatiraju, Jake Bruce, Peter Buchlovsky, David Budden, Trevor Cai, Aidan Clark, Ivo Danihelka, Antoine Dedieu, Clau- dio Fantacci, Jonathan Godwin, Chris Jones, Ross Hemsley, Tom Hennigan, Matteo Hessel, Shaobo Hou, Steven Kapturowski, Thomas Keck, Iurii Kemaev, Michael King, Markus Kunesch, Lena Martens, Hamza Merzic, Vladimir Mikulik, Tamara Norman, George Papamakarios, John Quan, Roman Ring, Francisco Ruiz, Alvaro Sanchez, Laurent Sartran, Rosalia Schneider, Eren Sezener, Stephen Spencer, Srivatsan Srinivasan, Milovs Stanojevic, Wojeiech Stokowiec, Luyu Wang, Guangyao Zhou, and Fabio Viola. The DeepMind JAX Ecosystem, 2020. URL http:
//github. com/ google- delle--deeinddd
Franck Djeumou, Cyrus Neary, Eric Goubault, Sylvie Putot, and Ufuk Topeu. Taylan- neural ordinary differential equations: Toward fast training and evaluation of neural odes. arXiv preprint ar arXiv:2201.0577 2022.
Michael B Elowitz and Stanislas Leibler. A synthetic oscillatory network of transcriptional regula- tors. Nature, 403(6767) 3338338, 2000.
Chris Finlay, Jörn-Henrik Jacobsen, Levon Nurbekyan, and Adam Oberman. How to train your neural ode: the world of jacobian and kinetic regularization. In International conference machine learning, pp. 31 1164. PMLR, 2020.
Timothy S Gardner, Charles R Cantor, and James J Collins. Construction of genetic toggle switch in escherichia coli. Nature, 303(6767) 2000.
Amir Gholami, Kurt Keutzer, and George Biros. Anode: Unconditionally accurate memory-efficient gradients for neural odes. arXiv preprint arXiv:1902.10200 2019.
Arnab Ghosh, Harkirat Behl, Emilien Dupont, Philip Torr, and Vinay Namboodiri. Steer: Simple temporal regularization for neural ode. Advances Neural Information Processing Systems, 33: 14831-14843, 2020.
Benjamin Gompertz. Xxiv. on the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. in letter to francis baily, esq. frs &c. Philosophical transactions the Royal Society of London, (115):513-583, 1825.
Emanuele Guidotti and David Ardia. Covid-19 data hub. Journal of Open Source Software, 5(51): 2376, 2020. doi: 10.21105/joss.023777 023766
Wolfgang Hackbusch. Variational formulation. Elliptic Differential Equations: Theory and Numer- ical Treatment, pp. 159-180, 2017.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp- 770-778, 2016.
Jacob Kelly, Jesse Bettencourt, Matthew Johnson, and David K Duvenaud. Learning differential equations that are easy to solve. Advances in Neural Information Processing Systems, 33:4370- 4380, 2020.
Patrick Kidger. On neural differential equations. arXiv preprint arXiv:2202.02433. 2022.
Patrick Kidger and Cristian Garcia. Equinox: neural networks in JAX via callable PyTrees and filtered transformations. Differentiable Programming workshop at Neural Information Processing Systems 2021, 2021.
Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equa- tions for irregular time series. Advances in Neural Information Processing Systems, 33:6696- 6707, 2020.
Patrick Kidger, Ricky TQ Chen, and Terry J Lyons. " hey, that's not an ode": Faster ode adjoints via seminorms. In ICML., pp. 5443-5452, 2021.
12

---

Published as a conference paper at ICLR 2025
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
Sharon Kingsland. Modeling nature. University of Chicago Press, 1995.
Ian Knowles and Robert Renka. Methods for numerical differentiation of noisy data. Electron, J. Differ. Equ, 21:235-2466 2014.
Mikko Lehtimäki, Lassi Paunonen, and Marja-Leena Linne. Accelerating neural odes using model order reduction. IEEE Transactions on Neural Networks and Learning Systems, 2022.
Xin Li, Jingdong Zhang, Qunxi Zhu, Chengli Zhao, Xue Zhang, Xiaojun Duan, and Wei Lin. From fourier to neural odes: Flow matching for modeling complex systems. arXiv preprint arXiv:2405. 1542, 2024.
Yaron Lipman, Ricky TQ Chen, Heli Ben- Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint arXiv:2210.02747, 2022.
Ion Matei, Maksym Zhenirovskyy, Johan De Kleer, and John Maxwell. Sensitivity-free gradient descent algorithms. Journal of Muchine Learning Reseurch, 24(300):1-26, 2023.
Ho Huu Nghia Nguyen, Tan Nguyen, Huyen Vo, Stanley Osher, and Thieu Vo. Improving neural ordinary differential equations with nesterov's accelerated gradient method. Advances in Neural Information Processing Systems, 35:7712-7722,, 2022.
Alexander Norcliffe and Mare Peter Deisenroth. Faster training of neural odes using sauB-legendre quadrature. arXiv preprint 2rXiv:2308.10644 2023.
Derek Onken and Lars Ruthotto. Discretize-optimize vs. optimize-discretize for time-series regres- sion and continuous normalizing flows. arXiv preprint arXiv:2005. 13420, 2020.
Avik Pal, Yingbo Ma, Viral Shah, and Christopher Rackauckas. Opening the blackbox: Accel- erating neural differential equations by regularizing internal solver heuristics. In International Conference on Machine Learning, pp. 3325-8335. PMLR, 2021.
William H Press. Numerical recipes 3rd edition: The art of scientific computing. Cambridge uni- versity press, 2007.
Zhaozhi Qian, Krzysztof Kacprzyk, and Mihaela van der Schaar. D-code: Discovering closed-form odes from observed trajectories. In International Conference on Learning Representations, 2022.
Vishaal Ram and Laura Schaposnik. A modified age-structured sir model for covid-19 type viruses. Scientific reports, 11(1):1519,, 2021.
Yulia Rubanova, Ricky TQ Chen, and David K Duvenaud. Latent ordinary differential equations for irreuulary-samppled time series. Advances in neural information processing systems, 32, 2019.
EE Kov. Self-oscillations in glycolysis simple kinetic model. European Journal of Bio- chemistry, 4(1):79-86, 1968.
Hirokazu Takaya and CovsirPhy Development Team. Covsirphy version 3.1.2: Python library for covid-19 analysis with phase-dependent sir-derived ode models. https://gi com/ lisphi coo/vovviii// 9-sir, 2020-2024 Accessed: 2024-09-19.
James M Varah. A spline least squares method for numerical parameter estimation in differential equations. SIAM Journal on Scientific and Statistical Computing, 3(1):28 46, 1982.
Pauli Virtanen, Ralf Gommers, Travis F. Oliphant, Matt Haberland, Tyler Reddy, David Courna- peau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. Nel- son, Eric Jones, Robert Kern, Eric Larson, C Carey, ilhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mul- bregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272, 2020. doi: 10. 1038/03//411 592-019-0686--22
13

---

Published as a conference paper at ICLR 2025
Song Wen, Hao Wang, and Dimitris Metaxas. Social ode: Multi-agent trajectory forecasting with neural ordinary differential equations. In European Conference on Computer Vision, pp. 217-233. Springer, 2022.
Hedi Xia, Vai Suliafu, Hangjie Ji, Tan Nguyen, Andrea Bertozzi, Stanley Osher, and Bao Wang. Heavy ball neural ordinary differential equations. Advances in Neural Information Processing 2021.
Systems,
Jun Yin, Jiali Li, Iftekhar Karimi, and Xiaonan Wang. Generalized reactor neural ODE for dy- namic reaction process modeling with physical interpretability. Chemical Engineering Journal, 132:13 2023.
Qunxi Zhu, Xin Li, and Wei Lin. Leveraging neural differential equations and adaptive delayed feedback to detect unstable periodic orbits based irregularly sampled time series. Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(3), 2023.
14

---

Published as a conference paper at ICLR 2025
A LIMITATIONS
limitation of our approach is its reliance on polynomial-based approximation, which may not ef- fectively approximate complex trajectories intricate dynamical systems, such as chaotic systems, or may fail when the sampling is extremely sparse. However, the proposed VF-NODEs can still cap- ture underlying patterns of chaotic systems, as shown in Appendix Additionally, computing the numerical integrals requires the trajectory points are available, meaning VF-NODEs are currently applicable only time series tasks. In future work, we aim explore new techniques address these challenges, such as coordinate gradient descent (Matei et al., 2023).
B BACKGROUND: TRAINING OF NODES
B.1 FOUNDATIONAL. TRAINING FRAMEWORKS OF NODES
In current training frameworks of NODEs, the outcomes are directly computed using numerical ODE solvers for the forward pass. Consider scalar-value loss function L(-), the input the result of an ODE solver:
[(x((l)) = (a0(1o)+ 1.x fo((, C [(ODSSoove(fo 41),
(10)
Then for the backward pass, compute gradients 0c/00, there are two major approaches: discretize-then--piiite and sprimizzethenn-isereeiz. (Kidger, 2022; Onken & Ruthotto, 2020). In the wwweeizzeeteen-ptiine approach, operations of ODE solvers are directly backpropagated. In contrast, the optimize-then-disereee approach, also known the adjoint sensitivity method, intro- duces additional backward-in time adjoint ODEs and calculates the gradients with an integral:
da (1, ar)
ac
do
@ x)
dt = a
a
dt.
(11)
Ox
00
00
where a 8C is the adjoint state. Other training frameworks are typically modifications of these Par
two approaches.
B.2 ERROR ACCUMULATION OF NUMERICAL ODE SOLVERS
In this subsection, we discuss the issue error accumulation in numerical ODE solvers. ODE solvers widely used in the context of NODEs be expressed
En+1 + hg(tn, 2n),
where h the step size, and g(: denotes the updating formula. In essence, numerical ODE solvers the state the current time step predict the state at the next time step, which can be viewed as a form of autoregression.
To demonstrate enror accumulation, let consider the simple ODE i The analytical solution of this ODE is xt) = (0) exp(l). We compute the mean absolute percentage error (MAPE) between the analytical solution and the numerical solution obtained using L.SODA in sc scipy py (Virtanen al., 2020). The MAPE is shown in Fig. 3. As observed in the figure, the error increases over time. For untrained NODEs, such error accumulation can lead to poor performance during training.
C PROOF OF THEOREM
According (Qian al., 2022), we present the proof of Theorem 2, which is refined for neural networks.
Proof of Theorem 2. Given the definition of c:
c(2, fo. (1) (' m(t))e(t) dt S' folt, (0)))
(12)
15

---

Published as a conference paper at ICLR 2025
MAPE between Numerical Solution and Ground Truth
0.005
0.004
0.003
MAPE
0.002
Figure 3: The MAPE between the numerical solution and ground truth (analytic solution) for x.
For the first term in Eq. (12), we can express it as:
S x(t) dt So dØt()) =
''' a(t)
(13)
=- S (t)) dt.
Substituting Eq.(13) into (12) we obtain:
c(x, for oz) = [" folt.z(t))o / a(t))r(t) dt
[" (o(()=(()-<(())))))
+(1)) o))))
(14)
f. (fo(t,z(t)) - f(t, dt.
Based on Parseval's identity, for each element c(x, for Φ2), we can obtain:
E/ (Ifo,(t,x(t)) fs x(t))) (t) dt
Σ (z.
lim L Em
fo. fo, 0x) = lim
L 1=1
1 [Jas(1,((1)) - [,(1,
(15)
= II[es - f3) (t, (1,z)11
Thus, we have proven Theorem 2.
D GENERAL NUMERICAL INTEGRATION TECHNIQUES DO NOT WORK FOR OSCILLATORY INTEGRALS
For integrals over tabular data, Newton-Cotes formula-base methods (Press, 2007), such as the trapezoidal method and Simpson's 's method, are commonly used. These methods estimate the original integral by:
Sd hett) sin(wt)
K
Σ
sin)sin(ut) dt
Pk (t) dt,
(16)
1x
where Pk (t) the approximate n-th-order polynomial of h(t) sin(w)) over the time interval (1, (1]. This approach can compute the integral accurately when the original integrand is exact
16

---

Published as a conference paper at ICLR 2025
n-th-order polynomial in [&,, fx-1] However, based on Taylor expansion, such polynomials only ensure accurate computation of an integral within enror bound of ((g(at)) for any integrand g(t) (Deaño et al., 2017). For the oscillatory integrand lt(t) sin(wt), as w increases, the (n + 1)-th derivative of the integrand also increases rapidly and can be pretty large. In this scenario, polyno- mials may struggle to accurately approximate h(t) sin(wi), causing general numerical integration techniques fail in calculating such integrals.
E SPL.NE-BASED APPROXIMATION
In this section, we delve into the application of natural cubic splines for constructing closed-form approximations from sampled data. Generally, two techniques are employed for this purpose: (1) natural cubic spline interpolation, and (2) natural cubic spline regression.
Let's begin with natural cubic splines. Given sampled data (h(t+)) R= our objective to construct cubic polynomial (t) within each interval [1x, tx+1]:
3
2 (t)
akmm (t tx = ak,u + at, - tx) + 0x,2(t tx)² + ax,3(t tx)², m==
Σ
(17)
where Ox.m (k O,..., K, m = 0, 1,2,3) are the coefficients of this polynomial. To ensure smooth approximation in each interval, the following conditions must be met:
1. Continuity: The spline must be continuous at each point: q (tx+1) = gk+1 (fx+1)-
2. Smoothness:
The first and second derivatives of the spline must be continuous each point: & (tt++)) = 4 (f*+1), and of (fx+1) = t41 (t&+1).
3. Natural boundary condition: The second derivative at the endpoints should be zero:
(((t)) - =0 0.
Now, let's examine the difference between spline interpolation and spline regression. The crucial distinction lies in their treatment of the sampled data points. Spline interpolation constructs a piece- wise cubic polynomial that exactly passes through all the given data points:
q (t&) h(tx)) 4 h(tx).
However, not adept at handling noisy data as it strictly adheres to the given points. In contrast, spline regression constructs piecewise cubic polynomial that best fits the data points in least- squares sense while maintaining smoothness. This achieved by solving the following optimization problem:
K
(1- (b(4x) - ((wx)) (1- J. (1)) di
min ok.m
(18)
k==
where (t) = (t), E [t,, 11), and À [o, 1] is hyperparameter used to control smoothness. When À = 1, this formulation reduces to spline interpolation. Spline regression is better suited for handling noisy data. The process of solving these coefficients dx.m can be found in (De Boor, 1978).
In our method, we initially employ spline regression to construct precise closed-form approxima- tions of trajectories from noisy data, thereby obtaining estimated trajectories z(t). Subsequently, we utilize spline interpolation to construct accurate closed-for approximations vector fields (t, ((t)).
F ALGORITHM OF VF-NODE
The detailed algorithm is presented in Algorithm 1. First, employ natural cubic spline regression to construct an analytical approximation of trajectories from noisy and partially observed data æ, represented as the coefficients of the spline. Next, we estimate the values in trajectories based on the spline to remove noise and fill in missing values. Then, we compute the vector fields f based on estimated trajectories and utilize natural cubic spline regression again to build analytical approximation of f. Finally, using Filon's method, we compute oscillatory integrals based on the analytical approximations of and f.
17

---

Published as a conference paper at ICLR 2025
Algorithm 1: VF Loss
Data: A trajectory ((te, ((x))) k=o- where = ly 4 < ... < lx = T, and
x(tx) = [x) (tx), x2((t) (tx)...... Ta(ta)] E (R (x
Input: Smoothing coefficient [0, 1]. the number of basis functions L, and neural network
fo(- x())
Output: The Vososs
/+Perform spline zegression on to get spline coefficients
an.m SplineReegressonn(A. {(l,, a(t))( e Rd, == K - 1, 0,1,
/MMaee estimations of the trajectory x
for k = 0 K do
if < K then
d(4x) = 04,0 + Em=i a2, (1: - 14 "1 ax,0
else
L ((tx)= = Em ak 1,m - -K-1)""
/>Evaluate the vector fields ftx) = fe(t,, 2(t+)), & 0.......f K
*/
/*Perform spline interpolation on f to get spline coefficients
*/
bymm = (0)
ER, 0,.., K-1, 1, m 0,1,2,3
/*Compute oscillatory integrals for 1:
5e a()) dt = / - Ei akm pl+t (t - tx n cos
alt dt T
5 (t) dt = EK En bam f (t t sin nēt dt
T
f(x, 0- (t)) dt l (()) return Σ/= |c(æ, f, 21)1
18

---

Published as a conference paper at ICLR 2025
G ADDITIONAL EXPERIMENTS
G.1 ADDITIONAL RESULTS: AVERAGE TRAINING TIME ON OTHER DYNAMICAL SYSTEMS
In addition to the average training time per epoch on the glycolytic model presented in Fig. 2, we also provide additional timing results for the repressilator model and the age-structured SIR model, as shown in Fig. 4 and Fig. 5, respectively. For some methods, the uncertainty may not be visible due their stability. Across different cases, the proposed VF-NODEs achieve acceleration factors ranging from 10 to 1000 times compared to the baselines.
Avg Training Time Epoch (5)
VF-NODEs
Dis-Opt
Opt-Dis
Seminorm
10
10-1
VF-NODES
TayNODE
Latent
Latent
NCDE
Figure 4: Average training time per epoch (second) for each method the repressilator model. Our method can achieve 10 to 1000 times faster than the baselines. Due to the high stability of the training speed for methods, the uncertainty is negligible and not clearly visible in the figure.
Avg Training Time Epacc (s)
VF-NODEs
10a
Dis-Opt
Opt-Dis
Seminorm
10
10
VF-NODEs
Vanilla NODE
TayNODE
Latent ODE 1
ODE-RNN
Latent ODE
NCDE
Figure 5: Average training time per epoch (second) for each method on the age-structured SIR model. Our method achieve 10 to 1000 times faster than the baselines. Due to the high stability of the training speed for some methods, the uncertainty is negligible and not clearly visible in the figure.
G.2 ADDITIONAL EXPERIMENTS ON SIMULATION FOR DYNAMICAL SYSTEMS
Besides four dynamical systems mentioned in Section 5.1, also test additional dynamical systems: (1) the Gompertz model (Gompertz, 1825) and (2) Lotka- Volterra equations (Kingsland, 1995). The Testing MSE for method and other baselines presented in Table 6. These ad- ditional experimental results also demonstrate the outstanding performance of the proposed VF- NODEs. Detailed training settings are presented in Appendix I.3.
G.3 TEMPORAL EFFECT OF CHEMOTHERAPY ON TUMOR VOLUME
We also use VF-NODES to model the temporal effect of chemotherapy on tumor volume. Using the same data preprocessing (Qian et al., 2022), the results are shown in Table 7. NCDEs failed due to the setting fiffrax (Kidger, 2022). Note that our method does not outperform all baselines due to the inherent limitations of the standard NODE architecture, which assumes that time series data must be the solution of an ODE system. In contrast, ODE-RNNs, which combine ODEs with RNNs, not bound by this constraint. However, it is important to note that the training of VF-NODEs can be almost 180 times faster than that of ODE-RNNs Detailed training settings are provided in A Appendix I.3.
19

---

Published as a conference paper at ICLR 2025
Table 6: Testing MSE (meantstandard deviation) on two additional dynamical systems with 80% observed data (r 0.8). Lower values indicate better performance. Here refers x10# Latent ODE refers to Latent ODE with an RNN encoder. Latent ODE refers to Latent ODE with an ODE-RNN encoder. The best results are highlighted in bold black, and the gecond-bess results are highlighted in bold purple.
Compertz
Lotka- Vollerra
Interpolation
Extrapolation
Interpolation
Extrapolation
Vanilla NODE
(3.18--0)=(9.266-888)
(1.76e-07)=:8.74e1100
(I2Ie=000==55:9-001)
(2.1I=+01)=(588e--
TayNODE
(5.15e-05)=(1002--03
(4.12c-06)H2..64o-00
(150e+02))46822440000
(1644c+02))((2855-0000
Latent ODE
(4.54e-00)e(3.566-00) (2991-00))(3.87--3)
(4.62++02)++92222+00)
(1669c+03)4(7116ee))))
ODE-RNN
(1.52=-05)=(14422000 (2.04e-06))(1.82eeee)
(1997=+01)=(2.50=+001
(1.32++2)=((1.444)) (.884+02)
Latent ODE
(4.376-03)=:47788-055) (2.26e-03)*{2..) (2--0))
(2.39=+03)=(1855c+03
NCDE
(3.50e-02) (5.372-04))(2.22-000) 10==0))
(3.94e+00)+((555=+000
(3. (e144+0)(((247-000
ResNet Flow
(181c-02)==(L886600) (5.63e=00))(7779<000000
(2.57=+02
2299c=+02)((3666--22
GRU Flow
(2.6600)))
(495c-02))33626-022
(3. (.60==02)==(115--00))
VF-NODE (Ours)
(2.76--0))
(1.78c-0))
(1000=+00)((311c443))
(1366=+0))44666+00))
Table 7: Testing RMSE (meanzstandard deviation) of modeling the temporal effect of chemotherapy on tumor volume. Lower values indicate better performance. Here ein refers to x x10#=. The best results are highlighted in bold black, and the second-best results highlighted in bold purple.
RMSE
Vanilla NODE
(1.899-03))+99.66-06))
TayNODE
(2.31e-03)+(1.430-05))
Latent ODE (RNN Enc.)
(1.26-0)))
ODE-RNN
(1.590-03HH(443e-044)
Latent ODE (ODE-RNN Enc.)
(9.022-03)+((4442-033)
NCDE
(nan)+(naa)
ResNet Flow
(1.88--03)+((4446-00)
GRU Flow
(3.146-03))(3.666-0))
VF-NODE (Ours)
(1.87--33)+(3.71e-55)
20

---

Published as a conference paper at ICLR 2025
G.4 ADDITIONAL EXPERIMENTS ON ABLATION STUDIES
In this subsection, we first evaluate the sensitivity of the proposed VI-NODES to two hyper- parameters: the number of basis functions L and spline regression hyper-parameter À mentioned in Appendix E. Subsequently, we assess the performance of VF-NODEs under varying sampling and noise conditions.
The number of basis functions L. We set the number of basis function L to test the impact of L on the performance and training speed of VF-NODEs by setting L 50, 60, 70, 80, 90, 100, 110, as shown in Fig. 6 and Fig. 7, respectively. We can see that when L increases, the testing MSE decreases generally. This can be attributed to Filon's 's method and spline regression, which enable us compute oscillatory integrals precisely from noisy and partially observed data. Additionally, Fig. demonstrate that training speed of VF-NODES is not sensitive to L.. Considering the computational cost, set L = 80 in all experiments.
Glycolytic
Toggle
0.0003
Interp
Extrap
0.0002
MSC
0.0002
0.0001
Interp
0.0001
Extrap
0.0000
80
90
100
110
50
60
70
160
110
# Basis Functions
Basis Functions
Figure 6: The impact of the number of basis functions in the VF loss on model performance
Glycolytic
Toggle
Average Time (5)
0.025
0.025
0.020
0.020
0.015
50
60
70
80
90
100
110
50
60
70
80
90
100
110
Basis Functions
Basis Functions
Figure 7: The impact of the number of basis functions T. in the VF loss training speed of VF- NODEs.
As discussed in Appendix controls the smoothness of the splines. Typically. À values greater than 0.9 are suitable for most cases. To evaluate the sensitiv- ity of the proposed VF-NODES to we test VF-NODEs with 0.9, 0.99, 0.9999, and 0.99999. The experimental results are presented in Table 8. These results indicate that our method is robust to variations in X. In practice, À can be selected based on the noise level of the sampled trajecto- ries, with higher noise levels generally requiring smaller values of X, making it straightforward parameter to adjust.
Spline Regression Elyper-paramete
Table 8: Testing MSE (mcantstandard deviation) from the ablation study on the sensitivity of VF- NODEs to X. Lower values indicate better performance. Here ein refers to x
=10
Glycolytic
Toggle
A
Interpolation
Extrapolation
Interpolation
Extrapolation
0.9
(4.20e-03)+(344ee00)
(4.81e-03)4(6003e-0))
(5 (2) 04)
(7.5e-006
0.99
(1.14e-03)(2.000044) (2. (..16-00)) (1.76-003)((.266-04))
(8.296-00))
0.999
(7.276-09)+(65990004) (1.79e-03) (7.86-04))66.76-004)
(1.526-05)499.322-00
0.9999
(5.322-09))(44476004) (8.94e-03)x(11096-000 (1.45-00
(3.40c-02)#(4.80eee
0.99999
(6.35--05)1(2.6e-066)
(1.65-04))
(1669--04)(66006-055))
(4.79--0))+55.244-08
Effect of sampling. There are three key settings related to sampling: (1) the number of data points per trajectory, (2) the random dropping rate, 1 -r, and (3) sampling interval T. First, we fix 0.8,
21

---

Published as a conference paper at ICLR 2025
T = 10 and sample 50, 100, 300, 500, 700, and 1000 points for each trajectory, respectively. The experimental results are presented Table 9. Next, we fix the number of points at 100 and test with 0.4, 0.6, and 0.8, presented in Table 10. Finally, vary the sampling interval, setting T = 10, 20, and 50, while considering 100, 200, 500 sampled points respectively, as reported in Table 11. These experiments demonstrate that our method is robust across different sampling settings, particularly on sparsely sampled data. This can be attributed to spline regression's 's ability perform effectively on sparse datasets (Knowles Renka, 2014).
Table Testing MSE (meantstandard deviation) from the ablation study on the number of sampling points for each trajectory. Lower values indicate better performance. Here ein refers to x 10
Glycolytic
Toggle
# Point
Interpolation
Extrapolation
Interpolation
Extrapolation
50
(7.08e-04)=(3366--000
(1.56e-03))(8.333e00
(4.56e-00)+(84456-0))
(6.35e-05)4(2.686-06) (1.63--0)))43305--055))
(1.69--04)466009--055
(4. (4.79-007)66..24-000)
300
(4.940)))))(1.7400055 (1.200-04))((3722-005))
(2009c-44)=(722440000
(1.736-07)=((46655-00))
500
(4.39-0))) (1001=04))(3111e0055)
(2.022-04)=(1.0e=004))
(2.21e-07)=(70076-08)
700
(4.46e-05)4(1339-000 (8.39-0)5)
(1 (1.25-04))(5.155-055))
(1.05=007)<(5.15-008)
Table 10: Testing MSE (meandstandard deviation) from the ablation study for the random dropping Lower values indicate better performance. Here ečn refers to x 10
rate 1
Glycolytic
Toggle
r
Interpolation
Extrapolation
Interpolation
Extrapolation
0.4
(1006-0))))((2226--055
(1466c-04)+(14422--)))
(5.84c-04)t(1 19cc000)) (1.5lc-06))(4224000))
0.6
(1772-004)
(6.75--04) (7.236-77)+(23160007)
0.8
(6.35600)))((.668-00)) (16.3--04))(3.05-4455)
(1660-504))((600660555)
(4.79:-07)+(55544088))
Table 11: Testing MSE (meanzstandard deviation) from the ablation study on the sampling interval for each trajectory. Lower values indicate better performance. Here ein refers to x10
Glycolytic
Toggle
T
Interpolation
Extrapolation
Interpolation
Extrapolation
10
(L.63e-0)) 68635
(1.69--04)+(600060055)) (44799-07)
20
(4. (400))00)))4
(8.22-))) (5.555)))) (2.288-07)+(4255--08)) (3.2200))) (2334--33)+00119-00)) (388--04) (1.90--07)+(8.222-088)
50
Effect of noise. We set the relative noise level az 0.01, 0.03, 0.05, 0.07, 0.1, 0.3, and 0.5 respectively. The experimental results are shown in Fig. 8. By setting the smoothing hyperparameter 0.99 in spline regression (as discussed in Appendix E), VF-NODEs demonstrate robustness to different noise settings.
Glycolytic
Toggle
10-1
Interp
10-1
Interp
10
Extrap
Extrap
MSE
MSE
10-3
10-3.
10 5
10
10
10" Relative Noise Level 0%
Relative Noise Level de
Figure 8: Testing MSE on the effect of relative noise level ar
22

---

Published as a conference paper at ICLR 2025
H DISCUSSION: TRAINING NODES BASED ON NATURAL CUBIC SPLINES
The proposed VE-NODES utilize natural cubic splines numerically compute integrals in the VF loss. Beyond this, some other training methods for NODEs also leverage splines. One such method gradient matching (Varah, 1982: Li al., 2024). where the derivatives of sampled trajectories are estimated numerically using splines or other spectral methods. The DNN-based vector fields are then trained to match these estimated derivatives. In the second method, given the sampled trajectories (=< the DNN-based vector fields can be approximated as
ft(t, ((t)) x 0 +
by.m (t. 1, t
[1x, [x+1].
(19)
m1=1
Using this approximation, the solution of NODEs at the sampled time points can be computed as
le+1
((t++1) = x(xx)= +
folT, (()) dr
flatl
23 (((x))))
Jtn
bs.0 +
ba mx
dr
(20)
m=m1
(1++1
(x) m+
= =(tx)
bk.m
m+1
m7
The NODE can then be trained by minimizing the MSE between (tx) and ((t)). All these methods are ODE-solver-free free, similar to the proposed VF-NODES.
We refer to these methods Grad-Matching NODEs and Spline-Integ NODEs, respectively. We compared the proposed VF-NODEs with these methods, with experimental results presented in Ta- bles 12 and 13. While the training speed of these methods is similar (approximately 10-2 second per epoch), we observe that their performance significantly worse than that of VF-NODEs.
For Grad-Matching NODEs, this performance gap be attributed to their sensitivity to noisy and sparse sampling. In contrast, VF-NODES rely solely on numerical integrals, representing an improvement and generalization of gradient matching methods Brunel et al. (2014).
For Spline-Integ NODEs, the discrepancy arises due to the autoregressive nature as described in Eq. (20). Errors accumulate across iterations because each step depends the approximation from the previous one. VF-NODEs, by relying on global numerical integrals, avoid this issue of error accumulation, resulting more robust and accurate training.
Testing MSE (meantsstaddard deviation) for interpolation task of NODEs trained with splines on dynamical systems with 80% observed data (r = 0.8). Lower values indicate better performance, Here ein refers to x10#". The best results are highlighted in bold black.
Table 12:
Glycolytic
Toggle
Repressilator
AgeSIR
Grad-Matching NODE
(I. (8.19:00)) (7.55-00)))
(9.80-00)++(2200-000 (9.50e-03)4(2.00))
Spline-Integ NODE
(1.34403))(1.4440333)
(9.56-0032(7.766003)) (1.000022)110000022 (3.48-0211((.44-003
VF NODE (Ours)
(6.355-05))((2666-006) (1.00-00))=(60000055)) (1.922-0))+(2222-04)) (7.39--03)466.71ee04)
Table 13: Testing MSE (meantstandard deviation) for extrapolation task of NODEs trained with splines on dynamical systems with 80% observed data (r 0.8). Lower values indicate better performance. Here etn refers to x x10#". The best results are highlighted in bold black.
Glycolytic
Toggle
Repressilator
AgeSIR
Grud-Matching NODE
(6.322-04))(11866-00) (4:436-06)((222-066))
(3.77-00))a(4.90000)
(..16600))
Spline-Integ NODE
(5.44--0)) (5.226-05)+3.555-05))
66005--01)2(7550-000
(1
VF NODE (Ours)
(L.65e-0)))(.05--005)
(4.746-07)+(52446088)
(1.23--01=+(.488002)
(2.376-02):(L661ee033)
23

---

Published as a conference paper at ICLR 2025
I EXPERIMENTAL SETTINGS
All experiments in this work are implemented using jax (Bradbury et al., 2018). Specifically, the implementation of neural differential equation models is based on equinox (Kidger Garcia, 2021) and di ffrax (Kidger, 2022). To optimize models, we use the optax (DeepMind et al., 2020). All the experiments are implemented on the same server, equipped with 4 A5000 GPUs with 24GB graphics memory.
1.1 HYPERPARAMETERS FOR MODELS
To evaluate the performance of VF-NODEs, compare them with two categories of baseline mod- els: NODE-based models and Neural Flows. To maintain consistent parameters across all models, employed the following hyperparameters:
VF-NODEs: the vector field is parameterized as 4-layer MLP with 128 hidden units per layer. For the VF loss:
The number of basis functions L is set to 80 for all tasks.
the smoothing hyperparameter À is set 0.99999 for most tasks, while for the simu- lation of the age-structured SIR model, is set to 0.9999.
Vanilla NODES (Chen et al., 2018): the vector field is parameterized as 4-layer MLP with 128 hidden units per layer.
TayNODEs (Kelly et al., 2020):
The neural architecture of TayNODEs the same as that of Vanilla NODEs.
- The 5-th order derivative regularization is used to match the order Dopri5 solver. is set to 0.001.
Latent ODEs with RNN encoders (Chen et al., 2018):
RNN encoder: I-layer GRU with 25 hidden units.
NODE Decoder: the vector field is parameterized 4-layer MLP with 124 hidden units per layer.
-The latent size is set to 4.
ODE-RNNs Rubanova al. (2019): 1-layer GRU with hidden units 4-layer MLP with 127 hidden units per layer for NODE.
Latent ODEs with ODE-RNN encoders (Rubanova et al., 2019):
ODE-RNN Encoder: 1-layer GRU with 4 hidden units and 4-layer MLP with 84 hidden units per layer for NODE.
NODE Decoder: A 4-layer MLP with hidden units per layer for NODE.
Latent size set to 4.
NCDEs (Kidger et al., 2020): The hidden size is set to 4. A 4-layer MLP with 89 hidden units per layer is used for the vector field.
ResNet Flow (Biloš al., 2021): 4-layer MLP with 128 hidden units per layer.
GRU Flow (Biloš al., 2021): Three 4-layer MLP with 74 hidden units per layer.
In most experiments, we used the Dopri5 solver for NODE-based models. However, because training ODE-RNNs is computationally intensive, we employed the Midpoint solver for these models to reduce computational complexity.
I.2 SETTINGS FOR DYNAMICAL SYSTEMS
In this subsection, we provide detailed information on the dynamical systems used in this work, including their trajectories and the simulation results of VF-NODEs compared to the best baseline, as reported Tables and 2. For this study, we selected six dynamical systems. The specific parameters and trajectories for each system detailed below.
24

---

Published as a conference paper at ICLR 2025
Glycolytic oscillator Kov, 1968). The glycolytie oscillator is fundamental system in bio- chemistry that models the glycolysis process. can be expressed as
i 01 0221 2 21022
(21)
お = -2
32 + 0g.01 + 21
where 01 0.75, 82 = 83 = 0.1, and æ1 (0), x2 [0.1, 1.1].
3.0
Ground Truth
Ground Truth
VF-NODE
2.0
VF-NODE
Best Baseline
Best Baseline
2.0
1.5
1.5
1.0
1.0
0.5
0.5
10
15
15
20
Figure 9: The trajectory plots of the glycolytic model with simulation results of VF-NODEs and the best baseline (ODE-RNN).
Genetic toggle switch (Gardner et al., 2000). The genetic toggle switch key mechanism in genetic engineering and synthetic biology for controlling genes. It can be expressed as
01
21:
+ 22
(22)
02
42
02 T2:
1 + x
where d1 = 02 = 4, A1 = 122 = 3, and x (0), (0) E 4.0).
5.0
Ground Truth
Ground Truth
VF-NODE
0.4
VF-NODE
4.5
Best Baseline
Best Baseline
4.0
0.3
0.2
3.0
2.5
0.1
10
15
20
5
10
15
20
Figure 10: The trajectory plots of the toggle switch model with simulation results of VF-NODEs and the best baseline (Vanilla NODE).
Repressilator (Elowitz & Leibler, 2000). The repressilator genetic regulatory network. It can be expressed as
a P3 5: = -3(p: m,), where P; are three repressor- pppressor-protein concentrations, and m; are corresponding mRNA concentra- tions. We
ri Th +
+ ao: lacI, tetR, cl,
(23)
j = cl, lacI, tetR,
10, og 10-5, -,, 3 3, and (0))),((0)) E [0,
25

---

Published as a conference paper at ICLR 2025
15
Ground Truth VF-NODE
Ground Truth VF-NODE
Ground Truth
12.5
VF-NODE
Baseline
10.0
Baseline
10
Baseline
7.5
10
5
100
15
0
5
100
12.5
Ground Truth
Ground Truth VF-NODE
Ground Truth
VF-NODE
10.0
10.0
Best Baseline
Best Baseline
10.0
Best Baseline
&
Fore
5.0
5.0
5.0
2.5
5
10
10
10
Figure 11: The trajectory plots of the repressilator model with simulation results of VF-NODEs and the best baseline (Vanilla NODE).
Age-structured SIR model (Ram & Schaposnik, 2021). The age-structured SIR model a variant of the standard SIR model that considers the effects of different age groups. It can be expressed
5; = -B S N
MiyIj,
j=1
S
(24)
1, = 8
Mi
7141
N
==1
R2 = 7/4
where S, and R; (i . 9) denote the numbers of susceptible, infected, and removed individuals, respectively, for the age groups 0-9, 10-19, 70-79, and 80+. The age-contact matrix M is parameterized
19.2
4.8
3.0
7.1
3.7
3.1
2.3
1.4
1.47
4.8
42.4
6.4
5.4
7.5
5.0
1.8
1.7
1.7
3.0
6.4
20.7
9.2
7.1
6.3
2.0
0.9
0.9
7.1
5.4
9.2
16.9
10.1
6.8
3.4
1.5
1.5
M
7
7.5
7.1
10.1
13.1
7.4
2.6
2.1
2.1
3.1
6.3
6.8
7.4
10.4
3.5
1.8
1.8
2.3
1.8
2.0
3.4
2.6
3.5
7.5
3.2
3.2
1.4
1.7
0.9
1.5
2.1
1.8
3.2
7.2
7.2
1.4
1.7
0.9
1.5
2.1
1.8
3.2
7.2
7.2
In addition, we set 3 0.8, 7 = 0.5, and (,(0),,((0), (O) E [0.1, 10.1]
(i
1,...,9).
Gompertz model (Gompertz, 1825). The Gompertz model widely applied in medical research and tumor growth analysis as kind of growth model. It can be expressed as
t 0 log(02x),
(25)
where 0, 02 1.5, and (0) [0.1,
Lotka- Voltera equations (Kingsland, 1995). The Voltera equations are used model the interactions between the predator and prey populations time, capturing how the population sizes of each species affect the other. This system can be expressed as
0.T
1
(26)
i
, - 7/
26

---

Published as a conference paper at ICLR 2025
Ground Truth
Ground Truth VF NODE Best Baseline
VF-NODE
Ground Truth
6
Best Baseline
VF IF-D
Best Baseline
10
20
15
20
5
10 t
20
5.0
Ground Truth
Ground Truth
20
5
2.5
VF-NODE
Best Baseline
5
VF-NODE
Best Baseline
&
15
VF-NODE
Ground Truth
Best Baseline
0.0
10
10
15
20
10
20
1.0
Ground Truth
Ground Truth
12.5
VF NODE
VF NODE
0.5
Best Baseline
-22
Best Baseline
VF-NODE
0.0
5
10
15
20
10
15
20
5
10
15
20
t
Ground Truth
Ground Truth
15-
07
VF-NODE
Baseline
5
VF-NODE
Baseline
10-
VF-NODE
Ground Truth
Best Baseline
10
15
t
10
Ground Truth
Ground Truth VF-NODE
5
VF-NODE
20
Ground Truth
67
Best Baseline
Best Baseline
VF-NODE
Best Baseline
0
01
10-
10
15
20
10
15
20
0
10
15
20
41
Ground Truth
Ground Truth
VF NODE
VF -F-NODE
Ground
68
Best Baseline
2.5
Best Baseline
VF-NODE
Best Baseline
0.0
5
10
20
15
20
5
10 t
15
20
10
5.0
Ground
Ground
15
VF-NODE
VF-NODE
Ground Truth
5
2.5
Baseline
=
Best Baseline
& 10
VF-NODE
5
Best Baseline
Ground Truth
Ground Truth
5
VF NODE
VF-NODI
Ground Truth
Best Baseline
Best Baseline
VF-NODE
10
Best Baseline
5
10
15
20
0
10
15
20
0
5
10
15
20
Ground Truth
Ground Truth
15
68
VF-NODE
Best Baseline
VF-NODE
Best Baseline
10
VF NODE
Ground Truth
5-
Best Baseline
0
10 t
20
15
20
5
10 t
20
Figure 12: The trajectory plots of the age-structured SIR model with simulation results of VF- NODEs and the best baseline (Vanilla NODE).
27

---

Published as a conference paper at ICLR 2025
Ground Truth
0.8
Best Baseline
0.7
0.6
0.5
0.4
0.3
0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
Figure 13: The trajectory plots of the gompertz model with simulation results of VF-NODEs and the best baseline (Vanilla NODE).
Ground Truth VF-NODE Baseline
14
Ground Truth
VF-NODE Baseline
80
50
30
5.0
10.0
15.0
20.0
0.0
5.0
10.0
15.0
20.0
Figure 14: The trajectory plots of the Lotka- Volterra equations with simulation results of VF- NODEs and the best baseline (Vanilla NODE).
28

---

Published as a conference paper at ICLR 2025
where 1.0, 3 0.3, 8 = 0.1, 3.0, and (0) [10.0, 20.0), z2( E [5.0, 10.0].
To generate datasets from these dynamical systems, the Dopri5 solver in di f ffrax is used.
1.3 TRAINING SETTINGS
For all tasks, employed the Adam optimizer (Kingma & Ba, 2014), loading all training data in single epoch. The testing MSE losses of NODE-based baselines reported in all tables were evaluated based on models trained using the optimize- then-discretize approach, following (Kidger, 2022).
For all tasks except modeling the temporal effect of chemotherapy on tumor volume, we set the number of training epochs 5,000; for the tumor volume modeling task, we used 300 epochs.
For all tasks except the simulation of the age-structured SIR model, we utilized the wew.innonnccceeeeeeee
from optax as the learning rate scheduler, with an initial leam- ing rate of 0.001. The scheduler parameters were set as follows: transition steps equal the number of epochs, peak va I uu at 0.01, pet start at 0.2, div factor 100, and at 1,000. For the simulation of the age-structured SIR model, we employed
final Shaalldiv_fcttte the cosine wen.eddeca.sccheeeee
e as the learning rate scheduler, also with an initial learning rate of equal to the number of epochs, alpha
0.001. The parameters were set as follows: at 0.01, and exponent at 1.0.
J DISCUSSION: APPLYING VF-NODES TO CHAOTIC SYSTEMS
In this section, we discuss the capability of VF-NODEs to model chaotic systems. Since making long-term precise predictions for chaotic systems is nearly impossible, employ the unstable periodic orbit (UPO) detection using adaptive delayed feedback (ADF) technique (Zhu et al., 2023) evaluate whether VF-NODES capture the underlying patterns of chaotic systems. Consider the Lorenz system, defined as:
t oly x))
ý x(p - =) - y: xy A=z,
(27)
where we set 10, 8/3, and p 20, and initial conditions x00), =(), E 101 10). Then, following (Zhu et al., 2023), we apply control term c(t) (t))((t) to where e(t) = yt ait-pla p(t)) yt).. The variables r, y. 2, p. and y are determined following controlled system:
å o x),
ý = x(p z) + c(t),
1 xy Az,
(28)
P = ri (v(() - r(t p(t)),
4 = relu(t) - y(t p(t))
where we set p(0) 70) = 0.1, = 0.2 and = 20. In this setup, p(t) reflects the underlying pattern the chaotic system. The mean value p(t) on an interval with the minimum variance can be viewed a UPO for Eq. (27).
These results demonstrate that VF-NODES can effectively capture underlying patterns of chaotic systems, significantly outperforming Vanilla NODEs in these challenging scenarios. This improve- ment is largely due to VF-NODEs mitigating autoregression during the training process, thereby avoiding error accumulation. In contrast, Vanilla NODEs rely heavily ODE solvers, where error accumulation during training can significantly impact the results.
29

---

Published as a conference paper at ICLR 2025
0.6
0.4
a
Ground Truth VF-NODE Vanilla NODE
0.2
0.0
2.5
5.0
7.5
10.0
Figure 15: The comparison of ptt) between the VF-NODE and ground truth system in Eq. (28). The detected UPO for the ground truth is 0.6108, for VF-NODE (0.5963 + 0.0002), and for Vanilla NODE (0.5010 + 0.0150).
20
Ground
50
20
VI-NDDE
10
NCDE
30
-10
10
10
20
5
10
15
20
10
15
20
Figure 16: The trajectory plots of the Lorenz system with simulation results of VF-NODEs and Vanilla NODEs. The testing MSE error is (49.7407 + 18.5300) for VF-NODEs and for (54.4111 2.5761) Vanilla NODEs.
30


Reviews:
[{"review_id": "PiNSoqeQPV", "reviewer": "Reviewer_1H9x", "summary": "The authors present a new method for training Neural ODEs (NODEs) based on a variational formulation of an ODE loss. The resulting method used spline regression and to interpolate noisy data which allows for a computation of a variational loss. The advantage of the method is many few function evaluations during training.", "questions": "In the COVID-19 dataset example is the data gathered from real world observations? Or is the data generated from some parametric model of COVID-19 spread?\n\nDo the authors have a sense of the limits of the spline regression / interpolation? On what sorts of trajectories it might fail? \n\nIt would be very helpful if the authors could provide plots of training trajectories the ODEs in consideration so the reader could assess the noise levels / complexity of the trajectories.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The strengths of the paper are:\n- the proposed method introduces a reasonable approach to training a NODE based on a variational formulation of the loss\n- The training method requires many fewer function evaluations, allowing for much faster training than existing methods.\n- The method maintains or at times outperms competing methods in terms of accuracy.", "weaknesses": "The weaknesses of the paper are:\n- Many of the examples seem somewhat toy. It is unclear how the method might extend to much noisier / less structured ODEs with more complex dynamics. In particular one would expect the  spline regression / interpolation to eventually fail on very long, non-smooth or noisy trajectories.\nFrom the paper it is difficult to get a sense of how much these attributes are present in the given benchmarks. \n- Implementation of the method is not straightforward from a practitioner's point of view.\n- The method introduces a number of hyper parameters required for accurate interpolation which must be chosen/tuned.\n- The benefit of reduced training, while convenient, is not broadly important for many of the problems considered.", "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": "", "extra_scores": {}}, {"review_id": "OrDWcrTJX0", "reviewer": "Reviewer_Ycyv", "summary": "This paper proposes a new way to train NeuralODEs that is faster than current approaches. Current training strategies require numerically solving the differential equation which is very computationally expensive due to the high number of function evaluations at each time step. The approach introduced here relies on the variational formulation as a surrogate objective. However, the variational formulation still requires the value of the vector field and underlying time series at each time step, so the authors approximate this with cubic splines. To compute oscillatory integrals, the authors incorporate the Filon method. Finally, the authors demonstrate the performance of their method and show impressive computational time gains. Some prediction performance gains are also observed, due to the non-accumulation of errors in their training approach.", "questions": "As stated above, I would like the authors to discuss the limitations of using an explicit interpolation of the dynamics before training.\n\n4.2 Step 4, the fact of interpolating the vector field using the values only at the observation points seems strange to me as the interpolation will not coincide between observations. That is, between observations, with $\\hat{f}$ the interpolation of the vector field:\n$ \\hat{f}(t) \\neq f(\\hat{x}(t))$. Can the authors give more details about this discrepancy and motivate why it makes sense ? \n\nDuring training, your method requires computing cubic splines coefficients for each time series every time (at least for the $b$ coefficients). Can the author elaborate on the compuational cost of such a procedure ? I would also like the authors to explain how the gradient with respect to $\\theta$ can still flow from the computation of the $b$ coefficients - how is this end-to-end differentiable ?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Speeding up the training of neural ODE methods is a very important challenge with potentially high impact. \n\nThe method proposed shows clear computational time improvements as well as some prediction performance gains. \n\nThe method makes sense, and the paper is well written and easy to read.", "weaknesses": "One significant weakness that I see with this approach is that it requires defining the trajectories of the time series a priori (here with cubic splines). As such, the model with learn a vector field that agrees with the cubic spline interpolation and will not figure out alternative dynamics. If this is correct, this is a signicant limitation as one important application of NeuralODEs is for data imputation. I would like the authors to discuss that limitation more clearly in the paper and/or to argue against the reasoning above. \n\nGiven that you use cubic splines to interpolate the dynamics (and vector field), another baseline can now be considered. That is, you can now use the same cubic spline acceleration to directly integrate Equation 1 and compute the MSE at the observation points. Is such an approach reasonable ? This should be considered as an additional baseline, that doesn't use the variational approach but the cubic spline interpolation.", "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": "", "extra_scores": {}}, {"review_id": "B2KIjpmq9I", "reviewer": "Reviewer_BidZ", "summary": "This paper introduces VF-NODE, a new method for accelerating the training of Neural ODEs by using a variational formulation that evaluates vector fields only at observed data points, reducing function evaluations (NFEs) and eliminating autoregression, thereby minimizing error accumulation. The approach integrates Filon’s method to handle oscillations in the loss function and uses natural cubic spline regression to manage noisy or incomplete data. Experiments show VF-NODE is 10x to 1000x faster than traditional methods while maintaining or improving accuracy.", "questions": "1. The author utilized natural spline interpolation to fit the orbit **$x$** as well as the vector field $f_\\theta$ in this process. Then, why not directly adopt the strategy of gradient flow matching for training (see literature [1]). This approach seems more straightforward and bypasses the high computational costs associated with traditional numerical integration.\n\n2. As the author did not present the trajectory plots in experiments, I am concerned about the performance boundaries of this method in handling time series data. Is this method only effective on data with simple behaviors, or can it still outperform traditional NODE methods for more complex systems (such as chaotic systems) ?\n\n3. When using natural spline interpolation, is there an overfitting issue, such as the Runge phenomenon?\n\n4. This method directly inputs the data from the original system, thus it cannot train through modeling latent variables (like latent ODE method), which may result in a loss of the method's flexibility.\n\n[1] Li X, Zhang J, Zhu Q, et al. From Fourier to Neural ODEs: Flow matching for modeling complex systems[J]. arXiv preprint arXiv:2405.11542, 2024.\n\nIf the authors can answer the above questions well, I would be happy to consider raising the score.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This paper introduces a novel approach using variational formulation (VF) to greatly accelerate Neural ODE training, reducing function evaluations and improving accuracy. By integrating Filon's method with cubic spline regression for handling oscillatory integrals, the method achieves 10x to 1000x faster training across various dynamical systems with higher or competitive accuracy, demonstrating originality, technical rigor, and real-world applicability.", "weaknesses": "The author accelerates the training of the original NODEs method by adopting the VF-NODEs approach. However, the necessity and rationality of transforming the training problem into the optimization problem of Equation 6 still require further elaboration.  See Questions for details.", "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": "", "extra_scores": {}}, {"review_id": "gBCwxyBbDj", "reviewer": "Reviewer_wbux", "summary": "**Edit:** I have increased my score and confidence following the author responses.\n\n\nThis paper proposes a new method for speeding up the training of Neural ODEs. The method removes the need for an ODE solve during the forward and backward pass of a Neural ODE.\n\nThe theory behind the approach is to rewrite integration by parts:\n\n$\\int_0^T x\\dot{\\phi} dt = [x\\phi]_{0}^{T} - \\int_0^T \\dot{x}\\phi dt$ \n\nIf $\\phi(0)=0$ and $\\phi(T)=0$, then $\\int_0^T x\\dot{\\phi} dt + \\int_0^T \\dot{x}\\phi dt = 0$.\n\nThe aim is to approximate $\\dot{x}$ using a neural network $f_\\theta(x, t)$. And so rather than solving the ODE for $f_\\theta$ and applying MSE on the observations, the method attempts to minimise $\\int_0^T x\\dot{\\phi_l} dt + \\int_0^T f_\\theta \\phi_l dt$, for a set of orthonormal basis function $\\phi_l$ that are zero at the boundaries of the solve.\n\nA fourier basis is used as a natural choice for $\\phi_l$, and cubic splines are used to approximate $x(t)$ in the integrals to give them analytical solutions.\n\nExperiments demonstrate good performance measured by MSE, and a speed up in training time.", "questions": "- Is it possible to include two papers in the related work: 1) STEER: Simple Temporal Regularization For Neural\nODEs and 2) Interpolation Technique to Speed Up Gradients Propagation in Neural ODEs?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The writing is strong\n- The solution is novel and interesting, it's particularly interesting to see an approach which does not solve the ODE in the forward pass\n- The theory behind the method is strong\n- The experiments demonstrate good MSE\n- The ablations on the basis functions demonstrate using a Fourier basis is a good choice", "weaknesses": "- I found the evaluation confusing. As far as I can tell the only place where any wall clock timing is carried out is in Figure 2, on one dataset, the glycolytic model. All other tables seem to be about MSE. However, the main claims and messages of the paper are about accelerating neural ODEs. If the main claim is that this approach speeds up training then this should be shown more across all datasets. Are these numbers available, have I missed them?\n- Figure 2 shows no uncertainties on the time taken to train.\n- There is an error in line 137-138. If there is regularisation applied during training, the training time can also be reduced, since later iterations are made faster after regularisation is applied in earlier iterations of training.\n- Following on from this, the experiments would be made more convincing by testing some of these regularisation methods. For example minimising higher order derivatives. Currently the main baselines for speed are discretise-then-optimise, optimise-then-discretise and Seminorms, of which only Seminorms claim to speed up Neural ODEs. This point would not be as relevant if the main claims of the paper are about MSE rather than time to train.\n- It is not clear what ODE solver is used in the experiments, if adaptive solvers are used it is harder to make the argument that the forward solve is slow. The proposed method is linear in the number of observed points (to do cubic interpolation), however if there are many observed points this method will be slower compared to an adaptive solver if the trajectories are quite smooth. Have adaptive solvers been tested? The reverse is also true, if the dynamics are complex an adaptive solver can take smaller steps whereas the proposed method might suffer from accuracy issues as identified in the paper.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)