A Lagrangian generalization of time-reversible Born-Oppenheimer molecular dynamics Niklasson et al. [Phys. Rev. Lett. 97, 123001 (2006)10.1103/PhysRevLett.97.123001] is proposed. The formulation enables the application of higher-order symplectic or geometric integration schemes that are stable and energy conserving even under incomplete self-consistency convergence. It is demonstrated how the accuracy is improved by over an order of magnitude compared to previous formulations at the same level of computational cost. The proposed Lagrangian includes extended electronic degrees of freedom as auxiliary dynamical variables in addition to the nuclear coordinates and momenta. While the nuclear degrees of freedom propagate on the Born-Oppenheimer potential energy surface, the extended auxiliary electronic degrees of freedom evolve as a harmonic oscillator centered around the adiabatic propagation of the self-consistent ground state.