Homotopy-based training of NeuralODEs for accurate dynamics discovery
Authors: Joon-Hyuk Ko
Hankyul Koh
Nojun Park
Wonho Jhe
Venue: Submitted to ICLR 2023 (2023)
Keywords: neural ordinary differential equation
synchronization
homotopy
dynamical systems
physics

Abstract: Conceptually, Neural Ordinary Differential Equations (NeuralODEs) pose an attractive way to extract dynamical laws from time series data, as they are natural extensions of the traditional differential equation-based modeling paradigm of the physical sciences. In practice, NeuralODEs display long training times and suboptimal results, especially for longer duration data where they may fail to fit the data altogether. While methods have been proposed to stabilize NeuralODE training, many of these involve placing a strong constraint on the functional form the trained NeuralODE can take that the actual underlying governing equation does not guarantee satisfaction. In this work, we present a novel NeuralODE training algorithm that leverages tools from the chaos and mathematical optimization communities -- synchronization and homotopy optimization -- for a breakthrough in tackling the NeuralODE training obstacle. We demonstrate architectural changes are unnecessary for effective NeuralODE training. Compared to the conventional training methods, our algorithm achieves drastically lower loss values without any changes to the model architectures. Experiments on both simulated and real systems with complex temporal behaviors demonstrate NeuralODEs trained with our algorithm are able to accurately capture true long term behaviors and correctly extrapolate into the future.


Under review as a conference paper at ICLR 2023
HOMOTOPY-BASED TRAINING OF NEURALODES FOR ACCURATE DYNAMICS DISCOVERY
Anonymous authors Paper under double-blind review
ABSTRACT
Conceptually, Neural Ordinary Differential Equations (NeuralQDEs) pose an at- tractive way extract dynamical laws from time series data, as they are natural extensions of the traditional differential equation-based modeling paradigm of the physical sciences. In practice, NeuralODEs display long training times and sub- optimal results, especially for longer duration data where they may fail to fit the data altogether. While methods have been proposed to stabilize NeuralODE train- ing, many of these involve placing strong constraint on the functional form the trained NeuralODE can take that the actual underlying governing equation does not guarantee satisfaction. In this work, we present novel NeuralODF: training algorithm that leverages tools from the chaos and mathematical optimization com- munities synchronization and homotopy optimization for breakthrough in tackling the NeuralODE training obstacle. We demonstrate architectural changes unnecessary for effective NeuralODE training. Compared the conventional training methods, our algorithm achieves drastically lower loss values without any changes the model architectures. Experiments both simulated and real sys- tems with complex temporal behaviors demonstrate NeuralODEs trained with our algorithm are able accurately capture true long term behaviors and correctly extrapolate into the future.
INTRODUCTION
Predicting the evolution of a time varying system and discovering mathematical models that gov- em it is paramount to both deeper scientific understanding and potential engineering applications. The centuries-old paradigm to tackle this problem was to either ingeniously deduce empirical rules from experimental data, mathematically derive physical laws from first principles. However, the complexities of the systems interest have grown so much that these traditional approaches are now often insufficient. This has led to growing interest in using machine learning methods to infer dynamical laws from data.
One school of thought, such the seminal work of [Schmidt & Lipson 2009) Brunton et al. (2016), focuses on deducing the exact symbolic form of the governing equations from data using techniques such as genetic algorithm or sparse regression. While these methods are powerful in that they output mathematical equations that are directly human-interpretable, they require prior infor- mation on the possible terms that may enter the underlying equation. This hinders the application of symbolic approaches to scenarios where there is insufficient prior information on the possible candidate terms, or complex, nonlinear systems whose governing equations involve non-elementary functions.
On the other hand, neural network-bant -base methods, such as Raissi et al. (0018) leverage the univer- sal approximation capabilities of neural networks to model the underlying dynamics of the system without explicitly involving mathematical formulae, Of the various architectual designs in literature, Neural Ordinary Differential cause these seamlessly incorporate neural networks inside ordinary differential equations (ODES), thus bridging the expressibility and flexibility of neural networks with the de facto mathematical lan- guage the physical sciences. Subsequent works have expanded on this idea. including blending NeuralODEs with partial information on the form of the governing equation to produce "grey-box"
Chen et al. (2018) stand out in particular be-

---

Under review as a conference paper at ICLR 2023
dynamics model (Rackauckass et al. (021), and endowing Neur with mathematical structures that the system must satisty (Greydanus et al. 2019 Finzi et al. 2020p.
However, despite the conceptual elegance of NeuralODEs, training these models tend to result in long training times and sub-optimal results, a problem that is further exacerbated as the length of the training data grows (Ghosh al., 2020 Finlay
al., 2020). Different methods have been proposed to tackle the problem, but majority of these approaches to date involve placing either strong
(Choromanski al. 2020 Hasani al., 2021), ) or semi-strong constraints(Finlay et al. 2020 Kidger et al., 2021) to the functional form the NeuralODE can take something the underlying governing
equation does not guarantee satisfying.
Contributions. We introduce a novel training algorithm that does not require architectural con- straints to accurately train NeuralODEs on long time series data. As our algorithm inspired by ideas from the chaos and mathematical optimization literature, we provide a background survey on the ideas involved before providing both a general framework and specific implementation for our algorithm. Experiments on various systems of difficulties demonstrate that our method always out- performs conventional gradient-descent based training, with resulting trained NeuralODEs having both higher interpolation and extrapolation capabilities than their counterparts. Especially, for the relatively simple Lotka- Volterra system, we report a x x102 improvement for interpolation error and staggering x x10 improvement in extrapolation error, showcasing the power of our new approach.
2 BACKGROUND
2.1 NEURAL ORDINARY DIFFERENTIAL EQUATIONS
A NeuralODE (Chen al. 2018 is model of the form,
du
(t, u; 0), u(l (0)= 40
di
where ug E Rn is the initial condition or input given to the model, and U(...: R R is neural network with parameters 0 E " that governs the dynamics of the model state E R" over time t E R. The value of the model state at a given time can then be evaluated by numerically integrating [ starting from the initial conditions.
In this paper, we concern ourselves with the problem of training NeuralODEs on time series data. Specifically, given an monotonically increasing sequence of time points {t(9) N Ji=o and the corre- sponding vector-valued measurements {ū (1) E IR" Jizo, we wish to train a NeuralODE on the data to learn the underlying governing equation and forecast future data.
Conventionally, NeuralODE training starts with using an ordinary differential equation (ODE) solver numerically integrate equation Io obtain the model state u at given time points:
N = ODESolve du = (t, u: 0), J. uA
uu (1
(2)
dt
with u (0) being a shorthand for ue ( 0). Afterwards, the loss function ((0) R + computed according to
R is
N N+1 Σ [(0==>-\>'m
1
E'((000)) (0) -0))
((
where ((u, ú) the pairwise loss function. this paper, we adopt the widely used mean-squared error function ((u, ú) I|u à ///, but other metries such as the LI loss can be used (Finzi al. 2020
Kim et al. 2021).
Training performed by minimizing equation 3 via gradient descent. non-trivial aspect of this process is that computing VoL requires differentiating the ODESolve operation. This can be done by either directly backpropagating through the internals of the ODE solver algorithm which re- turns accurate gradients but is memory intensive or by the "adjoint method", which computes an auxiliary set of ODEs to obtain gradients at a low memory cost, but can yield inaccurate gradients. In this paper, we embrace recent advances in the field and use the "symplectic-adjon method", which brings the best of both worlds by having both low memory footprint and improved accuracy guarantees (Matsubara al. 2021)
2

---

which has an additional term in the a dynamics that Under review as a conference paper at ICLR 2023
2.2 SYNCHRONIZATION OF DYNAMICAL SYSTEMS
Consider the following two dynamical systems.
du
du
dt (t,,u;))) dt (1,
(4)
where u, i E I" are n-dimensional state vectors, and U, Ü R - R are the dynamics of the two systems with coefficients 0, O. The functional form of the first system U....: is known and wish fit this to data generated from the second system, (ū("), Ji= at time points ( Ti-o- We assume U U, that is, the functional form of the data generating dynamics the same as the equation to be fitted and the initial conditions of the two systems are identical.
As the coefficients 0, of the two systems are different, the resulting trajectories u(t), ú(t) will, in general, be independent of each other. For chaotic systems, the result is more drastic as the solutions to these systems are extremely sensitive to parameter variations (Figure 1). For the periodic system (left panel), we find that +50% parameter perturbation only changes the period and amplitude of the signal. In contrast, perturbation of the same relative magnitude chaotic system (middle panel) results in completely different behaviors for the resulting solution.
20
MXX
25
150
Lon: 25
5-
100
Lum::
25
-20
Parameter perturbation (96)
Figure 1: Dynamics of single parameter perturbed systems. The systems and the parameters used are described in Experiments (Left) Solutions for the periodic Lotka- Volterra equations with perturbed a parameter. (Middle) Solutions for the chaotic Lorenz system with perturbed parameter. (Right) MSE loss function landscape for the Lotka- Volterra equations (Up) and Lorenz system (Down) for different lengths of time series data for the loss calculation.
Determining the coefficients è of the known equation to find the unknown coefficients Ô un- derlying the data (+) ) proceeds analogously to the NeuralODE training described in section (. the equation to be fitted is numerically solved for given parameter guess, the loss function is calculated with respect to data, then minimized using gradient descent to yield the final answer. However, the independent evolution of two uncoupled dynamics can easily lead non-convex loss functions with sharp local minima, especially for longer, more irregular time series (Figure IJ right panel). This, in turn, leads unstable loss function minimization and finally resulting in inaccurate coefficients (Voss et al. 2004).. Furthermore, this problem of irregular loss function during train- ing also occurs recurrent neural network (RNN) and NeuralODE training as well (Doya 1993 Ribeiro et al. 2020).
Enter a slightly altered version of equation couples the two systems:
du
du ÜL, úz
Ü ū; (1, K ④),
(5)
dt
di
where the K = diag(k).- k..) the diagonal coupling matrix K. It was found that if the elements of K positive and are sufficiently large, the dynamics of and v synchronize with increasing time regardless
the parameter mismatch: meaning [|ū(t) - ú(t)ll2 - as - oo (Abarbanel) et al. 2009 Pecora & Carroll 2013). Figure iillustrates this phenomenon for both periodic (left panel) and chaotic (middle panel) systems, using = kI where k E R+ is the scalar coupling identity matrix.
strength and I is the
The new coupling term acts similarly to the proportional control term used in PID control, and actively drives the modified dynamics ūit) towards the reference trajectory u(l). Increasing the
3

---

Under review as a conference paper at ICLR 2023
k-C.0
20
0.0
2 KI
0-
1.0
100
0.0
0.3
25
20
20
40
10
Parameter perturbattin
Figure Dynamies of coupled systems with varying coupling strengths. The systems and the parameters used are identical to those of figure [ (Left) Results for the periodic I Lotka-Voltery equations. (Middle) Results for the chaotic Lorenz system. (Right) Loss function landscape of the coupled systems with different coupling strengths.
coupling strength correspond to turning up the proportional gain, and causes the two trajectories to converge even faster (Figure 2] left and middle panels).
Repeating the parameter estimation procedure, the loss function L between ū and á now becomes function of both 0 and K (= kT) (Figure 2] right panel). In the context of loss function calcula- tion, increasing the coupling strength k effectively acts as shortening the length of the data because increased coupling causes the two trajectories start synchronizing earlier time, after which the error between the two will diminish and contribute negligibly to the total loss. Therefore, synchro- nizing the model and the data generating equations leads to more favorable loss landscape, which gradient descent based optimizers will not struggle on.
However, directly minimizing this new loss function C(O, K) associated with the coupled system (equation 5) does not result in our desired goal of discovering the unknown data generating coeffi- cients O. This is because goal requires to fit the original uncoupled equation 44) on the data, and not the modified equation (5) that has an additional term. Therefore, a method that reaps the benefits of the original loss function C(0) is needed. That method is homotopy optimization, which we describe in the following section.
behaaed loss function using synchronization, but also results in minimizer of
2.3 HOMOTOPY OPTIMIZATION AND PARAMETER IDENTIFICATION
Consider the problem of finding the parameter vector 0* R that minimizes the function F(0): RT + F. While minimization algorithms, such as the commonly used gradient descent methods, can directly be applied to the target function, these tend to struggle if F is non-convex and has a complicated landscape riddled with local minima. To address such difficulties, the homotopy optimization method (Dunlavy & O' O'Leary 2005) introduces alternative objective function
{
9(0). is = F(O), if A =
H(O, A)
(6)
where ((0) : RT + R is an auxillary function whose minimum easily found, and H(0, A) R. x R + R is continuous function that smoothly interpolates between G and F as the homotopy parameter À is varied from to
The motivation behind this scheme is similar simulated annealing: starts out with relaxed version of the more complicated problem of interest, and finds series of approximate solutions while slowly morphing the relaxed problem back into its original non-trivial form. This allows the optimization process to not get stuck in spurious sharp minima and accurately converge to the minimum of interest.
To proceed with the method, one first selects the number of discrete steps for optimization, 8, as well as a series of positive decrement values for the homotopy parameter (A)( (19))-- that sum to 1. Afterwards, optimization starts with an initial À value of ( = 1, which gives H()))) = ( At each step, the objective function at the current iteration is minimized with respect to 0, using the
4

---

Under review as a conference paper at ICLR 2023
output from the previous step (*111) is the initial guess, and
10 (0) M(e, À ((h)) - ge(k) arg min (() (0)
Afterwards, À decremented to its next value, ( = A A and this iteration continues
until the final step s where ( = 0. H (0) - F(0). and the final minimizer (*(s) 0 the
sought-after solution to the original minimization problem F(0).
Application to synchronization. To combine homotopy optimization with synchronization, we slightly modify the coupling term of equation 5 by multiplying it with the homotopy parameter À:
K K
-AK û).
(8)
With this modification, applying homotopy optimation to the problem is straightforward. When = and the coupling matrix K has sufficiently large elements, synchronization occurs and the resulting loss function L(O, K) is well-behaved, serving the role of the auxillary function G in equation 6 When = 0, the coupled equation [3] reduces to the original equation and the corresponding loss function C(0) C(0.0 K) is the complicated loss function JF need to ultimately minimize. Therefore, starting with X = and successively decreasing its value to 0 in discrete steps, all the while optimizing for the coupled loss function L(O, XK) allows the well-behaved loss function landscape from synchronization while being able to properly uncover the system parameters (Vyasarayani et al. 2012: Schäfer al. 2019).
leverage
HOMOTOPY OPTIMIZATION FOR NEURALODE TRAINING
While the approaches in the previous sections were developed for the problem of determining the unknown coefficients of a differential equation from data, we propose to apply them to the problem of training NeuralODE on time series data. This conceptual leap based on two observations:
1. The previous formalism synchronization and homotopy optimization can directly be translated over to NeuralODE training by reinterpreting (..., 0) from 'a known equation whose coefficients 0 we wish to find' 'a NeuralODE with parameters 0 we wish learn the data generating equation Ü..,, 0) with.
2. While (ubsection 2.2 considered the case where the model equation differs only in coeffi- cients with the data generating equation (U Ú), has been discovered that two systems with different functional forms can also synchronize phenomena called "generalized synchronization" (Pecora & Carroll (0155). This motivates our coupling a NeuralODE to given time series data to facilitate training.
3.1 IMPLEMENTATION DETAILS
The previous sections describe a general methodology of using homotopy optimization for Neu- raIODE training. However, there are couple details that need addressing in order to create a specific implementation of our algorithm. We briefly discuss these points below.
Construction of the coupling term In the previous section, the coupling term KK( ú) in the right hand side of equation 5 was implicitly thought to be defined for all time t. However, Neu- ralODEs are trained on measurements ( === sampled at discrete time points {e(0) ===o which makes the control term undefined on any other time points. Leaving the control term partially defined is not option, most ODE solver algorithms require evaluating the right hand side of equation at intermediate time points as well. Therefore, one must extend the definition of the control term so that it is defined even at times where the measurements ()) are not supplied.
One solution to this problem is to apply the coupling term sequence of impulses, only applying it where it is defined, and keeping it zero at other times. While this an interesting avenue of interest especially so since this be shown to be mathematically equivalent to teacher
et al.] 2009), method used to mitigate instabilities during RNN training(Toomaran Toomarian & Barhen do not adopt this method in this paper. Instead, construct smoothing cubic splines
1992
5

---

Under review as a conference paper at ICLR 2023
from the data points to supply the measurement values at unseen time points that needed to evaluate the coupling term, with smoothing chosen over interpolation to make our approach robust to measurement noise in the data.
Scheduling the homotopy parameter While the homotopy parameter À can be decremented in various different ways, the most common approach is to use constant decrements: A 1/s, E (0, 8 - 1] where is the number of homotopy steps. In our study, we modify this slightly so that the decrements are decayed by constant ratio for each passing step (Figure 3) left panel). that is,
AX =
Σs =1 =1; D<NASL
k=
3.2
ALGORITHM OVERVIEW
Our final implementation of the homotopy training algorithm has five hyperparameters. Here, we briefly describe the effects and the tips for tuning each.
Number of homotopy steps (s) This determines how many relaxed problem the opti- mization process will pass through to get to the final solution. Similar to scheduling the temperature in simulated annealing, fewer steps results in the model becoming stuck in sharp local minima, and too many steps makes the optimization process unnecessarily long. We find using values in the range of 6-8 or slightly larger values for more complex systems yields satisfactory results.
Epochs per homotopy steps (Nepoch): This determines how long the model will train on given homotopy parameter value A. Too small, and the model lacks the time properly converge on the loss function; too large, and the model overfits on the simpler loss land- scape of + 0, resulting reduced final performance when O. We find for simpler monotonic or periodic systems, values of 100-150 work well; for more irregular systems, 200-250 are suitable.
Coupling strength (k) : This determines how trivial the auxillary function for the homo- topy optimization will be. Too small, and even the auxillary function will have jagged landscape; too large, and the initial auxillary function will become flat (Figure 2] left panel, k = 1.0, 1.5) resulting in very slow parameter updates. We find good choices of k tend to be comparable to the scale of the measurement values.
A decrement ratio (kx) This determines how the homotopy parameter À is decremented. Values close cause À decrease in nearly equal decrements, whereas smaller values cause a large decrease of À in the earlier parts of the training, followed by subtler decre- ments later on. We empirically find that KY values of 0.5-0.6 tends work well.
Learning rate (n) This is as same as conventional NeuralODE training. We found values in the range of 0.01-0.1 to be adequate for our experiments.
EXPERIMENTS
We evaluate the performance of our algorithm by training NeuralODEs on three different systems of varying difficulties selected from literature: the Lotka-Volterny system, the Lorenz equations and the double pendulum. In both homotopy and baseline experiments, the MSE loss was monitored during training, and the checkpoint corresponding to this MSE minimum was used for further anal- ysis. Each experiment was repeated three times with different random seeds, and we report the corresponding means and standard emrors solid lines and shaded bands wherever applicable. Ad- ditional details regarding the experiments can be found in the ??
4.1
SYSTEM
We start our experiments with the Lotka- Volterra system, which is a simplified model predator- prey dynamics given by the following equations,

---

Under review as a conference paper at ICLR 2023
dx di
dy
az - Bry,
di - -y + dxy
(10)
where xff), vt) denote the prey and predator populations as function of time, and a, 3, 7,8 are parameters regarding the interaction between the populations (Murray) (002). We follow the exper- imental design of Rackauckas et al. (202I), where the goal of the experiment is to fit the following hybrid model on the simulated data,
dæ dt
dy
or + 1;
= -y + U2(2, yū 02)
(11)
dt
where U1 y; 01). (2 y; 02) are neural networks employed to discover the corresponding "un- known" terms - Srz, dxy from the given data.
Train loss LIR.A1
20
Homatopy Vanils floor
7.5
Date Homotogy
Venill Extrap
Mean squared eneer
2.5
10-1
0.0
10a
10-2
State
cormyttyyyyy param.
10-1
2-
200
400 Fpoch
600
800
200
400
600
800
0
10
15
Epoch
Time
Figure 3: Lotka-Volterra dataset results. (Left) Train loss L(B, A) and homotopy parameter À during training. (Middle) Mean squared error as a function of training epochs for both homotopy-based and vanilla NeuralODE. Models corresponding to the minimum of this curved are used plot the prediction results on the right. (Right) Interpolation and extrapolation results for both the homotopy and vanilla trained models. The dashed gray line indicates the end of training data and the start of extrapolation.
Figure displays results for both vanilla and our homotopy- NeuralODE training. We find that not only does our training method reach the final MSE value of the vanilla Neural in about half the number of epochs, the best MSE error achieved is near to the noise floor of the dataset which is the lowest value the MSE can take without model overfitting. Inspecting Figure 3 (right), we find NeuralODEs trained with the homotopy method is able perfectly fit the training data in the interpolation regime as suggested by extremely low best MSE value of the middle panel. The high performance of the trained model continues far into the extrapolation regime (after the dashed gray line), where it near perfectly predicts the data across multiple future periods. As the final trained model was selected only using the training data and had not seen any data in the extrapolation regime, the perfect future predictions indicate that our training method accurately extracted the underlying dynamics of the data. These results are in stark contrast to those obtained by standard NeuralODE training, whose results overfit the y variable and neglect the ø variable almost entirely and has no predictive capabilities.
Training on different length trajectories. Simply training NeuralODEs with gradient-based optimizer gives drastically deteriorating results as the length of training trajectories increases. To illustrate this point, simulation trajectories of different time spans were generated: t [0, 3.1, [0,6 6.1, [0.9.1] A/ 0.1, and the same hybrid model (equation [T] was trained using both vanilla and homotopy methods with identical hyperparameter settings as before.
Figure 4 shows the training results for Lotka- Volterra train data of different lengths. In the case of interpolation, our method consistently trains models to loss values near the noise floor, regardless of the length of the training data. On the other hands, the conventional training method only manages to train NeuralODEs for very short data, and completely fails to learn the oscillatory dynamics that start arise at longer times. The differences are even more drastic for the case of extrapolation, where there is a x 107 MSE difference between our approach and the baseline for longer time series. In the case of short data, our result also gives high, albeit slightly lower loss value than the baseline.
7

---

Under review as a conference paper at ICLR 2023
10'
Date
Mean squared error
10
length: lwngth:
Noise floer
Humotopy (intep))
Varilla Griepp
10
Homotopy textrap))
Verilla (estraps)
Stzte
10-1
40
60
80
length
Time
Figure 4: Different train data length results. (Leftt) Interpolation and extrapolation MSE values as a function of training data length. Extrapolation MSE was calculated by computing prediction values for 50 points after the end of training data, then computing the MSE value against the noise-added true dynamics. Error bars indicate one standard error away from the mean over three runs. (Right) Interpolation and extrapolation results for models trained with our method. End of training data is marked by the colored dash lines, and extrapolations are performed for 50 points after the end of train data to match the results on the left panel.
Inspecting the right panel, this revealed be due lack of information in the training data given access to less than half the oscillation period, it would be impossible for the model to determine whether the true dynamics is oscillatory or not.
4.2 LORENZ EQUATIONS
For our next set of experiments, we turn to the Lorenz equations, given by
dx dt = -o - x),
dy
dz
dt = =)
zy Bz.
(12)
dt
The chaotic nature of these equations makes even the seemingly simple task of determining the un- known coefficients of the equations from the data extremely nontrivial. Hence, the Lorenz equations serve as the perfect stress test for training NeuralODE on the data.
The task of the experiment was to fit the "black-box" model of equation on the simulated and noise added trajectory of 31 points, and the resulting extrapolation performance of the trained models were inspected using another 31 points following the end of training data.
Data comatapy
Statex
Venill
Homptogy
Extrap.
40
Statey
10
10
-10
x
Figure 5: Lorenz system: (Left) State space trajectories for the train data and learned dynamics using the homotopy and vanilla training methods. (Right) Interpolation and extrapolation results on the Lorenz equations dataset for both homotopy and vanilla trained results.
Figure 5) agaain highlights the effectiveness of our homotopy training against the baseline approach. Without homotopy, the NeuralODE fails to learn any useful feature of the system dynamics, col- lapsing into a flat line after the first few points in the interpolation interval. In contrast, results from the homotopy training is able to accurately fit the training data in the interpolation interval, having learnt both the initial peak and the subsequent oscillations. Due to the complexity of the chaotic system, the model predictions does start to go oft-of-syne with the ground truth in the extrapolation regime, but it is evident the dynamics learned with homotopy is much richer and accurate than its baseline counterpart.
8

---

Under review as a conference paper at ICLR 2023
4.3 DOUBLE PENDULUM
For our final set of experiments, we used real-world measurement data of the dynamics of a dou- ble pendulum taken by Schmidt & Lipson (2009). The goal of the experiment was fit another "black-box" model on the training data, which consists of 100 time points, and monitor the model predictions for another 50 points afterwards.
State
10
20-
State 64g
State 622
-10
0.00
0.25
0.50
0.75
1.25
1.50
0.75 1.00 1.25 Timec
Time
Figure 6: Double pendulum: (Left) Diagram of the double pendulum with the state variables 0,, 02, Wi, w2 denoted. 1 and 2/ correspond to the lengths of the pendulum bars. (Right) Interpolation and extrapolation results on the double pendulum dataset for both homotopy and vanilla trained results.
From Figure 6 we see once again that our homotopy approach produces models that are much more accurate the interpolation regime than those by the conventional training method. For extrapola- tion, both homotopy and baseline results fail to accurately predict all four state variables. However, the homotopy result does much closely follow the rises and falls of the ground truth data (especially for 01) than its counterpart, demonstrating the effectiveness of our algorithm.
5 RELATED WORKS
Stabilizing NeuralODE training The necessity for improving NeuralODE training has been no- ticed by the community, and multiple different methods have been proposed. Works such as Choro- manski et al. (2020) or Hasani et al. (202T) use specific mathematical forms for their NeuralODEs that guarantee stability. However, for the purpose of discovering dynamical equations underlying the data, these approaches are inadequate as one cannot ascertain their hidden equations conform to the rigid mathematical structure of these approaches. Another class of methods, such as Finlay al. (2200) and Kelly et al. (2020) use softer form of regularization to boost training, by placing constraints on the Jacobian or high derivatives of the differential equation to be learned. However, these methods were also developed in the context of using NeuralODEs for input-output mapping. There, as long as the mapping is the same, the actual form of the dynamics underneath does not matter which is not the case in our problem setting. Ghosh et al.) (2020) introduces method that does not impose any constraint on the model form, but treats the integration time points as stochas- tic variables. Other methods such Kidger et al. (2021), ZZhangg et al. (2020), Matsubara et al. (2221), and Kim et al. (202T) study other aspects of NeuralODE training, including better gradient calculation using improved adjoint methods. Such methods are fully compatible with our proposed work, as evidenced by our use of the 'symplectic-adjoii" adjoint-method" method" of Matsubara et al. (2021) in our experiments.
DISCUSSION AND OUTLOOK
In this paper, we adapted the concepts of synchronization and homotopy optimization for the first time in the NeuralODE literature, and demonstrated that models trained with our proposed method are unrivaled in both interpolation and extrapolation accuracy across three different systems. As our training algorithm is first of its kind for NeuralODEs, there are diverse avenues of future search, including a more through investigation on the presented hyperparameter selection heuristics, using regularization terms in loss function decrease the coupling term opposed to manual scheduling (Abarbanel et al.) 2009) Schäfer et al. 2019). As NeuralODEs are continuous analogs of RNNs, work holds new prospect for stabilizing RNN training well.

---

Under review as a conference paper at ICLR 2023
REFERENCES
Henry D. L. Abarbanel, Daniel R. Creveling, Reza Farsian, and Mark Kostuk. Dynamical state and parameter estimation. SIAM Journal on Applied Dynamical Systems, 8: 8:1341-1381, 2009.
Steven L.. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences, 113:3932-3937, 2016.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K. Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems, volume 31, 2018.
Krzysztof M Choromanski, Jared Quincy Davis, Valerii Likhosherstov, Xingyou Song, Jean-Jacques Slotine, Jacob Varley, Honglak Lee, Adrian Weller, and Vikas Sindhwani. Ode to ode. In Advances in Neural Information Processing Systems, volume 33, 2020.
Kenji Doya. Bifurcations of recurrent neural networks in gradient descent learning. IEEE Transac- fions on neural networks, 1:218, 1993.
Daniel M. Dunlavy and Dianne O'Leary. Homotopy optimization methods for global optimiza- tion. Technical report, Sandia National Laboratories (SNL), Albuquerque, NM and Livermore, CA, 2005.
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, volume 37, 2020.
Mare Finzi, Ke Alexander Wang, and Andrew G. Wilson. Simplifying hamiltonian and lagrangian neural networks via explicit constraints. In Advances Neural Information Processing Systems, volume 33, 2020.
Arnab Ghosh, Harkirat Behl, Emilien Dupont, Philip Torr, and Vinay Namboodiri. Steer: Simple temporal regularization for neural ode. In Advances in Neural Information Processing Systems, volume 33, 2020.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In Advances in Neural Information Processing Systems, volume 32, 2019.
Ramin Hasani, Mathias Lechner, Alexander Amini, Daniela Rus, and Radu Grosu. Liquid time- constant networks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, 2021.
Jacob Kelly, Jesse Bettencourt, Matthew J Johnson, and David K Duvenaud. Learning differen- tial equations that are easy to solve. In Advances in Neural Information Processing Systems, volume 33, 2020.
Patrick Kidger, Ricky TQ Chen, and Terry J Lyons. " hey, that's not an ode": Faster ode adjoints via seminorms. In International conference machine learning. volume 38, 2021.
Suyong Kim, Weiqi Ji, Sili Deng, Yingbo Ma, and Rackauckas Christopher. Stiff neural ordinary differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31:093122, 2021.
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. arXiv preprint, arXiv: 1711.05101, 2017.
Takashi Matsubara, Yuto Miyatake, and Takaharu Yaguchi. Symplectic adjoint method for exact gradient of neural ode with minimal memory. In Advances in Neural Information Processing Systems, volume 34, 2021.
James D. Murray. Mathematical biology, volume 1. Springer, 2002.
Louis M. Pecora and Thomas L. Carroll. Synchronization of chaotic systems. Chaos: An Interdis- ciplinary Journal of Nonlinear Science, 25:097611, 2015.
10

---

Under review as a conference paper at ICLR 2023
John C. Quinn, Paul H. Bryant, Daniel R. Creveling, Sallee R. Klein, and Henry D. L. Abarbanel. Parameter and state estimation experimental chaotic systems using synchronization. Physical Review E, 80:016201, 2009.
Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner. Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal differential equations for scien- tific machine learning. arXiv preprint, arXiv:2001, .04385, 2021.
M. Raissi, P. Perdikaris, and G. Karniadakis. Multistep neural networks for data-driven discovery of nonlinear dynamical systems. arXiv preprint, arXiv: 1801.01236, 2018
António H Ribeiro, Koen Tiels, Luis A Aguirre, and Thomas Schön. Beyond exploding and vanish- ing gradients: analysing rnn training using attractors and smoothness. In International Conference on Artificial Intelligence and Statistics, volume 23, 2020.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. Science, 324:81-88 2009,
Kai Schäfer, Kathrin Flaßkamp, Jörg Fliege, and Christof Büskens. A combined homotopy- optimization approach to parameter identification for dynamical systems. In Proceedings Ap- plied Mathematics and Mechanics, volume 19, 2019.
Nikzad Benny Toomarian and Jacob Barhen. Learning trajectory using adjoint functions and teacher forcing. Neural networks, 5, 1992.
Henning U. Voss, Jens Timmer, and Jürgen Kurths. Nonlinear dynamical system identification from uncertain and indirect measurements. International Journal of Bifurcation and Chaos, 14:1905- 1933, 2004.
C. Vyasarayani, Thomas Uchida, and John McPhee. Single-shooting homotopy method for pa- rameter identification in dynamical systems. Physical Review E, 85:036201, 2012.
Juntang Zhuang, Nicha Dvornek, Xiaoxiao Li, Sekhar Tatikonda, Xenophon Papademetris, and James Duncan. Adaptive checkpoint adjoint method for gradient estimation in neural ode. In International Conference on Machine Learning, volume 37, 2020.
A EXPERIMENT DETAILS
In all experiments, the AdamW optimizer (Loshchilov & Hutter 2017) was used to minimize the re- spective loss functions for the vanilla and homotopy training. The number hidden layers was fixed to two for all experiments, with hyperbolic tangent activation for the input and hidden layers, and the identity function for the output layer. Hyperparameters were chosen by running hyperparameter sweeps and selecting the resulting best values.
A.1 LOTKA-VOLTERKY SYSTEM
Following the experimental design of Rackauckas al. (20TT) we numerically integrate equation TO in the time interval t [0, 6.1]. At 0.1 with the parameter values B, 7, 5 = 1.3, 0.9, 0.8, 1.8 and initial conditions (0), (0) 0.44249296, 4. 6.6280594 Continuing with the recipe, Gaussian random noise with zero mean and standard deviations with magnitude of 5% of the mean of each tra- jectory was added to both populations. The integration was performed using the dopri5 solver from the syyymppeecc-adjoiiiii package (Matsubara et al. 2021 was used with absolute tolerance of le-9 and a relative tolerance le-7.
The neural networks Ui (x, 0,) and U2(2, y/; 02) had 2 nodes for the input layer, 5 nodes for its two hidden layers, and node for the output layer. Deviating from Rackauckas et al. (202T), used hyperbolic tangent activation functions instead of the radial basis function activation that the previous authors used. This change was made to better adhere to the standard practices of the NeuralODE community.
11

---

Under review as a conference paper at ICLR 2023
A.2 LORENZ ATTRACTOR
To generate the data, we followed experimental settings of Vyasarayani et al. 2012). we used the parameter values o, p. 3 10, 28, 8/3 and the initial condition Ho, 90; =0 1.2, 2.1, 1.7 which, upon integration, gives rise to the well-known "butterfly attractor". Adhering to the paper, a Gaus- sian noise of mean 0 and standard deviation 0.25 was added to the simulated trajectories to model experimental noise. The ODE solver and the tolerance values were kept identical to the previous Lotka- Lotka-Volterra experiment.
The NeuralODE used for the experiment had 50 nodes for its two hidden layers, and 3 nodes for its input and output layers corresponding the three degrees of freedom of the state vector [z, y, 2]" of the system.
A.3 DOUBLE PENDULUM
The experimental data from Schmidt & Lipson (2099) consists of two trajectories of the double pendulum, captured using multiple cameras. The noise the data is subdued due to the LOESS smoothing performed by the original authors. For experiments, used the first 100 points of the first trajectory for training and the next 50 for evaluating the extrapolation capabilities of the trained model.
Similar to the previous experiment, we used a "black-box" NeuralODE with 50 nodes for the hidden layers, and the input and output layer nodes changed to 4 to reflect the four degrees of freedom of the system (Figure 6, left panel).
12


Reviews:
[{"review_id": "zEe-x73Tv5", "reviewer": "Reviewer_4PH8", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The manuscript proposes to use homotopy optimization between a stabilized version of a neural ODE and the original neural ODE to learn dynamical systems from time series data. The proposed method is tested on three synthetic problems, where it is shown to outperform a naive gradient descent training of the neural ODE.", "summary_of_the_review": "The manuscript addresses the important and timely problem of learning dynamical systems from time series data using neural ODEs. \nI appreciate the idea of using a homotopy to a sychronized system in order to leverage training of neural ODEs. \nCurrently, I have the a few main concerns regarding the manuscript in its current form:\n* The idea of homotopy optimization with a synchronized system is not new. However, prior works from the parameter identification literature are only mentioned at the very end of Subsection 2.3 and neither in the introduction nor in the related works section. In general I would describe the contribution of the paper rather the application of this approach in the context of neural ODE training as they fall in the abstract setting of parameter identification as far as I can see it (see for example Schäfer et. al. 2019). That being said, I don’t mind applying an idea from another field, but I think the contribution of the manuscript in comparison to prior works have to be sufficiently clear.\n* Being an empirical paper, I am currently not convinced that the experiments are comprehensive enough. In particular, the proposed method is only compare against one baseline, which is vanilla gradient descent.  \n* The overall presentation can be improved. \n\nOverall, I believe that the manuscript is not ready for publication in its current form. Nevertheless, I want to highlight that I find the idea of using homotopy optimization for neural ODE training intriguing and believe that with stronger experiments, this can be an interesting direction.\n", "strength_and_weaknesses": "**Strengths:**\n* The idea of using a homotopy to a problem with a more regular landscape to leverage difficulties in the training of neural ODEs is very attractive. Also the demonstration that such a problem can be constructed via synchronization is valuable.\n\n**Weaknesses:**\n* Apart from being an attractive idea, homotopy continuation with synchronization has been used for parameter identification in dynamical systems before. Hence, the contribution of the manuscript reduces to the application of an existing method for parameter identification to the special case of neural ODE training.\n* Being a purely empirical work, the experiments are not very strong. In particular, the proposed method is only compared against one baseline, which is the vanilla gradient descent training of the neural ODE.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "1: The contributions are neither significant nor novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "iTMC5LAAJZ", "reviewer": "Reviewer_3MeW", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a new training paradigm for Neural ODEs by annealing a coupling term between the learnt system and the target system. This coupling allows for a more efficient training of the Neural ODE by simplifying the loss landscape. The authors then showcase their approach on synthetic data (Lorenz, LV and double pendulum). ", "summary_of_the_review": "Interesting and promising idea for improving the training of neural odes. Nevertheless, the paper currently lacks depth in the investigation. In particular, theoretical motivation for the more behaved loss and/or more comprehensive experiments. In its current form, the paper is nor theoretically nor experimentally convincing but I fully believe there is the potential for it to be in the future.", "strength_and_weaknesses": "### Strengths\n\n- The coupling of the dynamics to improve the training is a very nice and promising idea. \n- The paper is easy to read and understand\n\n### Weaknesses\n\nIn general, I think this paper does not go deep enough either theoretically or practically, therefore missing the impact it could have.\n\n- I believe it would be very interesting to have a more in-depth investigation of the properties of the loss including the coupling. You show graphically that this coupling results in a smoother loss landscape but your work requires more theoretical foundations to motivate the approach in general.  You state that previous works have found that it synchronizes but you should do the extra step of deriving what it implies for the loss.\n\n- From an experimental point of view, it would be desirable to have some real world data as well. It's ok if the theoretical motivation is not strong but then you need a solid experimental setup. Additionally, it seems you are using only single realization of time series. How does your model transfer to the case of multiple realization of a time series (say different clinical trajectories for instance) ?\n\n- In terms of limitations, the fact that coupling is required at every time step is significant. I think it would be interesting to know when the spline interpolation fails, when the sampling rate because lower and lower.  In particular, the sampling rate used in your experiments seems quite high. \n\n- The mean square errors you get from the vanilla Neural ODEs in Figure 4 seem unrealistic and might point towards some regularization issue in the Neural ODE.  It's not clear from your experiments how you used regularization for training. I think it should be fair to have an extensive hyperparameters search to compare against your method. Also, why are the results so bad for the Lokta Voltera but seems to be acceptable for Lorenz and Double Pendulum which are arguably much more difficult systems ?\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "PGQCoGQTh_7", "reviewer": "Reviewer_89pi", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a new method for training Neural ODEs by employing techniques of synchronization and homotopy optimization. In the task of learning physical systems, the authors have shown that it is possible to learn appropriately for long duration data.", "summary_of_the_review": "The proposed method is novel, but its validity and effectiveness are somewhat unclear.", "strength_and_weaknesses": "S1. This paper is the first example of introducing synchronization and homotopy optimization for learning Neural ODEs.\nS2. The effectiveness of the proposed method is verified for multiple physical systems.\n\nW1. The proposed method is novel, however, it is somewhat unclear why the introduction of synchronization and homotopy optimization solves the problem. First of all, can you briefly explain why the learning becomes unstable for long duration data?　Then, can you explain how the proposed method solves the problem? An intuitive introduction would also be helpful.\nW2. As the authors said, It is a well-known fact that learning is unstable for long duration data. For practical use, a time-window is introduced and a shortened sequence is used as a mini-batch for training. What are the advantages of the proposed method over such techniques?\nW3. Can the proposed method be applied to machine learning tasks (regression, classification, etc.)? If so, it would be good to discuss its effectiveness.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "9yTK21i39g", "reviewer": "Reviewer_zjmg", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes a technique for matching the predicted trajectory to the trained ODE using a Coupling term. It also introduces the idea of Homotopy training, where we have the interpolation between a smooth, easy function and the target function that we want to optimise. Using these two techniques allow to better fit the Neural ODE to the target trajectory, including Lotka-Volterra models, Lorenz attractor and double pendulum. Using these techniques also lead to the better extrapolation of the predicted trajectory.", "summary_of_the_review": "The paper introduces an interesting technique to enforce the match between the ground-truth trajectory and the. Trained ODE, as well a simulated-annealing-like technique to start the optimisation from an easier function.\n\nHowever, the results of the baseline (Neural ODE) are not consistent with the results from the previous papers and external sources, suggesting that it was not trained properly or under-parameterised. I am giving the weak-reject, but willing to greatly increase the score if the authors re-train the Neural ODE model on Lotka-Volterra and Lorenz attractor and show the better fit than just predicting a constant function, or explain why their Neural ODE results might be inconsistent with the previous works.\n", "strength_and_weaknesses": "**Strengths**\n\nBoth techniques introduced in the paper come from the classic optimisation theory. Intuitively, the idea is similar to simulated annealing — first optimise an easier function that is already close to the ground-truth, and then increase the difficulty of the optimised task.\n\nThe paper demonstrates the approach on notoriously challenging systems, such as Lorenz attractor and double pendulum. Since Lorenz attractor is a chaotic system, it requires a very accurate reconstruction of the trajectory and is hard to extrapolate correctly. The coupling term ensures that the trained ODE can accurately fit the target trajectory. Figure 2 also provides an insight on how the optimised function changes because of the coupling term.\n\n**Weaknesses**\n\n1. My main concern is the evaluation of the Neural ODE baselines. It is surprising that Neural ODE has extremely poor fit on Lotka-Volterra model in Figure 3 and Lorenz attractor on Figure 5, but can fit the double pendulum reasonably well (Figure 6).  I trust that the coupling term and the homotopy would help the ODE training (mostly based on results in Figure 6), but it is hard to evaluate the benefits of the proposed approach, if the baseline is not well trained.\n\n\n    To my understanding, the training is done on one trajectory, so Neural ODE should be able to overfit the training data, or at least have a better fit than just a constant function. For example, in the original paper the model was able to fit the oscillating data: https://arxiv.org/abs/1907.03907. This is the example where Neural ODE is able to fit the trajectory from Lotka-Volterra model: https://gist.github.com/ChrisRackauckas/a531030dc7ea5c96179c0f5f25de9979\n\n    Figures 3 and 5 suggest that the Neural ODE baseline is either under-parameterized, or not trained long enough.\n\n\n    Can the authors provide a similar experiment to figure 3, where the Neural ODE is able to fit at least the training part of the trajectory, or explain why they think the results are not consistent with the results of the previous papers?\n\n\n2. The proposed approach uses splines to impute the input data to construct the coupling term. The splines can provide the reasonable interpolation on densely-sampled data with low amounts of noise, but might require a more careful consideration in case of sparse data or data with higher amount of noise (a common scenario in the real-world data). However, this is not a big limitation, because the splines are only used in the coupling term, which is set to zero in the end of the training.\n\n**Minor**\nSection 2.1: “Differentiating the ODESolve operation … can be done by the ”adjoint method”, which …. can yield inaccurate gradients.”\n\nCan the authors give examples or citations when the adjoint gradients can be inaccurate, but the ”symplectic-adjoint” integrators are accurate?", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "3UHIR8GfBY", "reviewer": "Reviewer_nRZJ", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors propose a novel training strategy for neural ordinary differential equations (NODEs) by introducing a coupling term between the true dynamics and NODEs. They formulate the coupling-based training framework by using the homotopy optimization, which optimizes a homotopy between the simple coupled loss landscape and the complicated uncoupled loss landscape. It results in stable and fast training of NODEs with respect to the targeted dynamical system. The authors validate their proposed learning framework for three standard benchmarks (Lotka-Volterra, Lorenz, and double pendulum) that frequently used in NODE literatures.", "summary_of_the_review": "I think the paper is interesting, as I mentioned above. I would like to vote to (weak) accept for this paper.", "strength_and_weaknesses": "Strength:\n1. The motivation of the introduced coupling term for synchronizing NODEs and true dynamics is clear and principled. The authors convincely explain the effectiveness of such a coupling term when estimating the optimal parameters of dynamical systems.\n\n2. The coupling-based synchronization cannot be used directly for training NODEs, because it gives an undesired model that depends on the unknown true dynamics. To overcome this issue, the authors introduce a homotopy between the coupled and uncoupled NODE loss functions. To me, the application of the homotopy optimization for training NODEs with the synchronization is novel.\n\n3. The paper is very well-written and easy-to-follow. \n\nWeakness:\n\n1. The proposed method has four additional hyper-parameters that should be tuned carefully. Although the authors provide rough guidelines for tuning such hyper-parameters, it is not entirely clear whether the proposed method is sufficiently robust and general.\n\n2. The experimental results on Lorenz and double pendulum are not sufficiently strong. There are some extrapolation errors. Considering that they are chaotic, it is acceptable (for me). For the double pendulum experiment, it might be interesting to apply the proposed homotopy optimization to Hamiltonian NODEs.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject