PaperlayerView exact revision on arXiv

Source-linked AI summary

A multi-resolution approximation for massive spatial datasets

Matthias Katzfuss

arXiv:1507.04789v2stat.MEstat.CO

TL;DR

Massive spatial datasets make exact Gaussian-process inference computationally infeasible because dense covariance operations scale poorly and parallelize difficultly. The paper proposes the multi-resolution approximation (M-RA), whose optimized multi-scale basis representation enables scalable inference and prediction; compared with the full-scale approximation, it provides better approximation at equal computational cost or similar approximation faster.

  • Problem

    Exact spatial-statistical inference for massive datasets is hindered by O(n3) time, O(n2) memory, and substantial communication overhead from dense covariance computations.

  • Method

    The M-RA approximates Gaussian processes with basis functions at multiple resolutions, optimally chosen for a given covariance function, and uses a block-sparse precision structure for scalable inference.

  • Results

    The M-RA compares favorably with the full-scale approximation, providing a better approximation at the same computational complexity and time as the 1-RA or a similar approximation at a fraction of the computation time.

  • Takeaways & Limitations

    The M-RA supports scalable and distributed computations for Gaussian-process inference and prediction on massive spatial datasets.

  • Takeaways & Limitations

    The model assumes a zero-mean process and known covariance structure up to parameters θ; in applications, the mean must generally be estimated and subtracted.

Abstract

from arXiv · show

Automated sensing instruments on satellites and aircraft have enabled the collection of massive amounts of high-resolution observations of spatial fields over large spatial regions. If these datasets can be efficiently exploited, they can provide new insights on a wide variety of issues. However, traditional spatial-statistical techniques such as kriging are not computationally feasible for big datasets. We propose a multi-resolution approximation (M-RA) of Gaussian processes observed at irregular locations in space. The M-RA process is specified as a linear combination of basis functions at multiple levels of spatial resolution, which can capture spatial structure from very fine to very large scales. The basis functions are automatically chosen to approximate a given covariance function, which can be nonstationary. All computations involving the M-RA, including parameter inference and prediction, are highly scalable for massive datasets. Crucially, the inference algorithms can also be parallelized to take full advantage of large distributed-memory computing environments. In comparisons using simulated data and a large satellite dataset, the M-RA outperforms a related state-of-the-art method.

1 Introduction

The paper addresses the computational difficulty of applying spatial statistics to massive, irregularly observed spatial datasets by proposing a scalable multi-resolution Gaussian-process approximation.

  • Massive satellite and aircraft observations could support insights into climate, agriculture, and weather if efficiently exploited.
  • Traditional kriging is infeasible for big datasets because it requires decomposing dense n×n matrices.
  • Existing approximations often impose restrictive covariance assumptions or omit fine-scale or large-scale dependence.
  • The M-RA represents Gaussian processes at irregular locations with basis functions across multiple spatial resolutions.
  • The paper develops the M-RA, scalable inference and prediction algorithms, computational-complexity analysis, and comparisons on simulated and real data.
  • Its multi-resolution block-sparse precision structure and small regional basis sets support repeated Sherman-Morrison-Woodbury updates and scalable inference.

2 Multi-resolution approximation (M-RA)

This section introduces the true Gaussian process, establishes preliminaries, and presents the multi-resolution approximation and its properties.

  • The paper first describes the true Gaussian process that the M-RA will approximate.
  • It then provides preliminaries before introducing the multi-resolution approximation.
  • The section also discusses properties of the M-RA in subsequent subsections.

2.1 The true Gaussian process

The paper models a continuous spatial field as a zero-mean Gaussian process with a parameterized covariance, then explains why exact inference becomes infeasible for massive datasets.

  • The spatial field is modeled as a zero-mean Gaussian process on a continuous, non-gridded domain with covariance function C0 parameterized by θ.
  • Observed values at n spatial locations follow an n-variate Gaussian distribution with covariance matrix C0(S, S).
  • Spatial-statistical inference estimates θ and predicts the field at unobserved locations while quantifying the associated posterior distributions.
  • Exact computations require multiple Cholesky decompositions with O(n3) time and O(n2) memory complexity, becoming infeasible around n = 10^5.
  • Dense covariance computations are difficult to parallelize because they require substantial communication overhead, motivating approximations or simplifying assumptions.

2.2 Domain partitioning and knots

The M-RA organizes the spatial domain recursively into finer regions and assigns multi-resolution knot sets to support block-independent approximations and basis construction.

  • The spatial domain is recursively partitioned into J regions and then smaller subregions through level M.
  • At resolution m, the block-independent process preserves covariance within each region and sets covariance to zero between regions.
  • Each resolution-M subregion contains r knots, with r << n, and the notation Q(m) collects all knots at resolution m.
  • For convenience, the knots at the finest resolution are taken to be the observation locations in each region.
  • Equal numbers of subregions and knots are not required; the partitioning and knots are treated as fixed and known for inference.

2.3 Definition of the multi-resolution approximation (M-RA)

The M-RA recursively approximates a Gaussian process by predictive-process remainders across spatial resolutions, yielding a multiresolution basis-function representation plus a finest-resolution remainder.

  • At each resolution, the M-RA approximates the remainder between the true process and lower-resolution terms independently within spatial regions.The remainder is represented through a predictive process before proceeding to finer resolutions.
  • Predictive-process components are conditional-expectation approximations of the successive remainder processes.The construction begins with a predictive-process approximation of the original process and repeats it for each remainder.
  • Figure 1 contrasts a 1-RA with a 3-RA using equal computational complexity, showing increasingly fine covariance deviations across resolutions.The toy example has n = 54 observations, Mr = 6 for both models, and uses an exponential covariance on D = [0, 1].
  • The resulting process combines basis functions from resolutions 0 through M−1 with a remainder process at resolution M.The basis-function form uses region-specific weight vectors and a finest-resolution residual.

2.4 Properties of the M-RA

The M-RA is a valid Gaussian-process approximation whose multiresolution structure can represent fine- and large-scale dependence while improving approximation quality over a same-complexity full-scale approximation.

  • The total number of M-RA basis functions can exceed n, allowing representation of variation at very small as well as larger spatial scales.This distinguishes the M-RA from low-rank approaches that use relatively few basis functions.
  • Predictive-process components form orthogonal contributions, while spatial basis functions are block-orthogonal when their regions are disjoint.Orthogonality is exact for the predictive process and remainder in the probabilistic decomposition.
  • Proposition 1 establishes that the M-RA defines a valid Gaussian process with a nonnegative definite covariance function.
  • The basis functions are automatically selected for a given covariance function without restrictions on that covariance, including highly nonstationary cases.Figure 2 illustrates adaptation to spatially varying smoothness and lower-resolution basis placement.
  • For sufficiently large r, the M-RA can better approximate dependence between intermediate subregions than the 1-RA at the same computational complexity.The two approximations agree within the finest subregions and across different coarsest-level regions; extra M-RA basis functions address other pairs.

2.5 More on the choice of knots and partitions

The recommended design balances resolution depth, partitioning, and within-region knots: use small M and J, large r, and adapt partitions or knot placement to the observation geometry and covariance behavior.

  • Choose M and J as small and r as large as computational resources allow, subject to rJ^M ≥ n.
  • Equal-area recursive partitions are suitable for approximately uniform observations, whereas markedly nonuniform locations may require more complex schemes for fast inference.
  • Knots may use equidistant grids, but placing more knots near region boundaries can reduce the dependence omitted between neighboring subregions.The ignored cross-interior covariance can be small when the boundary distance is large or a screening effect holds.
  • For the one-dimensional exponential covariance without nugget, boundary knots make the M-RA exact because the screening effect holds exactly.Without a screening effect or in higher dimensions, exactness generally does not hold.
  • More elaborate partition and knot strategies, including clustering or reversible-jump Markov chain Monte Carlo, can be adopted when needed.Boundary effects may also be alleviated by carrying out several M-RA models.

3 Inference

Inference for the M-RA uses recursively computed local quantities and small matrix operations, enabling scalable likelihood evaluation, prediction, and distributed computation.

  • Parameter inference: Inference computes prior quantities and posterior distributions for the M-RA weight vectors at each resolution level.For fixed θ, the covariance function, basis functions, and covariance matrices are fixed before inference.
  • Parameter inference: The M-RA log-likelihood is a sum of log-determinants and quadratic forms involving only r×r matrices.This supports likelihood-based methods including maximum likelihood, Markov chain Monte Carlo, and particle filtering.
  • Spatial prediction: Spatial prediction is performed after parameter inference by constructing the posterior predictive distribution at prediction locations.The posterior predictive distribution has the same form as the prior M-RA process.
  • Distributed computing: Distributed inference assigns subsets of data and knots to nodes, with computations at each resolution carried out in parallel.Only connected nodes communicate in the tree-like computational structure.
  • Computational complexity: The M-RA inference algorithm has O(nr^2M^2) time and O(nrM) memory complexity in the nondistributed setting.With many distributed nodes, overall time is O(M^3r^3) and per-node memory is O(Mr^2), assuming communication does not dominate computation.
  • Computational complexity: With M increasing as log n, M-RA computation is quasilinear in n on one computer and polylogarithmic in distributed settings.The M-RA can retain fixed r while increasing M with n, unlike the 1-RA comparison described here.

4 Numerical comparisons and illustrations

The numerical studies compare the M-RA with full-scale approximations on simulated datasets and a large satellite TPW dataset. The M-RA achieves scalable computation while maintaining strong covariance approximation and long-range prediction performance.

  • Experimental design: The experiments compare the M-RA with full-scale approximations on simulated and real spatial datasets.The simulations include roughly 2 million one-dimensional observations and over 3.2 million two-dimensional observations; the real-data analysis uses 271,014 TPW measurements.
  • Simulation study: The one-dimensional simulations evaluate likelihood scores and computation times under fixed-domain and increasing-domain sampling as sample size grows.Datasets range from roughly 2,000 to 2 million observations, using both equally spaced subsets over the full domain and prefixes of the full dataset.
  • Simulation study: As n increases, the M-RA log-scores approach those of the true process and the slow full-scale approximation, while the fast full-scale approximation becomes relatively worse.The true process is computationally infeasible at large n; the fast full-scale approximation has fixed r and linear scaling, whereas the slow version has cubic scaling.
  • Simulation study: At n = 1,966,080, the 2-RA and 4-RA were roughly 8.7 and 11.8 times faster than the fastest comparable 1-RA while achieving equal or greater log-score.None of the tested 1-RAs achieved a log-score as high as the 8-RA.
  • TPW analysis: In the TPW analysis, the 6-RA produced by far the best long-range predictions, whereas 1-RA and block approximations showed strong blocky artifacts.The first three methods had similar computation times, maximum loglikelihood values, and short-range predictions; evaluations used held-out random observations and 5°×5° test regions.

5 Conclusions and Future Work

The conclusions present the M-RA as a flexible Gaussian-process approximation with multi-resolutional sparsity that supports scalable and distributed computation. They report favorable comparisons with a state-of-the-art full-scale approximation and outline extensions to complex and spatio-temporal models.

  • Conclusions: The M-RA approximates Gaussian processes with any covariance function using basis functions across multiple spatial resolutions.Its basis functions are chosen optimally for the specified covariance function.
  • Conclusions: A multi-resolutional block-sparse precision matrix enables scalable inference and distributed computations.The M-RA is a valid Gaussian process in its own right, rather than only a covariance-matrix approximation.
  • Conclusions: The M-RA provides a better approximation at the same computational complexity and time as the 1-RA, or a similar approximation at a fraction of the computational time.The 1-RA is the full-scale approximation of Sang et al. (2011), represented as the M-RA with M = 1.
  • Future work: Embedding the M-RA in hierarchical models can support complex measurement processes, non-Gaussian data, and data fusion across instruments.These extensions rely on the M-RA being a valid Gaussian process.
  • Future work: Because the entire joint posterior predictive distribution can be stored and propagated, the M-RA could support Kalman-filter-type inference in massive spatio-temporal state-space models.The paper identifies this as a potential alternative to ensemble Kalman filtering in certain settings.

A Proofs

The proofs establish that the M-RA remains a valid Gaussian process with nonnegative-definite covariance functions and supports sparse precision-based computation and prediction.

  • Gaussian-process validity: The M-RA is a Gaussian process because each finite collection of its values is a linear combination of jointly Gaussian basis-function weights.This follows from the multivariate normal distribution of the complete weight vector.
  • Covariance properties: The covariance functions of the residual, component, and M-RA processes are nonnegative definite.The proof obtains this property iteratively using the law of total variance.
  • Precision structure: The relevant precision matrix is sparse and block-structured, with only one nonzero block at the next resolution level.This structure underlies the recursive computations used in the proofs.
  • Precision structure: A variant of the Sherman-Morrison-Woodbury formula supports recursive manipulation of the sparse precision structure.The proof applies the formula across resolution levels.
  • Prediction: Spatial prediction uses a multivariate-normal posterior predictive distribution characterized by its posterior mean and covariance matrix.The proof derives these quantities using standard conditional-normal identities, total expectation, and total covariance.
Loading 1507.04789v2…