Parallel Partial GP Model

• Parallel Partial GP models are scalable GP constructions that avoid full N×N kernel inversion by using techniques like truncated Mercer decompositions, inducing points, and block partitions.
• They leverage parallel computing across CPUs, GPUs, and distributed networks to reduce computational cost while maintaining reliable prediction and uncertainty estimation.
• Variants such as PPE and VPPE optimize multi-output GP regression by balancing shared input structure with partially decoupled outputs to achieve significant computational speedups.

Searching arXiv for recent and foundational papers on "parallel partial Gaussian process" and closely related GP formulations. In the cited literature, the expression parallel partial Gaussian process model does not denote a single universally standardized GP prior. Rather, it refers to a class of scalable Gaussian-process constructions in which exact full-covariance inference is replaced by a partial representation—such as a truncated kernel expansion, a support-set approximation, a partitioned mixture of experts, a reduced multi-output representation, or a parallel partial emulator—and the resulting computation is organized for parallel execution on CPUs, GPUs, or distributed machines. Across these variants, the common objective is to avoid direct formation or inversion of the full dense $N\times N$ kernel matrix while retaining GP-based prediction and uncertainty quantification (Carminati, 2024, Dai et al., 2014, Chen et al., 2014, Semochkina et al., 14 Feb 2025).

1. Terminological scope and conceptual identity

The term is best understood as an umbrella description for several related GP strategies rather than as a single canonical model family. In one line of work, it refers to a scalable sparse-variational framework whose objective is only partially factorized over datapoints, so that large data-dependent terms can be computed in parallel while a small inducing-point system remains globally coupled (Dai et al., 2014). In another, it appears as Parallel Partial Emulation (PPE), a formal multi-output emulator in which output locations are conditionally independent a priori but share common input-space regressors and correlation structure (Semochkina et al., 14 Feb 2025). In a third, it describes a decomposed or reduced GP in which the kernel is replaced by a finite Mercer expansion and posterior computation is accelerated with CUDA (Carminati, 2024).

This diversity of usage matters because “partial” is not used uniformly. In Mercer-based models, the partiality lies in a low-rank kernel decomposition. In support-set and inducing-point methods, it lies in keeping only a compact global representation while distributing datapoint-local computations. In PPE, it lies in partial decoupling across output coordinates while preserving shared dependence on the input domain. In hyperspectral imaging, a “parallel GP” may operate on the full spatial domain but only on a reduced set of components obtained by unmixing, which is again partial in the representational sense (Kalinin et al., 2020).

The literature also identifies nearby but distinct usages. “Physics-assisted Gaussian processes” concern PDE-informed losses rather than a parallel partial GP architecture, and the word “partial” there refers to partial differential equations, not to partial covariance structure (Zhang et al., 2022). Likewise, GP-based time-parallel ODE solvers and GP surrogate Bayesian inference with batch evaluation are parallel GP methods, but not parallel partial GP models in the sense of partitioned, decomposed, or partially coupled GP regression (Pentland et al., 2022, Järvenpää et al., 2019).

2. Mathematical motifs underlying partiality

The standard GP regression bottleneck is the $\mathcal{O}(N^3)$ solve involving the training covariance matrix. For training data $(X,\mathbf y)$ and test data $X_*$, the predictive posterior takes the familiar form

$$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$

with $K=k(X,X)$ (Carminati, 2024). Parallel partial GP constructions alter this computation by changing the covariance representation, the output structure, or both.

A first motif is kernel truncation via Mercer decomposition: $$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$ In matrix form,

$$k(X,X') \approx \mathbf\Phi_{(X)}\,\Lambda\,\mathbf\Phi_{(X')}^T.$$

This replaces the full covariance by $N\times n$ feature matrices and an $n\times n$ diagonal eigenvalue matrix. After Woodbury reduction, the expensive inversion is transferred to

$\mathcal{O}(N^3)$0

whose inversion costs $\mathcal{O}(N^3)$1 instead of $\mathcal{O}(N^3)$2 (Carminati, 2024).

A second motif is conditional block structure. In support-set approximations such as PITC and PIC, or in partitioned expert models, data are divided into subsets, and the full covariance is approximated by a support-conditioned block structure or a block-diagonal covariance conditional on a partition (Chen et al., 2014, Zhang et al., 2017). The computational effect is that one replaces one large dense inversion by multiple smaller inversions, local summaries, and a small amount of global synchronization.

A third motif is partial coupling of outputs. In the EELS/STEM formulation, multiple spatial loading maps are reconstructed jointly with covariance

$\mathcal{O}(N^3)$3

so outputs share a common spatial kernel while retaining component-specific residual variance and effective noise levels (Kalinin et al., 2020). In PPE, the matrix- or tensor-valued emulator instead assumes output-wise independent residuals through a diagonal coregionalization matrix,

$\mathcal{O}(N^3)$4

while retaining shared regressors $\mathcal{O}(N^3)$5 and shared input correlation $\mathcal{O}(N^3)$6 (Semochkina et al., 14 Feb 2025).

These motifs can coexist. A model may be partial because it is low-rank in the kernel, partial because it is only partially coupled across outputs, or partial because it distributes only data-local terms while retaining a small global latent system.

3. Principal model families in the literature

The major formulations can be organized by the locus of approximation and the locus of parallelism.

Family	Partial structure	Representative paper
Mercer-decomposed GP	Finite eigenfunction expansion $\mathcal{O}(N^3)$7	(Carminati, 2024)
Support-set / low-rank distributed GP	Local blocks plus global support set or low-rank factor	(Chen et al., 2014)
Sparse variational inducing-point framework	Datapoint-local sufficient statistics plus small inducing system	(Dai et al., 2014)
Partitioned mixture-of-experts GP	Block-diagonal covariance conditional on sampled partitions	(Zhang et al., 2017)
Reduced multi-output parallel GP	Joint GP on loading maps after spectral decomposition	(Kalinin et al., 2020)
Parallel Partial Emulator	Shared input regressors and kernel, diagonal output covariance	(Semochkina et al., 14 Feb 2025, Seidman et al., 26 Aug 2025)

In the Mercer-decomposed formulation, the squared exponential kernel is expanded in terms of eigenfunctions and eigenvalues. In one dimension,

$\mathcal{O}(N^3)$8

with eigenfunctions

$\mathcal{O}(N^3)$9

and eigenvalues

$(X,\mathbf y)$0

The multidimensional ARD extension has product eigenfunctions and product eigenvalues, but the basis size grows as $(X,\mathbf y)$1, which is the central computational limitation motivating GPU acceleration (Carminati, 2024).

In the support-set framework, pPITC and pPIC divide the training set into blocks $(X,\mathbf y)$2, compute local summaries such as

$(X,\mathbf y)$3

aggregate them into global summaries

$(X,\mathbf y)$4

and then distribute prediction over test partitions (Chen et al., 2014). The paper proves predictive equivalence between pPITC and centralized PITC, between pPIC and centralized PIC, and between pICF-based GP and centralized incomplete-Cholesky GP.

In the sparse variational inducing-point framework, a lower bound is written in terms of additive datapoint statistics such as

$(X,\mathbf y)$5

so the heavy work becomes embarrassingly parallel over datapoints and only a small $(X,\mathbf y)$6 matrix involving the inducing variables remains non-distributable (Dai et al., 2014).

In partitioned expert models, latent cluster assignments define $(X,\mathbf y)$7 independent GP experts conditional on a partition, yielding a block-diagonal covariance. Importance sampling over many partitions then restores flexibility by averaging over multiple block structures rather than committing to a single one (Zhang et al., 2017).

In PPE, the emulator at each output coordinate $(X,\mathbf y)$8 is

$(X,\mathbf y)$9

with independent residuals across output locations and a shared input kernel $X_*$0. The tensor-variate GP framework shows PPE to be a special case of a broader Kronecker-structured multi-output emulator, alongside the outer product emulator (OPE) (Semochkina et al., 14 Feb 2025).

4. Parallelization mechanisms and computational layout

Parallelism enters these models at different algorithmic levels. In the CUDA implementation of the Mercer-decomposed GP, the expensive stages are the construction of $X_*$1 and $X_*$2, the dense products $X_*$3, and the posterior mean and covariance terms. The implementation uses cuBLAS for GEMM/GEMV operations and cuSOLVER for LU factorization and solution of the reduced $X_*$4 linear system associated with $X_*$5; the CPU baseline is implemented in C++ with Eigen, which uses OpenMP for multithreading in GEMM-like routines (Carminati, 2024).

In distributed support-set methods, parallelism is coarse-grained. Each machine computes its local summaries from its own data block, a master aggregates these summaries, and prediction is then redistributed. The communication burden depends on the algorithm: pPITC and pPIC communicate summaries whose size depends on the support set $X_*$6, whereas pICF-based GP communicates quantities that depend on the number of prediction points $X_*$7 (Chen et al., 2014).

In the sparse variational framework, the computation has a clear map-reduce flavor. Workers independently accumulate local contributions to $X_*$8, $X_*$9, $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$0, and gradient terms, and only the small inducing-point algebra needs global coordination. GPU acceleration targets repeated evaluation of $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$1 and $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$2, which can dominate more than $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$3 of inference time in large unsupervised models (Dai et al., 2014).

In mixture-of-experts inference, parallelism occurs both across importance samples and across experts inside each partition sample. Conditional on a sampled partition, each expert fits its own local GP, with block inversions of size roughly $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$4. The only synchronization is the normalization of importance weights and the final averaging of predictions (Zhang et al., 2017).

The computational gains reported in the literature are therefore heterogeneous. For the CUDA Mercer model, with $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$5, varying retained eigenvalues $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$6, and $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$7, the GPU provides no major advantage at $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$8, begins to dominate as $$\boldsymbol\mu_* = m(X_*) + K_* (K+\sigma^2 I)^{-1}(\mathbf y - m(X)), \qquad \Sigma_* = K_{**} - K_*(K+\sigma^2 I)^{-1}K_*^T,$$9 and $K=k(X,X)$0 grow, and on the strongest reported setup—RTX 2080 Super versus Ryzen 5 3600—is about $K=k(X,X)$1 faster at $K=k(X,X)$2, with runtime about $K=k(X,X)$3 versus $K=k(X,X)$4 (Carminati, 2024). In the distributed GP framework, empirical evaluation on a cluster of 20 nodes shows pPITC, pPIC, and pICF-based GP to be significantly more time-efficient and scalable than centralized counterparts and exact/full GP, with pPITC and pPIC often $K=k(X,X)$5–$K=k(X,X)$6 orders of magnitude faster than full GP (Chen et al., 2014).

5. Multi-output reconstruction and parallel partial emulation

A distinct branch of the literature uses the term in connection with structured multi-output GP regression. In hyperspectral EELS/STEM analysis, the workflow is to decompose the data cube

$K=k(X,X)$7

using NMF or a related unmixing method, then reconstruct the spatial loading maps $K=k(X,X)$8 jointly with a GP that shares a common spatial covariance structure across components (Kalinin et al., 2020). The GP therefore acts neither on the full three-dimensional cube nor on each map independently; it operates on the full spatial domain but only on a reduced set of components.

This model is “parallel” because several maps are reconstructed simultaneously, and “partial” because the coupling is not fully unrestricted. Information is shared via the common spatial kernel structure, while output covariance $K=k(X,X)$9 and per-component effective noise levels $$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$0 allow variability in relative noise magnitude and morphology (Kalinin et al., 2020). Hyperparameters are trained jointly by maximizing the marginal log likelihood, and the evolution of kernel length scale and noise during optimization is used diagnostically to assess whether a physically meaningful scale has been learned.

The experimental EELS example uses a LaAlO$$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$1/SrTiO$$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$2 interface dataset of size $$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$3, decomposed with NMF using $$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$4 components. Independent GP reconstruction recovers atomic-scale contrast in the first four components but fails for components 5 and 6; the parallel GP reconstruction recovers atomic-scale features in these weaker components as well, indicating that higher-order spatial information can be reconstructed and spatially localized by borrowing structure from better-constrained components (Kalinin et al., 2020).

The Parallel Partial Emulator is a more formal construction. In the tensor-variate GP framework, PPE is the special case where output locations are conditionally independent a priori, but all output coordinates share the same regressors in the input domain and the same input correlation function. Its mean tensor has the form

$$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$5

and, for noninformative priors on regression coefficients, the generalized least-squares estimator of the PPE mean does not depend on the output covariance $$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$6 (Semochkina et al., 14 Feb 2025). This makes PPE parsimonious and often computationally attractive when output grids are large and explicit output-space dependence is weak, unknown, or undesirable.

A later extension, VPPE, incorporates the Scaled Vecchia approximation into PPE so that larger training sets become tractable. The original PPE has complexity

$$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$7

whereas VPPE reduces this to

$$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$8

for neighbor size $$k(\mathbf x,\mathbf x')=\sum_{i=1}^{\infty}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x'), \qquad k(\mathbf x,\mathbf x') \approx \sum_{i=1}^{n}\lambda_i \phi_i(\mathbf x)\phi_i(\mathbf x').$$9 (Seidman et al., 26 Aug 2025). VPPE preserves the same partial-emulation structure—shared correlation parameters across outputs, output-specific variances and trends, and Student-$$k(X,X') \approx \mathbf\Phi_{(X)}\,\Lambda\,\mathbf\Phi_{(X')}^T.$$0 predictive distributions—while replacing the exact covariance likelihood by a scaled Vecchia likelihood approximation.

6. Limitations, controversies, and nearby formulations

A recurrent limitation is that the partial structure can itself become expensive or statistically restrictive. In the Mercer-decomposed model, the multidimensional ARD basis is a tensor product with $$k(X,X') \approx \mathbf\Phi_{(X)}\,\Lambda\,\mathbf\Phi_{(X')}^T.$$1 and $$k(X,X') \approx \mathbf\Phi_{(X)}\,\Lambda\,\mathbf\Phi_{(X')}^T.$$2, so the number of terms grows exponentially in the input dimension $$k(X,X') \approx \mathbf\Phi_{(X)}\,\Lambda\,\mathbf\Phi_{(X')}^T.$$3; the reduced model is therefore still computationally heavy in high dimensions (Carminati, 2024). In the EELS formulation, performance depends strongly on the quality of the initial unmixing decomposition and on correct kernel constraints; forcing a length scale in $$k(X,X') \approx \mathbf\Phi_{(X)}\,\Lambda\,\mathbf\Phi_{(X')}^T.$$4 pixels when the true structure is atomic-scale can cause reconstruction failure (Kalinin et al., 2020).

Partition-based methods introduce another trade-off. A single fixed partition implies zero inter-block correlation and can be overconfident if cross-block uncertainty is ignored. The importance-sampled mixture-of-experts model addresses this by averaging over many partitions, but that introduces dependence on the proposal distribution and the possibility of weight degeneracy if the proposal is poor (Zhang et al., 2017).

Support-set and low-rank distributed methods are also approximation-sensitive. Larger support sets $$k(X,X') \approx \mathbf\Phi_{(X)}\,\Lambda\,\mathbf\Phi_{(X')}^T.$$5 improve pPITC and pPIC accuracy at higher cost; larger ICF rank $$k(X,X') \approx \mathbf\Phi_{(X)}\,\Lambda\,\mathbf\Phi_{(X')}^T.$$6 improves pICF-based GP accuracy but increases computation, and if $$k(X,X') \approx \mathbf\Phi_{(X)}\,\Lambda\,\mathbf\Phi_{(X')}^T.$$7 is too small, predictive variance quality can degrade and even become non-PSD in practice (Chen et al., 2014). In PPE, the lack of borrowing across output locations can be a weakness when the simulator output has meaningful spatial or temporal geometry. The tensor-variate comparison with OPE shows that OPE is preferable when output locations have strong dependence that can be modeled and exploited, whereas PPE is preferable when output dependence is weak, hard to justify, or the output dimension is so large that independent residuals are advantageous (Semochkina et al., 14 Feb 2025).

The literature also warns against terminological conflation. PAGP is a physics-assisted GP framework in which PDE residuals are added to GP training losses; it is not a parallel GP framework in the sense of parallel partial modeling (Zhang et al., 2022). GParareal is a time-parallel ODE solver that models fine-minus-coarse corrections with a GP emulator, but its partiality lies in time-slice decomposition of an IVP rather than in a partial GP covariance architecture (Pentland et al., 2022). GP-based surrogate Bayesian inference with noisy likelihood evaluations uses a hierarchical GP and batch-sequential design so that expensive simulations can be run in parallel, but it does not introduce a special partial GP construction (Järvenpää et al., 2019).

Taken together, these papers suggest that the enduring significance of the parallel partial Gaussian process idea is not a single invariant model form, but a recurring design principle: retain the probabilistic semantics of GP regression while replacing global dense covariance computation by a smaller, structured, or partially decoupled representation that admits effective parallel execution (Dai et al., 2014, Carminati, 2024, Semochkina et al., 14 Feb 2025).