Adaptive Resolution for Finite-Rank Gaussian Processes

Technical Summary: Adaptive Resolution for Finite-Rank Gaussian Processes

Problem Statement
Gaussian process (GP) regression is a standard tool for Bayesian nonparametric learning, yet exact implementation requires $O(n^2)$ memory and $O(n^3)$ computational complexity, rendering it infeasible for moderately large sample sizes $n$. While scalable approximations (sparse, low-rank, and distributed methods) have been developed, there is limited theoretical understanding of when these approximations retain the statistical guarantees (specifically posterior contraction rates) of their corresponding "parent" GP priors. A critical gap exists in finite-rank approximations: the resolution parameter $N$, which determines the dimension of the expansion, is often fixed. If $N$ is too small, the prior support may not contain the true function; if too large, computational costs become prohibitive. The paper addresses how to construct finite-rank GP priors that adaptively learn the resolution $N$ while preserving the optimal posterior contraction rates of the parent GP.

Methodology
The authors propose a general class of finite-rank GP priors, termed d-frGPloc, defined by a basis expansion with locally supported basis functions and dependent Gaussian coefficients:
$$f_N(x) = \sum_{j \in J_N} w_j \psi_j(x),$$
where the coefficient vector $w_N$ follows a dependent Gaussian prior $w_N | N, \theta \sim N(0, \Sigma_{N,\theta})$. The covariance matrix $\Sigma_{N,\theta}$ is induced from a parent GP $f \sim \text{GP}(0, K_\theta)$.

The paper studies two specific constructions within this framework:

1. SPDE-based Finite Element Approximation: Based on the stochastic partial differential equation (SPDE) representation of Matérn GPs (Lindgren et al., 2011). The coefficients are derived from a finite-element discretization of the SPDE, resulting in a sparse precision matrix structure. The bandwidth parameter $\kappa$ is treated as a deterministic sequence $\kappa_n$ dependent on sample size and smoothness, inherited from parent Matérn theory.
2. Lattice-based GP Interpolation (GPI): A regular-grid interpolation approach where coefficients on a grid are assigned a Gaussian prior induced from a parent squared-exponential (SE) GP. Here, both the bandwidth $\kappa$ and the resolution $N$ are assigned hierarchical priors.

Key Contributions
The primary theoretical contribution is the establishment of adaptive-resolution results. The authors demonstrate that by placing a suitable prior on the resolution parameter $N$, finite-rank approximations can inherit the posterior contraction rates of their parent GPs.

• Theoretical Results:

  • SPDE Construction (Theorem 3.3): Under the bandwidth scaling $\kappa_n$ derived from parent Matérn theory (Fang and Bhadra, 2025), placing a prior on $N$ that satisfies a mild lower-tail condition (e.g., polynomial or log-normal) ensures the finite-rank SPDE prior achieves the minimax posterior contraction rate $n^{-\alpha/(2\alpha+d)}$ (up to logarithmic factors) for true functions with smoothness $\alpha$. The proof relies on a spectral comparison between the continuous Neumann operator and its mass-lumped finite-element discretization, showing that the RKHS norms and small-ball probabilities match up to constants for sufficiently large $N$.
  • GPI Construction (Theorem 3.4): For the squared-exponential parent GP, placing hierarchical priors on both $\kappa$ and $N$ yields the minimax-optimal posterior contraction rate up to logarithmic factors. The proof combines the optimal prior concentration of the parent SE GP with interpolation error bounds, showing that the approximation error is negligible when $N$ is sufficiently large relative to $\kappa$.

• Computational Strategy:
  To address the computational challenge of varying dimension $N$, the authors develop a posterior sampler that avoids reversible-jump updates. The algorithm integrates out the coefficient vector $w_N$ to sample $(N, \kappa)$ from the marginal posterior, then samples $w_N$ conditionally from its Gaussian posterior.

  • For the GPI method, this exploits Toeplitz or block-Toeplitz structures.
  • For the SPDE method, the sparsity of the precision matrix allows for efficient sparse Cholesky factorization.
    The computational cost per iteration scales linearly with the sample size $n$ (due to local support) and cubically with the coefficient dimension $|J_N|$ (or linearly for 1D SPDE due to bandedness), making the approach scalable for large $n$.

Results
Numerical studies in one and two dimensions validate the theoretical findings:

• Performance Tracking: In 1D, the adaptive finite-rank procedures (both SPDE and GPI) closely track the predictive performance (measured by Average Mean Squared Error) of their corresponding exact parent GP baselines across different smoothness regimes.
• Adaptivity: In 2D, where exact parent GP computation is infeasible, the adaptive procedures outperform fixed-$N$ baselines. Fixed resolutions show high sensitivity to the chosen grid size (under-fitting if too small, over-fitting/computationally expensive if too large), whereas the adaptive method automatically adjusts the effective model complexity to the data without manual tuning.
• Posterior Behavior: The posterior distributions of $N$ and $\kappa$ adapt correctly: $N$ increases for rougher true functions and decreases for smoother ones, while $\kappa$ shifts according to the scale of the data.
• Scalability: The adaptive finite-rank methods remain computationally efficient as $n$ increases, avoiding the $O(n^3)$ bottleneck of exact GP inference while maintaining comparable accuracy.

Significance and Claims
The paper claims that finite-rank approximations need not sacrifice statistical optimality for computational tractability if the resolution is treated adaptively.

• Unification of Perspectives: The SPDE result connects the numerical approximation perspective (Sanz-Alonso and Yang, 2022) with statistical contraction theory (Fang and Bhadra, 2025), showing that finite-element discretizations can preserve parent GP contraction rates under appropriate priors on resolution.
• Adaptive Hierarchy: The work demonstrates that the finite-dimensional approximation dimension $|J_N| = (N+1)^d$ can be learned from data rather than fixed, preserving the parent-GP contraction rate while retaining the computational advantages of locally supported bases.
• Modesty: The authors note that while the adaptive-resolution hierarchy is theoretically attractive, it introduces computational challenges that are addressed via specific integration strategies. They acknowledge that extending these results to very high dimensions ($d \ge 5$) remains challenging due to the $O(N^{3d})$ cost, suggesting structured additive models as a potential future direction, and note that full adaptation over smoothness parameters for the SPDE construction remains an open problem.