Vecchia Parallel Partial Emulation (VPPE)

• Vecchia Parallel Partial Emulation (VPPE) is a scalable Gaussian process emulator that combines outputwise emulation with a Scaled Vecchia approximation to efficiently handle vector-valued outputs.
• It leverages parallel construction of emulators with shared input correlation structures, replacing dense Gaussian likelihood calculations with ordered nearest-neighbor factorizations.
• Empirical results show that VPPE achieves similar predictive accuracy to PPE while reducing computational time to a fraction of the original runtime.

Searching arXiv for papers on VPPE, PPE, and scaled Vecchia approximation. Vecchia Parallel Partial Emulation (VPPE) is a Gaussian-process emulator for computer models with vector-valued outputs and large numbers of training runs. It combines Parallel Partial Emulation (PPE), which models multidimensional simulator outputs through output-specific Gaussian processes with shared input-side correlation parameters, with the Scaled Vecchia approximation, which replaces dense Gaussian likelihood calculations by ordered nearest-neighbor conditional factorizations. In the reported computer experiments, VPPE yields comparable predictive accuracy to PPE at a fraction of the runtime (Seidman et al., 26 Aug 2025, Semochkina et al., 14 Feb 2025).

1. Conceptual lineage and definition

PPE was developed for simulators with matrix-valued or more general tensor-valued outputs. In the tensor-variate Gaussian-process framework, PPE is characterized by fitting “separate, shared-characteristic emulators at each output grid point,” with a common input-side regressor structure and common kernel hyperparameters, while emulators at different output locations are assumed conditionally independent a priori. In the matrix-output case, this places PPE alongside outer product emulators (OPE), but without the additional dependence structure that OPE assumes across output dimensions (Semochkina et al., 14 Feb 2025).

VPPE is the extension of PPE obtained by replacing the exact PPE marginal likelihood with a Scaled Vecchia approximation inside the RobustGaSP-style posterior used to estimate the shared range parameters. In this construction, PPE supplies the multidimensional-output model structure, Scaled Vecchia supplies a fast approximate Gaussian likelihood for large $n$, and RobustGaSP supplies the marginal-posterior estimation strategy for $\boldsymbol{\lambda}$, integrating out trend and variance parameters rather than maximizing a profile likelihood. The scalar-output analogue is termed VRGaSP (Seidman et al., 26 Aug 2025).

The designation “parallel partial” has a precise statistical meaning in this setting. “Parallel” refers to outputwise emulator construction with shared characteristics, and “partial” refers to the fact that PPE shares the input-side correlation structure but does not introduce a full cross-output covariance. VPPE preserves that statistical architecture and changes the computational regime in which it can be fitted (Semochkina et al., 14 Feb 2025).

2. Statistical structure inherited from PPE

For scalar output, the emulator is written as

$$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$

with mean

$\mu(x)=\mathbf h(x)\boldsymbol\beta,$

and covariance

$C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$

For a training design $x^D=\{x_1,\dots,x_n\}$, the response vector satisfies

$$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$

where $h$ is the basis matrix and $\mathbf R$ is the input correlation matrix (Seidman et al., 26 Aug 2025).

PPE extends this to vector-valued output $Y^\mathcal D\in\mathbb R^{n\times k}$ by assigning each output coordinate its own regression coefficients and variance, while all coordinates share the same input correlation matrix $\boldsymbol{\lambda}$0. Under common $\boldsymbol{\lambda}$1, the predictive distribution at $\boldsymbol{\lambda}$2 is

$\boldsymbol{\lambda}$3

with predictive mean

$\boldsymbol{\lambda}$4

A key PPE identity is the weighted-sum representation

$\boldsymbol{\lambda}$5

which is emphasized as preserving physical features such as conservation inherited from the simulator (Seidman et al., 26 Aug 2025).

This weighted-sum predictor is central to the practical appeal of PPE and therefore of VPPE. The approximation introduced by VPPE acts primarily in fitting the shared range parameters; prediction retains the PPE form once $\boldsymbol{\lambda}$6 has been estimated (Seidman et al., 26 Aug 2025).

3. Scaled Vecchia approximation inside the PPE posterior

The Vecchia approximation starts from the exact factorization

$\boldsymbol{\lambda}$7

and replaces the full conditioning history by a small conditioning set $\boldsymbol{\lambda}$8: $\boldsymbol{\lambda}$9 If the conditioning set size is $$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$0, the total cost is $$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$1, much smaller than $$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$2 when $$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$3 (Seidman et al., 26 Aug 2025).

VPPE uses the Scaled Vecchia approximation. Inputs are rescaled by initial range values,

$$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$4

and Scaled Vecchia performs maximin ordering on the scaled inputs together with nearest-neighbor conditioning in that scaled space. The motivation follows earlier scaled Vecchia emulation work, which interprets $$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$5 as input relevance and uses the transformed input geometry to make neighborhood selection more faithful in anisotropic computer experiments (Seidman et al., 26 Aug 2025, Katzfuss et al., 2020).

For scalar output, integrating the Vecchia likelihood over $$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$6 yields

$$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$7

For $$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$8-dimensional output, VPPE reuses the same $$y(\cdot)\sim \mathcal{GP}\big(\mu(\cdot), C(\cdot,\cdot)\big),$$9 and $\mu(x)=\mathbf h(x)\boldsymbol\beta,$0 across output coordinates and computes a separate $\mu(x)=\mathbf h(x)\boldsymbol\beta,$1 for each output dimension: $\mu(x)=\mathbf h(x)\boldsymbol\beta,$2 When $\mu(x)=\mathbf h(x)\boldsymbol\beta,$3, this approximation recovers the exact RobustGaSP marginal likelihood (Seidman et al., 26 Aug 2025).

In this sense, VPPE is not a new output-side stochastic specification. It is a replacement of the dense $\mu(x)=\mathbf h(x)\boldsymbol\beta,$4 Gaussian marginal calculations inside PPE by an ordered nearest-neighbor conditional approximation that is shared across the $\mu(x)=\mathbf h(x)\boldsymbol\beta,$5 output coordinates (Seidman et al., 26 Aug 2025).

4. Computational interpretation and parallel implementations

The computational bottleneck of PPE is the repeated use of the $\mu(x)=\mathbf h(x)\boldsymbol\beta,$6 correlation matrix $\mu(x)=\mathbf h(x)\boldsymbol\beta,$7 when estimating $\mu(x)=\mathbf h(x)\boldsymbol\beta,$8. The reported PPE fitting cost is

$\mu(x)=\mathbf h(x)\boldsymbol\beta,$9

where $C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$0 is the number of optimization iterations or posterior evaluations. VPPE replaces this with

$C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$1

by substituting many small $C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$2 matrix operations for a single dense $C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$3 factorization. The paper also notes that the conditional Gaussian terms $C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$4 can be computed in parallel (Seidman et al., 26 Aug 2025).

This computational view places VPPE within a broader Vecchia implementation literature. GPU-accelerated Vecchia likelihood evaluation for geospatial Gaussian processes has been shown to reduce the time to solution compared to ExaGeoStat by up to $C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$5, $C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$6, and $C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$7 on 32GB GV100, 80GB A100, and 80GB H100 GPUs, respectively, and to accommodate up to $C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$8M locations on a single NVIDIA GPU while maintaining application accuracy (Pan et al., 2024). A separate GPU study emphasizes that Vecchia approximation “can be calculated with embarrassingly parallel algorithms” and reports that the GpGpU implementation achieves faster runtimes and better predictive accuracy on datasets including a large satellite dataset with $C(\cdot,\cdot)=\sigma^2 c(\cdot,\cdot).$9 points (James et al., 2024). At the distributed end, the Scaled Block Vecchia (SBV) algorithm uses MPI and MAGMA to achieve near-linear scalability on up to 64 A100 and GH200 GPUs, handles 320M points, and is presented as the first distributed implementation of any Vecchia-based GP variant (Pan et al., 16 Apr 2025).

These results are not VPPE itself, but they define the computational ecosystem in which VPPE belongs: outputwise emulation over shared hyperparameters, combined with local conditional Gaussian calculations that are naturally batched, GPU-amenable, and, in distributed settings, compatible with MPI-style reductions.

5. Empirical behavior in computer-model emulation

The reported VPPE experiments cover a synthetic GP simulator, a hydrology model based on Richards’ equation, and a volcanic flow model. Across all three, the dominant pattern is near-PPE predictive accuracy with substantially lower fit time (Seidman et al., 26 Aug 2025).

Study	Setup	Reported result
Synthetic GP simulator	$x^D=\{x_1,\dots,x_n\}$0, $x^D=\{x_1,\dots,x_n\}$1, $x^D=\{x_1,\dots,x_n\}$2	PPE took 6100 s; VPPE took 176 s; minimum relative RMSEs were essentially identical at about 0.0199
Hydrology model	1696 training runs, 300 test runs	PPE fit time 488 s with RMSE $x^D=\{x_1,\dots,x_n\}$3; VPPE fit times 19.4–384 s with RMSE $x^D=\{x_1,\dots,x_n\}$4 to $x^D=\{x_1,\dots,x_n\}$5 for $x^D=\{x_1,\dots,x_n\}$6 to $x^D=\{x_1,\dots,x_n\}$7
Volcanic flow model	3863 training runs, $x^D=\{x_1,\dots,x_n\}$8, $x^D=\{x_1,\dots,x_n\}$9	PPE fit time 23.66 hours with RMSE 0.0952; VPPE fit time 3.34 hours with RMSE 0.0959

In the synthetic study with $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$0, the practical gain appears once $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$1 is a few hundred. For $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$2 and $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$3, VPPE matched PPE to four decimal places in relative RMSE while using less than $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$4 of the runtime in one of the reported comparisons (Seidman et al., 26 Aug 2025).

The hydrology study provides a more detailed picture of approximation behavior. Across 20 repetitions with 300 training runs and 1696 test runs, ANOVA showed that the method explained only $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$5 of RMSE variation; median RMSE was 0.00237 for PPE and 0.00236 for VPPE. In the larger 1696/300 split, smaller $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$6 slightly improved RMSE, and the paper interprets this as a case where using only the closest shapes may improve predictive fit. The same study also reports that nearest-neighbor prediction with $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$7 around 150–175 out of 300 considered was best, and that using too many neighbors could slightly degrade prediction (Seidman et al., 26 Aug 2025).

The volcanic flow model is the most extreme application. After restricting the TITAN2D output to a region of interest, the output dimension was $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$8, and the training design contained 3864 runs before holding one out for testing. In that case, randomly sampling $$y^D \mid \boldsymbol\beta,\sigma^2,\boldsymbol\lambda \sim MVN\big(h\boldsymbol\beta,\sigma^2\mathbf R\big),$$9 of output dimensions when estimating $h$0 did not noticeably alter final estimated range parameters in practice. The paper also reports that local prediction with $h$1 improved RMSE to 0.0949 and reduced prediction runtime to less than 5 seconds, compared with over 17 minutes using all training data (Seidman et al., 26 Aug 2025).

6. Scope, related Vecchia developments, and limitations

VPPE inherits PPE’s statistical assumptions. Output coordinates are modeled as independent conditional on shared input correlation parameters, dependence across outputs is handled indirectly via common $h$2 and the weighted-sum predictor, and one common range parameter vector $h$3 must be suitable across all outputs. This differs sharply from OPE, which assumes additional dependence structure across output dimensions and can provide more accurate predictions when advantage can be taken of correlation in the output dimensions (Semochkina et al., 14 Feb 2025).

Within the broader Vecchia literature, VPPE sits on top of several methodological themes. Scaled Vecchia factorization implies a sparse inverse-Cholesky structure and supports unbiased mini-batch gradients in large-scale GP regression (Cao et al., 2022). Prediction-focused Vecchia work shows that different ordering and conditioning choices have strong effects on uncertainty quantification and computational cost, and distinguishes between scalable joint prediction and highly parallel local marginal prediction (Katzfuss et al., 2018). A later unification result proves that “partial Cholesky + Vecchia = Vecchia,” meaning that a low-rank global component plus a sparse local residual correction can be interpreted exactly as a Vecchia approximation with an augmented sparsity pattern (Kaminetz et al., 5 Mar 2026). These developments clarify that VPPE is one instance of a larger family of sparse conditional GP emulators.

The limitations reported for VPPE are correspondingly specific. Approximation quality depends on the ordering, the neighbor size $h$4, and the prediction neighbor size $h$5. If $h$6 is too small, long-range dependence may be underrepresented. If outputs differ too much in their effective correlation over inputs, the shared $h$7 assumption may be restrictive. When $h$8 is enormous, subsampling output dimensions during estimation may be used, as in the volcanic example. The paper also notes that if one needs exact global conservation in prediction, local prediction with small $h$9 may weaken that property (Seidman et al., 26 Aug 2025).

Taken together, these features define VPPE as a computational extension of PPE rather than a replacement for its multi-output modeling assumptions. Its contribution is to make the PPE marginal posterior tractable for larger training sets by inserting a Scaled Vecchia approximation at the point where dense Gaussian-process algebra would otherwise dominate runtime (Seidman et al., 26 Aug 2025).