  Gaussian latent variable models are a key class of Bayesian hierarchical
models with applications in many fields. Performing Bayesian inference on such
models can be challenging as Markov chain Monte Carlo algorithms struggle with
the geometry of the resulting posterior distribution and can be prohibitively
slow. An alternative is to use a Laplace approximation to marginalize out the
latent Gaussian variables and then integrate out the remaining hyperparameters
using dynamic Hamiltonian Monte Carlo, a gradient-based Markov chain Monte
Carlo sampler. To implement this scheme efficiently, we derive a novel adjoint
method that propagates the minimal information needed to construct the gradient
of the approximate marginal likelihood. This strategy yields a scalable
differentiation method that is orders of magnitude faster than state of the art
differentiation techniques when the hyperparameters are high dimensional. We
prototype the method in the probabilistic programming framework Stan and test
the utility of the embedded Laplace approximation on several models, including
one where the dimension of the hyperparameter is $\sim$6,000. Depending on the
cases, the benefits can include an alleviation of the geometric pathologies
that frustrate Hamiltonian Monte Carlo and a dramatic speed-up.
