  We study the design and implementation of numerical methods to solve the
generalized Langevin equation (GLE) focusing on canonical sampling properties
of numerical integrators. For this purpose, we cast the GLE in an extended
phase space formulation and derive a family of splitting methods which
generalize existing Langevin dynamics integration methods. We show exponential
convergence in law and the validity of a central limit theorem for the Markov
chains obtained via these integration methods, and we show that the dynamics of
a suggested integration scheme is consistent with asymptotic limits of the
exact dynamics and can reproduce (in the short memory limit) a superconvergence
property for the analogous splitting of underdamped Langevin dynamics. We then
apply our proposed integration method to several model systems, including a
Bayesian inference problem. We demonstrate in numerical experiments that our
method outperforms other proposed GLE integration schemes in terms of the
accuracy of sampling. Moreover, using a parameterization of the memory kernel
in the GLE as proposed by Ceriotti et al [9], our experiments indicate that the
obtained GLE-based sampling scheme outperforms state-of-the-art sampling
schemes based on underdamped Langevin dynamics in terms of robustness and
efficiency.
