Many dynamic statistical models are characterized by a joint distribution with potentially high-dimensional nuisance parameters. To avoid the computational cost of numerically marginalizing over such parameters, we can instead use the Laplace approximation for the integral. This approximation is exact if and only if the distribution is Gaussian in the nuisance parameters; therefore, its effectiveness depends intimately on the true shape of the distribution. Here, we review a probabilistic numeric tool for assessing the quality of the Laplace approximation. This diagnostic tool recasts the approximation problem in a Bayesian framework, modelling the distribution with a Gaussian process prior and evaluating it at a coarse grid of points to obtain a posterior from which inferences can be made. We will discuss approaches for optimizing the tool, choosing evaluation points and other parameters to maximize our ability to detect non-Gaussianity. Our methods will be demonstrated on a variety of distributions, showing the effects of shape and dimension on the optimization. Finally, we will implement our methods using a GPU, highlighting performance advantages over the more conventional CPU.