  Most physical processes posses structural properties such as constant
energies, volumes, and other invariants over time. When learning models of such
dynamical systems, it is critical to respect these invariants to ensure
accurate predictions and physically meaningful behavior. Strikingly,
state-of-the-art methods in Gaussian process (GP) dynamics model learning are
not addressing this issue. On the other hand, classical numerical integrators
are specifically designed to preserve these crucial properties through time. We
propose to combine the advantages of GPs as function approximators with
structure preserving numerical integrators for dynamical systems, such as
Runge-Kutta methods. These integrators assume access to the ground truth
dynamics and require evaluations of intermediate and future time steps that are
unknown in a learning-based scenario. This makes direct inference of the GP
dynamics, with embedded numerical scheme, intractable. Our key technical
contribution is the evaluation of the implicitly defined Runge-Kutta transition
probability. In a nutshell, we introduce an implicit layer for GP regression,
which is embedded into a variational inference-based model learning scheme.
