Numerically solving an ODE system is a computationally intensive task, and often dominates the cost of one unnormalized posterior density evaluation in Bayesian ODE models. Computational requirements are amplified by the fact that MCMC inference typically requires a large number of these ODE solutions, and this number can be several orders of magnitude larger than the posterior effective sample size (ESS) due to the need to warm-up the sampler and high auto-correlation of intermediate draws. The ratio of ESS to the number of unnormalized posterior density evaluations can be drastically increased by using gradient-based MCMC methods such as Hamiltonian Monte Carlo (HMC) [1]}, and using the gradient information is essential in the case of high-dimensional parameter spaces or complex posterior geometries in order to achieve good sampling performance. In modern statistical software such as Stan [2]}, PyMC [3]}, and Turing.jl [4]}, gradients are computed using automatic differentiation [5]}, [6]}. This poses a challenge for ODE models, for which computing the gradient is computationally significantly more demanding and subtle than the plain likelihood evaluation.
