Neural Ordinary Differential Equations

Here is one of the most beautiful ideas in modern deep learning. You know how a Residual Network computes its output? Each layer takes the current hidden state and adds a small update: h(t+1) = h(t) + f(h(t), theta_t). The output of layer 6 is the input of layer 1, plus six small increments stacked on top of each other. Now squint at that formula. Does it remind you of anything from calculus?

It should. That is Euler's method — the simplest way to numerically solve a differential equation, with step size 1. Every ResNet you have ever trained was secretly doing numerical integration. The authors of this paper asked a question that seems obvious in hindsight: what if we took the step size to zero and the number of layers to infinity? Instead of a staircase of discrete layers, you get a smooth ramp — a continuous transformation governed by an ordinary differential equation.

This is not just a mathematical curiosity. It unlocks three practical superpowers. First, constant memory: you do not need to store intermediate layer activations for backprop because the adjoint method recomputes them on the fly. Second, adaptive computation: the ODE solver automatically decides how many steps to take based on how complex the dynamics are — simple inputs get fewer steps, hard inputs get more. Third, parameter efficiency: one set of parameters shared across continuous time, instead of separate weights per layer.

Let us be precise about what a Neural ODE actually computes. Instead of specifying a discrete sequence of hidden layers, we parameterize the derivative of the hidden state using a neural network f. Given an input h(0), the output h(T) is the solution to an initial value problem: start at h(0) and follow the velocity field defined by f for time T. The answer is whatever state you end up at.

Think of it this way: the neural network f does not transform data directly. Instead, it defines a velocity at every point in space and time. If you are at position h and the clock reads t, then f(h, t, theta) tells you which direction to move and how fast. The ODE solver follows these instructions, tracing out a smooth curve from h(0) to h(T). The same parameters theta define the velocity everywhere — unlike a ResNet where each layer has its own independent parameters.

Why does this matter? Because the dynamics function f is shared across all times, a Neural ODE is dramatically more parameter-efficient. A 6-block ResNet has 6 separate parameter sets. A Neural ODE with equivalent representational power has just one. And because the ODE solver chooses how many steps to take, the effective depth is adaptive — not fixed at architecture design time.

The ODE solver is the forward pass of a Neural ODE. It is not an approximation or a hack — it is a rigorous numerical method for computing the exact solution to the differential equation (up to solver tolerance). Different solvers offer different trade-offs between accuracy, speed, and stability, and choosing the right one is critical for both training and inference.

The paper uses Dormand-Prince (DOPRI5), a 5th-order adaptive step-size Runge-Kutta method that has been the workhorse of scientific computing for decades. Here is how it works: at each step, it computes two estimates of the next state — one 4th-order and one 5th-order. The difference between them estimates the local truncation error. If the error is too large, the step is rejected and retried with a smaller step size. If the error is tiny, the next step is made larger. This means the solver automatically takes big steps through smooth regions and tiny steps through turbulent ones.

Why not just use Euler's method with many small steps? Because Euler is 1st-order — halving the step size only halves the error. DOPRI5 is 5th-order — halving the step size reduces the error by a factor of 32. For the same accuracy, DOPRI5 needs far fewer function evaluations, which translates directly to fewer neural network forward passes and faster training.

Here is the problem that almost killed Neural ODEs before they were born. To train a neural network, you need gradients. Standard backpropagation computes gradients by storing all intermediate activations and walking backward through them. But an ODE solver might take hundreds of internal steps during a single forward pass. Storing all those intermediate states would require enormous memory — potentially more than the entire GPU.

The adjoint sensitivity method solves this elegantly. It is a technique from optimal control theory (1960s!) that computes exact gradients with O(1) memory — constant memory regardless of how many solver steps are taken. The idea is beautiful in its simplicity: instead of storing the forward trajectory and walking backward through it, we define a new ODE — the adjoint equation — that computes the gradient by integrating backward in time.

Define the adjoint state a(t) = dL/dh(t) — how the loss changes with the hidden state at time t. The key insight: a(t) satisfies its own ODE: da/dt = -a(t)^T (df/dh). We can integrate this backward from t=T to t=0, simultaneously recomputing the forward state h(t) on the fly and accumulating parameter gradients. No stored activations. No checkpointing. Just solve three coupled ODEs backward: the adjoint, the state, and the parameter gradient.

One of the paper's most impactful contributions was not just theoretical — it was a new way to build generative models. Traditional normalizing flows transform a simple distribution (like a Gaussian) into a complex one through a chain of carefully designed invertible layers. The catch? Each layer must be invertible with an efficiently computable Jacobian determinant. This severely restricts what architectures you can use — coupling layers, autoregressive transforms, and other constrained designs.

Continuous Normalizing Flows (CNFs) blow this restriction wide open. Since an ODE is automatically invertible (just reverse time), you can use any neural network architecture for the dynamics function f. No coupling layers. No autoregressive constraints. No architectural gymnastics to ensure invertibility. The ODE solver handles it all.

But wait — how do you compute the density? In discrete normalizing flows, you need the log-determinant of each layer's Jacobian, which costs O(D^3). The magic of CNFs comes from the instantaneous change of variables formula: the log-density changes at a rate given by just the trace of the Jacobian — O(D) instead of O(D^3). For D=1000 dimensions, that is a billion times cheaper. And the trace can be estimated stochastically using the Hutchinson estimator: just compute v^T (df/dz) v for a random vector v. One Jacobian-vector product instead of the full Jacobian matrix.

A Neural ODE learns a vector field over the hidden state space. At every point, the network defines a velocity: which direction and how fast the state should move. The ODE solver traces trajectories through this field, and the collection of all possible trajectories is called the phase portrait — one of the most powerful tools in dynamical systems theory.

Understanding the phase portrait reveals what the model has actually learned. You can see fixed points (where the dynamics stop — these are attractors or repellers), limit cycles (closed loops — periodic behavior), and separatrices (boundaries between regions with qualitatively different behavior). A classification Neural ODE, for instance, learns a vector field that routes different classes into different regions of space.

But here is a fundamental constraint that limits what Neural ODEs can represent: by the Picard-Lindelof uniqueness theorem, two ODE trajectories starting from different initial conditions cannot cross in the same phase space. If two data points start at different positions, their trajectories will never pass through the same point at the same time. This means Neural ODEs preserve the topology of the input — they cannot untangle interleaved spirals or separate concentric circles. This is a real limitation, and it drove the development of Augmented Neural ODEs (which add extra dimensions to allow effective crossing).

The paper introduces the Latent ODE model — perhaps the most practically impactful contribution for real-world applications. The idea: combine a Variational Autoencoder (VAE) with a Neural ODE decoder. An encoder processes observations (potentially at completely irregular times) to produce a distribution over latent initial states. The Neural ODE decoder then evolves this state forward to generate predictions at any desired time.

Why is this so powerful? Consider medical data. A patient visits the hospital at random times. Their blood pressure is measured at 9am, temperature at 11am, blood test at 2pm tomorrow, another visit next week. Traditional RNNs treat these as a fixed-step sequence, discretizing time into bins and padding with zeros. But a Latent ODE models the patient's state as evolving continuously. Between observations, the hidden state follows the learned dynamics — no padding, no binning, no information loss. When an observation arrives, the state is updated with a GRU-like jump.

The encoder is an ODE-RNN: between observations, the hidden state evolves continuously according to a learned ODE. At observation times, it receives a discrete update. This explicitly models the passage of time between events — unlike a standard RNN where the hidden state is frozen between observations. The decoder then samples z(0) from the encoded distribution and integrates forward, producing a smooth trajectory that can be queried at any time point.

Here is something no fixed-depth network can do: trade accuracy for speed at inference time by turning a single dial. In a Neural ODE, that dial is the solver's error tolerance. Tighter tolerances (e.g., 10^-5) mean more function evaluations (NFE), more computation, and higher accuracy. Looser tolerances (e.g., 10^-1) mean fewer evaluations, faster inference, and a controlled accuracy sacrifice.

The Number of Function Evaluations (NFE) is the effective depth of a Neural ODE. A ResNet-6 always has depth 6, whether the input is trivial or complex. A Neural ODE might use 8 evaluations for a simple input and 45 for a complex one — the adaptive solver automatically allocates more computation where it is needed. This is input-dependent adaptive computation without any architecture changes, conditional paths, or early-exit mechanisms. The solver just does it.

An interesting observation from the paper: NFE tends to increase during training as the model learns more complex dynamics. At initialization, the dynamics are nearly zero (random weights produce small outputs), so the solver needs few steps. As training progresses and the model learns to solve the task, the dynamics become more complex and the solver takes more steps. This can be problematic — later work on regularized Neural ODEs addresses this by penalizing the dynamics magnitude to keep NFE manageable.

Training Neural ODEs requires careful attention to solver choice, tolerance settings, and regularization. The paper demonstrates results on three very different tasks — image classification, density estimation, and time series modeling — each requiring different training strategies. Let us look at the concrete numbers and learn what makes Neural ODE training succeed or fail.

Neural ODEs are elegant but have real limitations. Understanding these limitations is important not just for using the models correctly, but because they drove a wave of follow-up work that significantly expanded the continuous-depth model family. Every limitation spawned at least one new paper with a clever solution.

Neural ODEs found their niche in domains where data is naturally continuous, where physics provides dynamical constraints, and where adaptive computation is valuable. They shine brightest when the underlying process you are modeling is actually governed by differential equations — which turns out to be surprisingly many real-world problems.

In scientific computing, Neural ODEs learn physical dynamics directly from data: planetary orbits, molecular interactions, fluid flows. Hamiltonian Neural Networks (a direct descendant) guarantee energy conservation by parameterizing the Hamiltonian rather than the vector field. In healthcare, Latent ODEs model patient disease trajectories from sparse clinical visits, predicting outcomes and treatment responses with proper uncertainty. In generative modeling, CNFs enabled the probability flow ODE that powers modern diffusion models.

The Neural ODEs paper did more than introduce a new model — it opened a new research field at the intersection of deep learning and differential equations. Before 2018, these two communities rarely talked to each other. After the paper, a flood of work connected neural networks to dynamical systems, optimal control theory, numerical analysis, and mathematical physics.