S-SOLVER: Numerically Stable Adaptive Step Size Solver for Neural ODEs
Authors: Eliska Kloberdanz
Wei Le
Venue: Submitted to ICLR 2023 (2023)
Keywords: neural ODEs
ODE solvers
numerical integration
numerical stability

Abstract: A neural ordinary differential equation (ODE) is a relation between an unknown function and its derivatives, where the ODE is parameterized by a neural network. Therefore, to obtain a solution to a neural ODE requires a solver that performs numerical integration. Dopri5 is one of the most popular neural ODE solvers and also the default solver in torchdiffeq, a PyTorch library of ODE solvers. It is an adaptive step size solver based on the Runge-Kutta (RK) numerical methods. These methods rely on estimation of the local truncation error to select and adjust integration step size, which determines the numerical stability of the solution. A step size that is too large leads to numerical instability, while a step size that is too small may cause the solver to take unnecessarily many steps, which is computationally expensive and may even cause rounding error build up. Therefore, accurate local truncation error estimation is paramount for choosing an appropriate step size to obtain an accurate, numerically stable, and fast solution to the ODE. In this paper we propose a novel local truncation error approximation that is the first to consider solutions of four different RK orders to obtain a more reliable error estimate. This leads to a novel solver S-SOLVER (Stable Solver), which is more numerically stable; and therefore accurate. We demonstrate S-SOLVER's competitive performance in experiments on image recognition with ODE-Net, learning hamiltonian dynamics with Symplectic ODE-Net, and continuous normalizing flows (CNF).

Under review as a conference paper at ICLR 2023
S-SOLVER: NUMERICALLY STABLE ADAPTIVE STEP SIZE SOLVER FOR NEURAL ODES
Anonymous authors Paper under double-blind review
ABSTRACT
A neural ordinary differential equation (ODE) is a relation between an unknown function and its derivatives, where the ODE is parameterized by neural network. Therefore, to obtain a solution to a neural ODE requires solver that performs numerical integration. Dopris one of the most popular neural ODE solvers and also the default solver in torchdiffeq, a PyTorch library of ODE solvers. It is an adaptive step size solver based on the Runge-Kutta (RK) numerical methods. These methods rely estimation of the local truncation error select and adjust integration step size, which determines the numerical stability of the solution. A step size that is too large leads to numerical instability, while step size that is too small may cause the solver to take unnecessarily many steps, which is compu- tationally expensive and may even cause rounding error build up. Therefore, ac- curate local truncation error estimation is paramount for choosing an appropriate step size to obtain an accurate, numerically stable, and fast solution the ODE, this paper we propose novel local truncation error approximation that is the first to consider solutions of four different RK orders to obtain more reliable error estimate. This leads to novel solver S-SOLVER (Stable Solver), which is more numerically stable; and therefore accurate. We demonstrate S-SOLVER's com- petitive performance in experiments image recognition with ODE-Net, learn- ing hamiltonian dynamics with Symplectic ODE-Net, and continuous normalizing flows (CNF).
INTRODUCTION
Neural ODEs are continuous depth deep learning models that combine neural networks and ODEs.
Since their first introduction in (Chen et al. 2018). they have been used in many applications such as: stochastic differential equations (Li et al. 2020). physically informed modeling (Sanchez-Gonzal
ct al.
2019: Zhong et al. 2020). free-form continuous generative models (Grathwohl et al.) 2019] Finlay et al. 2020 mean mean-field games (Ruthotto et al. 2020). and irregularly sampled time-series
(Rubanova et al. 2019).
Neural ODEs parameterize the derivative of the hidden state using neural network; and therefore, learn non-linear mappings via differential equations. A differential equation is a relation between an unknown function and its derivatives. Ordinary differential equations describe the change of only one variable (as opposed to multiple) with respect to time, i.e.: dr/dt = f(t,z). Typically, an ODE is formulated as an initial value problem (IVP), which has the following form. Given function derivative dr/ dtt a time interval solution to the IVP yields evaluated at time b. The method for approximating z=b) is nu- merical integration; therefore, all the various ODE solvers include different methods for performing integration.
(a. b) and an initial value (e.i.: at time t = a). the
Adaptive step size solvers are amongst the most popular solvers for neural ODEs. In fact, the default solver torchdiffey (a library of ODE solvers implemented in PyTch))) is Dopri5, the Dormand- Prince 5(4) embedded adaptive step size method of the Kunge-Kutte (RK) family. Adaptive step size RK solvers perform two approximations: one of order and another of p-p1 1 and compare them to obtain the local truncation eтror, which is used to determine the integration step size. Specifically, the error is used to make a decision whether to accept or reject the solution step under the current step size and decide how to modify the step size for the next step. A step size that is too large leads

---

Under review as a conference paper at ICLR 2023
to numerical instability, while step size that is too small may cause the solver to take unnecessarily many steps, which is computationally expensive and may even cause the rounding error to build up. Therefore, accurate local estimation is paramount for choosing appropriate step size to obtain an accurate, numerically stable, and fast solution to the ODE.
The local truncation error is defined as the difference between the exact and approximate solution obtained at a given time step. All currently available adaptive step neural ODE solvers rely on estimating the local error as the difference between order and p - 1 solutions, which assumes that the order solution is exact. This is not necessarily true and the solution far from the exact one, the local error estimate inaccurate, which results in the solver making poor decisions regarding its step size.
In this paper we propose novel local truncation error estimation that takes into account multiple orders of the RK method as opposed to just order and 1 to obtain more accurate estimate of the local truncation error that guides the integration step size. Specifically, we modify the local trun- cation error estimation of Dopri8, the Dormand- Prinncc 8(7) embedded adaptive step size method. Dopri8 calculates the local truncation error as the difference between its Sth and 7th order solution. Our modification computes this error as the average of the difference between both its 8th and 7th, and also 4th and 5th order solution. This leads to a new ODE solver, S-SOLVER (Stable Solver), a modified Dopri8 integrator with more accurate local truncation enror estimation that provides more reliable information for step size calculations; and therefore, more numerically stable solution. To best knowledge, S-SOLVER is the first solver that uses multiple solution orders estimate local truncation error for adjusting its step size.
BACKGROUND
2.1 NEURAL ORDINARY DIFFERENTIAL EQUATIONS
Traditional neural networks are defined as discrete models with discrete sequence hidden layers, where the depth of the network corresponds to the number of layers. Neural ODEs (Chen et al., 2018) are continuous depth deep learning models, which parameterize the derivative the hidden state using neural network. Specifically, they ODEs that are parameterized by a neural network, which has many benefits such as memory efficiency, adaptive computation, and parameter efficiency.
Neural ODEs inspired by the dynamic systems interpretation of residual and other networks Haber et al. (2018); Weinan (2017). These networks perform sequence of transformations to hidden state:
stale state, [stale, 0z),
which can be viewed as discretized forward Euler method applied to continuous transformation. Given this interpretation, the transformation to a hidden state can formulated as an ODE:
(
= (
(2)
where state(t = 0) the input layer and state(t = T) is the output layer. Therefore, the neural ODE is IVP:
9). for to st t subject to x(to) = Feg
(3)
where f(..,)) 0) is the deep neural network, 2t, is the input, and Zt, is the output.
Neural ODEs are trainable through loss minimization, but due to their continuous nature the op- timization process is slightly different from classical discrete deep learning models. The forward pass solves the ODE with ODE solver and the backward pass computes the gradients either by backpropagating through the ODE solver or with the adjoint method (Chen et al., 2018). In this work focus the forward pass, which outputs solution the ODE.
2

---

Under review as a conference paper at ICLR 2023
2.2 NEURAL ODE SOLVERS
Solving neural ODEs that we generalized in Equation 3 requires numerical integration that can be described as follows:
z(t1) = z(to)= [.d) f(t, ()d)
(4)
This equation can be solved with an ODE solver, which returns the value ar(1) that represents the solution at the end of the time interval that satisfies the initial condition x: 2((n) Zit-
There are different types of ODE solvers that use different methods and algorithms for performing numerical integration, and the Runge- (RK) of methods that are amongst the most popular (Seiler & Seiler 1989). The basic idea behind the RK integration methods is to re-write de and di in Equation as finite steps Az and At and multiply the equations by At. which provides a change in with respect to At. The finite time step At is called the step size (Seiler & Seiler, 1989) and is typically represented h. The the simplest RK method, the Euler method, illustrates this well:
F Ptat1 Zt, + (f, Itt..
(5)
RK methods leverage the differential equation for computing the slope k of the tangent line to the function f. The slope is then used approximate f at the next time step 1++ 1. As shown (Bogacki & Shampine 1989), this can be represented as:
1t,, 11 = De, +
(6)
1=1
where
b ((tnsar.)
i-1 agght) for RK stages i 2, .....
K = f(tn + cih, + h
Σ ==1
4:
Oij
Since the Euler method approximates the slope only to proceed from / to l 1, it can be expressed using the general RK method shown in Equation 6 as:
Tt. + ((bkk1)
(
The number of times that the slope k approximated between and t++ impacts the local truncation error the RK method Burden et al. (2015); Burrage & Burrage (000). The local truncation error the difference between the exact and approximated solution determines the order of the RK method [Burden et al. (015). The order of the RK method corresponds to the order of the local truncation error minus one. For example, the local truncation error for the Euler's method O(h²), resulting in first order numerical technique.
3 NUMERICAL STABILITY OF NEURAL ODE SOLVERS
When approximating the solution of IVP. there are two primary sources of error: the roundoff error and the truncation error (Abell Braselton 2014), which impact the numerical stability the ODE solver that yields the approximate solution.
3.1 NUMERICAL STABILITY
Numerical stability can be viewed property of an algorithm, which describes the sensitivity of solution to numerical errors (Higham 2002). An unstable numerical method produces large changes

---

Under review as a conference paper at ICLR 2023
in outputs in response to small changes in inputs (Jong 1977). which lead to unexpected outputs or errors. Numerical instability arises due rounding, and truncation errors (Higham 2002). Roundoff errors are caused by approximating real numbers with finite precision, while truncation errors caused by approximating a mathematical process. Many numerical methods (e.g.: Huler's method for solving differential equations)
be derived by taking finitely many terms of Taylor series. The terms omitted constitute the truncation error, which often depends on a parameter called the step size (Higham 2002). In this paper introduce novel truncation error estimation used for setting an adaptive step size that achieves numerically stable solution.
3.2 STABILITY DIFFERENT ODE SOLVERS
There are two notions of numerical stability of ODEs: zero-stability and absolute stability. Zero- stability implies that on a fixed time interval, small perturbations data yield bounded perturbations in the solution as the step size h approaches zero (LeVeque) 2005). Absolute stability, a stronger notion of stability, guarantees the same behavior, but for a fixed step size as the time interval approaches infinity. Generally, numerical method for for solving initial value ODE is numerically stable if "small changes or perturbations in the initial conditions produce correspondingly small changes in the subsequent approximations" (Burden al. 2013).
Different ODE solvers have different numerical stability. This can be demonstrated with a canonical example of an ODE that describes a swinging pendulum:
ml df((t) (ingsiin((()))))
d1
where (()) is the angle between the pendulum and vertical axis at a time 1, l is the length of the pendulum, m is the pendulum mass, and g represents gravity. of numerical stability of three different methods that can be used to solve this pendulum IVP.
the varying degrees
Explicit
Implicit
15
Trapzzoddal
Exact
10
-5
-10
Figure 1: Comparison of numerical stability various ODE solvers
3.3 ANALYSIS ON NUMERICAL STABILITY WITH RESPECT TO ODE SOLVER STEP SIZE
To illustrate the impact of step size on the numerical stability of ODE solvers, we provide numer- ical stability analysis of the Euler method and derive its stability condition. The stability condition pertaining explicit Euler method's 's step size can be derived using the test equation:
ky, (0) a.
(10)
Applying the forward Euler method this ODE yields:
To a Fi+1 24 h(kare) (1 hk)r;
(11)

---

Under review as a conference paper at ICLR 2023
Solving for 24417
2++1 (1 hk)a; (1+ hb) xo = (1 (1+h
(12)
The exact solution is:
vtt exp(kt)
(13)
The absolute error is the absolute difference between the exact and approximated solution:
lv(4;) .ı| |eep(iik)) (1 hk) exp(hk) (1 -b))'laa)
(14)
If A: > O. the problem is unstable. If k S and [1 hk| the forward Euler method will be stable. This condition is called the stability region and pertains to the notion of absolute stability. Specifically, the analysis of the stability region is useful for determining step size that ensure absolute stability.
3.3.1 ILLUSTRATIVE EXAMPLE
We demonstrate the practical application of the theoretical numerical stability analysis shown above with an illustrative example of an ODE dy dy/dt /t = -2.3 with an initial value of y(0) 1. Figure 2a compares the exact solution 2233 with approximate solutions obtained with the explicit Euler method with varying step sizes: = 1.0, 0.7, 0.1. The solution obtained with step size = 1.0 erratic and inaccurate, while the solution with the smallest step size h 0.1 yields stable solution that is very close to the exact one. The reason for that is that kh for = 1.0 and h = 0.7 are far away from the stability region represented the blue circle in Figure 2b] Therefore,
observe that the step size h has a significant impact on the numerical stability and accuracy of the ODE solution.
Approximate, Appraximate, h=0.7
Approxmate, h=0.1
Fract
0
0
Legend h 1.0
-1
0.1
22
11
(a) The importance of step size for accurate solutions (b) The stability region of the explicit Euler method
Figure 2: The relationship between numerical stability of an ODE solution and ODE solver step size
METHOD
Prior adaptive step size solvers approximate the local truncation error as the difference between order p- and p solution, where the p order solution is assumed to be the exact solution. This means that the error estimates are not exact, but only accurate to the leading order in h, i.e.: order (Press & Teukolsky 1992). S-SOLVER an adaptive step size solver with novel, more accurate local error estimation that is used for adjusting the solver step size to achieve an accurate numerically stable solution to neural ODEs.
5-SOLVER is based on Dopri8, the Dormand-Prine 8(7) embedded adaptive step size method, which is 8th order RK method that requires 13 function evaluations per integration step (Prince & Dormand 1981) as shown in Equation [5]
13
Rt = 2t, +
(15)
5

---

Under review as a conference paper at ICLR 2023
where
h1 f(tns sta)
h; f(hn + Gghtten
h Σ aiyky) for i= =2,..,
(16)
Dij
==1
The coefficients a, b, and c in Equations 15 and 16 are defined using the Butcher tableau provided in (Prince & Dormand 1981).
In contrast Dopri8, which calculates local error as the difference between the 7th and 8th order solution, S-SOLVER uses order 8, 5, and Specifically, given a neural ODE:
x(t1) = z(to) + 1 ((t,
(17)
suppose that the solver has progressed in integration to some time step f, and approximated xtt,) as ). make further progress, the solver needs to take a step forward and compute the value of at time step f,, + h, where h is the step size. Suppose that this is approximated as it.,. + h) and that the step's error is Zerror. S-SOLVER computes Terror as an average of the difference between
8th and 7th order solution and 5th and 4th order solution to obtain more reliable estimate:
(0 + Worder
5(nn + h) ardery) + (1 + h) order, ((n + h) horder))
Ferror
(18)
2
The 5th and 4th order solution is computed using a similar process, but only with 6 stages as follows:
Rt,,11 = It., + h
(19)
1=1
where
h1 = f(n: (r.))
k; [tn + cilt, Re, +
Σ j=1
aisky) for i = 22....
(20)
G agj.
j=1
where the coefficients a, b, and с are given in the Butcher tableau provided in (Lawrence 1986).
Given a pre-defined upper bound on relative error RTOL (Ie-7 default in torchdiffeq) and upper bound on absolute error ATOL (1e-9 default in torchliff"), the solver then computes an error ratio as follows:
Terror
(21)
scale
where scale defined as:
scale ATOL + RTOL max() h)).
(22)
If r < the step is accepted, otherwise rejected and the value x at time step t, + h is approximated again with smaller step size h.

---

Under review as a conference paper at ICLR 2023
5 EXPERIMENTS
We implement S-SOLVER new solver that is part of the torchdiffeq library httpps: 7/anonymous 4open. sci sctence/r//SSOOVER-ECC B ReadMe md) and perform expen- ments image recognition with ODE-Net, learning hamiltonian dynamics with Symplectic ODE- Net, and generating new distributions with continuous normalizing flows (CNF). We demonstrate the S-SOLVER accurate and numerically stable thanks to better local error estimation that deter- mines the step size, which in turn afects the numerical stability of the ODE solution as shown in section
5.1 STIFF NEURAL ODE AND ERROR MONITORING
We first validate S-SOLVER's s numerical stability on solving the following stiff neural ODE obtained from page 353 of Burden et al. (2015):
dy/di di = 5 exp(())) - 1)2 + for S1sl, subject to y(1 = 0) = -1
(23)
This ODE equation is stiff, which means that it is likely the error due approximation is amplified and becomes dominating in the solution calculations leading to numerically unstable solution
(Burden et 2015 Kim al., 2021).
As shown in Figure
3 the neural ODE solved with S-SOLVER yields a solution that is very close to
the exact solution:
y(t) = t exp(-5t)
(24)
100
0.75
0.00
-025
-030
noural DDC solution with 5-SDIVER
-100
0.0
6.6
Figure 3: Solution a stiff Neural ODE with S-SOLVER
In addition demonstrating that S-SOLVER solve stiff neural ODEs, which typically have numerical stability issues, we also examine the local error. Figure 4 shows the local error estimate produced by Dopri5 (default solver in torchdiffeq), S-SOLVER, and also comparison of the two, which suggests that Dopri5 underestimates the local error.
30
Depi foo
C-DOUCRR loal NEWY
S-SOEVER sot BTD
2.8
-0.5
9.5
-3.0
60
3.5
Tran Raraatio
T1all Baration
Figure 4: A comparison of the local error produced by Dopris and S-SOLVER
7

---

Under review as a conference paper at ICLR 2023
5.2 IMAGE RECOGNITION
The next set of experiments focuses on image recognition with ODE-Nets. We train an ODE-Net with S-SOLVER and compare its results with ODE-Net trained with DopriS (default solver in torchdiffeq) and also a classical ResNet on two datasets: MNIST and FASHION MNIST. Table shows that the highest test accuracy on both datasets is achieved with our ODE-Net with S-SOLVER. The test accuracy on MNIST beats prior SOTA results in (Chen et al. 2018). who report 0.42% test error, i.c.: 99.58% test accuracy. Using the same experiment settings as (Chen et al., 2018), thanks to S-SOLVER we push the test accuracy to 99.73%. Our results are also better compared to, for example, (Ghosh al., 2020) who report 98.3% test accuracy that is achieved with their proposed temporal regularization.
Table 1: Results for ODE-Net with S-SOLVER on image recognition tasks
MNIST
FASHION MNIST
train acc
test acc
loss
train
test acc
loss
ODE-Net with S-SOLVER
99.99%
99.73%
0.00013
97.58%
94.00 %
0.079816
ODE-Net with dopri5 ResNet
99. 99.98%
99.69%
0.04623
97.75%
93.72%
0.055531
99.96%
99.68%
7.2E-05
98.52%
93.94%
0.065489
5.3 LEARNING HAMILTONIAN DYNAMICS
We test S-SOLVER on Symplectic ODE-Net Zhong et al. 2020). which can learn Hamiltonian dynamics. Specifically, choose the problem of "acrobot" (Murray & Hauser] 2010 Suttonn & Barto 2005). which simulates physical system with two joints and two links, where the joint between the two links is actuated. Initially, the links are hanging downwards, and the goal is to swing the end of the lower link up to given height. In Figure 5 we show that the validation loss obtained with S-SOLVER is more stable than with Dopri5 (default solver in torchdiffeq) and therefore, preferable. We interpret this observation to be the result of S-SOLVER's more reliable local estimation that controls the step size, which in turn impacts the stability of the ODE solution.
Dopri
0.015
1o3s
0.010
0.005
0.000
10
epoch
Figure 5: Acrobot: The validation loss obtained by solving Symplectic ODE-Net with S-SOLVER is more stable than with Dopri5 solver (default solver in torchdiffeg)
5.4 CONTINUOUS NORMALIZING FLOWS
Continuous Normalizing Flows (CNF) are generative models introduced by (Chen et al. 2018) that leverage neural ODEs. CNFs are based normalizing flows (Rezende & Mohamed, 2015), which perform transformations of a simple probability distribution into more complex one by a sequence of invertible and differentiable mappings (Kobyzev et al. 2021)
8

---

Under review as a conference paper at ICLR 2023
We perform experiments with CNFs that use S-SOLVER and visualize how the model generates the Two Circles distribution from random noise Figure 6 Figure @ the evolution of the generated distribution (samples) and probability density (log probability) with respect the Two Circles distribution (target) over time from time-step 0.0 to 10.0. It can be observed that by the last time step, the random distribution has been transformed into the Two Circles distribution.
Log Probability
Log Probability
Log Probability
Log Probability
Log Probabillty
Samples
Samples
Samples
Samples
Samples
Target
Target
Target
Target
Target
Time Step 0.00
Time 2.50
Time Step 5.00
Time Step 7.50
Time Step
Figure Continuous normalizing flows for fitting the Two Circles distribution with S-SOLVER
RELATED WORK
This paper focuses ODE solvers for solving neural ODEs, which involves performing numeri- cal integration. Specifically, we focus on adaptive step size ODE solvers which have become the standard for solving neural ODEs. While more numerically stable ODE solver that is based on more accurate local truncation error estimation, there are several prior works that also study numerical integration neural ODEs. Zhu et al. (2022) perform numerical analysis of numerical integration in neural ODEs and propose IMDE, or inverse modified differential equations. Zhuang al. (2021) propose MALI, new numerical integrator that is memory-efficient Ghosh et al. (2020) introduce STEER, simple temporal regularization that randomly perturbs the numerical integration time limits. Pal et al. (2021) propose regular- ization method for adaptive ODE solvers that the internal cost heuristics. Yan et al. (2020) study the robustness of the Fuler method, which the simplest, but important neural ODE solver. Krishnapriyan et al. (2022) develop convergence test that can be used to select ODE solver that is suitable for particular task.
best knowledge,
are the first ones propose
7 CONCLUSION
In this paper we demonstrate the importance of appropriately choosing and adapting the step size in ODE solvers for obtaining numerically stable; and therefore, accurate solutions neural ODEs. To this end we propose S-SOLVER, a new neural ODE solver that is more numerically sta- ble thanks to more accurate local truncation error estimation that is based on comparing multiple approximations as opposed to just two, which has been the standard approach. We provide the- oretical analysis of the impact of solver step size on numerical stability and also perform practical experiments with S-SOLVER. We show that S-SOLVER can solve stiff neural ODE and that image recognition ODE-Nets learned with 5-SOLVER surpass the test accuracy of prior solvers as well as classical ResNets on MNIST and FASHION MNIST. In fact, S-SOLVER achieves new SOTA test accuracy on MNIST. We also show that the process of learning Hamiltonian dynamics with Symplectic ODE-Nets the acrobot example is more stable with S-SOLVER than with Dopri5, the solver used in prior neural ODE works. Finally, we also show that S-SOLVER works well for CNFs in experiment, where successfully learn data distribution from random noise.

---

Under review as a conference paper at ICLR 2023
REFERENCES
Martha L. Abell and James P. Brasclton. Introductory differential equations. 2014.
Przemyslaw Bogacki and Lawrence F. Shampine. A 3(2) pair of runge kutta formulas. Applied Mathematics Letters, 2:321-325, 1989.
Richard L Burden, Douglas Faires, and Annette M Burden. Numerical analysis. Cengage learning, 2015.
Kevin Burrage and Pamela M. Burrage. Order conditions of stochastic runge-kutta methods by b-series. SIAM Numer. Anal., 38:1626-16660 1646, 2000.
Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David Kristjanson Duvenaud. Neural ordi- nary differential equations. ArXiv, abs/1806.07366 2018.
Chris Finlay, Jörn-Henrik Jacobsen, Levon Nurbekyan, and Adam M. Oberman. How to train your neural ode: the world of jacobian and kinetic regularization. In ICML, 2020.
Arna Ghosh, Harkirat Singh Behl, Emilien Dupont, Philip S. Torr, and Vinay Namboodiri. Steer Simple temporal regularization for neural odes. ArXiv, abs/2006. 10711, 2020.
Will Grathwohl, Ricky T. Q. Chen, Jesse Bettencourt, Ilya Sutskever, and David Kristjanson Duve- naud. Ffjord: Free-form continuous dynamics for scalable reversible generative models. ArXiv, abs/1810.01367, 2019.
Eldad Haber, Lars Ruthotto, and Elliot Holtham. Learning across seales a multiscale method for convolution neural networks. ArXiv, abs/1703.00009, 2018.
Nicholas John Higham. Accuracy and stability of numerical algorithms, second edition. 2002.
L. D. Jong. Towards a formal definition of numerical stability. Numerische Mathematik, 28:211- 219, 1977.
Suyong Kim, Weiqi Ji, Sili Deng, Yingbo Ma, and Christopher Rackauckas. Stiff neural ordinary differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(9):0931 122,, 2021.
Ivan Kobyzev, Simon Prince, and Marcus A. Brubaker. Normalizing flows: An introduction and review current methods. IEEE Transactions on Puttern Analysis and Machine Intelligence, 43: 3964-3979, 2021.
Aditi S. Krishnapriyan, A. F. Queiruga, N. Benjamin Erichson, and Michael W. Mahoney. Learning continuous models for continuous physics. ArXiv, abs/2202.004944, 2022.
F Shampine Lawrence. Some practical runge-kutta formulas. Mathematics of Computation, 46: 135-150, 1986.
Randall J. LeVeque. Finite difference methods for differential equations. 2005.
Xuechen Li, Ting- Kam Leonard Wong, Ricky T. Q. Chen, and David Kristjanson Duvenaud. Scal- able gradients for stochastic differential equations. ArXiv, abs/2001.00 0012000,,,38,,,, 2020.
Richard M. Murray and John Hauser. A case study in approximate linearization: The acrobot exam- ple. 2010.
Avik Pal, Yingbo Ma, Viral B. Shah, and Chris Rackauckas. Opening the blackhox: Accelerating neural differential equations by regularizing internal solver heuristics. In ICML, 2021.
William H. Press and Saul A. Teukolsky. Adaptive stepsize runge-kutta integration. Computers in Physics, 6:188-191, 1992.
P. Prince and J. R. Dormand. High order embedded runge-kutta formulae. Journal of Computa- fional and Applied Mathematics, 7:67-75, 1981.
10

---

Under review as a conference paper at ICLR 2023
Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows. In ICML, 2015.
Yulia Rubanova, Tian Qi Chen, and David Kristjanson Duvenaud. Latent ordinary differential equa- tions for irregulary-saampl time series. NeurIPS, 2019.
Lars Ruthotto, Stanley J. Osher, Wuchen Li, Levon Nurbekyan, and Samy Wu Fung. machine learning framework for solving high-dimensional mean field game and mean field control prob- lems. Proceedings of the National Academy of Sciences, 117:9183- 9193, 2020.
Alvaro Banchez-Gonzaled Victor Bapst, Kyle Crammer, and Peter W. Battaglia. Hamiltonian graph networks with ode integrators. ArXiv, abs/1909. 12790, 2019.
Mary C. Seiler and Fritz A. Seiler. Numerical recipes in The of scientific computing. Risk Analysis, 9:415-416, 1989.
Richard S. Sutton and Andrew G. Barto. Reinforcement learning: An introduction. IEEE Transac- lions on Neural Networks, 16:285-286, 2005.
Weinan. A proposal on machine learning via dynamical systems. 2017.
Hanshu Yan, Jiawei Du, Vincent Yan Fu Tan, and Jiashi Feng. On robustness of neural ordinary differential equations. ArXiv, abs/1910.05513, 2020.
Yaofeng Desmond Zhong. Biswadip Dey, and Amit Chakraborty. Symplectic ode-net: Learning hamiltonian dynamics with control. ArXiv, abs/1999..2077, 2020.
Aiqing Zhu, Pengzhan Jin, Beibei Zhu, and Yifa Tang. On numerical integration in neural ordinary differential equations. ICML, 2022.
Juntang Zhuang, Nicha C. Dvomek, Sekhar C. Tatikonda, and James S. Duncan. Mali: memory efficient and reverse accurate integrator for neural odes. ArXiv, abs/2102.04665 2021.
11


Reviews:
[{"review_id": "k_zaUlVHKAu", "reviewer": "Reviewer_cFyQ", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This study proposed S-SOLVER, which is a numerical integrator with adaptive step size based on eighth-order Dormand-Prince method. It can solve even a stiff equation.\n", "summary_of_the_review": "The proposed method appears to be just a mixture of dopri5 and dopri8, with no theoretical support that it is superior to dopri8. In addition, the proposed method is related to numerical analysis and not directly to learning. Hence, I do not feel the need to publish it from ICLR.\n", "strength_and_weaknesses": "**Strength**\n\nIt worked better than dopri5 for some tasks of neural ODEs.\n\n**Weakness**\n\nThe main concern is whether the proposed method is suited for ICLR publication. The proposal is a higher-order adaptive stepping numerical solver, which is generally applicable to any ODEs and not specialized to neural ODEs. I think the present study is about numerical analysis and is out of scope of ICLR.\n\nThere is no theoretical support that the error estimate in (8) is more accurate than that of dopri8. It appears that (8) is merely an average of the estimates by dopri8 and dopri5.\n\nS-Solver uses the methods of orders 8,7,5,4 to estimate the local truncation error. The computation cost is simply increased from dopri8 or dopri5.\n\nThe coefficients of Runge-Kutta methods are not uniquely determined, but they are not provided explicitly. They are the same as those used for dopri8 and dopri5?\n\nFor any adaptive stepping methods, one can choose the tolerance to obtain arbitrary accuracy. So, it is unclear why the proposed method worked better than dopri5. Because the proposed method is based on dopri8, it should be compared with dopri8, not dopri5.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "1: The contributions are neither significant nor novel.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "CTbHarq3mPt", "reviewer": "Reviewer_BCpH", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors aim to propose a new ODE solver for Neural ODE, specifically targeted at stiff ODEs. To achieve this goal, the authors modified the Dopri5 method, and give an estimate of the local error as the mean of error estimate of 5-th order and 8-th order. The authors also validated their method in Neural ODE experiments.", "summary_of_the_review": "The proposed method averages error estimates of different orders, which is not meaningful in theory. Furthermore, the experimental validations are too simple.", "strength_and_weaknesses": "Strengths:\nTo my knowledge, solving stiff ODEs is not extensively studied in the field of Neural ODE. The authors are discussing an important yet not well studied problem. \nFurthermore, the paper is in general well-written and easy to follow.\n\nWeakness:\nI have a few questions and hope the authors could address them, listed below.\n\n1) On the theoretical groundings of the proposed method. I'm not sure if the authors have thought about this, but Runge-Kutta methods are derived by solving a system of equations, see wikipedia for detail https://en.wikipedia.org/wiki/Runge–Kutta_methods#Adaptive_Runge–Kutta_methods. The idea is that, a $p$-th order RK method should have a error on the order of $h^p$ where $h$ is the stepsize, and the system of equations are derived to satisfy this constraint. \nUnder this consideration, the 8-th and 7-th equations gives an error estimate of order $h^8$, and the difference between 5-th and 4-th equations give an error estimate of order $h^5$. I don't think it's meaning to average two errors of different orders.\n\n2) On experiments. The experimental validations are focused on toy examples rather than some harder problems. For example, Neural ODE has been applied to ImageNet classification (1000 classes), CNF on Cifar10 and ImageNet64, which are much more complicated than the experiments on MNIST and Swiss roll here.\n\n3) (Minor, some missing related works). There are some related works that the author could consider for reference, references [1] [2] consider the numerical issue with Neural ODE solvers and have achieved SOTA results on much harder experiments.\n\n[1] Zhuang, Juntang, et al. \"Adaptive checkpoint adjoint method for gradient estimation in neural ode.\" International Conference on Machine Learning. PMLR, 2020.\n[2] Matsubara, Takashi, Yuto Miyatake, and Takaharu Yaguchi. \"Symplectic adjoint method for exact gradient of neural ode with minimal memory.\" Advances in Neural Information Processing Systems 34 (2021): 20772-20784.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "5mmBtxppG-b", "reviewer": "Reviewer_H6ea", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The step size for the ODE solvers is crucial for training Neural ODEs. A step size that is too large leads to numerical instability, while a step size that is too small may cause the solver to take unnecessarily many steps, which is computationally expensive and may even cause rounding error build-up. Therefore, accurate local truncation error estimation is paramount for choosing an appropriate step size to obtain an accurate, numerically stable, and fast solution to the ODE. In the paper under review, the authors propose a local truncation error approximation that is the first to consider solutions of four different RK orders to obtain a more reliable error estimate. This leads to a novel solver, S-SOLVER (Stable Solver), which is more numerically stable; and, therefore, accurate. The authors test S-SOLVER’s performance in experiments on image recognition with ODE-Net, learning hamiltonian dynamics with Symplectic ODE-Net, and continuous normalizing flows (CNF).", "summary_of_the_review": "I feel sorry to provide a strong rejection for this paper as my concerns are substantial about this paper. First, the authors seem confused about the concept of numerical accuracy and numerical stability. Second, the proposed s-solver is not more accurate than the existing adaptive ODE solver but is computationally much more expensive. Third, the evaluation performed in this paper is very weak, and no details is provided in the experimental section to reproduce the numerical results.", "strength_and_weaknesses": "Strength: the paper is easy to read, and the context is most about basic numerical ODE solvers.\n\nWeakness: I have many concerns about this work. In particular, the novelty, accuracy, and computational efficiency of the proposed method compared to the standard adaptive ODE solver. Below, I list some detailed comments.\n\n1. The contribution of this paper is very marginal and seems problematic to me. As stated in the paper about the contribution of this paper \"We modify the local truncation error estimation of Dopri8, the Dormand-Prince 8(7) embedded adaptive step size method. Dopri8 calculates the local truncation error as the difference between its 8th and 7th-order solution. Our modification computes this error as the average of the difference between both its 8th and 7th, and also 4th and 5th order solution.\" --- This will make the resulting ODE solver has no better convergence order but computationally much more expensive than Dopri5.  In particular, in each ODE numerical integration, we need to evaluate the right hand side of the Neural ODE many more times than the standard methods that have the same order of accuracy at best.\n\n\n2. Many existing methods, which have been developed for several decades, use multiple different-order methods for local truncation error estimation. The authors may refer to the book “Hairer, E., Nørsett, S.P., Wanner, G., Solving Ordinary Differential Equations I, Nonstiff Problems, Second Revised Edition, Springer, 2008” and “Hairer, E., Wanner, G., Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems, Second Revised Edition, Springer, 2002”. As far as I remember, DOP853 belongs to such category.\n\n\n3. The authors seem confused about numerical stability and numerical accuracy. I do not see that s-solver is more stable. Usually higher-order explicit ODE solvers are less stable than the low-order explicit methods. Implicit ODE solvers are much more stable than explicit ODE solvers.\n\n\n4. The ODE in section 5.1 is not a stiff ODE. The dynamics of a stiff ODE should change abruptly. The following tutorial provides you an example of stiff ODE\nhttps://www.mathworks.com/company/newsletters/articles/stiff-differential-equations.html\n\nAgain, this example shows that the authors seem to have some misunderstanding of stiff ODEs.\n\n\n5. Figure 3 is meaningless; you should contrast s-solver with other ODE solvers. Also, y-label should be y(t) rather than f(t). Why the error between the numerical solution and the exact one is so large when t is around 0.2? I also want to see the behavior of numerical and exact solutions when t is much larger than 1. For instance, you may plot for t from 0 to 10.\n\n6. All currently available adaptive step neural ODE solvers rely on estimating the local error as the difference between order $p$ and $p − 1$ solutions, which assumes that the order $p$ solution is exact. This is not necessarily true, and if the $p$ solution is far from the exact one, the local error estimate is inaccurate, which results in the solver making poor decisions regarding its step size. - - - This statement seems wrong. \n\n\n7. The authors claim s-solver is stable; why is this the case? I would like to see some discuss of this. For instance, the authors may follow some existing discussions on the stability of ODE solvers. E.g., https://math.stackexchange.com/questions/2282814/stability-of-numerical-ode-solvers\n\n\n8. There are many issues with the experiments: 1) There is no detail of the experimental settings, e.g., hyperparameters from training the machine learning model, neural network architecture, etc. 2) The testing examples are rather toyed examples. 3) The authors did not report the computational time.\n\n\n9. The authors may test the impact of the error of dopri5 and s-solver on the performance of the trained machine learning model.\n\n\n\nMinor comments:\n\n1. The writing style is really not appreciated. The entire section 2, which takes slightly less than two pages, is reviewing existing neural ODE. The entire section 3, which takes around two pages, is simply adapted from some standard numerical ODE textbook. The paper contains really minimal new knowledge.\n\n2. There are some typos. E.g., in the abstract, \"build up\" should be replaced with \"build-up\".\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "1: The main claims of the paper are incorrect or not at all supported by theory or empirical results.", "technical_novelty_and_significance": "1: The contributions are neither significant nor novel.", "empirical_novelty_and_significance": "1: The contributions are neither significant nor novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject