Neural ODEs and Continuous-Depth Networks

A residual network with update $h_{t+1} = h_t + \frac{1}{L} f_\theta(h_t, t/L)$ for $t = 0, 1, \ldots, L-1$ is the Euler discretization of the initial value problem:

$\frac{dh}{dt} = f_\theta(h(t), t), \quad h(0) = x, \quad t \in [0, 1]$

with step size $\Delta t = 1/L$. In the limit $L \to \infty$, the discrete trajectory $\{h_0, h_1, \ldots, h_L\}$ converges to the continuous solution $h(t)$ (under Lipschitz conditions on $f_\theta$).

The output of the network is $h(1) = h(0) + \int_0^1 f_\theta(h(t), t) \, dt$.

Each ResNet layer adds a small correction to the hidden state. In the continuous limit, these corrections become a vector field that flows the input through a smooth trajectory. The network's "depth" becomes a continuous time variable. Deeper networks correspond to longer integration times, and the network learns the vector field $f_\theta$ that transforms inputs into useful representations.

This perspective explains why ResNets work: the skip connection $h_{t+1} = h_t + f(h_t)$ is not just a gradient-flow trick. It makes each layer an incremental transformation, and incremental transformations compose into smooth, invertible maps (under mild conditions). This is why ResNets can be very deep without the representation collapsing.

It also enables replacing the fixed $L$-layer architecture with an adaptive ODE solver that chooses its own step size. Easy inputs get fewer steps; hard inputs get more. This is adaptive computation depth without architectural changes.

The continuous limit requires the dynamics $f_\theta$ to be Lipschitz continuous. If $f_\theta$ is not Lipschitz (e.g., if it has sharp discontinuities), the ODE may not have a unique solution, and the convergence of Euler's method is not guaranteed. In practice, standard neural network architectures with smooth activations (tanh, softplus) satisfy this condition. ReLU technically does not (it is not differentiable at 0), but works in practice.

To compute $\frac{dL}{d\theta}$ for a neural ODE, define the adjoint state $a(t) = \frac{dL}{dh(t)}$. The adjoint satisfies a backward ODE:

$\frac{da}{dt} = -a(t)^\top \frac{\partial f}{\partial h}(h(t), t, \theta)$

integrated backwards from $t = T$ to $t = 0$ with initial condition $a(T) = \frac{dL}{dh(T)}$. The parameter gradient is:

$\frac{dL}{d\theta} = -\int_0^T a(t)^\top \frac{\partial f}{\partial \theta}(h(t), t, \theta) \, dt$

Memory cost: $O(1)$ in depth (constant, regardless of the number of ODE solver steps), compared to $O(L)$ for standard backpropagation through $L$ layers.

Standard backprop stores all intermediate activations $h_0, h_1, \ldots, h_L$ to compute gradients, costing $O(L)$ memory. The adjoint method avoids this by solving the adjoint ODE backwards in time, recomputing $h(t)$ on the fly by integrating the forward ODE backwards. This trades memory for compute: you solve two ODEs (forward and backward) instead of storing all activations.

This is the continuous analog of activation checkpointing, taken to its logical extreme.

Constant-memory backpropagation enables training of very deep (or continuous-depth) networks without running out of GPU memory. For standard ResNets with $L = 100$ layers, the memory saving is $\sim 100\times$. For neural ODEs where the solver may take thousands of steps, the saving is even larger.

The adjoint method is not new. It was developed in optimal control theory in the 1960s (Pontryagin's maximum principle). Neural ODEs brought it to deep learning.

The backward ODE recomputation of $h(t)$ introduces numerical error. If the forward and backward solvers use different discretizations or if the dynamics are chaotic, the recomputed $h(t)$ can diverge from the original, causing gradient inaccuracy. In practice, this is mitigated by using the same adaptive solver in both directions, but it remains a source of subtle bugs. Checkpointed approaches (solving forward, saving a few checkpoints, recomputing between them) offer a middle ground.