summaries: This paper presents a memory-efficient asynchronous leapfrog integrator for numerically solving neural ODEs, referred to as MALI. The method comes with a constant memory guarantee (like the adjoint method) and also guarantees reverse-time accuracy (like the adaptive checkpoint adjoint (ACA) method). The authors also give a rigorous theoretical analysis of MALI, and also discuss a "damped" version with an increased stability region. The method is evaluated on a variety of tasks which includes classification, dynamical modelling and generative modelling. There are typically two methods for estimating the gradients with respect to the loss for neural ODEs. The naive method directly backpropagates through the steps of the ODE solver leading to accurate gradients but very large memory cost. The adjoint method in contrast does not store the entire trajectory in memory, but has reverse trajectory errors (i.e. the numerical solution in the reverse direction will not be the inverse of the numerical solution in the forward direction). In this paper, the authors propose a method that is both reverse accurate and has low memory cost. This paper proposes a new algorithm for solving neural ODEs. Each numerical solver step of the neural ODE is implemented as an invertible neural network via a variant of the asynchronous leafprog integrator. While still computing an accurate gradient, this allows memory savings by discarding intermediate data from the numerical integration steps since it can be reconstructed using the inverse. A theoretical stability analysis is provided. The experimental results show that the algorithm achieves similar performance to previous methods (e.g. ACA) while using less memory. 
abstractText: Neural ordinary differential equations (Neural ODEs) are a new family of deeplearning models with continuous depth. However, the numerical estimation of the gradient in the continuous case is not well solved: existing implementations of the adjoint method suffer from inaccuracy in reverse-time trajectory, while the naive method and the adaptive checkpoint adjoint method (ACA) have a memory cost that grows with integration time. In this project, based on the asynchronous leapfrog (ALF) solver, we propose the Memory-efficient ALF Integrator (MALI), which has a constant memory cost w.r.t number of solver steps in integration similar to the adjoint method, and guarantees accuracy in reverse-time trajectory (hence accuracy in gradient estimation). We validate MALI in various tasks: on image recognition tasks, to our knowledge, MALI is the first to enable feasible training of a Neural ODE on ImageNet and outperform a well-tuned ResNet, while existing methods fail due to either heavy memory burden or inaccuracy; for time series modeling, MALI significantly outperforms the adjoint method; and for continuous generative models, MALI achieves new state-of-theart performance.We provide a pypi package: https://jzkay12.github. io/TorchDiffEqPack
authors.name: Juntang Zhuang
authors.name: Nicha C. Dvornek
authors.name: Sekhar Tatikonda
authors.name: James S. Duncan
references.authors: Jens Behrmann Will Grathwohl Ricky TQ Chen David Duvenaud Jörn-Henrik Jacobsen
references.title: Invertible residual networks
references.venue: In International Conference on Machine Learning,
references.authors: Ricky TQ Chen Yulia Rubanova Jesse Bettencourt David K Duvenaud
references.title: Neural ordinary differential equations
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Ricky TQ Chen Jens Behrmann David K Duvenaud Jörn-Henrik Jacobsen
references.title: Residual flows for invertible generative modeling
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Talgat Daulbaev Alexandr Katrutsa Larisa Markeeva Julia Gusak Andrzej Cichocki Ivan Oseledets
references.title: Interpolated adjoint method for neural odes
references.venue: arXiv preprint arXiv:2003.05271,
references.authors: Alfredo M Ozorio De Almeida
references.title: Hamiltonian systems: chaos and quantization
references.authors: Laurent Dinh Jascha Sohl-Dickstein Samy Bengio
references.title: Density estimation using real nvp
references.venue: arXiv preprint arXiv:1605.08803,
references.authors: Emilien Dupont Arnaud Doucet Yee Whye Teh
references.title: Augmented neural odes
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Chris Finlay Jörn-Henrik Jacobsen Levon Nurbekyan Adam M Oberman
references.title: How to train your neural ode: the world of jacobian and kinetic regularization
references.venue: In International Conference on Machine Learning,
references.authors: Ian J Goodfellow Jonathon Shlens Christian Szegedy
references.title: Explaining and harnessing adversarial examples
references.venue: arXiv preprint arXiv:1412.6572,
references.authors: Will Grathwohl Ricky TQ Chen Jesse Bettencourt Ilya Sutskever David Duvenaud
references.title: Ffjord: Free-form continuous dynamics for scalable reversible generative models
references.venue: arXiv preprint arXiv:1810.01367,
references.authors: Eldad Haber Lars Ruthotto
references.title: Stable architectures for deep neural networks
references.venue: Inverse Problems,
references.authors: YAN Hanshu DU Jiawei TAN Vincent FENG Jiashi
references.title: On robustness of neural ordinary differential equations
references.venue: In International Conference on Learning Representations,
references.authors: Kaiming He Xiangyu Zhang Shaoqing Ren Jian Sun
references.title: Deep residual learning for image recognition
references.venue: In Proceedings of the IEEE conference on computer vision and pattern recognition,
references.authors: Jonathan Ho Xi Chen Aravind Srinivas Yan Duan Pieter Abbeel
references.title: Flow++: Improving flowbased generative models with variational dequantization and architecture design
references.authors: Junteng Jia Austin R Benson
references.title: Neural jump stochastic differential equations
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Patrick Kidger Ricky T.Q. Chen Terry Lyons
references.title: Hey, that’s not an ODE”: Faster ODE Adjoints with 12 Lines of Code
references.authors: Patrick Kidger James Morrill James Foster Terry Lyons
references.title: Neural controlled differential equations for irregular time series
references.venue: arXiv preprint arXiv:2005.08926,
references.authors: Durk P Kingma Prafulla Dhariwal
references.title: Glow: Generative flow with invertible 1x1 convolutions
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Xuechen Li Ting-Kam Leonard Wong Ricky TQ Chen David Duvenaud
references.title: Scalable gradients for stochastic differential equations
references.venue: arXiv preprint arXiv:2001.01328,
references.authors: Kuang Liu
references.title: Train cifar10 with pytorch. 2017
references.venue: URL https://github.com/kuangliu/ pytorch-cifar
references.authors: Yiping Lu Aoxiao Zhong Quanzheng Li Bin Dong
references.title: Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations
references.venue: In International Conference on Machine Learning,
references.authors: Stefano Massaroli Michael Poli Jinkyoo Park Atsushi Yamashita Hajime Asama
references.title: Dissecting neural odes
references.venue: arXiv preprint arXiv:2002.08071,
references.authors: Ulrich Mutze
references.title: An asynchronous leapfrog method ii
references.venue: arXiv preprint arXiv:1311.6602,
references.authors: Razvan Pascanu Tomas Mikolov Yoshua Bengio
references.title: On the difficulty of training recurrent neural networks
references.venue: In International conference on machine learning,
references.authors: Lev Semenovich Pontryagin
references.title: Mathematical theory of optimal processes
references.authors: Alessio Quaglino Marco Gallieri Jonathan Masci Jan
references.title: Koutnı́k. Snode: Spectral discretization of neural odes for system identification
references.authors: Alejandro F Queiruga N Benjamin Erichson Dane Taylor Michael W Mahoney
references.title: Continuousin-depth neural networks
references.venue: arXiv preprint arXiv:2008.02389,
references.authors: Yulia Rubanova Ricky TQ Chen David K Duvenaud
references.title: Latent ordinary differential equations for irregularly-sampled time series
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Lars Ruthotto Eldad Haber
references.title: Deep neural networks motivated by partial differential equations
references.venue: Journal of Mathematical Imaging and Vision,
references.authors: Lars Ruthotto Stanley J Osher Wuchen Li Levon Nurbekyan Samy Wu Fung
references.title: A machine learning framework for solving high-dimensional mean field game and mean field control problems
references.venue: Proceedings of the National Academy of Sciences,
references.authors: Alvaro Sanchez-Gonzalez Victor Bapst Kyle Cranmer Peter Battaglia
references.title: Hamiltonian graph networks with ode integrators
references.venue: arXiv preprint arXiv:1909.12790,
references.authors: John R Silvester
references.title: Determinants of block matrices
references.venue: The Mathematical Gazette,
references.authors: Endre Süli David F Mayers
references.title: An introduction to numerical analysis
references.venue: Cambridge university press,
references.authors: Paulo Tabuada Bahman Gharesifard
references.title: Universal approximation power of deep neural networks via nonlinear control theory
references.venue: arXiv preprint arXiv:2007.06007,
references.authors: Loup Verlet
references.title: Computer” experiments” on classical fluids. i. Thermodynamical properties of Lennard-Jones molecules
references.venue: Physical review,
references.authors: Gerhard Wanner Ernst Hairer
references.title: Solving ordinary differential equations II
references.authors: Antoine Wehenkel Gilles Louppe
references.title: Unconstrained monotonic neural networks
references.venue: In Advances in Neural Information Processing Systems,
references.authors: E Weinan
references.title: A proposal on machine learning via dynamical systems
references.venue: Communications in Mathematics and Statistics,
references.authors: Haruo Yoshida
references.title: Construction of higher order symplectic integrators
references.venue: Physics letters A,
references.authors: Yaofeng Desmond Zhong Biswadip Dey Amit Chakraborty
references.title: Symplectic ode-net: Learning hamiltonian dynamics with control
references.venue: arXiv preprint arXiv:1909.12077,
references.authors: Juntang Zhuang Nicha Dvornek Xiaoxiao Li Sekhar Tatikonda Xenophon Papademetris James Duncan
references.title: Adaptive checkpoint adjoint method for gradient estimation in neural ode
references.venue: International Conference on Machine Learning,
references.authors: Zhuang
references.title: 2020) and an open-source repository
references.authors: However Queiruga
references.title: 2020) argues that many numerical discretizations fail to be meaningful dynamical systems, while our experiments demonstrate that our model is continuous hence invariant to discretization schemes. Adversarial robustness Besides the high accuracy and robustness to discretization scheme, another advantage of Neural ODE is the robustness to adversarial attack. The adversary robustness
sections.text: Neural ordinary differential equations (Neural ODEs) are a new family of deeplearning models with continuous depth. However, the numerical estimation of the gradient in the continuous case is not well solved: existing implementations of the adjoint method suffer from inaccuracy in reverse-time trajectory, while the naive method and the adaptive checkpoint adjoint method (ACA) have a memory cost that grows with integration time. In this project, based on the asynchronous leapfrog (ALF) solver, we propose the Memory-efficient ALF Integrator (MALI), which has a constant memory cost w.r.t number of solver steps in integration similar to the adjoint method, and guarantees accuracy in reverse-time trajectory (hence accuracy in gradient estimation). We validate MALI in various tasks: on image recognition tasks, to our knowledge, MALI is the first to enable feasible training of a Neural ODE on ImageNet and outperform a well-tuned ResNet, while existing methods fail due to either heavy memory burden or inaccuracy; for time series modeling, MALI significantly outperforms the adjoint method; and for continuous generative models, MALI achieves new state-of-theart performance.We provide a pypi package: https://jzkay12.github. io/TorchDiffEqPack
sections.heading: 1 INTRODUCTION
sections.text: Recent research builds the connection between continuous models and neural networks. The theory of dynamical systems has been applied to analyze the properties of neural networks or guide the design of networks (Weinan, 2017; Ruthotto & Haber, 2019; Lu et al., 2018). In these works, a residual block (He et al., 2016) is typically viewed as a one-step Euler discretization of an ODE; instead of directly analyzing the discretized neural network, it might be easier to analyze the ODE.
  Another direction is the neural ordinary differential equation (Neural ODE) (Chen et al., 2018), which takes a continuous depth instead of discretized depth. The dynamics of a Neural ODE is typically approximated by numerical integration with adaptive ODE solvers. Neural ODEs have been applied in irregularly sampled time-series (Rubanova et al., 2019), free-form continuous generative models (Grathwohl et al., 2018; Finlay et al., 2020), mean-field games (Ruthotto et al., 2020), stochastic differential equations (Li et al., 2020) and physically informed modeling (SanchezGonzalez et al., 2019; Zhong et al., 2019).
  Though the Neural ODE has been widely applied in practice, how to train it is not extensively studied. The naive method directly backpropagates through an ODE solver, but tracking a continuous trajectory requires a huge memory. Chen et al. (2018) proposed to use the adjoint method to determine the gradient in continuous cases, which achieves constant memory cost w.r.t integration time; however, as pointed out by Zhuang et al. (2020), the adjoint method suffers from numerical errors due to the inaccuracy in reverse-time trajectory. Zhuang et al. (2020) proposed the adaptive checkpoint adjoint (ACA) method to achieve accuracy in gradient estimation at a much smaller memory cost compared to the naive method, yet the memory consumption of ACA still grows linearly with integration time. Due to the non-constant memory cost, neither ACA nor naive method are suitable for large scale datasets (e.g. ImageNet) or high-dimensional Neural ODEs (e.g. FFJORD (Grathwohl et al., 2018)).
  In this project, we propose the Memory-efficient Asynchronous Leapfrog Integrator (MALI) to achieve advantages of both the adjoint method and ACA: constant memory cost w.r.t integration time and accuracy in reverse-time trajectory. MALI is based on the asynchronous leapfrog (ALF) integrator (Mutze, 2013). With the ALF integrator, each numerical step forward in time is reversible. Therefore, with MALI, we delete the trajectory and only keep the end-time states, hence achieve constant memory cost w.r.t integration time; using the reversibility, we can accurately reconstruct the trajectory from the end-time value, hence achieve accuracy in gradient. Our contributions are:
  1. We propose a new method (MALI) to solve Neural ODEs, which achieves constant memory cost w.r.t number of solver steps in integration and accuracy in gradient estimation. We provide theoretical analysis.
  2. We validate our method with extensive experiments: (a) for image classification tasks, MALI enables a Neural ODE to achieve better accuracy than a well-tuned ResNet with the same number of parameters; to our knowledge, MALI is the first method to enable training of Neural ODEs on a large-scale dataset such as ImageNet, while existing methods fail due to either heavy memory burden or inaccuracy. (b) In time-series modeling, MALI achieves comparable or better results than other methods. (c) For generative modeling, a FFJORD model trained with MALI achieves new state-of-the-art results on MNIST and Cifar10.
sections.heading: 2 PRELIMINARIES
sections.heading: 2.1 NUMERICAL INTEGRATION METHODS
sections.text: An ordinary differential equation (ODE) typically takes the form dz(t)
  dt = fθ(t, z(t)) s.t. z(t0) = x, t ∈ [t0, T ], Loss = L(z(T ), y) (1)
  where z(t) is the hidden state evolving with time, T is the end time, t0 is the start time (typically 0), x is the initial state. The derivative of z(t) w.r.t t is defined by a function f , and f is defined as a sequence of layers parameterized by θ. The loss function is L(z(T ), y), where y is the target variable. Eq. 1 is called the initial value problem (IVP) because only z(t0) is specified.
  Algorithm 1: Numerical Integration Input initial state x, start time t0, end time T , error tolerance etol, initial stepsize h. Initialize z(0) = x, t = t0 While t < T
  error est =∞ While error est > etol
  h← h×DecayFactor ẑ, error est = ψh(t, z)
  If error est < etol h← h× IncreaseFactor t← t+ h, z ← ẑ
  Notations We summarize the notations following Zhuang et al. (2020).
  • zi(ti)/z(τi): hidden state in forward/reverse time trajectory at time ti/τi. • ψh(ti, zi): the numerical solution at time ti + h, starting from (ti, zi) with a stepsize h. • Nf , Nz: Nf is the number of layers in f in Eq. 1, Nz is the dimension of z. • Nt/Nr: number of discretized points (outer iterations in Algo. 1) in forward / reverse integration.
  • m: average number of inner iterations in Algo. 1 to find an acceptable stepsize.
  Numerical Integration The algorithm for general adaptive-stepsize numerical ODE solvers is summarized in Algo. 1 (Wanner & Hairer, 1996). The solver repeatedly advances in time by a step, which is the outer loop in Algo. 1 (blue curve in Fig. 1). For each step, the solver decreases the stepsize until the estimate of error is lower than the tolerance, which is the inner loop in Algo. 1 (green curve in Fig. 1). For fixed-stepsize solvers, the inner loop is replaced with a single evaluation of ψh(t, z) using predefined stepsize h. Different methods typically use different ψ, for example different orders of the Runge-Kutta method (Runge, 1895).
sections.heading: 2.2 ANALYTICAL FORM OF GRADIENT IN CONTINUOUS CASE
sections.text: We first briefly introduce the analytical form of the gradient in the continuous case, then we compare different numerical implementations in the literature to estimate the gradient. The analytical form
  Figure 2: In backward-pass, the adjoint method reconstructs trajectory as a separate IVP. Naive, ACA and MALI track the forward-time trajectory, hence are accurate. ACA and MALI only backpropagate through the accepted step, while naive method backpropagates through the search process hence has deeper computation graphs.
  of the gradient in the continuous case is
  dL dθ = − ∫ 0 T a(t)> ∂f(z(t), t, θ) ∂θ dt (2)
  da(t) dt + (∂f(z(t), t, θ) ∂z(t) )> a(t) = 0 ∀t ∈ (0, T ), a(T ) = ∂L ∂z(T ) (3)
  where a(t) is the “adjoint state”. Detailed proof is given in (Pontryagin, 1962). In the next section we compare different numerical implementations of this analytical form.
sections.heading: 2.3 NUMERICAL IMPLEMENTATIONS IN THE LITERATURE FOR THE ANALYTICAL FORM
sections.text: We compare different numerical implementations of the analytical form in this section. The forwardpass and backward-pass of different methods are demonstrated in Fig. 1 and Fig. 2 respectively. Forward-pass is similar for different methods. The comparison of backward-pass among different methods are summarized in Table. 1. We explain methods in the literature below.
  Naive method The naive method saves all of the computation graph (including search for optimal stepsize, green curve in Fig. 2) in memory, and backpropagates through it. Hence the memory cost is NzNf ×Nt×m and depth of computation graph are Nf ×Nt×m, and the computation is doubled considering both forward and backward passes. Besides the large memory and computation, the deep computation graph might cause vanishing or exploding gradient (Pascanu et al., 2013).
  Adjoint method Note that we use “adjoint state equation” to refer to the analytical form in Eq. 2 and 3, while we use “adjoint method” to refer to the numerical implementation by Chen et al. (2018). As in Fig. 1 and 2, the adjoint method forgets forward-time trajectory (blue curve) to achieve memory cost NzNf which is constant to integration time; it takes the end-time state (derived from forward-time integration) as the initial state, and solves a separate IVP (red curve) in reverse-time. Theorem 2.1. (Zhuang et al., 2020) For an ODE solver of order p, the error of the reconstructed initial value by the adjoint method is ∑N−1 k=0 [ hp+1k DΦ T tk (zk)l(tk, zk) +
  (−hk)p+1DΦtkT (zk)l(tk, zk) ]
  + O(hp+1), where Φ is the ideal solution, DΦ is the Jacobian of Φ, l(t, z) and l(t, z) are the local error in forward-time and reverse-time integration respectively.
  Theorem 2.1 is stated as Theorem 3.2 in Zhuang et al. (2020); please see reference paper for detailed proof. To summarize, due to inevitable errors with numerical ODE solvers, the reverse-time trajectory (red curve, z(τ)) cannot match the forward-time trajectory (blue curve, z(t)) accurately. The error in z propagates to dLdθ by Eq. 2, hence affects the accuracy in gradient estimation.
  Adaptive checkpoint adjoint (ACA) To solve the inaccuracy of adjoint method, Zhuang et al. (2020) proposed ACA: ACA stores forward-time trajectory in memory for backward-pass, hence guarantees accuracy; ACA deletes the search process (green curve in Fig. 2), and only backpropagates through the accepted step (blue curve in Fig. 2), hence has a shallower computation graph (Nf × Nt for ACA vs Nf × Nt ×m for naive method). ACA only stores {z(ti)}Nti=1, and deletes the computation graph for {f ( z(ti), ti ) }Nti=1, hence the memory cost is Nz(Nf + Nt). Though the memory cost is much smaller than the naive method, it grows linearly with Nt, and can not handle very high dimensional models. In the following sections, we propose a method to overcome all these disadvantages of existing methods.
sections.heading: 3 METHODS
sections.heading: 3.1 ASYNCHRONOUS LEAPFROG INTEGRATOR
sections.text: In this section we give a brief introduction to the asynchronous leapfrog (ALF) method (Mutze, 2013), and we provide theoretical analysis which is missing in Mutze (2013). For general firstorder ODEs in the form of Eq. 1, the tuple (z, t) is sufficient for most ODE solvers to take a step numerically. For ALF, the required tuple is (z, v, t), where v is the “approximated derivative”. Most numerical ODE solvers such as the Runge-Kutta method (Runge, 1895) track state z evolving with time, while ALF tracks the “augmented state” (z, v). We explain the details of ALF as below.
  Algorithm 2: Forward of ψ in ALF Input (zin, vin, sin, h) where sin is current
  time, zin and vin are correponding values at time sin, h is stepsize.
  Forward s1 = sin + h/2 k1 = zin + vin × h/2 u1 = f(k1, s1) vout = vin + 2(u1 − vin) zout = k1 + vout × h/2 sout = s1 + h/2 Output (zout, vout, sout, h)
  Algorithm 3: ψ−1 (Inverse of ψ) in ALF Input (zout, vout, sout, h) where sout is current time, zout and vout are corresponding values at sout, h is stepsize. Inverse s1 = sout − h/2 k1 = zout − vout × h/2 u1 = f(k1, s1) vin = 2u1 − vout zin = k1 − vin × h/2 sin = s1 − h/2 Output (zin, vin, sin, h)
  Figure 3: With ALF method, given any tuple (zj , vj , tj) and discretized time points {ti}Nti=1, we can reconstruct the entire trajectory accurately due to the reversibility of ALF.
  Procedure of ALF Different ODE solvers have different ψ in Algo. 1, hence we only summarize ψ for ALF in Algo. 2. Note that for a complete algorithm of integration for ALF, we need to plug Algo. 2 into Algo. 1. The forward-pass is summarized in Algo. 2. Given stepsize h, with input (zin, vin, sin), a single step of ALF outputs (zout, vout, sout).
  As in Fig. 3, given (z0, v0, t0), the numerical forwardtime integration calls Algo. 2 iteratively:
  (zi, vi, ti, hi) = ψ(zi−1, vi−1, ti−1, hi)
  s.t. hi = ti − ti−1, i = 1, 2, ...Nt (4)
  Invertibility of ALF An interesting property of ALF is thatψ defines a bijective mapping; therefore, we can reconstruct (zin, vin, sin, h) from (zout, vout, sout, h), as demonstrated in Algo. 7. As in Fig. 3, we can reconstruct the entire trajectory given the state (zj , vj) at time tj , and the discretized time points {t0, ...tNt}. For example, given (zNt , vNt) and {ti} Nt i=0, the trajectory for Eq. 4 is reconstructed:
  (zi−1, vi−1, ti−1, hi) = ψ −1(zi, vi, ti, hi) s.t. hi = ti − ti−1, i = Nt, Nt − 1, ..., 1 (5)
  In the following sections, we will show the invertibility of ALF is the key to maintain accuracy at a constant memory cost to train Neural ODEs. Note that “inverse” refers to reconstructing the input from the output without computing the gradient, hence is different from “back-propagation”.
  Initial value For an initial value problem (IVP) such as Eq. 1, typically z0 = z(t0) is given while v0 is undetermined. We can construct v0 = f(z(t0), t0), so the initial augmented state is (z0, v0).
  Difference from midpoint integrator The midpoint integrator (Süli & Mayers, 2003) is similar to Algo. 2, except that it recomputes vin = f(zin, sin) for every step, while ALF directly uses the input vin. Therefore, the midpoint method does not have an explicit form of inverse.
  Local truncation error Theorem 3.1 indicates that the local truncation error of ALF is of order O(h3); this implies the global error is O(h2). Detailed proof is in Appendix A.3. Theorem 3.1. For a single step in ALF with stepsize h, the local truncation error of z isO(h3), and the local truncation error of v is O(h2).
  A-Stability The ALF solver has a limited stability region, but this can be solved with damping. The damped ALF replaces the update of vout in Algo. 2 with vout = vin + 2η(u1 − vin), where η is the “damping coefficient” between 0 and 1. We have the following theorem on its numerical stability. Theorem 3.2. For the damped ALF integrator with stepsize h, where σi is the i-th eigenvalue of the Jacobian ∂f∂z , then the solver is A-stable if
  ∣∣∣1 + η(hσi − 1)±√η[2hσi + η(hσi − 1)2]∣∣∣ < 1, ∀i Proof is in Appendix A.4 and A.5. Theorem 3.2 implies the following: when η = 1, the damped ALF reduces to ALF, and the stability region is empty; when 0 < η < 1, the stability region is nonempty. However, stability describes the behaviour when T goes to infinity; in practice we always use a bounded T and ALF performs well. Inverse of damped ALF is in Appendix A.5.
sections.heading: 3.2 MEMORY-EFFICIENT ALF INTEGRATOR (MALI) FOR GRADIENT ESTIMATION
sections.text: An ideal solver for Neural ODEs should achieve two goals: accuracy in gradient estimation and constant memory cost w.r.t integration time. Yet none of the existing methods can achieve both goals. We propose a method based on the ALF solver, which to our knowledge is the first method to achieve the two goals simultaneously.
  Algorithm 4: MALI to acheive accuracy at a constant memory cost w.r.t integration time Input Initial state z0, start time t0, end time T Forward
  Apply the numerical integration in Algo. 1, with the ψ function defined by Algo. 2. Delete computation graph on the fly, only keep end-time state (zNt , vNt) Keep accepted discretized time points {ti}Nti=0 (ignore process to search for optimal stepsize)
  Backward Initialize a(T ) = ∂L∂z(T ) by Eq. 3, initialize dL dθ = 0
  For i in {Nt, Nt − 1, ..., 2, 1}: Reconstruct (zi−1, vi−1) from (zi, vi) by Algo. 7 Local forward (zi, vi, ti, hi) = ψ(zi−1, vi−1, ti−1, hi) Local backward, get ∂f(zi−1,ti−1,θ)∂zi−1 and ∂f(zi−1,ti−1,θ) ∂θ
  Update a(t) and dLdθ by Eq. 2 and Eq. 3 discretized at time points ti−1 and ti Delete local computation graph
  Output the adjoint state a(t0) (gradient w.r.t input z0) and parameter gradient dLdθ
  Procedure of MALI Details of MALI are summarized in Algo. 4. For the forward-pass, we only keep the end-time state (zNt , vNt) and the accepted discretized time points (blue curves in Fig. 1 and 2). We ignore the search process for optimal stepsize (green curve in Fig. 1 and 2), and delete other variables to save memory. During the backward pass, we can reconstruct the forward-time trajectory as in Eq. 5, then calculate the gradient by numerical discretization of Eq. 2 and Eq. 3.
  Constant memory cost w.r.t number of solver steps in integration We delete the computation graph and only keep the end-time state to save memory. The memory cost is Nz(Nf + 1), where
  NzNf is due to evaluating f(z, t) and is irreducible for all methods. Compared with the adjoint method, MALI only requires extra Nz memory to record vNt , and also has a constant memory cost w.r.t time step Nt. The memory cost is Nz(Nf + 1).
  Accuracy Our method guarantees the accuracy of reverse-time trajectory (e.g. blue curve in Fig. 2 matches the blue curve in Fig. 1), because ALF is explicitly invertible for free-form f (see Algo. 7). Therefore, the gradient estimation in MALI is more accurate compared to the adjoint method.
  Computation cost Recall that on average it takes m steps to find an acceptable stepsize, whose error estimate is below tolerance. Therefore, the forward-pass with search process has computation burden Nz ×Nf ×Nt ×m. Note that we only reconstruct and backprop through the accepted step and ignore the search process, hence it takes another Nz ×Nf ×Nt × 2 computation. The overall computation burden is NzNf ×Nt × (m+ 2) as in Table 1. Shallow computation graph Similar to ACA, MALI only backpropagates through the accepted step (blue curve in Fig. 2) and ignores the search process (green curve in Fig. 2), hence the depth of computation graph is Nf ×Nt. The computation graph of MALI is much shallower than the naive method, hence is more robust to vanishing and exploding gradients (Pascanu et al., 2013).
  Summary The adjoint method suffers from inaccuracy in reverse-time trajectory, the naive method suffers from exploding or vanishing gradient caused by deep computation graph, and ACA finds a balance but the memory grows linearly with integration time. MALI achieves accuracy in reversetime trajectory, constant memory w.r.t integration time, and a shallow computation graph.
sections.heading: 4 EXPERIMENTS
sections.heading: 4.1 VALIDATION ON A TOY EXAMPLE
sections.text: We compare the performance of different methods on a toy example, defined as
  L(z(T )) = z(T )2 s.t. z(0) = z0, dz(t)/dt = αz(t) (6)
  The analytical solution is
  z(t) = z0e αt, L = z20e 2αT , dL/dz0 = 2z0e 2αT , dL/dα = 2Tz20e 2αT (7)
  We plot the amplitude of error between numerical solution and analytical solution varying with T (integrated under the same error tolerance, rtol = 10−5, atol = 10−6) in Fig 4. ACA and MALI have similar errors, both outperforming other methods. We also plot the memory consumption for
  different methods on a Neural ODE with the same input in Fig. 4. As the error tolerance decreases, the solver evaluates more steps, hence the naive method and ACA increase memory consumption, while MALI and the adjoint method have a constant memory cost. These results validate our analysis in Sec. 3.2 and Table 1, and shows MALI achieves accuracy at a constant memory cost.
sections.heading: 4.2 IMAGE RECOGNITION WITH NEURAL ODE
sections.text: We validate MALI on image recognition tasks using Cifar10 and ImageNet datasets. Similar to Zhuang et al. (2020), we modify a ResNet18 into its corresponding Neural ODE: the forward function is y = x+ fθ(x) and y = x+ ∫ T 0 fθ(z)dt for the residual block and Neural ODE respectively, where the same fθ is shared. We compare MALI with the naive method, adjoint method and ACA.
  Results on Cifar10 Results of 5 independent runs on Cifar10 are summarized in Fig. 5. MALI achieves comparable accuracy to ACA, and both significantly outperform the naive and the adjoint method. Furthermore, the training speed of MALI is similar to ACA, and both are almost two times faster than the adjoint memthod, and three times faster than the naive method. This validates our analysis on accuracy and computation burden in Table 1.
  Accuracy on ImageNet Due to the heavy memory burden caused by large images, the naive method and ACA are unable to train a Neural ODE on ImageNet with 4 GPUs; only MALI and the adjoint method are feasible due to the constant memory. We also compare the Neural ODE to a standard ResNet. As shown in Fig. 6, the accuracy of the Neural ODE trained with MALI closely follows ResNet, and significantly outperforms the adjoint method (top-1 validation: 70% v.s. 63%).
  Invariance to discretization scheme A continuous model should be invariant to discretization schemes (e.g. different types of ODE solvers) as long as the discretization is sufficiently accurate. We test the Neural ODE using different solvers without re-training; since ResNet is often viewed as
  a one-step Euler discretization of an ODE (Haber & Ruthotto, 2017), we perform similar experiments. As shown in Table 2, Neural ODE consistently achieves high accuracy (∼70%), while ResNet drops to random guessing (∼0.1%) because ResNet as a one-step Euler discretization fails to be a meaningful dynamical system (Queiruga et al., 2020).
  Robustness to adversarial attack Hanshu et al. (2019) demonstrated that Neural ODE is more robust to adversarial attack than ResNet on small-scale datasets such as Cifar10. We validate this result on the large-scale ImageNet dataset. The top-1 accuracy of Neural ODE and ResNet under FGSM attack (Goodfellow et al., 2014) are summarized in Table 3. For Neural ODE, due to its invariance to discretization scheme, we derive the gradient for attack using a certain solver (row in
  Table 4: Test MSE (×0.01) on Mujoco dataset (lower is better). Results marked with superscript numbers correspond to literature in the footnote.
  Percentage of training data RNN 1 RNN-GRU1 Latent-ODEAdjoint1 Naive2 ACA2 MALI 10% 2.451 1.972 0.471 0.362 0.312 0.35 20% 1.711 1.421 0.441 0.302 0.272 0.27 50% 0.791 0.751 0.401 0.292 0.262 0.26
  Table 5: Test ACC on Speech Command Dataset
  Method Accuracy (%) Adjoint3 92.8± 0.4
  SemiNorm3 92.9± 0.4 Naive 93.2± 0.2 ACA 93.2± 0.2 MALI 93.7± 0.3
  Table 3), and inference on the perturbed images using various solvers. For different combinations of solvers and perturbation amplitudes, Neural ODE consistently outperforms ResNet.
  Summary In image recognition tasks, we demonstrate Neural ODE is accurate, invariant to discretization scheme, and more robust to adversarial attack than ResNet. Note that detailed explanation on the robustness of Neural ODE is out of the scope for this paper, but to our knowledge, MALI is the first method to enable training of Neural ODE on large datasets due to constant memory cost.
sections.heading: 4.3 TIME-SERIES MODELING
sections.text: We apply MALI to latent-ODE (Rubanova et al., 2019) and Neural Controlled Differential Equation (Neural CDE) (Kidger et al., 2020a;b). Our experiment is based on the official implementation from the literature. We report the mean squared error (MSE) on the Mujoco test set in Table 4, which is generated from the “Hopper” model using DeepMind control suite (Tassa et al., 2018); for all experiments with different ratios of training data, MALI achieves similar MSE to ACA, and both outperform the adjoint and naive method. We report the test accuracy on the Speech Command dataset for Neural CDE in Table 5; MALI achieves a higher accuracy than competing methods.
sections.heading: 4.4 CONTINUOUS GENERATIVE MODELS
sections.text: We apply MALI on FFJORD (Grathwohl et al., 2018), a free-from continuous generative model, and compare with several variants in the literature (Finlay et al., 2020; Kidger et al., 2020a). Our experiment is based on the official implementaion of Finlay et al. (2020); for a fair comparison, we train with MALI, and test with the same solver as in the literature (Grathwohl et al., 2018; Finlay et al., 2020), the Dopri5 solver with rtol = atol = 10−5 from the torchdiffeq package (Chen et al., 2018). Bits per dim (BPD, lower is better) on validation set for various datasets are reported in Table 6. For continuous models, MALI consistently generates the lowest BPD, and outperforms the Vanilla FFJORD (trained with adjoint), RNODE (regularized FFJORD) and the SemiNorm Adjoint (Kidger et al., 2020a). Furthermore, FFJORD trained with MALI achieves comparable BPD to stateof-the-art discrete-layer flow models in the literature. Please see Sec. B.3 for generated samples.
sections.heading: 5 RELATED WORKS
sections.text: Besides ALF, the symplectic integrator (Verlet, 1967; Yoshida, 1990) is also able to reconstruct trajectory accurately, yet it’s typically restricted to second order Hamiltonian systems (De Almeida, 1990), and are unsuitable for general ODEs. Besides aforementioned methods, there are other methods for gradient estimation such as interpolated adjoint (Daulbaev et al., 2020) and spectral method (Quaglino et al., 2019), yet the implementations are involved and not publicly available. Other works focus on the theoretical properties of Neural ODEs (Dupont et al., 2019; Tabuada & Gharesifard, 2020; Massaroli et al., 2020). Neural ODE is recently applied to stochastic differential equation (Li et al., 2020), jump differential equation (Jia & Benson, 2019) and auto-regressive models (Wehenkel & Louppe, 2019).
sections.heading: 6 CONCLUSION
sections.text: Based on the asynchronous leapfrog integrator, we propose MALI to estimate the gradient for Neural ODEs. To our knowledge, our method is the first to achieve accuracy, fast speed and a constant memory cost. We provide comprehensive theoretical analysis on its properties. We validate MALI
  01. Rubanova et al. (2019); 2. Zhuang et al. (2020); 3. Kidger et al. (2020a); 4. Chen et al. (2018); 5. Finlay et al. (2020); 6. Dinh et al. (2016); 7. Behrmann et al. (2019); 8. Kingma & Dhariwal (2018); 9. Ho et al. (2019); 10. Chen et al. (2019)
  with extensive experiments, and achieved new state-of-the-art results in various tasks, including image recognition, continuous generative modeling, and time-series modeling.
sections.heading: 7 ACKNOWLEDGEMENT
sections.text: This research was funded by the National Institutes of Health (NINDS-R01NS035193)
sections.heading: CONTENTS (APPENDIX)
sections.heading: A Theoretical properties of ALF integrator 13
sections.text: A.1 Algorithm of ALF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 A.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 A.3 Local truncation error of ALF . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 A.4 Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 A.5 Damped ALF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
sections.heading: B Experimental Details 18
sections.text: B.1 Image Recognition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
  B.1.1 Experiment on Cifar10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 B.1.2 Experiments on ImageNet . . . . . . . . . . . . . . . . . . . . . . . . . . 18
  B.2 Time series modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 B.3 Continuous generative models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
  B.3.1 Training details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 B.3.2 Addtional results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
  B.4 Error in gradient estimation for toy examples when t < 1 . . . . . . . . . . . . . . 20 B.5 Results of damped MALI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
sections.heading: A THEORETICAL PROPERTIES OF ALF INTEGRATOR
sections.heading: A.1 ALGORITHM OF ALF
sections.text: For the ease of reading, we write the algorithm for ψ in ALF below, which is the same as Algo. 2 in the main paper, but uses slightly different notations for the ease of analysis.
  Algorithm 1: Forward of ψ in ALF Input (ẑin, v̂in, sin, h) = (ẑ0, v̂0, s0, h) where s0 is current time, ẑ0 and v̂0 are correponding values at time s0; stepsize h. Forward
  s1 = s0 + h/2 (1) ẑ1 = ẑ0 + v̂0 × h/2 (2) v̂1 = f(ẑ1, s1) (3) v̂2 = v̂1 + (v̂1 − v̂0) (4) ẑ2 = ẑ1 + v̂2 × h/2 (5) s2 = s1 + h/2 (6)
  Output (ẑout, v̂out, sout, h) = (ẑ2, v̂2, s2, h)
  For simplicity, we can re-write the forward of ALF as[ ẑ2
  v̂2
  ] = ẑ0 + hf(ẑ0 + h2 v̂0, s0 + h2 ) 2f(ẑ0 + h 2 v̂0, s0 + h 2 )− v̂0  (7) Similarly, the inverse of ALF can be written as[
  ẑ0 v̂0
  ] = ẑ2 − hf(ẑ2 − h2 v̂2, s2 − h2 ) 2f(ẑ2 − h2 v̂2, s2 − h 2 )− v̂2  (8)
sections.heading: A.2 PRELIMINARIES
sections.text: For an ODE of the form dz(t)
  dt = f(z(t), t) (9)
  We have:
  d2z(t)
  dt2 =
  d dt f(z(t), t) = ∂f(z(t), t) ∂t + ∂f(z(t), t) ∂z dz(t) dt (10)
  For the ease of notation, we re-write Eq. 10 as
  d2z(t)
  dt2 = ft + fzf (11)
  where ft and fz represents the partial derivative of f w.r.t t and z respectively.
sections.heading: A.3 LOCAL TRUNCATION ERROR OF ALF
sections.text: Theorem A.1 (Theorem 3.1 in the main paper). For a single step in ALF with stepsize h, the local truncation error of z is O(h3), and the local truncation errof of v is O(h2).
  Proof. Under the same notation as Algo. 1, denote the ground-truth state of z and v starting from (ẑ0, s0) as z̃ and ṽ respectively. Then the local truncation error is
  Lz = z̃(s0 + h)− ẑ2, Lv = ṽ(s0 + h)− v̂2 (12)
  We estimate Lz and Lv in terms of polynomial of h.
  Under mild assumptions that f is smooth up to 2nd order almost everywhere (this is typically satisfied with neural networks with bounded weights), hence Taylor expansion is meaningful for f . By Eq. 11, the Taylor expansion of z̃ around point (ẑ0, v̂0, s0) is
  z̃(s0 + h) = ẑ0 + h dz dt + h2 2 d2z dt2 +O(h3) (13)
  = ẑ0 + hf(ẑ0, s0) + h2
  2
  ( ft(ẑ0, s0) + fz(ẑ0, s0)f(ẑ0, s0) ) +O(h3) (14)
  Next, we analyze accuracy of the numerical approximation. For simplicity, we directly analyze Eq. 7 by performing Taylor Expansion on f .
  f(ẑ0 + h
  2 v̂0, s0 +
  h 2 ) = f(ẑ0, s0) + h 2 ft(ẑ0, s0) + hv̂0 2 fz(ẑ0, s0) +O(h 2) (15)
  ẑ2 = ẑ0 + hf(ẑ0 + h
  2 v̂0, s0 +
  h 2 ) (16)
  Plug Eq. 14, Eq. 15 and E.q. 16 into the definition of Lz , we get
  Lz = z̃(s0 + h)− ẑ2 (17) = [ ẑ0 + hf(ẑ0, s0) + h2
  2
  ( ft(ẑ0, s0) + fz(ẑ0, s0)f(ẑ0, s0) )] − [ ẑ0 + h ( f(ẑ0, s0) + h
  2 ft(ẑ0, s0) + hv̂0 2 fz(ẑ0, s0)
  )] +O(h3) (18)
  = h2
  2 fz(ẑ0, s0)
  ( f(ẑ0, s0)− v̂0 ) +O(h3) (19)
  Therefore, if ∣∣∣f(ẑ0, s0)− v̂0∣∣∣ is of order O(1), Lz is of order O(h2); if ∣∣∣f(ẑ0, s0)− v̂0∣∣∣ is of order O(h) or smaller, then Lz is of order O(h3). Specifically, at the start time of integration, we have∣∣∣f(ẑ0, s0)− v̂0 = 0∣∣∣, by induction, Lz at end time is O(h3). Next we analyze the local truncation error in v, denoted as Lv . Denote the ground truth as ṽ(t0 +h), we have
  ṽ(s0 + h) = f ( z̃(s0 + h), s0 + h ) (20)
  = f(ẑ0, s0) + hft(ẑ0, s0) + ( z̃(s0 + h)− ẑ0 ) fz(ẑ0, s0) +O(h 2) (21)
  Next we analyze the error in the numerical approximation. Plug Eq. 15 into Eq. 7,
  v̂2 = 2f(ẑ0 + h
  2 v̂0, s0 +
  h 2 )− v̂0 (22)
  = f(ẑ0, s0) + ( f(ẑ0, s0)− v̂0 ) + hft(ẑ0, s0) + hv̂0fz(ẑ0, s0) +O(h 2) (23)
  From Eq. 14, Eq. 21 and Eq. 23, we have
  Lv = ṽ(s0 + h)− v̂2 (24) = ( f(ẑ0, s0)− v̂0 ) + ( z̃(s0 + h)− ( ẑ0 + hv̂0 )) fz(ẑ0, s0) +O(h 2) (25)
  = ( f(ẑ0, s0)− v̂0 ) + h ( f(ẑ0, s0)− v̂0 ) fz(ẑ0, s0) +O(h 2) (26)
  The last equation is derived by plugging in Eq. 14. Note that Eq. 26 holds for every single step forward in time, and at the start time of integration, we have ∣∣f(ẑ0, s0) − v̂0∣∣ = 0 due to our initialization as in Sec. 3.1 of the main paper. Therefore, by induction, Lv is of order O(h2) for consecutive steps.
sections.heading: A.4 STABILITY ANALYSIS
sections.text: Lemma A.1.1. For a matrix of the form [ A B C D ] , if A,B,C,D are square matrices of the same
  shape, and CD = DC, then we have det [ A B C D ] = det(AD −BC)
  Proof. See (Silvester, 2000) for a detailed proof.
  Theorem A.2. For ALF integrator with stepsize h, if hσi is 0 or is imaginary with norm no larger than 1, where σi is the i-th eigenvalue of the Jacobian ∂f∂z , then the solver is on the critical boundary of A-stability; otherwise, the solver is not A-stable.
  Proof. A solver is A-stable is equivalent to the eigenvalue of the numerical forward has a norm below 1. We calculate the eigenvalue of ψ below.
  For the function defined by Eq. 7, the Jacobian is
  J =  ∂ẑ2∂z0 ∂ẑ2∂v̂0 ∂v̂2 ∂z0 ∂v̂2 ∂v̂0  = I + h∂f∂z h22 ∂f∂z 2× ∂f∂z h ∂f ∂z − I  (27) We determine the eigenvalue of J by solving the equation
  det(J − λI) = h∂f∂z + (1− λ)I h22 ∂f∂z 2× ∂f∂z h ∂f ∂z − (1 + λ)I  = 0 (28) It’s trivial to check J satisfies conditions for Lemma A.1.1.Therefore, we have
  det(J − λI) = det [( h ∂f
  ∂z + (1− λ)I
  )( h ∂f
  ∂z − (1 + λ)I
  ) − (h2
  2
  ∂f
  ∂z
  )( 2× ∂f
  ∂z
  )] (29)
  = det [ − 2λh∂f
  ∂z + (λ2 − 1)I
  ] (30)
  Suppose the eigen-decompostion of ∂f∂z can be written as
  ∂f ∂z = Λ σ1 σ2 ... σN Λ−1 (31) Note that I = ΛIλ−1, hence we have
  det(J − λI) = det Λ { − 2λh σ1 σ2 ... σN + (λ2 − 1)I}Λ−1 (32) =
  N∏ i=1 (λ2 − 2hσiλ− 1) (33)
  Hence the eigenvalues are
  λi± = hσi ± √ h2σ2i + 1 (34)
  A-stability requires |λi±| < 1,∀i, and has no solution. The critical boundary is |λi±| = 1, the solution is: hσi is 0 or on the imaginary line with norm no larger than 1.
sections.heading: A.5 DAMPED ALF
sections.text: Algorithm 2: Forward of ψ in Damped ALF (η ∈ (0, 1] ) Input (ẑin, v̂in, sin, h) = (ẑ0, v̂0, s0, h) where s0 is current time, ẑ0 and v̂0 are correponding values at time s0; stepsize h. Forward
  s1 = s0 + h/2 (35) ẑ1 = ẑ0 + v̂0 × h/2 (36) v̂1 = f(ẑ1, s1) (37) v̂2 = v̂0 + 2η(v̂1 − v̂0) (38) ẑ2 = ẑ1 + v̂2 × h/2 (39) s2 = s1 + h/2 (40)
  Output (ẑout, v̂out, sout, h) = (ẑ2, v̂2, s2, h)
  Algorithm 3: ψ−1 (Inverse of ψ) in Damped ALF (η ∈ (0, 1] ) Input (ẑout, v̂out, sout, h) where sout is current time, ẑout and v̂out are corresponding values at sout, h is stepsize. Inverse
  (ẑ2, v̂2, s2, h) = (ẑout, v̂out, sout, h) (41) s1 = s2 − h/2 (42) ẑ1 = z2 − v̂2 × h/2 (43) v̂1 = f(ẑ1, s1) (44) v̂0 = (v̂2 − 2ηv̂1)/(1− 2η) (45) ẑ0 = ẑ1 − v̂0 × h/2 (46) s0 = s1 − h/2 (47)
  Output (ẑin, v̂in, sin, h) = (ẑ0, v̂0, s0, h)
  The main difference between ALF and Damped ALF is marked in blue in Algo. 2. In ALF, the update of v̂2 is v̂2 = (v̂1 − v̂0) + v̂1 = 2(v̂1 − v̂0) + v̂0; while in Damped ALF, the update is scaled by a factor η between 0 and 1, so the update is v̂2 = 2η(v̂1 − v̂0) + v̂0. When η = 1, Damped ALF reduces to ALF.
  Similar to Sec. A.1, we can write the forward as For simplicity, we can re-write the forward of ALF as [
  ẑ2
  v̂2
  ] = ẑ0 + ηhf(ẑ0 + h2 v̂0, s0 + h2 ) + (1− η)hv̂0 2ηf(ẑ0 + h 2 v̂0, s0 + h 2 ) + (1− 2η)v̂0  (48) Similarly, the inverse of ALF can be written as
  [ ẑ0
  v̂0
  ] = ẑ2 − h 1−η1−2η v̂2 + h η1−2ηf(ẑ2 − h2 v̂2, s2 − h2 ) 1
  1−2η v̂2 − 2η 1−2ηf(ẑ2 − h 2 v̂2, s2 − h 2 )  (49) Theorem A.3. For a single step in Damped ALF with stepsize h, the local truncation error of z is O(h2), and the local truncation errof of v is O(h).
  Proof. The proof is similar to Thm. A.3. By similar calculations using the Taylor Expansion in Eq. 15 and Eq. 14, we have
  ẑ2 − z̃(s0 + h) = (1− η)hv̂0 + hη [ f(ẑ0, s0) + h
  2 ft(ẑ0, s0) + hv̂0 2 fz(ẑ0, s0) ] − h [ f(ẑ0, s0) + h
  2 ftẑ0, s0 +
  h 2 fz(ẑ0, s0)f(ẑ0, s0)
  ] +O(h2) (50)
  = (1− η)h ( v̂0 − f(ẑ0, s0) ) + η − 1
  2 h2ft(ẑ0, s0)
  + h2
  2
  ( ηv̂0 − f(ẑ0, s0) ) fz(ẑ0, s0) +O(h 2) (51)
  Using Eq. 21, Eq. 15 and Eq. 14, we have
  ṽ2 − v̂2 = (1− 2η)v̂0 + (2η − 1)f(ẑ0, s0) + (1− η)hft(ẑ0, s0) + ( z̃(s0 + h)− ẑ0 − ηhv̂0 ) fz(ẑ0, s0) +O(h 2) (52)
  = (2η − 1) [ f(ẑ0, s0)− ẑ0 ] + (1− η)hft(ẑ0, s0)
  + η [ hf(ẑ0, s0)− hv̂0 ] fz(ẑ0, s0) +O(h 2) (53)
  Note that when η = 1, Eq. 51 reduces to Eq. 19, and Eq. 53 reduces to Eq. 26. By initialization, we have |f(ẑ0, s0) − v̂0| = 0 at initial time, hence by induction, the local truncation error for z is O(h2); the local truncation error for v is O(h) when η < 1, and is O(h2) when η = 1.
  Theorem A.4 (Theorem 3.2 in the main paper). For Dampled ALF integrator with stepsize h, where σi is the i-th eigenvalue of the Jacobian ∂f∂z , then the solver is A-stable if ∣∣∣1 + η(hσ − 1)± √ η [ 2hσi + η(hσi − 1)2
  ]∣∣∣ < 1, ∀i. Proof. The Jacobian of the forward-pass of a single step damped ALF is
  J = I + ηh∂f∂z (1− η)hI + η h22 ∂f∂z 2η ∂f∂z ηh ∂f ∂z + (1− 2η)I  (54) when η = 1, J reduces to Eq. 27. We can determine the eigenvalue of J using similar techniques. Assume the eigenvalues for ∂f∂z are {σi}, then we have
  det(J − λI) = det (1− λ)I + ηh∂f∂z (1− η)hI + η h22 ∂f∂z 2η ∂f∂z ηh ∂f ∂z + (1− 2η − λ)I  (55) = det [( (1− λ)I + ηh∂f
  ∂z
  )( ηh ∂f
  ∂z + (1− 2η − λ)I ) − ( (1− η)hI + ηh 2
  2
  ∂f
  ∂z
  )( 2η ∂f
  ∂z
  )] (56)
  = N∏ i=1 [ 1 + η(hσi − 1)± √ η [ 2hσi + η(hσi − 1)2 ]] (57)
  when η < 1, it’s easy to check that ∣∣∣1+η(hσi−1)±√η[2hσi + η(hσi − 1)2]∣∣∣ < 1 has non-empty
  solutions for hσ.
  For a quick validation, we plot the region of A-stability on the imaginary plane for a single eigenvalue in Fig. 1. As η increases, the area of stability decreases. When η = 1, the system is no-where A-stable, and the boundary for A-stability is on the imaginary axis [−i, i] where i is the imaginary unit.
sections.heading: B EXPERIMENTAL DETAILS
sections.text: B.1 IMAGE RECOGNITION
sections.heading: B.1.1 EXPERIMENT ON CIFAR10
sections.text: We directly modify a ResNet18 into a Neural ODE, where the forward of a residual block (y = x + f(x)) and the forward of an ODE block (y = x + ∫ T 0 f(z, t)dt where T = 1) share the same parameterization f , hence they have the same number of parameters. Our experiment is based on the official implementation by Zhuang et al. (2020) and an open-source repository (Liu, 2017).
  All models are trained with SGD optimizer for 90 epochs, with an initial learning rate of 0.01, and decayed by a factor of 10 at 30th epoch and 60th epoch respectively. Training scheme is the same for all models (ResNet, Neural ODE trained with adjoint, naive, ACA and MALI). For ACA, we follow the settings in (Zhuang et al., 2020) and use the official implementation torch ACA 1, and use a Heun-Euler solver with rtol = 10−1, atol = 10−2 during training. For MALI, we use an adaptive version and set rtol = 10−1, atol = 10−2. For the naive and adjoint method, we use the default Dopri5 solver from the torchdiffeq2 package with rtol = atol = 10−5. We train all models for 5 independent runs, and report the mean and standard deviation across runs.
sections.heading: B.1.2 EXPERIMENTS ON IMAGENET
sections.text: Training scheme We conduct experiments on ImageNet with ResNet18 and Neural-ODE18. All models are trained on 4 GTX-1080Ti GPUs with a batchsize of 256. All models are trained for 80 epochs, with an initial learning rate of 0.1, and decayed by a factor of 10 at 30th and 60th epoch. Note that due to the large size input 256 × 256, the naive method and ACA requires a huge memory, and is infeasible to train. MALI and the adjoint method requires a constant memory hence is suitable for large-scale experiments. For both MALI and the adjoint menthod, we use a fixed stepsize of 0.25, and integrates from 0 to T = 1. As shown in Table. 2 in the main paper, a stepsize of 0.25 is sufficiently small to train a meaningful continuous model that is robust to discretization scheme.
  Invariance to discretization scheme To test the influence of discretization scheme, we test our Neural ODE with different solvers without re-training. For fixed-stepsize solvers, we tested various step sizes including {0.1, 0.15, 0.25, 0.5, 1.0}; for adaptive solvers, we set rtol=0.1, atol=0.01 for MALI and Heun-Euler method, and set rtol = 10−2, atol = 10−3 for RK23 solver, and set rtol = 10−4, atol = 10−5 for Dopri5 solver. As shown in Table. 2, Neural ODE trained with MALI is robust to discretization scheme, and MALI significantly outperforms the adjoint method in terms of accuracy (70% v.s. 63% top-1 accuracy on the validation dataset). An interesting finding is that when trained with MALI which is a second-order solver, and tested with higher-order solver (e.g.
  1https://github.com/juntang-zhuang/torch_ACA 2https://github.com/rtqichen/torchdiffeq
  RK4), our Neural ODE achieves 70.21% top-1 accuracy, which is higher than both the same solver during training (MALI, 69.59% accuracy) and the ResNet18 (70.09% accuracy).
  Furthermore, many papers claim ResNet to be an approximation for an ODE (Lu et al., 2018). However, Queiruga et al. (2020) argues that many numerical discretizations fail to be meaningful dynamical systems, while our experiments demonstrate that our model is continuous hence invariant to discretization schemes.
  Adversarial robustness Besides the high accuracy and robustness to discretization scheme, another advantage of Neural ODE is the robustness to adversarial attack. The adversary robustness of Neural ODE is extensively studied in (Hanshu et al., 2019), but not only validated on small-scale datasets such as Cifar10. To our knowledge, our method is the first to enable effectuve training of Neural ODE on large-scale datasets such as ImageNet and achieve a high accuracy, and we are the first to validate the robustness of Neural ODE on ImageNet. We use the advertorch 3 toolbox to perform adversarial attack. We test the performance of ResNet and Neural ODE under FGSM attack. To be more convincing, we conduct experiment on the pretrained ResNet18 provided by the official PyTorch website 4. Since Neural ODE is invariant to discretization scheme, it’s possible to derive the gradient for attack using one ODE solver, and inference on the perturbed image using another solver. As summarized in Table. 3, Neural ODE consistently achieves a higher accuracy than ResNet under the same attack.
sections.heading: B.2 TIME SERIES MODELING
sections.text: We conduct experiments on Latent-ODE models (Rubanova et al., 2019) and Neural CDE (controlled differential equation) (Kidger et al., 2020a). For all experiments, we use the official implementation, and only replace the solver with MALI. The latent-ODE model is trained on the Mujoco dataset processed with code provided by the official implementation, and we experiment with different ratios (10%,20%,50%) of training data as described in (Rubanova et al., 2019). All models are trained for 300 epochs with Adamax optimizer, with an initial learning rate of 0.01 and scaled by 0.999 for each epoch. For the Neural CDE model, for the naive method, ACA and MALI, we perform 5 independent runs and report the mean value and standard deviation; results for the adjoint and seminorm adjoint are from (Kidger et al., 2020a). For Neural CDE, we use MALI with ALF solver with a fixed stepsize of 0.25, and train the model for 100 epochs with an initial learning rate of 0.004.
  3https://github.com/BorealisAI/advertorch 4https://pytorch.org/docs/stable/torchvision/models.html
sections.heading: B.3 CONTINUOUS GENERATIVE MODELS
sections.heading: B.3.1 TRAINING DETAILS
sections.text: Our experiment is based on the official implementation of (Finlay et al., 2020), with the only difference in ODE solver. For a fair comparison, we only use MALI for training, and use Dopri5 solver from torchdiffeq package (Chen et al., 2018) with rtol = atol = 10−5. For MALI, we use adaptive ALF solver with rtol = 10−2, atol = 10−3, and use an initial stepsize of 0.25. Integration time is from 0 to 1.
  On MNIST and CIFAR dataset, we set the regularization coefficients for kinetic energy and Frobenius norm of the derivative function as 0.05. We train the model for 50 epochs with an initial learning rate of 0.001.
sections.heading: B.3.2 ADDTIONAL RESULTS
sections.text: We show generated examples on MNIST dataset in Fig. 3, results for Cifar10 dataset in Fig. 4, and results for ImageNet64 in Fig. 5.
  B.4 ERROR IN GRADIENT ESTIMATION FOR TOY EXAMPLES WHEN t < 1
  We plot the error in gradient estimation for the toy example defined by Eq.6 in the main paper in Fig. 6. Note that the integration time T is set as smaller than 1, while the main paper is larger than 20. We observe the same results, MALI and ACA generate smaller error than the adjoint and the naive method.
sections.heading: B.5 RESULTS OF DAMPED MALI
sections.text: For all experiments in the main paper, we set η = 1 and did not use damping. For completeness, we experimented with damped MALI using different values of η. As shown in Table. 7, MALI is robust to different η values.
title: MALI: A MEMORY EFFICIENT AND REVERSE ACCU- RATE INTEGRATOR FOR NEURAL ODES