  Among the single-trajectory Gaussian-based methods for solving the
time-dependent Schr\"{o}dinger equation, the variational Gaussian approximation
is the most accurate one. In contrast to Heller's original thawed Gaussian
approximation, it is symplectic, conserves energy exactly, and may partially
account for tunneling. However, the variational method is also much more
expensive. To improve its efficiency, we symmetrically compose the second-order
symplectic integrator of Faou and Lubich and obtain geometric integrators that
can achieve an arbitrary even order of convergence in the time step. We
demonstrate that the high-order integrators can speed up convergence
drastically compared to the second-order algorithm and, in contrast to the
popular fourth-order Runge-Kutta method, are time-reversible and conserve the
norm and the symplectic structure exactly, regardless of the time step. To show
that the method is not restricted to low-dimensional systems, we perform most
of the analysis on a non-separable twenty-dimensional model of coupled Morse
oscillators. We also show that the variational method may capture tunneling
and, in general, improves accuracy over the non-variational thawed Gaussian
approximation.
