Scaled Vecchia Approximation

• Scaled Vecchia approximation is a Gaussian process method that redefines neighbor selection by rescaling input coordinates with covariance parameters to capture anisotropic effects.
• It replaces full conditioning with small, geometry-aware predecessor sets, reducing the likelihood evaluation cost from O(n³) to near-linear time.
• Advanced implementations leverage batched GPU operations and distributed frameworks to efficiently handle large-scale problems in emulation and high-dimensional regression.

Scaled Vecchia approximation is a family of Gaussian-process (GP) approximations built on Vecchia’s ordered conditional factorization, in which scalability is achieved by replacing full conditioning with small predecessor sets and, in its most characteristic form, by carrying out ordering and neighbor selection in an anisotropically transformed input space. In the computer-emulation and high-dimensional regression literature, the defining modification is to rescale each input coordinate according to estimated covariance range or relevance parameters before applying maximin ordering and nearest-neighbor conditioning, so that the approximation follows the covariance geometry rather than raw Euclidean distance (Katzfuss et al., 2020, Cao et al., 2022). In GPU-oriented geospatial work, the same label has also been used for a high-performance implementation of the classical Vecchia approximation based on batched matrix computations at million-location scale (Pan et al., 2024).

1. Conditional factorization and sparse inverse-Cholesky structure

For a Gaussian process observed at locations $S_n=\{x_1,\dots,x_n\}$, the exact joint density factorizes as

$$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$

Vecchia approximation replaces each full conditioning set $h(i)$ by a much smaller subset $g(i)\subset h(i)$ of size at most $m$,

$$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$

or, in equivalent notation used elsewhere,

$$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$

This yields a proper Gaussian likelihood with a sparse inverse-Cholesky structure,

$\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$

where $U_{n,m}$ is sparse upper triangular (Kang et al., 2024, Katzfuss et al., 2020).

The computational consequence is the principal reason for the method’s adoption. Each conditional involves only a small $(m+1)\times(m+1)$ covariance block, so the full likelihood evaluation costs $$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$0, or, with grouping strategies, about $$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$1. Exact GP inference, by contrast, requires $$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$2 time and $$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$3 memory. The approximation therefore converts dense cubic-cost likelihood evaluation into near-linear or quasilinear computation when $$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$4 (Kang et al., 2024).

Scaled Vecchia preserves this conditional factorization. The modification is not a new covariance model; rather, it changes how the ordering and local conditioning neighborhoods are defined. In that sense, scaled Vecchia is best viewed as a geometry-aware specialization of the classical Vecchia construction (Pan et al., 2024, Katzfuss et al., 2020).

2. Anisotropic scaling and covariance-aware geometry

The distinctive feature of scaled Vecchia is that ordering and nearest-neighbor conditioning are carried out in a scaled input space rather than in the raw domain. For anisotropic covariance parameters $$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$5, the scaled coordinates are

$$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$6

and the corresponding distance is

$$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$7

The Vecchia approximation is then applied in the transformed space $$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$8, with maximin ordering and nearest-neighbor conditioning computed in that geometry (Katzfuss et al., 2020).

This rescaling is motivated by relevance or factor sparsity. If an input dimension $$p_n(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{h(i)}\right), \qquad h(i)=\{1,\dots,i-1\}.$$9 has a small range parameter $h(i)$0, then changes in that input strongly affect the response; if $h(i)$1 is large, the dimension is weakly relevant, and in the extreme $h(i)$2 it becomes effectively inactive. The quantity $h(i)$3 is therefore used as a measure of relevance (Katzfuss et al., 2020). In the variable-selection setting of VGPR, the same idea is written through the anisotropically scaled distance

$h(i)$4

so that the approximation itself adapts to the current relevance vector $h(i)$5 (Cao et al., 2022).

The practical effect is that the neighbor graph follows the covariance geometry implied by the active covariates. Under maximin ordering, each new point is placed as far as possible from previously ordered points, while each conditioning set is chosen as the nearest previously ordered neighbors in the scaled metric. This arrangement is intended to improve the screening effect, reduce the neighbor size $h(i)$6 needed for a given accuracy, and avoid poor neighborhoods caused by unscaled Euclidean geometry in anisotropic or high-dimensional problems (Kang et al., 2024, Katzfuss et al., 2020).

The same principle recurs in later formulations. VPPE applies maximin ordering and nearest-neighbor selection after the transformation

$h(i)$7

and the distributed Scaled Block Vecchia method rescales each coordinate as

$h(i)$8

before block construction and neighbor search (Seidman et al., 26 Aug 2025, Pan et al., 16 Apr 2025).

3. Parameter estimation, likelihood optimization, and prediction

Scaled Vecchia methods typically estimate covariance parameters by maximizing the Vecchia log-likelihood rather than the exact GP likelihood. In large-scale computer emulation, Fisher scoring is used after profiling out the mean parameters. With $h(i)$9 denoting the Vecchia log-likelihood, the iteration is

$g(i)\subset h(i)$0

where $g(i)\subset h(i)$1 is the gradient and $g(i)\subset h(i)$2 is the expected negative Hessian. Because the scaling parameters alter the transformed inputs, the ordering and neighborhood structure should in principle change during optimization; the practical scheme is to ignore the derivative of the ordering and neighborhoods with respect to $g(i)\subset h(i)$3 when computing gradients, but to update the ordering and neighbor structure occasionally at iterations $g(i)\subset h(i)$4 (Katzfuss et al., 2020).

In high-dimensional regression and variable selection, the scaled Vecchia likelihood is embedded in a penalized objective,

$g(i)\subset h(i)$5

with the likelihood written as

$g(i)\subset h(i)$6

VGPR traverses a regularization path from strong to weak penalization, adds candidate covariates according to the squared-relevance gradient, and removes irrelevant covariates through quadratic constrained coordinate descent. It also introduces Vecchia-based mini-batch subsampling,

$g(i)\subset h(i)$7

with unbiased gradient estimator

$g(i)\subset h(i)$8

relative to the scaled Vecchia objective (Cao et al., 2022).

Prediction is handled by applying Vecchia to the joint vector of training and prediction variables. In the emulation setting, training points are ordered before prediction points and neighbor sets $g(i)\subset h(i)$9 are chosen among previously ordered variables, yielding

$m$0

This produces joint Gaussian predictive distributions with sparse inverse-Cholesky structure and supports fast predictive means, variances, and joint simulation, with prediction cost $m$1 for fixed $m$2 (Katzfuss et al., 2020).

In Robust GaSP-based formulations such as VPPE, the same conditional structure is combined with marginalization over $m$3 and $m$4, producing a Vecchia marginal likelihood for the range parameters and a Vecchia marginal posterior for $m$5 (Seidman et al., 26 Aug 2025). This indicates that scaled Vecchia is compatible with both direct likelihood maximization and integrated marginal-posterior optimization.

4. Computational realizations: batched GPUs, blocks, and distributed systems

A major development of scaled Vecchia is the conversion of its many-small-matrix workload into hardware-efficient batched linear algebra. In GPU-accelerated geospatial inference, each location $m$6 produces a neighbor covariance block $m$7, a cross-covariance vector $m$8, and the neighboring observations $m$9. These are processed with batched covariance generation, batched Cholesky factorization batchPOTRF, batched triangular solves batchTRSV, and batched dot products. The implementation uses the KBLAS library and a strided layout so that batched matrices are stored contiguously in GPU memory (Pan et al., 2024).

The resulting performance gains are substantial. Compared with the exact MLE implementation in ExaGeoStat, the GPU Vecchia implementation reports speedups of up to about $$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$0 on NVIDIA GV100 32GB, about $$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$1 on NVIDIA A100 80GB, and about $$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$2 on NVIDIA H100 80GB. It also accommodates up to 1 million locations on a single 80GB A100 or H100 GPU while maintaining the necessary application accuracy (Pan et al., 2024).

Block formulations modify the conditional decomposition itself. Block Vecchia replaces the product of univariate conditionals with a product of multivariate block-level conditionals,

$$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$3

with blocks formed using K-means and processed by MAGMA batched kernels. Classical Vecchia is recovered when the block size equals 1. The method reports approximately $$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$4 speedup over classical Vecchia in some settings and the ability to handle problem sizes about $$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$5 larger (Pan et al., 2024).

Scaled Block Vecchia combines anisotropic input scaling with block Vecchia and distributed GPU execution. Its pipeline consists of anisotropic scaling, MPI-based redistribution, Random Anchor Clustering, filtered $$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$6-nearest-neighbor search, and batched GPU linear algebra using MAGMA. The method is described as the first distributed implementation of any Vecchia-based GP variant, achieves near-linear scalability on up to 64 A100 and GH200 GPUs, handles 320M points, and is demonstrated on a 50M point simulation from a respiratory disease model (Pan et al., 16 Apr 2025).

These realizations clarify that “scaled” may refer either to covariance-aware input scaling or to extreme-scale implementation. Both uses retain the same sparse conditional backbone.

5. Asymptotic theory, screening effect, and inferential equivalence

The main rigorous asymptotic theory currently available concerns Vecchia approximation itself. Under fixed-domain asymptotics, for a large class of covariance functions with boundary conditioning, consistency and asymptotic normality of maximum exact-likelihood estimators imply the same properties for maximum Vecchia-likelihood estimators, provided the conditioning-set size grows polylogarithmically,

$$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$7

Under the same scaling, the approximation becomes asymptotically exact in Kullback–Leibler divergence,

$$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$8

and exact GP prediction is asymptotically matched by Vecchia prediction,

$$\hat p_{n,m}(y;\theta)=\prod_{i=1}^n p_\theta\!\left(y_i \mid y_{g(i)}\right),$$9

The paper identifies this as a rigorous manifestation of the screening effect: only a polylogarithmic number of nearest neighbors is enough for asymptotically accurate prediction (Kang et al., 2024).

The same work shows that if

$$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$0

then for some $$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$1,

$$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$2

The inferential message is that quasilinear-complexity Vecchia inference is asymptotically equivalent to cubic-complexity exact GP inference in that setting (Kang et al., 2024).

A separate fixed-domain asymptotic analysis for Matérn models establishes analogous results for the microergodic parameter $$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$3 under Vecchia’s approximation. On $$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$4 with equally spaced points and neighbor sets of size $$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$5, $$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$6, the Vecchia-based MLE preserves the same asymptotic normal limit as the full-likelihood estimator for the microergodic combination (Zhang et al., 2021).

These results do not directly constitute a theorem for every anisotropically scaled variant. They are, however, the principal theoretical justification for replacing dense GP likelihoods by sparse conditional Vecchia likelihoods. The boundary-conditioned results are formally proved only in that setting, while extension to ordinary Matérn GPs without boundary conditioning is argued theoretically and numerically rather than fully proved (Kang et al., 2024).

6. Empirical behavior and application domains

In large computer-model emulation, scaled Vecchia has been reported to outperform standard Vecchia, simplified low-rank alternatives, and local approximate GP methods. In a 10-dimensional Matérn example with two relevant dimensions, SVecchia with $$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$7 was more accurate than Vecchia or LowRank with $$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$8, implying about a 1000-fold reduction in cost for the same accuracy level because the cost scales roughly like $$p(y) \approx \prod_{i=1}^n p(y_i \mid y_{c(i)}), \qquad |c(i)|=\min(m,i-1).$$9. On the satellite-drag simulator, SVecchia met a 1% relative error benchmark, was the most accurate method for all six species, and the best-performing laGP method could take up to about 200 times longer while still being less accurate (Katzfuss et al., 2020).

In high-dimensional GP regression and variable selection, the scaled Vecchia approximation is reported to be substantially more accurate than unscaled Vecchia and to outperform FIC, FITC, and PIC in KL-divergence experiments. The scaled version produced additional improvements of several orders of magnitude in approximation error compared with unscaled Vecchia when the neighborhoods were rebuilt using the relevant covariate scaling. The resulting VGPR method is described as handling $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$0 responses and $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$1 covariates, finishing in under forty minutes on a workstation in the application study (Cao et al., 2022).

For vector-valued computer-model outputs, VPPE embeds scaled Vecchia inside Parallel Partial Emulation. In the synthetic study with $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$2 and $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$3, the maximum VPPE runtime was 176 seconds compared with 6100 seconds for PPE, while the minimum relative RMSE was essentially the same for both methods, about $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$4. In the volcanic flow application with $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$5 runs and $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$6 output pixels, VPPE fit in 3.34 hours versus 23.66 hours for PPE, with RMSE $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$7 versus $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$8 (Seidman et al., 26 Aug 2025).

The approximation has also been used in large geospatial and spatio-temporal applications. In Argo temperature modeling, a constant-mean Matérn GP with ARD over 5D or 7D spatio-temporal inputs is made scalable by Vecchia approximation, using $\hat K_{n,m}=(U_{n,m}U_{n,m}^\top)^{-1},$9 for fitting and $U_{n,m}$0 for prediction, and the reported predictive gains include over 13% RMSE reduction at 300 dbar and 1500 dbar in some settings (Liu et al., 20 Oct 2025). In low-cost urban sensor monitoring, six ordering schemes, three distance metrics, and two sensor-conditioning rules were compared; the min-max distance algorithm and random ordering achieved, on average, 102% better results in log Kullback–Leibler divergence compared with four other orderings (Idir et al., 27 Nov 2025).

These studies support a consistent empirical pattern: when anisotropy or covariate relevance is pronounced, scaling the input geometry materially improves the conditional approximation, and when the implementation is matched to batched or distributed linear algebra, the resulting methods remain feasible at very large $U_{n,m}$1.

7. Generalizations, hybrids, and learned Vecchia kernels

Several later developments retain the scaled or local-conditional logic of Vecchia while changing the covariance representation. NeuVec models the deterministic quantities arising in the Vecchia factorization,

$U_{n,m}$2

and learns the kriging coefficients $U_{n,m}$3 and conditional standard deviations $U_{n,m}$4 directly. These quantities uniquely determine the covariance structure under the Vecchia model through

$U_{n,m}$5

The method is built around permutation-invariant and permutation-preserving architectures for conditioning sets and is claimed to achieve linear complexity in $U_{n,m}$6 and also linear complexity in $U_{n,m}$7 (Cao et al., 6 May 2026).

Hybrid models combine Vecchia with other scalable GP mechanisms. Vecchia-inducing-points full-scale approximations decompose the process into a global inducing-point component and a local residual approximated by Vecchia, interpolating between FITC and pure Vecchia as special cases. The method uses a correlation-based neighbor metric for the residual process and extends to non-Gaussian likelihoods through Laplace-based iterative methods that reduce computational costs by several orders of magnitudes compared to Cholesky-based computations in some settings (Gyger et al., 7 Jul 2025). Earlier Vecchia-Laplace methods already combined Vecchia’s sparse Gaussian algebra with Laplace approximations for generalized Gaussian processes, covering Bernoulli, Poisson, Gamma, and log-Gaussian Cox process models (Zilber et al., 2019).

A plausible implication is that scaled Vecchia approximation should be understood less as a single algorithm than as a design principle: preserve the sparse directed conditional factorization, but adapt the geometry, optimization, or hardware realization to the structure of the problem at hand. Within that principle, the recurring technical themes are covariance-aware local neighborhoods, sparse inverse-Cholesky representations, and computational regimes in which accuracy improves faster than the cost of increasing the neighborhood size.