PaperlayerView exact revision on arXiv

Source-linked AI summary

A Case Study Competition Among Methods for Analyzing Large Spatial Data

Matthew J. Heaton, Abhirup Datta, Andrew Finley, Reinhard Furrer, Rajarshi Guhaniyogi, Florian Gerber, Robert B. Gramacy, Dorit Hammerling, Matthias Katzfuss, Finn Lindgren, Douglas W. Nychka, Furong Sun, Andrew Zammit-Mangion

arXiv:1710.05013v2stat.ME

TL;DR

Large spatial datasets make traditional Gaussian-process analysis computationally infeasible. This paper reviews scalable alternatives and compares expert implementations in a blinded competition on simulated and real data, finding strong predictive accuracy across methods while noting important scope and comparison limits.

  • Problem

    Traditional Gaussian processes require O(N^3) operations and O(N^2) memory, making them computationally infeasible for large spatial datasets.

  • Method

    The paper reviews scalable spatial methods and compares implementations produced by expert research groups using blinded simulated and real datasets on a common computing environment.

  • Results

    All methods performed very well in predictive accuracy on both simulated and real data, while uncertainty-quantification performance varied across methods.

  • Takeaways & Limitations

    The competition provides a review, practitioner code, and a common-task framework for future comparisons of large-spatial-data methods.

  • Takeaways & Limitations

    The study is non-comprehensive because some existing and future large-spatial-data methods were not included.

Abstract

from arXiv · show

The Gaussian process is an indispensable tool for spatial data analysts. The onset of the "big data" era, however, has lead to the traditional Gaussian process being computationally infeasible for modern spatial data. As such, various alternatives to the full Gaussian process that are more amenable to handling big spatial data have been proposed. These modern methods often exploit low rank structures and/or multi-core and multi-threaded computing environments to facilitate computation. This study provides, first, an introductory overview of several methods for analyzing large spatial data. Second, this study describes the results of a predictive competition among the described methods as implemented by different groups with strong expertise in the methodology. Specifically, each research group was provided with two training datasets (one simulated and one observed) along with a set of prediction locations. Each group then wrote their own implementation of their method to produce predictions at the given location and each which was subsequently run on a common computing environment. The methods were then compared in terms of various predictive diagnostics. Supplementary materials regarding implementation details of the methods and code are available for this article online.

1. Introduction

Gaussian processes are powerful but become computationally infeasible for large spatial datasets, motivating scalable approximations and parallel methods. The paper reviews these methods and compares expert implementations through a blinded competition on simulated and real data.

  • Gaussian-process density evaluation requires O(N^3) operations and O(N^2) memory, limiting its use for moderately large datasets.
  • Existing alternatives include conditional approximations, pseudo-likelihoods, low-rank methods, lattice kriging, and stochastic partial differential equations.
  • Many approximations can oversmooth data, and some impose upper limits on dataset size, motivating parallelizable methods for modern computing platforms.
  • The study gives an overview of modern methods and conducts a blinded competition using one simulated and one real spatial dataset.
  • Each method was implemented by a research group with strong expertise, reducing concerns that differing familiarity would make comparisons unfair.

2. Overview of Methods for Analyzing Large Spatial Data

Fixed Rank Kriging and LatticeKrig use multiresolution basis representations to reduce large spatial computations. FRK models fine-scale residual variation, whereas LatticeKrig uses sparse grid-based dependence structures.

  • Fixed Rank Kriging: FRK models a spatial process as a mean function, spatial random effects, and a spatially uncorrelated fine-scale process.
  • Fixed Rank Kriging: FRK represents spatial random effects with K basis functions, reducing estimation and prediction inversions to K × K matrices when K ≪ N.
  • Fixed Rank Kriging: The FRK implementation uses three resolutions of bisquare basis functions and block-diagonal coefficient covariance.
  • LatticeKrig: LatticeKrig uses the FRK setup without the fine-scale process, but applies a different multiresolution structure and constraints.
  • LatticeKrig: LatticeKrig uses normalized Wendland polynomial basis functions across resolutions to reduce edge effects and better approximate stationary covariance.
  • LatticeKrig: LatticeKrig places basis functions on regular grids and uses sparse spatial autoregression or Markov random fields, allowing K to exceed N with limited additional cost.

2.3. Predictive Processes

Predictive processes approximate spatial random effects using values at dispersed knot locations and covariance-based basis functions. This yields substantial computational savings but can underestimate variance and require repeated basis recalculation.

  • Predictive processes replace spatial random effects with a covariance-based representation defined by well-dispersed knot locations.
  • The predictive-process basis functions are determined by the selected covariance function rather than separately specified.
  • The basis expansion systematically underestimates the marginal variance, so an additional structure is used to counterbalance that deficiency.
  • The resulting likelihood requires inversion and determinant calculations for a dense K × K matrix plus diagonal N × N matrices when K ≪ N.

2.4. Spatial Partitioning

Spatial partitioning divides the domain into subregions with conditional within-region dependence and cross-region independence, enabling parallel likelihood computation while using all observations. Covariance tapering instead induces sparsity through compact support, trading computational efficiency against parameter-estimation bias.

  • Spatial Partitioning: Spatial partitioning assumes conditional dependence within subregions and conditional independence across subregions.
  • Spatial Partitioning: Shared regression and basis parameters allow smoothing across regions while cross-region independence enables parallel likelihood computation.
  • Spatial Partitioning: Unlike divide and conquer, spatial partitioning uses all observations simultaneously rather than fitting separate subset models and pooling results.
  • Spatial Partitioning: Partition choices include equal-area, centroid-clustering, and hierarchical-clustering approaches based on spatial gradients.
  • Covariance Tapering: Covariance tapering multiplies the covariance function by a compactly supported covariance, producing sparse matrices for efficient linear-system solves.
  • Covariance Tapering: One-taper likelihoods produce biased parameter estimates, whereas two-taper estimating equations are unbiased but substantially less computationally efficient.

2.6. Multiresolution Approximations

The MRA represents spatial processes with compactly supported basis functions across recursive resolutions, adapting to covariance smoothness while enabling scalable computation. The section also contrasts this strategy with nearest-neighbor Gaussian-process constructions.

  • Multiresolution approximation: MRA combines multiresolution basis expansions with predictive-process choices that adapt basis functions and weights to the covariance’s spatial smoothness.Its basis supports become smaller at finer resolutions, increasing covariance sparsity.
  • MRA-block: MRA-block recursively partitions the domain into smaller subregions, using a small number of basis functions per region at each resolution.At the finest resolution, the approximation of the process, variance, and smoothness within each region is exact.
  • MRA-block: MRA-block inference assigns computational nodes to regions, allowing parallel processing with limited communication overhead and polylogarithmic computational complexity.The regional decomposition supports massively distributed computing.
  • Nearest neighbor Gaussian process: NNGP approximates the full Gaussian-process conditional specification by retaining only a fixed number of nearest neighbors, producing a sparse formulation for scalable likelihood evaluation.The resulting sparse matrix has at most m nonzero entries per row.
  • Nearest neighbor Gaussian process: The conjugate NNGP model fixes selected covariance parameters by grid-search cross-validation, offering a pragmatic alternative to fully Bayesian inference for massive spatial datasets.Predictive distributions at new locations can be obtained as t-distributions.

2.8. Stochastic PDEs

The SPDE approach uses the relationship between Matérn covariance fields and stochastic PDEs to construct sparse spatial models. Its implementation combines triangulated basis functions, sparse precision calculations, and periodic-embedding ideas for gridded data.

  • Stochastic PDE formulation: SPDE exploits the equivalence between Matérn covariance fields and stochastic PDEs, together with a Markov property for suitable smoothness parameters on two-dimensional domains.The spatial field is represented using basis functions that are piecewise linear on a domain triangulation.
  • Model formulation: The SPDE observation model combines covariates, sparse basis evaluations, latent coefficients, and diagonal observation noise in a joint Gaussian formulation.The design matrix is A = (H, X), where H evaluates basis functions and X contains covariates.
  • Computation: Posterior inference uses sparse precision equations and Cholesky factorization of Qz|y, with automatic reordering to preserve sparsity.The resulting computational and storage cost for posterior predictions and likelihood evaluation is O(K3/2).
  • Computation: Takahashi recursions recover the posterior covariance elements needed to compute predictive variances at new locations.Only covariance elements corresponding to the nonzero structure of Qz|y are obtained.
  • Periodic embedding: For regular stationary grids, spectral methods use approximately decorrelating DFTs and FFTs, while periodic embedding and imputation address edge effects and missing values.The method expands the domain and imputes missing observations so the expanded vector is periodic.

2.10. Metakriging

Spatial metakriging combines inference from data subsets into a metaposterior rather than relying on one full-data computation. It uses geometric-median aggregation and subset posterior predictive distributions for scalable prediction.

  • Subset inference: Spatial metakriging applies a spatial model independently to K data subsets in parallel, producing subset posteriors for the model parameters.The approach is algorithmic and is not tied to one specific spatial model.
  • Metaposterior construction: The method combines subset posteriors into a legitimate metaposterior by using their unique geometric median in a Banach space.The geometric median lies in the convex hull of the individual posteriors, ensuring a valid probability density.
  • Metaposterior construction: Weiszfeld’s iterative algorithm estimates the geometric-median weights because the defining equation has no apparent closed-form solution.The algorithm estimates αρ,k(y) for each subset posterior.
  • Prediction: Posterior prediction is approximated by combining samples from the subset posterior predictive distributions, yielding predictive medians and 95% posterior predictive intervals.The empirical metaposterior predictive distribution is formed from subset-level predictive samples.
  • Implementation: The study uses K = 30 exhaustive, mutually exclusive random subsets so each subset can reasonably represent the full spatial domain.The chosen partitioning scheme supports the article’s metakriging implementation.

2.11. Gapfill

Gapfill is a distribution-free, non-Gaussian-process method designed for low computational workload and parallel processing. It predicts each missing value from a sufficiently large local neighborhood of observed pixels.

  • Gapfill method: Gapfill is purely algorithmic, distribution-free, and not based on Gaussian processes, distinguishing it from the other methods presented.Its low computational workload makes it attractive for satellite imaging applications.
  • Prediction procedure: For a target location s0, gapfill selects a local neighborhood containing enough observed values, then bases prediction on that subset using sorting algorithms.This study requires a neighborhood spanning at least 5 × 5 pixels and containing at least 25 non-missing values.
  • Implementation: The gapfill implementation was designed for spatio-temporal data, so this study supplies shifted versions of the images to mimic multiple-image input.The algorithm is applied to 13 images consisting of the original and shifted images.
  • Local approximate Gaussian processes: Local approximate Gaussian processes use prediction-specific local designs selected through a greedy search, with calculations for distinct prediction locations parallelizable.A local design size of m = 50 typically works well, while global predictive covariance is unavailable.
  • Competition context: The competition compared participating research groups and their selected methods using common datasets and a shared evaluation setting.Table 1 identifies the participating groups and associated competitors.

3. The Competition

The competition compared expert implementations of large-spatial-data methods on simulated and observed datasets using common prediction, uncertainty, and computation outputs.

  • Participating teams implemented their associated methods after outreach to groups with strong methodological expertise.
  • The competition used one real and one simulated dataset on the same 500×300 spatial grid.The simulated data used a correctly specified covariance function, whereas the real data represented covariance misspecification.
  • The datasets were split into training and test sets using a realistic missing-data pattern from MODIS observations.Both training sets contained 105,569 observations; the MODIS and simulated test sets contained 42,740 and 44,431 observations, respectively.
  • Each group supplied test-location predictions, 95% prediction intervals or standard errors, per-iteration time, and total clock time.Teams were also instructed to use an exponential correlation function when applicable to reduce confounding factors.
  • Methods were compared using predictive accuracy, uncertainty quantification, and computation time diagnostics.The reported metrics included MAE, RMSE, CRPS, interval score, and prediction-interval coverage.

4. Competition Results

Across simulated and satellite data, competing methods generally achieved strong predictive accuracy, while uncertainty quantification showed the clearest differences and depended on the data and covariance setting.

  • 4.1. Results for Simulated Data: Simulated-data MAE ranged from 0.61 to 1.03, while RMSE ranged from 0.83 to 1.31 across methods.The paper describes even the worst RMSE as accurate relative to the simulated-data range of 19.89.
  • 4.1. Results for Simulated Data: On simulated data, LatticeKrig, LAGP, metakriging, MRA, periodic embedding, and NNGP achieved near-nominal 95% coverage.FRK, Gapfill, partitioning, and PP had lower coverage, whereas SPDE and tapering had higher coverage.
  • 4.1. Results for Simulated Data: FRK and Gapfill were less suited to exponential covariance fields, which may partially explain poorer prediction or coverage in the simulated setting.The comparison also notes that several other methods could specify or approximate the correct exponential correlation function.
  • 4.1. Results for Simulated Data: Simulated-data predictive surfaces showed a smeared Gapfill pattern and possible oversmoothing for FRK and PP.The Gapfill pattern followed image shifting used to create pseudo-datasets, while the oversmoothing interpretation is linked to prior arguments about low-rank methods.
  • 4.2. Results for Real Data: The largest satellite-data RMSE was 2.52 over a data range of 31.04, and predictive accuracy was described as very good.Under the competition setup, low scores reflected good data fit rather than necessarily the best approximation of the exponential covariance.
  • 4.2. Results for Real Data: On MODIS data, LatticeKrig, metakriging, MRA, NNGP, and periodic embedding again had near-nominal coverage with small interval scores and CRPS.SPDE and tapering also approached nominal coverage, while FRK, Gapfill, LAGP, partitioning, and predictive processes had too-small coverage and larger interval scores.
  • 4.2. Results for Real Data: MODIS predictive surfaces were generally smoother in the north-central region, where predictions were long-range and lacked nearby observations.The Gapfill shifting procedure again produced a smeared surface, which the authors hypothesize would not occur in its spatio-temporal setting.

5. Conclusions

The study finds that the compared methods generally perform well for prediction but differ in uncertainty quantification. Its conclusions are bounded by spatial-only, non-comprehensive comparisons across limited data scenarios and non-standardized implementations.

  • Comparison results: The methods performed very well in predictive accuracy, suggesting that each is suited to spatial prediction.The authors also report differences in uncertainty quantification and note that performance can change with dataset size, measurement error variance, and missingness.
  • Computational limitations: Runtime comparisons can be misleading because implementations differ across languages, libraries, optimization levels, and development maturity.Reported runtimes are characterized as off-the-shelf estimates rather than optimized comparisons.
  • Computational limitations: Results may partly reflect unequal effort among participating groups, including differences in coding time and prior cross-validation.The authors nevertheless regard the findings as informative about methods’ strengths and weaknesses.
  • Scope: The competition was non-comprehensive because other existing methods and future methods were not included.The authors consider the compared methods representative of common large-spatial-data methods at the time of writing.
  • Scope: The tested scenarios do not cover anisotropy, non-stationarity, varying spatial ranges, or varying signal-to-noise ratios.Further competitions could distinguish methods under these more complicated spatial settings.
  • Scope: The study’s results apply only to spatial data, while spatio-temporal data require further comparative research.The authors note that spatio-temporal data are often larger and more complex, although several methods have extensions to that setting.

Supplement A: Computational Details

The supplementary materials provide implementation details, datasets, and code for reproducing the reported analyses.

  • Supplement B provides implementation details for some methods applied to the training datasets.
  • Supplement B provides method code, datasets, and code to reproduce results such as Tables 2 and 3.
Loading 1710.05013v2…