  Physical systems can often be described via a continuous-time dynamical
system. In practice, the true system is often unknown and has to be learned
from measurement data. Since data is typically collected in discrete time, e.g.
by sensors, most methods in Gaussian process (GP) dynamics model learning are
trained on one-step ahead predictions. This can become problematic in several
scenarios, e.g. if measurements are provided at irregularly-sampled time steps
or physical system properties have to be conserved. Thus, we aim for a GP model
of the true continuous-time dynamics. Higher-order numerical integrators
provide the necessary tools to address this problem by discretizing the
dynamics function with arbitrary accuracy. Many higher-order integrators
require dynamics evaluations at intermediate time steps making exact GP
inference intractable. In previous work, this problem is often tackled by
approximating the GP posterior with variational inference. However, exact GP
inference is preferable in many scenarios, e.g. due to its mathematical
guarantees. In order to make direct inference tractable, we propose to leverage
multistep and Taylor integrators. We demonstrate how to derive flexible
inference schemes for these types of integrators. Further, we derive tailored
sampling schemes that allow to draw consistent dynamics functions from the
learned posterior. This is crucial to sample consistent predictions from the
dynamics model. We demonstrate empirically and theoretically that our approach
yields an accurate representation of the continuous-time system.
