Neural ordinary differential equation

A neural ordinary differential equation, commonly abbreviated as neural ODE, is a differential equation whose vector field is parameterized by a neural network. It represents the evolution of a hidden state as a continuous function of an independent variable, which is usually interpreted as time or network depth. The standard initial-value formulation is

[ \frac{d\mathbf{z}(t)}{dt}

f_{\theta}!\left(\mathbf{z}(t),t\right), \qquad \mathbf{z}(t_0)=\mathbf{z}_0, ]

where (\mathbf{z}(t)) is the state, (f_{\theta}) is a neural network with parameters (\theta), and (\mathbf{z}_0) is an initial condition. The output at (t_1) is determined by integrating the learned vector field:

[ \mathbf{z}(t_1)

\mathbf{z}(t_0) + \int_{t_0}^{t_1} f_{\theta}!\left(\mathbf{z}(t),t\right),dt. ]

Neural ODEs are continuous-depth models rather than networks containing a fixed sequence of discrete transformations. Their evaluation therefore depends on a numerical ordinary differential equation solver, which approximates the state trajectory to a prescribed numerical tolerance.

Historical development

The mathematical foundations of neural ODEs combine earlier work on numerical integration, dynamical systems, optimal control, and trainable differential equations. Leonhard Euler's discretization of continuous dynamics established the direct relationship between differential equations and iterative state updates, while Carl Runge and Martin Kutta contributed the higher-order integration schemes now used in many numerical solvers. Lev Pontryagin and his collaborators developed the adjoint equations of optimal control, providing a mathematical basis for differentiating objectives constrained by differential equations.

The modern continuous-depth formulation was introduced in 2018 by Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, David Duvenaud, and You Watanabe. Their formulation treated the transformation performed by a residual network as the numerical discretization of a continuous hidden-state trajectory. It also applied an adjoint-state calculation to obtain parameter gradients without retaining every internal state produced by the forward integration.

Subsequent work connected neural ODEs to continuous-time generative models and irregularly sampled temporal data. Yulia Rubanova, Ricky T. Q. Chen, and David Duvenaud developed the latent ODE framework, in which a probabilistic latent state follows learned continuous dynamics between observations. Will Grathwohl, Ricky T. Q. Chen, Jesse Bettencourt, Ilya Sutskever, and David Duvenaud developed the free-form continuous dynamics framework for scalable continuous normalizing flows.

Relation to residual networks

A residual neural network updates a hidden state through transformations of the form

[ \mathbf{z}_{k+1}

\mathbf{z}k + h f{\theta_k}(\mathbf{z}_k), ]

where (k) indexes the network layer and (h) is a step size. This expression has the same form as a forward Euler discretization of an ordinary differential equation:

[ \frac{\mathbf{z}_{k+1}-\mathbf{z}_k}{h}

f_{\theta_k}(\mathbf{z}_k). ]

Under a continuous-depth interpretation, the layer index becomes a continuous coordinate (t), while the residual transformation becomes the derivative (d\mathbf{z}/dt). Parameter sharing across depth corresponds to an autonomous vector field (f_{\theta}(\mathbf{z})). Allowing explicit dependence on (t) yields a nonautonomous system whose dynamics can change over the integration interval.

This correspondence does not make every residual network identical to a neural ODE. A finite residual network implements a particular discrete composition, whereas a neural ODE defines a continuous initial-value problem whose computed output depends on the chosen integration method and its error control. The correspondence instead identifies the residual update as one numerical approximation to continuous hidden-state dynamics.

Forward evaluation

Forward evaluation consists of solving the initial-value problem from (t_0) to (t_1). An adaptive solver selects internal evaluation points according to an estimate of local truncation error. Consequently, the number of evaluations of (f_{\theta}) is determined by the learned dynamics, the integration interval, the numerical method, and the specified tolerances rather than by a fixed architectural depth.

For sufficiently regular vector fields, the Picard–Lindelöf theorem establishes local existence and uniqueness of the state trajectory. A neural network using smooth or locally Lipschitz activation functions generally defines a locally Lipschitz vector field with respect to the state. Global existence additionally depends on the growth of the vector field, since rapidly increasing dynamics can produce finite-time divergence or trajectories that are expensive to resolve numerically.

The numerical computation is part of the realized model. Distinct solvers approximate the same mathematical trajectory with different truncation errors, and adaptive algorithms may evaluate the vector field at different times. Training and inference therefore involve both approximation error from statistical estimation and discretization error from numerical integration.

Differentiation and the adjoint state

Let a scalar loss depend on the terminal state:

[ L = \ell!\left(\mathbf{z}(t_1)\right). ]

The adjoint state is defined as

[ \mathbf{a}(t)

\frac{\partial L}{\partial \mathbf{z}(t)}. ]

For a smooth vector field, it satisfies the reverse-time differential equation

[ \frac{d\mathbf{a}(t)}{dt}

\mathbf{a}(t)^{\mathsf T} \frac{\partial f_{\theta}}{\partial \mathbf{z}} \left(\mathbf{z}(t),t\right), ]

with terminal condition

[ \mathbf{a}(t_1)

\frac{\partial \ell}{\partial \mathbf{z}(t_1)}. ]

The parameter gradient is

[ \frac{dL}{d\theta}

\int_{t_1}^{t_0} \mathbf{a}(t)^{\mathsf T} \frac{\partial f_{\theta}}{\partial \theta} \left(\mathbf{z}(t),t\right),dt. ]

This construction is the continuous analogue of backpropagation. It expresses reverse-mode differentiation as another differential equation and permits an implementation that reconstructs states during reverse integration instead of storing the complete forward solver trace.

The continuous adjoint and direct differentiation through a numerical solver are not numerically identical. Direct differentiation computes derivatives of the discretized algorithm, whereas the classical continuous adjoint first differentiates the differential equation and then discretizes the resulting adjoint system. Reconstruction error, adaptive step selection, and nonreversible numerical integration can cause the two gradients to differ. Checkpointed adjoint methods retain selected forward states and recompute intermediate segments, placing memory use and numerical consistency between full trajectory storage and unrestricted reverse reconstruction.

Expressivity and topology

A neural ODE generates a continuous flow map when its vector field produces unique solutions. For any fixed integration interval, such a flow continuously deforms the initial state space without allowing distinct trajectories to intersect at the same time. This property imposes topological restrictions that do not apply to arbitrary feedforward networks.

In one dimension, for example, a unique continuous flow preserves the ordering of initial states. A transformation that reverses this ordering cannot be represented directly by an ordinary one-dimensional neural ODE over a finite interval. More generally, a flow cannot realize a terminal map that requires trajectories to cross or that changes topology through discontinuous identification.

Augmented neural ODEs address this restriction by embedding the state in a higher-dimensional space before integration. Trajectories that would intersect in the observed coordinates can remain disjoint in the augmented dimensions. The terminal state is then projected back into the output space, allowing the projected transformation to represent mappings unavailable to the unaugmented flow.

Continuous normalizing flows

A continuous normalizing flow uses an invertible neural ODE to transform a probability distribution. If (\mathbf{z}(t)) follows the learned dynamics, its log-density changes according to the instantaneous change-of-variables equation:

[ \frac{d}{dt}\log p!\left(\mathbf{z}(t)\right)

\operatorname{tr} \left( \frac{\partial f_{\theta}}{\partial \mathbf{z}} \right). ]

Integration gives

[ \log p!\left(\mathbf{z}(t_1)\right)

\log p!\left(\mathbf{z}(t_0)\right)

\int_{t_0}^{t_1} \operatorname{tr} \left( \frac{\partial f_{\theta}}{\partial \mathbf{z}} \right)dt. ]

This relation replaces the determinant of a finite-dimensional transformation Jacobian with the trace of the vector-field Jacobian accumulated along the trajectory. Exact trace computation scales with state dimension, while stochastic trace estimators approximate it through Jacobian–vector products. The resulting model remains an instance of a normalizing flow because the unique ODE flow defines an invertible transformation between the initial and terminal states.

Time-series models

Neural ODEs provide a continuous latent-state representation for observations recorded at irregular times. An encoder maps available observations to a distribution over an initial latent state, after which the learned differential equation propagates that state to each observation time. A decoder or observation model then relates the latent trajectory to measured data.

This structure separates continuous evolution from discrete observation events. In latent ODE models, uncertainty in the initial state induces uncertainty over the entire trajectory. In ODE recurrent networks, recurrent updates incorporate observations while ODE integration evolves the hidden state between them. Related neural controlled differential equations replace autonomous latent evolution with dynamics driven by a continuous interpolation of the observed input path.

Numerical characteristics

The computational cost of a neural ODE is commonly measured by the number of evaluations of its vector field. Smooth trajectories with moderate curvature can be integrated using comparatively large steps, whereas oscillatory or stiff dynamics require smaller steps or methods designed for stiff equations. Since the vector field is learned, optimization can alter not only the represented transformation but also the difficulty of numerical integration.

Adaptive solvers introduce data-dependent computation because different initial conditions can require different internal step sequences. The solver's error tolerance controls numerical accuracy rather than statistical generalization, and tighter tolerances usually increase the number of vector-field evaluations. Discontinuities in the vector field or in external forcing require event handling or hybrid formulations because a standard smooth ODE does not by itself encode instantaneous state jumps.

The continuous-depth description also does not imply an exact continuous computation in an implementation. Every finite computation uses a discretization, and the resulting approximation inherits the stability properties of the selected numerical method. The distinction between the differential equation and its solver remains central when interpreting trained dynamics, gradient calculations, and reproducibility across software systems.