  The article introduces a method to learn dynamical systems that are governed
by Euler--Lagrange equations from data. The method is based on Gaussian process
regression and identifies continuous or discrete Lagrangians and is, therefore,
structure preserving by design. A rigorous proof of convergence as the distance
between observation data points converges to zero is provided. Next to
convergence guarantees, the method allows for quantification of model
uncertainty, which can provide a basis of adaptive sampling techniques. We
provide efficient uncertainty quantification of any observable that is linear
in the Lagrangian, including of Hamiltonian functions (energy) and symplectic
structures, which is of interest in the context of system identification. The
article overcomes major practical and theoretical difficulties related to the
ill-posedness of the identification task of (discrete) Lagrangians through a
careful design of geometric regularisation strategies and through an exploit of
a relation to convex minimisation problems in reproducing kernel Hilbert
spaces.
