Vecchia–Laplace approximations of generalized Gaussian processes for big non-Gaussian spatial data

Abstract

Generalized Gaussian processes (GGPs) are highly flexible models that combine latent GPs with potentially non-Gaussian likelihoods from the exponential family. GGPs can be used in a variety of settings, including GP classification, nonparametric count regression, modeling non-Gaussian spatial data, and analyzing point patterns. However, inference for GGPs can be analytically intractable, and large datasets pose computational challenges due to the inversion of the GP covariance matrix. A Vecchia–Laplace approximation for GGPs is proposed, which combines a Laplace approximation to the non-Gaussian likelihood with a computationally efficient Vecchia approximation to the GP, resulting in simple, general, scalable, and accurate methodology. Numerical studies and comparisons on simulated and real spatial data are provided. The methods are implemented in a freely available R package.

Introduction

Dependent non-Gaussian data are ubiquitous in time series, geospatial applications, and more generally in nonparametric regression and classification. A popular model for such data is obtained by combining a latent Gaussian process (GP) with conditionally independent, potentially non-Gaussian likelihoods from the exponential family. This is traditionally referred to as a spatial generalized linear mixed model (SGLMM) in the spatial statistics literature (Diggle et al., 1998), but the same model has more recently also been referred to as a generalized GP (GGP) (e.g., Chan and Dong, 2011); we will use the latter, more concise term throughout. GGPs are highly flexible, interpretable, and allow for natural, probabilistic uncertainty quantification. However, inference for GGPs can be analytically intractable, and large datasets pose additional computational challenges due to the inversion of the GP covariance matrix.

Popular techniques to numerically perform the intractable marginalization necessary for inference are, in order of increasing speed: Markov chain Monte Carlo (MCMC), expectation propagation, variational methods, and Laplace approximations. See Shang and Chan (2013) for a recent review of deterministic techniques and Filippone and Girolami (2014) for a comparison of MCMC and expectation propagation. Rue et al. (2009) argue that variational methods and expectation propagation suffer from underestimated and overestimated posterior variances, respectively. Here, we consider the Laplace approximation (e.g., Tierney and Kadane, 1986, Williams and Barber, 1998), a classic technique for integral evaluation based on second-order Taylor expansion. Bonat and Ribeiro (2016) show numerically that the Laplace approximation can be a practical and accurate method for GGP inference.

It has long been recognized that the computational cost for GPs is cubic in the data size. Much work has been doneon GP approximations that address this problem in the context of Gaussian noise (as recently reviewed in Heaton et al., 2019). Low-rank approaches (e.g., Higdon, 1998, Quiñonero-Candela and Rasmussen, 2005, Banerjee et al., 2008, Cressie and Johannesson, 2008, Katzfuss and Cressie, 2011) have limitations in the presence of fine-scale structure (Stein, 2014), but they have proved popular due to their simplicity. Approximations relying on sparsity in covariance matrices (Furrer et al., 2006, Kaufman et al., 2008) by definition can only capture local, short-range dependence and cannot guarantee low computation cost. Approaches that take advantage of Toeplitz or Kronecker structure (e.g., Dietrich and Newsam, 1997, Flaxman et al., 2015, Guinness and Fuentes, 2017) can be extremely efficient but are not as generally applicable. Recently, methods relying on sparsity in precision matrices (Rue and Held, 2005, Lindgren et al., 2011, Nychka et al., 2015) have gained popularity due to their accuracy and flexibility. In particular, a class of highly promising GP approximations (Vecchia, 1988, Stein et al., 2004, Datta et al., 2016, Guinness, 2018, Katzfuss and Guinness, 2019, Katzfuss et al., 2020a) rely on ordered conditional independence that can guarantee linear scaling in the data size while resolving dependence at all scales.

There are also a number of existing methods for large non-Gaussian datasets modeled using GGPs. A popular approach is to combine a low-rank GP with an approximation of the non-Gaussian likelihood, as the dimension reduction inherent in the low-rank approximation carries through even to the non-Gaussian case. Sengupta and Cressie (2013) estimate parameters using an expectation–maximization algorithm with low-rank and Laplace approximations. Sheth et al. (2015) use variational inference to obtain the posterior and select a set of conditioning points for their low-rank approximation. Some methods of dimension reduction, including random sketching (e.g., Yang et al., 2017) and projection, offer theoretical guarantees and can be combined with MCMC methods for the analysis of non-Gaussian data (e.g., Hughes and Haran, 2013, Guan and Haran, 2018), but are still subject to the limitations of low-rank methods. Nickisch et al. (2018) develop state–space models for one-dimensional non-Gaussian time series, which can perform inference in linear time and memory using a set of knots in time, a form of low-rank approximation. Alternate priors (e.g., log gamma priors for count data, Bradley et al., 2018) are an interesting but specialized approach to completely avoid the intractability issues with GGPs.

Similar to what we shall propose, some authors have combined a sparse-precision approach with a non-Gaussian approximation. A prominent example is Lindgren et al. (2011), in which an integrated nested Laplace approximation (INLA) is combined with a sparse-precision approximation of the GP using its representation as the solution to a stochastic partial differential equation. Datta et al. (2016) proposed to apply the GP approximation of Vecchia (1988) to a latent GP, but did not provide an explicit algorithm for large non-Gaussian data. While both Lindgren et al. (2011) and Datta et al. (2016) limit the number of nonzero entries per row or column in the precision matrix to a small constant, the computational complexity for decomposing this sparse $n\times n$ matrix is not linear in $n$, but rather $\mathcal{O}\left(n^{3∕2}\right)$ in two dimensions (Lipton et al., 1979 Thm. 6), and at least $\mathcal{O}\left(n^2\right)$ in higher dimensions. In the Gaussian setting, this scaling problem can be overcome by applying a Vecchia approximation to the observed data (Vecchia, 1988) or to the joint distribution of the observed data and the latent GP (Katzfuss and Guinness, 2019).

To handle both scaling and intractability, we propose a Vecchia–Laplace (VL) approximation for GGPs. The posterior mode necessary for the Laplace approximation is found using the Newton–Raphson algorithm, which can be viewed as iterative GP inference based on Gaussian pseudo-data. At each iteration of our VL algorithm, the joint Gaussian distribution of the pseudo-data and the latent GP realizations is approximated using the general Vecchia framework (Katzfuss and Guinness, 2019, Katzfuss et al., 2020a). By modeling the joint distribution of pseudo-data and GP realizations at each iteration, our VL approach can ensure sparsity and guarantee linear scaling in $n$ for any dimension, overcoming the scaling issues of the sparse-matrix approaches mentioned above.

To our knowledge, we provide the first explicit algorithm extending and applying the highly promising class of general-Vecchia GP approximations to large non-Gaussian data. We believe it to be a useful addition to the literature due to its speed, simplicity, guaranteed numerical performance, and wide applicability (e.g., binary, count, right-skewed, and point-pattern data). In addition, as shown in Katzfuss and Guinness (2019), the general Vecchia approximation includes many popular GP approximations (e.g., Vecchia, 1988, Snelson and Ghahramani, 2007, Finley et al., 2009, Sang et al., 2011, Datta et al., 2016, Katzfuss, 2017, Katzfuss and Gong, 2019) as special cases, and so our VL methodology also directly provides extensions of these GP approximations to non-Gaussian data.

The remainder of this document is organized as follows. In Section 2, we review the Laplace approximation and general Vecchia. In Section 3, we introduce and examine our VL method, including parameter inference and predictions at unobserved locations. In Sections 4 Simulations and comparisons, 5 Application to satellite data, we study and compare the performance of VL on simulated and real data, respectively. Some details are left to the appendix. A separate Supplementary Material document contains Sections S1–S6 with additional derivations, simulations, and discussion. The methods and algorithms proposed here are implemented in the R package GPvecchia (Katzfuss et al., 2020) with sensible default settings, so that a wide audience of practitioners can immediately use the code with little background knowledge. Our results and figures can be reproduced using the code and data at https://github.com/katzfuss-group/GPvecchia-Laplace.