Vecchia Gaussian Processes: on probabilistic and statistical properties

Gaussian Processes (GPs) are widely used to model dependencies in spatial statistics and machine learning. However, exact inference is computationally intractable for GP regression, with a time complexity of $O(n^3)$. The Vecchia approximation scales up computation by introducing sparsity into the spatial dependency structure, represented by a directed acyclic graph (DAG). Despite its practical popularity, this approach lacks rigorous theoretical foundations, and the choice of DAG structure remains an open problem. In this paper, we systematically study the Vecchia approximation of the popular, isotropic Matérn GP as standalone stochastic process and uncover key probabilistic and statistical properties. We propose selecting parent sets as norming sets with fixed cardinality in the Vecchia approximation. On the probabilistic side, we show that the conditional distributions of Matérn GPs, as well as their Vecchia approximations, can be characterized by polynomial interpolations. This enables us to establish several results on small ball probabilities and the Reproducing Kernel Hilbert Spaces (RKHSs) of Vecchia GPs. Building on these probabilistic results, we prove that in the nonparametric regression model, the corresponding posterior contracts around the truth at the optimal minimax rate, both under oracle rescaling and hierarchical tuning of the prior. We illustrate the theoretical findings through numerical experiments on synthetic datasets. Our core algorithms are implemented in C++ with an R interface.

💡 Research Summary

This paper provides a rigorous theoretical investigation of Vecchia approximations for Gaussian Processes (GPs), focusing on isotropic Matérn kernels, and treats the resulting Vecchia GP as an independent stochastic process rather than merely an approximation of the “mother” GP. The authors address two central methodological questions: (1) what is the minimal size of parent sets required for optimal statistical inference, and (2) how should these parent sets be selected.

The key probabilistic contribution is the characterization of conditional distributions in a Vecchia GP via local polynomial interpolation. For a Matérn GP with smoothness parameter α, the conditional mean of a point given its parent set equals a polynomial of degree ⌊α⌋ fitted on the parents, while the conditional variance equals the squared interpolation error. This links the Vecchia construction to classical Gaussian interpolation results but strengthens them to uniform convergence rates. The authors introduce the notion of a norming set: a finite collection of points that controls the supremum norm of any polynomial of degree ≤⌊α⌋ on the whole domain. By choosing parent sets that are norming sets with a fixed cardinality roughly ⌊α⌋ + d·⌊α⌋, they guarantee that the interpolation error decays as O(h^α), where h denotes the mesh size.

Because the Vecchia GP is defined through conditional distributions, standard tools based on marginal covariance spectra are unavailable. The paper develops new techniques to bound small‑ball probabilities directly from the conditional structure, thereby characterizing the reproducing kernel Hilbert space (RKHS) of the Vecchia GP. The RKHS is shown to have the same smoothness order as the original Matérn RKHS, despite the non‑stationary and sparsely structured covariance of the Vecchia process.

On the statistical side, the authors consider the non‑parametric regression model Y_i = f_0(X_i) + ε_i with ε_i ∼ N(0,σ²). They place a Vecchia GP prior on f and study posterior contraction rates. By appropriately scaling the time and space hyperparameters (τ and s) and using the fixed‑size norming parent sets, they prove that the posterior contracts around the true Hölder‑α function at the minimax optimal rate n^{‑α/(2α+d)}. Crucially, this optimal rate is achieved without requiring the parent sets to grow with n; a constant‑size parent set suffices. To handle unknown smoothness, they embed τ and s in a hierarchical two‑level prior, showing that the resulting adaptive posterior also attains the minimax rate under mild conditions.

The paper supplies constructive algorithms for building norming parent sets. For data on regular grids, an explicit construction is given; for irregular designs, a greedy algorithm is proposed that seeks points satisfying the norming property. Computationally, the Vecchia GP retains O(n) time and memory complexity because each conditional depends on a bounded number of parents. The authors implement the core algorithms in C++ with an R interface and conduct synthetic experiments. Results demonstrate that the proposed fixed‑size norming parent sets outperform traditional nearest‑neighbor selections in terms of predictive accuracy while preserving linear computational cost.

In summary, this work establishes a solid probabilistic foundation for Vecchia GPs, derives their RKHS and small‑ball behavior, and proves that, when used as priors in Bayesian non‑parametric regression, they achieve optimal statistical performance with minimal computational overhead. The findings resolve longstanding questions about parent‑set size and selection, and open avenues for extending Vecchia approximations to broader kernel families, hierarchical models, and online settings.