Global Convergence in Neural ODEs: Impact of Activation Functions
Authors: Tianxiang Gao
Siyuan Sun
Hailiang Liu
Hongyang Gao
Venue: ICLR 2025 Oral (2025)
Keywords: Neural ODEs
Gradient Descent
Neural Tangent Kernel (NTK)
Primary area: learning theory

Abstract: Neural Ordinary Differential Equations (ODEs) have been successful in various applications due to their continuous nature and parameter-sharing efficiency. However, these unique characteristics also introduce challenges in training, particularly with respect to gradient computation accuracy and convergence analysis. In this paper, we address these challenges by investigating the impact of activation functions. We demonstrate that the properties of activation functions—specifically smoothness and nonlinearity—are critical to the training dynamics. Smooth activation functions guarantee globally unique solutions for both forward and backward ODEs, while sufficient nonlinearity is essential for maintaining the spectral properties of the Neural Tangent Kernel (NTK) during training. Together, these properties enable us to establish the global convergence of Neural ODEs under gradient descent in overparameterized regimes. Our theoretical findings are validated by numerical experiments, which not only support our analysis but also provide practical guidelines for scaling Neural ODEs, potentially leading to faster training and improved performance in real-world applications.

Published as a conference paper at ICLR 2025
GLOBAL CONVERGENCE IN NEURAL ODES: IMPACT OF ACTIVATION FUNCTIONS
Tianxiang Gao DePaul University tgao9edepaulledd edu
Styuan Sun
Hailiang Liu
lowa State University sxs 44730 assate.edd
lowa State University h1 uu lastate... edu
Hongyang Gao Iowa State University hygaoliastate.
ABSTRACT
Neural Ordinary Differential Equations (ODEs) have been successful in vari- applications due to their continuous nature and parameter-sharing efficiency. However, these unique characteristics also introduce challenges training, par- ticularly with respect to gradient computation accuracy and convergence analysis. In this paper, we address these challenges by investigating the impact of activation functions. We demonstrate that the properties of activation functions caartions--pecfcccalll smoothness and nonlincarity- are critical to the training dynamics. Smooth acti- vation functions guarantee globally unique solutions for both forward and back- ward ODEs, while sufficient nonlinearity is essential for maintaining the spectral properties of the Neural Tangent Kernel (NTK) during training. Together, these properties enable us to establish the global convergence of Neural ODEs under gradient descent overparameterized regimes. Our theoretical findings are val- idated by numerical experiments, which not only support analysis but also provide practical guidelines for scaling Neural ODEs, potentially leading to faster training and improved performance in real-world applications.
INTRODUCTION
In recent years, deep neural networks have achieved remarkable success across wide range of appli- cations. Among these advancements, Neural Ordinary Differential Equations (ODEs) (Chen al.] 2018b) stand out due to their continuous nature and parameter efficiency through shared parame- ters. Unlike conventional neural networks with discrete layers, Neural ÓDEs model the evolution of hidden states continuous-time differential equation, allowing them to better capture dynamic systems. This parameter serraietr--shring mechanism ensures consistent dynamics throughout the continu- transformation and reduces the number of parameters, improving both memory efficiency and computational complexity. These unique properties make Neural ODEs particularly effective not only for traditional machine learning tasks like image classification (Chen al. (0188) and natural language processing (Rubanova et al. 2019), but also for more complex tasks involving continu- ous processes, such as time series analysis (Kidger al., 2020), reinforcement learning (Du al. 2020], and diffusion models (Song et al., 2020). However, while these features offer flexibility and efficiency, they also introduce sigmificant challenges during training, particularly in gradient compu- tation and convergence analysis.
One of the key challenges in training Neural ODEs is accurately computing gradients. Unlike tradi- tional networks, where backpropagation can be computed through a discrete chain of layers, Neural ODEs require solving forward and backward ODEs using numerical solvers. These solvers intro- duce numerical errors, which can lead to inaccurate gradients and slow convergence even sub- optimal model performance (Rodriguez et al.] 2022). Moreover, ensuring the well well-posedness of Theorem, solutions
ODE solutions during training is nontrivial. According the may not always exist or may only exist locally, which can cause training divergence or significant numerical errors (Gholami et al. 2019 Ott et al. 2020 Sander al. 2022). Even with advanced

---

Published as a conference paper at ICLR 2025
solvers (Zhuang et al., 2020a b; Matsubara et al. 2021 Ko et al. 2023). it remains open problem whether simple first-order methods, such as stochastic gradient descent (SGD), can reliably train Neural ODEs to convergence. While discretizing Neural ODEs finite-depth networks offers po- tential solution, it results in a deeper computation graph (Zhuang et al. 2020a b), raising questions about whether the gradients computed in this manner truly match those of the continuous model.
Another essential challenge lies in analyzing the training dynamics of Neural ODEs. The optimiza- tion problem in training neural networks is inherently nonconvex, making theoretical analysis diffi- cult. Recent work by Jacot et al. (2018) has shown that the training dynamics of overparameterized networks can be understood through the lens of the Neural Tangent Kernel (NTK), which converges to a deterministic limit as network width increases. This convergence has enabled researchers to methods overparameterized regimes,
establish global convergence guarantees for gradient-basee provided the NTK remains strictly positive definite (SPD) (Du et al. 2019a Allen-Zhu et al.) 2019) Nguyen| 2021, Gao et al.) 2021). The analysis of the NTK's 's strict positive definiteness began with Daniely et al. 2016), who introduced the concept of dual activation for two-layer networks, later extended to deeper, finite networks (Jacot et al. 2018: Du et al., 2019a). However, these results are limited networks with discrete layers, raising the question of whether the same properties hold for continuous models like Neural ODEs.
In this paper, address these challenges by exploring the impact of activation functions on training Neural ODEs. We show that activation function ity-play critical roles in determining the well-posedness of ODE solutions and the spectral prop- erties of the NTK. Through our analysis, we demonstrate that smooth activation functions lead to globally unique solutions for both forward and backward ODEs, ensuring the stability of the training process. Additionally, we extend existing results on the NTK from discrete-layered neural networks to continuous models, demonstrating that the NTK for Neural ODEs is well-defined. Importantly we find that a higher degree of nonlinearity in the activation function not only helps maintain the SPD property of the limiting NTK, but also practically speeds up Neural ODE convergence.
smoothness and nonlinear-
1. We investigate the significance of the smoothness of activation functions for the well-posedness, of forward and backward ODEs in Neural ODEs. Using random matrix theory, we demonstrate the existence of globally unique solutions. Additionally, we show that no additional regularity is needed forward and backward ODEs are combined in a weakly coupled ODE system.
2. We propose new mathematical framework for studying continuous models from the approxi- mation theory perspective. By using a sequence of finite-depth neural networks to approximate Neural ODEs, show that key properties like activation and gradient propagation are preserved as depth approaches infinity. This allows us to apply the Moore-Osgood theorem from functional analysis to prove that the NTK of Neural ODEs well-defined.
3. Unfortunately, the SPD property of the NTK may not hold at infinite depth, even we can show every finite-depth approximation satisfies it. To address this, we conduct a fine-grained analysis and derive an integral form for the limiting NTK of Neural ODEs. This form reveals that the NTK remains SPD if the activation function is non-polynomial Leveraging this integral representation provides valuable insights into continuous models and may inspire further research.
4.
We conduct a series of numerical experiments to support our theoretical findings. Beyond validat- ing our analysis, these experiments also provide practical guidelines for training Neural ODEs. We show that activation function smoothness and nonlinearity accelerate convergence and im- prove performance (see Figure 7 Conversely, improper ODE scaling leads to damping from accumulated numerical errors (see Figure while adaptive solvers struggle with efficiency in large-scale Neural ODEs (see Figure causing instability and high computational overhead.
2
PRELIMINARIES
2.1
NEURAL ODES
In this paper, we consider a simple Neural ODE (a; 0) defined as follows dv o ((r)),
fl: 0)
(1)
Vn
'The general form of Neural ODEs is discussed in Appendix
2

---

Published as a conference paper at ICLR 2025
where he R is the hidden state that satisfies the following ordinary differential equation
hg C...U//VVa and h, a. W Wheh)//T, AI E [0,
where é is the activation function a F is input, U Raxd W E R"x", and E R" are learnable parameters. These parameters. denoted by 0 == vec v), are randomly initialized
(Glorot & Bengio 2010 He et al. 2015) from standard Gaussian distribution:
Uzj
Wej,
Vi ltd. 10,
with variance hyperparameters Ju. Jwr Ov > Due to the continuous nature of Neural ODEs, computing gradients through standard backpropagation is not feasible. Instead, method (Chen et al. 2018b) to compute gradients by solving the backward ODE:
the adjoint
Ay = di (hz))v/ Vn, and A = di (h+)) ( Wt////n
(4)
where A is the adjoint state. When both forward and backward ODEs well-posed, we can compute the gradients of So with respect to (w.r.t.) the parameters 0 as follows
vof(2; 0) -- (ht)) Vwf(x; ['" On (he) dtt Vuf(w; = Ca An Vn
(5)
Vn
Further details on these derivations are provided in Appendix B In Section 3 we demonstrate that if é is Lipschitz continuous, the forward and backward ODEs have globally unique solutions. Ad-
ditionally, in S we prove this well poll-posssssss holds throughout the entire training process.
2.2 NEURAL TANGENT KERNEL
Given a training dataset {(0;, y:)] the objective is to learn a parameter 0 that minimizes the empirical loss:
L( Σ 1 (f(x; x1)² --- -1)²
221lu yll?.
(6)
where (u,,... uN) is the prediction vector with u, f(z;; 0), and y (v1.... MNN is the output vector. Gradient descent with a learning rate 17 > 0 is used to minimize the loss:
g++1
0² L
Following (Du et al. 2019a) and (Jacot et al. 2018). under some required conditions, the evolution of the prediction vector u" can be approximated as follows:
k+1 y ~ (I - hH) (u* y),
(8)
where H ER
ERR NxN is Gram matrix defined through the NTK (Jacot et al. 2018):
K K(x, ē; 9) f( 0), Vof(=;)))
The training dynamics Eq. 8 is governed by the spectral property of the NTK Gram matrix H If there exists a strictly positive constant Ao > 0 s.t. Amin (H*) 2 An for all k. then the residual ² -.--1/ 3/ decreases to zero exponentially, provided the learning rate 1j > 0 sufficiently small (Allen- Zhu et al. 2019: Nguyen 2021). As the parameters 0² are updated during training, the NTK Ke changes over time, making its spectral analysis challenging. Fortunately, previous research (Yang) 2020 Jacot et al. 2018) shows that the time-varying Ko converges to deterministic limiting NTK Rxc as the network width n - x. By leveraging this result and the concept of dual activation (Daniely et al. (2016), can study the spectral properties of Kg during training through Ko0 using perturbation theory.
However, these prior results apply only finite-depth neural networks, while Neural ODEs infinite-depth networks due to their continuous nature. Moreover, as prior analyses are based on induction techniques, there is no guarantee that these essential properties will also hold as depth tends to infinity. In Section 4.1 and Section 4.2 we introduce a new framework to study Neural ODEs as infinite-depth networks. We demonstrate that the smoothness of the Φ is crucial to ensure these essential properties hold in Neural ODEs. Additionally, to study the spectral properties of the limiting NTK Koc. we provide a fine-grained analysis by expressing the limiting NTK of Neural ODEs in an integral form, which we conclude that the nonlinearity of activation plays a critical role ensuring the strict positive definiteness of K
To simplify the analysis, assume the activation threshold value for ø and further (00) 0.
3

---

Published as a conference paper at ICLR 2025
3 WELL-POSEDNESS OF NEURAL ODES AND ITS GRADIENTS
As continuous models, Neural ODEs pose significant challenge in accurately computing gradi- ents. In this section, we explore the challenges associated with two methods for computing gra- dients: optimize--hern-disreliizz and discrelize-then-opiiiize (Gholami al. 2019 Finlay et al,) 2020 Onken et al. 2021). Through our exploration, we emphasize the essential role of smoothness activation functions to guarantee the possedossss of Neural ODEs and their gradients.
3.1 OPTIMIZE-HHENDDSSEEE
METHOD
As discussed in Section Neural ODEs require numerical ODE solvers to solve the forward and backward ODEs Eq. 2 and Eq. (4) to compute the gradients, employing a method known as the pptimize-then-dixrreee
method. When solving ODEs, ensuring their well-posedness is of primary concern. In Proposition]] we demonstrate that o is Lipschitz continuous, the forward and back- ward ODEs have globally unique solutions, thus ensuring the well-posedness. The detailed proofs are provided in Appendix
!
Proposition 1. For any given T > O, if the activation function o is Lipschitz continuous with Lipsehitz constant L1. then the forward ODE Eq. (2) and the backward ODE Eq. (4) have unique solutions h2 and X for all E [0. T) and x E Rd almost surely over random initialization Eq. (3)
In addition, (x) (x; 0))/Oh, is the solution to the backward ODE.
Although Neural ODEs and their gradients are well-defined under these conditions, this does not guarantee that gradients can be computed accurately by solving the ODEs numerically. One primary issue is that the magnitudes of the hidden state h, and adjoint state X can grow exponentially over the time horizon T. leading to accumulated numerical errors. This issue is illustrated in Figure [ where long-time horizon leads to damping in the early stages of training. To mitigate this problem, suggest scaling the dynamics by setting dw (1/T), which ensures that the magnitudes of h2 and A, remain mild and independent of T. This scaling stabilizes the norms, thereby allowing numerical solvers to produce much more accurate gradient estimates.
Additionally, calculating gradients using Eq. (5) requires storing the values of h and A at every time step & E [0, T', which can consume a significant amount of memory in practice. To address this, Chen et al. (2018)) propose solving an augmented backward ODE (defined in Appendix (36 for our setup), which combines an additional gradient state with both the backward ODE and the reversed forward ODE. This approach eliminates the need for storing intermediate states. However, since the hidden state h4 is no longer constant in the augmented ODE, additional regularization conditions on the dynamics are typically required to ensure the stability of the solution. Fortunately, we demonstrate that such regularization conditions are unnecessary for Neural ODEs because the Lipschitz continuity of ensures that ha is well defined for all t j0, T). Therefore, the augmented ODE approach can be used without the need for additional regularization. A detailed discussion of this result can be found in Appendix C
3.2 THSCRETIZE-TTEN--OPTTI
METHOD
Alternatively, we can discretize the ODE using Euler's method, treating the continuous ODE as finite-depth Residual Network (ResNet) f= 0 with shared parameters across layers:
(z 0) = T o(h)(x)), Vn
(10a)
h¹ = h f-1 + ow Wo(h--))
E .1,
(10b)
h Tu Ux, Vd
(10c)
where k T/L represents the time step. The gradient can then be estimated using backpropagation through the finite depth ResNet (x; 0), referred to as the discretize tiscrriizeetheenoopiiize approach.
Here the superscript in fo indicates it has I. time steps, while actually fb has totally L 2 layers.
4

---

Published as a conference paper at ICLR 2025
As finite-depth ResNet, the gradients of f' 0) are always well defined. However, it remains open question whether the gradients of the finite-depth approximation f"(x; 0) converge to the gra- dients of the continuous Neural ODE f(; 0) as the depth L + In Proposition we demonstrate that the smoothness of & ensures this convergence. Thus, in the limit of infinite depth (or infinites- imally small time steps), both the optimize- theniizz-tee--dicretiize and discretize-then-opiiii yield the same gradients, provided that the activation function is sufficiently smooth. The detailed proofs are provided Appendix E
methods
Proposition Given a E R. if the activation function & and derivative o' are L1- and Ly- Lipschitz continuous, respectively, the following inequalities hold a.s. over random initialization:
- S CCL 4l (......... L.:
(11)
where C > O is a constant depending only on L1, L2, T, ov ow Ou and ||==]]-
To further validate our theoretical findings, we conduct experiments that compare training efficiency and gradient accuracy with and without the Lipschitz continuity of ø'. These experiment results, il- lustrated Figure 7 demonstrate the necessity of Lipschitz continuity for ensuring smooth gradient computation and achieving faster convergence during training.
4 NNGP AND NTK FOR NEURAL ODES
Understanding how activation and gradient propagate through neural networks is crucial for analyz- ing their training dynamics and generalization abilities, as emphasized in previous studies (Glorol & Bengio 2010 Poole et al. 2016;; Schoenholz et al.) [017). The frameworks of Neural Network Gaussian Processes (NNGP) (Lee al.) 2018) and Neural Tangent Kernels (NTK) (Jacot et al.] 2018), grounded in mean-field theory, provide powerful analytical tools to study these dynamics. In this section, we establish theoretical results that demonstrate the well-defined nature of NNGP and NTK for Neural ODEs and explore their properties with respect to information propagation.
4.1 NNGP: FORWARD PROPAGATION OF INPUTS
Previous work has shown that in the infinite- width limit, randomly initialized fiiia neural net- works converge Gaussian processes with mean zero and recursively defined covariance functions, known as NNGP kernels (Neal 2012 Lee et al. [0088] Daniely et al. 2016) Yang. 2019 Gao et al., 2023 Gao 2024). When approximating Neural ODEs using sequence of finite-depth ResNets fø establish the NNGP for f/ Detailed proofs provided in Appendix D
Proposition 3. Suppose & is Ln-Lipschitz continuous. Then, as width n. - co, the finite-depth neural network fo defined Eq. 10) converges distribution to a centered Gaussian Process with a covariance function or NNGP kernel EL+1 = C1+1,5+1 defined recursively:
(x, .)) = 50, 0- (0.
Vk: E (0, 1, ... L L 1}
(12)
d CC. (xx, ②) Eplu (o
ve,k E (1,2,.... +1}
(13)
where k T/L and (u², w²) are centered Gaussian random variables with covariance
E(u'at) = 00.0 x) x Σ Chd (t,
Ve,k (0, 1. L).
(14)
1-=1
Although Proposition 3 shows that sequence of Gaussian processes (GPs) can be derived for the sequence of f/ this does not necessarily mean that the Neural ODE So will also converge to a Gaus- sian Process peffnire-wiitth--he---eee
oc. The challenge lies in the difference between two convergence patterns: seithe-depttt-he--eee
These often lead to different limits. For example, consider the simple double sequence an := n/(e + n). This double sequence demon- strates how taking different convergence paths-first first in width, then in depth, or vice versa-can yield different results. Specifically, we observe that
and
lim ( lim an.e) 1, while
lim ( lim an,z) = 0. n (-000
E-000 n-000
5

---

Published as a conference paper at ICLR 2025
This phenomenon has been noted in several recent studies (Hayou & Yang 2023: Yang al. 2024
Gao et al. 2023 Gao 2024) across various neural network architectures. Specifically. for commonly used neural networks, the two convergence patterss often do not coincide, leading different limits for infinite-depth networks. Hence, the NNGP correspondence does not generally hold for infinite- depth neural networks. For Neural ODEs, the seflnite-deell--hee----eee relevant, as Neural ODEs are equivalent to infinite-depth neural networks from the standpoint of numerical discretization. However, the continuous nature and parameter sharing in Neural ODEs present unique challenges that make previous mathematical tools inapplicable directly.
convergence pattern is more
Fortunately, if the activation function & is sufficiently smooth, we can show that these two lim- its commute, and both convergence patterns share the same limit. One crucial intermediate result involves proving that the double sequence (d(h²), <b(h")) /n converges in depth I. uniformly in width (almost surely). This uniform convergence ensures that as depth increases, the behavior of the system remains stable regardless of width, which is crucial for showing that the limits commute and establishing the well-posedness of the NNGP for Neural ODEs. The proof relies Euler's convergence theorem and is provided Appendix D
Lemma 1. Let á be Ln-Lipschitz continuous. (o(h"), ( (hh"))/ n satisfies
For any x., sd-1, the double sequence
Kobh') (((h)),((h)) (b(hT), (hT))) /n S CL"!
(15)
where C is a constant depending solely on L1. Tw. ow and T.
By combining Lemma I and Proposition 3 with Moore- Osgood theorem, as stated in Theorem g of Appendix A we can establish the NNGP correspondence for the Neural ODE fe-
Theorem 1. Suppose L Li-Lisschitz continuous 6. As width n + co, the Neural ODE fo defined in Eq. 7 converges in distribution to centered Gaussian process with covariance function Σ defined as the limit of L given in Proposition 3
7 the covariance func- (xxx, CL-1 This polynomial rate of convergence preserves the geometry of the input space (Yang & Schoenholz 2017|. This stands in contrast to classical feedforward networks, where the input space geometry often collapses unless the variance hyperparameters are set precisely on the edge of chaos (Poole el al. 2016
Remark 1. Thanks to the uniform convergence result established in Lemmal
tion L converges E' with a rate of (E'(2,
4.2 NTK: BACKPROPAGATION OF GRADIENTS
While NNGP governs the forward propagation of inputs, the NTK (Jacot et al. 2018 governs the backward propagation of gradients. Understanding both is key to comprehending the full dynamics of Neural ODEs during training. As defined for Neural ODEs in Eq. (9) we can also define the NTK
K 0 for the finite-depth network J// in Eq. (10) as follows:
K- (z. h; 0) == (oo (z; (z; 0), (2;
(16)
In the same infinite-width limit, as highlighted in previous works (Jacot et al. 2018 Yang 2020),
the NTK K3 converges deterministic kernel KL that remains constant throughout training. Notably, this deterministic limiting NTK KL (and Koo defined in Theorem 2 dynamics of Neural ODEs under gradient descent.
governs the training
Below are the results for our setup, with proofs provided in Appendix E
Proposition 4. Suppose & is L, Li-Lisschitz continuous. Then, as the network width n + oo, the NTK KJ converges almost surely to a deterministic limiting kernel: VL 2 0
K/4 (xx x) C.+11++ (x, + 2
(x, (x, (2,
(z, (1) (17)
t,k=1
where C.. are defined in Proposition and D are defined recursively:
D (x, x) = =o, (o'")'(a")))
Vk E 10.1, ,L}.
(18)
6+1,4+1
(x. x) = x? Σ
(",
Ve,ke E (,,2, I., I 1).
(19)

---

Published as a conference paper at ICLR 2025
The same problem of different convergence patterns converging different limits, observed in the NNGP kernel Σ", also arises when computing the NTK of Neural ODEs. While the Lipschitz conti- nuity of enables well-posed forward propagation of inputs in Neural ODEs, additional regularity is required for backward propagation of gradients. Specifically, Lipschitz continuity of ø is suffi cient to ensure uniform convergence in the NTK. With e and d' both being Lipschitz continuous, we can obtain a uniform convergence result similar to Lemma] |
Lemma 2. If & is Li-Lipschitz continuous and o' is Lg-Lipschitz continuous, then the following inequality holds almost surely:
Kg (r, .) - Ko( x) S CI Va, x E
(20)
where C' > 0 is a constant dependent only on the constants de, Tw, Gu. L1, L2, and T.
Combining Lemma 2 with Proposition 4 and Moore-Osgood Theorem 8. we can interchange the limits L and n in the double sequence Kg (æx, a)) and show that the NTK Kg of Neural ODE con- verges to deterministic limiting kernel.
Theorem 2. Suppose & is L1 -ipischchtzz continuous and o' is L2- Ly-Licschitz continuous. As the net- work width + 0C, the NTK Ke converges almost surely to deterministic limiting kernel:
Ke + K as n + 00,
(21)
where Kx is the limit of the NTK K defined in Proposition 4 as depth L 00.
Remark 2. Using the uniform convergence from Lemma we observe that |K (x, Koo (x, ~ CL-1, This polynomial convergence not only guarantees that gradients neither explode nor vanish L- - 0o (Yang & Schoenholz] 2017; Schoenholz et al.) [2077) but also im- plies that the limiting NTK, Kx., has an integral form, as suggested by Eq. (17 and provided in Appendix E.5) This integral form provides key insight for studying the spectral properties of the NTK Kx directly, without relying on the inductive techniques used in previous works.
5 GLOBAL CONVERGENCE ANALYSIS FOR NEURAL ODES
As discussed in Eq. (8), the dynamics of the residual r² y under gradient descent can be char- acterized using the NTK Ke. In the infinite-width limit, shown in Theorem [ this time-varying kernel Ke converges to a deterministic limiting kernel Kx.. provided the activation function á is sufficiently smooth. Therefore, in this section, we establish the global convergence of Neural ODEs under gradient descent by examining the spectral property of the NTK Kg and its limit K
The limiting NTK Kx a deterministic kernel function, and its spectral properties are key to under- standing global convergence. Previous studies (Jacot et al. 2018 Nguyen 2021) have highlighted that the strictly positive definiteness (SPD) of the NNGP kernel Σ' is suffiiennt guarantee the SPD property of Kx.. Since Σ' is a component of Koo defined in Eq. (9), demonstrating the SPD property of Σ* is critical for proving convergence.
However, prior analyses have relied on inductive proofs for finite-depth neural networks, which are not directly applicable infinite-depth and continuous networks like Neural ODEs. That is because, depth increases, information propagation can become trivial (i.c., gradients vanishing or exploding), potentially diminishing the SPD property at the infinite-depth limit (Poole et al. 2016 Schoenholz al.) 2017] Hayou & Yang 2023). Fortunately, results in Section demonstrated stable information propagation in both forward and backward directions, regardless of the choice of ou, aw. and ou. This allows to retain the SPD property of the NTK of Neural ODEs as the depth L approaches infinity.
Specifically, recall from Theorem 1 that
express Σ" as: [x .) = ([)(u)
where (u, 7) are centered Gaussian random variables with covariance S* 2) defined by
S* (a, σ) lim (0,0( L
σ) ( + K Σ cdick (,(
F,==1
with λ = T/L. This expression S* can be interpreted as a double f whose explicit form is included in Appendix E.5) By leveraging the results from Section 4 derive key properties of S* in Lemma[3] These properties serve as a fundamental basis for analyzing the SPD properties of the NNGP and NTK.
7

---

Published as a conference paper at ICLR 2025
Lemma 3. For any x, B. E Si-1, we have
(xx B) is well defined, and < St x)) (a, x) < 0c.
(xx 2 (x, 3) and the equality holds if and only if x = a.
Lemma 3 implies (z æ) = ( for all S This allows to study the SPD property of the NNGP kernel Σ" using its Hermitian expansion from the perspective dual activation (Daneely et al. 2016). Detailed analysis and proofs are provided in Appendix F Additionally, 5 (ær, x)) (" Ø) for all x * æ implies that the pathology known as the loss of imput dependence, observed in other large-depth networks such as feedforward (Poole et al. 2016). ResNet (Hayou & Yang. 2023), and RNN (Gao et al. 2023), does not occur here. This stability results from combination of several factors, including skip connections, scaling K, and smoothness and nonlinearity of ø. With stable information propagation in Neural ODEs, we can use the nonlinearity of Φ to show that the NNGP kernel Σn and the limiting NTK Ko are SPD.
Proposition 5. If é is Lipschitz, nonlinear but non-polynomial, then the NNGP kernel Σ* is SPD. Corollary 1. Suppose & and o' are Lipschitz continuous. If ④ is nonlinear but aon-polynomial, then the limiting NTK K is SPD.
With these results, we can establish the global convergence of Neural ODEs under gradient descent with appropriate assumptions about the activation function & and the training data.
Assumption 1. Let (æt, yi Vi be a training set. Assume
Training set: æ, E gd-1 and æ; + æz for all + j: (l).
Smoothness: 6 and ø' are L1- and Læ-Lipschizz continuous, respectively.
Nonlinearity: & is nonlinear and non-polynomial.
Under [ssumption I the parameters 8 remain close to their initialization 8 This proximity ensures that the Neural ODE and its gradients are well-posed not only initialization, as proved in Proposition T but also throughout the entire training. This consistency in parameter updates enables us to prove that the NTK Ke retains SPD during training, ensuring that the training errors of Neural ODEs consistently decrease to zero at a linear rate. Detailed analysis sis and proofs are provided in Appendix G
employ inductive proofs show that in the overparameterized regime,
Theorem 3. Suppose Assumption 7 holds and the learning rate Tj is chosen such that < 7 S x/ Then for any 8 >> O, there exists natural number such that for all widths n 2 ns the following results hold with probability at least 8 over random initialization Eg. 3
1. The parameters H stay in neighborhood of 0", i.e.
100h -001 S CI CI L(00)/Xo.
(22)
2. The loss function L(Ox) consistently decreases to zero at an exponential rate, i.e.,
L(Ox) S (1
17
L(Oo),
(23)
16
where do = Amin (Koo) > O, and the constant C > 0 only depends on L1, L2. ou ow ou and T.
EXPERIMENTS
In this section, validate our theoretical findings through several experiments on Neural ODEs. We focus on the approximation errors between the continuous Neural ODE and its finite-depth ResNet approximations, the NTK behavior, and the empirical convergence properties of Neural ODEs under gradient descent. Further experimental details, including additional experiments on smooth vs. non- smooth activations and scaling for long-horizon stability, be found in the appendix.
Gradient and Output Approximation by Finite-Depth ResNet. As established Section 3, Neu- ral ODE outputs and gradients be approximated by finite-depth ResNet, with error rate that
8

---

Published as a conference paper at ICLR 2025
Oucuct snfe Laylag Sooke Sett
Gradien fomece lag Sap Spaat Sohol
accommeeent
STC Labieasene
(a)
(d)
Figure 1: Analysis of Neural ODE output, gradient differences, and NTK convergence. (a) Output differences between Neural ODE and finite-depth ResNet across different widths using Softplus activation. (b) Gradient differences for Neural ODE and ResNet models under Softplus activation. (c) NTK convergence behavior across different widths, showing the NTK approximation converging to the limiting NTK as width increases. (d) NTK convergence behavior log-log scale, further emphasizing the rapid convergence at larger widths.
decays as I/L, where L is the depth of the ResNet. We empirically verify this by measuring the output and gradient differences between the continuous Neural ODE and its finite-depth approxi- mation at initialization. Both the Neural ODE and ResNet were initialized with the same random weights and evaluated on the MNIST dataset, with ResNet depths L ranging from 10 to 1,000. We used Softplus activation ensure smoothness. Figure [b)) demonstrates that the approximation error for both outputs and gradients decreases as 1/L, with convergence being uniform across dif- ferent widths, consistent with theoretical results. These findings confirm that smooth activation functions lead to well-posed ODE solutions, with accurate approximations by finite-depth networks.
NTK Approximation Error from Finite-Depth ResNet. As discussed in Section 4, the NTK of Neural ODEs can be approximated by the NTK of a finite-depth ResNet, with the approximation error decaying as 1/L. This follows from the fact that the NTK is the inner product of gradients, and as shown in Proposition 2 and Figure [a)-(b), the gradient difference between Neural ODEs and finite-depth ResNets also decays as 1/L. By applying the triangle inequality the gradient differences, it is straightforward to conclude that the NTK approximation error inherits the same I/I, decay rate. Given this reasoning and the page limit, we skip this experiment, but refer readers to Section 4 and Theorem 3 for a detailed theoretical analysis.
NTK Convergence to Deterministic Limiting NTK. In Theorem we prove that as the width of Neural ODEs tends to infinity, the NTK converges to deterministic limiting NTK. While no theoretical convergence rate is provided, we conducted experiments to empirically investigate this convergence. We evaluated Neural ODE models with increasing widths, ranging from 10 to 1, 000, and computed the NTK for each width. These NTKs were then compared to an approximate limiting NTK derived from random matrix theory. As shown in Figure [c)- the NTK converges to the limiting NTK as the width increases. The empirical convergence rate falls between 1/m and 1 1//m, with a tendency closer 1/vm when plotted on logarithmic scale. This indicates that Neural ODEs exhibit rapid convergence to their limiting NTK, validating the theoretical analysis.
NTK's SPD and Global Convergence. In Proposition and Corollary I we established that the NTK of Neural ODEs is SPD when the activation function is nonlinear but not polynomial, which guarantees global convergence under gradient descent. Specifically, the NTK's smallest eigenvalue remains positive, ensuring the well-conditioning of the model during training. Additionally, we showed that the model parameters remain close to their initial values during training, further sup- porting the global convergence claim.
To empirically verify these results, we conducted experiments with Neural ODE models of varying 5idths-500,, 1000, 2000, and while monitoring both the NTK's 's smallest eigenvalue and the distance of the model parameters from their initial values over 100 epochs. Softplus was used as the activation function to ensure smoothness and non-polynoiial nonlincarity. At each epoch, computed the smallest eigenvalue of the NTK and the Euclidean distance between the current and initial parameter values.
The Smallest Eigenvalue of NTK: As shown in Figure [aa), we observed that as the width of the Neural ODE increases, the smallest eigenvalue of the NTK becomes larger. For widths greater than the number of training samples (i.e., which is 1000 in our experiments), the smallest eigenvalue remains strictly positive throughout the training process, confirming the NTK's striet positive defi-

---

Published as a conference paper at ICLR 2025
niteness and ensuring that the model is well-conditioned for gradient descent. However, for widths smaller than the number of training samples, the smallest eigenvalue becomes negative, indicating poor conditioning at smaller widths.
Parameter Distance: The results also confirm that the parameter distance remains stable as train- ing progresses, staying within manageable bound of O (1), shown in Figure ( As the width increases, the parameter distance grows, but the growth remains stable and does not devi- significantly. This supports the theoretical result that the parameters do not stray far from their initialization, ensuring stable training and global convergence.
Train and Test Loss: Finally, as depicted in Figure (c)- we observed that larger widths lead to faster convergence of the gradient descent. Models with larger widths (2000 and 4000) exhibited lower test losses and faster convergence rates compared to smaller widths (500 and 1000). This behavior demonstrates that larger widths allow the model to generalize better and converge more efficiently during training.
Smaliie Ogeersate by Valris
Persmene Deeance by Witth
hai Laoas lag as Nuts
Las 1NC 280
Naciona: was L00
The wiim 1962
Taltten Wiln (a0n
1818
Newm Detanee toen 1890 maxi Luitanos toch +94) Cacht
(a)
(b)
(c)
(d)
Figure 2: Empirical results of Neural ODEs with varying widths: (a) NTK smallest eigenvalue grows and stabilizes the width increases, with negative values for widths below the training size. (b) Parameter distances stay stable and bounded within O (1). (c) Licle train and test losses show faster convergence for larger widths. (d) Log-scale losses further confirm improved generalization for wider models.
In the appendix, we present supplementary experiments that validate and extend our findings. Without proper scaling (c.g., Tw ~ 1/T). Neural ODEs exhibit carly-stage damping during training over long-time horizons (see Figure 3.. Smooth activations like Softplus converge faster than non-smooth ones like ReLU, likely due to more accurate gradient computation (see Figure 77. Additionally, while non-polynomial nonlinearity is sufficient for an SPD NTK, experiments show that quadratic activations also yield SPD NTKs, though with slower convergence (see Figure 88. These results highlight the importance of activation functions and model design for Neural ODE performance. We also include convergence analysis diverse datasets, such as CIFAR-10, AG News. and Daily Climate, as well as additional activations like GELU, further demonstrating the generalizability of our findings.
Additional Experimental Results.
7 CONCLUSIONS
In this paper, we examined the crucial role of activation functions in the training dynamics of Neu- ral ODEs. Our findings demonstrate that the choice of activation function significantly impacts the dynamics, stability, and global convergence of the Neural ODE models under gradient descent. Specifically, we found that using smooth activations like Softplus ensures that the forward and back- ward dynamics in Neural ODEs are well-posed, allowing for accurate approximation by finite-depth ResNets. As a result, the NTK of Neural ODEs converges to deterministic limiting NTK that gov- emn the model's training dynamics. Additionally, demonstrated that when using nonlinear but non-polynomial activations, the NTK remains SPD, ensuring well-conditioned training and global convergence. Through extensive experiments, verified that suitable activation functions, Neural ODEs exhibit stable parameter behavior, rapid NTK convergence, and faster optimization, particu- larly at larger widths. These findings highlight the importance of selecting activation functions with appropriate smoothness and nonlinearity to ensure the robustness and scalability of Neural ODEs, establishing them as a powerful approach for continuous time deep learning.
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song. A convergence theory for deep learning via over- parameterization. International conference on machine learning, pp. 242-252 PMLR, 2019.
Sanjeev Arora, Simon S. Du, Wei Hu, Zhiyuan Li, Ruslan Salakhutdinov, and Ruosong Wang. On exact computation with an infinitely wide neural net. In Thirty-third Conference on Neural Information Processing Systems, 2019.
Zhi-Dong Bai and Yong- Qua Yin. Limit of the smallest eigenvalue of large dimensional sample covariance matrix. In Advances In Statistics. World Scientific, 2008.
Jinghui Chen, Dongruo Zhou, Yiqi Tang, Ziyan Yang, Yuan Cao, and Quanquan Gu. Closing the generalization gap adaptive gradient methods training deep neural networks. arXiv preprint arXiv:1806.06777 2018a.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 201 8b.
Amit Daniely, Roy Frostig, and Yoram Singer. Toward deeper understanding of neural networks: The power of initialization and a dual view on expressivity. Advances in neural information processing systems, 2016.
Jianzhun Du, Joseph Futoma, and Finale Doshi-Velez Model-based reinforcement learning for semi-markov decision processes with neural odes. Advances in Neural Information Processing Systems, 33:1 19805-19816,, 2020.
Simon Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai Gradient descent finds global minima of deep neural networks. International conference on machine learning, pp. 1675- 1685. PMLR, 2019a.
Simon Du, Kangeheng Hou, Russ R Salakhutdinov, Barnabas Poczos, Ruosong Wang, and Keyulu Xu. Graph neural tangent kernel: Fusing graph neural networks with graph kemnlss. Advances in neural information processing systems, 32, 2019b.
Chris Finlay, neural ode: the world of jacobian and kinetic regularization. In International conference on machine learning, pp. 3164-3164. PMLR, 2020.
Jacobsen, Levon Nurbekyan, and Adam Oberman. How to train your
Tianxiang Gao. Mastering infinite depths: Optimization and generalization in deeper neural net- works. PhD thesis, Iowa State University, 2024.
Tianxiang Gao and Hongyang Gao. Gradient descent optimizes infinite-depth relu implicit networks with linear widths. arXiv preprint arXiv:2205.07-663 2022a.
Tianxiang Gao and Hongyang Gao. On the optimization and generalization of overparameterizze implicit neural networks. arXiv preprint arXiv:2209. 2022b.
Tianxiang Gao, Hailiang Liu, Jia Liu, Hridesh Rajan, and Hongyang Gao. A global convergence theory for deep relu implicit networks via over-paramettrization.. International Conference on Learning Representations, 2021.
Tianxiang Gao, Xiaokai Huo, Hailiang Liu, and Hongyang Gao. Wide neural networks as gaussian processes: Lessons from deep equilibrium models. Advances in Neural Information Processing Systems, 18:54918-549511 2023.
Amir Gholami, Kurt Keutzer, and George Biros. Anode: unconditionally accurate memory-efficient gradients for neural odes. In Proceedings of the 28th International Joint Conference Artificial Intelligence, 2019.
Amir Gholaminejad, Kurt Keutzer, and George Biros. Anode: Unconditionally accurate memory- efficient gradients for neural odes. In International Joint Conferences on Artificial Intelligence, 2019.
11

---

Published as a conference paper at ICLR 2025
Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp. 249-256. JMLR Workshop and Conference Proceedings, 2010.
Tilmann Gneiting. Strictly and non-strictly positive definite functions spheres. Bernoulli, 2013.
Will Grathwohl, Ricky TQ Chen, Jesse Bettencourt, Ilya Sutskever, and David Duvenaud. Ffjord: Free-form continuous dynamics for scalable reversible generative models. arXiv preprint arXiv:1810.013377 2018.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. Advances in neural information processing systems, 32, 2019.
Soufiane Hayou and Greg Yang. Width and depth limits commute residual networks. Interna- tional Conference on Machine Learning, pp. 12700-12722 PMLR, 2023.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In Proceedings of the IEEE international conference on computer vision, pp. 1026- 2015.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp- 770-778, 201
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Identity mappings in deep residual net- works. In Computer Vision-ECCY 2016: 14th European Conference, Amsterdam, The Nether- lands, October 11-14, 2016, Proceedings, Part IV 14, pp. 630-645. Springer, 2016b.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840-683 6851, 2020.
Hochreiter. Long short-term memory. Neural Computation MIT-Press, 1997.
Jiri Hron, Yasaman Bahri, Jascha Sohl-Dickstein, and Roman Novak. Infinite attention: Nngp and ntk for deep attention networks. In International Conference Machine Learning, pp. 4376- 4386. PMLR, 2020.
Sergey loffe. Batch normalization: Accelerating deep network training by reducing internal covari- ate shift. arXiv preprint arXiv:150..03167, 2015.
Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kemnl: Convergence and gen- cralization in neural networks. Advances in neural information processing systems, 2018.
George Em Karniadakis, Ioannis Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang- Physics-informed machine learning. Nature Reviews Physics, 3(6):422-44) 440, 2021.
Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equa- tions for irregular time series. Advances in Neural Information Processing Systems, 33:6696- 6707, 2020.
Joon-Hyuk Ko, Hankyul Koh, Nojun Park, and Wonho Jhe. Homotopy-based training of neuralodes for accurate dynamics discovery. Advances in Neural Information Processing Systems, 36:64725- 64752, 2023.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86( 86(11):2278-2324, 1998.
Jachoon Lee, Jascha Sohl-dickstein, Jeffrey Pennington, Roman Novak, Sam Schoenholz, and
Yasaman Bahri. Deep neural networks gaussian processes. In International Confer-
ence on Learning Representations, 2018. URL httpps://openeeveew.... htwpppewwwwwwwweeeeeeeeeeeeeeeee net/forum?id
B 11EA-M-Z
Xuechen Li. Ting- Kam Leonard Wong, Ricky TQ Chen, and David Duvenaud. Scalable gradients for stochastic differential equations. In International Conference on Artificial Intelligence and Statisties, pp. 3880-3882, PMLR, 2020.
12

---

Published as a conference paper at ICLR 2025
Stefano Massaroli, Michael Poli, Jinkyoo Park, Atsushi Yamashita, and Hajime Asama. Dissecting neural odes. Advances in Neural Information Processing Systems,
2020.
Takashi Matsubara, Yuto Miyatake, and Takaharu Yaguchi. Symplectic adjoint method for exact gra- dient of neural ode with minimal memory. Advances in Neural Information Processing Systems, 14:20772-20784 2021.
Radford M Neal. Bayesian learning for neural networks, volume 118. Springer Science & Business Media, 2012.
Quynh Nguyen. On the proof of global convergence of gradient descent for deep relu networks with linear widths. In International Conference on Machine Learning. pp.
PMLR, 2021.
Derek Onken, Samy Wu Fung, Xingjian Li, and Lars Ruthotto. Ot-flow: Fast and accurate continu- normalizing flows via optimal transport. In Proceedings of the AAAI Conference Artificial Intelligence, volume 35, pp. 2223-923,, 2021.
Katharina Ott, Prateek Katiyar, Philipp Hennig, and Michael Tiemann. Resnet after all: Neural odes and their numerical solution. In International Conference Learning Representations, 2020.
Sunghyun Park, Kangyeol Kim, Junsoo Lee, Jaegul Choo. Joonscok Lee, Sookyung Kim, and Ed- ward Choi. Vid-ode: Continuous-timg video generation with neural ordinary differential equa- tion. In Proceedings of the AAAI Conference Artificial Intelligence, volume 35, pp. 2412-2422, 2021.
Ben Poole, Subhaneil Lahiri, Maithra Raghu, Jascha Sohl-Dickstein, and Surya Ganguli. Exponen- tial expressivity in deep neural networks through transient chaos. Advances in neural information processing systems, 29, 2016.
Maziar Raissi, Paris Perdikaris, and George Karniadakis. Physics-informes neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Ivan Dario Jimenez Rodriguez, Aaron Ames, and Yisong Yue. Lyanet: A Iyapunov framework for training neural odes. In International Conference on Machine Learning, pp. 18687-18703. PMLR, 2022.
Yulia Rubanova, Ricky TQ Chen, and David K Duvenaud. Latent ordinary differential equations for irreuularyy-amppled time series. Advances in neural information processing systems, 32, 2019.
Michael Sander, Pierre Ablin, and Gabriel Peyré. Do residual neural networks discretize neural ordinary differential equations? Advances in Neural Information Processing Systems, 2022.
Samuel S. Schoenholz, Justin Gilmer, Surya Ganguli, and Jascha Sohl-Dickstein. Deep information propagation. In International Conference Learning Representations, 2017. URL https:
7/ //openreview. wttp:/rwwwwww..e/ffoooooee
W1 UN 9qg
Aleksei Sholokhov, Yuying Liu, Hassan Mansour, and Saleh Nabi. Physics-informed neural ode (pinode): embedding physies into models using collocation points. Scientific Reports, 13(1): 10166, 2023.
Justin Sirignano and Konstantinos Spiliopoulos. Dgm: A deep learning algorithm for solving partial differential equations. Journal of computational physics, 375:1339-1364, 2018.
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In Interna- tional Conference on Learning Representations, 2020.
A Vaswani. Attention is all you need. Advances in Neural Information Processing Systems, 2017.
Roman Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv preprint 101 2010.
Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449:1 10768, 2022.
13

---

Published as a conference paper at ICLR 2025
Ge Yang and Samuel Schoenholz. Mean field residual networks: On the edge of chaos. Advances in neural information processing systems, 30, 2017.
Greg Yang. Wide feedforward or recurrent neural networks of any architecture are gaussian pro- cesses. Advances in Neural Information Processing Systems, 2019.
Greg Yang. Tensor programs ii: Neural tangent kernel for any architecture. arXiv preprint arXiv:2006. 14548, 2020.
Greg Yang. Dingli Yu, Chen Zhu, and Soufiane Hayou. Tensor programs VI: Feature learning in infinite depth neural networks. In The Twelfth International Conference on Learning Representa- fions, 2024. URL http::://ppeeeeeeeee wttp:///oooruum///111//pppppp
Juntang Zhuang, Nicha Dvonek, Xiaoxiao Li, Sekhar Tatikonda, Xenophon Papademetris, and James Duncan, Adaptive checkpoint adjoint method for gradient estimation in neural ode. In International Conference Machine Learning, pp. 11639-11649, PMLR, 2020a.
Juntang Zhuang, Nicha C Dvornek, James Duncan, et Mali: memory efficient and reverse accurate integrator for neural odes. In International Conference on Learning Representations, 2020b.
Difan Zou, Yuan Cao, Dongruo Zhou, and Quanquan Gu. Gradient descent optimizes over- parameterized deep relu networks. Machine learning, 109:467-492, 2020.
14

---

Published as a conference paper at ICLR 2025
A USEFUL MATHEMATICAL RESULTS
Theorem (Bai-Yin law, see Vershynin 2010): Bai & Yin 2008)) Let be an x n random ma- trix whose entries and independent copies of random variable with mean, unit variance, and finite fourth moment. Suppose that N and n grow to infinity while the aspect ratio n/N converges constant in [0, 1]. Then
Smin ( ()) = / - vn + o Smax VN + vn + (Vv), almost surely.
(24)
theorem). Let f (a, bj x R + R be a function. If f is contimoous in the first argument and Lipschitz continuous with coefficient in the second argument, then the ODE
Theorem
x(1) - f(t,
(25)
possesses unique solution on (a - E, + E) for each possible initial value æ<a) = IO E R", where E < L/L.
Theorem 6 (Peano Existence Theorem). If the function is continuous in neighborhood of (to, zo). then the ODE equation 23 has at least one solution defined in a neighborhood of to.
Theorem (Convergence for Euler's Method). Let In be the result of applying Euler's method to the ordinary differential equation defined as follows
i = f(x, t), E ti), and x(0) = zo.
(26)
the solution has bounded second derivative and f is L-Lipschitz continuous in x, then the global truncation error is bounded by
hM
|x((t) =xall S
(c -1),
(27)
2L
where h is the time step, and M is upper bound on the second derivative of x on the given interval.
Lemma 4 (Gronwall's inequality). Let T= (a, b) for interval such that a < b co. Let u, a, be real-valued continuous functions that satisfies the integral inequality
u(t) S a(t) + 1 B(s)u(s)ds, 1
(28)
Then
u(t) S a(1++ f a()))(s) exp B())a)) V.ET.
(29)
If. in addition, all) is non-decreusing, then
u(t) S a(t) exp (f, ((3)ds) Vte.. I.
Theorem Moore-Osgood Theorem). If lim anm by. uniformly in m, and lim On.m = Cn n-+00 for each n, then both lim bi and lim Cn exists and equal the double limit, i.e.,
m-
lim ( lim an,m) = lim ( lim an,m) lim On.m 1-000
(30)
mo 400
n0000 m00
B DERIVATION OF GRADIENT THROUGH ADJOINT METHOD
In this section, we provide the detailed derivation of the adjoint method to compute the gradients. We first recall the forward ODE follows:
h,
Ow Wo(he),
Vn
Ou
hn
Ux.
Vd
15

---

Published as a conference paper at ICLR 2025
To compute the gradients, introduce the Lagrange function
((h,0,,,))= fo(z)+ + ['"
Fw (e W ((h) - i) dt + nT (v h(0))]
folx)+
x
Vn
where h is extra variable that are independent from 0 and (2, Lagrangian multipliers. Observe that with h we have
0(,,0, fo(x). v)..
Thus, the derivatives of L w.r.t. is equal to gradients of fo w.r.t 6.
Now, we consider a variation (8h, 50) at point (h, 0). Then the correspondence variation of L is given by
8L 0, = Jo ((h))((((TT)) + dv u' (( + + T vn
Ju (6U)x - -8h(0)
vn
Vd
A
(BW))(h)
Tw W Wdiag ( (h))ōh 8h dl
Vn
Vn
ov TT (lag(d((((T)))) (ov) ((h(T))
+ nr
(6U) -8h())))(()))
vn
vn
Vd
Toho + !'" TT + ' -
A
(6W))())) + Vn
ow
au
W Wiag((o diag(d((h))))
Vn
a (80) To((T))
+ ao
dida(d(' (d(()) A(T)7 5h((T
nn
Vn
T
Ou (BU) + (A(0) 14) 81((0)
+ 1
+
Ow Vn"
x WWdiagg (hi)) öhdt + S0
Vn" x (oWW)o())d
where we use integration by parts in the second equality. Then we choose (A, 11) such that
/1 =(()),
dv Vn
X(T)=
diag( (o' (h(T)))u,
X()) vn
diag( ((((())) X(I.
Then the variation of & becomes
5c(h, 0, A, = d.
(8U)xx /'' aw ow ((W()))((h))d. ou vill ((((T))
Vn
vn"
Thus, we obtain the gradients of fo as
vofo(z) vv
o(h(T))
Vn"
Vwff(()) = ['' ow
x T
vn
VuJe(w) Tu 100)xT.
C WELL POSEDNESS OF NEURAL ODES AND ITS GRADIENTS
To show the existence and uniqueness, first recall the Picard-Lindelit theorem follows.
16

---

Published as a conference paper at ICLR 2025
C.1 FORWARD ODE IS WELL-POSED
As we assume the activation function is Lipschitz continuous, we can immediately obtain the local result that the hidden state hy exists near the initial time. Lemma (Local solution). If the outictotonfuunctionn & is L Liiisschitzz continuous, then hú uniquely exists for all |t| S E, where ε <1/owl1.
we know || ~ v a.s. Accordingly, we can show the mapping f: x + to W ((x) is Lipschitz continuous:
Proof. By using Bai-Yin law 4
1 -
=1 =1 @oWll Φ(x) Vn
Tu W Wo(z)]l
Vn
Sa (owllo(x) - x(=)]l
5a L1 || =||-
Hence is ow L--iishhst continuous a.s. As to = 0, it follows from Picard-Lindelöf theorem that unique h2 exists locally for all [t] S e, where e < 1/ow,11
Lemma (Global solution). For any given T > if o is L1 Lipschhtt continuous, then h2 uniquely exists for all itl < T.
Proof. We have shown unique ht exists locally. Specifically, let (x be the solution flow from initial condition to the solution at I. For any ho. chose 1/. Then the solution h1 == oz (zo) is well-defined based on the local solution result. As the dynamics is the same and the Lipschitz coefficient is uniform, we have h2 (h1) is also well-defined. By repeating this process for any finite steps N, we have ly = @e is well-defined. Hence, as T < oo, there exist N such that eN 2 T. Therefore, h, is well-defined for all [t] T and the desired result is obtained.
Then the result for global solution simply implies that result in Proposition
C.2 BACKWARD ODE IS WELL POSED
Recall the backward ODE as follows
\T diag diag( (bbTr,,
(31)
Vn
A, Tw diag( (ho))) A,-
(32)
Vn
Observe that if hy is well defined in t E [0, then the dynamics of Xe becomes a linear dynamics. Hence, with similar argument, we can easily show the corresponding VIP of X is well posed. Lemma 7. Given T, if the activation function ó L--isscciizz continuous, then de is uniquely determined for all [t) T and X, = / is the solution.
Proof. It follows Lemma and that h, is well defined for all t E [0, a.s. By Theorem 5 suffices to show g x vn diag ø' ((h.)) is Lipschitz continuous:
(1a(x) - -g(g=](===10- = -" diaglo' (ho)) (x (x)] S Lw z||,
Vn
where we use the fact [W]| ~ /n a.s. by Theorem 4 and |o'l < L1. Hence, the mapping g is Ly Lipschitz continuous. follows from Theorem that X; uniquely exist for E [T - T' + for < 1/o,,Ii Then with similar argument, we can show the existence of local solution can be extended to global solution since ò is uniformly Lipschitz continuous. Therefore, X is well defined for all [O.
Additionally, we can show A(t)= (ff(()) is solution. Specifically, the differential of fo given by ah(t)
dfo du do (((((T))
diag (((((
Vn
17

---

Published as a conference paper at ICLR 2025
Then we have Moreover, for any > it follows the chain rule that
fa(2) 1 diag( (h(T)))
diag(o((((T)))
(33)
Oh(T)) vn
@f((r) ah(t +) 0f)(z) ah(t)
ah(t)
Sh(t +E)
where we have
t+5 1 Wo(h(s))ds.
h(t +) hh()) +
(34)
Then we have
afa Bb((+c)
f (h(())
d
(ff((a)
lim 6-+01
dt
@h(t)
8
0fa Ob((+e)
8h(+++) a
ah(() ah(++e)
lim
6-+0+
afa 8h((+e) Bh(t)
(h(t) + tu +o(h())==+(() Of ah(t+e)
lim
5-+01
8fa (b((+e)
(1+
vn diag( diag( +0 +0((2)
81e (N([+s)
lim
diag(ø' ( Ofo
Vn
Qh(t)*
Thus, we can see 0ff(z)/Oh() is the solution to backward ODEs.
C.3 AUGMENTED BACKWARD ODE WELL POSED UNDER THE SAME REGULARITY
It can be seen from equation 5 that to compute gradient of fe w.r.t. W. need to evaluate values of ht and X for every t E [0. 7- However, Chen et al. (0188) suggested to solve an augmented backward ODE. As result, there is no need to store the intermediate values of h and )
We first recall the gradients of fo w.r.t. 0 in a vectorization form:
0.fe(z) de o(h(T))
(35a)
Vn
{" Tw
Owfe(z)
(6(h,) & Ax)dt
(35b)
Vn
fo Vd ((0))
(35c)
The according augmented backward ODE is given by
he
Wo(h)) diaglo' (h) TT (h,,) 3 Ar
Ow
Vt E [0,2]
(36)
Vn
[
where g: _=*xx and the initial condition hr and AT combined with gr
Once this augmented backward ODE is solved, the gradients of fe(x) w.r.t. W can be obtained by
V w fox) =0
=g(T) + f gidt
=9(T) + f
0
- vn
Ow
((h)) @ Aedt
T
(6(h.) & Xe] dt, = =1" ow
In
18

---

Published as a conference paper at ICLR 2025
where the fact 9r 0. Unlike h, is known in equation 4 h4 is unknown state in the augmented backward ODE equation 36 Hence, it follows from Theorem that extra smoothness is generally required to ensure the well posedness such ø' probably need to be Lipschitz continuous. However, the dynamics of h4 is decoupled from the dynamics of X and Gt- Hence, one can solve he first (in the backward manner), then solve the dynamics system for X and Jt- In this manner, he is still an known states. Hence, one can use the same regularity condition in Proposition to show existence of unique solutions for & [0, 27. Therefore, additional smoothness needed solve the augmented backward ODE.
D NNGP CORRESPONDENCE FOR NEURAL ODES
In this section, we establish the NNGP correspondence for Neural ODEs. It follows from the Euler method that Neural ODE can be approximated by finite-depth neural network fá equation 10 From the asymptotic perspective, Neural ODEs is equivalent to an infinite-depth ResNet with shared parameters in its all hidden layers and a special depth-dependent scaling hyperparameter T/ /T.
D.1 FINITE-DEPTH NEURAL NETWORKS AS GAUSSIAN PROCESSES
As the finite-depth neural network If can be considered an approximation to the Neural ODE Se. first study its signal propagation by establishing the NNGP correspondence for f/.
We define vectors gf E R"
gd (z) = Uz, o (x) T W (h ((hh11),
(37)
... Vl E ... ...
(38)
The vectors gf are G-vars in Tensor program Yang 2019). Tensor program is an representation of the neural network computations that only involves linear and element wise nonlinear operations. In the Yang 2019), the authors claim that computation using G-vars is equivalent to another com- putation that using corresponding list of oneediimnsional Gaussian variables in the infinite-width limit, long the computation only involves controllable nonlinear functions. The corresponding definitions and Theorems are reformulated as follows.
Definition 1. (Yang 2019 Simplified version of Definition 5.3) real-valued functionn vp :
R
R
called controllable if there exists some absolute constants C, c > O such that [+(x)] S
Ce Eizilel
Theorem 9. (Yang 2019 Theorem 5.4) Consider a NETSOR program that has forward computation for a given finite- neural network. Suppose the Gaussian random initialization and controllable activation functions for the given neural network. For any controllable sp RMM -> IR, as width + oo, any finite collection of G-vars Ia with size M satisfies
Σ E. 2al a3 Es((2", ... M).
(39)
where {=0.... .... - are Gaussian random variables whose mean and covariance are computed by the corresponding NETSOR Program
Notably, controllable functions are not necessarily smooth, although smooth functions can be eas- ily shown to be controllable. Moreover, controllable functions. as defined in (Yang) 2019 Defi- nition 5.3), can grow faster than exponential but remain I,¹ and i.2-interrable with respect to the Gaussian measure. However, the simplified definition presented here encompasses almost most functions encountered in practice. Moreover, the vectors of or G-vars are not necessary to encode the same input Φ. Hence, g (x) and g are two different G-vars in Tensor program. However, Theorem still holds for any finite collection of G-vars, even they have the different input encoded. Therefore, by utilizing Theorem 9 we can show as Gaussian Process weakly and the result is stated in Proposition 3 and the associated Tensor program for f/ is provided in Algorithm
the finite-depth network f/ tends to
19

---

Published as a conference paper at ICLR 2025
Algorithm ResNet of Forward Computation on Input x
Input: Ux/Vd: G(n)
Input: W n(n,,)
Input: v: G(n)
1: 10 = Un/Vdl vd: G(nz)
2: for E [L] do
3: x' o(h') H(n)
4: g = Wx //V G(n)
5: h ht-1 +..g. g. G(n)
6: end for
7: xt o(h²) H(n)
Output: o' 2t///0 /vn
In the rest of this subsection, we will provide rigours proof to show the NNGP correspondence for fg through induction. For simplicity, the proof assume only input is given, while the result
for multiple inputs similar. Additionally, we also assume a, = ow = oa = since their values not significant in the proof long their values strictly positive.
BASIC CASE L =0
As L O, we have (z) x" ((h))//)/))) Hence, we don't have the hidden layers. Based on the random initialization equation we have
Lid
9k
(x..))
Let 8 be the smallest o-algebra generated by g.. By condition on 8 we have
18 ~ NO. 100//))))
where 60 = =) It follows from the law of large that
0² Σ (0(h))²
o? 12
- 101
n Σ ((9h)) o. Ed(2")² = (x,,x)
where 00 ~ (0, (x,, x)). As the limit is deterministic, the conditional and unconditional distri- bution converge to the same limit. Therefore, have
08 + (P(O,E')
where
Σ (x, x) = ( (()
where we use 20 ~ O to denote centered Gaussian random variable(s) whose (co)variances can be computed using covariance function n
GENERAL CASE I.
Now consider f = =TTohh)//m Here we have h² = h1-1 + 8g and g" Wo(h'-1), where 3 T As W is used before, let BL-1 be the smallest o-algebra generated by (g0 ... g6 1-1. Then have
v 2 ... ,.-.....111111 -- g = Wo(h
equivalently
g
W 0
where G (L-1) and + (L-
20

---

Published as a conference paper at ICLR 2025
We can obtain the conditional distribution of W by solving the following optimization problem
min W 2
W s.t. G = We
The Lagrange function is given by
1
L(
/ | P + (V.G W
Then
VWL(W,V) - VV = W =V²
As G W we have
G W=
- V = G(D")) W = G(D""))0'
Thus, we have
WB = W + WII" = -( // + W (In - (
where II = In - 1 W is c.i. of W, and o' ( ) 4z
Since g" Wo(hh-1), we have the conditional distribution of of as follows
NG ( ( "T== o, ||T /n). where Gk. denotes the k-th row of matrix G and e  for simplicity.
gk
independent
As Lipschitz continuous activation is controllable function, it follows from Theorem and inductive hypothesis that
(0' /n- p(h))(()))
(6 + Bgt +... Bok) - BBg + .... 85%)
0 Ex(2 + 32 + (2'))(2) + B2² +... B22)
=')(('))((u),,
where we define another Gaussian random variable w' to simplify the notation: w' = 20 Bz +... + 3z
Therefore, we have
x /n a Ep("')((w)) (4 0)////n ( o) / # Es((''d(uttt)
(tTT),//n = (
For 8 E (0,1, ... L - 1}, let U = fu', ,u'} be collection of w'. We define UKk [(++1)((++1] as
I(U Ubj E(u', u²) Ep((')d(()) Wi 10,1, O.J E 10,1,... .k].
Therefore, we have
(6' ($ 0/n1 (0*)/n) + E(UL22 Q(U--200--00 ( wt1)
Moreover, observe that
1
| /1
2 (I, ΦΦi)Φ
1
0" 6 " "0"((""0'0"""
n
n
=0 T h/n n (0 (
~E(u 4-1 n -1) I(uh-) 1 2 21--,t (nb-1)
21

---

Published as a conference paper at ICLR 2025
Therefore, for any controllable function d/, it follows from Theorem that
1
(gk. gh. ... .k) EE [1(2 20, .. ...,,,)))) [2]]
n 1=1
where
(or(="(z), 0, ve21 (ov(2((x), = E oo (u (x)) @ (u x (a))
ve,k 1
Let B be the smallest o-algebra generated by {0o. ... gL By condition on B", we have (x))BB NN 110
(40)
where
1
1
)hp
71
k=1
I
2
+
3
9k
k=1
I
E
②
0 +
8
= Eld(un)) = =--+=- (x,
1=1
Thus, we obtain
f/ -P
where
(x, (x, R) = =(((())))(( o (u
D.2 NEURAL ODES AS GAUSSIAN PROCESSES
In this subsection, we prove Neural ODEs tends to a Gaussian process as the width - oo. As the output parameter is independent from all previous weights, by conditioning on the previous hidden layers, the Neural ODEs becomes Gaussian random variable with covariance || (x) i.e,
[o(z) N
where we denote (T(xx) := ((hr(x)) to simplify the notation, hT is the exact solution from the forward ODE, and B is the smallest olalebrra generated by previous hidden layers. Here we also assume og = aw ou = a as their values are not important in the proof as long as they are strictly positive.
follows from convergence analysis of Euler's method, stated in Theorem that (x(x) + ør (x), as - 20,
where we denote (x()) = ((h
Thus, the focus analysis becomes to study the convergence of this double sequence
On := (0*(z), /n.
By leveraging the convergence result for Euler's method Theorem 7 we can show the double sequence Ln converges as L + oo and this convergence is uniform in n a.s.
Lemma 8. o is Lipischhit continuous, then the following inequalities hold for every x E Sd-1 a.s.:
||hal| S Cole
VI E (0,T)
(41)
and
A
U - -1--0 S
(e Bta 1) IVV.
(42)
2B
where A = CoLIT and B = Ca for some absolute constant C 0.
22

---

Published as a conference paper at ICLR 2025
Proof. Recall from Lemma 5 that the mapping f 1 H Vn Wo(2) is - continuous.
Observe that
d(h) do (((())/
W Wdiag (o' (h(t))) dh(tt
W Wdiag (6' (h(t)) h(t)dt
Wdiag (h(t)) dt.
Then we have
d
h
-
W Wdiag (d(()(()) W (o(h()) Vn
dl
and
|ill S c L
where we the fact ||]]| CVn from Theorem for some absolute constant C > 0 and is L1 Lipschitz continuous.
Then
M(1) 1(() + f(hate hnds
implies
[л(!II ((1A0)) / 11 - Wo(h(s))l
Vn
+|1+00)+++ ( CoLillb(s)lds.
By using the Gronwall's inequality, we have
[パ(())| < ||<(O)] exp Coladd) = [|((0)]le
Additionally, as ||0|| C.n almost surely and ||>||| 1, have ||А(O)|| - CVn,, and obtain
[I(tt]] & CoLi Vt E [0,
Therefore, we obtain
[M(0)) S Co Li Colst 1I [o,
By the Euler's convergence theorem stated in Theorem
we have
A
T
|| h(te)ll 2B where A = Co Lo2 and B = CaLi.
(e Bt
1)
L Vn.
Lemma 9, Suppose L1 Lipschitt continuous activation à and he(x)) is the exact solution with input
x. Given L, we have
1 (()("())).,(h(cc - (o(t,, ((h,, (+))) SGL-, Vk, l [L] n
(43)
where t = k3 and C, > some constant that does not dependent on n and L. Therefore, the double sequence (0(h(()) ((h((())) /n converges w.r.t. L. and uniformly w.r.t. n.
23

---

Published as a conference paper at ICLR 2025
Proof. For simplicity, assume the activation function is 1-Lipschitz continuous, i.e., L1 1. For S I L. we denote of = (x)) ø' = =ohh' (ø)), o(t) = ((he(z)), and 6(t) = (hh
where (y(x)) is the exact solution the ordinary differential equation that encodes imput z. Then consider
- -0-,0--3(13)) 1 - -
(ot ;0¹) //nnnnnn(k300 (o(k()())
where B T/L is the time step.
Note that
| = |1 + T w (o0()lll s |141|
T Co 11 (1 CoTyL)II'1
n
where use the fact that Φ is 1-Lipschitz continuous and |W]| CVn a.s. Repeat this argument times and we have
|| (1 (oT
Therefore, we obtain
1 [+(())] - (1 CoT/L) CoTT/LL ( CoTel/
where we also use ||U|| < C Ti a.s. and | =
Moreover, we have
((BB)] < // - h(E3)) S Ci CiviL
where C > 0 a constant that does not dependent on and L.
Therefore, we obtain
1 / (0(33))))) /n| CIvR - C
Hence, (6(h((x)) ((h /n converges w.r.t. L and uniformly in
Combining Lemma with Moore-Osger theorem, stated in Theorem 8 the double sequence an,L. = (0(h" (z)), (h"(()) /n has both iterated limits that are equal to the double limit, i.e.
lim (6(h((z)), (br (i))) /n
lim lim (0(1(x)) (1))))(0(100 (+))) /n
a L
lim lim (0(h"(x)), ((h²(a)) /n
L n-000
lim E (x,T)
-x
As Σ" is a deterministic function, the conditioned and unconditioned distributions of fe(x) are equal in the limit: they are centered Gaussian random variables with covariance E' 22(22,2) 2).. This complete the proof of Theorem
E NTK FOR NEURAL ODE
In this section, derive the neural tangent kernel (NTK) for Neural ODEs and provide sufficient
condition to show when the NTK well defined for Neural ODEs. Under our exploration, the smoothness of the activation function play significant role to study the NTK of Neural ODEs. For example, additionally smoothness is required to ensure the uniqueness and existence of the adjoint state D, in the backward ODE equation or augmented backward ODE equation 36
24

---

Published as a conference paper at ICLR 2025
E.1 CONVERGENCE ANALYSIS EULER'S METHOD FOR BACKWARD ODE
Similar to the forward ODE, we can also discretize the backward ODE as follows:
x X -
Ow
diaglø' (hw,)) ( AI E [1,2, .. ,
(44)
vn
where B = T and he is the solution from the forward ODE equation 2and b == Be. Additionally, we can further discretize h4 and obtain
X x+1 X 3 ow diag(ø' (a ( v 11,, L].
(45)
vn
As L 0o or CC, we have h he, and X + At.- By utilizzing the result from convergence analysis of Euler's method equation| 7 we obtain the convergence rate, that indicates this conver- gence is uniform in width n, if the activation function is smooth. This result serves as fundamental result to ensure the NTK for Neural ODE is well defined and allow us to study the training dynamics of Neural ODEs under gradient-based methods.
Lemma 10. If o and o' are L1- and L2 Li-isccchtzz continuous, then the following inequalities hold for every Sd-1 a.s.:
|| CaLie (OLi(T)) [0
(46)
and
1 Xx1l S T,
_ -
C1
(47)
C2
where = Co0 C2 Ca In CoLn for some absolute constant >
0.
Proof.
For the mapping f (A,t) + -i (1((h))))))))) WTT
TM X, we consider
di =d
diag (o' (h(t)) W >
d +'(h(t)) @ w²>)
[de'(he)) @ TV
-di(h) dh o WT>
o²(h.) o W hall
(h.) O WT @ We(hi)dl
diag (ø" (h+) diag (WTT>) W (ha)ddt
where @ denotes element-wise product and we denote W = oW/ vū. Thus, we have
(A, -diag (he) diag
(TT)) W Wo(h.)
Let Wy be the k-th column of W. As in this case consider À as fixed, w. follows Gaussian distribution with zero mean and variance a² o211A11//n. We obtain the inequality
||0. 4)I < |0"1 S/ ||||| [| [1o(h))] s C2, |/.al//vm.
where we use the assumption |Ø"| Le and C' > 0 is some absolute constant. Observe that
| ll < PArll + J. |illts s CoLi + S. CoLi
Then it follows from the Gronwall's inequality that
|all < CaL, exp ( [CCoLdds) < COL
25

---

Published as a conference paper at ICLR 2025
Combining the above bound of d with equation 41 we have
[0,f(A, f)ll - CLL COLTT == C1
Note that Ci > 0 is independent from L and n.
With argument alike Theorem 7 we can obtain the global truncation error for X In Proposition and Proposition 2 we have shown the uniqueness and existence of hy and X for all f [0, To study the convergence of A to Ate. it is equivalent to apply Euler's s method to numerically solve A in the reverse order from truncation errors for X
é T. Hence, we will assume Ag is know and provide the global
Note that
||AC+1 x(te+1)ll - (4")) (xcte) 3r(te) -
82 x(te))
92 -Ci,
S11A - A(4)) | B diag of o' (11)) (o((dd((())))
Ci,
where 8 = T/L and we use X(te) = |10.f(te)l & C₁. Additionally, the triangle incquality implies that
||diaglo' (1 [/")) -
(/ X(6)))
<i|diaglo' (h (h²)) ( -xee)l + || (diag]] oo (hh)] -diaglo (beeo((((tt))))>>>>> x(te]) STa||W/11 >,,1 11 + T22lth 14z || A
SC2 (IAA - >1, + 11e - hstel),
where the constant C2 = Cola + Co CoLT, Hence, we have
|X -eee S 11 X(te)l + 3C2 (IIX- -helllell-oo + 11 - heell) + B
Denote E =1AA- dtell + 11h - -heel.. then we have
|X-- - dtell s S (1 BC,)E B1.
By the induction, have
(1 (1+ -1 (1+ BC2) -
E < =1+
+ BC
Since E 0 and B =3=17 -T/L, we have
T C1 e
C2(T Ca(T-tt)-1) -1)
Ex <
L C2
This completes the proof.
Additionally, as we have X₂ = =ffo/Oh, is the solution to the backward ODE, We have
0.e
aft
v [1,2, .. L], -
(48)
Oh.
ahd
where Co > 0 is some constant that is not dependent and T..
26

---

Published as a conference paper at ICLR 2025
E.2 GRADIENT ALIGNMENTS
By using Lemma 8 and TO, we can show the gradients obtained from the optimize-thhen-disrreee and discrete-then-optmmize as the depth + Do. Observe that for any æ. we have
11(hh) - ((h(T))l & CL VASCL
Of
T 1 0s
1w fell
((h,)) dtt
Σ t=1
Vn The ((h
vn Oly,
L
0s (hh)
as-
Th o(h o(h-))) dt
vn
0h,
t=1
4/-1
05
054
05L Oh
||4||| +
111 - 11- -1-11112 1101
Vn
0h,
dh
E=1
C
VnL dt
16-1
L.
<0
Σ
1-2 = =CL
1=1
IVolh Vulell |=|||a| doll S CT
Hence, combining the three results together prove Propostion 2
E.3
NTK FOR -DETT NEURAL NETWORKS
For Neural ODE define equation I
its NTK is given by
Ko x = (VoJo(z), Volo(T))
(49)
As we have shown in Proposition and Ko (z, .) is well defined for every x, ã as observed in Jacot et al. (2018). in the infinitewidthh limit, it converges to an explicit deterministic kernel K called limiting NTK. Hence, we will show that Koo is well-defined and provides its explicit form.
V [e) is well defined for every æ E gd-1 (a.s). Hence, While Ko is random and varies during the training,
3d gd
Recall that we use finite-depth neural network ff defined in equation 10 that approximates Neural ODE fo. As result, we can also approximate the NTK Ko using Ky defined as follows
(x, (Volo (x),, (x), (z))
(50)
We denote KL be the limit of K/- as width + xc. In this subsection, we provide the explicit form for K. For the convergence analysis, we leverage the Master Theorem introduced in (Yang, 2020 Theorem 7.2). This result similar to Theorem but it consider the backward information propagation, and it is reformed as follows.
Theorem 10. (Yang 2020 Theorem 7.2) Consider NETSOR program that has both forward and backward computation for a given finite-depth neural network. Suppose the Gaussian random ini- tialization and controllable activation functions for the given neural network. For any controllable 3) : RRM. M -> IR, as width n - xc, any finite collection of G-vars Ja with size M satisfies
-. M 2-g (0.0M......M11
(51)
where {2²,... (."..... 2 are Gaussian random variables whose mean and covariance are computed by the corresponding NETSOR
27

---

Published as a conference paper at ICLR 2025
Algorithm ResNet of Forward and Backward Computation on Input x
Input: vt//vd:0 G(n)
Input: W n(n,,)
Input: v: G(n)
1: 10 = Un/Vd: G(n)
2: for E f1.2,... 2..2...... D) do
3: x' o(h H(n)
4: g¹ = Wx /
5:
h= h1-1 + G(n)
6: end for
7: xt o(h²) H(n)
8: dx- u G(n)
9: dh dr _ Q d' :(n)
10: for LE E [L, ...
dg dh'" H(n)
12: dx W dg /
/ G(n)
dh dh + dh 1) o dx H(n)
14: end for
Output: ||-| 12/n + Ele (d x /n + xd dhhx Ty /d
As result, this type of Tensor program called NESTOR and includes additional G-vals from the backward information propagation. In our setup, to compute the gradients of fg defined in equation I0, the following new G-vals are introduced
dg b+ diag ø (h"
Vn
dg' diag of (h* W 211,2, L].
vn
and the associated NESTOR T is given in Algorithm 2 In the rest of this subsection, provide rigorous proof to show the convergence K3 KL as stated in Proposition 4
Without loss of generality, we assume on = = 1 and o./Vd = 1. As vec(u, W,U), we have
(z. (V.I (z), (z)) (V Vw (w) (Vuth (z), vuft
Hence, we will show the convergence of each term. To simplify the notation, we abbreviate f f (x) and f= 1z())
CONVERGENCE OF (V.,,V V.J)
By using simple calculus, we have
vof= o(h))(h+\/v
(52)
VnLf @ 0' d'(h')//r.
(53)
By Theorem 10
we have
(V..Vo.)) 1 T o(h2) 35 Bo( CL+1L++((,, x).
where d' 0 + k _
z is a centered Gaussian random variable, of are centered Gaussian random variables defined in Proposition the convergence result follows Theorem 9
28

---

Published as a conference paper at ICLR 2025
CONVERGENCE OF f. Fwf)
To show the convergence, we first rewrite the forward propagation suggested by the Tensor program: for all & E f1,2, ()
g Wx
Vn
he 11" 1--1 + tkg
x² o(h').
By using the chain rule, we obtain
Vwf= =
vn E (v ( (I--1)T
1
(54)
Then, the quantity can be written follows:
1. (Wwf. Σ
"ig',d"))
( /n
f,k=1
where dz is denoted the gradient of f w.r.t. vector z occurred in the forward propagation.
follows from Theorem[9th
1) (x (x 1 (,(ht1)))((h-1)) " (x)
(55)
n
n
Moreover, we have
1 Vi
da
W dg
and
dg' kdh = x(dh 1+1 + dx 2 o'()))) dg 6+1
|dx' o o'h'))
Repeat this recursive relation, and we obtain
dg
dx c s'hh')
(56)
Yang (2220), is equivalent consider the coordinates in dx are asymptoti- cally i.i.d. following some centered Gaussian random variables which satisfies:
By Theorem 10
k
E d f-1 2da
=x
zan' zar d'(")oo(x))
t,k
x
E
Zzz 7dx
E(o'(u')o'''))
Hence, we obtain
2+1,k+1
cC (x, = x2
Σ
Di (1
(57)
iデ=
As result, have
(Wwl o. 2
Σ
DL
(58)
t,k=1
29

---

Published as a conference paper at ICLR 2025
CONVERGENCE OF (Vuf. Vuf)
As Un h = Eド= Uzjzz implies
/ = 6g,
Observe that
Of 0f
(oo..uu) Σ
vi Uij
ohg of
oh of aU ong
oRi 0ND
1/
8
0f of
Oh
Ohg
0hg
ohg
aU 1j aU:
a
0f 0f
Σ a.3
Σ
8a, 1z 10,3110
0he ohg 4
0f 0s
Σ
Sa.3a
0ho 0hg
0
Σ 0f 0f TTa
Dhu O
4 D (r, 3 C x),
where (z.0 B) = =x==xx TT.
Putting everything together yields
(Vol.Vof) = (V... v.J) (Ww/ + (Vel.
0. C (x,
(LM,,,) (x, +
f,k=1
Hence, we obtain Kg (zz 3) converges a.s. to KL (I, 3) defined follows
(x, CLth (c, + Σ CL (z. CO. (2,)))))
E.4 NTK FOR NEURAL ODES
In the previous subsection, we have shown the NTK KJ- converges to deterministic limiting NTK KL as the width n - oc. In this subsection, in the same limit, we will show the NTK Ke of Neural ODE fa defined in equation | converges to the limiting NTK Koc.
Similar to the NNGP kernel Σ*, the NTK Koc can be considered as the limit of a double sequence:
(20 F) lim (VoJo, VoJo) lim lim (VoJe, Volg) L
We have shown lim
K& (x, x) in the previous subsection. Hence, the con- vergence of Koo is equivalent to show the two indices, i.e., depth and width, are interchangeable. Fortunately, if the activation function ø is sufficiently smooth, the two indices are indeed swappable and the NTK K is well defined.
120000
Based on Moore- Moore-Osgood Theorem stated in Theorem a double sequence has well defined iterated limits that are equal to the double limit if the double sequence converges in one index and uniformly in the other. Hence, we will show the NTK K!- as the double sequence converges in depth and uniformly with respect to the width n.
30

---

Published as a conference paper at ICLR 2025
Proof.
Without loss of generality, will assume o, dw 1 and u./Vd= 1. Observe that Ko(o, = (V.Jolæ), V.Jo(T)) + (Vw/o(x) VwJo(F)) + (Volll())VV..(((()
Hence, the rest proof is establish the convergence rate for each term in the summation.
Note that
(V.j(()), v.Jo(E)))
(o(h"(z)) - (o(h(z,T)), o(h(a,
(0(h
(MM (M(xTT)) - (0(1)(()) (),,,, T))) (1* (s)]ll M(, 7) L
n n |1 M1,,, T)1l
L
CV- VnL-
n
=CL-
where we use Lipschitz continuous of é and Lemma
Next, we can first show // and VV Wwff f'|l are upper bounded by some constants as long as T < co. Observe that
Vw =11 {'" Vn
T
<|| Co(T-())
Vn
SCoTe
where we use Lemma 8 and 10
Similarly, we have
L T 0f-
=1
Σ
vn The
f=1
T L 1
0fL
||||
L
Vn
the
1=1
T
1
(1 + (l+oT/L)=- (10t))=-..((-///1111 (1+OT/L)*- Covn
I
Vn
TT
where we have the facts
||~`|| s (1+
(59)
0ʃL the
< (1 + oT/2)4-91a5 /
(60)
for all é (0, 1,..,.1.. I.).
31

---

Published as a conference paper at ICLR 2025
Additionally, it follows from Lemma and
10 that
I W (x) - VwJo(x)
T 05 (bo))dl
- T 1 0sc
Σ T. vn ohe 611111)
ahe
8=1
L
0f
0sL.
(h1)
th 6( (1(-11) dt
Oh,
1
4
05 054
852
Σ
||bll| +
Oh 11 - -1 11at
0h,
Oh
1=1
C L
VnL- de
L
<c EcEne
1, = -
1=1
Hence, we obtain
(w (WW-HH)) (V
(Vw//()) -10 (W
SCL-1
equivalently
(wwr)(o)-wf()))-NNWoo)) Vwl((E) CCL-1
(61)
Next, observe that
Vul - (V
(x) ( ()o())) (3)
(34(1)) ) (ax, 0f((z) 0f((T)
ah(z.0)' ah(i.0)
Then we have
(z) 03)
0f((z) 0f{(3)
)
()) 0l0(3)
OM(x,0) 0) (h(r,
(z)
(f)
(.)
+ (0f'(z)
(.((z) 0f((z)
(x) Oh*(T)
Oh(F,0)
(x)
Jh(r, Oh(0 7))
(z)(0)
|1 (f)
(z)
+|| 0f-(x)
0fe(z)
ff)
(x)
(T) Oh(2.0)
(z)
0h(x, 0)
@h(0,
<CL-1,
where we use the Lipschitz smoothness of ø' and Lemma 10 Therefore, we have
Vo) (Vu Volo(T)) SCL!
Then putting everything together yields
(V)s"(>), vol"(z)) - (Volo()) Volo(T)) S CL.
(62)
Therefore, it converges uniformly in T. and uniformly in n.
32

---

Published as a conference paper at ICLR 2025
Combining Lemma 2 with Proposition 4 and Moore-Oste Osgood Theorem 8 the double sequence (x 2) and obtain the desired result
switch L and n in
Koo (z, σ)
lim Kø(x, Φ)
20000
lim (Vefe, Vefo) n=000
lim lim (Volg: Volo) n-000 1-0
lim lim (Volg, Vole 40000 n
lim K (x) +000
E.5 INTEGRAL FORM OF NNGP AND NTK
In this subsection, we provide the explicit form of the NNGP and NTK of Neural ODEs as the limits of tL and K. It follows from Proposition 5 and Lemma 8 that
<0 60.t d L [0T] "x B) = -o
(63)
Vt, E [0,T],
(64)
where (ur, Ux) are centered Gaussian random variables with covariance
1/")
B(u, ú.) = 0000
(x) a)d') ds'
(65)
Hence, the NNGP kernel of Neural ODE is given by
E'(x, x)) = DT. (v. = = "
(66)
For the NTK of Neural ODEs, we have
(*'(x, S: J. K" (" (a, (),
ds',
(67)
where $2''6' (x. == Ede(ue) )dl (ūw). As a result, the NTK of Neural ODE is given by
(x, >'' (x, (-') 1'" (z) (2, 0 (x, (x.
(68)
F STRICT POSITIVE DEFINITENESS OF NEURAL ODE'S NTK
In this subsection, we will prove the NTK Koo of Neural ODEs are strictly positive definite, We first recall the definition strict positive definite for kernel function.
Definition 2. A kernel function k X X - R is strictly positive definite (SPD) if. for any finite set of distinct points Ti,".. ,NN E X, the symmetric matrix = [k(z,, 2; 25)) is strictly positive definite, i.e., c Ke > for all nonzero vector
Recall that
Ko(2, (VV,,)) V.f(()) + ( Vwle(a)) + (Vofe(a), Vulo(z))
In Theorem we have shown that Ko(x, Ø) + Kx (z, 2 as n + xo, provided @ is sufficient smooth, and
(V.oo()), v.b(()) + Σn (x,
Hence, to show K-. is SPD, it is sufficient to show Σ' is SPD.
Moreover, follows from Theorem that lim Σ (x, σ) Σ' (2, We first show Σ is SPD.
L-+00
33

---

Published as a conference paper at ICLR 2025
F.1 DUAL ACTIVATION AND SPD OF DEPTH NETWORK'S NNGP KERNEL
We first provide the result for finite-depth network Sb defined by 10 where the depth < 00. Proposition 6. Suppose Φ is L, Lippschit continuous. If o is non-polynomial nonlinear, then Σ is SPD on gd-1 for S I. xc.
The proof is based the concept of dual activation and Hermitian expansion. Here brief intro-
duction is provided as follows. For details, we refer readers to Appendices from (Gao et al. 2021
Daniely al. 2016).
Let ~ (N,, 1) and f R + R be a real-valued function. We can define an inner product using expectation:
(f.g)== == (
Thus, we can further define Hilbert space of functions H. that is, f E H. if and only if
|| = (f.f) E= (f)) 00.
Apply Gram-Schmidt process the polynomial functions (1, x², ,..,,,. w.r.t. to the inner product defined before, and we obtain } the (normalized) Hermite polynomials that is orthonor- mal basis to the Hilbert space H:
d
hn(x)) - (-1)
dr
The dual activation :=:-1,1] IR of an activation function & is defined by
(pp) == ((,))),,,
1 p] p 1
where N, is multidimensional Gaussian distribution with mean 0 and covariance matrix
Then the dual kernel K& is defined over the unit sphere Sd-1; for every pair x, F E Sd-1 the dual kernel K5 gd-1 x gd-1 - IR is defined by
Ka(x, (((TT)
If function @ E H. only can obtain an expansion of ø by using the orthonormal basis of Hermitian polynomials but also an expansion the dual activation ø by using the same Hermitian coefficients. As a consequence, the corresponding dual kernel K can be shown to be strict positive definite by using the Hermitian expansion.
Lemma 11. Daniely al. 2016 Lemma 12) If H. then the Hermitian expansion is given by
6(x) = Σ (n))(())
(69)
R=0
00
8(p) Σ amp""
(70)
n=0
where a, : (Ras&) 6) is the Hermite coefficients.
Theorem 11. (Jacot et al. 2018 Theorem G)eiiing
2013
Theorem 1) For function x Sd + R defined by
[-1,]] + R with f(p) =000= 2 bnp", the kernel K₁
S
K,(x,0) f(xTx)
is strictly positive define for any d 2 if and only if the coefficients by for infinitely many even and odd integer n.
Now, with these results, we are ready prove the SPD of L Lemma 12. d is nonlinear and non-polynomial, then E) is SPD.
34

---

Published as a conference paper at ICLR 2025
Proof. We first show E is SPD. As 52 (x. T) o (x, x) and we have
E' (x, = o. (u, [d())((v)]
where
o2 [ G
(x,
(5,,0)
By the notion of dual activation, have
Σ (x o²(
where (xxx) = ((out/va)
Clearly, Lipschitz continuous since Φ is. Then E H. and let the expansion of j in Hermite polynomials {hn J to be given as Σ -==0 On where On (n, hn) are the Hermitian coefficients. Then we can write as (p(p) =n= "ypp" and we have
(x 3) = oni("" = o dx"
n=0
Note that is non-polynomial if and only if & is non-polynomial. As we assume ø is non- polynomial, we have μ non-polynomial, hence there are infinitely many number of nonzero o, in the expansion. That indicates b, := a >> O for infinitely many even and odd numbers. As o > 0, we have 5 strictly positive definite.
Next, we can show if Y is SPD, then y++1 is also SPD for all L 2 1. Lemma 13. Suppose nonlinear non-polynomial 6. Given I < co, then
1. El (x, x² Ei Ch (2, is SPD for all é E (1,2, LL+11 1),
2. E is also SPD.
Proof. As are working with finite-depth network f/ it is fine to assume to simplify the notations. Then Σ and Σ have the recurrent relation, stated in Proposition and so C and C...K, By Theorem 12 we have C is SPD. Additionally, we have
(1 F) = Ed(u (o('') = Es(( (o(")))(()) C*.l,
where we use the fact E[2"=*] 50 (ar Thus, C.. is SPD for all l. Recall that Blu'u* (x, + E Σ Cld (rr, F). Using this relation, we can write
Elu' u' = 00,0(1, .) + (1, 2) + 2 Σ (2, Σ (2,
1.==
As C is SPD for all i, the symmetry of Oh implies E[u is SPD. Now, assume the contrary, i.e., gf+1 cé+1,8+1 is not SPD. Then there exists distinct {or,... NNN and nonzero E RRN such that
2 0 a.ay ara,C 1 (Ii, = Σ ara Blo(()))))) = E Σ acoll() )
We must have E, a.o(u)) 0. As we already show at = (uf.... uxx) R is non-degenerate Gaussian random variables, nonlinearity of implies a = 0, which contradicts a * 0. Hence, E(+1 = 1+1,6+ is SPD.
35

---

Published as a conference paper at ICLR 2025
F.2 STRICT POSITIVE DEFINITENESS OF NEURAL ODE'S 'S NNGP KERNEL
Observe that the previous result uses induction to show _ is SPD. However, the strict positive definiteness of Σ- might diminish L -.. To address this, conduct fine-grained analysis of the properties of _ and demonstrate that these properties persist when L + oo. Consequently, E retains these essential properties. which are crucial for proving that E" is SPD.
Recall from Theorem that
(x,, = Flo("")((T"))
where (u", a²) are centered Gaussian random variables with covariance S (I, defined as the limit of (1, (1) i.e.
(I, = (,, (x +
Σ
(x, as L 1-00 co.
(71)
t,k=1
Based on the proof of Theorem 1 given as follows.
well defined. Some essential properties of s and S are
Lemma 14. Suppose L co. For any I Sd-1, we have
1. 5t (z, x) = (t,
2. 8th (2, 2 5t (x, 2) and the equality holds if and only if x = E
Proof. As L < oc, we can assume k 1 for simplicity. To prove the result, make the inductive hypothesis that Ce. x) = Ce2 (z, 3) for all E, k S L. Then observe that
s (x, ;) = Sz ;)) + Σ
ch.++(xx + ( (x,
Em1
Using the inductive hypothesis, for any E (1,2, + we have C x) Eò(u ( = (o(aa=))((a)) ChLtt
where (ut-1, u are centered Gaussian random variables with covariance
L-1
2-1
= C0.0 (x, Σ
= (x0(,)))
Σ
Eluu Blut-tu4
1==11
i==1
This shows 87 (2,, = -T (0,0) and also (C, C (3, for all t,k < 1
Next, using C (Z,(Ixx)) = C (Z, .). we have
$²(2,x) - -² (z, 2) 2 |I -|-- 2112
1
/ 1o (x) - (
where the function (0(x) == k. Ele olut). This indicates (I, 2 5t 2) and the equality holds if and only if x.
Corollary 2. For any x gd have
1. 0 < (x, St (τ. x) < 20
2. S* (x 2 (x, and the equality holds if and only if
Proof. Observe that
where g(x) lim (x) for -
By Lemma 17 we obtain [g*(2))
L-+00
O (1) uniform in L. Hence, g(x) is well defined and [g(x) O (1). Therefore, we obtain
S S* Φ) = ((1). Additionally. it follows from the relation of S² x) - 5t that
1 S² (x. (z, 2 11 I
lg(z)) g(3))
which allows obtain the second result.
36

---

Published as a conference paper at ICLR 2025
Now, we ready prove the SPD of E" Lemma 15. If & is nonlinear and non-polynomial, then E¹ is SPD.
Proof. As @ is Lipsehitz continuous, we can use Hermitian expansion rewrite E*:
(x = Fra.s))-sraaslouu)
(1, T=0
where we use (u, ū) (x, œ) to denote centered Gaussian random variables with covariance
computed using kernel 5 (z. F) an is the Hermitian coefficients of function vb()) == Φl p(vSjp S/2) with
So := 5* (xx x²), and we also use the facts So = 5* (z, x) S (æ, σ) for all x, 2 and So = ((1) from Corollary
Suppose we are given any finite distinet
{ N
from gd-1 and nonzero e R. Observe that
N
00
N
Σ q.
(x,)"(m,))))=
= Σ
a
Σ
ij=1
N
2
Σ
[ x3 Eg(z.)g(z,))),
a
0
where we use (*x, = xT + (g) By using fundamental properties for positive definite matrices from linear algebra, have
Σ +==1
Cety TTT - =o(r))o()))" + *Bg(X)g(Xi))0n
N
2c (x"(XX]] c Σ
Go xrTT,]";
where @ is Hadamard product. Then we obtain
N
N Σ
Σ
N
n n=U
Gc 2,,"
0
1,1=1
GG
o [o) xz 18)"
i,==
Ge CCFBLuaaaaaa//ssss(00000
Ge Gc.B GC Solo Sol
+F
Σ
Gej B (
GG E (z;
1,1=1
where we use the definitions of Hermitian coefficients an and gb. By Lemma 12 Σ' is SPD and so E* is also SPD.
As corollary result, we have the NTK of Neural ODE is also SPD. Corollary 3. Suppose & and d' are nonlinear Lipschitz continuous. If o is non-polynomial, then the NTK K of Neural ODE is SPD.
37

---

Published as a conference paper at ICLR 2025
Theorem 12. Let { M be training set. Assume
1. 2, E gd- lusl S 1, and z; + 2 for all i i j.
2. the activation & is L1 -iiischhht nonlinear contimuous, but non-polynomial,
3. its derivative &' are Ly-Lipschitz nonlinear continuous,
4. and we choose the learning rate 1) S 1/
For any ō > O, there exists a natural number n& such that for all n. 2 nh the parameter 0 stays in neighborhood of 00, iee.
|00h - -0 S Cl yll/
(72)
and the loss function L(Os) consistently decrease to zero an exponential rate, i.e.
L(O%)= - (1 ndo
&
L(00).
(73)
16
where C' > 0 is some constant only depends L1, L2, Twe Tu and T.
Proof. Given distinct {x. we consider the limiting NTK matrix He E NxN defined as H = (z), As é is non-polynomial, we have Å == Amin {H > O. Let 0g denote the parameters at initialization and (0) E IR NxN be the corresponding NTK computed by Oo at initialization.
Then for any 6e O, there exists a natural number ng such that with probability at least (1-5g) over random initialization Amin (H) 2o/2 for all that with probability at least do). the initial residual || S o.v 2NN log N/ō for all n 2 n Therefore, for any > O, we choose 5n = 6/2, and it follows from Lemma 16 that, with probability least (1 - 6) over random initialization, have
By Theorem we have (0) converges to H the width
no- By Lemma 19 there exists another natural number ny such
||e vl
|| w
1 s C|| wll/>o,
and
|/ wl1 s (1
17/00
1n00 yll.
16
for
2 max fno, n1. No log [N)
G GLOBAL CONVERGENCE OF NEURAL ODES
In this section, we provide the convergence analysis of Neural ODEs defined equation gradient descent.
under
As we use square loss, the loss function given by
(()) (fo( 2.
(74)
By using the vectorization form equation 35 and chain rule, the gradients are given by
aL(0)
N
Te
((r(()))(( - ve).
(75)
Ov
Vn
OL(
N
[' Ow ((((())))) & A((ra)d (So(z:) - y1).
(76)
OW
OL(0)
[ri x (xx)) (fo(x.) -.).
(77)
aU
38

---

Published as a conference paper at ICLR 2025
Consider the gradient descent
0k+1 0h 0l((k)
(78)
00
Assume the inductive hypothesis: For all i S k, there exist some constants Quy Awy On. that
such
1. |0 ll [| || S CVT.
2. Iu vll $ 110))) 5/1,
where C > 0 is a constant and on min vn $0)
Without loss generality, we assume o, = 1, Fw = o, u./
and L1 = L2 1.
Observe that
050
||
-1 [(rr)l -
1
|
|
Ov
vn
Note that
05a
5 vn /-
0
OW
s vn {- T
a
|11
|||||| (e) dt
|10|1 ||l||
(oT)
| w
/n vn
Observe that
f au
||a|| vn
< |oll & |||||
exp (OT][WII/VW)
By using the inductive hypothesis, obtain
< Ce CoT
|||||
(79)
0v
Ofe
S (oT)Ce
(80)
Ofe
aU || < Ce
Coo
(81)
Then we obtain
OL(0')
l1okt1
11 Sn Σ1
Dr
i=0
S1
Ce |X|l| yl1
SnCc |X|
(1 naa))) y/1
CoT vll/ag
Note that the RHS is an constant after initialization. If we assume ||ж||| = 1 and |g/ = 1, then we need ensure
Ce COT
X X1ll vll/ad s CVn.
(82)
39

---

Published as a conference paper at ICLR 2025
And as a result, we have Similarly, we have
||
s ||| - -011 + ||00|| S CVT.
aL(O')
| - 1 S1
OW
1=0
(TT) CaT
X|llla vll
k Σ
Sn(oT)Co COT |X|
(1 naa)ll vl
1=0
S(OT)CE
|X|l| - yll//0
Then need to ensure
(oT) | [x]])
- yl1/a - CVT.
(83)
Then we obtain
| S || - | + | S CVn.
Observe that
OL(09
- 1
2=0
Ce
| - yll
k CoT [ Σ (1 - m yl
1=
SCeCoT|x
-yll/a
Hence, we obtain
|| S || - - + |1 s Cvn
Next, observe that
k-1
k+1
1 + -u² - y)
=
Ou
(ak+1 -Oh) + (ut - -2)
00
Out
(ut y) + (u* v)
00
1| 00
I 1
(u+ 1 v)
30
00
Out
]
( Ou
Pū
I 1
(u* - y) +
(ut
y)
0e
00
Oe
00
where ă = v0) and 0 is an interpolation in between @ and o++1
Note that
0f
of
1
1
1
=1 ((hT)) ((rr)!
dv
Du
Vn
vn
1
<
nn ||/r- - hr
C
Vn |10 @lle
CCT
40

---

Published as a conference paper at ICLR 2025
where we use the Lemma and the inductive hypotheses. Similarly, note that
0f
0f
0
((ha) 20 Xx - (((h)) s0 Aedell
aW
aW
Vn
!'
T
(11b: helll/ell |i 1112 Jell)
vn
1"
110 0||0000||>|1
co((T
Vn 0
C
((oT)
110 oll Co ||r|
vn
and
0f
of
C -
$
|
Xoll
10
0|e
0
c
CoT
OU
OU
vn
Hence, we have
0f
00
of
00
=
Of
dv
0f
du
OW
0f
OW
0f
+
dU
Of
OU
0f
$ (oT)
10 -
vn
C
B|ee |||||
Then
Ouk I1 00
0ū
0g || s (oT)
vn 1111s 110 01|00oT|110 S (oT) C
C
-110 gk+z 1et21|e |X|
where we use the fact 0 = atk
+ (1 - abgtt for some a E [0.1].
Observe that
OL(Ok)
Out
||@k+1 8411 7/l
7/ll
(u* -ull
w S ((dT)cccoo yll.
00
Hence, we obtain
0ū
C
S nooor)
_e CoT
| 1
vll.
00
Vn
and
Ou k
Ou
S(oT) Ya
- -
CoT |X|
00
00
k-1
<(oT)
Σ
||0tt1
8'11
1
C
SnooT) Vn
CoT
1X
11 -v--vl vll
Sn(oTT/ Vn
CoT 1
Σ (1 - "1"" 1/l
C
S(oT)2
COT | vll/a8
Vn
500/2,
where we use the assumption
Vn >
CoT yll/ag-
(84)
41

---

Published as a conference paper at ICLR 2025
It follows from the Weyl's inequality that
Ouu
Ouk
Ou
d min
min
2 an/2.
00
08
09
0g
and
Out
Ou
Amin
2 o
00
00
Therefore, we obtain
yll + ( C
|| vl1 - [1 n///44 111t
||| 1
Vn v
vll
1 - -0 17 (OT,3
C
COT
|
w1l|
11 - 2ll
Vn
1 a - n(TT)s
CoT co
1X 2/1
11 - vl
vn
nog/8 /8 1| vll
where we assume
Vn 2 8Cn((T)) CoT 1x vll/o0-
This finishes proving Lemma 16
Lemma 16. Assume o and o' are L1- and L Lipschitt continuous and Nu == Amin (Ko.) > 0. Suppose choose the width || yll2/>8) and the learning rate n 1 Then the parameters H stays in the neighborhood of 60, i.e.,
- 11 || w |IU U S C||x yll/Aa.
(85)
and the residual |lu yll consistently decreases, i.e.,
l - -111 S (1 (1- 11A0
1 ull.
(86)
where C > 0 is some constant only depends on Li, Lzz au ou aw and T.
Lemma 17. Given 0 we have
{ at vn (T
|
|hell s ||
exp
(87)
|||||
|w]l
|ell <
exp
(88)
Vn
Vn
for all L [O,
Proof. Observe that
h4 = ho +
w
(hh))ds
vr
and so
vn ||W1] //'' |/hallds
hell < |
Then it follows from the Gronwall's inequality that
al /
[w]) Wi [0,T].
|
<
||||||
exp
(89)
42

---

Published as a conference paper at ICLR 2025
Similarly, we have
de = xr 1-' a diag e (he)) Asds
implies
A S |rll + vn Lill S. I|Alllds.
By the Gronwall's inequality, we obtain
/ a/IWIl/ Vids
| < ||| exp
>| exp Will/// f)}
By \r diaglo' (hr)]e, obtain the final result.
Lemma 18. Given 0, 8, we have
|
|1, hell
5
01
0 WW+WWW//W ||
(90)
[
|||||
|1 Xell
<|10
81|
c (OTT-0([]1) (/V /vn
(91)
for all t e [0,
Proof. Observer that
he he-h hs = (ho - ho) + - ( [Wo(ha) W [hhh)] ds
Then we have
l|ke hell 511h0 - holl + Vn {" [II1 - Will + |||l| hall] ds
</1 holl Vn |W W|ll / |Ua|| exp aas /I/////// ds + In |W|| /
1/,
hes lds
Using the bound of ||hg]l. we have
|U>||||- Wil (' exp (asl ds
Vn |U>||||W -WIl (% G
|U]
|| - W| ((-tIW//M ||
|
Then by Grownwall's inequality, obtain
|1||
|| h.lll s [1ho holl + I|W WI
- - ||||| -1ll111/
[W]
||U||
S (I + I WIl)
||]
Then we obtain the result.
Observe that
A 7, (Dr 2x)) + In 1'
X, 3.-dian --diigl diaglo'(h,.)
X,
43

---

Published as a conference paper at ICLR 2025
Then obtain
S. (IWV -- + || III 2xl1) ds
|| Xell S Ar
Arll +
By using the bound of
| AzlI. we obtain
|ell
o -
|M Jas
σ
W
WI|
exp
vn
vn
/n
|||||
|ill
( (1)
IW-
W||
Vn
vn
vn
|||||
W WIl
[
Then by Grownwall's inequality, have
|1a|| ( (elt--o)wwlVv -1) cOCT-OIIIIV
D - xe1l 1 | - - +
W
WIl
S 1
|oll
(110 - el| + I WIl)
(T (T- +IWII//
[|
Lemma 19. Given ó > 0, there exists a natural number ns such that for all n 2 ni, with probability at least 5 over random initialization, we have
||l|| S / 2N log(N/o),
(92)
where 0² Ex, St-1
Proof. Fix denote u == fo(x) T. ((t((z))// By Theorem distribution to centered Gaussian random variable with variance o² 8 > O, we have there exists ns such that 2 nı implies
we have converges in
E (), Hence, given
P(u 8) P(= 2 e) < 5/2.
where = ~ o2). Then we have
P(u 2 < 8/2 + P(= 2 8) < 6/2 + p2 s
where the last inequality is due to o 2 log(2/0). Similarly, we obtain two two-tailed bound, ie.,
P(|u| 2 2) 6.
Now, denote u = fx(X) E R. as vector. We have
N
P([|u|| 2 cu) =(||u||| + P(
|u 2 2 =
N P([u+| 2 2 cg/N)
N
= Σ
P(|u,l 2 co/VN)
1=1
-=1
<,,
where we use the fact PIEN 2 8)
< MN P(æ, 2 e/N) and eg == 2N log(N/6).
H ADDITIONAL EXPERIMENTS
In this appendix, we provide supplementary experiments that complement the results in the main paper. These experiments explore the impact of different activation functions, scaling for long time horizons, and the behavior of Neural ODEs when approximated by Gaussian processes. Addition- ally, examine the behavior of the NTK when using polynomial activations.
44

---

Published as a conference paper at ICLR 2025
H.1 SCALING FOR LONG-TIME HORIZONS
As discussed in Proposition I and Proposition 2 smooth activations ensure that the forward and backward dynamics of Neural ODEs have globally unique solutions. However, extending the time range or working with long-time horizons in the dynamics can introduce difficulties for numerical solvers, leading to higher numerical errors. To understand how Neural ODEs behave over extended time horizons, investigated their behavior at initialization as the time horizon increases, focusing on how output magnitudes and variance are affected. The objective was to understand how extending the time horizon impacts the model's outputs and the subsequent training process.
At initialization, as the time horizon T increases, the output magnitudes grow larger, resulting in increased variance, as shown in Figure a). This increased variance negatively impacts the training of Neural ODEs with gradient descent, leading to damping in the early stages of training as illus- trated in Figure ()). We observed that increasing the width of neural networks reduced the output variance, as shown in Figure Additionally, for long-time horizons, where T is large, we suggest scaling the dynamies by setting the weight variance Tw 1/'T'. This approach effectively mitigates the growth in output magnitudes, as shown in Figure 3(c) c), and reduces early-stage damping during training. as illustrated in Figure dd).
Taltand Bez
and Lestee Wirth (Soaled)
with
Leki I 201
wet With 1906
(a)
(b)
(d)
Figure 3: Effects of increasing time horizons on Neural ODE outputs and training. (a) Output variance increases as the time horizon T' becomes large at initialization.
(b) This leads to damping
during the early stages of training with gradient descent. (c) Scaling the dynamics by setting the
weight variance ow ~ 1/T reduces the output variance. (d) This scaling also mitigates early-stage damping, improving training stability.
H.2 GAUSSIAN PROCESS APPROXIMATION
In Section 4, we established that Neural ODEs tend toward a Gaussian Process (GP) as their width increases, demonstrated in Theorem 1 The associated NNGP kernel of this Gaussian process is non-degenerate, as stated in Lemma To empiricaly verify these theoretical findings, we conducted series of experiments.
First, we fixed an input æ and initialized 10,000 random Neural ODEs. We then plotted the output histograms for various network widths and fitted the distributions with Gaussian model, as shown in Figure Additionally, we ran statistical tests to confirm whether the output distributions followed a Gaussian distribution. The Kolmogorov-Sminng Sminnov (KS) test statistics and p-values indicated that as long as the width exceeds 100, the outputs closely follow Gaussian distribution.
Next, we analyzed the independence of the output neurons by plotting pairwise outputs across two coordinates. According to Theorem [ the output neurons should become independent as the width increases. Figure [ confirms this: while the diagonal plots show Gaussian bell shapes, the off- diagonal plots resemble random ball shapes, indicating that the neurons are uncorrelated and, there- fore, independent as the width increases.
Finally, we investigated whether Neural ODEs preserve the structure of input data at the output. We constructed a matrix X of 10 samples and calculated the input covariance matrix X Then, we initialized 10,000 Neural ODEs with random weights and evaluated them on the input X, computing the output covariance matrix. As shown in Figure 6 the output covariance matrix retained the cor- relation patterns of the input matrix but with reduced magnitudes, indicating that Neural ODEs act as structure-preatt preservvng smoothers, reducing the spread of the data while maintaining its underlying relationships.
45

---

Published as a conference paper at ICLR 2025
Histogram of Outputs (Width=10)
Histegram of Outputs (Width:50)
Histogram Outputs (Width=100)
Goussin ftt Hategram
Gaussion Tt Vistogram
Gausdann ft: etttogmm
De-stt
Dersaty
01
Durpot Value
Outpor Value
Output Value
(a) KSS Statistic: 0.0730,
(b) KS Statistic: 0.0200, P-value: 0.0366.
(c) KS Statistic: 0.0130,
P-value 0.0000.
P-value: 0.3609.
Histogram of Outputs (W0thh2000
Histogram of Cutputs (Width=500)
Histogram of Ourputs (W0dth-1000)
Coussian Hitoggm
Halogram
Denafty
21
Daippl Vallo
Cutpuz Value
Clutpul Matuee
(d) KS Statistic: 0,0085,
(e) KS Statistic: 0.0079,
(a) KS Statistic: 0.0069, P-value: 0.9688.
P-value: 0.8600.
P-value: 0.9084.
Figure 4: Gaussian fit of the sample distribution from 10,000 randomly initialized Neural ODEs across widths 10, 50, 100, 200, 500, and 1000. The corresponding KS statistics and p-values are displayed, showing improved Gaussian fit width increases.
hfrplr opee
Oupprs
Pepple OCC
Figure 5: Pairplots of output neurons given the same input data, showing that output neurons become independent as network width increases.
Additionally, we computed the smallest eigenvalues of the output covariance matrix. By sampling 500 examples from the MNIST training set and computing the covariance matrix for 10,000 inde- pendently initialized Neural ODEs, we found that the smallest eigenvalues became strictly positive when the width exceeded the number of samples. This result indicates that the NNGP and NTK kernels are strictly positive definite, aligning with our theoretical findings.
H.3 SMOOTH vs. NON-SMOOTH ACTIVATIONS
To compare the performance of smooth and non-smooth activation functions, we evaluated both Softplus (smooth) and ReLU (non-smooth) across different widths, measuring their training and test losses, parameter distances, and NTK least least eigenvalues. The results highlight several key differences between these activations.
46

---

Published as a conference paper at ICLR 2025
Covariance Matrix
Outout Ceveriance Mabris Witth 200
Smallest Eigewvalue vs. Wisth Neural CDE
000000
800011
Bjetnoace
520050
000
(c) Least Eigenvalues
(a) Input Covariance
Output Covariance
Comparison of input and output covariance matrices. (a) Input covariance matrix. (b) Output covariance matrix from Neural ODEs, showing similar structure but reduced variance. (c)
Figure 6:
Least cigenvalues of the covariance matrices, confirming positive definiteness width increases.
Training and Test Loss Behavior: As illustrated in Figure (a)- the log-scale plots show that Softplus consistently converges faster than ReLU. ReLU experiences slower convergence, particu- larly during the early stages of training, while Softplus benefits from smoother and faster optimiza- tion, especially at larger widths.
We also measured the distance between the model parameters and their initial values throughout training. As shown in Figure (i), parameters remained relatively closer to their initialization for Softplus, while the parameter distance for ReLU was significantly higher. This suggests that Sofiplus's ' smoother nature results in more stable parameter updates during training, contributing to its faster convergence.
Parameter Distance:
NTK Least Eigenvalues: Regarding the NTK least eigenvalues, the results shown in Figure 7 dd) indicates both activations exhibited strictly positive eigenvalues, with ReLU's slightly larger than Softplus's. s. However, despite this, Softplus converged more rapidly based on both training and test loss results. We hypothesize that Softplus's smoothness allows for more accurate gradient computa- tion by the numerical solver, leading more efficient optimization compared ReLU, which may suffer from less precise gradient calculations due to its non-smooth nature.
In summary, Softplus not only converges faster in terms of loss, but also leads more stable parameter behavior during training, despite the slightly smaller NTK least eigenvalues compared to ReLU. These findings suggest that the smoother nature of Softplus provides significant advantages for Neural ODE training.
Toir Lep Compertice Setalu et ESU L9g stele
Texl Loen Comparber: Solte SE'S 20g Soilt
Schylus
Tan Salana, A0a-
Golake Viee
Tum Nolin 1905
lan Solaks Vien
From Aelute
Sclaie I
Paten (esaacc 2611 Vicn 3003 Penan
Diuia furuine Vie
(a)
(c)
(d)
Figure 7: Comparison of Smooth (Softplus) and Non-Smooth (ReLU) Activation Functions across Neural ODE Training. (a) Log-scale plot of training loss shows faster convergence for Softplus compared to ReLU, especially at larger widths. (b) Log-scale plot of test loss similarly demon- strates faster convergence for Softplus, with ReLU exhibiting slower progress in the early stages. (c) Parameter distance from initialization, showing that Softplus keeps parameters closer to their initial values, suggesting more stable updates. (d) NTK least eigenvalues for both activations, where ReLU's eigenvalues are slightly larger, though Softplus achieves faster overall convergence.
47

---

Published as a conference paper at ICLR 2025
Trass Low Dorparscy SeN VE
Thit Loss Compation Sctprst Quad
Puram: Dar Compestiers Setplas Dasd
MIX bug Conpurser: Soltphe Quod
Bal Bulghio
Buluko A 1330
Bain Out Waan 100
Hat cocd vee 2009
Uik Que Die
Dusd vict
tade Veo
Mil juft 2018
Poas cadaa ca autinia
Wit Duai Wudiv1300
Perse (leaart Daed was 1030
(a)
(c)
(d)
Figure 8: Comparison of Neural ODEs with Softplus (non-polynomial) and Quadratic (polynomial) activation functions. (a) Training loss for Softplus and Quadratic activations, showing slower con- vergence for the Quadratic case. (b) Test loss comparison, further illustrating the slower convergence for Quadratic activations. (c) Parameter differences from initialization, showing that Quadratic ac- tivations lead to slightly larger parameter drift compared to Softplus.
(d) NTK least eigenvalues, where both activations show strictly positive eigenvalues, with Quadratic's s being slightly larger than Softplus's.
H.4 POLYNOMIAL ACTIVATIONS FOR NTK AND GLOBAL CONVERGENCE
In this experiment, we tested quadratic activation functions to assess their impact on NTK behavior and convergence. While previous results suggested that nonlinearity but non-polynomiality is sufficient condition for the strict positive definiteness (SPD) of the NTK, our experiments reveal that it is not necessary condition.
We observed that the NTK of Neural ODEs using quadratic activation is also strictly positive definite, with the smallest eigenvalue slightly higher than that of Softplus, as shown in Figure d().. Despite this, the quadratic Neural ODE converged much more slowly than Softplus, as illustrated in the training and test losses (Figure [b))
In terms of parameter behavior (Figure (())), the parameter differences for the quadratic activation were slightly larger than those for Softplus, meaning the parameters drifted further from their initial values. However, these differences remained within the same order of magnitude, indicating that the model still satisfies the conditions for global convergence, even though it does meet the sufficient condition of being non-polynomial.
In summary, while quadratic activation functions result in strictly positive definite NTKs similar to non-polynomial activations, they lead to slower convergence and slightly less stable parameter be- havior compared smoother activations like Softplus. This suggests that while pon-polynomiality is not strictly necessary for SPD and convergence, smoother activations may offer practical benefits for faster and more stable training.
H.5 CONVERGENCE ANALYSIS ON DIVERSE DATASETS
In the main paper, we focused on the convergence properties of Neural ODEs using different activa- tion functions on the MNIST dataset. To ensure that these findings generalize across different types of data and tasks, we extended our experiments to three additional datasets: CIFAR-10 (image clas- sification), AG News (text classification), and Daily Climate (time series forecasting). This section details the performance of three key activation functions-Soffi datasets, highlighting their effects on convergence speed, stability, and generalization.
ReLU,, and GELU on these
For each dataset, we trained Neural ODE models with different widths (i.e., 500, 1000, 2000, 3000) using Softplus, ReLU, and GELU activations. We monitored the training loss and test loss, compar- ing how different activations influence convergence behavior across datasets. The optimizer used was gradient descent with a learning rate of 0.1, and models were trained for 100 epochs.
For CIFAR- 10, the results showed minimal differences between the activation functions.
Softplus, ReLU, and GELU all exhibited similar convergence patterns, with larger widths leading faster convergence across the board.
48

---

Published as a conference paper at ICLR 2025
Larger widths consistently resulted in lower training and test losses, but the specific choice of activation did not have a significant impact on the overall performance or convergence speed.
These results suggest that for CIFAR-10, the activation function choice is less critical, particularly when the network is sufficiently wide, Le., see Figure
mash 2000
(a) Softplus
(b) ReLU
(c) GELU
Figure 9: Training and loss behavior for CIFAR-10 across different activations: (a) Softplus, (b) ReLU, and (c) GELU. All activations show similar convergence patterns, with larger widths leading to faster convergence.
For AG News, we observed distinct convergence patterns across the activation functions:
Softplus converged the fastest, followed by RcLU, with GELU converging the slowest. Despite GELU being smooth activation, its derivative differs significantly compared the other activa- tions, which may explain the slower convergence rate.
All three activations shared the same trend: larger widths led faster convergence and lower test losses. However, the differences between activation functions were more pronounced at smaller widths, where GELU lagged behind (Figure I0).
This suggests that while GELU's smoothness offers theoretical benefits, in practice its derivative may cause slower optimization dynamics, particularly for text-based tasks like AG News.
(a) Softplus
(b) ReLU
(c) GELU
Figure 10: Training and test loss behavior for AG News across different activations: (a) Softplus, (b) ReLU, and (c) GELU. Softplus converges fastest, while GELU lags due to its derivative behavior.
tsaliaatio
Epacss
(a) Softplus
(b) ReLU
GELU
Figure 1I: Training and test loss behavior for Daily Climate time series forecasting: (a) Softplus, (b) ReLU, and (c) GELU, All activations show similar convergence, with larger widths leading to faster results.
49

---

Published as a conference paper at ICLR 2025
H.6 NON-SMOOTH ACTIVATION FUNCTIONS: COMPARING "DISCRETIZE-THER THEN-OHTINNEEEE OPTIMIZE" AND "OPTIMIZE- THEN- DISCRETIZE"
(a)
(c)
(d)
Figure 12: Comparison of Neural ODEs and ResNets with Softplus (smooth) and ReLU (non- smooth) activations. (a. b) Softplus: Output and gradient differences in log-log scale, both showing L/L convergence. (c, d) ReLU: Output difference shows 1/L convergence, while gradient differ- ence remains constant as depth L increases.
(
(c)
(d)
Figure 13: Training comparison of Neural ODEs and ResNets with ReLU (non-smooth) activation. (a) Training loss decreases consistently and aligns closely for both Neural ODEs and ResNets. (b) Test loss shows similar trend, consistently decreasing for both models. (c) Output difference remains consistently small throughout training. (d) Gradient difference oscillates during training.
In this experiment, we investigate the impact non-smooth activation functions, specifically ReLU, on the performance of Neural ODEs and their ResNet approximations under the DDiscretize-Then frameworks. While the output differences between the
Optimize" " and "Optimize- two frameworks decrease as the depth L increases, our results reveal that the backward gradients fail to converge due to the non-smooth nature of ReLU's derivative.
Smooth Activation Functions (Softplus). For smooth activation functions like Softplus, both the output difference and gradient difference between the two frameworks decrease at rate of 1/L as the depth L increases. This behavior aligns with Proposition and is illustrated in Figure 12 ab))).
Non-Smooth Activation Functions (ReLU). In contrast, for ReLU, the output difference still decreases at a rate of 1/L, as shown in Figure (2.c). However, the gradient difference fails to converge, as illustrated in Figure d2. Initially, the gradient difference reduces as depth increases, but eventually stagnates at a fixed error level. Increasing the network width does not resolve this issue. Notably, the largest gradient difference is observed at width 500, whereas smaller errors achieved for both smaller and larger widths, such as width 200 and 1000. These results confirm that the lack of continuous derivative in ReLU introduces inconsistencies in gradient computations between the two frameworks.
Training Dynamics. Despite this mismatch in gradient computation, we did not observe signif- icant differences in the training dynamics between Neural ODEs and ResNets. We trained Neural ODEs and their finite-depth ResNet approximations (fixed depth 200, as further depth increases did not reduce errors, as shown in Figure d)) on a subset of MNIST. As illustrated in Figure [3]] both models exhibit similar training and test losses. While the output differences remain consis- tently small during training, the gradient differences oscillate, as shown in Figure [3] ResNets, as finite-depth networks, are known to exhibit global convergence guarantees under gradient descent in overparameterized regimes (Du et al. 2019a), so their convergence is unsurprising. What is unex- pected, however, is the near-identical training dynamics between Neural ODEs and ResNets despite the gradient mismatch caused by ReLU's moo-smoothness. Our hypothesis is that while gradient differences oscillate during training, they remain within small magnitudes because MNIST is a sim- ple dataset and ReL.U's derivative almost continuous everywhere except the origin. This partial
50

---

Published as a conference paper at ICLR 2025
smoothness may mitigate the adverse effects of the gradient mismatch. However, anticipate that in realistic applications involving more complex datasets, these differences could lead to divergent training trajectories and dynamics for Neural ODEs and ResNets using non-smooth activations.
H.7 SENSITIVITY OF "OPTIMIZE- THENT-HHENN---E-EEE TO ODE SOLVERS
(a)
(b)
Figure 14: Sensitivity of "Optimize- peptimize-then-eee dieenddiceetize" to ODE solvers. (a) Training and test losses decrease consistently for all three solvers at width 500. (b) Training and test losses decrease consis- tently for all three solvers at width 2000. (c) Time taken by each ODE solver across different widths, highlighting the scalability advantage of fixed-step solvers.
In this subsection, we investigate the impact of different numerical ODE solvers on the accuracy of gradient computation and overall training dynamics the "Optimize-then-dieeeriie framework. The solvers considered in experiments are Euler, rk4, and dopri5.
As illustrated Figure [4(a) (b), the choice of ODE solver does not significantly affect the accuracy of gradient computation or the overall training dynamics in our specific setting. This is consistent with the theoretical guarantees established in Proposition [ where we demonstrated that the ODE dynamics in Equations Eq. (2) and Eq. 4 possess globally unique solutions under the smoothness conditions on activation functions. Given the relatively simple nature the system studied, the numerical errors introduced by the solvers appear be negligible in this context. However, this observation may not generalize to more complex systems or practical applications where numerical errors can be influenced by other factors such as stiffness or stability in the dynamics, which are beyond the scope of this paper.
An interesting observation from our experiments is the computational efficiency of the solvers. While adaptive solvers like dopriS provide high accuracy. they require significantly more computa- tion time the neural network width increases. In contrast, fixed-step methods such Euler and rk4 scale more efficiently with width, making them preferable in scenarios where computational cost a concern. This is illustrated Figure [4(c), where we compare the time taken by the solvers across different widths.
I DISCUSSION ON GENERAL DYNAMIC FORM NEURAL ODES
In this section, we discuss extending results from the specific form equation and equation 2 to a more general dynamic formulation. Specifically, we first consider generalized nonlinear trans- formation:
Tw w f(h,, Vn
E [0,7],
where the original nonlinear activation function @ equation replaced by a general nonlinear mapping f I8" x IR + RR, defined as f : (h, + f(h t). This generalization introduces explicit time dependence, transforming the system from an autonomous to a non-autonomous system. Non- autonomous systems are prevalent in applications such as diffusion models (Song et al., 2020) and physics-informe informed neural networks (PINNs) (Sholokhov et al. 2023). The function f can represent another shallow neural network al. 1998). gating mechanisms (Hochreiter 1997). attention mechanisms (Vaswani 2017), or batch normalization (Toffe] 2015).
complex operations, such as convolution layers (LeCun
51

---

Published as a conference paper at ICLR 2025
For this generalized form, the backward dynamics take the form:
X Ow J (hu, t) W X
Vn
where = 0f/Oh R"** is the Jacobian matrix of f with respect to h. By Theorem 5 the forward ODE has unique global solutions f is continuous in and Lipschitz continuous in h, with Lipschitz constant independent of I. This generalizes the continuity requirement for the activation function @ to f. Additionally, if the Jacobian matrix J globally bounded, the backward ODE also admits a unique global solution. Since f is Lipschitz continuous in h, the boundedness of J is naturally satisfied (on a compact set). Therefore, appropriate smoothness conditions on f ensure well well-posed forward and backward dynamics with unique solutions.
Using Euler's method, we discretize the forward and backward dynamics as follows:
h+ h² +
Tw Vn
W f(h bmy)s
x x vn
J(h bmy W X
where x L. Ensuring convergence of (h' X²) to (he Ae) is critical to aligning the gradients ceemmiie-then-dieeeee
methods. As discussed in Proposition 2 additional smoothness of the backward ODE is required for gradient equivalence. By Theorem 7 the mapping i - J ((,,, t)) must be continuous in ł, which implies that J is Lipschitz continuous in h with Lipschitz constant independent of The smoothness J with respect to f. Specifically, bounding
obtained from the ftiscretize-then then-optimize" and
h can be guaranteed by imposing second-order regularity conditions the Jacobian tensor /ah under suitable norms, such as the operator norm or Frobenius norm, ensures the required regularity. Although OJ /J represents higher-order tensor, these regularity conditions allow the gradient consistency results from Proposition to extend seamlessly to this generalized formulation.
Theorem 7 provides not only convergence guarantees but also a uniform convergence rate un- der globally uniform smoothness conditions. Consequently, by Theorem 8 the iterated limits in Lemma D and Lemma 2 converge to the same double limit. As result, the NNGP and NTK of the generalized Neural ODE remain well-defined. If the limiting NNGP or NTK strictly positive definite (SPD). global convergence under gradient descent can also be established.
Finally, we discuss extending the dynamics to a post-activation formulation:
h f a
Wh Whe,
where the linear transformation h + Cx Wh is applied before the nonlinear mapping f. The anal- ysis remains analogous because the linear transformation is globally 1-Lipschitz continuous under Theorem 4 However, we focus primarily on the pre-activation form, as it consistently achieves superior empirical performance compared the post-activation formulation (He et al. 2016b).
.nn
J RELATED WORKS
Neural Ordinary Differential Equations (Neural ODEs) (Chen al., 2018b) introduced a continuous- depth framework for modeling dynamics by replacing discrete-layer transformations with parame- terized differential equations. This innovative framework has since inspired extensive research, leading to both theoretical advancements and practical applications.
Neural ODEs are distinguished by their continuous-time representation and memory efficiency through parameter sharing, setting them apart from traditional architectures like ResNet (He et al.) 2016a). Building on this foundation, several extensions have been proposed to address more complex systems. Notable examples include Neural Stochastic Differential Equations (SDEs) for stochastic dynamics (Li et al., 2020), Neural Partial Differential Equations (PDEs) for spatiotem- poral systems (Sirignano & Spiliopoulos 2018 Raissi et al.) 2019). Neural Controlled Differential Equations (CDEs) for irregular time-series data (Kidger et al. (020)). and Neural Variational and Hamiltonian Systems for capturing conserved quantities in physical dynamics (Greydanus et al.
52

---

Published as a conference paper at ICLR 2025
2019). These advanced formulations have broadened the applicability of Neural ODEs to diverse domains, such as time-series modeling (Rubanova et al. 2019), computer vision (Chen et al. 2018b]
Park et al.] 202I). and reinforcement learning (Du et al.] 2020). In generative modeling, Neu- ral SDEs underpin approaches like FFJORD (Grthhwhll et al.] 2018), score-based methods (Song et al. 2020), and diffusion models (Ho et al.) 2020). Similarly, in physics-informed machine learn- ing, Neural PDEs and thysics-Informed Neural Networks (PINNs) have proven critical for solving physical systems while incorporating domain-specife knowledge (Sholokhov et al. 2023 Karni- addaii et al., 2021 Raissi et al., 2019). However, while these features offer flexibility and efficiency, they also introduce significant challenges during training.
A key challenge in training Neural ODEs lies in gradient computation. The original adjoint method introduced by Chen et al. (20T8b) computes gradients with minimal memory overhead. However, this approach can suffer from numerical instabilities, as observed in Gholamineeja et al. (2019). To address these issues, advanced methods have been developed. For instance, [Zhuang et al. (2020a) integrates adjoint techniques with checkpointing to balance memory usage and computational cost, while Matsubara et al. (202T) employs symplectic integrators to preserve ODE structure, ensuring stability in long- long-time horizons and oscillatory systems. Finlay et al. (2020) regularizes the Jaco- bian norm of the dynamics to improve stability and generalization. Ko et al. (2023) introduces a homotopy-based approach, starting with simplified dynamics and gradually transitioning to target approach, where (aug-
dynamics. These methods generally follow an septimize-then-dieeee mented) backward ODEs are solved numerically to compute gradients. Conversely, the "discretize- then-optimize" approach, which discretizes the forward ODE into finite-depth network for gradient computation via backpropagation, has been explored by Massaroli et al. (2020). However, as noted in Zhuang et al. (2020a b). this method often results in deeper computational graphs, raising con- cerns about gradient accuracy.
To address the challenge in gradient computation, several theoretical studies have been conducted, focusing on well-posedness and stability. For instance, [Gholaminejad et al. (2019) highlighted sig- nificant numerical instabilities when using ReL.U activations in Neural ODEs. Meanwhile, [Ro-] driguez et al. (022)) investigated the stability of Neural ODEs through a Lyapunov framework de- rived from control theory. Despite these advancements, none of these works address when and how and "optimize- wewpniiiitthenniicretttee ents. Moreover, the question of whether simple first-order optimization methods, such stochastic gradient descent, can reliably train Neural ODEs to convergence remains unexplored.
the
methods can yield equivalent gradi-
Another essential challenge lies in analyzing the training dynamics of Neural ODEs due to the in- herent nonconvexity of neural network optimization. A significant breakthrough in this area came from the Neural Tangent Kernel (NTK) framework introduced by [aacot et al. (2018), which demon- strated that the NTK governs the training dynamics of feedforward networks (FFNs) under gradient descent and converges tates global convergence guarantees for gradient- based optimization in overparameterized regimes, provided the NTK remains strictly positive definite (SPD) (Du et al.) 2019a] Allen-Zhu et al. 2019] Nguyen| 2021). The strict positive definiteness of the NTK has been extensively studied, beginning with dual activation analysis for two-layer networks (Daniely et al., 2016) and later extended to finite-depth FFNs (Jacot et al.) 2018 Du et al. 2019a). Recent work has further applied NTK the- ory to diverse architectures, including convolutional neural networks (CNNs) (Arora et al. 2019). recurrent neural networks (RNNs) (Yang) 2020). transformers (Hron et al. 2020). physics-in formedd neural networks (PINNs) (Wang et al. 2022). and graph neural networks (GNNs) (Du et al. 2019b). NTK analysis has also been explored for various optimization methods, such as stochastic gradient descent (SGD) (Zou et al.) (020) and adaptive gradient algorithms (Chen et al. 2018a). A few recent works start studying large-depth neural networks Gao & Gao (2022b a); Gao 2022).. However, ap- plying NTK theory to continuous-deye depth models like Neural ODEs and determining whether similar SPD and convergence properties hold remains an open and active area of research.
deterministic limit as network width increases. This convergence facili-
Despite significant advancements, challenges persist in understanding the training dynamics of Neu- ral ODEs and ensuring gradient consistency between the busucretzze-heenoppiniei and "optimize- then-discretize" approaches. Our work addresses these gaps by:
1. Gradient Equivalence: Establishing conditions under which the gradients computed by the two methods are equivalent, as demonstrated in Proposition] and Proposition 2 emphasizing the role of smooth activations.
53

---

Published as a conference paper at ICLR 2025
2. NTK Analysis:
Providing rigorous conditions for the well-definedness of the Neural ODE NTK, demonstrating its strict positive definiteness (SPD) under suitable activation function properties, stated in Theorem 2 and Corollary !I
3. Global Convergence: rameterized Neural ODEs, bridging the gap between discrete and continuous-depth models, as outlined in Theorem 3
Extending global convergence guarantees for gradient descent in overpa-
54


Reviews:
[{"review_id": "l6eijrHrVx", "reviewer": "Reviewer_C5fC", "summary": "This paper presents several theoretical results about Neural ODEs, built around sufficient regularity of activations:\n- It proves that \"optimize then discretize\" setting, built upon the use of classical ODE solvers, is well-posed under the assumption of Lipschitz-continuous activations (proposition 1).\n- It proves that \"discretize-then-optimize\" setting, built as an Euler scheme discretization of the ODE by growing the depth of a residual neural network, is also well posed in the sense that the solution of the discretization scheme converges to the actual solution of the ODE in the infinite depth limit, if the activations and there derivatives are lipschitz-continuous (proposition 2). It has to be noted that authors gives explicit convergence rates.\n- It shows that under the \"discretize-then-optimize\" setting with fixed depth L, the associated residual network converges in infinite width limit to a fixed kernel with explicit form, under the assumption of Lipschitz continuity of activations (proposition 3).\n- It shows that under the assumption of Lipschitz continuity of activations, the infinite depth and width limits commute, implying that the neural ODE converges to a Gaussian process with covariance given as the infinite depth limit of the covariance of proposition 3 (theorem 1).\n- Similar results are shown for the convergence of the NTK (proposition 4 and theorem 2)\n- In section 5, authors show that the resulting limiting covariance of the gaussian process as well as NTK are Semi-positive definite (lemma 3, proposition 5, corrolary 1), resulting in derivation of explicit convergences rates (theorem 3).\nFinally empirical evidences of the theoretical results are presented in section 6.", "questions": "I don't think it's quite right (or at least I don't understand exactly why, but I'm not a specialist in Neural ODEs) that the title of the paper focuses on the impact of activations and their regularities, whereas most of your results use these as hypotheses for theorem proving. In my opinion, you should assume a more theoretical title (or present more results comparing activations, but it seems to me that this is a slightly different job than what you present, which is amply sufficient from a theoretical point of view).", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- precise, rigorously demonstrated theoretical results\n- clear presentation\n- extension of fundamental NTK results to neural ODEs, using classical ODE analytical tools (Picard theorem, Groenwald theorem etc).", "weaknesses": "- Typo : the index of the last h in equation 1 should probably be T instead of t\n- Figure 1 : a,b,c should rather be in log-log scale to better appreciate the slope comparison with 1/L.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "K1Lb9aogZu", "reviewer": "Reviewer_E13D", "summary": "This paper investigates the role of activation functions in the training dynamics of Neural ODEs, particularly focusing on the convergence and stability properties of two key methods: \"Optimize-then-discretize\" and \"Discretize-then-optimize.\" The authors establish theoretical guarantees for the convergence of both methods under certain conditions on the activation functions, and they provide numerical experiments to support their findings.", "questions": "1.\nCould you provide a numerical experiment that specifically highlights the convergence problem of the \"Discretize-then-optimize\" method when using an activation function that does not meet the Lipschitz continuity requirements? This would strengthen the paper by providing empirical support for your theoretical analysis.\n\n2.\nHow does the choice of numerical ODE solver in the \"Optimize-then-discretize\" method impact the accuracy of the gradient computation and the overall training dynamics? Exploring the sensitivity of the method to different solvers would be valuable.\n\n3.\nWhat are the practical implications of your findings for the choice of activation functions in real-world applications of Neural ODEs? Providing guidance on activation function selection based on the desired trade-off between theoretical guarantees and empirical performance would be beneficial.\n\nAddressing these points would significantly improve the paper and provide a more comprehensive understanding of the impact of activation functions on Neural ODE training.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Theoretical Analysis: The paper offers a rigorous mathematical analysis of the well-posedness of Neural ODE solutions and their gradients, focusing on the impact of activation function smoothness. The authors clearly demonstrate the importance of Lipschitz continuity for both the activation function and its derivative to ensure stable and accurate gradient computation.\n\nNTK Analysis: The paper extends the concept of the Neural Tangent Kernel (NTK) to the continuous setting of Neural ODEs. This is a valuable contribution, as it provides insights into the training dynamics of these models. The authors show that the NTK converges to a deterministic limit under certain conditions, which helps explain the convergence properties of gradient descent in Neural ODE training.\n\nGlobal Convergence Proof: The authors establish the global convergence of Neural ODEs under gradient descent in the overparameterized regime, under assumptions about the activation function's nonlinearity and non-polynomial nature. This theoretical result provides a strong foundation for understanding the optimization behavior of Neural ODEs.\n\nEmpirical Validation: The paper presents several numerical experiments to validate the theoretical findings. These experiments demonstrate the convergence of Neural ODE outputs and gradients to their finite-depth ResNet approximations, the behavior of the NTK, and the convergence properties of Neural ODEs under gradient descent.", "weaknesses": "Lack of Demonstration for \"Discretize-then-optimize\" Convergence Issue: While the paper theoretically analyzes the convergence problem of \"Discretize-then-optimize\" with non-smooth activation functions, it does not provide a dedicated numerical experiment to showcase this issue. This would be valuable for empirically confirming the theoretical claims and for understanding the practical implications of using non-smooth activations with this method.\n\nThe lack of a numerical demonstration of the \"Discretize-then-optimize\" convergence problem weakens the overall impact of the paper. Without empirical evidence, the theoretical claims regarding the convergence issue remain unsubstantiated, and the practical implications of using non-smooth activation functions with this method remain unclear.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "yRXK9XBGPD", "reviewer": "Reviewer_Zcba", "summary": "The authors studied activation functions in the training dynamics of neural ODEs. They showed that the choice of activation function has significant impacts on the dynamic and asymptotic behavior. Additionally, they propose a new mathematical framework to study continuous models from the perspective of approximation theory. Finally, the authors conduct experiments to empirically validate their theoretical results.", "questions": "No questions.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-written. The ideas presented are novel and represent a significant contribution in the context of neural ODEs. The main strength is the theoretical foundation of the work.", "weaknesses": "The proofs are generally well-written and clear, but since they are in the appendix, the flow of the paper's main narrative can sometimes be interrupted.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ewI9Ckitew", "reviewer": "Reviewer_Hphg", "summary": "This paper studies the convergence of gradient descent for training neural ODEs from the neural tangent kernel perspectives, with experimental results validating the theoretical claims.", "questions": "See weakness.", "limitations": "", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "(1) The paper extends NTK theory from general regression problems with ReLU nets to Neural ODEs with ResNet architectures. The analysis and results are novel for Neural ODEs.\n\n(2) The empirical results validate the theorems.", "weaknesses": "(1) Some prior work has extended NTK theory to analyze CNNs (for image classification tasks), PINNs, Transformers, and various gradient-based algorithms beyond vanilla gradient descent. However, related references are missing. \n\n(2) SGD is more commonly used in practice, making a convergence analysis of SGD particularly valuable and interesting.\n\n(3) The results are specific to ODEs of the form given in Equation (2). Theorem 2 does not hold for PINNs if the underlying PDE is nonlinear. How would Theorem 2 apply if the ODE were not limited to this special case?\n\n(4) This paper claims that the nonlinear and non-polynomial nature of the activation function is crucial for convergence. However, I believe the convergence results (Theorem 3) may also hold with polynomial activation functions, provided that Condition (1) of Assumption 1 is further strengthened to ensure linear independence.\n\n(5) From the title, the topic of the paper is exploring the smoothness of activation functions in neural ODEs. However, after reading your paper, your main contribution is developing the global convergence of GD in training neural ODEs. The title is quite strange...", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Oral)