Iterative Methods for Vecchia-Laplace Approximations for Latent Gaussian
  Process Models
Pascal Kündig Fabio Sigrist
Latent Gaussian process (GP) models are flexible probabilistic non-parametric
function models. Vecchia approximations are accurate approximations for GPs to
overcome computational bottlenecks for large data, and the Laplace
approximation is a fast method with asymptotic convergence guarantees to
approximate marginal likelihoods and posterior predictive distributions for
non-Gaussian likelihoods. Unfortunately, the computational complexity of
combined Vecchia-Laplace approximations grows faster than linearly in the
sample size when used in combination with direct solver methods such as the
Cholesky decomposition. Computations with Vecchia-Laplace approximations thus
become prohibitively slow precisely when the approximations are usually the
most accurate, i.e., on large data sets. In this article, we present several
iterative methods for inference with Vecchia-Laplace approximations which make
computations considerably faster compared to Cholesky-based calculations. We
analyze our proposed methods theoretically and in experiments with simulated
and real-world data. In particular, we obtain a speed-up of an order of
magnitude compared to Cholesky-based inference and a threefold increase in
prediction accuracy in terms of the continuous ranked probability score
compared to a state-of-the-art method on a large satellite data set. All
methods are implemented in a free C++ software library with high-level Python
and R packages.