PaperlayerView exact revision on arXiv

Source-linked AI summary

A general framework for Vecchia approximations of Gaussian processes

Matthias Katzfuss, Joseph Guinness

arXiv:1708.06302v5stat.MEstat.CO

TL;DR

Exact GP inference is too costly for large datasets, especially for noisy spatial observations. The paper generalizes Vecchia approximation over latent and observed variables, unifies existing methods, and uses DAG-based sparsity analysis to develop SGV. SGV can improve approximation accuracy over standard Vecchia while retaining linear computational complexity, although prediction requires additional choices and the latent approach can scale quadratically in two-dimensional space.

  • Problem

    Exact GP inference has quadratic memory and cubic time costs, making it infeasible for large datasets; noisy spatial-data approximations lack a unified framework for comparison.

  • Method

    The paper applies Vecchia approximation to a vector containing latent GP realizations and noisy observations, analyzes the resulting DAGs and inference-matrix sparsity, and proposes SGV.

  • Results

    SGV can dramatically improve approximation accuracy over standard Vecchia while maintaining linear computational complexity; latent Vecchia can scale quadratically in two-dimensional space.

  • Takeaways & Limitations

    The authors recommend SGV with nugget or noise, standard Vecchia when noise is absent or nearly absent, and maxmin ordering in two-dimensional space.

  • Takeaways & Limitations

    Prediction requires addressing complex computational, ordering, and conditioning issues, and the framework assumes a GP with a specified covariance structure and additive independent Gaussian noise.

Abstract

from arXiv · show

Gaussian processes (GPs) are commonly used as models for functions, time series, and spatial fields, but they are computationally infeasible for large datasets. Focusing on the typical setting of modeling data as a GP plus an additive noise term, we propose a generalization of the Vecchia (1988) approach as a framework for GP approximations. We show that our general Vecchia approach contains many popular existing GP approximations as special cases, allowing for comparisons among the different methods within a unified framework. Representing the models by directed acyclic graphs, we determine the sparsity of the matrices necessary for inference, which leads to new insights regarding the computational properties. Based on these results, we propose a novel sparse general Vecchia approximation, which ensures computational feasibility for large spatial datasets but can lead to considerable improvements in approximation accuracy over Vecchia's original approach. We provide several theoretical results and conduct numerical comparisons. We conclude with guidelines for the use of Vecchia approximations in spatial statistics.

1 Introduction

The paper addresses the computational infeasibility of exact GP inference for large datasets by developing a general Vecchia framework for noisy spatial data. It unifies existing approximations, analyzes their computational structure through DAGs, and motivates a new sparse approximation.

  • Exact GP inference requires quadratic memory and cubic time in the number of observations, limiting use for datasets with tens of thousands or more observations.
  • Existing GP approximations use sparse covariance, sparse precision, low-rank, and algorithmic constructions to reduce computational demands.
  • Vecchia approximation replaces a high-dimensional joint distribution with univariate conditionals on small subsets of previous observations, producing a sparse precision Cholesky factor.
  • The generalized framework applies Vecchia approximation to latent variables and noisy observations, contains several popular GP approximations as special cases, and supports efficient likelihood computation.
  • DAG representations are used to study inference-matrix sparsity and computational complexity, motivating a sparse general Vecchia approximation and comparisons among methods.

2 A general Vecchia approach

The general Vecchia approach approximates a noisy GP by applying Vecchia factorization to an interleaved vector of latent variables and observations. Its approximation is determined by locations, grouping, ordering, conditioning, and the division between latent and observed conditioning variables.

  • 2.1 Noisy observations of a Gaussian process: The model assumes a zero-mean GP on a continuous domain, partitioning locations into subvectors whose process values form the latent vector y.
  • 2.1 Noisy observations of a Gaussian process: Observed vectors equal latent process values plus independent Gaussian nugget noise, with only a subset of observation vectors potentially observed.
  • 2.2 Review of Vecchia’s approximation: Exact likelihood evaluation becomes infeasible because the observed density incurs O(n_z^2) computational cost for n_z individual observations.
  • 2.2 Review of Vecchia’s approximation: Vecchia approximation replaces each observed history with a smaller conditioning index vector, commonly selecting nearby observations.
  • 2.3 The general Vecchia framework: General Vecchia applies the approximation to x = y ∪ z_o, interweaving latent vectors and their corresponding observed vectors in the ordering.
  • 2.3 The general Vecchia framework: For latent y_i, q_y(i) selects latent conditioners and q_z(i) selects observed conditioners; latent conditioning may improve accuracy but increase computational expense.
  • 2.3 The general Vecchia framework: The approximation is specified by location selection, partitioning, ordering, conditioning indices, and partitioning each conditioning set into latent and observed variables.
  • 2.4 Connections to directed acyclic graphs: Representing variables as DAG vertices makes conditioning relationships explicit and supports analysis of the sparsity of matrices needed for inference.

3 Existing methods as special cases

Many established GP approximations arise from particular choices within the general Vecchia framework. These choices differ in conditioning variables, grouping, ordering, spatial partitioning, and whether latent variables are introduced.

  • Standard Vecchia and extensions: Standard Vecchia uses singleton vectors, coordinate ordering, and conditioning only on previous observations.
  • Standard Vecchia and extensions: The NNGP conditions only on latent variables, while observed locations may condition only on variables in a latent knot set.
  • Independent blocks: Independent-block approximations use empty conditioning sets, treating spatial subvectors independently with parallelizable and inexpensive computations.
  • Latent vector autoregressive models: Latent vector autoregressive and state-space models condition on the latest m sets of latent variables and typically use Kalman filtering and smoothing.
  • Predictive-process approximations: The MPP and FSA-block use a common knot vector for conditioning, while FSA-block additionally groups remaining variables by spatial region.
  • Multiresolution approximation: The MRA iteratively partitions the domain into subregions and conditions on latent variables associated with containing regions; FSA and MPP are special cases of it.
  • Composite likelihood: Composite likelihood methods may use neighboring pairs or more general conditioning sets, but are not generally guaranteed to become exact or define valid joint densities.

4 Inference and computations

The framework represents general Vecchia approximations with sparse matrix factors whose structure determines inference cost. Analytical integration over latent variables enables likelihood-based inference for very large datasets, with complexity characterized through DAG-induced sparsity.

  • Matrix representations: The general Vecchia density has a sparse upper-triangular reverse-Cholesky factor U for the precision matrix.U = rchol(bC−1).
  • Integrated likelihood: Analytical integration over latent variables yields an observed-data likelihood using W = U_Y U_Y′, V = rchol(W), and transformed observations.This avoids numerical integration over a high-dimensional latent vector.
  • Inference: The integrated likelihood can be evaluated quickly for parameter values, supporting frequentist and Bayesian inference on very large datasets.The paper reports similar SGV parameter estimates to exact maximum likelihood estimates in a spatial-range simulation.
  • Interpretation: The approximation defines a valid probability model, and its approximation error disappears when m = n − 1.Inference is exact under the approximating model itself.
  • Prediction: Prediction requires additional ordering, conditioning, and fast-summary strategies, which are referred to separate work.The posterior y|z is available, but efficient prediction is not developed here.
  • Sparsity and DAGs: The sparsity of U, W, and V is determined by parent relations and paths in the Vecchia DAG.For V, nonzeros depend on paths through latent variables with observed descendants.
  • Computational complexity: U has O(nmr) nonzero elements and can be computed in O(nm^3r^2) time under bounded conditioning size m.The stated scaling assumes equal-order block sizes and conditioning vectors containing at most m subsets.
  • Computational complexity: W has at most O(nrm^2) nonzero elements, while reverse-Cholesky cost depends on squared per-column nonzero counts in V.The DAG sparsity structure therefore determines inference complexity for each Vecchia specification.

5 Sparse general Vecchia (SGV) approximation

SGV partitions conditioning variables between latent and observed components while enforcing a DAG-compatible structure. It improves approximation accuracy over standard Vecchia while retaining its favorable computational scaling, unlike latent Vecchia in spatial settings.

  • Methods: Standard Vecchia conditions only on observed vectors, latent Vecchia only on latent vectors, and SGV partitions each conditioning vector between both.SGV imposes an additional ordering constraint on jointly latent-conditioned indices.
  • SGV construction: SGV ensures a perfect-graph DAG and can include as much latent conditioning as allowed by its compatibility rule.Standard Vecchia is one special case of SGV.
  • Approximation accuracy: For fixed grouping, ordering, and conditioning indices, joint-distribution accuracy is ordered latent Vecchia, then SGV, then standard Vecchia.The ordering concerns the approximation to the joint distribution of x and does not guarantee the same ordering for observed-data KL divergence.
  • Computational complexity: Latent Vecchia can require O(n^3−2/dm^2/d) time to compute V on d-dimensional grids.The complexity approaches cubic in n as spatial dimension increases.
  • Computational complexity: SGV guarantees at most mr off-diagonal elements per V column and computes V in O(nm^2r^2) time.Its overall computational complexity is the same as standard Vecchia.
  • Summary: SGV improves approximation accuracy over standard Vecchia while retaining linear computational complexity in n.Latent Vecchia improves accuracy but can become infeasible for large n.

6 Ordering and conditioning

Ordering and conditioning choices affect Vecchia approximation quality and computational behavior, especially beyond one spatial dimension. The paper recommends adapting conditioning size m to smoothness, ordering, locations, and available resources.

  • Ordering: Ordering has negligible effect in some studies but can materially affect Vecchia approximation quality in others.Coordinate ordering is common, while higher-dimensional spatial orderings are not uniquely obvious.
  • Conditioning size: The required conditioning size m depends on the covariance smoothness, spatial dimension, ordering, observation locations, and other factors.The paper treats m as a key design choice for a given ordering.
  • One-dimensional guidance: For a one-dimensional Matérn process with smoothness ν = 0.5, left-to-right latent conditioning with m = 1 gives an exact approximation.The process is Markov of order 1 in this case.
  • Practical guidance: The authors suggest starting with small m and increasing it using warm-started parameter estimates until convergence or resource limits are reached.This recommendation applies when the necessary m is unknown.
  • Conditioning schemes: Nearest-neighbor conditioning is common, while same-conditioning-set schemes may benefit from placing a coarse spatial grid early in the ordering.More elaborate conditioning schemes are also possible.

7 Numerical study

The numerical study compares Vecchia variants across one- and two-dimensional spatial settings, emphasizing approximation accuracy, ordering and conditioning choices, sparsity, and computation. SGV improves substantially over standard Vecchia in several settings while retaining linear scaling and outperforming competing linear-scaling methods after accounting for complexity.

  • One-dimensional domain: Latent Vecchia and equivalent SGV performed much better than standard Vecchia on a one-dimensional Matérn GP with coordinate ordering.The setting used nz = 100 and λ = 0.9; the results also numerically supported approximate screening for latent Vecchia when m > ν.
  • Two-dimensional domain: In two dimensions, KL divergences generally ordered latent Vecchia best, SGV next, and standard Vecchia worst.The comparison used Matérn covariances on the unit square with nz = 6,400 and λ = 0.9; methods were equivalent for SNR = ∞.
  • Ordering and conditioning: Maxmin ordering often produced tremendous improvements over coordinate ordering, except for standard Vecchia, whose results were similar under both orderings.This ordering effect was observed in the two-dimensional experiments.
  • Ordering and conditioning: For very smooth covariances, conditioning on the first m variables under maxmin ordering strongly improved SGV over nearest-neighbor conditioning.The comparison used nz = 400, m = 16, and λ ≈ 2; SGV and latent Vecchia were equivalent for first-m conditioning.
  • Computational feasibility: SGV kept nonzero elements per column in V at or below m, yielding linear scaling in n, whereas latent Vecchia's denser V required roughly O(n^2) computation.For nz around 100,000, latent Vecchia with MMD ordering was roughly two orders of magnitude slower than SGV with reverse ordering.
  • Large datasets: SGV outperformed standard Vecchia, MRA, and independent blocks in large-n comparisons even after adjusting for differences in computational complexity.The four methods all scaled linearly in this comparison, which used maxmin ordering where applicable and SGV with m = 40 for the reference loglikelihood.

8 Conclusions and guidelines

The paper develops a general Vecchia framework and proposes SGV, which can improve approximation accuracy while retaining linear computational complexity. It recommends SGV with nugget or noise and standard Vecchia when noise is negligible.

  • The general framework applies Vecchia approximations to latent GP realizations and noisy observations, encompassing several commonly used GP approximations.
  • SGV can dramatically improve approximation accuracy over standard Vecchia while maintaining linear computational complexity.
  • Latent Vecchia can scale quadratically in data size in two-dimensional space.
  • SGV is recommended when spatial data include nugget or noise, whereas standard Vecchia is recommended when the noise term is zero or almost zero.
  • In two-dimensional space, maxmin ordering is recommended, and the conditioning size m can be increased until inference converges or computational resources are exhausted.
  • The methods are implemented in the R package GPvecchia, with extensions for prediction and further computational details available elsewhere.

A Vector notation

The paper establishes notation for ordered vectors, subvectors, unions, intersections, matrices, and directed acyclic graphs. It also defines d-separation and its conditional-independence implication for multivariate normal variables.

  • Vectors are ordered lists of same-type elements equipped with union and intersection operations; integer vectors use non-bold lowercase notation, while real-valued vectors use bold lowercase notation.
  • Subvectoring uses index vectors and preserves the ordering specified by those indices.
  • The index function returns an element’s position within a vector and acts as an inverse of the union operation.
  • Matrices are treated as two-dimensional vectors with context-defined blocks, and functions are vectorized over location vectors.
  • In a directed graph, vertices connected by directed edges have parent-child relationships determined by edge direction.
  • For multivariate normal variables, d-separation of two sets by a third implies their conditional independence given that third set.

C Same conditioning sets (SCS)

SCS assigns a common conditioning set to variables, yielding sparsity and computational advantages. Its accuracy in two-dimensional space may require conditioning sets whose size grows with the number of spatial observations.

  • SCS gives every latent variable the same conditioning vector, a strategy used by MPP and FSA with different conditioning-set sizes.
  • Latent Vecchia automatically follows the SGV rules for SCS, guaranteeing the corresponding sparsity structure.
  • SCS provides exact marginal distributions for the latent variables and preserves the sparsity structure in V^-1, supporting fast prediction and Kalman-filter-type inference.
  • Table 1 reports simulated range-parameter estimation and 95% confidence intervals for MSE differences relative to the method in the last row.
  • SCS and MRA may require r1 = O(√nz) for accurate two-dimensional approximations, resulting in time complexity O(n^3/2).

D Comparison to composite likelihood

The paper compares SGV with full-conditional and pairwise-block composite likelihood methods for estimating a spatial covariance range parameter. In the reported simulations, SGV performs better than the tested composite-likelihood settings, including a large-block PBL configuration.

  • The simulations compare SGV with full-conditional likelihood and pairwise-block likelihood for maximum-likelihood estimation of the unknown range α.
  • The FCL study uses 300 datasets on a 30 × 30 grid with true range α = 10 and compares exact likelihood, FCL, and SGV with m = 10, 15, and 20.
  • The PBL study uses 300 datasets on a 100 × 100 grid with α = 30 and evaluates several PBL and SGV settings.
  • Even SGV with m = 20 performed better than PBL with 100 blocks of size 100 each.

E Illustration of Bayesian inference

The Bayesian-inference illustration evaluates SGV-based posterior and predictive distributions on simulated 30×30 and 100×100 grid settings. On held-out locations, the resulting posterior predictive intervals were well calibrated.

  • SGV with m = 30 was used for Bayesian inference under the settings described in Appendix D.
  • On the 30×30 grid, the study compared exact and SGV-implied posteriors for α and predictive distributions at an unobserved center point.
  • For the 100×100 grid, 9,000 randomly selected locations were observed and the remaining GP realizations served as test data.
  • 79.8% of posterior 80% intervals covered the true values at 1,000 held-out test locations, indicating well-calibrated predictive distributions.

F Proofs

The proofs establish how conditioning choices and graph structure determine approximation accuracy and matrix sparsity. They show that SGV retains sparse computation, with complexity governed by n, m, and r.

  • The approximate covariance precision is factored through a nonsingular upper-triangular matrix U, yielding bC−1 = UU′.
  • Conditional-independence and d-separation arguments characterize zeros in the precision and Cholesky factors by blocking paths in the associated graph.
  • Adding variables to Vecchia conditioning vectors cannot increase KL divergence from the true model, supporting the conditioning relations used in Proposition 4.
  • Under lexicographic ordering on a regular grid, graph paths imply nonzero entries in V across a neighborhood whose width scales as O(m1/d).
  • For SGV, V has at most mr off-diagonal nonzero elements per column, which determines its sparse structure.
  • SGV requires O(nm2r2) time for W-related computations, O(nm3r2) time for U, and O(nmr) memory for U.
Loading 1708.06302v5…