---
title: Vecchia Parallel Partial Emulation (VPPE)
url: https://www.emergentmind.com/topics/vecchia-parallel-partial-emulation-vppe
type: topic
---

# Vecchia Parallel Partial Emulation (VPPE)

Searching arXiv for recent 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 [2508.19144][2502.10319].

## 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 [2502.10319].

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 [2508.19144].

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 [2502.10319].

## 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 [2508.19144].

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 \(\mathbf R(\boldsymbol\lambda)\). Under common \(\boldsymbol\lambda\), the predictive distribution at \(x^*\) is
\[
\mathbf y(x^*)\mid x^D,Y^\mathcal D,\boldsymbol\lambda \sim t_{n-q}\!\left( \hat{\mathbf y}(x^*), \mathrm{diag}(\hat{\boldsymbol\sigma}^2)c^{**} \right),
\]
with predictive mean
\[
\hat{\mathbf y}^\top(x^*) = \mathbf h(x^*)\hat{\boldsymbol\beta} + \mathbf r^\top(x^*)\mathbf R^{-1}\big(Y^\mathcal D-h\hat{\boldsymbol\beta}\big).
\]
A key PPE identity is the weighted-sum representation
\[
\hat{\mathbf y}(x^*) = \boldsymbol\omega(x^*)Y^\mathcal D,
\]
which is emphasized as preserving physical features such as conservation inherited from the simulator [2508.19144].

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 \(\hat{\boldsymbol\lambda}\) has been estimated [2508.19144].

## 3. Scaled Vecchia approximation inside the PPE posterior

The Vecchia approximation starts from the exact factorization
\[
p(\mathbf y)=\prod_{i=1}^n p(y_i\mid y_1,\dots,y_{i-1}),
\]
and replaces the full conditioning history by a small conditioning set \(b(i)\subset\{1,\dots,i-1\}\):
\[
\hat p(\mathbf y)=\prod_{i=1}^n p(y_i\mid \mathbf y_{b(i)}).
\]
If the conditioning set size is \(m\), the total cost is \(\mathcal O(nm^3)\), much smaller than \(\mathcal O(n^3)\) when \(m\ll n\) [2508.19144].

VPPE uses the Scaled Vecchia approximation. Inputs are rescaled by initial range values,
\[
\tilde x = (x_1/\lambda_1,\dots,x_p/\lambda_p)^\top,
\]
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 \(1/\lambda_l\) as input relevance and uses the transformed input geometry to make neighborhood selection more faithful in anisotropic computer experiments [2508.19144][2005.00386].

For scalar output, integrating the Vecchia likelihood over \((\boldsymbol\beta,\sigma^2)\) yields
\[
\mathcal L_m(\boldsymbol\lambda\mid x^D,y^D) \propto \left(\prod_{i=1}^n \omega_{im}^{-1/2}\right) |\tilde{\Sigma}|^{-1/2} (\tilde S^2)^{-(n-q)/2}.
\]
For \(k\)-dimensional output, VPPE reuses the same \(\omega_{im}\) and \(\tilde\Sigma\) across output coordinates and computes a separate \(\tilde S_l^2\) for each output dimension:
\[
\mathcal L_m(\boldsymbol\lambda\mid x^D,Y^D) \propto \left( \left(\prod_{i=1}^n \omega_{im}^{-1/2}\right) |\tilde\Sigma|^{-1/2} \right)^k \prod_{l=1}^k (\tilde S_l^2)^{-(n-q)/2}.
\]
When \(m=n-1\), this approximation recovers the exact RobustGaSP marginal likelihood [2508.19144].

In this sense, VPPE is not a new output-side stochastic specification. It is a replacement of the dense \(n\times n\) Gaussian marginal calculations inside PPE by an ordered nearest-neighbor conditional approximation that is shared across the \(k\) output coordinates [2508.19144].

## 4. Computational interpretation and parallel implementations

The computational bottleneck of PPE is the repeated use of the \(n\times n\) correlation matrix \(\mathbf R\) when estimating \(\boldsymbol\lambda\). The reported PPE fitting cost is
\[
\mathcal O(t n^3) + \mathcal O(t n^2 k),
\]
where \(t\) is the number of optimization iterations or posterior evaluations. VPPE replaces this with
\[
\mathcal O(nm^3)+\mathcal O(m^2k),
\]
by substituting many small \(m\times m\) matrix operations for a single dense \(n\times n\) factorization. The paper also notes that the conditional Gaussian terms \(p(y_i\mid y_{b(i)})\) can be computed in parallel [2508.19144].

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 \(700\times\), \(833\times\), and \(1380\times\) on 32GB GV100, 80GB A100, and 80GB H100 GPUs, respectively, and to accommodate up to \(1\)M locations on a single NVIDIA GPU while maintaining application accuracy [2403.07412]. 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 \(n>10^6\) points [2407.02740]. 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 [2504.12004].

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 [2508.19144].

| Study | Setup | Reported result |
|---|---|---|
| Synthetic GP simulator | \(n=4000\), \(k=100\), \(m=30\) | 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 \(6.91\times10^{-4}\); VPPE fit times 19.4–384 s with RMSE \(6.37\times10^{-4}\) to \(6.86\times10^{-4}\) for \(m=5\) to \(100\) |
| Volcanic flow model | 3863 training runs, \(k=130{,}262\), \(m=20\) | 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 \(m=30\), the practical gain appears once \(n\) is a few hundred. For \(k=100\) and \(n=4000\), VPPE matched PPE to four decimal places in relative RMSE while using less than \(3\%\) of the runtime in one of the reported comparisons [2508.19144].

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 \(0.03\%\) of RMSE variation; median RMSE was 0.00237 for PPE and 0.00236 for VPPE. In the larger 1696/300 split, smaller \(m\) 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 \(m_{\text{pred}}\) around 150–175 out of 300 considered was best, and that using too many neighbors could slightly degrade prediction [2508.19144].

The volcanic flow model is the most extreme application. After restricting the TITAN2D output to a region of interest, the output dimension was \(k=130{,}262\), and the training design contained 3864 runs before holding one out for testing. In that case, randomly sampling \(10\%\) of output dimensions when estimating \(\boldsymbol\lambda\) did not noticeably alter final estimated range parameters in practice. The paper also reports that local prediction with \(m_{\text{pred}}=200\) improved RMSE to 0.0949 and reduced prediction runtime to less than 5 seconds, compared with over 17 minutes using all training data [2508.19144].

## 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 \(\boldsymbol\lambda\) and the weighted-sum predictor, and one common range parameter vector \(\boldsymbol\lambda\) 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 [2502.10319].

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 [2202.12981]. 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 [1805.03309]. 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 [2603.05709]. 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 \(m\), and the prediction neighbor size \(m_{\text{pred}}\). If \(m\) is too small, long-range dependence may be underrepresented. If outputs differ too much in their effective correlation over inputs, the shared \(\boldsymbol\lambda\) assumption may be restrictive. When \(k\) 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 \(m_{\text{pred}}\) may weaken that property [2508.19144].

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 [2508.19144].

Source: https://www.emergentmind.com/topics/vecchia-parallel-partial-emulation-vppe