Prob-GParareal: Probabilistic Parallel-in-Time Solver

• Prob-GParareal is a probabilistic parallel-in-time solver that uses Gaussian process posterior laws to propagate uncertainty across differential equation trajectories.
• It replaces the deterministic correction in GParareal with a random sampling approach, yielding empirical distributions over trajectories and supporting probabilistic initial conditions.
• The scalable Prob-nnGParareal variant leverages a nearest-neighbors GP strategy to reduce computational cost while maintaining accuracy in high-dimensional systems.

Prob-GParareal is a probabilistic numerical Parallel-in-Time solver for differential equations that extends GParareal by replacing the deterministic Gaussian-process correction with a random correction and propagating that uncertainty across time. In place of a single corrected trajectory, it returns empirical distributions over trajectories, supports probabilistic initial conditions, and retains compatibility with classical numerical solvers. In the formulation introduced in 2025, Prob-GParareal is presented as the first probabilistic extension of GParareal, while Prob-nnGParareal adapts the nearest-neighbors GP strategy to improve scalability on larger problems (Gattiglio et al., 4 Sep 2025, Pentland et al., 2022, Gattiglio et al., 2024).

1. Algorithmic lineage and defining construction

Classical Parareal splits the time interval into subintervals, advances a cheap coarse solver $G$ sequentially, and corrects it with an expensive fine solver $F$ computed in parallel. In the notation used for Prob-GParareal, its update is

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$

with convergence checked by

$$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$

GParareal retains the same coarse/fine structure but replaces the explicit Parareal correction by a GP posterior mean for the correction function $f_c=(F-G)$, yielding

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$

Prob-GParareal keeps this learned-correction viewpoint but makes the correction random: $${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$ with

$$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$

This makes the iteration a Markov process in time, driven by a GP-induced conditional law rather than a deterministic residual (Gattiglio et al., 4 Sep 2025).

The defining shift from GParareal is therefore not the use of Gaussian processes per se, but the use of the full GP posterior law rather than only the posterior mean. GParareal already models the correction function $x\mapsto(\mathcal F-\mathcal G)(x)$ and can use all accumulated correction data, or even legacy data from prior runs, to accelerate convergence (Pentland et al., 2022). Prob-GParareal adds uncertainty quantification to that architecture by sampling the correction and propagating the resulting distribution forward in time (Gattiglio et al., 4 Sep 2025).

2. Gaussian-process model for the correction

The correction function is modeled coordinatewise. For each component $s=1,\ldots,d$,

$F$0

The method uses $F$1 independent scalar-output GPs, one per coordinate, rather than a full multi-output GP, for computational reasons. The training set at iteration $F$2 is

$F$3

In Prob-GParareal, these data are collected using sample summaries, specifically the sample mean $F$4 of the propagated samples: $F$5 For a test point $F$6, the posterior mean and variance of the $F$7-th GP are

$F$8

$F$9

Here $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$0 contains the GP inputs, $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$1 the outputs, and $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$2 is a nugget/jitter regularization parameter (Gattiglio et al., 4 Sep 2025).

The principal kernel in the experiments is the Gaussian or RBF kernel,

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$3

and the paper also discusses the Matérn family,

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$4

The hyperparameters $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$5 are tuned independently for each coordinate by maximizing the marginal log-likelihood,

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$6

Within the method’s interpretation, the GP posterior covariance is epistemic information about the learned discrepancy model, and Prob-GParareal uses it directly rather than discarding it as GParareal does (Gattiglio et al., 4 Sep 2025).

3. Uncertainty propagation, empirical forecasts, and stopping rules

Conditioned on the previous state, the probabilistic update is Gaussian: $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$7 Unconditionally, the law is a continuous Gaussian mixture,

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$8

which is generally intractable. The implementation therefore uses ancestral sampling: for each interval $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \left(\mathscr{F}-\mathscr{G}\right)({u}_{i-1,k-1}),$$9 and iteration $$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$0, it draws $$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$1 sample trajectories $$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$2 by sampling from the conditional Gaussian along the time-marching recursion (Gattiglio et al., 4 Sep 2025).

This sample-based propagation changes both the solver output and the stopping logic. Instead of comparing point values between iterations, Prob-GParareal can compare empirical distributions using a Wasserstein distance,

$$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$3

The algorithm can also stop by an external exit condition such as a runtime or variance budget. This allows early stopping while retaining a probabilistic forecast whose spread reflects remaining numerical uncertainty (Gattiglio et al., 4 Sep 2025).

Forecast quality is assessed with multivariate probabilistic scores. The energy score is

$$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$4

and the variogram score is

$$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$5

with $$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$6 and $$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$7. The paper also reports MAD and MSE, and uses

$$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$8

to connect probabilistic error to squared Wasserstein distance (Gattiglio et al., 4 Sep 2025).

Probabilistic initial conditions are handled directly. For deterministic data, the initial law is a Dirac mass,

$$\max_{1\leq i\leq N-1}\|{u}_{i,k}-{u}_{i,k-1}\|_\infty < \epsilon.$$9

whereas random initial conditions are represented by draws from $f_c=(F-G)$0; the experiments study Gaussian perturbations

$f_c=(F-G)$1

The reported behavior is that $f_c=(F-G)$2 sets a lower bound on the solution uncertainty, larger $f_c=(F-G)$3 increases the spread of the forecast, and too-large initial uncertainty can slow convergence or cause nonconvergence, especially for chaotic systems (Gattiglio et al., 4 Sep 2025).

4. Error analysis, variance control, and computational complexity

The theoretical analysis assumes that the fine solver is exact,

$f_c=(F-G)$4

that the coarse solver has local truncation error of order $f_c=(F-G)$5,

$f_c=(F-G)$6

that the correction function is Lipschitz,

$f_c=(F-G)$7

that $f_c=(F-G)$8 is Lipschitz,

$f_c=(F-G)$9

and that the GP posterior variance decays with fill distance in one of three regimes: differentiability, Sobolev norm-equivalence, or infinite smoothness. The local fill distance is

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$0

Under these assumptions, the main convergence theorem states

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$1

where $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$2 depends on the regime under consideration and $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$3 depends on fill distance and GP approximation quantities (Gattiglio et al., 4 Sep 2025).

For the infinitely smooth case, the coefficient is given as

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$4

A corollary states that if $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$5, then

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$6

The analysis also provides a variance recursion for the maximum coordinate variance

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$7

namely

$${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$8

and a bound comparing the mean of the probabilistic iterate $${u}_{i,k} = \mathscr{G}({u}_{i-1,k}) + \widehat{f}_{\rm GPara}({u}_{i-1,k}), \qquad \widehat{f}^{(s)}_{\rm GPara}({u}_{i-1,k})=\boldsymbol{\mu}^{(s)}_{\mathcal{D}_k}({u}_{i-1,k}).$$9 with the deterministic GParareal iterate $${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$0,

$${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$1

These results formalize the dependence of forecast quality on GP approximation quality and on the regularity and Lipschitz properties of the coarse/fine discrepancy (Gattiglio et al., 4 Sep 2025).

The model cost reflects both GP training and sampling. For the correction-model component at iteration $${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$2,

$${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$3

and the total cost is approximated by

$${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$4

The distance-computation term $${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$5 is typically $${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$6 for Wasserstein-$${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$7, but $${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$8 for $${U}_{i, k}=G\left({U}_{i-1, k}\right)+\boldsymbol{Z}_{i, k}, \quad i=1,\ldots,N, \quad k \geq 1,$$9 under Gaussian assumptions. The paper notes that with sufficient parallelism the extra cost from the $$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$0 trajectories can be fully parallelized, making the total wall-clock cost comparable to deterministic GParareal (Gattiglio et al., 4 Sep 2025).

5. Numerical behavior and the Prob-nnGParareal variant

The reported ODE benchmarks are FitzHugh–Nagumo, Rössler, Hopf bifurcation, double pendulum, and Lorenz. The summary table in the source reports good calibration and sharpness across all systems, with small metrics on converged runs. FitzHugh–Nagumo shows MAD $$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$1 and very small VS/ES/MSE; Hopf bifurcation and double pendulum show small but nonzero errors; Lorenz and Rössler show larger ES and growing uncertainty over time, while the forecast remains meaningful. A normalized bias comparison against Parareal and GParareal indicates that the mean accuracy of Prob-GParareal is comparable to deterministic solvers, but chaotic systems display larger errors, as expected (Gattiglio et al., 4 Sep 2025).

The temporal behavior of uncertainty is system-dependent. The reported standard deviation grows steadily for Lorenz and Rössler, grows more mildly for the double pendulum, and stabilizes for FitzHugh–Nagumo and Hopf. For the Rössler system, long-horizon propagation of uncertain initial conditions produces a 95% interval from Prob-GParareal that matches the empirical interval obtained from many fine-solver simulations. The early-stopping experiments show that stopping at about $$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$2 of convergence often preserves good forecast scores for FitzHugh–Nagumo, Lorenz, and the double pendulum, with reduced runtime (Gattiglio et al., 4 Sep 2025).

Prob-nnGParareal addresses the cubic GP cost by adapting the nearest-neighbors strategy developed for nnGParareal. Instead of conditioning on all $$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$3 points, each local GP is trained on only the $$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$4 nearest neighbors of the prediction point, with $$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$5, and a kd-tree is used for efficient neighbor search. In the probabilistic formulation, this changes the local model cost to $$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$6 per interval for training and $$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$7 per draw for sampling (Gattiglio et al., 4 Sep 2025, Gattiglio et al., 2024).

The PDE demonstration for Prob-nnGParareal uses the viscous Burgers equation

$$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$8

with $$\boldsymbol{Z}_{i,k}|{U}_{i-1, k}={u}_{i-1,k} \sim \mathcal{N}_d(\boldsymbol{\mu}_{\mathcal{D}_k}({u}_{i-1, k}), \Sigma_{\mathcal{D}_k}({u}_{i-1,k})).$$9, $x\mapsto(\mathcal F-\mathcal G)(x)$0, finite-difference spatial discretization, $x\mapsto(\mathcal F-\mathcal G)(x)$1, $x\mapsto(\mathcal F-\mathcal G)(x)$2, $x\mapsto(\mathcal F-\mathcal G)(x)$3, $x\mapsto(\mathcal F-\mathcal G)(x)$4, fine solver RK8, and coarse solver RK1. The reported outcome is that Prob-nnGParareal converges in fewer iterations than nnGParareal, runtime drops to about 33 minutes, nnGParareal requires roughly 1h 3m, GParareal about 1h 39m, and the fine solver alone about 13h 5m; bias and MSE remain comparable to the deterministic variants, while uncertainty is slightly larger than in Prob-GParareal because of the local GP approximation (Gattiglio et al., 4 Sep 2025).

6. Position within the broader Parareal literature

Prob-GParareal belongs to the family of GP-enhanced Parareal methods, but its role is distinct. GParareal already replaced the explicit Parareal defect by GP emulation of the correction function and showed that one can converge in fewer iterations than Parareal, sometimes succeed where Parareal fails, and incorporate archives of legacy solutions from prior runs (Pentland et al., 2022). nnGParareal then made that idea scalable by replacing the global GP with a nearest-neighbors GP whose model complexity is reduced from cubic to log-linear in the sample size, improving performance on higher-dimensional systems and larger processor counts (Gattiglio et al., 2024). Prob-GParareal takes the next step by treating the learned correction as a random variable and propagating the resulting uncertainty across the time decomposition (Gattiglio et al., 4 Sep 2025).

A related but different probabilistic precursor is the micro-macro Parareal scheme for stochastic differential equations with time-scale separation, where the fine model is a weighted Euler–Maruyama Monte Carlo ensemble, the coarse model is a one-dimensional Fokker–Planck equation, and the macroscopic correction acts on probability densities through update operators such as the quantile iterator. That method addresses probability distributions directly, but its state representation and correction mechanism are distributional rather than GP-emulative (Lelièvre et al., 2019). Prob-GParareal instead keeps the coarse/fine Parareal architecture of GParareal and injects uncertainty through a conditional Gaussian correction law (Gattiglio et al., 4 Sep 2025).

The surrounding deterministic Parareal literature clarifies what Prob-GParareal inherits and what it does not yet resolve. Robust convergence for parabolic problems with backward Euler coarse propagators and arbitrarily high-order fine propagators has been analyzed under strong-stability conditions such as $x\mapsto(\mathcal F-\mathcal G)(x)$5, with linear convergence factors near $x\mapsto(\mathcal F-\mathcal G)(x)$6 for suitable coarse/fine ratios (Yang et al., 2021). Parareal has also been adapted to differential-algebraic systems arising in eddy current simulation by separating differential and algebraic components and exploiting the automatic projection effect of implicit Euler (Schöps et al., 2017), and to coupled elliptic-parabolic saddle-point systems such as quasi-static Biot poroelasticity through monolithic, fixed-stress, and multirate fixed-stress coarse propagators with explicit contractivity conditions (Jimenez-Ciga et al., 21 Jan 2026). This suggests that extending Prob-GParareal beyond the ODE and Burgers settings studied so far would likely require the same kind of structure-aware treatment of stiffness, DAE consistency, and coarse-propagator contractivity, but the published Prob-GParareal formulation itself is centered on uncertainty-aware GP correction rather than on those deterministic stability mechanisms (Gattiglio et al., 4 Sep 2025).