PaperlayerView exact revision on arXiv

Source-linked AI summary

Hierarchical Nearest-Neighbor Gaussian Process Models for Large Geostatistical Datasets

Abhirup Datta, Sudipto Banerjee, Andrew O. Finley, Alan E. Gelfand

arXiv:1406.7343v2stat.ME

TL;DR

Large geostatistical datasets make conventional spatial-process inference computationally prohibitive. This paper develops valid, scalable NNGP models with sparse precision structure and embeds them in hierarchical Bayesian inference. The NNGP matches full-GP inferential performance in the reported setting and significantly outperforms competing low-rank processes while enabling analysis at massive scale.

  • Problem

    Conventional Gaussian-process inference requires computations and storage that become prohibitive for large datasets, while scalable alternatives do not provide a fully process-based framework for broad inference.

  • Method

    The paper constructs NNGPs from directed-acyclic-graph neighbor sets and uses them as sparsity-inducing priors in hierarchical Bayesian spatial models.

  • Results

    The NNGP performs at par with the full GP, significantly outperforms competing low-rank processes, and shows very robust inference across reference-set choices.

  • Takeaways & Limitations

    NNGPs provide a highly scalable model for fully process-based inference, including estimation, prediction, and interpolation at massive geostatistical scale.

  • Takeaways & Limitations

    The resulting NNGP may be a poor approximation of the full Gaussian process under some reference-set choices.

Abstract

from arXiv · show

Spatial process models for analyzing geostatistical data entail computations that become prohibitive as the number of spatial locations become large. This manuscript develops a class of highly scalable Nearest Neighbor Gaussian Process (NNGP) models to provide fully model-based inference for large geostatistical datasets. We establish that the NNGP is a well-defined spatial process providing legitimate finite-dimensional Gaussian densities with sparse precision matrices. We embed the NNGP as a sparsity-inducing prior within a rich hierarchical modeling framework and outline how computationally efficient Markov chain Monte Carlo (MCMC) algorithms can be executed without storing or decomposing large matrices. The floating point operations (flops) per iteration of this algorithm is linear in the number of spatial locations, thereby rendering substantial scalability. We illustrate the computational and inferential benefits of the NNGP over competing methods using simulation studies and also analyze forest biomass from a massive United States Forest Inventory dataset at a scale that precludes alternative dimension-reducing methods.

1 Introduction

Large geostatistical datasets make conventional Gaussian-process inference computationally prohibitive, while existing scalable approximations have important inferential or computational limitations. The paper introduces NNGP models as a valid, fully process-based framework for scalable estimation, prediction, and interpolation.

  • Computational challenge: ∼n3 flops and n2 storage for covariance operations become prohibitive when the number of locations is large and covariance structure is unavailable.Conventional model fitting requires the inverse and determinant of the covariance matrix.
  • Existing approximations: Low-rank models reduce fitting cost to O(nr2 + r3) ≈ O(nr2), but large datasets may require large r, making computation exorbitant.They construct the process on a lower-dimensional subspace using r << n knots or centers.
  • Existing approximations: Low-rank models can perform poorly when neighboring observations are strongly correlated and spatial signal dominates noise, while bias adjustment increases computational burden.These limitations motivate alternatives based on sparsity.
  • Existing approximations: Covariance tapering is effective for parameter estimation and kriging but has not been fully developed for general inference on residual or latent processes.Other sparse likelihood methods may also lack a corresponding process for arbitrary new locations, limiting process-based prediction and uncertainty assessment.
  • NNGP contribution: The NNGP constructs valid spatial processes from neighbor sets in directed acyclic graphs, with finite-dimensional Gaussian realizations and closed-form sparse precision matrices.It extends finite-dimensional sparse probability models to a valid process over uncountable spatial domains.
  • NNGP contribution: Embedded as a Bayesian sparsity-inducing prior, the NNGP supports joint inference for model parameters and spatial processes at observed and unobserved locations.The framework includes spatially varying regression and is demonstrated on forest biomass at a scale where competing dimension-reduction methods are unimplementable.

2 Nearest-Neighbor Gaussian Process

The NNGP replaces large Gaussian-process conditioning sets with small neighbor sets organized through an acyclic graph, yielding a valid spatial process with sparse precision. Its construction extends finite-dimensional Gaussian densities consistently across the spatial domain.

  • Construction: The NNGP approximates each conditional distribution using at most m neighbors, with m much smaller than the reference-set size k.Neighbor sets are chosen from previously ordered locations to ensure an acyclic graph.
  • Validity: If the neighbor graph is acyclic, the resulting approximation is a proper multivariate Gaussian density.The covariance differs from that of the parent Gaussian process.
  • Sparsity: The resulting Gaussian distribution has a sparse precision matrix with at most km(m + 1)q^2/2 non-zero entries.Sparsity follows from the bounded neighbor-set size.
  • Practical choices: Inference is generally more sensitive to the number of neighbors than to the particular ordering of spatial locations.The authors report that their simulations found inference extremely robust to ordering, while selecting an optimal neighbor subset is difficult.
  • Spatial-process extension: Finite-dimensional NNGP densities satisfy Kolmogorov consistency criteria and therefore define a valid spatial process over the domain.The process can be derived from any parent spatial process and fixed reference set.
  • Properties: The NNGP is a proper, nondegenerate, sparsity-inducing Gaussian-process prior that remains usable when the reference set is as large as the dataset.The reference set may be larger than the observed dataset without undermining the construction.

3 Bayesian estimation and implementation

The hierarchical NNGP model replaces a customary Gaussian-process prior for spatial effects and supports Bayesian estimation, prediction, and MCMC without large matrix operations. Its computational cost is linear in the number of locations for fixed neighbor size, while small neighbor sets can closely reproduce full-geostatistical inference.

  • Hierarchical model: The hierarchical model uses an NNGP prior for spatially varying effects within a regression model containing fixed predictors and measurement error.The framework accommodates multivariate responses and spatially varying coefficients.
  • Bayesian estimation: The Gibbs sampler updates spatial effects and covariance parameters using local conditional distributions and sparse quantities.The updates avoid storing or factorizing any n × n matrices.
  • Computational complexity: The total per-iteration cost is O((n + k)m^3), which is linear in the total number of locations for fixed m.A full Gaussian process requires O(n^3) flops for updating spatial effects in each iteration.
  • Empirical performance: With m ≈ 10, NNGP inference is reported to be almost indistinguishable from full geostatistical models.The comparison is based on simulation results and appendix experiments.
  • Memory requirements: NNGP storage uses small m × m neighbor matrices for each location rather than the full n × n distance matrix.This reduces storage demands for very large datasets.
  • Reference-set choice: Unlike low-rank models, increasing the reference-set size raises NNGP cost only linearly and does not impose the same knot-size constraint.This permits large or flexible reference sets, including dense grids.
  • Reference-set choice: Choosing S = T avoids sampling additional reference-set effects and produces inference almost indistinguishable from using a grid reference set.The authors describe this as a legitimate choice because the NNGP is valid for any fixed reference set.

4 Alternate NNGP models and algorithms

The paper develops alternative NNGP formulations for response modeling, marginalized likelihoods, block updates, and spatiotemporal or non-Gaussian data. These alternatives trade computational convenience against residual-surface inference and can have structure-dependent costs or convergence limitations.

  • Block updates: Sequential updates can sometimes converge slowly, motivating block-update alternatives for wS.The stated computational efficiency does not eliminate this convergence issue.
  • Block updates: Block updates of wS can efficiently produce Cholesky factors and facilitate Gibbs sampling.This alternative addresses dependence among sequentially updated elements of wS.
  • Response modeling: An NNGP prior directly on the response can avoid full conditionals for the latent spatial effects.This produces a Bayesian analogue of earlier neighbor-based likelihood approaches but prevents inference on the residual surface w(s).
  • Response modeling: Modeling the latent process w(s), rather than only the response, preserves access to residual spatial contours useful for identifying unexplained patterns.The response-only formulation may not recover w when nugget effects are present.
  • Marginalized models: Marginalizing w yields a likelihood with covariance Σy = Z ˜CSZ′ + Dn and substantially reduces the number of Gibbs-sampler variables.The nugget effect from the parent model is retained, but conjugacy for covariance parameters can be lost.
  • Marginalized models: The marginalized sampler’s cost depends on the sparse structure of ˜C_S^-1 and can substantially exceed the linear cost of the unmarginalized model.The appropriate fitting algorithm should therefore be selected according to the dataset’s sparsity structure.
  • Extensions: For spatiotemporal models, independent NNGPs can replace independent Gaussian-process priors at each time point while preserving computational tractability.The framework supports spatial interpolation at discrete time points and can be incorporated into dynamic models.
  • Extensions: NNGP spatial effects can also be embedded in generalized linear models for binary, count, and other non-Gaussian responses.The spatial effects enter through a link function and retain structured dependence.

5 Illustrations

Simulation and forest-biomass analyses assess NNGP models against full Gaussian Process and predictive-process alternatives. NNGP models achieve similar or better inference and prediction while offering substantial computational advantages, and spatial effects improve biomass modeling.

  • Simulation experiment: The simulation compares Full GP, NNGP models across m values, and a 64-knot Gaussian Predictive Process using held-out observations.Models were trained on 2,000 of 2,500 locations, with 500 observations withheld for prediction.
  • Simulation experiment: NNGP models provide clear computational advantages over Full GP and both inferential and computational advantages over the predictive-process model.The comparison includes computing times for 25,000-iteration chains and fit metrics such as DIC and GPD.
  • Simulation experiment: NNGP predictive performance becomes comparable to the Full GP at approximately m = 10, with negligible differences between reference-set configurations.Both NNGP configurations outperform the 64-knot Predictive Process when the number of knots is small.
  • Simulation experiment: All models achieve appropriate 95% credible interval coverage, while NNGP spatial surfaces closely approximate the true and Full GP surfaces.The 64-knot predictive-process model greatly smooths over small-scale spatial patterns.
  • Forest biomass data analysis: The forest-biomass analysis uses spatial models to represent dependence that NDVI alone does not adequately capture.Non-spatial, SVI, and SVC models have PMSE values of 0.52, 0.41, and 0.42, respectively.
  • Forest biomass data analysis: Allowing the NDVI coefficient to vary spatially improves fit over the SVI model and reveals stronger positive relationships in the Pacific Northwest and northern California.Near-zero coefficients in western New England indicate that NDVI is less effective there at distinguishing forest-structure differences.

6 Summary and conclusions

The NNGP is presented as a scalable spatial-process model whose inference remains robust with modest neighbor sizes and suitable reference sets, while retaining a broad modeling scope. Its performance is strong for several covariance settings, although numerical and design limitations remain.

  • Contributions: NNGP models provide scalable spatial inference and are framed as models rather than likelihood approximations.The framework is presented as applicable to large geostatistical datasets and as a unifying approach for scalable modeling.
  • Computational and inferential performance: Inference is very robust to the reference set S, and modest values of m, typically much less than 20, usually suffice.Larger reference sets may be needed for larger datasets without thwarting computation.
  • Limitations: Using observed locations as the reference set can create poor approximations across large gaps because nearby points may receive very different neighbor sets.Simulations report a very flat NNGP covariance field in such gaps, although the full GP also provides weak information there.
  • Covariance settings: NNGP estimation and kriging closely emulate true Matérn GP models, including cases with slowly decaying covariances.The paper relates this performance to Matérn screening conditions and prediction based on a few neighbors.
  • Covariance settings: For wave covariances, NNGP estimates can remain close to true parameters and kriged surfaces can resemble the true surface, while numerical stability varies by covariance.The full GP can be unstable for some wave covariances; NNGP also encounters instability for the cardinal sine covariance.
  • Future scope: The NNGP is extensible to multivariate, discretized spatiotemporal, network, and regionally aggregated spatial settings.The paper identifies additional scope for jointly modeling space and time through spatiotemporal covariance functions.

B Properties of ˜C

The NNGP covariance representation is built from sparse conditional-regression blocks, yielding a sparse precision structure whose nonzero pattern is controlled by the neighbor size.

  • Sparse construction: Conditional Gaussian coefficients define sparse block rows with at most m + 1 nonzero blocks.The construction uses neighbor sets of size m for each location.
  • Sparse construction: The matrix B_S is sparse and lower triangular with ones on its diagonal.This structure follows directly from the block construction based on ordered neighbor sets.
  • Sparsity pattern: Nonzero covariance-precision blocks arise only when locations share a relevant neighbor-set relationship.Each neighbor set contributes at most m(m + 1)/2 potentially nonzero block pairs.
  • Sparsity pattern: For m much smaller than k, the resulting precision matrix remains sparse relative to the number of locations.The sparsity argument is based on the bounded size of each neighbor set.

ordering of locations

A simulation on a long, skinny domain examines whether coordinate ordering affects NNGP inference. Across three orderings, parameter and spatial-effect estimates remain highly consistent with full-GP inference.

  • Simulation design: The experiment uses n = 2500 locations in a long, skinny domain to expose sensitivity from unequal x- and y-axis scales.The locations are ordered by x, by y, or by f(x, y) = x + y.
  • Parameter inference: Point estimates and 95% credible intervals from all three NNGP models are extremely consistent with the full Gaussian process model.The comparison is summarized in Table 3.
  • Spatial effects: The impact of alternative ordering on posterior estimates of the spatial residual surface is negligible.Figure 5 displays surfaces for the full model and NNGP alternatives with S = T and m = 10.
  • Spatial effects: Differences between full-GP and NNGP spatial-effect estimates are negligible relative to differences between true effects and full-GP estimates.This result holds for orderings by x, y, and x + y coordinates.

D Kolmogorov Consistency for NNGP

The NNGP construction is shown to define proper finite-dimensional densities and to be invariant to permutations of locations under a fixed reference set.

  • Proper densities: The NNGP conditional specification defines a proper density over any finite set of locations.The proof integrates conditional factors over an ordering induced by an acyclic directed graph.
  • Permutation and consistency: With S fixed, the NNGP density for a finite set V depends only on neighbor sets for locations outside S and is invariant to location permutations.The consistency argument extends the density by integrating out locations not in the reference set.

E Properties of NNGP

The NNGP defines valid Gaussian finite-dimensional distributions while inducing sparse covariance and precision structures. Its inference closely resembles full-GP inference, including for small neighbor sets.

  • Sparsity: BU is sparse because each row has at most m non-zero entries.This sparsity follows from conditioning each location on at most m neighbors.
  • Validity: All finite-dimensional realizations of the NNGP process are Gaussian because its component densities are Gaussian.The construction combines Gaussian densities for the reference set and conditional Gaussian distributions outside it.
  • Covariance: For locations outside S, the NNGP covariance combines a conditional variance term with neighbor-mediated covariance.The covariance is expressed using Fv1, Bv1, and the covariance among neighbor locations.
  • Approximation: Increasing m improves the approximation, with the NNGP becoming identical to the full GP when m equals the sample size.The paper treats NNGP as an independent model while noting this limiting relationship.
  • Inference: NNGP inference closely emulates full-GP inference, and parameter credible intervals remain nearly identical even for small m.The simulation compares full GP with NNGP settings including m = 10 and m = 100.

G Simulation experiment: Data with gaps

When data contain large gaps, an NNGP using data locations as its reference set can poorly approximate the full GP as a spatial process in gap regions. A grid reference set improves covariance approximation, while kriging and model-fitting results remain similar to the full GP.

  • Gap effects: Large gaps can make the NNGP covariance nearly uncorrelated for nearby locations inside the gap.The example evaluates a point surrounded by data and another point centered in the gap.
  • Gap effects: A reference set with large gaps can yield a poor full-GP process approximation in some domain regions.The issue arises because locations outside S are correlated through their neighbor sets, which may be far away.
  • Reference-set choice: A uniformly distributed grid reference set produces NNGP covariance functions that closely resemble the true GP covariance.The grid used 14 × 7 points over [0, 3]×[0, 1], with a size similar to the original sample.
  • Kriging: Full-GP and NNGP kriging means and variances are very close even for data with gaps.The comparison uses domain-wide kriging surfaces for the two models.
  • Model comparison: Full-GP and NNGP models produce very similar parameter estimates and kriging results for locations with gaps.The analysis compares posterior surfaces, variance surfaces, and model-fitting results.
  • Scope: Neither full GP nor NNGP with S = T provides enough information for locations inside large gaps.Thus, poor process approximation by NNGP in gaps does not necessarily translate into different model-fitting performance.

H Simulation experiment: Slow decaying covariance

The slow-decaying covariance experiment tests whether nearest-neighbor conditioning can recover inference when distant observations remain correlated. Across simulated parameter settings, NNGP posterior inference closely matches full-GP inference.

  • Limitation: A potential limitation is that m-nearest neighbors may miss information from distant observations under very flat-tailed covariance functions.Such functions can keep distant observations significantly correlated with a given location, affecting covariance-parameter information.
  • Simulation design: NNGP is evaluated for marginal variances σ2 from 0.05 to 0.5 and effective ranges from 0.1 to 1, with nugget variance τ2 = 0.1.The simulations vary σ2 and 3/φ while holding τ2 fixed.
  • Results: NNGP and full-GP posterior samples look identical across all simulated parameter choices.Figure 13 compares confidence intervals for the variance and effective range parameters.
  • Results: NNGP inference remains similar to full-GP inference for slowly decaying covariance functions.The result supports the chosen neighbor sets across the tested parameter range.

I Simulation experiment: Wave covariance function

The wave-covariance experiment examines NNGP performance when correlation does not decrease monotonically with distance. NNGP closely approximates the wave GP in parameter estimation and kriging, while avoiding unstable full-GP computation.

  • Motivation: The experiment tests NNGP with a two-dimensional damped cosine covariance function whose correlations do not monotonically decrease with distance.This setting probes performance beyond the commonly used Matérn covariance functions.
  • Neighbor selection: Nearest-neighbor selection yields lower KL divergence than the alternate selection scheme.The alternate scheme's KL numbers are always higher in the comparison.
  • Approximation: For damped cosine covariance, KL divergence is small when m ≥25 across sample sizes and neighbor-selection schemes.Small KL values indicate close approximation to the true damped cosine GP.
  • Computation: The full GP could not be fitted because of computational instability in the large wave covariance matrix.NNGP fitting remained feasible because it avoids inverting large matrices.
Loading 1406.7343v2…