PaperlayerView exact revision on arXiv

Source-linked AI summary

Vecchia approximations of Gaussian-process predictions

Matthias Katzfuss, Joseph Guinness, Wenlong Gong, Daniel Zilber

arXiv:1805.03309v4stat.MEstat.CO

TL;DR

Direct Gaussian-process prediction is computationally infeasible for many large datasets, despite its flexibility and uncertainty quantification. This paper studies spatial prediction through a general Vecchia framework and its existing and novel variants. It establishes linear-complexity methods and gives recommendations tailored to dimension, accuracy, and computational trade-offs.

  • Problem

    Direct GP computation is cubic in data size, limiting its use for large datasets, while scalable joint prediction at observed and unobserved locations requires suitable approximations.

    Source: §1 · pp. 1–2

  • Method

    The paper extends a general Vecchia approximation framework to spatial prediction and studies choices of ordering and conditioning for marginal and joint predictive distributions.

    Source: §1 · pp. 2–3, §3.1 · p. 5

  • Results

    The new methods guarantee linear time and memory complexity for fixed conditioning-set size, and the paper recommends LF-auto in one dimension and generally RF-full in two dimensions.

    Source: §1 · p. 3, §7 · p. 23

  • Takeaways & Limitations

    Method choice can trade accuracy against computational cost: RF-stand has computational advantages when prediction locations are much fewer than observations, while LF-full can be very accurate but does not scale linearly.

    Source: §7 · p. 23

    LF-ind assumes conditional independence of prediction locations given observed latent values, which can make joint predictions inaccurate, and neither LF-ind nor LF-full guarantees linear computational complexity.

    Source: §4.2 · p. 10, §4.3.3 · p. 14

Abstract

from arXiv · show

Gaussian processes (GPs) are highly flexible function estimators used for geospatial analysis, nonparametric regression, and machine learning, but they are computationally infeasible for large datasets. Vecchia approximations of GPs have been used to enable fast evaluation of the likelihood for parameter inference. Here, we study Vecchia approximations of spatial predictions at observed and unobserved locations, including obtaining joint predictive distributions at large sets of locations. We consider a general Vecchia framework for GP predictions, which contains some novel and some existing special cases. We study the accuracy and computational properties of these approaches theoretically and numerically, proving that our new methods exhibit linear computational complexity in the total number of spatial locations. We show that certain choices within the framework can have a strong effect on uncertainty quantification and computational cost, which leads to specific recommendations on which methods are most suitable for various settings. We also apply our methods to a satellite dataset of chlorophyll fluorescence, showing that the new methods are faster or more accurate than existing methods, and reduce unrealistic artifacts in prediction maps.

1 Introduction

Gaussian processes support flexible spatial prediction with uncertainty quantification, but their cubic computational cost makes direct use infeasible for many large datasets. This paper extends a general Vecchia framework to spatial prediction and studies accuracy, computation, and method choices for observed and unobserved locations.

1 Introduction

  • Direct GP computation scales cubically in data size, while sparsity-based approaches may not scale linearly in higher dimensions and local methods do not scale well to joint predictions at many locations.Low-rank approaches are described as poorly suited to capturing fine-scale dependence.Source: §1 · pp. 1–2
  • Prior Vecchia prediction approaches include one-at-a-time prediction, conditioning on observed variables, and domain-partition methods, with cited limitations including simulation error, noisy-data approximation accuracy, and partition-boundary artifacts.The paper compares its methods with Bayesian inference and prediction approaches and discusses these connections.Source: §1 · p. 2
  • The paper extends the general Vecchia framework to marginal and joint spatial predictions at observed and unobserved locations, addressing ordering and conditioning choices needed for accurate joint prediction.It also studies fast computation of summaries of joint predictive distributions.Source: §1 · pp. 2–3
  • The new methods guarantee linear time and memory complexity in the number of data points and predictions for fixed conditioning-set size.The framework supports both frequentist and Bayesian inference on hyperparameters, while the study focuses on predictions conditional on hyperparameters.Source: §1 · p. 3
  • Joint predictive distributions are needed to quantify uncertainty in spatial averages, totals, and other downstream analyses of multiple predictions.The paper also considers latent-process prediction at observed locations for spatial smoothing and Vecchia-Laplace approximations of non-Gaussian spatial data.Source: §1 · p. 3

2 Exact Gaussian-process prediction

The section formulates exact GP prediction for noisy observations and distinguishes observed locations from unobserved prediction locations. It defines the conditional predictive distribution and notes the computational burden of obtaining joint uncertainty summaries.

2 Exact Gaussian-process prediction

  • The model treats responses as independent noisy versions of latent GP values, with covariance C = K + D and diagonal noise variances in D.The locations are divided into observed indices o and unobserved prediction indices p.Source: §2 · pp. 3–4
  • For fixed hyperparameters, prediction uses the conditional distribution of the latent process given observed responses; predicting noisy responses adds their noise covariance to the predictive covariance.The paper describes both point-mass hyperparameter estimation and Bayesian parameter averaging as routes to prediction.Source: §2 · p. 4
  • Joint predictive distributions are needed for uncertainty in linear combinations and for posterior simulations, whereas individual point predictions use marginal distributions.The section states that obtaining the full covariance matrix and even its diagonal can be computationally expensive.Source: §2 · p. 4
  • For large numbers of observations or locations, exact GP prediction is computationally infeasible, so approximations or simplifying assumptions are necessary.The passage identifies the issue for large nO or n.Source: §2 · p. 4

3 The general Vecchia framework for GP prediction

The general Vecchia framework approximates the joint distribution by restricting each variable's conditioning set, then derives predictive quantities from sparse precision-factor representations. Ordering and conditioning choices shape both approximation accuracy and computational cost.

3.1 Definition of the framework

  • The framework applies Vecchia factorization to x = zo ∪ y, conditioning each variable on a subset of preceding variables that may be selected by spatial proximity.Using all preceding variables recovers the exact distribution; bounding conditioning-set size by m can reduce computation to decompositions of m × m matrices.Source: §3.1 · pp. 4–5
  • Approximation accuracy depends on variable ordering and conditioning sets, while computational efficiency depends in part on sparsity of the posterior precision matrix and its Cholesky factor.Conditioning on latent rather than response variables can be more accurate but also more computationally expensive.Source: §3.1 · p. 5

3.2 Matrix representations

  • The matrix representation uses U as an upper-lower Cholesky factor of the joint precision matrix and W as the latent-variable posterior precision matrix.U has at most m off-diagonal nonzeros per column and can be computed in O(nm^3) time.Source: §3.2 · pp. 5–6

3.3 General Vecchia predictions

  • General Vecchia prediction approximates f(y|zo) from the joint approximation, with posterior mean and marginal variances computed using sparse-factor operations.Selected inversion computes covariance entries for pairs corresponding to nonzero entries of W, including the prediction variances.Source: §3.3 · p. 6
  • The framework also supports variances of linear combinations and conditional simulation, but only moderate numbers of linear combinations are generally computationally feasible.Computing these prediction tasks requires V = rchol(W), whose factorization cost depends on its sparsity.Source: §3.3 · p. 6

4 Specific methods and their properties

The paper compares response-first and latent-first Vecchia orderings and conditioning schemes for spatial GP prediction. Their accuracy and computational costs depend on these choices, with some methods guaranteeing linear complexity and others trading accuracy against cost.

4 Specific methods and their properties

  • The framework compares response-first and latent-first orderings and their special cases, providing a systematic basis for evaluating Vecchia methods for GP prediction.The methods are summarized in Table 1 and illustrated as directed acyclic graphs in Figure 1.Source: §4 · p. 7, Table 1 · p. 7, Figure 1 · p. 8

4.1 Response-first ordering

  • RF-full conditions on eligible previously ordered variables, whereas RF-stand restricts prediction variables to conditioning on prediction and observed-response variables.RF-stand can reduce prediction cost when only unobserved locations are needed, at some loss of accuracy relative to RF-full.Source: §4.1 · p. 9

4.2 Latent-first ordering

  • Among latent-first methods, LF-auto uses the preceding m latent variables, LF-full allows broader conditioning, and LF-ind excludes prediction locations from conditioning.LF-full is usually the most accurate in the examples considered but can be computationally demanding; MRA uses domain partitioning to ensure sparsity and linear complexity.Source: §4.2 · pp. 10–11

4.3.1 Zero noise

  • With zero nugget, RF-full, RF-stand, and LF-full become equivalent, as do LF-ind and RF-ind.Source: §4.3.1 · p. 11

4.3.2 Approximation accuracy

  • Approximation accuracy generally improves as conditioning-set size m increases, but RF-ind and LF-ind converge to exact marginals rather than the exact joint predictive distribution.LF-auto is exact for a one-dimensional exponential-covariance GP when m ≥ 1 and nearly exact for Matérn smoothness ν when m > ν.Source: §4.3.2 · p. 13

4.3.3 Computational complexity

  • For response-first methods, V is a submatrix of U; LF-auto and MRA also avoid fill-in through their special conditioning structures.These structures bound the off-diagonal nonzeros per column of V by m, supporting the linear-complexity result for fixed m.Source: §4.3.3 · pp. 13–14
  • LF-auto, MRA, and all response-first methods have linear prediction complexity in the number of locations, while LF-full and LF-ind have no such guarantee.For LF-full and LF-ind, Cholesky fill-in can make computation grow quadratically with the number of observations on regular two-dimensional grids.Source: §4.3.3 · p. 14

4.3.4 Consistent framework

  • Prediction locations can be appended after likelihood inference without changing the distribution of existing variables, and existing matrices can be reused or augmented.The prediction locations therefore need not be known when training the model.Source: §4.3.4 · p. 15

5 Simulation study

The simulation study compares Vecchia prediction methods on one- and two-dimensional problems, timing their computations and evaluating large-data and benchmark settings. Accuracy depends on the method and conditioning scheme, with trade-offs between joint prediction quality and computational scalability.

5 Simulation study

  • The simulation comparison used random observation locations and an equidistant prediction grid, with known hyperparameters to isolate approximation accuracy.It evaluated joint and marginal KL divergences for predictions and observations across repeated simulations.Source: §5 · p. 15

5.1 Numerical comparison in 1-D

  • LF-auto was exact for ν = 0.5 at any m ≥ 1 and was much more accurate than other approaches for prediction locations when ν = 1.5.RF-ind and LF-ind plateaued in joint prediction accuracy as m increased, while improving on marginal measures.Source: §5.1 · p. 16, Figure 4 · p. 17

5.2 Numerical comparison in 2-D

  • LF-full performed best in two dimensions but was not computationally scalable, while RF-full performed well on both joint and marginal accuracy measures.RF-ind and LF-ind did not converge to zero joint prediction KL divergence, though they were more competitive on marginal measures.Source: §5.2 · p. 16, Figure 5 · p. 18

5.3 Timing comparison

  • On a unit square, U time grew roughly linearly with n; for LF-ind and LF-full, V time using a fill-reducing permutation grew roughly as O(n^(3/2) and failed for large nO due to memory limitations.For response-first methods, V computation was negligible relative to U.Source: §5.3 · p. 16

5.4 Comparison for large n

  • RF-full was much more accurate than the other scalable response-first methods across the large-data settings considered.The comparison used approximate KL divergences for the joint predictive distribution.Source: §5.4 · p. 17, Figure 7 · p. 19

5.5 Heaton comparison

  • RF-full and RF-stand outperformed the other methods on the Heaton comparison; only MRA achieved comparable scores, at much larger conditioning-set size and computation time.RF-ind’s joint log score was much worse than RF-full’s and RF-stand’s.Source: §5.5 · p. 19, Table 2 · p. 20

6 Application to satellite data

The satellite application compares scalable Vecchia prediction methods on dense, swath-patterned chlorophyll fluorescence data. RF-full performs best across the reported settings, while RF-ind produces noisier maps with visible streaks and discontinuities.

6 Application to satellite data

  • The study analyzed 245,236 observations collected over the contiguous United States from August 1–31, 2018, and restricted modeling to a spatial process because temporal change appeared limited.Predictions were computed on a 24,407-location grid over the contiguous United States and an 18,576-location grid over Texas.Source: §6 · pp. 20–21
  • The residual field used an isotropic Matérn covariance with estimated variance 0.1097, range 100.8 km, smoothness 0.0982, and noise variance 0.1869.The model included 50 Gaussian basis functions with 637 km range to capture long-range variability.Source: §6 · pp. 20–21
  • RF-ind predictions were noisier and showed streaks and discontinuities near the Texas swaths, which remained visible at m = 60.The authors conjecture that nearby prediction locations can condition on different observation sets when they are nearly equidistant from separate, densely sampled swaths.Source: §6 · p. 21, Figure 9 · p. 22
  • Cross-validation assessed short-range prediction on random held-out sets and long-range prediction on held-out swaths using RMSE and total log score.Each experiment used 10 test sets; the average held-out swath contained 3,493 locations.Source: §6 · pp. 21–22
  • RF-full performed best in all settings; for random tests, RF-stand and local kriging needed about m = 20 to match RF-full at m = 10, while log scores showed stronger uncertainty-quantification differences.RF-ind was not competitive for long-range swath predictions; computation scales cubically in m.Source: §6 · p. 22, Figure 10 · p. 23

7 Conclusions

The conclusions recommend different Vecchia prediction methods according to domain dimension, accuracy needs, and computational constraints. The appendices provide vector-index notation and formulas for computing U.

7 Conclusions

  • The paper extends Vecchia approximations to GP prediction in a general framework that includes existing and novel approaches, with theoretical and numerical study of accuracy and computational properties.The methods can also be extended to unknown hyperparameters.Source: §7 · p. 23
  • For two-dimensional domains, RF-full is generally recommended because it scales linearly and performs well across accuracy measures; other methods trade accuracy, joint-distribution fidelity, or linear scaling for computational advantages.For one-dimensional domains, LF-auto is recommended without qualification; RF-stand can help when prediction locations are much fewer than observations, while RF-ind ignores joint predictive dependence.Source: §7 · p. 23
  • The methods are implemented in the R package GPvecchia, but their properties in dimensions above two warrant further investigation.The conclusion notes follow-up work on higher-dimensional emulation and non-Gaussian spatial data.Source: §7 · p. 23

A Vector and indexing notation

  • Vector-index notation defines subvectors by ordered index vectors, and unions or intersections preserve a specified ordering.For example, an index vector (4, 1, 2) selects elements in that order.Source: A Vector and indexing notation · p. 24

B Computing U

  • In computing U, g(i) denotes the indices of variables on which x_i conditions, while C denotes covariance under the true model.The passage specifies covariance for latent and observed variables using K and the observation-noise variance.Source: B Computing U · p. 24

C Prediction with unknown parameters

The appendix describes likelihood-based parameter inference for Vecchia prediction and a numerical issue in Cholesky-based computation. It contrasts plug-in and Bayesian parameter handling and explains how rounding can create numerical nonzeros.

C Prediction with unknown parameters

  • The Vecchia likelihood can have low evaluation cost, with linear cost in n guaranteed under conditions on the conditioning sets g(i).The likelihood integrates the approximate joint density over the latent process.Source: C Prediction with unknown parameters · p. 24
  • Likelihood-based inference can estimate unknown parameters by maximizing the Vecchia likelihood, then compute predictive summaries at the estimate, but this plug-in approach ignores parameter-estimation uncertainty.The passage also describes Bayesian Metropolis–Hastings sampling, which yields a Gaussian-mixture prediction over parameter draws.Source: C Prediction with unknown parameters · pp. 24–25

D Numerical nonzeros in V

  • Cholesky rounding errors can introduce numerical nonzeros in V where its theoretical value is zero; extracting V as a submatrix of U avoids them.For some j < i, W_ji can be nonzero even when V_ji is zero, with the equality involving the Cholesky sum potentially failing numerically.Source: D Numerical nonzeros in V · p. 25

D.1 Implications for selected inverse

Computing selected inverses from a copied submatrix of V may not yield exact posterior variances. Exact variances are available through a separate calculation whose cost depends on the prediction-set size.

D.1 Implications for selected inverse

  • The selected inverse from copied-submatrix V is not guaranteed to give exact posterior variances unless V is zero-padded, which adds cost.In experiments, its additional approximation error was negligible relative to Vecchia error and went to zero as m became large enough for the Vecchia approximation to be accurate.Source: §D.1 · p. 25
  • Exact variances for RF-full and RF-stand require O(nmnP) time, so overall complexity is not increased when nP = O(m^2).The passage gives the exact variance calculation as a diagonal of products involving V^-1.Source: §D.1 · p. 25

E Proofs

The appendix proofs use conditional-independence and conditioning-set arguments to establish propositions. Supplementary results show block Cholesky forms improve sparsity, while reverse ordering raises computation times for LF methods.

E Proofs

  • Proposition 1's proof uses zero-mean, unit-variance standardized residuals and conditioning-set inclusion to establish its parts.The proof applies the law of total variance and also treats conditional densities for subvectors.Source: E Proofs · p. 26
  • Proposition 2's proof uses conditional independence given yj and conditioning-set inclusion; response-first methods have nested sets as m increases.For a one-dimensional GP with exponential covariance, left-to-right conditioning on all preceding variables recovers the exact density.Source: E Proofs · p. 26

S1 Sparsity and computation time

  • Block forms (5) or (7) avoid numerical nonzeros and increase V sparsity except for RF-ind, whose W and V are diagonal.Figure S1 contrasts the block forms with brute-force Cholesky under reverse row-column ordering.Source: S1 Sparsity and computation time · p. 30, Figure S1 · p. 30
  • Reverse ordering of Woo without a fill-reducing permutation algorithm leads to higher computation times for LF-ind and LF-full.Time for computing U is similar across methods; for RF methods, computing V is negligible using (5).Source: S1 Sparsity and computation time · p. 30

S2 Vecchia predictions based on interweaved ordering

Interweaved predictions (IWP) extend the SGV likelihood approximation to GP prediction by interweaving observed responses and latent variables. The resulting precision-factor sparsity yields linear computational complexity.

S2 Vecchia predictions based on interweaved ordering

  • IWP interweaves observed responses yo and latent variables zo under OP ordering, followed by prediction variables yp.For observed locations, its conditioning rules follow SGV; it tries to maximize latent conditioning while using zj where needed for linear complexity.Source: S2 Vecchia predictions based on… · p. 30
  • IWP is a natural extension of the SGV likelihood approximation, and the two approaches provide a consistent framework for inference and prediction.Its implied likelihood for zo is equivalent to that of SGV.Source: S2 Vecchia predictions based on… · pp. 30–31
  • For IWP, V has at most m off-diagonal nonzero elements in each column.The sparsity result follows from Vji being zero unless j=i or j is in qy(i).Source: S2 Vecchia predictions based on… · p. 31
  • This sparsity guarantees linear computational complexity in n for prediction using IWP.For RF methods, computing V is negligible using (5); for LF methods, V is computed using (7).Source: S2 Vecchia predictions based on… · p. 31, Figure S2 · p. 31

S3 Numerical comparison to IWP and MRA

The comparison considers RF-full, MRA, and IWP, which all scale linearly in n. IWP is most accurate for small m, but its accuracy improves little or can decrease at intermediate m.

S3 Numerical comparison to IWP and MRA

  • RF-full, MRA, and IWP all scale linearly in n in the comparison.The comparison uses a unit square with nO=nP=4,900 and effective range 0.15.Source: S3 Numerical comparison to IWP and MRA · p. 31
  • The MRA comparison uses default GPvecchia tuning, which could be improved, and does not adjust for MRA's O(nm^2) implementation complexity.That complexity is a factor of m lower than the complexity of RF-full and IWP.Source: S3 Numerical comparison to IWP and MRA · p. 31
  • IWP is the most accurate method for small m and is exact at m=n−1, but its accuracy improves little and sometimes decreases for intermediate m around 10 to 50.For some observed-location pairs, sparsity rules can require conditioning on zj instead of yj, which may make the approximation to their joint distribution poor.Source: S3 Numerical comparison to IWP and MRA · pp. 31–32

S4 Satellite data: supplementary plots

The supplementary satellite-data figures provide additional prediction and uncertainty plots. They include the original SIF data, residuals after trend removal, and maps over the contiguous USA and Texas.

S4 Satellite data: supplementary plots

Loading 1805.03309v4…