Homotopy-based training of NeuralODEs for accurate dynamics discovery
Authors: Joon-Hyuk Ko
Hankyul Koh
Nojun Park
Wonho Jhe
Venue: NeurIPS 2023 poster (2023)
Keywords: neural ordinary differential equations
synchronization
homotopy optimization
loss landscape
dynamical systems
TL;DR: We present a new NeuralODE training method based on synchronization and homotopy optimization that boosts the training of NeuralODEs without introducing any architectural constraints.

Abstract: Neural Ordinary Differential Equations (NeuralODEs) present an attractive way to extract dynamical laws from time series data, as they bridge neural networks with the differential equation-based modeling paradigm of the physical sciences. However, these models often display long training times and suboptimal results, especially for longer duration data. While a common strategy in the literature imposes strong constraints to the NeuralODE architecture to inherently promote stable model dynamics, such methods are ill-suited for dynamics discovery as the unknown governing equation is not guaranteed to satisfy the assumed constraints. In this paper, we develop a new training method for NeuralODEs, based on synchronization and homotopy optimization, that does not require changes to the model architecture. We show that synchronizing the model dynamics and the training data tames the originally irregular loss landscape, which homotopy optimization can then leverage to enhance training. Through benchmark experiments, we demonstrate our method achieves competitive or better training loss while often requiring less than half the number of training epochs compared to other model-agnostic techniques. Furthermore, models trained with our method display better extrapolation capabilities, highlighting the effectiveness of our method.

Homotopy-based training of NeuralODEs for accurate dynamics discovery
Joon-Hyuk Ko
Hankyul Koht
Department of Physics & Astronomy Seoul National University
Department of Physics & Astronomy Seoul National University
Seoul, 08826, South Korea
Seoul, 08826, South Korea
jhko7250smu. ac kr
physics1130es ac
Nojun Park
Wonho Jhe
Department of Physics Massachusetts Institute of Technology
Department of Physics & Astronomy Seoul National University Scoul, 08826, South Korea whjhe@snu. ac kr
MA, 02142, United States baj115260mit.. edu
Abstract
Neural Ordinary Differential Equations (NeuralODEs) present attractive way to extract dynamical laws from time series data, as they bridge neural networks with the differential equation- based modeling paradigm of the physical sciences. However, these models often display long training times and suboptimal results, especially for longer duration data. While common strategy in the literature imposes strong constraints to the NeuralODE architecture to inherently promote stable model dynamics, such methods the unknown governing equation is not guaranteed to satisfy the assumed con- straints. this paper, we develop new training method for NeuralODEs, based on synchronization and homotopy optimization, that does not require changes to the model architecture. We show that synchronizing the model dynamics and the training data tames the originally irregular loss landscape, which homotopy opti- mization can then leverage enhance training. Through benchmark experiments, we demonstrate our method achieves competitive or better training loss while often requiring less than half the number of training epochs compared to other model-agnostic techniques. Furthermore, models trained with our method display better extrapolation capabilities, highlighting the effectiveness of our method.
ill-suited for dynamics discovery
Introduction
Predicting the evolution of time varying system and discovering mathematical models that govern it is paramount to both deeper scientific understanding and potential engincering applications. The either ingeniously deduce empirical rules
tackle this problem from experimental data, or mathematically derive physical laws from first principles. However, the complexities of systems of 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.
centuries-old paradigm
One line of research stems from the dynamical systems literature, where the problem of interest was determine the unknown coefficients of a given differential equation from data. While the task seems rather simple, it was found that naive optimization cannot properly recover the desired coefficients
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

due to irregularities in the optimization landscape, which worsens with increasing nonlinearity and
data length 664 50 56
A parallel development occurred in the study of re- current neural networks (RNNs), which neural network analogues of discrete maps. It was quickly found that training these models can be quite unsta- ble and results in exploding or vanishing gradients. While earlier works in the field resulted in techniques such as teacher forcing or gradient clipping that are agnostic the RNN architecture [14 [5), more re- cent approaches have been directed towards circum- venting the problem entirely by incorporating specific mathematical forms into the model architecture to constrain the eigenvalues of the model Jacobian and 42). While
stabilize RNN dynamics 44 29
these latter methods can effective, for the purpose of discovering dynamical equations underlying the data, these approaches are inadequate as one cannot ascertain the hidden equations conform to the prede- termined rigid mathematical structures especially so if the dynamics to be predicted is incompatible with the constrained model architecture [57 41
Figure 1: Optimization trajectories for our homotopy method (left) and vanilla gradient descent (right). While convention training meanders on the pathological loss landscape, our method provides a series of relaxed land- scapes that effectively guide the optimizer to the loss minimum.
In this work, we focus on Neural Ordinary Differen- tial Equations (NeuralODEs) [[]] These models are powerful tools in modeling natural phenomena, bridging the expressibility and flexibility of neural networks with the de facto mathematical lan- guage the physical sciences. This has led to researchers amalgamating NeuralODEs with partial
information about the governing equation to produce "grey-box" dynamics models [55 71]. and
endowing NeuralODEs with mathematical structures that the system must satisfy [22] 18 Despite
their conceptual clegance, training NeuralODEs tend to result in long training times and sub-optimal
results, a problem that is further exacerbated as the length the training data grows (20 I7 Similar
to the case of RNNs, methods have been proposed to tackle the problem involve placing either
strong 12, 25 or semi-strong constraints [7](3)) to the functional form the NeuralODE can take something the underlying governing equation does not guarantee satisfying.
Contributions We introduce method to accurately train NeuralODEs long time series data by leveraging tools from the chaos and mathematical optimization literature: synchronization and homotopy optimization
| 15 Through loss landscape analyses, show that longer training data lead to more irregular loss landscapes, which in turn lead to deteriorating performances of the trained models. We show that synchronizing NeuralODE dynamics with the data can smooth the loss landscape, on which homotopy optimization can be applied to enhance training. Through performance benchmarks, we show that not only does our model far outperform conventional gradient descent training, it also surpasses multiple shooting in terms of convergence speed and extrapolation accuracy.
2 Preliminaries
2.1 Neural ordinary differential equations
NeuralODE 11 is model of the form,
du
u(t,u:)) u(to) Ug.
(I)
where Un RD is the initial condition or input given to the model, and U...: 0) R RD is
a neural network with parameters Tm that governs the dynamics of the model state u(l) E ".. While this model can be used as input-output mapping we are interested in the problem of training
In this case, NeuralODEs become the continuous analog of residual networks

---

NeuralODEs on time series data to learn the underlying governing equation and forecast future dynamics.
Given an monotonically increasing sequence of time points and the corresponding measurements {Elo) ( I Ji NeuralODE training starts with using an ordinary differential equation (ODE) solver
to numerically integrate Equation
(I) to obtain the model state given time points:
(0) 11 (0)
Ut, u: 0)tt ODESolve (U(+, u: 0); (()), up
(2)
g(0)
Afterwards, gradient descent is performed the loss function [(0) -N+1 Ei n1 (0) arrive at the optimal parameters 0"
(0)
(())
2.2 Homotopy optimization method
In topology, homotopy refers to the continuous deformation of one function into another. Motivated by this concept, to find the minimizer 0" I of complicated function F(0), the homotopy
optimization method [67 68 3 [5] introduces an alternative objective function
{
((0), if F(O). if
A=1
[( A)
x=
where G()) an auxillary function whose minimum is easily found, and C( A) R R R smoothly interpolates between G and J as the homotopy parameter À varies from 1 to O. While commonly used construction is the convex homotopy ((O, A) problem-specific forms can be employed.
()) AG more
Homotopy optimization shines in plicated landscape riddled with local minima. Similarly to simulated annealing [33], one starts out with relaxed version of the more complicated problem at hand and finds series of approximate solutions while slowly morphing the relaxed problem back into its original non-trivial form. This allows the optimization process not get stuck in spurious sharp valleys and accurately converge the minimum of interest.
where the target function F Fs non-convex and has com-
2.3 Synchronization dynamical systems
Here, we briefly summarize the mathematical phenomenon of synchronization, described in the 49 later sections, we will show that synchronization can used mitigate the difficulties arising in NeuralODE training. Consider two states ú, E RR, each evolving via du
literature [5 28
dú dt
(1, ú; 0)
U(, u: 0) - K K(u ú).
dt
Here, U U. + R" are the dynamics of the two systems parameterized by 0, 0, respectively, ), is the diagonal coupling matrix between the two systems. The two systems are said to synchronize if the error between û(t) and u(t) vanishes with increasing time.
and K diag(ki, kn,
Theorem 1 (Synchronization). Assuming U- = Ú O = ô and u(t) in the vicinity of int), the elements K can be chosen that the two states synchronize: that is, for error dynamics {(t) u(t) - û(t), [|E())]l2 = ||u(t) - ú(t)]l2 + as t + ac.
Proof. By subtracting the two state equations, the error, up first order terms, can be shown to follow:
dÉ(t)
du(t) dû(t) Tau(t, O)
E(t) = AE(().
(5)
dt dt dt
du
u=
Increasing the elements of K causes the matrix A to become negative definite which in turn, causes the largest Lyapunov exponent of the system,
|KE()))
lim t
log
(6)
[|E(O)]
It is also possible use other types losses, such as the L. loss used in 18 In the literature, these methods are also referred to as "continuation methods"

---

to also become negative as well. Then ||((t)]22 | - ú(t) - 0 as + oo, which means the dynamics of u and synchronize with increasing time regardless of the parameter mismatch.
Remark. Our assumption of = U. that the functional form of the data dynamics is the same as the equation be fitted can also be thought hold for NeuralODE training. This is because if the model is sufficiently expressive, universal approximation theorem will allow the model to match the data equation for some choice of model parameter.
2.3.1 Constructing homotopy with synchronization
To combine homotopy optimization with synchronization, slightly modify the coupling term of Equation (4) by multiplying it with the homotopy parameter A:
du dt U(,, u: XK ú).
With this modification, homotopy optimization is introduced to the problem as follows. When A and the coupling matrix K has sufficiently large elements, synchronization occurs and the NeualODE the data trajectory û(t) for large t. As will confirm in Section 4 this causes the resulting loss function C(O, 1) to become well-behaved, serving the role of the auxillary function G in Equation (3). When À = O, the coupled Equation (4) reduces to the
prediction u(t) starts to converge
original model form of Equation (.. Accordingly, the corresponding loss function L()) = L(0, 0) is the complicated loss function
need ultimately minimize. Therefore, starting with À 1 and successively decreasing its value to O, all the while optimizing for the coupled loss function L(O, A) allows one leverage the well-behaved loss function landscape from synchronization while being
able to properly uncover the system parameters 66 58
3 Related works
Synchronization in ODE parameter estimation After the discovery of synchronization in chaotic
systems by Pecora & Carroll 47 48), methods to leverage this phenomenon to identify unknown ODE parameters were first proposed by Parlitz 44. and later refined by Abarbanel et al. [2 [ Parameter identification in these earlier works would proceed by either designing an auxillary dynamics for the ODE parameters [4 [4), including regularization term on the coupling strength in the loss function [53] [.. The idea of applying homotopy optimization to this problem was pioneered by
Vyasarayani al. 655 66 and further expanded upon by later works [58 63 52
Homotopy methods in machine learning As homotopy methods can outperform gradient descent for complex optimization problems, it has been used in diverse fields, including Gaussian homotopy
for non-convex optimization (27) and training feedforward networks 13 24] 10]. Curriculum learning [6] is another example of homotopy optimization, with the each task in the curriculum corresponding to particular value of the homotopy parameter as it is decreased to O. Not to be confused with the method discussed in this paper is the homotopy perturbation technique (26) used in the dynamical systems literature, which is a method to obtain analytical solutions to nonlinear differential equations.
Improving the training of RNNs and NeuralODEs As previously mentioned, many recent approaches to improve RNN and NeuralODE training involve constraining the model expressivity. The limitations of this strategy is investigated in [57], where it is demonstrated that such "stabilized" models cannot learn oscillatory or highly complex data that does match the built-in restrictions.
recent development with similar spirit ours is Mikhaeil et al. [41], where the authors developed an architecture- agnosst training method for RNNs by revisiting the classical RNN training strategy of teacher forcing [60] and refining it via the chaos theory concept of Lyapunov exponents. Other aspects of NeuralODE training, such as better gradient calculation, have been studied to improve NeuralODE performance 30 72
39 31 Such methods are fully compatible with our proposed work, and we showease an example in Appendix G by changing the gradient caleulation scheme to the "symplectic-adjoii waitettiiicciiiii-eemtttooo from 39].

---

4 Loss landscapes in learning ODEs
4.1 Complicated loss landscapes in ODE parameter estimation
NeuralODEs suffer from the same training difficulties of ODE fitting and RNNs. Here we illustrate that the underlying cause for this difficulty is also identical by analyzing the loss landscape of NeuralODEs for increasing lengths of training data. The data were generated from the periodic Lotka- Volterra system and the landscapes were visualized by calculating the two dominant eigenvectors of
the loss Hessian, then perturbing the model parameter along those directions (9). We also provide
analogous analyses for ODE parameter estimation on the same system in Appendix B to further emphasize the similarities between the problem domains.
COAN
CIRT
CI
Quuling cherge
Figure 2: Irregular loss landscape in NeuralODE training. (Upper left) Loss landscape and (Lower left) eigenvalue spectrum of the loss Hessian for increasing lengths of train data. (Upper right) Loss landscape and (Lower right) eigenvalue spectrum of the loss Hessian with increasing coupling strength. For clarity, loss values were clipped above at 4000. All Hessian related information was calculated using the PyHessian package (69).
From Figure 2 (upper left), we find that the loss landscape is well-behaved for short data, but becomes pathological as the data lengthens. In fact, the landscape for t [0, 9.1] features stiff cliff, which is known to cause exploding gradients in RNN training 14 45) I Indeed, this picture is also confirmed by alternative principal component analysis- based visualization of the loss landscape and the optimization path [35 (Figure ! right). We also observe that the overall magnitude of the loss increases with data length. This is direct consequence of model dynamics utt) being independent of the data dynamics (t) and therefore tending to diverge away with time.
To further confirm our findings and compensate for the fact that our visualizations are low dimensional projections of the actual loss landscape, we computed the eigenvalue spectrum of the loss function hessian for increasing lengths of training data (Figure 2 lower left). The computed spectrums show that hessian eigenvalues spread away from zero with increasing length of training data, signalling that the loss morphs from a relatively flat to more corrugated landscape. We also observe large positive eigenvalue outlier peaks, whose values also rapidly increase with the length of training data. As such large outliers have been report to slow down training [19), this also provides additional insights as to why NeuralODE training stagnates on long time series data.
4.2
Loss landscape of synchronized dynamics
With coupling, the loss function L between the model predictions and the data now becomes a function of both the model parameters @ and the coupling strength k. Increasing the coupling strength k effectively acts as shortening the data because increased coupling causes the two trajectories to start synchronizing earlier in time, after which the error between the two will diminish and contribute negligibly to the total loss. Therefore, synchronizing the model and the data generating equations leads to a more favorable loss landscape that gradient descent based optimizers will not struggle on.
4 'For more details, see Section 6
For details on the visualization method, see Section B of the supplementary of the cited reference.
5

---

Once again, we confirm this point by visualizing the loss landscape and the Hessian eigenvalue spectrum for long time series dataset while increasing the strength of the coupling. From Figure 2 (upper right) we find that with increasing coupling strength, the steep cliff and the flat basin of the original loss surface are tamed into a well-behaved slope. This is reflected in the eigenvalue spectrum as well (Figure 2 lower right), with the all of the spectrum gathering near zero as the coupling is turned up. In fact, when the coupling strength is increased too much, we find that all the eigenvalues bunch together at zero, meaning that the loss function becomes completely flat. Qualitatively, this corresponds to the case where the coupling completely dominates the synchronized dynamics and forces the model predictions onto the data regardless of the model parameters.
Homotopy optimization for NeuralODE training
Having established that applying synchronization between the model and data can tame the irregular loss landscape during training. we now turn towards adapting the amnchronzzaion-homooyg approach of Section 2.3. to NeuralODE setting. In doing so, there couple details that needs addressing arrive at specific implementations. We discuss such points below.
5.1 Constructing the coupling term
The coupling term for synchronization - AK((u- û) in Equation depends on the data dynamics ú, whose values are only available at discrete timestamps {t() Ji=o This is at a disaccord with the ODE solvers used during training, which require evaluating the model at intermediate time points as well. Therefore, to construct a coupling term defined for arbitrary time, we used a cubic smoothing spline (21) of the data points to supply the ú values unobserved times. Smoothing was chosen over interpolation to increase the robustness of our approach to measurement noise.
Note that cubic spline is not the only solution even defining the coupling term as a sequence of impulses is possible, only applying where it is defined, and keeping it zero at other times. This is an interesting avenue of future research, as this scheme can be shown 53 to be mathematically equivalent to teacher forcing [60].
5.2 A little more on the coupling matrix
In general, the coupling matrix K can be any n x matrix, provided it renders the Jacobian of the eтror dynamics (matrix A in Equation (6)) to be negative definite along the data trajectory ú(t). However, naively following this approach leads to n² hyperparameters just for the coupling matrix, leads which can quickly get out of hand for high dimensional systems. To circumvent this issue, in this work have used a coupling matrix of the form K kI, where I, the n x n identity matrix, to arrive at single scalar hyperparameter characterizing coupling strength.
5.3 Scheduling the homotopy parameter
The homotopy parameter can be decreased in various different ways, the simplest approach being Water K²/ "" nd where AXO) denotes the j-th decrement and Sstep is number of homotopy steps. This was inspired from our observations of the synchronized loss landscape (Figure right), where the overall shape change as well as the width of the eigenvalue spectrum seemed be a superlinear function of the coupling strength.
constant steps. In our study, we instead used power law decrements of the form (A\ 1==
5.4 Hyperparameter overview
Our final implementation of the homotopy training algorithm has five hyperparameters, as described below.
Coupling strength (k) This determines how trivial the auxillary function for the homotopy opti- mization will be. Too small, and even the auxillary function will have jagged landscape; too large, and the initial auxillary function will become flat (Figures and TT left panel, = 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.

---

Algorithm 1 Homotopy-based NeuralODE training
Input: Data {((0) [
12=0+ untrained model 0(...;0)
Hyperparameters: Coupling strength(k), homotopy parameter decrement ratio(s), number of homotopy steps(s), epochs per homotopy step(Mepoch). learning rate(n)
i(t)= CwwiccSmothinngpllie.. ((t(1) (i0 (o)))
> Construct data interpolant
x-1
Initialize homotopy parameter
for in 0 ..t.estpp
do
for epoch in 0.. ...oooooch Nepoch do
[u(O]N i=0 = ODESolve (ú U(t. u: 0) - AAI - a(t))) > Solve coupled dynamics
[ū()]N 1=0 = ODE (n U(t. 0))
D Also solve uncoupled dynamics
MSE N+1
E,
(6) (0) ú())
if MSE MS Ebest then
MSet + MSE
>> Model performance is gauged by uncoupled dynamics error > Checkpoint best model parameters
0" t
end if
(()) N+I Ei ( 5() θ + O-rttmizeerteppp
Loss function is calculated using coupled dynamics
VeL)
end for
f- A- K²/ / ===0
Metepa-1
Decrease homotopy parameter with power law decrements
end for Output: Trained model parameters A 0
Homotopy parameter decrement ratio (x) This determines how the homotopy parameter is decremented for each step. Values close to cause À to decrease in nearly equal decrements, whereas smaller values cause large decrease of À in the earlier parts of the training, followed by subtler decrements later on. We empirically find that x values near 0.6 tend to work well.
Number of homotopy steps (s) This determines how many relaxed problem the optimization 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 step (Nepoch) This determines how long the model will train given homotopy parameter value A. Too small, and the model lacks the time to properly converge on the loss function; too large, and the model overfits on the simpler loss landscape of X & O, resulting in reduced final performance when À = O. We find for simpler monotonic or periodic systems, values of 100-150 work well; for more irregular systems, 200-300 are suitable.
Learning rate (n) This is as same as in conventional NeuralODE training. We found values in the range of 0.002-0.1 to be adequate for our experiments.
6 Experiments
To benchmark our homotopy training algorithm against existing methods, we trained multiple NeuralODE models across time-series datasets of differing complexities. The performance metrics we used in our experiments were mean squared етrоr on the training dataset and on model rollouts of 100% the training data length, which from now on we refer as interpolation and extrapolation MSE, respectively. For each training method, the hyperparameters were determined through grid search performed prior to the main experiments. All experiments were repeated three times with different random seeds, and we report the corresponding means and standard errors whenever applicable. For conciseness, we outline the datasets, models and the baseline methods used in the subsequent sections, and refer interested readers to Appendices F and F for extensive details on data generation, model configuration and hyperparameter selection. We provide all of the code for this paper in https:/ / / /gi thub. com http://hhho//////eeaaaaaaataaaaa New.raaaeEeeEeootooo...
7

---

6.1 Datasets and model architectures
Lotka- Voltera system The Lotka-Voltena system is a simplified model of prey and predator populations [3. The system dynamics are periodic, meaning even models trained on short trajectories
should have strong extrapolation properties, provided the training data contains a full period of
oscillation. We considered two types of models for this system: "black-box" model of Equation (
and "gray-box" model used in (54) given by dz ar + Ui (2 yi 01), dy -Y + (2, ys 02).
Double pendulum For this system, we used the real-world measurement data from Schmidt & Lipson [59]. While the double pendulum system can exhibit chaotic behavior, we used the portion of the data corresponding to the aperiodic regime. We trained two types of models on this data: a black-box model and model with second-order structure. Our motivation for including the latter model stems from the ubiquity of second-order systems in physics and from 23 which reports these models can match the performance of Hamiltonian neural networks.
Lorenz system The Lorenz system displays highly chaotic behavior, and serves stress test for time series prediction tasks. For this system, we only employed a black-box model.
6.2 Baselines
As the focus of this paper is to effectively train NeuralODEs without employing any specialized model design or constraining model expressivity, this was also the criteria by which we chose the baseline training method for our experiments. Two training methods were chosen as our baseline: vanilla gradient descent, which refers to the standard method of training NeuralODEs, and multiple-shooting method. As the latter also stems from the ODE parameter estimation literature, and its inner working relevant to the discussion of our results, briefly describe the multiple-shootin algorithm before discussing our results.
Multiple shooting The muliple-shooting algorithm 62 g circumvents the difficulty of fitting models on long sequences by dividing the training data into shorter segments and solving multiple initial value problems. To keep the predictions for different segments in syne, the loss function is augmented with a continuity enforcing term. Unlike synchronization and homotopy, which we employ in NeuralODE training for the first time, works applying multiple shooting to NeuralODEs are scant but do exist [54] 38 61 As the three previous references differ in their implementations, we chose the scheme [54] for its simplicity and implemented a PyTorch adaptation to use for our experiments. Further details about our multiple-shotingg implementation can be found in Appendix D
7 Results
7.1
Training NeuralODEs time series data of increasing lengths
In the previous sections, we found that longer training data leads to a more irregular loss landscape. To confirm that this does translate to deteriorated optimization, we trained NeuralODEs on Lotka- Volterra system trajectories of different lengths. Inspecting the results (Figure 3 first panel) we find while all training methods show low interpolation MSE for shortest data, this due to overfitting, as evidenced by the diverging extrapolation MSE As the length of the data grows and encompasses a full oscillation period, we find that the model performances from homotopy and multiple-shootine stabilize for both interpolation and extrapolation. This is clearly not the case for the vanilla gradient descent, whose model performances degrade exponentially with as data length increases.
Model capacity is not the problem The reason NeuralODEs struggle on long sequences can easily be attributed to insufficient model capacity. We demonstrate this is not the case by repeating the previous experiment for the mid-length data while decreasing the model size used (Figure 3 second panel). In the case of vanilla gradient descent, the model size is directly correlated to its performance,
WWiile not explicitly mentioned the paper, DiffEqFlux. j1 library has an example tutorial using this technique to train NeuralODEs. Indeed, plots of the model prediction trajectory in Appendix G.3 also confirm this fact.
8

---

Interpstetion MISE
100
10F
120
10-1
2014
2012
Extrapotetton MSEE
Homotopy Rhattahoet Vaniila
10
10
704
100
20
10-1
20Th
1000
2012
length
Nodes
Sampiig period
Noise strengtth (%
Figure 3: Results for the Voltera system. (First) Interpolation and extrapolation errors for increasing data length. (Second) Errors for the mid-length data with decreasing model capacity. (Third) Errors for fixed time interval but with increasing data sparsity. (Fourth) Errors for mid-length data with increasing noise.
lending the claim some support. However, we find that both homotopy and multiple shooting methods can effectively train much smaller models to match the accuracy of a larger model trained with vanilla gradient descent. This signifies that the cause of the problem not the capacity itself, but rather the inability of training algorithms to fully capitalize on the model capacity.
Robustness against data degradation We performed additional set of experiments evaluate the performances of different training methods against increasing data sparsity and noise. This is particularly of interest for our homotopy training as we employ smoothing cubic splines (Section 5. whose reliability can be adversely affected by the data quality. From the results (Figure 3 third, fourth panels) we find this not the case and that that our algorithm performs well even for degraded data. One unexpected result is the drastic underperformance of the multiple shooting method with increased sampling period, which is side effect from the fact we held the time interval fixed for this experiment. Due to the space constraint, we continue this discussion in Appendix G.3
7.2 Performance benchmarks the training algorithms
Having confirmed that both the homotopy and multiple-shootin methods successfully mitigate the NeuralODE training problem, we now present comprehensive benchmark in Figure
Hunology
4330
Multahost
Vani
20
1030
530
Abrk D0yy die Socved Diiie Morck
dark
Litba Vollora
Deaale Ariidar
Lobto Walloa
Lotlel Wallere
Dabie Feasiker
Figure 4: Performance benchmarks. (Left) Interpolation and (Center) extrapolation MSEs for each training method. (Right) Epochs where the minimum interpolation MSE was achieved. For all plots, the bar and the error bars denote the mean and standard errors across three runs.
From the results, it is clear both our homotopy method and multiple shooting can both boost the interpolation performance of NeuralODEs, with our method being slightly more performant. Inspecting the extrapolation performance, however, we do find a clear performance gap between our method and multiple shooting, with the latter even performing worse than vanilla gradient descent for the double pendulum dataset (Figure 4 center, Figure [] left). This underperformance of multiple shooting for extrapolation tasks can attributed to the nature of algorithm: during training time, multiple shooting integrates NeuralODEs on shorter, subdivided time intervals, whereas during inference time the NeuralODB needs be integrated for the entire data time points and then some more for extrapolation. This disparity between training and inference render models trained by

---

multiple shooting to be less likely to compensate for error accumulation in long term predictions, resulting in poor extrapolation performance.
Another stark difference between our homotopy training and the other methods is in the number of epochs required to arrive minimum interpolation MSE (Figure + right panel). While the exact amount of speedup is variable and seems to depend on the dataset difficulty, we find a 90%- 40% reduction compared vanilla training, highlighting the effectiveness of our approach. One unexpected observation from our experiments was that the epoch at which minimum interpolation MSE epoch was achieved did not always correspond the homotopy parameter À being zero (Appendix G.2). While seemingly counterintuitive, this is in line with the literature, where nonzero homotopy parameter has resulted in better generalization performances [36 24).
What makes homotopy training effective: inspecting training dynamics To obtain a comple- mentary insight to the NeuralODE training dynamics and relate the effectiveness of our method with the properties of the loss landscape, we calculated the trace of the loss hessian during training (Figure right panel). Small, positive trace values indicate more favorable training dynamics [69], and we find that our method results in the smallest trace, followed by multiple shooting then vanilla gradient descent. The fact that the hessian trace maintains low value throughout the entire homotopy training, even as the loss surface is morphed back to its unsynchronized, irregular state, also confirms that our homotopy schedule works well and managed to gradually guide the optimizer into a well performing loss minimum.
Stat 9
Hezsion trace
0.0
0.5
1.0
1.5
2.0
0.0
0.5
1.0
1.5
2.0
Stato us.
101
Stere Mt
20
Porricooy
101
-hhost
-10
-20
0.0
0.5
1.0
1.5
2.0
0.0
0.5
1.0
1.5
2.0
1000
2000
3000
4000
Time
Time C
Train epochs
Figure 5: Second-order NeuralODE training results for the double pendulum dataset. (Left) Predicted trajectories. Lines and bands correspond to mean and standard errors for three training runs and the dashed line indicate the extrapolation start time, (Right) Estimation of the hessian trace during training. The homotopy curve stops early due the difference in the number of maximum training curves used.
8 Limitations and discussion
In this paper, we adapted the concepts of synchronization and homotopy optimization for the first time in the NeuralODE literature stabilize NeuralODE training on long time-series data. Through loss landscape analysis. we established that the cause of the training problem is identical between differential equations. RNNs and NeuralODEs and that synchronization can effectively tame the irregular landscape. We demonstrated that our homotopy method can successfully train NeuralODF models with strong interpolation and extrapolation accuracies.
Note that our work serves as an initial confirmation that the methods developed in the ODE parameter estimation literature can be adapted be competitive in the NeuralODE setting. As such, there are multiple opportunities for future research, including theoretical studies on the convergence properties of our algorithm and further testing out ideas in the synchronization literature, such alternative parametrization for the coupling matrix 228, or incorporating the homotopy parameter into the loss function regularization term [ 588. Extending our algorithm to high-dimensional systems is another important direction, which we believe will be possible through the use of latent ODE models 70. Nevertheless, our work hold merit in that lays the groundwork upon which ideas from ODE parameter estimation, control theory, and RNN/NeuralODE training can be assimilated, which we believe will facilliate further research connecting these domains.
Author Contributions
J. Ko and H. Koh contributed equally to the paper.
10

---

Acknowledgments
This work was supported by grants from the National Research Foundation of Korea (No. 2016R1A3B190866) to W. Jhe.
References
[1] H. D. I. Abarbanel, D. Creveling. R. Farsian, and M. Kostuk. Dynamical State and Parameter Estimation. SIAM Journal on Applied Dynamical Systems, 8(4):1341-131,, 2009.
[2] H. D. L. Abarbanel, D. R. Creveling, and J. M. Jeanne. Estimation of parameters in nonlinear systems using balanced synchronization. Physical Review E, 77(1):016208, 2008.
[3] E. L.. Allgower and K. Georg. Numerical Continuation Methods: An Introduction. Springer Berlin Heidelberg, Berlin, Heidelberg, 1990.
[4] M. Arjovsky, A. Shah, and Y. Bengio. Unitary Evolution Recurrent Neural Networks. In Proceedings of The 33rd International Conference Machine Learning, pages 1120-1128.
PMLR, 2016.
[5] E. Baake, M. Baake, H. G. Bock, and K. M. Briggs. Fitting ordinary differential equations to chaotic data. Physical Review A, 45(3):5524-5529 1992.
[6] Y. Bengio,, J. Louradour, R. Collobert, and J. Weston. Curriculum learning. In Proceedings of the 26th Annual International Conference on Machine Learning, ICML '09, pages 41 41-48, New York, NY, USA, 2009. Association for Computing Machinery.
[7] L. Biewald. Experiment tracking with weights and biases, 2020. Software available from wandb.com.
[8] H. G. Bock. Recent Advances in Parameter Identification Techniques for O.D.E. In P. Deuflhard and E. Hairer, editors, Numerical Treatment of Inverse Problems in Differential and Integral Equations:
Proceedings of an International Workshop. Heidelberg, Fed. Rep. of Germany, August 30- September 3, 1982, Progress in Scientific Computing, pages 95 121, Birkhäuser, Boston, MA, 1983.
[9] B. Chang, M. Chen, E. Haber, and E. H. Chi. AntisymmetricRNN: A Dynamical System View on Recurrent Neural Networks. International Conference on Learning Representations, volume 7, 2019.
[10] Q. Chen and W. Hao. A homotopy training algorithm for fully connected neural networks. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 475(2231):201994) 200622 Nov. 2019.
[11] R. Q. Chen, Y. Rubanova, Bettencourt, and D. K. Duvenaud. Neural Ordinary Differen- tial Equations. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
[12] K. M. Choromanski, J. Q. Davis, V. Likhosherstov, X. Song, J.-J. Slotine, J. Varley, H. Lee, A. Weller, and V. Sindhwani. Ode to an ODE. Advances in Neural Information Processing Systems, volume 33, pages 3338-3350. Curran Associates, Inc., 2020.
[13] Chow, L. Udpa, and Udpa. Homotopy continuation methods for neural networks. In 1991., IEEE International Sympoisum on Circuits and Systems, volume 5, pages 2483-2486, 1991.
[14] K. Doya. Bifurcations of recurrent neural networks in gradient descent learning. IEEE Transactions on neural networks, 1(75):218-229, 1993.
[15] D. M. Dunlavy and D. P. O'LLery Homotopy optimization methods for global optimization. Technical Report SAND2005-7495, Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA (United States), 2005.
[16] W. Falcon and The PyTorch Lightning team. PyTorch Lightning, Mar. 2019.

---

[17] C. F.Filla, J.-H. Jacobsen, L. Nurbekyan, and A. Oberman. How to Train Your Neural ODE: The World of Jacobian and Kinetic Regularization. In Proceedings of the 37th International Conference Machine Learning. volume 37, pages 3154-3164. PMLR, 2020.
[18] M. Finzi, K. A. Wang, and A. G. Wilson. Simplifying Hamiltonian and Lagrangian Neural Networks via Explicit Constraints. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 13880-13889, Curran Associates, Inc., 2020.
[19] B. Ghorbani, S. Krishnan, and Y. Xiao. An Investigation into Neural Net Optimization via Hessian Eigenvalue Density. In Proceedings of the 36th International Conference on Machine Learning, pages 2232-222 2241. PMLR, 2019.
[20] A. Ghosh, H. Behl, E. Dupont, P. Torr, and V. Namboodiri. STEER SSimple Temporal Regularization For Neural ODE, In Advances in Neural Information Processing Systems, volume 33, pages 14831-14843, Curran Associates, Inc., 2020.
[21] P. J. Green, B. W. Silverman, P. J. Green, and B. W. Silverman. Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach. Number 58 Monographs on Statistics and Applied Probability. Chapman & Hall /CRC, Boca Raton, Fla., 2000.
[22] S. Greydanus, M. Dzamba, and J. Yosinski. Hamiltonian Neural Networks. In H. Wallach, H. Larochelle, A. Beygelzimer, A'Alché-Bue, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
[23] N. Gruver, M. A. Finzi, S. D. Stanton, and A. G. Wilson. Deconstructing the Inductive Biases of Hamiltonian Neural Networks. In International Conference on Learning Representations. volume 9, 2021.
[24] C. Gulcehre, M. Moczulski, F. Visin, and Y. Bengio. Mollifying Networks. In International Conference on Learning Representations, volume 5, 2017.
[25] R. Hasani, M. Lechner, A. Amini, D. Rus, and R. Grosu. Liquid Time-constant Networks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 7657-7666, 2021.
[26] J.-H. He. Homotopy perturbation technique. Computer Methods in Applied Mechanics and Engineering, 172335)):57-222, 1999.
[27] H. Iwakiri, Y. Wang, S. Ito, and A. Takeda. Single Loop Gaussian Homotopy Method for Non-convex Optimization. Advances in Neural Information Processing Systems, 35:7065-7076. 2022.
[28] G. A. Johnson, D. J. Mar, T. L. Carroll, and L. M. Pecora. Synchronization and imposed bifurcations in the presence of large parameter mismatch. Physical Review Letters, 80(18 3):3966 3959, 1998.
[29] S. Kanai, Y. Fujiwara, and S. Iwamura. Preventing Gradient Explosions in Gated Recurrent Units. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
[30] P. Kidger, R. T. Q. Chen, and T. J. Lyons. "Hey, that's not an ODE": Faster ODE Adjoints via Seminorms. In Proceedings of the 38th International Conference Machine Learning, pages 5443-5452. PMLR, 2021.
[31] S. Kim, W. Ji, S. Deng, Y. Ma, and C. Rackauckas. Stiff neural ordinary differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science. 31(9):093122, 2021.
[32] A. Krishnapriyan, A. Gholami, S. Zhe, R. Kirby, and M. W. Mahoney. Characterizing possi- ble failure modes in physics-ifformed neural networks. In Advances in Neural Information Processing Systems, volume 34, pages 26548 -665000 Curran Associates, Inc., 2021.
[33] J. M. Laarhoven and F. H. Aarts. Simulated Annealing: Theory and Applications. Springer Netherlands, Dordrecht, 1987.
12

---

[34] M. Lezcano-Casaal and D. Martínez-Rubee Cheap Orthogonal Constraints in Neural Networks: Simple Parametrization of the Orthogonal and Unitary Group. In Proceedings of the 36th International Conference on Machine Learning, pages 7799-3803. PMLR, 2019.
[35] H. Li, Z. Xu, G. Taylor, C. Studer, and Goldstein. Visualizing the Loss Landscape of Neural Nets. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
[36] X. Lin, Z. Yang, X. Zhang, and Q. Zhang. Continuation path learning for homotopy optimization. In Proceedings of the 40th International Conference on Machine Learning, pages 21288 21311. PMLR, 2023.
[37] I. Loshchilov and F. Hutter. Decoupled Weight Decay Regularization. In International Confer- ence Learning Representations, volume 2019.
[38] S. Massaroli, M. Poli, S. Sonoda, T. Suzuki, J. Park, A. Yamashita, and H. Asama. Differentiable Multiple Shooting Layers. In Advances in Neural Information Processing Systems, volume 34, pages -6542-16544. Curran Associates, Inc., 2021.
[39] T. Matsubara, Y. Miyatake, and T. Yaguchi. Symplectic Adjoint Method for Exact Gradient of Neural ODE with Minimal Memory. In Advances in Neural Information Processing Systems, volume 34, pages 20772-20784. Curran Associates, Inc., 2021.
[40] A. Maybhate and R. E. Amritkar. Use of synchronization and adaptive control in parameter estimation from time series. Physical Review E, 59(11)284429)) 293, 1999.
[41] J. Mikhaeil, Z. Monfared, and D. Durstewitz. On the difficulty of learning chaotic dynamics with RNNs. Advances in Neural Information Processing Systems, 35:1 1192-1 2022.
[42] J. Miller and M. Hardt. Stable Recurrent Models. In International Conference on Learning Representations, volume 7, 2019.
[43] J. D. Murray. Mathematical Biology. 1: An Introduction. Number 17 in Interdisciplinary Applied Mathematics. Springer- GmbH, Berlin Heidelberg, softcover reprint of the hardcover 3rd edition 2002, corrected second printing edition, 2004.
[44] U. Parlitz. Estimating Model Parameters from Time Series by Autosynchronization. Physical Review Letters, (118):1232-133: 1996.
[45] R. Pascanu, T. Mikolov, and Y. Bengio. On the difficulty of training recurrent neural networks. In Proceedings of the 30th International Conference Machine Learning, pages 1310-1318. PMLR, 2013.
[46] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang. Z. VeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L.. Fang, J. Bai, and S. Chintala. Pytorch: An imperative style, high-performance deep learning library. arXiv preprint 1rXiv:19/2.01703. 2019.
[47] L. M. Pecora and T. L. Carroll. Synchronization in chaotic systems. Physical Review Letters, 64(8):821-824, 1990.
[48] L. M. Pecora and T. L. Carroll. Driving systems with chaotic signals. Physical Review A. 44(4):2374-2383, 1991.
[49] L. M. Pecora and T. L. Carroll. Synchronization of chaotic systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(9):097611, 2015.
[50] M. Peifer and J. Timmer. Parameter estimation in ordinary differential equations for biochemical processes using the method of multiple shooting. IET Systems Biology, 1(2):78-88, 2007.
[51] H. Peng, E. J. Ding, M. Ding, and W. Yang. Synchronizing hyperchaos with a scalar transmitted signal. Physical Review Letters, 7606))904-907, 1996.
13

---

[52] V. G. Polisetty, S. K. Varanasi, and P. Jampana. Stochastic state-feedback control using homotopy optimization and particle filtering, International Journal of Dynamics and Control, 10(3):94-9555 2022.
[53] J. C. Quinn, P. H. Bryant, D. R. Creveling, S. R. Klein, and H. D. I. Abarbanel. Parameter and state estimation of experimental chaotic systems using synchronization. Physical Review E. 80(1):016201, 2009.
[54] C. Rackauckas, M. Innes, Y. Ma, J. Bettencourt, L. White, and V. Dixit. DiffEqFlux.ji A Julia Library for Neural Differential Equations, 2019.
[55] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, Skinner, A. Ramadhan, and A. Edelman. Universal Differential Equations for Scientific Machine Learning, 2021.
[56] J. O. Ramsay, G. Hooker, D. Campbell, and J. Cao. Parameter estimation for differential equations: A generalized smoothing approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(5):741-796, 2007.
[57] A. H. Ribeiro, K. Tiels, L. A. Aguirre, and T. Schön. Beyond exploding and vanishing gradients: Analysing RNN training using attractors and smoothness. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Stutistics, pages 2370-2380. PMLR, 2020.
[58] K. Schäfer, K. FlaBkamp, J. Fliege, and C. Büskens. A Combined CoomoooyyOpttimizann Approach to Parameter Identification for Dynamical Systems. PAMM, 19(1), 2019.
[59] M. Schmidt and H. Lipson. Distilling Free-Form Natural Laws from Experimental Data. Science, 324((923):81-35 2009.
[60] N. B. Toomarian and J. Barhen, Learning a trajectory using adjoint functions and teacher forcing, Neural Networks, 5(3):473-444, 1992.
[61] E. M. Turan and J. Jäschke. Multiple shooting for training neural differential equations on time series. IEEE Control Systems Letters, 6: 6:1897-1902, 2022.
[62] B. van Domselaar and P. Hemker. Nonlinear parameter estimation in initial value problems. Technical Report MC-NW-75-18, Amsterdam Univ. Math. Cent. Num. Sci., 1975.
[63] S. K. Varanasi, P. Jampana, and C. P. Vyasarayani. Minimum attention stochastic control with homotopy optimization. International Journal of Dynamics and Control, 9(1 266,246,,, 2021,
[64] H. U. Voss, J. Timmer, J.d Kurths. Nonlinear Dynamical System Identification from Uncertain and Indirect Measurements. International Journal of Bifurcation and Chaos, 19056))19 2004.
[65] C. P. Vyasarayani, T. Uchida, A. Carvalho, and J. McPhee. Parameter identification in dynamic systems using the homotopy optimization approach. Multibody System Dynamics, 26(4):411- 424, 2011.
[66] C. P. Vyasarayani, T. Uchida, and J. McPhee. Single-shooting homotopy method for parameter identification in dynamical systems. Physical Review E, (5(3))033001,211
[67] L. T. Watson. Globally convergent homotopy methods: A tutorial. Applied Muthemutics and Computation, 31:369-396, May 1989.
[68] L. T. Watson and R. Haftka. Modern homotopy methods in optimization. Computer Methods in Applied Mechanics and Engineering, 74(3):289-305, 1989.
[69] 7. Yao, A. Gholami, K. Keutzer, and M. W. Mahoney. PyHessian: Neural Networks Through the Lens of the Hessian. In 2020 IEEE International Conference on Big Data (Big Data), pages 581-590, IEEE, 2020.
[70] C. Yildiz, M. Heinonen, and H. Lahdesmaki. ODE2VAE: Deep generative second order ODEs with Bayesian neural networks. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
14

---

[71] Y. Yin, Le Guen, J. Dona, E. de Bézenac, I. Ayed, N. Thome, and P. Gallinari. Augmenting physical models with deep networks for complex dynamics forecasting. Journal of Statistical Mechanics: Theory and Experiment, 2021(12):124012, 2021.
[72] J. Zhuang, N. Dvomek, X. Li, S. Tatikonda, X. Papademetris, and J. Duncan. Adaptive Checkpoint Adjoint Method for Gradient Estimation in Neural ODE. In Proceedings of the 37th International Conference Machine Learning, pages 11639-11649 11649, PMLR, 2020.
15

---

A Additional figures
In this section, we present some of the figures that were excluded from the main text for brevity.
Figure O below highlights the efficiency of our training approach with respect to the number of training epochs by plotting the trajectory predictions of the second-order models during training.
Ho Homolooyy Ecoch 32
Epoch
Epeci
Epech - 9303bct
shont.
Spect
Spact 3285-baes
Figure 6: Second-order model predictions for the double pendulum data during training. Blue and orange lines indicate data and model predictions respectively. (Top) Results from our homotopy approach. (Middle) Results from multiple shooting. (Bottom) Results from vanilla gradient descent.
Figures to show model prediction trajectories corresponding One intriguing point is that even though both the gray-box model (Lotka- Volterra system, Figure and the second-order model (double pendulum, Figure contain more prior information about their respective systems, this does not necessarily translate to hetter predictive capabilities if an effective training method is used. This suggests that introducing physics into the learning problem can obfuscate optimization, something that has been reported for physics-informe neural networks [32]. also highlights the effectiveness of our homotopy training algorithm, our method can properly train such difficul-to-optimize models.
benchmark results of Figure 4
Black-box
Gray-box
10
10
Homotopy
Multisshoot
Vanill
Estrupi statt Cine
State
10
Timee
Time
Figure 7: Trajectory predictions for the Lotka- Lotka-Volterra system from models trained with different methods. Dashed vertical line indicates the start of extrapolation.
16

---

Black-box
Second-order
0:
Hometopy
Mullisholl
Varilla
1.0
0.0
State
0.0
0.0
-2.
-2.
1.0
2.0
5.0
State ut:
2.0
13
uz
20
20
1.0
2.0
0.0
100
1.3
Time
Time
Figure 8: Trajectory predictions for the double pendulum from models trained with different methods. Dashed vertical line indicates the start of extrapolation.
20
Homotopy
40
15
Multi-shoot
State
10
Vanilla
Extrap. start time
State y
10-
30
10-
-10
2
Figure 9: Trajectory predictions for the Lorenz system from models trained with different methods. Dashed vertical line indicates the start of extrapolation.
B Loss landscape and synchronization in ODE parameter estimation
Here, we illustrate the similarity between the ODE parameter estimation problem and the NeuralODE training problem by repeating analyses analogous to those of Section in the setting of estimating the unknown coefficients of given ordinary differential equation.
In the left and middle panels of Figure 10 we show the trajectories of the Lotka- Volteraa and the Lorenz systems where single parameter was perturbed. We can easily observe that the original trajectory and the perturbed trajectories evolve independently in time. Furthermore, while trajectories from the simpler Lotka-Volterra Volterra system retain the same periodic behavior and differ only in their amplitudes and phases, trajectories from the Lorenz system display much different characteristics with increasing time. This translates over to the shape of the loss landscape (Figure [0] right) with the loss for the Lotka- Volterra system only becoming steeper with increasing data length, whereas the loss for the Lorenz system displays more and more local minima.
Once synchronization is introduced to the perturbed trajectory, its dynamics tends to converge with the reference trajectory as time increases, with the rate of convergence increasing for larger values of the coupling strength & (Figure 11 left and middle panels). In terms of the loss landscape (right panel), larger coupling does result in smoother landscape, with excessive amounts of resulting in a flat landscape that is also detrimental to effective learning. These results are a direct parallel to our observations Figure and provide additional justifications as to why techniques developed in the field of ODE parameter estimation also work well on NeuralODE training.
17

---

10-
20
+50%
-50%
5
25
100
Lenu 50
Len.: 100
25
-40
-20
20
40
Parameter perturbation
Figure 10: Dynamics of single parameter perturbed systems. The systems and the parameters used are described Section 6 (Left) Solutions for the periodic Lotka- Volterra equations with perturbed parameter. (Middle) Solutions for the chaotic Lorenz system with perturbed 3 parameter. (Right) MSE loss function landscape for the Lotka- Volterra equations (right, upper) and Lorenz system (right, lower) for different lengths of time series data for the loss calculation.
20
k.0.0
10
= 0.0
0.1
2-0
0.1
k= 1.0
100
Vmn
5.0
0.3
50
0.0
-40
-20
40
perturbation
Figure 11: Dynamics of coupled systems with varying coupling strengths. The systems and the parameters used are identical to those of Figure [O (Left) Results for the periodic Lotka-Voltern equations. (Middle) Results for the chaotic Lorenz system. (Right) Loss function landscape of the coupled systems with different coupling strengths.
C Further description of the homotopy optimization procedure
To implement homotopy optimization, usually selects the number of discrete steps for optimiza- tion, as well as series of positive decrement values for the homotopy parameter (*))))) that sum to 1. Afterwards, optimization starts with an initial X value of ( 1, which gives () (0) (()). At each step. the objective function at the current iteration is minimized with respect 0, using the output from the previous step ot (k-1) is the initial guess:
((² (0) H(0, A = (H)) + (k(k) argmin H ())
(8)
Afterwards, decremented to its next value, (k+ = ( A and this iteration continues until the final step s where ( = 0, H(()(()) = F(0). and the final minimizer (*))) = is the sought-after solution to the original problem F(0).
D Multiple shooting algorithm
Our implementation of the multiple shooting algorithm for training NeuralODEs closely mirrors the example code provided in the DiffEqFlux..ii package [54) of the Julia programming language ecosystem.
Given time series data (t,, ú, oN multiple shooting method partitions the data time points into m overlapping segments:
[to (10)
(m
(m
tx]] +((-11
,,), i 1.... 1.
10
18

---

During training, the NeuralODE is solved for each time segment, resulting in m segmented trajecto- ries:
(m (/-
Taking into account the overlapping time point between segments, these trajectories can then all be concatenated produce the full trajectory prediction:
1| [uo.... ....
From this, the data loss is defined identically to conventional NeuralODE training: Latta (f) N+I Σ, (0) ú
However, due to the segmented manner in which the above trajectory is generated, optimizing only over this quantity will result in a model that could generate good piecewise predictions, but is unable to generate proper, continuous global trajectory. Therefore, to ensure the model can generate smooth trajectories, continuity constraints achieved differs across the literature, and our implementation following that of [54) introduces regularization term in the loss:
(f-1) u ug
(4) (i 1 m 1) must be enforced. How this is
Ccontimuity (0) 1
Σ
(1-1)
2 u
Finally, the train loss for the multiple shooting method then defined as
C( 3) = Ldara (0) + B Lcomtinsity (O)
and is minimized using gradient descent, where B is hyperparameter that tunes the relative impor- tance the two terms during training.
E Experiment details
Our experiments were implemented in PyTorch 46], using the pytorch--ightning [16] frame- work. The wandb library (7] was used to keep track of all experiments as well as to perform hyperparameter sweeps. all experiments, the AdamW optimizer (37] was used to minimize the respective loss functions for the vanilla and homotopy training. Each experiment was repeated using random seed values of 10, 20, and 30 compute the prediction means and standard errors.
E.1 Voltera system
The Lotka- Volterra system is a simplified model of predator-prey dynamics given by,
dx = ox Bry, di
dy
Ty + day
where the parameters a, 3, 3,7,6. characterize interactions between the populations.
Data preparation. Following the experimental design of Rackauckas et al. 55 we numeri- cally integrate Equation 99) in the time interval E [0, 6.1], At 0.1 with the parameter values a, 0,,87,68 1.3, 0.9, 0.8. T.8 and initial conditions zr(),, (() 0.442-49296, 4.6280594. Continu- ing with the recipe, Gaussian random noise with zero mean and standard deviations with magnitude of 5% of the mean of each trajectory was added to both states. For both data generation and Neu- raIODE prediction, integration was performed using the adaptive step size dopri5 15 solver from the torchdiffeg package [1I] with an absolute tolerance of lc-9 and a relative tolerance of le-7.
For the experiments of Figure 3 training data were generated in similar manner, but with specific parameters varied to match the corresponding experiments. Data varying train data length (Figure 3 left panel) were generated using time spans of E [0. 3.1). [0, 6.1], [0.9.1] respectively with shared parameter values of At 0.1 and noise amplitude of 5% of the mean. Data with differing sampling periods (Figure 3 third panel) were generated with a fixed time span / [0, 6. noise amplitude of 5% of the mean and sampling periods of A/ 0.1, 0.3, 0.5. Data with different noise amplitudes (Figure 3 fourth panel) were generated using a time span & [0, 6.1). sampling period A/ 0.1 and noise amplitudes of 5%, 10%, 20%, 50% of the mean.
19

---

Model architecture. We used two types of models for this dataset, black-box NeuralODE of Equation ( and a gray-box NeuralODE used in Rackauckas et al. [54] that incopporaees partial information about the underlying equation, given by: dy
dx dt oz + ( yr 01), dt - (2(,, y:; 02).
(10)
In our default setting, the black-box NeuralODE had 3 layers with [2. 32, 2] nodes respectively. The gray-box NeuralODE had 4 layers with [2, 20, 20, 2] nodes, with each node in the output layer corresponding to the output of U1. U2 of Equation (1O). Following the results from [I], a non- saturating gelu activation was used. For the model capacity experiment (Figure second panel), the number of nodes in the hidden layer was changed accordingly.
used for all layers except for the final layer, where identity activation
E.2 Double pendulum
The double pendulum system is a canonical example in classical mechanies, and has four degrees of freedom 01, 82, w1, 42 corresponding to the angles and angular velocities of the two pendulums with respect to the vertical. The governing equation for the system can be derived using Lagrangian formulation of classical mechanics and is given by:
d01
d02
dt W1,
W2
(11)
di
dwi
mylww sin A cos A + m3l2w3 sin A + m2g sin 02 cos A (m1 + m2)g sin 01
dt
(mi + m2)4 m2l1 cos Ao
(12)
dwz m2lzis sin A0 cos A0 + (m1 m2 1- sin A0 + g sii 0₁ cos A0 g sin s2)
(13)
(m1 + m2 - m2l2 cos? AO
where A0 = 02 01, m., m2 are the masses of each rod, and y is the gravitational acceleration.
= 4
Data preparation. While simulated trajectories for the double pendu- lum can be generated using the equations above, we instead used the experimental data from Schmidt & Lipson [59]. This consists of two trajectories from the double pendulum, captured using multiple cam- eras. The noise in the data is subdued due to the LOESS smoothing performed by the original authors. For our experiments, we used the first 100 points of the first trajectory for training and the next 100 to evaluate the extrapolation capabilities of the trained model.
4
02
02
1/
Model architecture. Two types of models were used for this dataset: black-box NeuralODE (Equation (T) and NeuralODE with second- order structure 23), which for this system, takes the form:
0x
d0; Wi, dt
dw,
U. (0,,, 02 W1, W2; 0:); i= 1,2.
Figure 12: Diagram of a double pendulum.
dt
The black-box NeuralODE had 4 layers with [4, 50, 50, 4] nodes, with the input and output node numbers corresponding to the degrees of freedom of the system. The second-order model also had 4 layers with [4. 50. 50. 2] nodes. Note that there are now 2 output nodes instead of 4 since incorporating second-ordd structure requires the neural network to only model the derivatives of Wi and w2. For both models, gelu activations were used for all layers except the last, which used identity activation. Identical to the previous dataset, the NeuralODEs were integrated using the dopri5 solver from the torchdiffeg package [I] with an absolute tolerance of le-9 and relative tolerance of le-7,
E.3
Lorenz system
The Lorenz system is given by the equations
dx dt = oly - xx,
dy
dz
x(p 2) y,
zy B2.
(14)
dt
dt
20

---

Data preparation. To generate the data, we followed experimental settings of Vyasarayani et al. [66] and used parameter values o, p. 10, 28, 8/3 and the initial condition zo. No. 40 1.2, 2.1, 1.7. These conditions, upon integration, gives rise to the well-known "butterfly attractor". The training data was generated in the interval E [0, 3. and still adhering the paper, a Gaussian noise of mean O and standard deviation 0.25 was added to the simulated trajectories to emulate experimental noise. For data generation and NeuralODE prediction, the adaptive step size dopri5 solver from the forchdiffed eq package [T] was used with an absolute tolerance of le-9 and relative tolerance of le-7. ODE solver and the tolerance values were kept identical to the previous Lotka-Von Volterra experiment.
Model architecture. We used a single NeuralODE for this dataset a black-box model having layers with [3, 50, 50, 3] nodes. The number of nodes for the input and output layers correspond to the three degrees of freedom of the state vector u = [x, y, 2] of the system. The activations for the model was kept identical to the previous experiments gelu activation for all layers, except for the last layer which used identity activation.
F Hyperparameter selection
For all combinations of models and datasets used, hyperparameters for the optimization were chosen by running sweeps prior to the experiment with a fixed random seed of 10, then selecting hyperparameter values that resulted in the lowest mean squared error value during training.
Note that for the experiment of Figure 3 the same hyperparameters used in Figure I was used for all experiments, due to the sheer cost of sweeping for the hyperparameters for every change of independent variables.
Vanilla gradient descent. Vanilla gradient descent has learning rate as its only hyperparameter. Due to finite computation budget, we set the maximum training epochs to 4000 for all sweeps and experiments. We list the values used in the hyperparameter sweep well as the final selected values for each dataset in Table [ We found that larger learning rates lead to numerical underflow in the adaptive ODE solver while training the more difficult Lorenz and double datasets; hence, the sweep values for the learning rate parameters were taken to be lower than those for the Lotka- Volterra dataset.
Dataset
Model
Hyperparameter
Sweep values
Selected value
Black-box
0.05
Lotka-Volle Volterra
Learning rate
0.005, 0.01, 0.02, 0.05, 0.1
Gray-box
0.02
Black-box
0.02
Double pendulum
Learning rate
0.002, 0.000, 0.01, 0.02
Second-order
0.05
Lorenz
Black-box
Learning rate
0.002, 0.005, 0.01, 0.02
0.005
Table 1: Hyperparameter sweep and selected values for vanilla gradient descent
Multiple shooting. The multiple shooting method has three hyperparameters: learning rate, number of segments, and continuity penalty. For the majority of experiments, we set number of segments to 5, and performed hyperparameter sweep the remaining two parameters.
Homotopy optimization. Our homotopy- based NeuralODE training has five hyperparameters: learning rate (17), coupling strength (k), number of homotopy steps (s), epochs per homotopy steps (Nepoch), and homotopy parameter decrement ratio (x).
While sweeping over all five hyperparameters would yield the best possible results, searching such high dimensional space can pose a computation burden. In practice, we found that sweeping over only the first two hyperparameters and fixing the rest at predecided values still returned models that functioned much better than their vanilla counterparts. We list the predecided hyperparameter values for each dataset in Table 3 Note that we increased the epochs per homotopy step from 100 the Lotka- Volterra model to 300 for the Lorenz system and double pendulum datasets to account for the increased difficulty of the problem.
21

---

Dataset
Model
Hyperparameter
Sweep values
Selected value
Learning rate Continuity penality
0.005, 0.01, 0.02, 0.05
0.05
Black-box
0.005, 0.002, 0.01
0.01
VotterVolterra
Learning rate
0.005, 0,01, 0.02, 0.05
0.05
Gray-box
Continuity penalty
0,005, 0.002, 0.01
0.002
Learning rate
0.005, 0,01, 0.02, 0.05
0.02
Black-box
Continuity penalty
0.005, 0.002, 0.01
0.002
Double pendulum
Learning rate
0.005, 0.01, 0.02, 0.05
0.02
Second-order
Continuity penalty
0.005, 0.002, 0.01
0.002
Learning rate
0.005, 0.01, 0.02, 0.05
0.02
Lorenz
Black-box
Continuity penalty
0.005, 0.002, 0.01
0.01
Table 2: Hyperparameter sweep and selected values for multiple shooting
Dataset
Model
Hyperparameter
Predecided value
Number of homotopy steps
6
Black-box
Lotka-Volteraa
Epochs per step
100
Gray-box
decay ratio
0.6
Number of homotopy steps
7
Black-box
Epochs per step
300
decay ratio
0.6
Double pendulum
Number of homotopy steps
6
Second-order
Epochs per step decay ratio
300
0.6
Number of homotopy steps
8
Lorenz
Black-box
Epochs per step
300
decay ratio
0.6
Table 3: Predecided hyperparameter values for homotopy optimization
Table shows the sweeped hyperparameters well the final chosen values for each dataset. The total number of training epochs is given by multiplying the number of homotopy steps with epochs per homotopy steps. This resulted in 600 epochs for the Lotka- Volterra dataset, and 1800 epochs for both the Lorenz system and the double pendulum datasets.
Dataset
Model
Hyperparameter Learning rate Control strength
Sweep values 0.01, 0.02, 0.05
Selected value
0.05
Black-box
3, 4,
6
Lotka-Volters
Learning rate Control strength
0.01, 0.02, 0.05
0.02
Gray-box
2, 5,
4
Learning rate
0.01, 0.02, 0.05
0.05
Black-box
Control strength
3, 4, 6, 7, 8, 9, 10
10
Double pendulum
Learning rate
0.01, 0.02, 0.05
0.02
Second-order
Control strength
3, 4, 5, 6, 7, 8, 9, 10
10
Learning rate
0.005, 0,01, 0.02
0.02
Lorenz
Black-box
Control strength
6. 7, 7, 9, 10
6
Table 4: Hyperparameter sweep and selected values for homotopy optimization
22

---

G Further results
G.1 Table of the benchmark results
Here, we present the benchmark results of Figure 4 in a table format as an alternative data representa- tion.
Table Experimental results for various NeuralODE models trained with three different methods. The values for best train epochs correspond to random seed values of 10, 20, and 30, except for runs marked with where the backup random seed was used due to unstable training the original random seed. MSE values are reported with the mean and standard error from three different runs.
Daiaset
Lotka Lotka-Volterra
Double Pendulum
Lorenz System
Type
Black
Gray
Black Box
Best train epochs
Baseline
(3932, 3818, 3134) (3644, 1482, 1359)
(3999, 3980", 3992) (391 2127, 1198)
(3998, 3731, 3633)
(2890, 3914, 3964)
(3960. 3996, 3958) (3710, 3666*, 3976)
Multi- Shoot.
(3438, 3454, 3646)
(3286, 2207, 3831)
Homotopy
(299, 208, 227)
(291, 265, 216*)
(1201, 2087, 1502)
(1200, 600, 900)
(2399, 2399, 2377*)
Mcan Squared Error (x 10-?)
Baseline (interp.)
1.76+0.13
8..3
166
1.03=0.12
259:+6.05
Multi- Shoot.(interp.)
0.95::0.22
0.74:10.02
0.81:10.33
0.28=0003
12.9:11003
Homotopy (interp.)
0.78:=0005
0.72:#0.01
0,98+0008
0,18=0.02
18.6:t1.75
Bascline (extrap.)
2.73+00117
6254:5022
944.9+159.7
1619#172.99
603+ 1.36
Multi-, (extrap)
1.42::0.01
262.2:422.4
2676::391.7
2365:1105.1
138:437.0
Homotopy (extrap.)
1.22::0099
7.656=24.577
1031:1003.2
428.5=$1.90
92.8++2.59
G.2 Example of the training curves
Here, include some of the training curves for experiments. Figure 13 displays the averaged training curves for the black-box model trained on the VolterVolerrra dataset. From the right panel, it clearly be seen that homotopy method optimizes the MSE much rapidly than the other methods, arriving at the noise floor in less than 500 epochs. The abrupt jumps in the MSE curve for homotopy is due the discontinous change in the train loss that occurs every the the homotopy parameter is adjusted. To make this connection clearer, we show the train loss as well as the homotopy parameter for our method in the left panel of the same figure.
10a
Nokse floor
Homotopy
shooting
Train loss Lie.
10a
107
Vanilla
Mean squared error
10-1
1.2
10-1
1.0
10-1
0.5
10"2
0.0
200
400
600
0
500 10000
15500
2000
2500
3000
3500
4000
Epoch
Epech
Figure 13: Training curves for the black-box model on the Lotka- Volterra dataset, corresponding to the benchmark results of Figure (Left) Train loss L(0, X) and the homotopy parameter X. (Right) Mean squared error as a function training epochs. Note that the curves for homotopy stops early due to the choice of the algorithm hyperparameters (see Appendix for further details).
We present analogous results for the second-order model trained on the double pendulum dataset in Figure [4] Once again, we find that our homotopy method converges much quicker than its competitors.
23

---

Noise floor Homotopy
-shooting
Train L(B.A)
100
Vanilla
Mean squared
10"
10
10"3
10-3
1.0
10-2
0.5
0.0
10-2
300
600
900
1200
1500
1800
500
1000
1500
2000
2500
3000
4000 3500
Epoch
Epoch
Figure 14: Train loss, homotopy parameter, and mean squared error as a function of epochs for the results of Figures 45 and 6
The layout of the plots, well their interpretations are identical to
that Figure
G.3 Further discussion on the Lotka- system results
In this section, we present additional results regarding our experiment of Figure 3 and continue our discussion.
Increasing data length. Figure 15 depicts the predicted trajectories for differently trained models. For the shortest data, it can be seen that all models did overfit on the training data and failed to capture periodic nature the system. This justifies our attributing the large extrapolation error of models in Figure to overfitting. However, we emphasize that this overfitting is not due to failings rather due to the insufficient information the training data, as the
in the training algorithms, models cannot be expected to learn periodicity without being given at least single period worth of data. As the data increased, we find that both homotopy and multiple shooting properly capture the dynamics of the system, whereas the vanilla method was unable properly learn the system for the longest data.
Homotopy
Multishoot
State
5
State
10
15
10
15
10
15
Time 2
Time t
Time t
Figure 15: Predicted trajectories for the models with increasing train data length, t=(3.1. 6.1. 9.1).
Corresponds to the first panel Figure The red dashed line indicates the start of extrapolation.
Reducing model capacity. In Figure 16 we show the prediction trajectories corresponding to the model capacity experiment. clear from the results that both the homotopy and multiple shooting methods can produce models that accurately portray the dynamics, regardless of the reduction in the model size. In contrast, predictions from vanilla training deteriorate as the model size decreases, with the smallest model inaccurately estimating both the amplitude and the phase of the oscillations.

---

7.5
Data
7.5
7.5
Homotopy
State
5.0
Murtishoot
5.0
5.0
Vanilla
2.5
2.5
2.5
0.0
0.0
0.0
6-
6
6-
State
4
4-
2
2
0
0
0
0
5
10
0
5
10
5
10
Time
Time
Time
Figure 16: Training curves for each number of nodes, nodes=(8, 16, 32). Corresponds to the second
panel of Figure
Increasing sampling period. The model trajectories for decreasing sampling period in the training data is shown below in Figure 17] Intriguingly, we find that the multiple shooting method results suffer greatly with the increased data sparsity. This is side effect of our experiment setting, where we set the time interval constant while increasing the sampling period. To elaborate, choice of fixed time interval causes the number of training data points decrease as the period increases. However, as multiple shooting trains by subdividing the training data, it struggles small datasets because this leads to each segment containing even less data points that do not convey much information about the data dynamics.
Another interesting observation that can be made is that vanilla training gives better results when sparser data is used, which is also confirmed by the error values in the third panel of Figure This suggests an alternate training improvement strategy for training NeuralODEs on long time series data by training on subsampled data, then if necessary, gradually anneal the data intervals back to its original form as training proceeds. Of course, the effectiveness of such a scheme will need to be tested for more complex systems, which is outside the scope of this paper.
M
7.5
7.5
Homotopy
States
5.0
Multishoot Vanilla
5
2.5
0.0
0.0
State
4-
2.5
2
2
5
5
10
Time
Time
Time
Figure 17: Predicted trajectories for the models with increasing sampling periods, dt (0.1, 0.3, 0.5). Corresponds to the third panel of Figure 3
As our homotopy training method utilized a cubic smoothing spline to supply the coupling term, we also inspected the quality of this interpolant as the data period was increased. From Figure [8] we find that the cubic spline relatively resistant to the increased data sparsity, which may be one of the reasons why our homotopy method is very robust against this mode of data degradation.
25

---

10
10
10
Data
Ground truth
Cubic spline
5
6-
6
State
4
2
0
2
6
2
Time
Time
Time
Figure 18: Cubic spline interpolation results for data with increasing sampling periods, dt=(0.1, 0.3,
0.5). Corresponds to the third panel of Figure
Increasing data noise. Finally, Figure 19 shows the model trajectories as the noise the training data is increased. Here, we do find that while all methods return deteriorating predictions as the noise increased, our homotopy method found be most robust to noise, followed by multiple shooting, and finally vanilla gradient descent. A part of our algorithm's robustness to noise be explained by the use of the cubic smoothing splines, As our use of synchronization couples the NeuralODE dynamics to the cubic spline trajectory during training, one could argue that the denoising effect of cubic splines is thus directly transferred to the model.
Data
Homotopy
Multishoot
5
Vanilla
State
2-
2
10
10
5
10
5
10
Time
Time t
Time t
Time t
Figure 19; Predicted trajectories for the models with increasing noise, noise=(0.05, 0.1, 0.2, 0.5).
Corresponds to the fourth panel Figure
However, this does not seem to be the full picture. Figure 20 shows the cubic smoothing splines used during modeling training for each the noise levels. As we held the amount of smoothing constant throughout our experiments, we see that for larger noise amplitudes, the spline fails to reject all noise, and displays irregular high frequency behaviors. On the other hand, the corresponding trained trajectories for the homotopy method (Figure 19 third and fourth pancls) do not mirror such irregular oscillations, indicating that the homotopy optimization procedure itself also has an intrinsic robustness to noise.
26

---

10
10
10
10
Data
Ground truth Cubic Subicespline
State
4
4
6
2
6
Time
Time
Time
Time
Figure 20: Cubic spline interpolation results for data with increasing noise amplitude, Noise=(0.055.
0.1, 0.2, 0.5). Corresponds to the fourth panel of Figure
G.4 Comparison to non-deep learning algorithm: SINDY
In this section, we briefly compare our algorithm to a more traditional symbolic regression- -base method. We choose the well-known SINDy algorithm, which is readily available in the pysindy package. This algorithm takes as input, time series measurements of the states the dynamical system in question, well dictionary of candidate terms that could constitute the governing equation. Afterwards, the time derivative of the states is estimated numerically, and a sparse regression is performed to determine which terms exist in the governing equation, well what their coefficient values are.
Compared to neural network-based approaches, the unique feature of SINDy is that its results are given in the form of analytical expressions. Therefore if the governing equation for the data is simple (e.g. consists of low-order polynomial terms), or if the user has sufficient prior information about the system and can construct a very small dictionary that contains all the terms that appear in the true equation, SINDy can accurately recover the ground truth equation.
On the other hand, if the system high dimensional so that the space of possible candidate terms start to grow drastically, or if the governing equation has arbitrary. complex form that cannot be casily guessed, SINDy either does not fit the data properly, or even if does, recovers inaccurate equations that are difficult to assign meaning to. These correspond to situations where NeuralODEs, and hence our method for effectively training them, shines.
To illustrate this point, we choose the double pendulum dataset and compare the predicted trajectories from SINDy and homotopy algorithm. As one can see from Equation [3). the governing equation of this system has an extremely complicated form that cannot be guessed easily. To reflect this difficulty, we chose the basis of the dictionary of the candidate terms to be the state variable themselves, their sines and cosines, and the higher harmonics of the sines and cosines.
1.5
7.5
20
Data
5.0
1.0
10
0.5
State 8;
0.0
0.0
State
State
-0.5
1.0
-5.0
-5
100
-1.0
-10.0
-15
-2.0
Figure 21: Comparison between our homotopy method and SINDy on the double pendulum dataset.
27

---

From Figure 21 we find that SINDy is able to fit the data to modest accuracy. While this may seem weak compared to the result from our homotopy method, the accuracy of the SINDy prediction is likely to improve if the candidate dictionary is updated to better reflect the actual form of Equation (T3). However, this amounts to supplying the algorithm with excessive amounts of prior information about the system, which is rarely available in practice. This is also constrasted with NeuralODE training with our homotopy method- which does not require any additional information about the system to produce the results above.
G.5 Changing the gradient calculation scheme
As we commented in Section 3 methods that focus on other aspects of NeuralODE training. such as alternate gradient calculation schemes, are fully compatible with our homotopy training method. To demonstrate this point, we used the symplectic-adii torch-symplecticc a--ajdoonnn library as substitute for the direct backpropagation used for gradient calculation our experiments.
method [39] provided in the
Table Benchmark results with the gradient calculation scheme changed the symplectic-adjoint method.
Dataset
Lotka-Volters.
Double Pendulum
Lorenz System
Model Type
Black Box
Gray Box
Black Box
Second Order
Black Box
Best train epochs
Baseline Multi-, Shoot, Homotopy
(3999, 3999, 3853) (3982, 3968, 3670) (299, 208, 228)
(3997, 3999, 3992)
(3909, 3384, 3998)
(3983, 3972)
3972)
(3994, 3999, 3990) (3977, 3418, 3954)
(3901, 3958, 3996)
(3805, 3913, 3820)
(3988, 3970. 3942)
(898, 291)
(1201, 1799, 1774)
(900, 1799,
(1799, 918, 1799)
Mean Squared Emor (x 10 2)
Baseline (interp.)
2.47++0.13
187+133
1.69-0.20
1.02+-0.12
184+-60.1
Multi-, Shoot.(interp.) )
1.20::0.02
0.95=0.02
0.76=0.01
0.61
29.1+15.8
Homotopy (interp.)
0.77=1004
1.68
1.08=0.12
0.25::0.03
33,325,43
Comparing the above results to our benchmark results in FigureMaan and Table 5 we see that the overall results remain the same our homotopy method is able to train the models effectively, in much smaller number of epochs.
28


Reviews:
[{"review_id": "4bKrGbRT4G", "reviewer": "Reviewer_9EGS", "summary": "This paper presents a new training method for Neural Ordinary Differential Equations (NeuralODEs) that aims to improve their performance in extracting dynamical laws from time series data. The proposed method is based on synchronization and homotopy optimization, which does not require changes to the model architecture. The authors demonstrate that their method achieves competitive or better training loss while often requiring less than half the number of training epochs compared to other techniques. Furthermore, models trained with their method display better extrapolation capabilities.", "questions": "see above", "limitations": "see above", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "[+] The paper is well-organized and clearly written, which makes it easy to follow and understand the proposed method. The introduction to the relevant background is sufficient, and the presentation of the experimental part is well-structured.\n\n[+] The proposed method is well motivated, and some visualizations play a significant role in understanding the problems with traditional methods and the performance of the newly proposed method.\n\n[+] The newly proposed method is very easy to integrate with models. In fact, the new method does not change the model but only changes the training process. It is simple to use and has excellent results.\n\n[-] Considering that there are five main hyperparameters, if there are targeted sensitivity experiments in the experiments to demonstrate how the selection of hyperparameters affects the experimental results to some extent, the experimental results might be more convincing.\n\n[-] The few systems involved in the experimental part are low-dimensional ODE systems. Can the proposed algorithm be applied to higher-dimensional problems? What new challenges might be encountered?\n\n[?] What is the relationship between the K in the algorithm framework and the hyperparameters: Coupling strength (k)?\n\n[?] It would be better if there is an analysis of computational complexity or a comparison of computation time with the baseline algorithm.\n", "weaknesses": "see above", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the detailed comments and suggestions. We have written our responses to some of the common questions in the global comment above, and have listed responses to the rest of the reviewer's questions below.\n\n> Considering that there are five main hyperparameters, if there are targeted sensitivity experiments in the experiments to demonstrate how the selection of hyperparameters affects the experimental results to some extent, the experimental results might be more convincing.\n\nWhile we agree with the reviewer, due to limited computational resources, we were only able to sweep a subset of the full hyperparameter space while running our experiments. Still, to provide some more information, we have included training curves for our preliminary runs where we individually varied some of the hyperparameters in Figure 3 of the pdf attached in the global comments. In case of the homotopy parameter decrement ration $\\kappa$ and the number of homotopy steps, a suboptimal choice causes the train loss to jump in the later stages of the optimization due to the loss landscape changing too drastically after a homotopy step. In case of the learning rate, larger values cause the train loss to decrease rapidly in the initial stages, but also likely to stick in a bad local minimum later on during training. For the coupling strength, larger values cause the training loss to start from a smaller value and facilitate stable optimization. However, using too large values cause the coupled train loss to act completely independent of the true mean squared error due to insensitivity of the model parameters to the train loss.\n\n> It would be better if there is an analysis of computational complexity or a comparison of computation time with the baseline algorithm.\n\nFirst regarding the computational complexity, our homotopy algorithm has almost identical complexity as the vanilla training of NeuralODEs. This is because our algorithm only changes the differential equation to be solved (by adding on a coupling term to the model). The major factors of computational complexity are how the differential equation is solved and how the gradients are calculated (direct backpropagation or adjoint method-based techniques) - which were kept identical for both the homotopy and vanilla methods in our experiments.\n\nOne difference between homotopy and vanilla training arises from the calculation of the cubic spline interpolant. However, this is calculated only once prior to the training and the resulting interpolations are cached for use during training. Therefore this only adds a negligible constant factor to the total algorithm. \n\nWhile this was not a direction explored in our paper, in the case of the multiple shooting, the computational complexity can be made different if one follows the approach of [1]. In this implementation, all of the segments for multiple shooting are solved in parallel, bringing down the computational complexity by a factor equal to the number of multiple shooting segments used. However, we would like to clarify that this does not always guarantee shorter training time there is a competition between each epoch of the multiple shooting method taking shorter time and our homotopy method requiring much less total epochs to converge. The use of adaptive solvers (which we describe below) further obfuscates this problem, and a further investigation is required for a final verdict.\n\nRegarding the computation time, we want to preface this point by saying that we did not account for such analysis at the time of our experiments, and thus our implementations are unfortunately not streamlined for speed. In the case of the homotopy method, some part of the cubic spline coefficients were recomputed each epoch and our multiple-shooting implementation was not parallelized in time. Furthermore, we had used multiple machines to run some of the longer experiments in parallel, which made it difficult to provide a comprehensive computational time benchmark. We apologize for this inconvenience and will improve on this point in future works.\n\nWith this, we present the computation time benchmarks for the Lotka-Volterra training results in Table 1 of the attached pdf in the global comment. We see that while our homotopy method does take longer time per epoch than the other methods, the time reached to reach best epoch is much shorter due to the much less number of epochs required to train the model. A point worth stressing is that since we used 5 segments for the multiple shooting training, the expected speed up from parallelizing the multiple shooting algorithm is about 5x. However, our homotopy method takes less than this time to converge, indicating that the our method is competitive with even a optimized,  parallelized version of the multiple-shooting algorithm.\n\nFinally, we want to comment on one factor, other than the algorithm time complexity, that can greatly affect the computation time. This is due to the use of the adaptive ODE solver for solving the NeuralODE. As adaptive solvers evaluate the ODE variable amount of times in order to reach given tolerance values, a method that has a better per-epoch time complexity can show worse computation times if the solver evaluates the ODE more times for that method. This number of function evaluations is itself dependent on the form of the differential equation - i.e. the parameters of the NeuralODE. Therefore, a training method that guides the NeuralODE through a more favorable region of the parameter space will result in lesser number of function evaluations throughout training, which can result in reduced total training time even if the method’s per-epoch time complexity is worse than its competitors.\n\n[1] M. Stefano et. al, Differentiable Multiple Shooting Layers. NeurIPS 34 (2021)."}, {"review_id": "0g4ZDaaPtN", "reviewer": "Reviewer_jYQA", "summary": "The paper proposes a homotopy-based training method for NeuralODE models, in particular for the case with cases with long sequence of training data.  Comprehensive experimental results demonstrated the effectiveness of the proposed method. ", "questions": "1. I question the statement around Line 196, on the usage of cubic splines and its potential consequences. Cubic spline is a relatively high-order method that may be inaccurate when the underlying dynamic system is **stiff** (or there are widely spread time constants). Will the cubic spline interpolation be able to handle responses from _stiff_ systems? \n\n2. I question the convergence of the homotopy method when constant $\\Delta \\lambda$ is used. The introduction of homotopy to the original ODE effectively increases the order of the ODE by 1, and transform it into a DAE.  The authors should address whether there is the potential of non-convergence, and when it happens, how to address it. ", "limitations": "No negative societal impact perceived.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. The proposed method is based on homotopy method, which is a classic global optimization method, with sound theoretical foundation.\n2. Comprehensive experimental results. Good results when compared with the baseline methods\n3. The paper explains key concepts and results well. ", "weaknesses": "1.  The argument on the loss function landscape lacks theoretical insights. The connection is established by empirical results. \n2.  Although the overall idea of using homotopy is good, some details are not properly addressed. See the questions and limitation part for more details. \n3. A few key references related to homotopy methods are not cited. The references by JH He really put homotopy methods in the mainstream of scientific computing: https://scholar.google.com/citations?user=tzM7c2cAAAAJ&hl=en\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the comments and the very interesting questions. We list our answers below, and have grouped the common questions in the global comment.\n\n> A few key references related to homotopy methods are not cited. The references by JH He really put homotopy methods in the mainstream of scientific computing: https://scholar.google.com/citations?user=tzM7c2cAAAAJ&hl=en\n\nThank you for kindly pointing us to key references regarding the homotopy method. We will update the manuscript to include these references in the related works section.\n\n> I question the statement around Line 196, on the usage of cubic splines and its potential consequences. Cubic spline is a relatively high-order method that may be inaccurate when the underlying dynamic system is **stiff** (or there are widely spread time constants). Will the cubic spline interpolation be able to handle responses from *stiff* systems?\n\nWe thank the reviewer for bring up an interesting discussion point. To address the reviewer’s comment, we created simulation data from the Robertson’s chemical reaction system and constructed a cubic spline interpolant from it. The results are displayed in Figure 2 of the pdf linked in the global comment. The results show that in this case, cubic spline interpolation does return satisfying results, except for a slight ripple after the first peak in the y variable.\n\nStill, we do agree that it is possible that on some time series datasets, such as those from more pathological stiff systems or data with large amounts of noise, cubic spline interpolation may return sub-par results. However, our experiences seem to indicate that our training method is relatively robust to the quality of the interpolant. \n\nTo support this point, we would like to guide the reviewer to Figure 19 and Figure 20 of our appendix (please look at the appendix.pdf file in the zip file of our Supplementary materials as it is the updated version compared to the one at the end of the main paper pdf), where we plot trained NeuralODE predictions, as well as the cubic spline interpolants used in the training for increasing noise in the data. We find that when the noise amplitude becomes large, the cubic spline interpolant displays high frequency artifacts (Figure 20, rightmost panel). However, the corresponding model predictions (Figure 19, rightmost panel) show no signs of this artifact carrying over, and still outperforms all other baseline methods.\n\n> I question the convergence of the homotopy method when constant $\\Delta \\lambda$ is used. The introduction of homotopy to the original ODE effectively increases the order of the ODE by 1, and transform it into a DAE. The authors should address whether there is the potential of non-convergence, and when it happens, how to address it. \n\nWe thank the reviewer for bringing up another important point. We agree with the reviewer that our homotopy training method does not guarantee convergence to the global minimum. The reason for this is because homotopy optimization method, while effective, is still a local optimization method. One possible resolution to this problem is to use the HOPE algorithm from [2], which is a combination of homotopy optimization and particle swarm optimization. By employing multiple \"particles\" to explore the parameter space, the algorithm has a much higher chance to converge to the global optimum at the cost of increased computation.\n\nFinally, regarding the reviewer’s comment about the homotopy turning the original ODE into a DAE of 1 order higher, we are less familiar with DAEs and are unsure if this is indeed the case. The reason behind our confusion is that in the original ODE, the states evolve with time t, whereas the homotopy parameter changes with the training epochs. Therefore, we are not sure if our situation can be cast into the typical form of DAEs, given by\n\n$$\n\\frac{d\\bf u}{dt}={\\bf f}(t, \\bf u, v);\\quad 0={\\bf g}(t, \\bf u, v)\n$$\n\n\nWe would love to discuss more about this point and learn from the reviewer.\n\n[1] H. H. Robertson, The solution of a set of reaction rate equations, *Numerical Analysis: An Introduction*, pp. 178–182, Academic Press, London (1966).\n\n[2] D. M. Dunlavy and D. P. O'Leary, *Sandia Reports* SAND2005-7495 (2005)."}, {"review_id": "BegSgvU87f", "reviewer": "Reviewer_bK9U", "summary": "The authors present a new method for training Neural Ordinary Differential Equations (NeuralODEs), a modeling approach that merges neural networks with the paradigm of differential equations from physical sciences. While NeuralODEs offer significant potential for extracting dynamic laws from time series data, they often suffer from long training times and subpar results. In this paper, the authors propose a training method built on synchronization and homotopy optimization. Their results demonstrate that the new method achieves competitive or even better training loss, often requiring fewer epochs compared to traditional model-agnostic techniques, and also exhibits better extrapolation capabilities.", "questions": "1. How does the efficacy of homotopy and multiple-shooting methods change with the size of the NeuralODE model? Is it possible to find a balance between model size and training method for optimal performance?\n2. What are the effects of data sparsity and noise on the performance of different training methods for NeuralODEs, particularly the homotopy training method?\n3. How can this methodology be extended to high-dimensional systems or systems with partially observed states?\n4. What are the practical implications and potential applications of this training method in various industries or scientific disciplines?", "limitations": "The paper mainly uses synthetic datasets for testing the performance of training algorithms. The findings might not generalize to real-world datasets, which can exhibit unique complexities and noise structures.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper addresses a notable challenge in the field of NeuralODEs: the long training times and subpar results. Their proposed solution is both innovative and logical, introducing a novel concept in the NeuralODE literature.\n2. The authors demonstrate the advantages of their method convincingly, showing that it improves the training loss and requires fewer epochs for training.", "weaknesses": "1. The method's extension to high-dimensional systems or partially observed states is yet to be explored, limiting its current applicability.\n2. The paper could have benefited from more real-world application examples, demonstrating the utility of the method beyond the context of benchmark experiments.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the comments as well as the detailed questions. We have listed the responses to the questions below and grouped the responses to some of the more common questions in the global comment above. \n\n> How does the efficacy of homotopy and multiple-shooting methods change with the size of the NeuralODE model? Is it possible to find a balance between model size and training method for optimal performance?\n\nWhile more experiments comparing the two methods are needed for a final confirmation, we found in our experiments that the relative efficacy of the two methods were not very sensitive to the model size. Both methods returned orders of magnitude improvements compared to the vanilla training albeit homotopy method resulting in better extrapolation performance than the multiple-shooting method.\n\nHowever, continuing the discussion, we do think there is a different factor that can change the relative efficacy of the two methods -  the length of the training data. This is because the multiple-shooting method can be implemented to compute all of the shooting segments parallel in time [1] - resulting in a time complexity proportional to the length of one segment - whereas the homotopy method scales with the length of the total data.\n\nTherefore, for short to long time series data, homotopy training method tends to be more effective than multiple shooting because even though a single epoch of the homotopy method can take longer than multiple shooting, its faster convergence and improved extrapolation qualities makes homotopy the more attractive option. For very long sequences however, multiple-shooting method can be more economical as the reduced time complexity starts to outshine the improved convergence speed of the homotopy method. This is a very interesting direction of research that we would like to pursue in the future.\n\nAnother potentially interesting direction is to use the homotopy method and the multiple-shooting method at the same time - that is, subdivide the data into segments and apply synchronization on each - as the two frameworks are compatible with each other. This has the potential to reap the benefits of the two methods at the cost of increased number of parameters (the sum of the hyperparameters of the two methods). \n\n> What are the effects of data sparsity and noise on the performance of different training methods for NeuralODEs, particularly the homotopy training method? \n\nThank you for your question. Regarding the performances of the different training methods with respect to data sparsity and noise, we would like to refer the reviewer to the last paragraph of Section 7.1 (corresponding to third and fourth columns of Figure 3) and Section G.3 in the appendix with emphasis on Figures 17-20 (please look at the appendix.pdf file in the zip file of our Supplementary materials as it is the updated version compared to the one at the end of the main paper pdf). \n\nFirst, regarding noise, we find that the vanilla method is very sensitive with the model predictions degrading rapidly with increasing noise. In comparison, both the homotopy and the multiple-shooting methods are much more robust, although for very large noise amplitudes our homotopy method prevails. An interesting point here is that large amounts of noise can introduce high frequency artifacts to the cubic spline interpolant used in homotopy training (Figure 20, rightmost column) but the corresponding training results (Figure 19, rightmost column) are not affected by these artifacts.\n\nAs for data sparsity, we find that multiple-shooting method suffers the most as data sparsity increases. This is intuitive because as multiple-shooting method subdivides the time series data, it works well as large, dense data regime where each shooting segment is given enough training data elucidate the system dynamics. Another unexpected phenomenon we observed that the vanilla method tends to work better as the data sparsity increases. This hints at another potential training improvement strategy where one first trains a NeuralODE with an undersampled version of the data, then subsequently training the model on denser and denser version of the data to arrive at the final training result. Finally, we find that our homotopy method is quite robust to data sparsity - model predictions do degrade as data sparsity is increased, but still the training model predictions are competitive with those from vanilla training. As the homotopy method produces such results at an accelerated training time, we believe homotopy method is still better in this situation.\n\n\n> What are the practical implications and potential applications of this training method in various industries or scientific disciplines?\n\nOur homotopy training method is valuable in multiple scientific disciplines where the goal is to discover governing equations underlying data or need to create a surrogate dynamics model that provides reliable forecasts of the data. \n\nOne example is [2], where the authors introduced a gray-box model to predict the temperature of the ocean as a function of depth. The neural network component in this hybrid model serves to account for the missing physics that are not captured by traditional physics based modeling. \n\nTo obtain a proper surrogate model, an effective training method is necessary and this is where our homotopy method shines. Especially, our results on the gray-box models for the Lotka-Volterra system (main paper, Figure 4) show that even though significant prior knowledge on the governing equation was added to the NeuralODE model, vanilla training was not able to fully capitalize on this and instead returned models the performed significantly worse than their black-box counterparts. Our homotopy method can help resolve this problem, thus aiding in creating proper surrogate models.\n\n[1] M. Stefano et. al, Differentiable Multiple Shooting Layers. *NeurIPS* **34** (2021).\n\n[2] A. Ramadhan et. al. *arXiv:2010.12559* (2023)."}, {"review_id": "OnLS5Bz6hR", "reviewer": "Reviewer_84Va", "summary": "Training neural ODE models on long sequences of data from a dynamical system is difficult. The authors argue empirically that this is due to a poorly conditioned loss landscape leading to difficulties in optimization. To rectify this difficulty the authors use tools from the literature on synchronization (which gives conditions for when two dynamical systems will converge based on an error correction coupling term. They provide a homotopy between the synchronized dynamical system and the original one and argue that the optimization dynamics will be nice on the sychronized system and progressively less nice along the homotopy to the original system. However by gradually training the parameters and moving along the homotopy the method “ratchets up the difficulty” in a way that makes the ultimate optimization problem easier. \n\nThey demonstrate the success of the system on a collection of commonly used dynamical systems taken from the neural ODE literature when compared to vanilla SGD and a multiple-shooting algorithm that handles the longer time horizons by breaking them up. The authors show that their method archives better accuracy, robustness, and convergence than the baselines on these problems.", "questions": "How is K chosen? I couldn’t find it in the main paper and the idea is intuitive even if K=I but I’d like to know?\n\nWould this method work for control systems as well? Seems like the counterfactual nature of the problem could be interesting. \n", "limitations": "This seems pretty general and seems to widely apply to the breadth of neural ODE architectures and dynamical systems.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The work combines the idea of synchonization with the standard Neural ODE training to improve the performance on long-horizon data. Since “long-horizon” data is not that long when training Neural ODEs, this is a welcome development. The idea of synchronization is intuitive and theoretically justified as well as empirically demonstrated to be useful through the loss landscape Hessian. The empirical results are not broad but they are deep; the authors study many variations of the problem on a small number of dynamical systems (3) against a small number of baselines (2). However, the method performs well under noise, sparsity, and long data on these and with better test performance and more quickly. \n", "weaknesses": "The theoretical explanation leaves a little bit lacking. I would love to know more about the properties of the dynamical system as it is more or less synchronized. If there were traces of such in an image I think it would be intuitive. It took me a while to work out that as \\lambda -> \\infty any setting of the parameters is optimal and I think it might be good to explicitly say / visualize that. Any more detailed theoretical statements about this method would also be welcomed. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the interest in our paper and the positive comments. We have listed the answers to the reviewer's question below, and have grouped the answers to some of the more commonly occurred questions from multiple reviewers in the global response above.\n\n> The theoretical explanation leaves a little bit lacking. I would love to know more about the properties of the dynamical system as it is more or less synchronized. If there were traces of such in an image I think it would be intuitive.\n\nWe had originally tried to illustrate the effect of coupling strength on the dynamics trajectories through Figure 10 in Section B of the Appendix  (please look at the appendix.pdf file in the zip file of our Supplementary materials as it is the updated version compared to the one at the end of the main paper pdf). In the figure, we perturb the coefficients of the Lotka-Volterra and Lorenz systems, couple them to a reference trajectory, then plot the resulting dynamics trajectory as well as the loss landscape as a function of coupling strength. We hope this can provide some additional information about our work.\n\nWe would also love to add additional visualizations the reviewer finds would be useful and look forward to the reviewer’s comments. \n\n> It took me a while to work out that as $\\lambda \\rightarrow \\infty$ any setting of the parameters is optimal and I think it might be good to explicitly say / visualize that.\n\nIn our paper, we wanted to convey this point through the upper right panel of Figure 2, where we drew the loss landscape for increasing values of the coupling strength. We will update this part to explicitly mention that the flat landscape for the large value of k means that any setting of the parameters is equally optimal. Also to better visualize the flatness of the landscape, we will adjust the y-axis limits as well.\n\n\n> How is K chosen? I couldn’t find it in the main paper and the idea is intuitive even if K=I but I’d like to know?\nWe apologize for the confusion. In our work, the matrix $\\bf K$ and the scalar algorithm hyperparameter $k$ are related to each other by ${\\bf K}=k\\bf I$, where $\\bf I$ is the identity matrix.\n\nWe would like to expound on this by clarifying that this is not the only way to construct the coupling matrix – any choice of $\\bf K$ can achieve synchronization provided that the matrix $\\frac{\\partial \\bf U}{\\partial \\bf u}-\\lambda\\bf K$ is always negative definite along the trajectory ${\\bf u}(t)$. However, naively constructing $\\bf K$ in a different manner results in additional hyperparameters corresponding to each of the independent elements of $\\bf K$. Therefore, we decided to reduce $\\bf K$ into a single scalar value to minimize the number of hyperparameters in our algorithm.\n\nWith that being said, we do think that exploring different parametrization of $\\bf K$ is an important future direction of our research. In systems where the order of magnitude of the different degrees of freedom differ significantly, it can be more intuitive to use ${\\bf K}=diag(k_1, k_2 , \\dots)​$.  Also in the chaos synchronization literature, [1] reports that $K = BV^T$ where $B, V$ are constant vectors can be tuned to achieve synchronization when the parameters of the coupled systems differ greatly.\n\n> Would this method work for control systems as well? Seems like the counterfactual nature of the problem could be interesting.\n\nWe thank the reviewer for the very interesting discussion point. We believe that our method can work on control systems, but the added complexity of the problem will require additional experimentation and theoretical studies to accurately confirm our hypothesis.\n\nTo elaborate, applying our homotopy method to control system is straightforward since everything can be kept the same and only the NeuralODE slightly modified to accept the control signal ${\\bf v}(t)$:\n\n$$\n\\frac{d\\bf u}{dt}={\\bf U}(t,\\bf u, v;\\theta)-K(u-\\hat u)\n$$\n\nwhere $\\bf \\hat u$ is the interpolant constructed from measured time series data ${\\bf u}_{data}(t_i)$.\n\nAs for whether this scheme will actually result in improved training, we can take some hints from the synchronization literature. In [2], synchronization and homotopy optimization were applied to estimate the coefficients of ODE systems that contain time-dependent forcing terms. Since the such forcing terms can be viewed as control terms, this study shows that it is indeed possible to use the homotopy method to discover control systems and outline the possibility of using our homotopy method to train NeuralODEs on such systems as well.\n\n[1] G. A. Johnson et. al, *Phys. Rev. Lett.* **80**, 18 (1998).\n\n[2] R. Manikantan et. al, *Mathematics* **8**, 7 (2020).\n\n"}, {"review_id": "jl7kzP4AXY", "reviewer": "Reviewer_U4Km", "summary": "This paper proposes a novel method for training NeuralODEs, based on synchronization and homotopy optimization. They show that the addition of the synchronization module can smooth the loss landscape, on which homotopy optimization can be applied to enhance training. The new training method achieves competitive results in convergence speed and interpolation and extrapolation accuracy when compared with other baseline methods, especially for long training data. In addition,they demonstrate the robustness of the method experimentally.\n", "questions": "1. Is it possible to add some theoretical discussions on the smoothness of the loss landscape？\n2. Could you provide experimental results on more highly complex real-world datasets?\n", "limitations": "The authors have adequately addressed the limitations.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The method of homotopy optimization has been introduced into the training of various neural networks, but this paper combines the idea of homotopy optimization with synchronization, and proposes a method that is very suitable for training NeuralODE on dynamical system-related tasks.\n2. The new training method can effectively improve the efficiency of network training and has the potential to alleviate a series of problems caused by irregular loss landscape.\n3. Compared with other methods, this method has no restrictions on the model structure or parameters at all, so it is more general.\n4. The experiments are valid.", "weaknesses": "1. The explanations in this paper are all experimental and lack rigorous theoretical support. For example, whether it is possible to analyze if the hessian trace of loss can be strictly bounded by the traditional method after adding coupling items for homotopy optimization.\n2. The authors only test the proposed method in low-dimensional toy models.\n3. It is mentioned in the paper that there are many stabilized models that improve training by limiting the expressivity of the model, but these methods are not compared experimentally with this new method to verify the advantages of this new method.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the valuable comments and the constructive feedback. We have added the responses to both of the reviewer's questions in the global response above. We look forward to additional discussions with the reviewer.\n\nWe also add a bit more theoretical description about the synchronization below:\n\nLet's say that our system under consideration is $$d\\mathbf{x}(t)/dt=\\mathbf{F}(\\mathbf{x}(t)),$$\nwhere $\\mathbf{x} \\in \\mathbb{R}^m$ is an m-dimensional vector. We take the input vector to be transformed as a scalar variable in the form $u(t)=\\mathbf{K}^T\\mathbf{x}(t)=K_1x_1(t)+K_2x_2(t)+...+K_mx_m(t)$, with $\\mathbf{K}$ a constant column vector.\n\nThe output subsystem is then written as $$d\\mathbf{y}(t)/dt=\\mathbf{F}(\\mathbf{y}(t))-\\mathbf{D}(v(t)-u(t))$$, where $v(t)=\\mathbf{K}^T\\mathbf{y}(t)$ and $\\mathbf{D}=(D_1, D_2, ..., D_m)^T$ is an m-dimensional constant vector.\n\nObserve that if $\\mathbf{y}(t)=\\mathbf{x}(t)$ is plugged into Eq. (2), the equation is satisfied, meaning that synchronization of the systems is possible for the combined system, Eqs. (1) and (2). It can be guaranteed that if we evaluate the largest Lyapunov exponent for the subsystem Eq. (2) with respect to the trajectory $\\mathbf{y}(t)=\\mathbf{x}(t)$.\n\nConsider infinitesimal deviations of $\\mathbf{y}(t)\\$ from $\\mathbf{x}(t)$,\n$$\n\\mathbf{y}(t)=\\mathbf{x}(t)+\\delta\\mathbf{y}(t).\n$$\nFrom (2),\n$$\nd\\delta\\mathbf{y}(t)/dt=[\\partial\\mathbf{F}(\\mathbf{y})/\\partial\\mathbf{y}\\mid_{\\mathbf{y}=\\mathbf{x}}-\\mathbf{DK}^T]\\delta\\mathbf{y}.\n$$\n\nHere, if you see the $[\\partial\\mathbf{F}(\\mathbf{y})/\\partial\\mathbf{y}\\mid_{\\mathbf{y}=\\mathbf{x}}-\\mathbf{DK}^T]$ as a matrix A and for $\\delta\\mathbf{y}$ as an error vector, then the Eq.(3) can be rewritten as:\n\n\\begin{equation}\nd\\mathbf{e}/dt=\\mathbf{A}\\mathbf{e}\n\\end{equation\n\nThe solution of this differential equation is given by\n\n\\begin{equation}\n\\mathbf{e}(t)=e^{\\mathbf{A} t}\\mathbf{e}(0)\n\\end{equation}\n\nwhere $\\mathbf{e}(0)$ is the initial value of $\\mathbf{e}(t)$.\n\nWe can always make this to be converged by choosing the proper value of $\\mathbf{K}$.\n\nConnecting this result to the hessian of the loss is a bit tricky at the moment as the matrix $\\bf A$ depends on the parameters as well, and we are continuing to research on this front at the moment.\n\nWe suspect it might be the case that it is impossible to state a global smoothness claim for the loss landscape. However, we believe by drawing on the discussion of homotopy from topology, it can be possible to weakly assert the local smoothness of the loss manifold.\n\nThe main general idea in topology is to study spaces which can be continuously deformed into one another. This idea is given mathematical substance by the introduction of homeomorphisms.\n\nIf we take two topological spaces $T_1$ and $T_2$ then a map $\\alpha$ from $T_1$ to $T_2$:\n\n\\begin{equation}\n\\alpha : T_1 \\rightarrow T_2.\n\\end{equation}\n\nis called a homeomorphism if it is both continuous and has an inverse which is also continuous. This notion, homotopy, is inspired, as is the notion of homeomorphism, by the more informal or intuitive notion of deformation.\n\nWhile homeomorphism generates equivalence classes whose members are topological spaces, homotopy generates equivalence classes whose members are continuous maps. Take two continuous maps $\\alpha_1$ and $\\alpha_2$ from a space $X$ to a space $Y$:\n\n\\begin{equation}\n\\begin{split}\n\\alpha_1 : X \\rightarrow Y, \\\\\n\\alpha_2 : X \\rightarrow Y.\n\\end{split}\n\\end{equation}\n\nThen the map $\\alpha_1$ is aid to be homotopic to the map $\\alpha_2$ if $\\alpha_1$ can be deformed into $\\alpha_2$, in precise mathematical terms:\n\n\\begin{equation}\nF:X \\times [0, 1] \\rightarrow Y,\\quad F \\ continuous.\n\\end{equation}\n\nand $F$ satisfies\n\n\\begin{equation}\n\\begin{split}\nF(x, 0)=\\alpha_1(x), \\\\\\\\\nF(x, 1)=\\alpha_2(x).\n\\end{split}\n\\end{equation}\n\nIn other words as the real variable $t$ in $F(x, t)$ varies continuously from 0 to 1 in the unit interval $[0, 1]$; the map $\\alpha_1$ is deformed continuously into the map $\\alpha_2$. Homotopy is an equivalence relation and it divides the space of continuous maps from $X$ to $Y$, which is written $C(X, Y)$, into equivalence classes. Since a homeomorphism is also a continuous map, these equivalence classes are unchanged under homeomorphism of $X$ or $Y$. "}]

Author rebuttal:
1. Can our algorithm be applied to potentially partially observed high dimensional data?

We believe that our algorithm can be applied to higher-dimensional problems, provided that a couple of challenges due to the increased dimensionality are addressed. Depending on the nature of the data, we believe there are two possible strategies.
    
The first strategy is to use a latent NeuralODE model where a NeuralODE is used to model the latent state dynamics ${\bf z}(t)$ via:
$\frac{d\bf z}{dt}={\bf U_1}(t, {\bf z};\bf\theta_1)$ with the latent state mapped to the observables ${\bf u}(t)$ by ${\bf u}(t)={\bf U_2}({\bf z}(t); \bf\theta_2)$. Depending on the available prior knowledge, ${\bf U_2}$ can either be a fixed equation, a parametrized equation, or a neural network. 

Coupling between the dynamics and the data can achieved by the following augmentation: 
$\frac{d\bf z}{dt}={\bf U_1}(t, {\bf z};{\bf\theta_1})-{\bf K}({\bf u}-{\bf \hat u})$ where $\bf \hat u$ is the interpolated version of the observed data. Such type of coupling has been reported in the ODE synchronization literature, an example being [1]. We expect this scheme to perform in the NeuralODE setting as well, provided that the observable $\bf u$ contains sufficient information about the latent state $\bf z$.

The second strategy to tackle observations from high-dimensional systems is to perform a phase space reconstruction using the observable time series, then train a simple NeuralODE on the reconstructed phase space data. In this case, our training method remains the same since we have access to all of the states of the reconstructed phase space. This idea was explored in [2], where the authors used temperature time series data from a weather station in Jena, Germany to train recurrent neural networks. As temperature corresponds to a scalar observable from an unknown high dimensional dynamical system, the data was detrended, denoised, then subjected to phase space reconstruction using time delay embedding.


2. Training on real world data

We thank the reviewer for the suggestion. To address the reviewer’s concerns and to continue the discussion in our previous response, we trained a NeuralODE on a sub-sampled version of this phase space reconstructed temperature dataset (to reduce computational burden) using our homotopy-based method. We present the results in Figure 1 of the attached pdf. 

We find that our model can predict the reconstructed dynamics well in the interpolation interval, but fails in the extrapolation regime. This is likely due to insufficient portion of the high dimensional phase space being sampled by the training data, and should be improved with additional training data.

We find that our model can predict the reconstructed dynamics well in the interpolation interval, but fails in the extrapolation regime. This is likely due to insufficient portion of the high dimensional phase space being sampled by the training data, and should be improved with additional training data. Still, this result shows the feasibility of using our method on real world data, and we plan to continue exploring in this area.

3. Relationship between K and k
We apologize for the confusion. In our work, the matrix $\bf K$ and the scalar algorithm hyperparameter $k$ are related to each other by ${\bf K}=k\bf I$, where $\bf I$ is the identity matrix.

4. 

Meta-review: This is a solid submission. Novelty and relevance to the NeurIPS community have been acknowledged by all reviewers. There were multiple suggestions regarding he experimental validation of the method and we strongly encourage the authors to incorporate the feedback and especially the promised improvements in the final version of the paper.
Decision: Accept (poster)