Probabilistic Numerics: Bayesian Computation

• Probabilistic numerics is a field that treats numerical computation as statistical inference by modeling discretization error as epistemic uncertainty.
• It employs Bayesian frameworks, like Gaussian processes, to quantify error in methods such as integration, ODE/PDE solvers, and linear algebra.
• This approach enables adaptive refinement, coherent uncertainty propagation through computation pipelines, and improved decision-making.

Probabilistic numerics is a research field that treats numerical computation as a problem of statistical inference, systematically modeling discretization and algorithmic error as epistemic uncertainty. This paradigm posits a prior distribution on the latent mathematical object of interest (e.g., a function, field, trajectory, or system state), uses information obtained during numerical computation as observational data, and computes a posterior distribution over the quantity of interest. The resulting probabilistic numerical method (PNM) yields both a point estimate (typically the posterior mean) and a full measure of uncertainty, which can be propagated through computational pipelines or explicitly used in decision-making, adaptive refinement, and uncertainty quantification.

1. Conceptual and Mathematical Foundations

Emphasizing epistemic uncertainty, probabilistic numerics models the error incurred by finite computation as a probability distribution over the output, in sharp contrast to traditional deterministic schemes that provide only point estimates or conservative worst-case bounds (Oates et al., 2019). The general framework is rooted in Bayesian inference:

• Prior: A distribution (often Gaussian process in function spaces) encodes structural assumptions (smoothness, regularity).
• Likelihood: The information derived from the numerical computation (e.g., pointwise evaluation, matrix-vector products, PDE residuals) is treated as a (usually linear) observation model.
• Posterior: Conditioning the prior on this information yields the posterior distribution, which quantifies the remaining epistemic uncertainty.

The pushforward of the posterior through the quantity-of-interest map $T$ yields a measure over outcomes; for example, the distribution over $\int_a^b f(x)\,dx$ when inferring integrals (Briol et al., 2015), or over the solution $u(x)$ to a PDE when using GP-based solvers (Bilionis, 2016, Cockayne et al., 2016).

This perspective also underpins the concept of Bayesian PNMs—methods that are, up to disintegration, equivalent to Bayesian inversion of an observation model. General conditions for existence and convergence for Bayesian PNMs, even in non-Gaussian and nonlinear settings, are established using measure-theoretic tools (e.g., regular conditional probabilities) (Cockayne et al., 2017).

2. Historical Development and Key Advances

Probabilistic numerics originated in foundational ideas from Sul’din (Wiener process priors for quadrature; 1959), Larkin (fully Bayesian error modelling for integration in the 1970s), and was formally advanced by O’Hagan, Skilling, and subsequent groups who linked average-case analysis, information-based complexity, and Bayesian decision theory (Oates et al., 2019). Parallel developments formalized the equivalence between Bayes-optimal numerical estimators and average-case optimal rules under specified priors (Oates et al., 2019).

Recent decades have witnessed explosive growth, driven by advances in Gaussian processes, Bayesian inverse problems, kernel quadrature, probabilistic ODE/PDE solvers, and the formalization of uncertainty quantification (UQ) in computational pipelines. This has yielded a technology stack encompassing linear algebra, integration, ODE/PDE simulation, optimization, and inverse problems (Hennig et al., 2015, Wenger et al., 2021).

3. Core Methodologies Across Domains

3.1 Bayesian Quadrature and Probabilistic Integration

Bayesian quadrature places a GP prior over the integrand $f$, computes a posterior GP on $f$ given evaluations, and pushes this posterior through the integration operator to obtain a Gaussian law for the integral (Briol et al., 2015):

• Posterior mean: a kernel quadrature rule
• Posterior variance: a principled measure of numerical uncertainty

Rates of posterior contraction coincide with optimal estimation rates in corresponding reproducing kernel Hilbert spaces (RKHS), and the probabilistic error shrinks under classical sample designs (IID, MCMC, QMC, higher-order nets). The approach enables coherent calibration and propagation of error in statistical computation pipelines (Briol et al., 2015).

3.2 Probabilistic Solvers for ODEs and PDEs

In ODEs, Gauss-Markov priors (e.g., integrated Wiener process, Ornstein-Uhlenbeck process) yield filter-based solvers where each time-step assimilates the ODE residual via Bayesian filtering, maintaining a trajectory distribution (Bosch et al., 2023, Fay et al., 2023). Probabilistic exponential integrators encode stiff linear dynamics exactly in the prior, yielding L-stability and calibrated covariance in stiff systems (Bosch et al., 2023).

For PDEs, a widely used construction assigns a GP prior to the unknown solution field $u$, builds a likelihood via collocation (enforcing the residuals of the PDE and boundary conditions at a finite set of points), and conditions on this information (Bilionis, 2016, Cockayne et al., 2016, Thakur et al., 13 Aug 2025). The resulting posterior GP quantifies epistemic uncertainty induced by discretization:

• Posterior mean: classical collocation estimator (e.g., meshless symmetric collocation)
• Posterior covariance: sharp, pointwise error bars for discretization error
• Adaptive refinement: highest posterior variance regions guide sample placement, yielding mesh-free and $h$-adaptive solvers (Thakur et al., 13 Aug 2025)

Matrix decompositions, stochastic dual descent, and clustering-based active learning enable scalability to high-dimensional PDEs beyond the cubic bottleneck of classical GP regression (Thakur et al., 13 Aug 2025).

3.3 Black-Box and Composite PNMs

Black-box PNMs (BBPN) construct posteriors over the limiting value of a convergent sequence of numerical approximations (e.g., Richardson extrapolation), generalizing classical methods to non-linear or adaptive solvers where only the final output at each resolution is accessible (Teymur et al., 2021). This approach yields not only acceleration in convergence (order boosting) but also calibrated UQ for broad numerical tasks, including nonlinear eigenproblems and PDE simulation.

PNMs can be composed modularly: for instance, solutions of probabilistic PDE solvers can be embedded within Bayesian inverse problem pipelines, automatically propagating solver error through inference (Cockayne et al., 2016, Cockayne et al., 2017). The composition properties are formally captured via measure-theoretic pushforward and disintegration.

4. Optimality, Experimental Design, and Adaptive Methods

Classical decision-theoretic optimality criteria (Bayesian average-case analysis) yield point-optimal numerical estimators but can be insufficient for uncertainty quantification, as they do not penalize posterior spread (Oates et al., 2019). For UQ-driven design, a distinct Bayesian probabilistic numerics (BPN) criterion targets minimization not just of estimator error but of posterior dispersion around the true state. In general, BPN-optimal information (e.g., node locations, discretization points) differs from classical experimental designs unless loss is strictly quadratic in the quantity of interest.

Sequential experimentation—choosing data adaptively by maximizing posterior variance reduction (active learning)—is enabled in probabilistic numerics by the closed-form expression for posterior risk (Bilionis, 2016, Thakur et al., 13 Aug 2025). Adaptive batch sizes and integration error-aware strategies are directly driven by the GP posterior variance in Bayesian quadrature for active learning (Adachi et al., 2023), providing dynamic cost-precision trade-offs unavailable to fixed heuristics.

5. Uncertainty Quantification, Propagation, and Applications

PNMs offer rigorous, statistical UQ for arbitrary quantities derived from numerical computation. Posterior variances are interpretable as epistemic uncertainty due to limited discretization, finite iteration, or incomplete data. This property is critical in several contexts:

• Sequential pipelines: Uncertainty can be coherently propagated through multiple computational stages (e.g., PDE solve $\rightarrow$ inverse problem $\rightarrow$ decision-making), preserving overall coverage and credibility (Cockayne et al., 2017, Wenger et al., 2021).
• Design of experiments: Bayesian design criteria can be directly applied to minimize the propagation of discretization error in statistical inferences, e.g., in PDE-constrained inverse problems (Cockayne et al., 2016).
• Monte Carlo estimators: For random-time-step integrators, the mean-square convergence rate of the Monte Carlo estimator is decoupled from the number of samples at appropriate noise tuning, making UQ feasible with minimal sampling overhead (Abdulle et al., 2018).
• Fast and reliable calibration: Posterior credible intervals automatically contract as grid density or data increase, ensuring frequentist coverage in a variety of scientific contexts (thermodynamic integration, complex computer models, and high-dimensional Bayesian computation) (Briol et al., 2015).

6. Computational Realizations and Software Infrastructure

State-of-the-art implementation of PNMs is exemplified by the ProbNum library (Wenger et al., 2021), which provides modular, composable probabilistic solvers for linear algebra, ODEs, Bayesian quadrature, and more:

• Architecture: Modular separation of prior models, information operators, belief updates, and adaptive policies.
• API: High-level interfaces compatible with standard scientific Python APIs.
• Uncertainty propagation: Native chaining and composition for composite pipelines.
• Customization: User-definable priors, kernels, and model parameters, with empirical Bayes/hyperparameter optimization routines.

Many PNM algorithms recover the convergence properties of classical solvers and are implemented with numerically stable filtering/backward-smoothing for trajectory inference and covariance control.

7. Open Challenges and Research Frontiers

While PNM methodology is now well-developed for quadrature, linear problems, ODE/PDE discretization with scalable matrix solvers, several difficult research questions remain (Oates et al., 2019):

• Prior calibration and automatic model selection: Determining appropriate kernel hyperparameters, error structure, and regularization in high dimensions, with robust empirical Bayes strategies.
• Scalability and computational efficiency: Reducing cubic scaling of GP solvers, especially for high-dimensional PDEs, via stochastic inference, sparsity, and surrogate models (Thakur et al., 13 Aug 2025).
• Adaptive and sequential design: Formal integration of fully Bayesian experimental design for adaptive mesh refinement, node placement, and information acquisition.
• Industrial and scientific use-cases: Widespread adoption in areas with dominant numerical uncertainty, such as high-dimensional Bayesian inverse problems, optimal experimental design, and complex systems simulation pipelines.
• Software ecosystems and education: Development and dissemination of robust, tested PNMs, integration with probabilistic programming, and creation of curricular resources for broader scientific communities.

The probabilistic numerics paradigm unifies mathematical-statistical inference and numerical computation, endowing scientific computing with rigorous, compositional, and adaptive uncertainty quantification that can be carried throughout computational workflows (Hennig et al., 2015, Oates et al., 2019, Wenger et al., 2021).