PaperlayerView exact revision on arXiv

Source-linked AI summary

Bayesian Probabilistic Numerical Methods

Jon Cockayne, Chris Oates, Tim Sullivan, Mark Girolami

arXiv:1702.03673v2stat.MEmath.NAmath.STstat.CO

TL;DR

Probabilistic numerics lacked a general statistical foundation and a clear account of what makes a method Bayesian. This paper formulates Bayesian PNMs as Bayesian inverse problems, establishes well-definedness conditions, develops approximations and composition theory, and applies the framework to nonlinear, non-Gaussian, and industrial settings.

  • Problem

    Probabilistic numerical methods lacked common theoretical foundations, including general conditions for well-definedness and a clear definition of Bayesian status.

  • Method

    The paper casts Bayesian PNMs as Bayesian inverse problems, uses disintegration to define their outputs, develops numerical approximations, and analyzes composed PNM pipelines.

  • Results

    The framework establishes that Bayesian PNMs are well-defined under mild conditions, with essentially unique disintegrations, and extends to nonlinear, non-Gaussian models and composed pipelines.

  • Takeaways & Limitations

    Bayesian PNMs provide an extensible framework for meaningful uncertainty propagation across numerical-method pipelines when constituent methods and priors are coherent.

  • Takeaways & Limitations

    The paper leaves sufficient conditions for posterior concentration on the data-generating quantity of interest, including identifiability and concentration rate, for further work.

Abstract

from arXiv · show

The emergent field of probabilistic numerics has thus far lacked clear statistical principals. This paper establishes Bayesian probabilistic numerical methods as those which can be cast as solutions to certain inverse problems within the Bayesian framework. This allows us to establish general conditions under which Bayesian probabilistic numerical methods are well-defined, encompassing both non-linear and non-Gaussian models. For general computation, a numerical approximation scheme is proposed and its asymptotic convergence established. The theoretical development is then extended to pipelines of computation, wherein probabilistic numerical methods are composed to solve more challenging numerical tasks. The contribution highlights an important research frontier at the interface of numerical analysis and uncertainty quantification, with a challenging industrial application presented.

1. Introduction

Probabilistic numerics addresses discretisation error, but existing methods lack a common statistical and Bayesian foundation. The paper establishes foundations for Bayesian methods and their composition in computational pipelines.

  • Motivation: Finite computational budgets force high- and infinite-dimensional problems to be discretised, producing approximations with consequential error.Declining processor improvements and increasing inference complexity make discretisation error increasingly important.
  • Open problem: Probabilistic numerical methods lack a general statistical foundation, making their assumptions and relationships difficult to compare.Existing approaches may model uncertainty in the unknown solution or in discretisation error, without a framework relating them.
  • Contributions: The paper defines Bayesian probabilistic numerical methods explicitly and clarifies how existing non-Bayesian methods relate to them.This framework builds on Bayesian inverse-problem theory and distinguishes among existing approaches.
  • Contributions: It establishes conditions for well-defined PNMs beyond conjugate Gaussian models and develops an asymptotically convergent numerical approximation scheme.The scope includes nonlinear and non-Gaussian models, while efficiency improvements are not the stated aim.
  • Contributions: The paper develops theory for composing PNMs, enabling statistical analysis of how individual discretisations contribute to pipeline error.The contribution introduces conditions under which ANOVA can be justified for composed numerical methods.

2. Probabilistic Numerical Methods

Probabilistic numerical methods represent uncertainty about numerical quantities through belief distributions updated from finite information. The paper formalises Bayesian PNMs using disintegration, addressing why ordinary Bayes’ theorem is problematic for noiseless infinite-dimensional information.

  • General framework: The probabilistic formulation treats the unknown state as fixed while using a random variable to represent epistemic uncertainty.This distinguishes formal uncertainty about the state from randomness in the underlying numerical problem.
  • General framework: A probabilistic numerical method maps finite information about an unknown state to a belief distribution over a quantity of interest.Its information operator records computational information, while its belief update operator produces the output distribution.
  • Examples: The distributional output of a PNM can share its mean with a classical numerical method while additionally representing uncertainty about the exact quantity or solution.This distinction is illustrated for both integration and partial differential equation problems.
  • Bayesian PNMs: For noiseless PNM information, the likelihood is supported on a prior-null set, so standard Bayes’ theorem is generally ill-defined outside conjugate Gaussian infinite-dimensional settings.This motivates a distinct Bayesian PNM definition based on conditioning through disintegration.
  • Bayesian PNMs: Under Radon prior and countably generated information-space conditions, a disintegration exists and is essentially unique, making Bayesian PNMs well-defined.The stated conditions include common settings such as separable complete metric spaces and Euclidean information spaces.

3. Decision-Theoretic Treatment

The paper frames PNM evaluation through Bayes risk, allowing probabilistic and classical numerical methods to be compared and information operators to be optimized. Under squared-norm loss, posterior means are Bayes rules, while optimal information for Bayesian PNMs coincides with average-case optimal information under the theorem's conditions.

  • Loss and Risk: Bayes risk evaluates the expected loss of a PNM's belief update after information about the unknown state is observed.The framework treats the belief update as a randomized decision rule and averages its risk over states drawn from a belief distribution.
  • Loss and Risk: A contraction rate describes how a PNM's Bayes risk decreases as its information operator supplies increasingly many pieces of information.This definition supports comparisons with classical numerical methods and average-case analysis.
  • Bayes Decision Rules: In general, Bayesian PNMs are not Bayes rules because their additional uncertainty increases Bayes risk relative to an appropriate classical numerical method.The paper nevertheless identifies a connection between Bayesian PNMs and Bayes rules through posterior means.
  • Bayes Decision Rules: Under squared-norm loss, returning the posterior mean is a Bayes rule for estimating the quantity of interest.This explains why several Bayesian PNMs are centered on classical numerical methods through their mean functions.
  • Connection to Average Case Analysis: Under squared-norm loss, optimal information for a Bayesian PNM and average-case optimal information are identical.This result allows average-case analysis results to be used to construct optimal Bayesian PNMs, although the paper also notes that the two notions need not coincide in general.
  • Connection to Average Case Analysis: In the integration example, the resulting Bayesian PNM is centered on the average-case optimal method but carries twice the Bayes risk of that method.The example illustrates that optimal information and Bayes-rule optimality are distinct properties.

4. Numerical Disintegration

The paper constructs approximate Bayesian posteriors for numerical computation by replacing exact conditioning with a bandwidth-controlled approximation and finite-dimensional projection. Under continuity and approximation assumptions, these approximations converge to the exact disintegration, while practical sampling uses Monte Carlo methods and introduces a prior-truncation limitation for non-Gaussian models.

  • Approximate disintegration: The approximate posterior µaδ weights states whose information A(x) is close to observed a, generalising approximate Bayesian computation.The bandwidth parameter δ controls the closeness requirement, while φ permits near-matching states rather than only exact matches.
  • Approximate disintegration: Convergence is studied using probability metrics dF that encompass total variation and Wasserstein distance, although useful theory depends on the choice of function space F.Total variation can be uninformative here because distinct conditional distributions may have distance one, and dF may only be a pseudometric.
  • Convergence: Under Hölder continuity of a↦µa and regularity assumptions on the observation density, µaδ converges weakly to µa as δ decreases.The result applies to infinite-dimensional X but relies on the stronger Hölder continuity assumption.
  • Convergence: Finite-dimensional approximations µaδ,N converge when the basis approximation error satisfies CδΨ(N)=o(1), combining bandwidth and projection limits.The total approximation error is bounded using the triangle inequality, with separate contributions from numerical conditioning and finite-dimensional projection.
  • Practical computation: For non-Gaussian priors, experiments truncate the series prior to XN, but the resulting prior perturbation is not analysed because weak convergence does not generally imply total-variation convergence.The paper defers theoretical analysis of this approximation to future work.
  • Practical computation: The approximate distributions can be sampled with standard Monte Carlo methods, including sequential Monte Carlo and parallel tempering.The construction separates sampling from the potentially substantial cost of obtaining the information a itself.

5. Computational Pipelines and PNM

The paper represents composed probabilistic numerical methods as directed graphical pipelines and establishes when their combined output retains a Bayesian interpretation. The key condition is coherence of the prior, expressed through conditional independences in the pipeline’s dependence graph.

  • Motivation: The pipeline framework targets composed numerical procedures because accumulated discretisation error can have a non-trivial structure that is difficult to analyse method by method.Probabilistic formulations can support analysis of variance to identify contributions from individual discretisations.
  • Pipeline representation: A computational pipeline is a directed acyclic bipartite graph connecting information nodes and method nodes, with one terminal principal QoI.Method nodes have labelled information inputs and produce QoIs that can feed later methods.
  • Pipeline representation: Pipelines propagate deterministic or statistical information between compatible PNMs, allowing methods to consume outputs from multiple parent methods.Compatibility requires consistent information spaces and operators across shared information nodes and method connections.
  • Bayesian composition: Redundant representations of the same information are permitted by compatibility but produce non-Bayesian pipelines.The framework does not require the converse of the shared-information condition, so information may be duplicated without being recycled.
  • Bayesian composition: Theorem 5.9 states that a computation of compatible Bayesian PNMs is Bayesian for the principal QoI when the prior is coherent for the pipeline.Coherence is a sufficient and verifiable condition for the composed computation’s Bayesian interpretation.
  • Bayesian composition: For the distributed-integration example, a Wiener-process prior satisfies the required conditional independences, justifying parallel computation of the component integrals from a Bayesian standpoint.The relevant independence follows because the process after the shared state is conditionally independent of its earlier history given that state.

6. Numerical Experiments

The experiments apply Bayesian PNMs to a linear PDE, a nonlinear ODE with multiple solutions, and an industrial process-monitoring pipeline. They demonstrate non-Gaussian beliefs, multimodal posteriors, prior-sensitive inference, and uncertainty propagation through composed numerical methods.

  • Experimental scope: Three experiments address a linear Poisson PDE, a nonlinear Painlevé ODE, and industrial process monitoring with a composed PNM pipeline.The experiments use non-Gaussian belief distributions to extend beyond previous work.
  • 6.1. Poisson Equation: The Poisson experiment specifies beliefs with tensor products of normalised Chebyshev polynomials and constrains the solution using interior, Dirichlet, and Neumann information locations.The spatial domain is two-dimensional, and the finite-element model solution uses a fine mesh.
  • 6.2. The Painlevé ODE: The Painlevé experiment represents solutions with Gaussian or Cauchy coefficient priors whose scales decay geometrically to include both known solutions in prior support.The basis uses normalised Chebyshev polynomials, with γ_i = αβ^-i, α = 8, and β = 1.5.
  • 6.2. The Painlevé ODE: At n = 15, the Bayes factor for the Cauchy prior over the Gaussian prior is 20.26, indicating strong evidence for the Cauchy prior at that discretisation.For n = 17 and δ = 1, the posterior is strongly multimodal with skewed coefficient correlations.
  • 6.2. The Painlevé ODE: As n increases to 22, posterior mass for each prior settles on the positive Painlevé solution, while smaller n can produce multimodal posteriors.At n = 12, the Gaussian-prior posterior becomes trimodal; the positive solution has smaller L2-norm.
  • 6.3. Application to Industrial Process Monitoring: The industrial pipeline uses P = 500 particles and n = 119 design points to recover rotating, mixing conductivity structure across the first eight post-injection frames.Integrated posterior standard deviation is compared between the uncertainty-propagating pipeline and a static approach without propagated uncertainty.

7. Discussion

The discussion frames Bayesian PNMs as a statistically coherent foundation for composing computational pipelines, while identifying unresolved questions about consistency, cost, priors, and extensions.

  • Bayesian PNMs are connected to Bayesian inverse problems to establish when their outputs are well-defined and meaningful.
  • Prior coherence can be established locally through a Markov condition, supporting valid Bayesian interpretations for composed Bayesian PNMs.
  • Approximations or non-Bayesian PNMs may reduce computational costs, but their role in PNM pipelines remains unclear.Numerical disintegration is identified as a potentially substantial offline cost.
  • The framework focuses on accessing information A(x), while numerical disintegration was the predominant computational cost in the reported applications.The authors note that information-access costs may dominate in applications such as climate models, and Monte Carlo disintegration methods could be improved.
  • The framework assumes a supplied belief distribution µ, while robust inference under prior misspecification and coherent multi-expert priors remain open concerns.Generalised Bayes risk with an assessment distribution different from µ is identified as an unconsidered robustness direction.
  • Further work is needed on consistent estimation, randomized information operators, and adaptive Bayesian PNMs.The paper highlights identifiability and concentration rates, random acquisition or measurement error, and adaptation beyond the canonical setting.

A. Proofs

The proofs establish risk identities for Bayesian PNMs and Bayes acts, then show that the pipeline computation agrees with the corresponding posterior distribution under coherence.

  • The proof expresses the joint density of X and A(X) in two equivalent ways as a basis for subsequent arguments.
  • Under integrability and density assumptions, the proof bounds the relevant expression for sufficiently small δ and completes the argument.The reasoning uses disintegration, Fubini’s theorem, and local density geometry.
  • With squared-error loss, the Bayesian PNM’s Bayes risk is related to the Bayes act through the Hilbert-space inner-product structure.The proof compares posterior-distribution output with the point-valued posterior mean.
  • The proof also uses centered variables under µ_a, with posterior mean b(a), to compare Bayesian PNM and Bayes-act risks.
  • The proof identifies the algorithmic pipeline output with the posterior distribution (Q_n)#µ_a on Q_n under the coherence condition.The argument uses conditional-density factorization and does not require existence of a joint density for validity.

S1. Philosophical Status of the Belief Distribution

This section examines what the belief distribution µ means in a PNM, surveying historical interpretations, recent literature, and philosophical objections with counter-arguments.

  • The section treats the semantic status of the belief distribution µ as a central issue in probabilistic numerical methods.
  • It surveys historical work and more recent literature concerning belief distributions in PNMs.
  • It also highlights philosophical objections to belief distributions and discusses corresponding counter-arguments.

S1.1. Historical Precedent

The historical record links probabilistic models of numerical objects to interpolation, rounding-error propagation, Gaussian measures, and the statistical interpretation of numerical analysis.

  • Poincaré used stochastic models for interpolation, effectively constructing a Gaussian measure over a finite-dimensional Hilbert-space basis.Later work connected spline interpolation formulae with means of Gaussian measures over Hilbert spaces.
  • Early ODE work modeled rounding error with uniform perturbations and propagated the resulting numerical uncertainty through subsequent computation.This line of work focused on rounding error rather than time-discretisation error.
  • Larkin’s agenda treated numerical analysis as statistical estimation using specific information and assumed general properties of solution functions.The permissible class could be a Hilbert space or be simply related to one.
  • Diaconis emphasized that finite information determines some properties of a function while leaving others, such as monotonicity or convexity, unknown.The associated Gaussian belief distributions encoded basic features such as smoothness, and some classical methods were Bayes rules.
  • Theoretical analysis of computation under finite information later proceeded largely in applied mathematics without an explicit statistical context.

S1.2. Contemporary Outlook

Information-based complexity motivates optimal approximation from finite information, while Bayesian probabilistic numerics connects this perspective to priors, losses, and information operators.

  • Information-based complexity studies computation from finite information, motivating optimal approximations under worst-case or average-case criteria.
  • The paper’s probabilistic numerical methods align particularly with average-case analysis and focus on that literature.
  • Earlier average-case work studied optimal numerical integration, function approximation, covariance-kernel classes, and spline optimality.
  • Game-theoretic analysis proposed a Gaussian prior for probabilistic meshless methods, offering one route toward default or objective priors for PNMs.
  • The paper frames practical PNM design around selecting information operators and Bayesian beliefs within this broader foundation.

S1.3. Paradise Lost?

Numerical algorithms contain multiple discretisation-error sources, creating a modeling choice about which errors probabilistic numerical methods should represent.

  • Finite element wave-equation solvers combine spatial and temporal discretisation, numerical quadrature, and finite-precision arithmetic.
  • PNM implementations omitted rounding error because floating-point error was insignificant relative to discretisation error in the examples.
  • The paper treats choosing which discretisation errors to model as a pragmatic simplification of the real world.

S2. Existence of Non-Randomised Bayes Rule

The existence result shows that whenever admissible belief-update operators exist, at least one classical numerical method is also admissible.

  • Sufficient conditions for Fubini’s theorem are assumed while establishing existence of non-randomised Bayes rules.
  • If B(A) is non-empty, it contains a classical numerical method B(µ, a) = δ ◦b(a), with b(a) a Bayes act for each a.
  • The proof represents any non-classical operator as a distribution over classical operators and uses its risk to obtain a classical operator with no worse risk.

S3. Optimal Information: A Counterexample

This material contrasts Bayesian PNM and average-case optimal information, then develops SMC-based numerical disintegration with convergence conditions and explicit failure boundaries.

  • Optimal Information: A Counterexample: In the counterexample, ACA optimal information can use S = {♠, ♦} or S = {♠, ♦, ♥}, whereas PNM optimal information is S = {♠, ♦}.
  • Optimal Information: A Counterexample: Bayesian PNM optimal information need not coincide with average-case optimal information, as shown by a discrete 0-1-loss counterexample.
  • Sequential Monte Carlo Algorithms for Numerical Disintegration: Monte Carlo approximation of µaδ is a rare-event problem because most posterior mass may lie in a set with small prior probability.
  • Sequential Monte Carlo Algorithms for Numerical Disintegration: SMC-ND uses particle initialization, weighting, propagation, and resampling to produce an empirical approximation, with related tempering approaches including splitting and parallel tempering.
  • Sequential Monte Carlo Algorithms for Numerical Disintegration: Uniform convergence is established under positivity and transition-kernel assumptions, with bounds independent of δ and extensions left for future work.
  • Case B: Indicator Function: For indicator weighting, the algorithm can fail when no particle enters the next tolerance set, so successive δi values must be sufficiently close to control failure probability.
Loading 1702.03673v2…