Bayesian Probabilistic Numerical Methods
Jon Cockayne∗

Chris Oates†

Tim Sullivan‡

Mark Girolami§

arXiv:1702.03673v2 [stat.ME] 7 Jul 2017

July 10, 2017

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
Numerical computation underpins almost all of modern scientific and industrial research and
development. The impact of a finite computational budget is that problems whose solutions are
high- or infinite-dimensional, such as the solution of differential equations, must be discretised in
order to be solved. The result is an approximation to the object of interest. The declining rate
of processor improvement as physical limits are reached is in contrast to the surge in complexity
of modern inference problems, and as a result the error incurred by discretisation is attracting
increased interest (e.g. Capistrán et al., 2016).
The situation is epitomised in modern climate models, where use of single-precision arithmetic
has been explored to permit finer temporal resolution. However, when computing in singleprecision, a detailed time discretisation can increase total error, due to the increased number
of single precision computations, and in practice some form of ad-hoc trade-off is sought (Harvey and Verseghy, 2015). It has been argued that statistical considerations can permit more
principled error control strategies for such models (Hennig et al., 2015).
Numerical methods are designed to mitigate discretisation errors of all forms (Press et al.,
2007). Nonetheless, the introduction of error is unavoidable and it is the role of the numerical
analyst to provide control of this error (Oberkampf and Roy, 2013). The central theoretical
results of numerical analysis have in general not been obtained through statistical considerations.
∗

University of Warwick, j.cockayne@warwick.ac.uk
Newcastle University and Alan Turing Institute, chris.oates@ncl.ac.uk
‡
Free University of Berlin and Zuse Institute Berlin, sullivan@zib.de
§
Imperial College London and Alan Turing Institute, m.girolami@imperial.ac.uk
†

1

More recently, the connection of discretisation error to statistics was noted as far back as
Henrici (1963); Hull and Swenson (1966), who argued that discretisation error can be modelled
using a series of independent random perturbations to standard numerical methods. However,
numerical analysts have cast doubt on this approach, since discretisation error can be highly
structured; see Kahan (1996) and Higham (2002, Section 2.8). To address these objections, the
field of probabilistic numerics has emerged with the aim to properly quantify the uncertainty
introduced through discretisation in numerical methods.
The foundations of probabilistic numerics were laid in the 1970s and 1980s, where an important shift in emphasis occurred from the descriptive statistical models of the 1960s to the use
of formal inference modalities that generalise across classes of numerical tasks. In a remarkable
series of papers, Larkin (1969, 1970, 1972); Kuelbs et al. (1972); Larkin (1974, 1979a,b), Mike
Larkin presented now classical results in probabilistic numerics, in particular establishing the
correspondence between Gaussian measures on Hilbert spaces and optimal numerical methods.
Re-discovered and re-emphasised on a number of occasions, the role for statisticians in this new
outlook was clearly captured in Kadane and Wasilkowski (1985):
Statistics can be thought of as a set of tools used in making decisions and inferences
in the face of uncertainty. Algorithms typically operate in such an environment.
Perhaps then, statisticians might join the teams of scholars addressing algorithmic
issues.
The 1980s culminated in development of Bayesian optimisation methods (Mockus, 1989; Törn
and Žilinskas, 1989), as well as the relation of smoothing splines to Bayesian estimation (Kimeldorf and Wahba, 1970b; Diaconis and Freedman, 1983).
The modern notion of a probabilistic numerical method (henceforth PNM) was described in
Hennig et al. (2015); these are algorithms whose output is a distribution over an unknown,
deterministic quantity of interest, such as the numerical value of an integral. Recent research in
this field includes PNMs for numerical linear algebra (Hennig, 2015; Bartels and Hennig, 2016),
numerical solution of ordinary differential equations (ODEs; Schober et al., 2014; Kersting and
Hennig, 2016; Schober et al., 2016; Conrad et al., 2016; Chkrebtii et al., 2016), numerical
solution of partial differential equations (PDEs; Owhadi, 2015; Cockayne et al., 2016; Conrad
et al., 2016) and numerical integration (O’Hagan, 1991; Briol et al., 2016).
Open Problems Despite numerous recent successes and achievements, there is currently no
general statistical foundation for PNMs, due to the infinite-dimensional nature of the problems
being solved. For instance, at present it is not clear under what conditions a PNM is welldefined, except for in the standard conjugate Gaussian framework considered in (Larkin, 1972).
This limits the extent to which domain-specific knowledge, such as boundedness of an integrand
or monotonicity of a solution to a differential equation, can be encoded in PNMs. In contrast,
classical numerical methods often exploit such information to achieve substantial reduction in
discretisation error. For instance, finite element methods for solution of PDEs proceed based
on a mesh that is designed to be more refined in areas of the domain where greater variation of
the solution is anticipated (Strang and Fix, 1973).
Furthermore, although PNMs have been proposed for many standard numerical tasks (see
Section 2.6.1), the lack of common theoretical foundations makes comparison of these methods
difficult. Again taking PDEs as an example, Cockayne et al. (2016) placed a probability distribution on the unknown solution of the PDE, whereas Conrad et al. (2016) placed a probability
distribution on the unknown discretisation error of a numerical method. The uncertainty modelled in each case is fundamentally different, but at present there is no framework in which to

2

articulate the relationship between the two approaches. Furthermore, though PNMs are often
reported as being “Bayesian” there is no clear definition of what this ought to entail.
A more profound consequence of the lack of common foundation occurs when we seek to
compose multiple PNMs. For example, multi-physics cardiac models involve coupled ODEs and
PDEs which must each be discretised and approximately solved to estimate a clinical quantity of
interest (Niederer et al., 2011). The composition of successive discretisations leads to non-trivial
error propagation and accumulation that could be quantified, in a statistical sense, with PNMs.
However, proper composition of multiple PNMs for solutions of ODEs and PDEs requires that
these PNMs share common statistical foundations that ensure coherence of the overall statistical
output. These foundations remain to be established.
Contributions The main contribution of this paper is to establish rigorous foundations for
PNMs:
The first contribution is to argue for an explicit definition of a “Bayesian” PNM. Our framework generalises the seminal work of Larkin (1972) and builds on the modern and popular
mathematical framework of Stuart (2010). This illuminates subtle distinctions among existing
methods and clarifies the sense in which non-Bayesian methods are approximations to Bayesian
PNMs.
The second contribution is to establish when PNMs are well-defined outside of the conjugate
Gaussian context. For exploration of non-linear, non-Gaussian models, a numerical approximation scheme is developed and shown to asymptotically approach the posterior distribution of
interest. Our aim here is not to develop new or more computationally efficient PNMs, but to
understand when such development can be well-defined.
The third contribution is to discuss pipelines of composed PNMs. This is a critical area of
development for probabilistic numerics; in isolation, the error of a numerical method can often
be studied and understood, but when composed into a pipeline the resulting error structure may
be non-trivial and its analysis becomes more difficult. The real power of probabilistic numerics
lies in its application to pipelines of numerical methods, where the probabilistic formulation
permits analysis of variance (ANOVA) to understand the contribution of each discretisation to
the overall numerical error. This paper introduces conditions under which a composition of
PNMs can be considered to provide meaningful output, so that ANOVA can be justified.
Structure of the Paper In Section 2 we argue for an explicit definition of Bayesian PNM
and establish when such methods are well-defined. Section 3 establishes connections to other
related fields, in particular with relation to evaluating the performance of PNMs. In Section 4 we
develop useful numerical approximations to the output of Bayesian PNMs. Section 5 develops
the theory of composition for multiple PNMs. Finally, in Section 6 we present applications of
the techniques discussed in this paper.
All proofs can be found in either the Appendix or the Electronic Supplement.

2. Probabilistic Numerical Methods
The aim of this section is to provide rigorous statistical foundations for PNMs.

2.1. Notation
For a measurable space (X , ΣX ), the shorthand PX will be used to denote the set of all distributions on (X , ΣX ). For µ, µ0 ∈ PX we write µ  µ0 when µ is absolutely continuous with respect

3

to µ. The notation δ(x) will be used to denote a Dirac measure on x ∈ X , so that δ(x) ∈ PX . Let
1[S] denote the indicator function of an event S ∈ ΣX . For a measurable function
f :X →R
R
and a distribution µ ∈ PX , we will on occasion use the notation µ(f ) = f (x)µ(dx) and
kf k∞ = supx∈X |f (x)|. The point-wise product of two functions f and g is denoted f · g. For a
function or operator T , T# denotes the associated push-forward operator1 that acts on measures
on the domain of T . Let ⊥
⊥ denote conditional
independence. The subset `p ⊂ R∞ is defined
P∞
to consist of sequences (ui ) for which i=1 |ui |p is convergent. C(0, 1) will be used to denote
the set of continuous functions on (0, 1).

2.2. Definition of a PNM
To first build intuition, consider numerical approximation of the Lebesgue integral
Z
x(t)ν(dt)
for some integrable function x : D → R, with respect to a measure ν on D. Here we may
directly interrogate the integrand x(t) at any t ∈ D, but unless D is finite we cannot evaluate
x at all t ∈ D with a finite computational budget. Nonetheless, there are many algorithms for
approximation of this integral based on information {x(ti )}ni=1 at some collection of locations
{ti }ni=1 .
To see the abstract structure of this problem, assume the state variable x exists in a measurable space (X , ΣX ). Information about x is provided through an information operator
A : X → A whose range is a measurable space (A, ΣA ). Thus, for the Lebesgue integration
problem, the information operator is


x(t1 )


A(x) =  ...  = a ∈ A.
(2.1)
x(tn )

The space X , in this case a space of functions, can be high- or infinite-dimensional, but the space
A of information is assumed to be finite-dimensional in accordance with our finite computational
budget. In this paper we make explicit a quantity of interest (QoI) Q(x), defined by a map
Q : X → Q into a measurable space (Q, ΣQ ). This captures that x itself may not be the object
of interest for the numerical
problem; for the Lebesgue integration illustration, the QoI is not
R
x itself but Q(x) = x(t)ν(dt).
The standard approach to such computational problems is to construct an algorithm which,
when applied, produces some approximation q̂(a) of Q(x) based on the information a, whose
theoretical convergence order can be studied. A successful algorithm will often tailor the information operator A to the QoI Q. For example, classical Gaussian cubature specifies sigma
points {t∗i }ni=1 at which the integrand must be evaluated, based on exact integration of certain
polynomial test functions.
The probabilistic numerical approach, instead, begins with the introduction of a random
variable X on (X , ΣX ). The true state X = x is fixed but unknown; the randomness is used
an abstract device used to represent epistemic uncertainty about x prior to evolution of the
information operator (Hennig et al., 2015). This is now formalised:
1

Recall that, for measurable T : X → A, the pushforward T# µ of a distribution µ ∈ PX is defined as T# µ(A) =
µ(T −1 (A)) for all A ∈ ΣA .

4

Definition 2.1 (Belief Distribution). An element µ ∈ PX is a belief distribution 2 for x if it
carries the formal semantics of belief about the true, unknown state variable x.
Thus we may consider µ to be the law of X. The construction of an appropriate belief distribution µ for a specific numerical task is not the focus of this research and has been considered
in detail in previous work; see the Electronic Supplement for an overview of this material.
Rather we consider the problem of how one updates the belief distribution µ in response to the
information A(x) = a obtained about the unknown x. Generic approaches to update belief distributions, which generalise Bayesian inference beyond the unique update demanded in Bayes
theorem, were formalised in Bissiri et al. (2016); de Carvalho et al. (2017).
Definition 2.2 (Probabilistic Numerical Method). Let (X , ΣX ), (A, ΣA ) and (Q, ΣQ ) be measurable spaces and let A : X → A, Q : X → Q and B : PX × A → PQ where A and Q are
measurable functions. The pair M = (A, B) is called a probabilistic numerical method for estimation of a quantity of interest Q. The map A is called an information operator, and the map
B is called a belief update operator.
The output of a PNM is a distribution B(µ, a) ∈ PQ . This holds the formal status of a belief
distribution for the value of Q(x), based on both the initial belief µ about the value of x and
the information a that are input to the PNM.
An objection sometimes raised to this construction is that x itself is not random. We emphasise that this work does not propose that x should be considered as such; the random variable
X is a formal statistical device used to represent epistemic uncertainty (Kadane, 2011; Lindley,
2014). Thus, there is no distinction from traditional statistics, in which x represents a fixed but
unknown parameter and X encodes epistemic uncertainty about this parameter.
Before presenting specific instances of this general framework, we comment on the potential
analogy between A and the likelihood function, and between B and Bayes’ theorem. Whilst
intuitively correct, the mathematical developments in this paper are not well-suited to these
terms; in Section 2.5 we show that Bayes formula is not well-defined, as the posterior distribution
is not absolutely continuous with respect to the prior.
To strengthen intuition we now give specific examples of established PNMs:
Example 2.3 (Probabilistic Integration). Consider the numerical integration problem earlier
discussed. Take D ⊆ Rd , X a separable Banach space of real-valued functions on D, and ΣX
the Borel σ-algebra for X . The space (X , ΣX ) is endowed with a Gaussian belief distribution
µ ∈ PX . Given information A(x) = a, define µa to be the restriction of µ to those functions
which interpolate x at the points {ti }ni=1 ; that µa is again Gaussian follows from linearity
of the
R
information operator (see Bogachev, 1998, for details). The QoI Q remains Q(x) = x(t)ν(dt).
This problem was first considered by Larkin (1972). The belief update operator proposed
therein, and later considered in Diaconis (1988); O’Hagan (1991) and others, was B(µ, a) =
Q# µa . Since Gaussians are closed under linear projection, the PNM output B(µ, a) is a univariate Gaussian whose mean and variance can be expressed in closed-form for certain choices of
Gaussian covariance function and reference measure ν on D. Specifically, if µ has mean function
m : X → R and covariance function k : X × X → R, then
B(µ, a) = N(z > K−1 (a − m̄), z0 − z > K−1 z)
2

(2.2)

Two remarks are in order: First, we have avoided the use of “prior” as this abstract framework encompasses
both Bayesian and non-Bayesian PNMs (to be defined). Second, the use of “belief” differs to the set-valued
belief functions in Dempster–Shafer theory, which do not require that µ(E) + µ(E c ) = 1 (Shafer, 1976).

5

R
where m̄, z ∈ Rn are defined
as m̄i = m(ti ), zi = k(t, ti )ν(dt), K ∈ Rn×n is defined as
RR
Ki,j = k(ti , tj ) and z0 =
k(t, t0 )(ν × ν)d(t × t0 ) ∈ R. This method was extensively studied in
Briol et al. (2016), who provided a listing of (ν, k) combinations for which z and z0 possess a
closed-form.
An interesting fact is that the mean of B(µ, a) coincides with classical cubature rules for
different choices of µ and A (Diaconis, 1988; Särkkä et al., 2016). In Section 3 we will show
that this is a typical feature of PNMs. The crucial distinction between PNMs and classical
numerical methods is the distributional nature of B(µ, a), which carries the formal semantics
of belief about the QoI. The full distribution B(µ, a) was examined in Briol et al. (2016), who
established contraction to the exact value of the integral under certain smoothness conditions
on the Gaussian covariance function and on the integrand. See also Kanagawa et al. (2016);
Karvonen and Särkkä (2017).
Example 2.4 (Probabilistic Meshless Method). As a canonical example of a PDE, take the
following elliptic problem with Dirichlet boundary conditions
−∇ · (κ∇x) = f

in D

x=g

on ∂D

(2.3)

where we assume D ⊂ Rd and κ : D → Rd×d is a known coefficient. Let X be a separable
Banach space of appropriately differentiable real-valued functions and take ΣX to be the Borel
σ-algebra for X . In contrast to the first illustration, the QoI here is Q(x) = x, as the goal is to
make inferences about the solution of the PDE itself.
Such problems were considered in Cockayne et al. (2016) wherein µ was restricted to be a
Gaussian distribution on X . The information operator was constructed by choosing finite sets
of locations T1 = {t1,1 , . . . , t1,n1 } ⊂ D and T2 = {t2,1 , . . . , t2,n2 } ⊂ ∂D at which the system
defined in Eq. (2.3) was evaluated, so that




−∇ · (κ(t1,1 )∇x(t1,1 ))
f (t1,1 )


 .. 
..


 . 
.




−∇ · (κ(t1,n1 )∇x(t1,n1 ))
f (t1,n1 )



 .
A(x) = 
a=


x(t2,1 )


 g(t2,1 ) 


 .. 
..


 . 
.
x(t2,n2 )
g(t2,n2 )

The belief update operator was chosen to be B(µ, a) = µa , where µa is the restriction of µ to
those functions for which A(x) = a is satisfied. In the setting of a linear system of PDEs such
as that in Eq. (2.3), the distribution B(µ, a) is again Gaussian (Bogachev, 1998). Full details
are provided in Cockayne et al. (2016).
As in the previous example, we note that the mean of B(µ, a) coincides with the numerical solution to the PDE provided by a classical method (the symmetric collocation method;
Fasshauer, 1999). The full distribution B(µ, a) provides uncertainty quantification for the unknown exact solution and can again be shown to contract to the exact solution under certain
smoothness conditions (Cockayne et al., 2016). This method was further analysed for a specific
choice of covariance operator in the belief distribution µ, in an impressive contribution from
Owhadi (2017).

6

2.2.1. Classical Numerical Methods
Standard numerical methods fit into the above framework, as can be seen by taking
B(µ, a) = δ ◦ b(a)

(2.4)

independent of the distribution µ, where a function b : A → Q gives the output of some classical
numerical method for solving the problem of interest. Here δ : Q → PQ maps b(a) ∈ Q to a
Dirac measure centred on b(a). Thus, information in a ∈ A is used to construct a point estimate
b(a) ∈ Q for the QoI.
The formal language of probabilities is not used in classical numerical analysis to describe
numerical error. However, in many cases the classical and probabilistic analyses are mathematically equivalent. For instance, there is an equivalence between the standard deviation of
B(µ, a) for probabilistic integration and the worst-case error for numerical cubature rules from
numerical analysis (Novak and Woźniakowski, 2010). The explanation for this phenomenon will
be given in Section 3.

2.3. Bayesian PNMs
Having defined a PNM, we now state the central definition of this paper, that is of a Bayesian
PNM. Define µa to be the conditional distribution of the random variable X, given the event
A(X) = a. For now we assume that this can be defined without ambiguity and reserve a more
technical treatment of conditional probabilities for Section 2.5.
In this work we followed Larkin (1972) and cast the problem of determining x in Eq. (2.1)
as a problem of Bayesian inversion, a framework now popular in applied mathematics and
uncertainty quantification research (Stuart, 2010). However, in a standard Bayesian inverse
problem the observed quantity a is assumed to be corrupted with measurement error, which is
described by a “likelihood”. This leads, under mild assumptions, to general versions of Bayes’
theorem (see Stuart, 2010, Section 2.2)
For PNM, however, the information is not corrupted with measurement error. As a result,
the support of the likelihood is a null set under the prior, making the standard approaches
to such problems, including Bayes’ theorem, ill-defined outside of the conjugate Gaussian case
when unknowns are infinite-dimensional. This necessitates a new definition:
Definition 2.5 (Bayesian Probabilistic Numerical Method). A probabilistic numerical method
M = (A, B) is said to be Bayesian 3 for a quantity of interest Q if, for all µ ∈ PX , the output
B(µ, a) = Q# µa ,

for A# µ-almost-all a ∈ A.

That is, a PNM is Bayesian if the output of the PNM is the push-forward of the conditional
distribution µa through Q. This definition is familiar from the examples in Section 2.2, which
are both examples of Bayesian PNMs.
For Bayesian PNMs we adopt the traditional terminology in which µ is the prior for x and
the output Q# µa the posterior for Q(x). Note that, for fixed A and µ, the Bayesian choice of
belief update operator B (if it exists) is uniquely defined.
It is emphasised that the class of Bayesian PNMs is a subclass of all PNMs; examples of nonBayesian PNMs are provided in Section 2.6.1. Our analysis is focussed on Bayesian PNMs due to
3

The use of “Bayesian” contrasts with Bissiri et al. (2016), for whom all belief update operators represent
Bayesian learning algorithms to some greater or lesser extent. An alternative term could be “lossless”, since
all the information in a is conditioned upon in µa .

7

their appealing Bayesian interpretation and ease of generalisation to pipelines of computation in
Section 5. For non-Bayesian PNMs, careful definition and analysis of the belief update operator
is necessary to enable proper interpretation of the uncertainty quantification being provided.
In particular, the analysis of non-Bayesian PNMs may present considerable challenges in the
context of computational pipelines, whereas for Bayesian PNMs this is shown in Section 5 to
be straight-forward.

2.4. Model Evidence
A cornerstone of the Bayesian framework is the model evidence, or marginal likelihood (MacKay,
1992). Let A ⊆ Rn be equipped with the Lebesgue reference measure λ, such that A# µ admits
a density pA = dA# µ/dλ. Then the model evidence pA (a), based on the information that
A(x) = a, can be used as the basis for Bayesian model comparison. In particular, two prior
distributions µ, µ̃, can be compared through the Bayes factor
BF :=

dA# µ̃
p̃A (a)
=
(a),
pA (a)
dA# µ

(2.5)

where p̃A = dA# µ̃/dλ. Here the second expression is independent of the choice of reference
measure λ and is thus valid for general A. The model evidence has been explored in connection
with the design of Bayesian PNM. For the integration and PDE examples 2.3 and 2.4, the model
evidence has a closed form and was investigated in Briol et al. (2016); Cockayne et al. (2016).
In Section 6 we investigate the model evidence in the context of non-linear ODEs and PDEs
for which it must be approximated.

2.5. The Disintegration Theorem
The purpose of this section is to formalise µa and to determine conditions under which µa exists
and is well-defined. From Definition 2.5, the output of a Bayesian PNM is B(µ, a) = Q# µa . If
µa exists, the pushforward Q# µa exists as Q is assumed to be measurable; thus, in this section,
we focus on the rigorous definition of µa .
Unlike many problems of Bayesian inversion, proceeding by an analogue of Bayes’ theorem
is not possible. Let X a = {x ∈ X : A(x) = a}. Then we observe that, if it is measurable, X a
may be a set of zero measure under µ. Standard techniques for infinite-dimensional Bayesian
inversion rely on constructing a posterior distribution based on its Radon–Nikodým derivative
with respect to the prior (Stuart, 2010). However, when µa 6 µ no Radon–Nikodým derivative
exists and we must turn to other approaches to establish when a Bayesian PNM is well-defined.
Conditioning on null sets is technical and was formalised in the celebrated construction of
measure-theoretic probability by Kolmogorov (1933). The central challenge is to establish
uniqueness of conditional probabilities. For this work we exploit the disintegration theorem to
ensure our constructions are well-defined. The definition below is due to Dellacherie and Meyer
(1978, p.78), and a statistical introduction to disintegration can be found in Chang and Pollard
(1997).
Definition 2.6 (Disintegration). For µ ∈ PX , a collection {µa }a∈A ⊂ PX is a disintegration of
µ with respect to the (measurable) map A : X → A if:
1 (Concentration:) µa (X \ X a ) = 0 for A# µ-almost all a ∈ A;
and for each measurable f : X → [0, ∞) it holds that

8

2 (Measurability:) a 7→ µa (f ) is measurable;
R
3 (Conditioning:) µ(f ) = µa (f )A# µ(da).

The concept of disintegration extends the usual concept of conditioning of random variables
to the case where X a is a null set, in a way closely related to regular conditional distributions
(Kolmogorov, 1933). Existence of disintegrations is guaranteed under general weak conditions:
Theorem 2.7 (Disintegration Theorem; Thm. 1 of Chang and Pollard (1997)). Let X be a
metric space, ΣX be the Borel σ-algebra and µ ∈ PX be Radon. Let ΣA be countably generated
and contain all singletons {a} for a ∈ A. Then there exists a disintegration {µa }a∈A of µ with
respect to A. Moreover, if {ν a }a∈A is another such disintegration, then {a ∈ A : µa 6= ν a } is a
A# µ null set.
The requirement that µ is Radon is weak and is implied when X is a Radon space, which
encompasses, for example, separable complete metric spaces. The requirement that ΣA is
countably generated is also weak and includes the standard case where A = Rn with the Borel
σ-algebra. From Theorem 2.7 it follows that {µa }a∈A exists and is essentially unique for all of
the examples considered in this paper. Thus, under mild conditions, we have established that
Bayesian PNMs are well-defined, in that an essentially unique disintegration {µa }a∈A exists. It
is noted that a variational definition of µa has been posited as an alternative approach, for when
the existence of a disintegration is difficult to establish (p3 of Garcia Trillos and Sanz-Alonso,
2017).

2.6. Prior Construction
The Gaussian distribution is popular as a prior in the PNM literature for its tractability, both in
the fact that finite-dimensional distributions take a closed-form and that an explicit conditioning
formula exists. More general priors, such as Besov priors (Dashti et al., 2012) and Cauchy priors
(Sullivan, 2016) are less easily accessed. In this section we summarise a common construction
for these prior distributions, designed to ensure that a disintegration will exist.
Let {φi }∞
i=0 denote an orthogonal Schauder basis for X , assumed to be a separable Banach
space in this section. Then any x ∈ X can be represented through an expansion
x = x0 +

∞
X

ui φi

(2.6)

i=0

for some fixed element x0 ∈ X and a sequence u ∈ R∞ . Construction of measures µ ∈ PX
is then reduced to construction of almost-surely convergent measures on R∞ and studying the
pushforward of such measures into X . In particular, this will ensure that µ ∈ PX is Radon
(as X is a separable complete metric space), a key requirement for existence of a disintegration
{µa }a∈A .
To this end it is common to split u into a stochastic and deterministic component; let ξ ∈ R∞
represent an i.i.d sequence of random variables, and γ ∈ `p for some p ∈ (1, ∞). Then with ui =
γi ξi , for the prior distribution to be well-posed we require that almost-surely u ∈ `1 . Different
choices of (ξ, γ) give rise to different distributions on X . For instance, ξi ∼ Uniform(−1, 1),
γ ∈ `1 is termed a uniform prior and ξi ∼ N (0, 1) gives a Gaussian prior, where γ determines
the regularity of the covariance operator C (Bogachev, 1998). The choice of ξi ∼ Cauchy(0, 1)
gives a Cauchy prior in the sense of Sullivan (2016); here we require γ ∈ `1 ∩ ` log ` for X a
separable Banach space, or γ ∈ `2 for when X is a Hilbert space.
A range of prior specifications will be explored in Section 6, including non-Gaussian prior
distributions for numerical solution of nonlinear ODEs.

9

2.6.1. Dichotomy of Existing PNMs
This section concludes with an overview of existing PNMs with respect to our definition of a
Bayesian PNM. This serves to clarify some subtle distinctions in existing literature, as well as
to highlight the generality of our framework. To maintain brevity we have summarised our
findings in Table 1.

3. Decision-Theoretic Treatment
Next we assess the performance of PNMs from a decision-theoretic perspective (Berger, 1985)
and explore connections to average-case analysis of classical numerical methods (Ritter, 2000).
Note that the treatment here is agnostic to whether the PNM in question is Bayesian, and
also encompasses classical numerical methods. Throughout, the existence of a disintegration
{µa }a∈A will be assumed.

3.1. Loss and Risk
Consider a generic loss function L : Q × Q → R where L(q † , q) describes the loss incurred when
the true QoI q † = Q(x) is estimated with q ∈ Q. Integrability of L is assumed.
The belief update operator B returns a distribution over Q which can be cast as a randomised
decision rule for estimation of q † . For randomised decision rules, the risk function r : Q×PQ → R
is defined as
Z
†
r(q , ν) = L(q † , q)ν(dq) .
The average risk of the PNM M = (A, B) with respect to µ ∈ PX is defined as
Z
R(µ, M ) = r(Q(x), B(µ, A(x)))µ(dx).

(3.1)

Here a state x ∼ µ is drawn at random and the risk of the PNM output B(µ, A(x)) is computed.
We follow the convention of terming R(µ, M ) the Bayes risk of the PNM, though the usual
objection that a frequentist expectation enters into the definition of the Bayes risk could be
raised.
Next, we consider a sequence A(n) of information operators indexed such that A(n) (x) is
n-dimensional (i.e. n pieces of information are provided about x).
Definition 3.1 (Contraction). A sequence M (n) = (A(n) , B (n) ) of PNMs is said to contract at
a rate rn under a belief distribution µ if R(µ, M (n) ) = O(rn ).
This definition allows for comparison of classical and probabilistic numerical methods (Kadane
and Wasilkowski, 1983; Diaconis, 1988). In each case an important goal is to determine methods
M (n) that contract as quickly as possible for a given distribution µ that defines the Bayes risk.
This is the approach taken in average-case analysis (ACA; Ritter, 2000) and will be discussed
in Section 3.4. For Examples 2.3 and 2.4 of Bayesian PNMs, Briol et al. (2016) and Cockayne
et al. (2016) established rates of contraction for particular prior distributions µ; we refer the
reader to those papers for details.

10

Method
Integrator

QoI Q(x)
R
x(t)ν(dt)

Information A(x)
{x(ti )}n
i=1

Non- (or Approximate) Bayesian PNMs
Osborne et al. (2012b,a); Gunter et al. (2014)

{ti }n
i=1 s.t. ti ∼ x
{(ti , x1 (ti ))}n
i=1 s.t. ti ∼ x2
{x(ti )}n
i=1
{∇x(ti )}n
i=1
{(x(ti ), ∇x(ti )}n
i=1

Kong et al. (2003); Tan (2004); Kong et al. (2007)

Optimiser

R
R f (t)x(dt)
x1 (t)x2 (dt)
arg min x(t)

Oates et al. (2016a)
Bayesian Optimisation (Mockus, 1989)
Hennig and Kiefel (2013)
Probabilistic Line Search (Mahsereci and Hennig,
2015)
Probabilistic Bisection Algorithm (Horstein,
1963)

{I[tmin < ti ]}n
i=1
Linear Solver

x

ODE Solver

x

−1

b

{I[tmin < ti ] + error}n
i=1
{xti }n
i=1

Waeber et al. (2013)

{∇x(ti )}n
i=1

(Skilling, 1992)
Filtering Methods for IVPs (Schober et al., 2014;
Chkrebtii et al., 2016; Kersting and Hennig, 2016;
Teymur et al., 2016; Schober et al., 2016)
Finite Difference Methods (John and Wu, 2017)
Hull and Swenson (1966); Mosbach and Turner
(2009)
Stochastic Euler (Krebs, 2016)
Chkrebtii et al. (2016); Raissi et al. (2017)

11

∇x + rounding error
PDE Solver

x(tend )
x

{∇x(ti )}n
i=1
{Dx(ti )}n
i=1
Dx + discretisation error

Bayesian PNMs
Bayesian Quadrature (Larkin, 1974; Diaconis,
1988; O’Hagan, 1991)

Probabilistic Linear Solvers (Hennig, 2015; Bartels and Hennig, 2016)

Probabilistic Meshless Methods (Owhadi, 2015,
2017; Cockayne et al., 2016; Raissi et al., 2016)

Conrad et al. (2016)

Table 1: Comparison of several existing Probabilistic Numerical Methods (PNMs).

3.2. Bayes Decision Rules
A (possibly randomised) decision rule is said to be a Bayes rule if it achieves the minimum
Bayes risk among all decision rules. In the context of (not necessarily Bayesian) PNMs, let
M = (A, B) and let


0
B(A) = B : R(µ, (A, B)) = inf0 R(µ, (A, B )) .
B

That is, for fixed A, B(A) is the set of all belief update operators that achieve minimum Bayes
risk.
This raises the natural question of which belief update operators yield Bayes rules. Although
the definition of a Bayes rule applies generically to both probabilistic and deterministic numerical methods, it can be shown4 that if B(A) is non-empty, then there exists a B ∈ B(A) which
takes the form of a classical numerical method, as expressed in Eq. (2.4). Thus in general,
Bayesian PNMs do not constitute Bayes rules, as the extra uncertainty inflates the Bayes risk,
so that such methods are not optimal.
Nonetheless, there is a natural connection between Bayesian PNMs and Bayes rules, as exposed in Kadane and Wasilkowski (1983):
Theorem 3.2. Let M = (A, B) be a Bayesian probabilistic numerical method for the QoI Q.
Let (Q, h·, ·iQ ) be an inner-product space and let the loss function L have the form L(q † , q) =
kq † − qk2Q , where k · kQ is the norm induced by the inner product. Then the decision rule that
returns the mean of the distribution B(µ, a) is a Bayes rule for estimation of q † .
This well-known fact from Bayesian decision theory5 is interesting in light of recent research in
constructing PNMs whose mean functions correspond to classical numerical methods (Schober
et al., 2014; Hennig, 2015; Särkkä et al., 2016; Teymur et al., 2016; Schober et al., 2016).
Theorem 3.2 explains the results in Examples 2.3 and 2.4, in which both instances of Bayesian
PNMs were demonstrated to be centred on an established classical method.

3.3. Optimal Information
The previous section considered selection of the belief update operator B, but not of the information operator A. The choice of A determines the Bayes risk for a PNM, which leads to a
problem of experimental design to minimise that risk.
The theoretical study of optimal information is the focus of the information complexity literature (Traub et al., 1988; Novak and Woźniakowski, 2010), while other fields such as quasi-Monte
Carlo (QMC, Dick and Pillichshammer, 2010) attempt to develop asymptotically optimal information operators for specific numerical tasks, such as the choice of evaluation points for
numerical approximation of integrals in the case of QMC. Here we characterise optimal information for Bayesian PNMs.
Consider the choice of A from a fixed subset Λ of the set of all possible information operators.
To build intuition, for the task of numerical integration, Λ could represent all possible choices of
locations {ti }ni=1 where the integrand is evaluated. For Bayesian PNM, one can ask for optimal
information:

Aµ ∈ arg inf R(µ, M ) s.t. M = (A, B), B = Q# µA
A∈Λ

4
5

The proof is included in the Electronic Supplement.
This is the fact that the Bayes act is the posterior mean under squared-error loss (Berger, 1985).

12

where we have made explicit the fact that the optimal information depends on the choice of prior
µ. Next we characterise Aµ , while an explicit example of optimal information for a Bayesian
PNM is detailed in Example 3.4.

3.4. Connection to Average Case Analysis
The decision theoretic framework in Section 3.1 is closely related to average-case analysis (ACA)
of classical numerical methods (Ritter, 2000). In ACA the performance of a classical numerical
method b : A → Q is studied in terms of the Bayes risk R(µ, M ) given in Eq. (3.1), for the PNM
M = (A, B) with belief operator B(µ, a) = δ ◦ b(a) as in Eq. (2.4). ACA is concerned with the
study of optimal information:


A∗µ ∈ arg inf inf R(µ, M ) s.t. M = (A, B), B = δ ◦ b .
A∈Λ

b

In general there is no reason to expect Aµ and A∗µ to coincide, since Bayesian PNM are not Bayes
rules6 . Indeed, an explicit example where Aµ 6= A∗µ is presented in Appendix S3. However, we
can establish sufficient conditions under which optimal information for a Bayesian PNM is the
same as optimal information for ACA:
Theorem 3.3. Let (Q, h·, ·iQ ) be an inner product space and the loss function L have the form
L(q † , q) = kq † − qk2Q where k · kQ is the norm induced by the inner product. Then the optimal
information Aµ for a Bayesian PNM and A∗µ for ACA are identical.
It is emphasised that this result is not a trivial consequence of the correspondance between
Bayes rules and worst case optimal methods, as exposed in Kadane and Wasilkowski (1983). To
the best of our knowledge, information-based complexity research has studied A∗µ but not Aµ .
Theorem 3.3 establishes that, for the squared norm loss, we can extract results on optimal
average case information from the ACA literature and use them to construct optimal Bayesian
PNMs. An example is provided next.
Example 3.4 (Optimal Information for Probabilistic Integration). To illustrate optimal information for Bayesian PNMs, we revisit the first worked example of ACA, due to Sul0 din (1959,
1960). Set X = {x ∈ C(0, 1) : x(0) = 0} and take the belief distribution µ to be induced
from the Weiner process on X , i.e. a Gaussian
process with mean 0 and covariance function
R1
k(t, t0 ) = min(t, t0 ). Our QoI is Q(x) = 0 x(t)dt and the loss function is L(q, q 0 ) = (q − q 0 )2 .
Consider standard information A(x) = (x(t1 ), . . . , x(tn )) for n fixed knots 0 ≤ t1 < · · · <
tn ≤ 1. Our aim is to determine knots ti that represent optimal information for a Bayesian
PNM with respect to µ and L.
Motivated by Theorem 3.3 we first solve the optimal information problem for ACA and
then derive
Pn the associated PNM. It will be sufficient to restrict attention to linear methods
b(a) = i=1 wi x(ti ) with wi ∈ R. This allows a closed-form expression for the average error:

 X
n
n
X
1
1 2
R(µ, (A, δ ◦ b)) = − 2
wi ti − ti +
wi wj min(ti , tj ).
3
2
i=1

(3.2)

i,j=1

Standard calculus can be used to minimise Eq. (3.2) over both the weights {wi }ni=1 and the
locations {ti }ni=1 ; the full calculation can be found in Chapter 2, Section 3.3 of Ritter (2000).
6

The distribution Q# µa will in general not be supported on the set of Bayes acts.

13

The result is an ACA optimal method
n

b(A(x)) =

2 X ∗
x(ti ),
2n + 1

t∗i =

i=1

2i
2n + 1

which is recognised as the trapezium rule with equally spaced knots. The associated contraction
rate rn is n−1 (Lee and Wasilkowski, 1986).
From Theorem 3.3 we have that ACA optimal information is also optimal information for the
Bayesian PNM. Thus the optimal Bayesian PNM M = (A, B) for the belief distribution µ is
uniquely determined:
 ∗ 
!
x(t1 )
n
2 X
1
 .. 
A(x) =  .  ,
B(µ, a) = N
ai ,
.
2n + 1
3(2n + 1)2
∗
i=1
x(tn )

Note how the PNM is centred on the ACA optimal method. However the PNM itself is not a
Bayes rule; it in fact carries twice the Bayes risk as the ACA method.
This illustration can be generalised. It is known that for µ induced from the Weiner process
on ∂ s x, Q a linear functional and φ a loss function that is convex and symmetric, equi-spaced
evaluation points are essentially optimal information, the Bayes rule is the natural spline of
degree 2s + 1, and the contraction rate rn is essentially n−(s+1) ; see Lee and Wasilkowski (1986)
for a complete treatment.
This completes our performance assessment for PNMs; next we turn to computational matters.

4. Numerical Disintegration
In this section we discuss algorithms to access the output from a Bayesian PNM. The approach
considered in this paper is to form an explicit approximation to µa that can be sampled. The
construction of a sampling scheme can exploit sophisticated Monte Carlo methods and allow
probing B(µ, a) at a computational cost that is de-coupled from the potentially substantial cost
of obtaining the information a itself.
The construction of an approximation to µa is non-trivial on a technical level. As shown in
Section 2.5, under weak conditions on the space X and the operator A, the disintegration µa
is well-defined for A# µ-almost all a ∈ A. The approach considered in this work is based on
sampling from an approximate distribution µaδ which converges in an appropriate sense to µa
in the δ ↓ 0 limit. This follows in a similar spirit to Ackerman et al. (2017).

4.1. Sequential Approximation of a Disintegration
Suppose that A is an open subset of Rn and that the distribution A# µ ∈ PA , admits a continuous and positive density pA with respect to Lebesgue measure on A. Further endow A with
the structure of a Hilbert space, with norm k · kA .
Let φ : R+ → R+ denote a decreasing function, to be specified, that is continuous at 0, with
φ(0) = 1 and limr→∞ φ(r) = 0. Consider


1
kA(x) − akA
µaδ (dx) := a φ
µ(dx)
Zδ
δ

14

where the normalisation constant
Zδa

:=

Z

φ



kã − akA
δ



pA (dã)

is non-zero since pA is bounded away from 0 on a neighbourhood of a ∈ A and φ is bounded
away from 0 on a sufficiently small interval [0, γ]. Our aim is to approximate µa with µaδ
for small bandwidth parameter δ. The construction, which can be considered a mathematical
generalisation of approximate Bayesian computation (Del Moral et al., 2012), ensures that
µaδ  µ. The role of φ is to admit states x ∈ X for which A(x) is close to a but not necessarily
equal. It is assumed to be sufficiently regular:
R
Assumption 4.1. There exists α > 0 such that Cφα := rα+n−1 φ(r)dr < ∞.
To discuss the convergence of µaδ to µa we must first select a metric on PX . Let F be a
normed space of (measurable) functions f : X → R with norm k·kF . For measures ν, ν 0 ∈ PX ,
define
dF (ν, ν 0 ) = sup |ν(f ) − ν 0 (f )|.
kf kF ≤1

This formulation encompasses many common probability metrics such as the total variation
distance and Wasserstein distance (Müller, 1997). However, not all spaces of functions F lead
0
to useful theory. In particular the total variation distance between µa and µa for a 6= a0 will be
one in general. Furthermore depending on the choice of F, dF may be merely a pseudometric7 .
Sufficient conditions for weak convergence with respect to F are now established:

Assumption 4.2. The map a 7→ µa is almost everywhere α-Hölder continuous in dF , i.e.
0

dF (µa , µa ) ≤ Cµα ka − a0 kαA
for some constant Cµα > 0 and for A# µ almost all a, a0 ∈ A.
Sufficient conditions for Assumption 4.2 are discussed in Ackerman et al. (2017), but are
somewhat technical.
Theorem 4.3. Let C̄φα := Cφα /Cφ0 . Then, for δ > 0 sufficiently small,
dF (µaδ , µa ) ≤ Cµα (1 + C̄φα )δ α
for A# µ almost all a ∈ A.

This result justifies the approximation of µa by µaδ when the QoI can be well-approximated by
integrals with respect to F. This result is stronger than that of earlier work, such as Pfanzagl
(1979), in that it holds for infinite-dimensional X , though it also relies upon the stronger Hölder
continuity assumption.
The specific form for φ is not fundamental, but can impact upon rate constants. For the
n
choice φ(r) = 1[r < 1] we have C̄φα = α+n
, which can be bounded independent of the dimension
n of A. On the other hand, for φ(r) = exp(− 12 r2 ) it can be shown that, for α ∈ N,
C̄φα =

(α + n − 1)!!
(n − 1)!!

(4.1)

so that the constant C̄φα might not be bounded. In general this necessitates effective Monte
Carlo methods that are able to sample from the regime where δ can be extremely small, in
order to control the overall approximation error.
7

For a pseudometric, dF (x, y) = 0 =⇒ x = y need not hold.

15

4.2. Computation for Series Priors
The series representation of µ in Eq. (2.6) of Section 2.6 is infinite-dimensional and thus cannot,
in general, be instantiated. To this end, define XN = x0 + span{φ0 , . . . , φN } and define the
associated projection operator PN : X → XN as
!
∞
N
X
X
PN x0 +
ui φi := x0 +
ui φi .
i=0

i=0

A natural approach is to compute with the modified information operator A ◦ PN instead of
A. This has the effect of updating the distribution of the first N + 1 coefficients and leaving
the tail unchanged, to produce an output µaδ,N . Then computation performed in the Bayesian
update step is finite-dimensional, whilst instantiation of the posterior itself remains infinitedimensional. A “likelihood-informed” choice of basis {φi } in such problems was considered in
Cui et al. (2016).
Inspired by this approach, we next considered convergence of the output µaδ,N to µaδ in the
limit N → ∞. In this section it is additionally required that φ be everywhere continuous with
φ > 0. Let ϕ = − log φ, so that ϕ is a continuous bijection of R+ to itself. The following are
also assumed:
Assumption 4.4. For each R > 0, it holds that |ϕ(r) − ϕ(r0 )| ≤ CR |r − r0 | for some constant
CR and all r, r0 < R.
Assumption 4.5. kA(x) − A ◦ PN (x)kA ≤ exp(m(kxkX ))Ψ(N ) for all x ∈ X , where m is
measurable and satisfies EX∼µ [exp(2m(kXkX ))] < ∞ and Ψ(N ) vanishes as N is increased.
Assumption 4.6. supx∈X kA(x)kA < ∞.
Assumption 4.7. kf k∞ ≤ CF kf kF for some constant CF and all f ∈ F.
Assumption 4.4 holds for the case ϕ(r) = 12 r2 with constant CR = R. Assumption 4.5 is
standard in the inverse problem literature; for instance it is shown to hold for certain series
priors in Theorem 3.4 of Cotter et al. (2010). Assumption 4.6 is, in essence, a compactness
assumption, in that it is implied by compactness of the state space X when A is linear. In
this sense it is a strong assumption; however it can be enforced in our experiments, where X is
unbounded, through a threshold map
(
A(x)
if kA(x)kA ≤ λmax ,
Ã(x) :=
A(x)
λmax kA(x)kA if kA(x)kA > λmax ,
where λmax is a large pre-defined constant. Assumption 4.7 places a restriction on the probability
metric dF in which our result is stated.
The following theorem has its proof in the Electronic Supplement:
Theorem 4.8. For some constant Cδ , dependent on δ, it holds that dF (µaδ,N , µaδ ) ≤ Cδ Ψ(N ).
An immediate consequence of Theorems 4.3 and 4.8 is that the total approximation error can
be bounded by applying the triangle inequality:
dF (µa , µaδ,N ) ≤ Cµα (1 + C̄φα )δ α + Cδ Ψ(N ).
In particular, we have convergence of µaδ,N to µa in the δ ↓ 0 limit provided that the number of
basis functions satisfies Cδ Ψ(N ) = o(1).

16

The approximate posterior µaδ,N analysed above can be sampledP
when µ is Gaussian, since the
first N + 1 coefficients can be handled with MCMC and the tail ∞
i=N +1 ui φi , being Gaussian,
can be sampled. However, when µ is non-Gaussian the tail is not recognised in a form that can
be sampled. For the experiments in Section 6, in which both Gaussian and non-Gaussian priors
µ are considered, the series in Eq. (2.6) was truncated at level N + 1, with the resultant prior
denoted µN . The associated posterior was then entirely supported on the finite-dimensional
subspace XN ; this is mathematically equivalent to working with the projected output PN µaδ,N .
Analysis of prior truncation, as opposed to modification of the information operator just reported, is known to be difficult. Indeed, while µN converges to µ weakly, it does not do so in
total variation, and this deficiency generally transfers to the associated posteriors. In general
the impact of prior perturbation is a subtle topic — see e.g. Owhadi et al. (2015) and the
references therein — and we therefore defer theoretical analysis of this approximation to future
work.

4.3. Monte Carlo Methods for Numerical Disintegration
The previous sections established a sequence of well-defined distributions µaδ (or µaδ,N for nonGaussian models) which converge (in a specific weak sense) to the exact disintegration µa .
From construction, µaδ  µa and this is sufficient to allow standard Monte Carlo methods to
be used. The construction of Monte Carlo methods is de-coupled from the core material in the
main text and the main methodological considerations are well-documented (e.g. Girolami and
Calderhead, 2011).
For the experiments reported in subsequent sections two approaches were explored; a Sequential Monte Carlo (SMC) method (Doucet et al., 2001) and a parallel tempering method (Geyer,
1991). This provided a transparent sampling scheme, whose non-asymptotic approximation
error can be theoretically understood. In particular, they provide robust estimators of model
evidence that can be used for Bayesian model comparison. Full details of the Monte-Carlo
methods used for this work, along with associated theoretical analysis for the SMC method, are
contained in Section S4.1 of the Electronic Supplement.

5. Computational Pipelines and PNM
The last theoretical development in this paper concerns composition of several PNMs. Most
analysis of numerical methods focuses on the error incurred by an individual method. However,
real-world computational procedures typically rely on the composition of several numerical
methods. The manner in which accumulated discretisation error affects computational output
may be highly non-trivial (Roy, 2010; Anderson, 2011; Babuška and Söderlind, 2016). An
extreme example occurs when one of the numerical methods in a pipeline is charged with
integration of a chaotic dynamical system (Strogatz, 2014).
In recent work, Chkrebtii et al. (2016), Conrad et al. (2016) and Cockayne et al. (2016) each
used PNMs within a broader statistical procedure to estimate unknown parameters in systems
of ODEs and PDEs. The probabilistic description of discretisation error was incorporated into
the data-likelihood, resulting in posterior distributions for parameters with inflated uncertainty
to properly account for the inferential impact of discretisation error. However, beyond these
limited works, no examination of the composition of PNMs has been performed. In particular,
the question of which PNMs can be composed, and when the output of such a composition
is meaningful, has not been addressed. This is important; for instance, if the output of a

17

composition of PNMs is to be used for analysis of variance to elucidate the main sources of
discretisation error, then it is important that such output is meaningful.
This section defines a pipeline as an abstract graphical object that may be combined with a
collection of compatible PNMs. It is proven that when compatible Bayesian PNMs are employed
in the pipeline, the distributional output of the pipeline carries a Bayesian interpretation under
an explicit conditional independence condition on the prior µ.
To build intuition, for the simple case where two Bayesian PNMs are composed in series,
our results provide conditions for when, informally, the output B2 (B1 (µ, a1 ), a2 ) corresponds to
a single Bayesian procedure B(µ, (a1 , a2 )). To reduce the notational and technical burden, in
this section we will not provide rigorous measure theoretic details; however we note that those
details broadly follow the same pattern as in Section 2.5.

5.1. Computational Pipelines
To analyse pipelines of PNMs, we consider n such methods M1 , . . . , Mn , where each method
Mi = (Ai , Bi ) is defined on a common8 state space X and targets a QoI Qi ∈ Qi . A pipeline
will be represented as a directed graphical model, wherein the QoIs Qi from parent methods constitute information operators for child methods. It may be that a method will take
quantities from multiple parents as input. To allow for this, we suppose that the information operator Ai : X → Ai can be decomposed into components Ai,j : X → Ai,j such that
Ai = (Ai,1 , . . . , Ai,m(i) ) and Ai = Ai,1 × · · · × Ai,m(i) . Thus, each component Ai,j can be
thought of as the QoI output by one of the parents of the method Mi .
Without loss of generality we designate the nth QoI Qn to be the principal QoI. That is, the
purpose of the computational pipeline is to estimate Qn . The case of multiple principal QoI is
a simple extension not described herein. Nodes with no immediate children are called terminal
nodes, while nodes with no immediate parents are called source nodes. We denote by A the set
of all source nodes.
Definition 5.1 (Pipeline). A pipeline P is a directed acyclic graph defined as follows:
• Nodes are of two kinds: Information nodes are depicted by , and method nodes are
depicted by .
• The graph is bipartite, so that edges connect a method node to an information node or
vice-versa. That is, edges are of the form  →  or  → .
• There are n method nodes, each with a unique label in {1, . . . , n}.
• The method node labelled i has m(i) parents and one child. Its in-edges are assigned a
unique label in {1, . . . , m(i)}.
• There is a unique terminal node and it is the child of method node n. This represents the
principal QoI Qn .
Example 5.2 (Distributed Integration). Recall the numerical integration problem of Example
3.4 and, as a thought experiment, consider partitioning the domain of integration in order to
distribute computation:
Z
Z
Z
1

8

|

0

0.5

x(t)dt =
{z } | 0
(c)

1

x(t)dt +
x(t)dt
0.5
{z
} | {z }

(a)

(5.1)

(b)

This is without loss of generality, since X can be taken as the union of all state spaces required by the individual
methods.

18

B1 (µ, ·)

x(t1 ), . . . , x(tm−1 )
x(tm )
x(tm+1 ), . . . , x(t2m )

B2 (µ, ·)

R 0.5
0

x(t)dt
B3 (µ, ·)

R1

0.5 x(t)dt

R1
0

x(t)dt

Figure 1: An intuitive representation of Example 5.2.
1

1

1

2
1
2

3
2

2

Figure 2: The pipeline P corresponding to Figure 1.
To keep presentation simple we consider an integral over [0, 1] with 2m + 1 equidistant knots
ti = i/2m. Let M1 be a Bayesian PNM for estimating Q1 (x) = (a) and M2 be a Bayesian PNM
for estimating Q2 (x) = (b).
In terms of our notational convention, we divide the information operator into four components; Ai,j , for i, j ∈ {1, 2}. A1,1 and A2,2 contain the information unique to M1 and M2 .
Specifically




x(t1 )
x(tm+1 )




..
..
A1,1 (x) = 
A2,2 (x) = 
,
.
.
.
x(t2m )

x(tm−1 )

A1,2 and A2,1 contain the information that is shared between the two methods; that is A1,2 =
A2,1 = {x(tm )}. To complete the specification we need a third PNM for estimation of Q3 (x) =
(c) which we denote M3 and which combines the outputs of M1 and M2 by simply adding them
together. Formally this has information operator A3 (x) = (A3,1 (x), A3,2 (x)) where A3,1 (x) =
(a) and A3,2 (x) = (b). Its belief update operator is given by:
B3 (µ, (a3,1 , a3,2 )) = δ(a3,1 + a3,2 )
An intuitive graphical representation of this set-up is shown in Figure 1. The pipeline P itself,
which is identical to Figure 1 but with additional node and edge labels, is shown in Figure 2.
In general, the method node labelled i is taken to represent the method Mi . The in-edge to
this node labelled j is taken to represent the information provided by the relationship Ai,j (x) =
ai,j . Here ai,j can either be deterministic information provided to the pipeline, or statistical
information derived from the output of another PNM. To make this formal and to “match
the input-output spaces” we next define what it means for the collection of methods Mi to be
compatible with the pipeline P . Informally, this describes the conditions that must be satisfied
for method nodes in a pipeline to be able to connect to each other.
Definition 5.3 (Compatible). The collection (M1 , . . . , Mn ) of PNMs is compatible with the
pipeline P if the following two requirements are satisfied:

19

(i) (Method nodes which share an information node must have consistent information spaces
and information operators.) For a motif
i

j0

i0

j

we have that Ai,i0 = Aj,j 0 and Ai,i0 = Aj,j 0 .
(ii) (The space Qi for the output of a previous method must be consistent with the information
space of the next method.) For a motif
j0

i

j

we have that Qi = Aj,j 0 .
Note that we do not require the converse of (i) at this stage; that is, the same information can
be represented by more than one node in the pipeline. This permits redundancy in the pipeline,
in that information is not recycled. It will transpire that pipelines with such redundancy are
non-Bayesian.
The role of the pipeline P is to specify the order in which information, either deterministic
of statistical, is propagated through the collection of PNMs. This is illustrated next:
Example 5.4 (Propagation of Information). For the pipeline in Figure 2, the propagation of
information proceeds as follows::
1. The source nodes, representing A(x) = {A1,1 (x), A1,2 (x) = A2,1 (x), A2,2 (x)} are evaluated
as {a1,1 , a1,2 = a2,1 , a2,2 }. This represents all the information on x that is provided.
2. The distributions
µ(1) := B1 (µ, (a1,1 , a1,2 ))
µ(2) := B2 (µ, (a2,1 , a2,2 ))
are computed.
3. The push-forward distribution
µ(3) := (B3 )# (µ, µ(1) × µ(2) )
is computed.
Here µ(1) × µ(2) is defined on the Cartesian product ΣA3,1 × ΣA3,2 with independent components
µ(1) and µ(2) . The notation (B3 )# refers to the push-forward of the function B3 (µ, ·) over its
second argument. The distribution µ(3) is the output of the pipeline and is a distribution over
the principal QoI Q3 (x).
The procedure in Example 5.4 can be formalised, but to keep the presentation and notation
succinct, we leave this implicit:

20

1

4

2
3

6
5

Figure 3: Dependence graph G(P ) corresponding to the pipeline P in Figure 2. The nodes are
indexed with a topological ordering (shown).
Definition 5.5 (Computation). For a collection (M1 , . . . , Mn ) of PNMs that are compatible
with a pipeline P , the computation P (M1 , . . . , Mn ) is defined as the PNM with information
operator A and belief update operator B that takes µ and A(x) = a as input and returns the
distribution µ(n) as its output B(µ, a), obtained through the procedure outlined in Example
5.4.
That is, the computation P (M1 , . . . , Mn ) is a PNM for the principal QoI Qn . Note that this definition includes a classical numerical work-flow just as a PNM encompasses a standard numerical
method.

5.2. Bayesian Computational Pipelines
Noting that P (M1 , . . . , Mn ) is itself a PNM, there is a natural definition for when such a
computation can be called Bayesian:
Definition 5.6 (Bayesian Computation). Denote by (A, B) the information and belief operators
associated with the computation P (M1 , . . . , Mn ) and let {µa }a∈A be a disintegration of µ with
respect to the information operator A. The computation P (M1 , . . . , Mn ) is said to be Bayesian
for the QoI Qn if
B(µ, a) = (Qn )# µa for A# µ-almost-all a ∈ A.
This is clearly an appealing property; the output of a Bayesian computation can be interpreted
as a posterior distribution over the QoI Qn (x) given the prior µ and the information A(x). Or,
more informally, the “pipeline is lossless with information”. However, at face value it seems
difficult to verify whether a given computation P (M1 , . . . , Mn ) is Bayesian, since it depends
on both the individual PNMs Mi and the pipeline P that combines them. Our next aim is to
establish verifiable sufficient conditions, for which we require another definition:
Definition 5.7 (Dependence Graph). The dependence graph of a pipeline P is the directed
acyclic graph G(P ) obtained by taking the pipeline P , removing the method nodes and replacing
all  →  →  motifs with direct edges  → .
The dependency graph for Example 5.2 is shown in Figure 3.
For a computation P (M1 , . . . , Mn ), each of the J distinct nodes in G(P ) can be associated
with a random variable Yj where either Yj = Ak,l (X) for some k, l, when the node is a source,
or otherwise Yj = Qk (X), for some k. Randomness here is understood to be due to X ∼ µ, so
that the distribution of the {Yj }Jj=1 is a function of µ. The convention used here is that the Yj
are indexed according to a topological ordering on G(P ), which has the properties that (i) the
source nodes correspond to indices 1, . . . , I, and (ii) the final random variable is YJ = Qn (X).

21

Definition 5.8 (Coherence). Consider a computation P (M1 , . . . , Mn ). Denote by π(j) ⊆
{1, . . . , j − 1} the parent set of node j in the dependence graph G(P ). Then we say that
µ ∈ PX is coherent for the computation P (M1 , . . . , Mn ) if the implied joint distribution of the
random variables Y1 , . . . , YJ satisfies:
Yj ⊥
⊥ Y{1,...,j−1}\π(j) | Yπ(j)
for all j = I + 1, . . . , J.
Note that this is weaker than the Markov condition for directed acyclic graphs (see Lauritzen,
1991), since we do not insist that the variables represented by the source nodes are independent.
It is emphasised that, for a given µ ∈ PX , the coherence condition can in general be checked
and verified.
The following result provides sufficient and verifiable conditions which ensure that a computation composed of individual Bayesian PNMs is a Bayesian computation:
Theorem 5.9. Let M1 , . . . , Mn be Bayesian PNMs and let µ ∈ PX be coherent for the computation P (M1 , . . . , Mn ). Then it holds that the computation P (M1 , . . . , Mn ) is Bayesian for the
QoI Qn .
Conversely, if non-Bayesian PNM are combined then the computation P (M1 , . . . , Mn ) need
not be Bayesian in general.
Example 5.10 (Example 5.2, continued). The random variables Yi in this example are:
Z 0.5
Z 1
m−1
2m
Y1 = {X(ti )}i=1 , Y2 = X(tm ), Y3 = {X(ti )}i=m+1 , Y4 =
X(t)dt, Y5 =
X(t)dt.
0

0.5

From G(P ) in Figure 3, coherence condition in Definition 5.8 requires that the non-trivial
conditional independences Y4 ⊥
⊥ Y3 | {Y1 , Y2 } and Y5 ⊥⊥ Y1 | {Y2 , Y3 } hold. Thus the distribution µ is coherent for the computation
if and only if, for X ∼R µ, the assoR 0.5 P (M1 , M2 , M3 ) 2m
1
⊥
ciated information variables satisfy 0 X(t)dt ⊥⊥ {X(ti )}i=m+1 |{X(ti )}m
i=1 and 0.5 X(t)dt ⊥
2m .
{X(ti )}m−1
|{X(t
)}
i i=m
i=1
The distribution µ induced by the Weiner process on x in Example 3.4 satisfies these conditions. Indeed, under µ the stochastic process {x(t) : t > tm } is conditionally independent of its
history {x(t) : t < tm } given the current state x(tm ). Thus for this choice of µ, from Theorem
5.9 we have that P (M1 , M2 , M3 ) is Bayesian and parallel computation of (a) and (b) in Eq. (5.1)
can be justified from a Bayesian statistical standpoint.
However, for the alternative of belief distributions induced by the Weiner process on ∂ s x, this
condition is not satisfied and the computation P (M1 , M2 , M3 ) is not Bayesian. To turn this
into a Bayesian procedure for these alternative belief distributions it would be required that
A1,2 (x) provides information about the derivatives ∂ k x(tm ) for all orders k ≤ s.

5.3. Monte Carlo Methods for Probabilistic Computation
The most direct approach to access µ(n) is to sample from each Bayesian PNM and treat the
output samples as inputs to subsequent PNM. This is sometimes known as ancestral sampling in
the Bayesian network literature (e.g. Paige and Wood, 2015), and is illustrated in the following
example:
Example 5.11 (Ancestral Sampling for PNM). For Example 5.2, ancestral sampling proceeds
as follows:

22

1. Draw initial samples
q1 ∼ B1 (µ, (a1,1 , a1,2 ))
q2 ∼ B2 (µ, (a2,1 , a2,2 ))

2. Draw a final sample
q3 ∼ B3 (µ, (q1 , q2 ))
Then q3 is a draw from µ(3) .
Ancestral sampling requires that PNM outputs can be sampled. Such sampling methods were
discussed in Section 4.3. For a more general approach, sequential Monte Carlo methods can be
used to propagate a collection of particles through the pipeline P , similar to work on SMC for
general graphical models (Briers et al., 2005; Ihler and McAllester, 2009; Lienart et al., 2015;
Lindsten et al., 2017; Paige and Wood, 2015).

6. Numerical Experiments
In this final section of the paper we present three numerical experiments. The first is a linear
PDE, the second is a nonlinear ODE and the third is an application to a problem in industrial
process monitoring, described by a pipeline of PNM. In each case we experiment with nonGaussian belief distributions and, in doing so, go beyond previous work.

6.1. Poisson Equation
Our first illustration is an instance of the Poisson equation, a linear PDE with mixed DirichletNeumann boundary conditions:
−∇2 x(t) = 0

x(t) = t1

x(t) = 1 − t1

∂x/∂t2 = 0

t ∈ (0, 1)2

(6.1)

t1 ∈ [0, 1]

t2 = 0

(6.2)

t2 = 1

(6.3)

t2 ∈ (0, 1)

t1 = 0, 1

(6.4)

t1 ∈ [0, 1]

A model solution to this system, generated with a finite-element method on a fine mesh, is
shown in Figure 4.
As the spatial domain for this problem is two-dimensional, the basis used for specification of
the belief distribution is more complex. Here tensor products of orthogonal polynomials have
been used: φi (t) = Cj (2t1 − 1)Ck (2t2 − 1), i + j ≤ Nc . The polynomials Ci were chosen to
be normalised Chebyshev polynomials of the first kind. Prior specification then follows the
formulation given in Section 2.6, where the remaining parameters were chosen to be x0 ≡ 1,
and γi = α(i + 1)−2 . The random variables ξ were taken to be either Gaussian or Cauchy and
the polynomial basis was truncated to N = 45 terms, corresponding to a maximum polynomial
degree of NC = 8. For both priors the parameter α was set to α = 1. Note that closed-form
expressions are available for analysis under the Gaussian prior (Cockayne et al., 2016) but, to
simplify interpretation of empirical results, were not exploited. Mathematical background on
Cauchy priors can be found in Sullivan (2016).
The information operator was defined by a set of locations ti ∈ [0, 1]2 , i = 0, . . . , Nt , where

either the interior condition or one of the boundary conditions was enforced. Denote by tI,i
23

1.0
0.8

t2

0.6
0.4
0.2
0.0
0.0

0.2

0.4

t1

0.6

0.8

1.0

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

Figure 4: Model solution x(t), t = (t1 , t2 ), generated by application of a finite element method
based on a triangular mesh of 50 × 50 elements.


the set of interior points, tD,j the set of Dirichlet boundary points and tN,k the set of
Neumann boundary points, where i = 1, . . . , NI , j = 1, . . . , ND and k = 1, . . . , NN , with
n = NI + ND + NN . Then, the information operator is given by the concatenation of the
conditions defined above:
A(x) = [AI (x)> , AD (x)> , AN (x)> ]> ,





 ∂
N,1 )
−∇2 x(tI,1 )
x(tD,1 )
∂t1 x(t






..
..
..
D
N
AI (x) = 
 , A (x) = 
 , A (x) = 

.
.
.
2
I,N
D,N
∂
N,N
I
D
N
−∇ x(t
)
x(t
)
)
∂t1 x(t

The Bayesian PNM output was approximated by numerical disintegration and sampled with
a Monte Carlo method whose description is reserved for the Electronic Supplement. In Figure 5
the mean and pointwise standard-deviations of the posterior distributions are plotted for Gaussian and Cauchy priors with n = 16. There is little qualitative difference between the posterior
distributions for the Gaussian and Cauchy priors. The mean functions match closely to the
mean function from the model solution, as given in Figure 4. The posterior variance is lowest
near to the Dirichlet boundaries where the solution is known, and peaks where the Neumann
condition is imposed. This is to be expected, as evaluations of the Neumann boundary condition
provide less information about the solution itself.
Next, the posterior distribution of the spectrum {ui } was investigated. In Figure 6 the
posterior distribution over these coefficients is plotted and it is seen that the correlation structure
between coefficients is non-trivial, c.f. the joint distribution between u0 and u3 .
Last, in Figure 7 convergence of the posterior distribution is plotted as the number of design
points is varied, for n = 16, 25, 36. In each case a Gaussian prior was used. As expected, the
standard deviation in the posterior distribution is seen to decrease as the number of design
points is increased. At n = 36, the shape of the region of highest uncertainty changes markedly,
with the most uncertain region lying between the Dirichlet boundary and the first evaluation
points on the Neumann boundary. This is likely due to the fact that the number of evaluation
points is approaching the size of the polynomial basis; when the number of points equals the
size of the basis the system is completely determined for a linear model. Thus, we need N  n
in order for discretisation error to be quantified.

24

0.8

t2

0.6
0.4
0.2
0.0
0.0

0.2

0.4

t1

0.6

0.8

1.0

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.8
0.6

t2

1.0

0.4
0.2
0.0
0.0

0.2

0.4

t1

0.6

0.8

0.0250
0.0225
0.0200
0.0175
0.0150
0.0125
0.0100
0.0075
0.0050
0.0025
0.0000

(a) Gaussian prior
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0.8

t2

0.6
0.4
0.2
0.0
0.0

0.2

0.4

t1

0.6

0.8

0.8
0.6

t2

1.0

0.4
0.2
0.0
0.0

1.0

0.2

0.4

t1

0.6

0.8

0.0250
0.0225
0.0200
0.0175
0.0150
0.0125
0.0100
0.0075
0.0050
0.0025
0.0000

(b) Cauchy prior

Figure 5: Posterior distributions for the solution x of the Poisson equation, with n = 16 and
different choices of prior distribution. Left: Posterior mean. Design points for the
interior, Dirichlet and Neumann boundary conditions are indicated by green dots,
green squares and green crosses, respectively. Right: Posterior standard deviation.

25

u0

8.5
6.2
3.8
1.5

u1

0.0
0.0
-0.0
-0.0

u2

0.0
0.0
-0.0
-0.0

u3

0.0
0.0
-0.0
-0.0

u4

-0.8
-0.8
-0.8
-0.8

u5

0.0
0.0
-0.0
-0.0

1.5 1.6 1.7 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.1 -0.8 -0.8 -0.8 -0.0 0.0 0.0
Figure 6: Posterior distributions for the first six coefficients of the spectrum for the solution
x of the Poisson equation, obtained with Monte Carlo methods and numerical disintegration, based on δ = 0.0008, n = 16. (NB: The posterior is Gaussian and can
be obtained in closed-form, but we opted to additionally illustrate the Monte Carlo
method.)

0.4
0.2
0.0
0.0

0.2

0.4

t1

0.6

0.8

(a) n = 16

0.6
0.4
0.2
0.0
0.0

0.2

0.4

t1

0.6

0.8

(b) n = 25

0.0250
0.0225
0.0200
0.0175
0.0150
0.0125
0.0100
0.0075
0.0050
0.0025
0.0000

0.8
0.6

t2

t2

0.6

0.8

t2

0.8

0.0250
0.0225
0.0200
0.0175
0.0150
0.0125
0.0100
0.0075
0.0050
0.0025
0.0000

0.4
0.2
0.0
0.0

0.2

0.4

t1

0.6

0.8

0.0250
0.0225
0.0200
0.0175
0.0150
0.0125
0.0100
0.0075
0.0050
0.0025
0.0000

(c) n = 36

Figure 7: Heat map of the point-wise standard deviation for the solution x to the Poisson
equation as the number n of design points is varied. In each case a Gaussian prior
has been used.

26

4

10 -2
10 -4

1

un

x(t)

2

0

10 -6
10 -8

1

10 -10

2

10 -12

3

Positive
Negative

10 0

3

0

2

4

t

6

8

10

10 -14

0

10

20

30

40

i

50

60

70

80

Figure 8: Two distinct solutions for the Painlevé ODE. The spectral plot on the right shows the
true coefficients {ui }, as determined by a model solver (the MatLab package chebfun).

6.2. The Painlevé ODE
In this section a Bayesian PNM is developed to solve a nonlinear ODE based on Painlevé’s first
transcendental
x00 = x2 − t,

t−1/2 x(t) → 1

t ∈ [0, ∞)x(0)

=0

as t → ∞ .

To permit computation, the right-boundary condition
was relaxed by truncating the domain to
√
[0, 10] and using the modified condition x(10) = 10.
Two distinct solutions are known, illustrated in Figure 8 (left). These model solutions were
obtained using the deflation technique described in Farrell et al. (2015). The spectrum plot
in Figure 8 (right) represents the coefficients {ui } obtained when each solution is represented
over a basis of normalised Chebyshev polynomials. As those polynomials are orthonormal with
respect to the L2 -inner-product, the slower decay for the negative solution compared to the
positive solution is equivalent to the negative solution having a larger L2 -norm. This explains
the preference that optimisation-based numerical solvers have for returning the positive solution
in general, and also explains some of the results now presented.
Such systems for which multiple solutions exist have been studied before in the context of
PNM, both in Chkrebtii et al. (2016) and in Cockayne et al. (2016). It was noted in both papers
that existence of multiple solutions can present a substantial challenge to classical numerical
methods.
To build a Bayesian PNM, a prior µ for this problem was defined by using a series expansion
as in Eq. (2.6). The basis functions were φi (t) = Ci ( 21 (t − 5)) where the Ci were normalised
Chebyshev polynomials of the first kind. Both Gaussian and Cauchy priors were considered
by taking ui := γi ξi , where ξi were taken to be either standard Gaussian or standard Cauchy
and in in each case x0 (t) ≡ 0. In accordance with the exponential convergence rate for spectral
methods when the solution to the system is a smooth function, the sequence of scale parameters
was set to γi = αβ −i , where α = 8 and β = 1.5. These values were chosen by inspection of the
true spectra (obtained with Matlab’s “chebfun” package) to ensure that both solutions were in
the support of the prior.

27

The information operator A was defined by the choice of locations {tj }, j = 1, . . . , m, which
determine the locations at which the posterior will be constrained. Analysis for several values
of m was performed. In each case t1 = 0, tm = 10 and the remaining tj were equally spaced on
[0, 10]. To be explicit, the information operator was
 00

x (t1 ) − (x(t1 ))2


..


.
 00

A(x) = x (tm ) − (x(tm ))2 




x(0)
x(10)

with the last two√elements enforcing the boundary conditions. Thus our information was a =
[−t1 , . . . , −tm , 0, 10], which is n = m + 2 dimensional.
The Bayesian PNM output B(µ, a) was approximated via numerical disintegration with the
first N = 40 terms of the series representation used. This was sampled with Monte Carlo
methods, the details of which are reserved for the Electronic Supplement.
Results for a selection of bandwidths δ, with n = 17, are shown in Figure 9. Note that a
strong preference for the positive solution is expressed at the smallest δ, with mass around both
solutions at larger δ. For the Gaussian prior, some mass remained around the negative solution
at the smallest δ, while this was not so for the Cauchy prior. This reflects the fact that, for
a collection of independent univariate Cauchy random variables, one element is likely to be
significantly larger in magnitude than the others, which favours faster decay for the remaining
elements.
Using the calculation described in Section S4.4, model evidence was computed for both the
Gaussian and the Cauchy prior at n = 15. The Bayes factor for the Cauchy, compared to the
Gaussian prior, was found to be 20.26, which constitutes strong evidence in favour of a Cauchy
prior for this problem at the given level of discretisation.
In Figure 10 the posterior distributions for first six coefficients ui at n = 17 and δ = 1 are
plotted. Strong multimodality is clear, as well as skewed correlation structure between the
coefficients. Illustration of such posteriors for smaller δ is difficult as the posteriors become
extremely peaked.
Figure 11 displays convergence of the posterior distributions as n is increased. Of particular
interest is that for n = 12, the posterior distribution based on a Gaussian prior becomes
trimodal. For each prior, the posterior mass settles on the positive solution to the system at
n = 22. This is in accordance with the fact that this solution has smaller L2 -norm. This perhaps
reflects the fact that, while in the limiting case both solutions should have an equal likelihood,
the curvature of the likelihood at each mode may differ. Prior truncation may also be influential;
in Figure 12 the log-likelihood of the negative solution increases at a slower rate than that of
the positive solution. Thus, while in the setting of an infinite prior series neither solution should
be preferred, in practice truncation might bias one solution over the other. Lastly, it is clear
that the parameters α and β may also have a significant effect on which solution is preferred.
Further theoretical work will be required to understand many of the phenomena that we have
just described.
Of particular interest is how a preference for the negative solution could be encoded into a
PNM. Owing to the flexible specification the information operator, there is considerable choice in
this matter. An elegant approach is the introduction of additional, inequality-based information
x0 (0) ≤ 0 .

28

(6.5)

(a) Gaussian Prior

(b) Cauchy Prior.

Figure 9: Posterior samples for the Painlevé system for n = 17. Blue and green dashed lines
represent the positive and negative solutions determined with chebfun. Grey lines are
samples from an approximation to the posterior provided by numerical disintegration
(bandwidth parameter δ).

29

0.8
0.6
0.4
0.1
3.5
2.8
2.1
1.4
-0.3
-0.8
-1.2
-1.7

u0

u1

u2

u3

0.7
0.2
-0.2
-0.7

u4

1.6
0.9
0.2
-0.6

u5

0.7
0.2
-0.2
-0.7

2.3 3.0 3.7 1.4 2.5 3.5 -1.7 -1.0 -0.3 -0.7 0.0 0.7 -0.6 0.5 1.6 -0.7 0.0 0.7
Figure 10: Posterior distributions for the first six coefficients obtained with numerical disintegration (bandwidth parameter δ = 1), at n = 17. Vertical dashed lines on the diagonal
plots indicate the value of the coefficients for the positive (blue) and negative (green)
solutions determined with chebfun.

30

Figure 11: Convergence for the numerical disintegration scheme as n is increased. Left: Gaussian prior. Right: Cauchy prior. In all cases δ = 10−4 .

31

Positive
Negative

108

106

− log φ



kA(x)−ak
δ



1010

104

102

0

10

20

30

40

50

60

70

N

Figure 12: Negative-log-likelihoods for the point-estimates of coefficients for the postive and
negative solutions given by chebfun, as the truncation level N is varied. The fact
that the likelihood for the positive solution decreases more rapidly than that of the
negative solution suggests indicates that the posterior may have a preference for that
solution over the other, though the level N = 40 has been selected in an attempt to
minimise the impact.
Such information can be difficult to incorporate in standard numerical algorithms, but is of
interest in many physical problems (Kinderlehrer and Stampacchia, 2000). For Bayesian PNM
we can extend the information operator to include 1[x0 (0) ≤ 0]. Posterior distributions for the
Gaussian prior at n = 17 are shown in Figure 13. Note that posterior mass has settled close to
the negative solution. This highlights the simplicity with which Bayesian PNMs can encode a
preference for a particular solution when a multiplicity of solutions exist.

6.3. Application to Industrial Process Monitoring
This final application illustrates how statistical models for discretisation error can be propagated
through a pipeline of computation to model how these errors are accumulated.
Hydrocyclones are machines used to separate solid particles from a liquid in which they are
suspended, or two liquids of different densities, using centrifugal forces. High pressure fluid is
injected into the top of a tank to create a vortex. The induced centrifugal force causes denser
material to move to the wall of the tank while lighter material concentrates in the centre,
where it can be extracted. They have widespread applications, including in areas such as
environmental engineering and the petrochemical industry (Sripriya et al., 2007). An illustration
of the operation is given in Figure 14.
To ensure the materials are well-separated the hydrocyclone must be moitored to allow adjustment of the input flow-rate. This is also important for safe operation, owing to the high pressures
involved (Bradley, 2013). However, direct monitoring is impossible owing to the opaque walls of
the equipment and the high interior pressure. For this purpose electrical impedance tomography
(EIT) has been proposed to allow monitoring of the contents (Gutierrez et al., 2000).
EIT is a technique which allows recovery of an interior conductivity field based upon measurements of voltage obtained from applying a stimulating current on the boundary. It is suited
to this problem, as the two materials in the hydrocyclone will generally be of different conductivities. In its simplified form due to Calderón (1980), EIT is described by a linear partial

32

Figure 13: Posterior distribution at n = 17, based on a Gaussian prior, with the negative
boundary condition given by Eqn. (6.5) enforced. Left: δ = 0.99. Right: δ = 0.0001.

overflow
input flow
less dense

more dense
underflow

(a) Hydrocyclone tank schematic

(b) Cross-section (top of tank)

Figure 14: A schematic description of hydrocyclone equipment. (a) The tank is cone-shaped
with overflow and underflow pipes positioned to extract the separated contents. (b)
Fluid, a mixture to be separated, is injected at high pressure at the top of the tank
to create a vortex. Under correct operation, denser materials are directed toward
the centre of the tank and less-dense materials are forced to the peripheries of the
tank.

33

differential equation similar to that in Section 6.1, but with modified boundary conditions to
incorporate the stimulating currents and measured voltages:
−∇ · (a(t)∇x(t)) = 0

∂x
ce
a(t) (t) =
0
∂n

t∈D
t = te
e
t ∈ ∂D \ {te }N
e=1

(6.6)

where D denotes the domain, modelling the hydrocyclone tank, e indexes the stimulating electrodes, te ∈ ∂D are the corresponding locations of the electrodes on ∂D, a is the unknown
∂
conductivity field to be determined and ∂n
denotes the derivative with respect to the outward
1
pointing normal vector. The electrode t is referred to as the reference electrode. The vector
c = (c1 , . . . , cNe ) denotes the stimulation current pattern. Several stimulation patterns were
considered, denoted cj , j = 1, . . . , Nj .
The experimental data described in West et al. (2005) were considered. In the experiment, a
cylindrical perspex tank was used with a single ring of eight electrodes. Translation invariance
in the vertical direction means that the contents are effectively a single 2D region and electrical
conductivity can be modelled as a 2D field. At the start of the experiment, a mixing impeller was
used to create a rotational flow. This was then removed and, after a few seconds, concentrated
potassium chloride solution was carefully injected into the tap water initially filling the tank.
Data, denoted yτ , were collected at regular time intervals by application of several stimulation
patterns c1 , . . . , cM .
To formulate the statistical problem, consider parameterising the conductivity field as a(τ, t),
where τ ∈ [0, T ] is a temporal index while t ∈ D is the spatial coordinate and D is the circular
domain representing the perspex tank in the experiment. A log-Gaussian prior was placed over
the conductivity field so that log ais a Gaussian
process with separable covariance function

kt−t0 k2
0
0
0
ka ((τ, t), (τ , t )) := λ min(τ, τ ) exp − 2`2
where ` is a length-scale parameter representing
the anticipated spatial variation of the conductivity field and λ is a parameter controlling the
amplitude of the field. Here ` was fixed to ` = 0.3, while λ = 10−3 . The problem of estimating
a based on data can be well-posed in the Bayesian framework (Dunlop and Stuart, 2016). Full
details of this experiment can be found in the accompanying report Oates et al. (2017).
Our aim is to use a PNM to account for the effect of discretisation on inferences that are
made on the conductivity

field. For fixed τ , a Gaussian prior was posited for x, with covariance
kt−t0 k2
0
kx (t, t ) := exp − 2`2
where `x was fixed to `x = 0.3. The associated Bayesian PNM, a
x
probabilistic meshless method (PMM), was described in Example 2.4.
The statistical inference procedure is formulated in a pipeline of computations in Figure 15.
It is assumed that the desired outcome is to monitor the contents of the tank while the current
contents are being mixed. This suggests a particle filter approach where a PMM Mτ is employed
to handle the intractable likelihood p(yτ |aτ ) that involves the exact solution of a PDE. The
distribution of aτ given y1 , . . . , yτ is denoted πτ an the computation P (M1 , . . . , Mτ ) is Bayesian
only if the particle approximation error due to the use of a particle filter is overlooked.
To briefly illustrate the method, Figure 16 presents posterior means for the field a(τ, ·), for
each post-injection time point τ = 1, . . . , 8. These are based on a particle approximation of size
P = 500, with method nodes based upon a Bayesian PNM, as in Example 2.4, with n = 119
design points. The high conductivity region representing the potassium chloride solution can be
seen rotating through the domain in the frames after injection, with its conductivity reducing
as it mixes with the water. The full posterior distribution over the conductivity field is inflated
as a result of explicitly modelling the discretisation error; an extensive analysis of these results
will be reported in the upcoming Oates et al. (2017).

34

yτ
2
1

...

...

τ

dπτ −1
dπ0

dπτ
dπ0

Figure 15: Pipeline for hydrocyclone application: The method node (black) represents the use
of PMM solvers, which are incorporated into the likelihood for evolving the particles
according to a Markov transition kernel.

τ =1

1.0

τ =2

1.0

0.5

τ =3

1.0

0.5

τ =4

1.0

0.5

0.5
138.1

0.0

0.0

0.0

0.0

123.3

−0.5

−0.5

−0.5

−0.5

108.4
93.6

−1.0
−1.0

−0.5

0.0

0.5

1.0

τ =5

1.0

−1.0
−1.0

−0.5

0.0

0.5

1.0

τ =6

1.0

−1.0
−1.0

−0.5

0.0

0.5

1.0

τ =7

1.0

−1.0
−1.0

−0.5

0.5

0.5

0.5

0.0

0.0

0.0

0.0

−0.5

−0.5

−0.5

−0.5

0.5

1.0

τ =8

1.0

0.5

0.0

78.7
63.9
49.0
34.1
19.3

−1.0
−1.0

−0.5

0.0

0.5

1.0

−1.0
−1.0

−0.5

0.0

0.5

1.0

−1.0
−1.0

−0.5

0.0

0.5

1.0

−1.0
−1.0

4.4

−0.5

0.0

0.5

1.0

Figure 16: Mean conductivity fields recovered in the hydrocyclone experiment, for the first 8
frames post-injection.

35

9

Pipeline
Static

30

8

Pipeline - Static

28

24
22

R

D

σ(t) dt

26

20
18

7
6
5
4
3

16

2

14
1

2

3

4

5

6

7

8

1

τ

2

3

4

5

6

7

8

τ

Figure 17: Left: Integrated standard-deviation over the domain, for the first 8 frames postinjection, for both the pipeline and the static approaches described in the text.
Right: The difference between these two quantities.
R
In Figure 17, the integrated standard-deviation D σ(t) dt is shown for τ = 1, . . . , 8 for both
the “pipeline”, as described above, and a “static” approach in which no uncertainty was propagated. In this static approach a symmetric collocation PDE solver9 was used to solve the
forward problem, and a separate Bayesian inversion problem was solved at each time point.
The parameters of the symmetric collocation solver were identical to those used in the PMM.
In the left panel we observe some structural periodicity, present in both the pipeline and the
static approach. We speculate that this may be due to the rotation of the medium causing
the area of high conductivity to periodically reach an area of the domain, relative to the 8
sensors, in which it is particularly easy to recover. With this periodicity subtracted in the
right panel, there was a clear increase in posterior uncertainty in the pipeline compared to
the static approach, which is depicted. Temporal regularisation would usually be expected to
reduce uncertainty; thus, the fact that the overall uncertainty increased with τ , relative to the
static formulation, demonstrates that we have quantified and propagated uncertainty due to
successive discretisation of the PDE at each time point.

7. Discussion
This paper has established statistical foundations for PNMs and investigated the Bayesian case
in detail. Through connection to Bayesian inverse problems (Stuart, 2010), we have established
when Bayesian PNM can be well-defined and when the output can be considered meaningful.
The presentation touched on several important issues and a brief discussion of the most salient
points is now provided.
Bayesian vs Non-Bayesian PNMs The decision to focus on Bayesian PNMs was motivated by
the observation that the output of a pipeline of PNMs can only be guaranteed to admit a valid
Bayesian interpretation if the constituent PNMs are each Bayesian and the prior distribution is
coherent. Indeed, Theorem 5.9 demonstrated that prior coherence can be established at a local
level, essentially via a local Markov condition, so that Bayesian PNMs provide a extensible
modelling framework as required to solve more challenging numerical tasks. These results
support a research strategy that focuses on Bayesian PNMs, so that error can be propagated
in a manner that is meaningful.
9

Recall that the PMM has a corresponding symmetric collocation solution to the PDE as its mean function.

36

On the other hand, there are pragmatic reasons why either approximations to Bayesian
PNMs, or indeed, non-Bayesian PNMs might be useful. The predominant reason would be to
circumvent the off-line computational costs that can be associated with Bayesian PNMs, such
as the use of numerical disintegration developed in this research. Recent research efforts, such
as Schober et al. (2014, 2016) and Kersting and Hennig (2016) for the solution of ODEs, have
aimed for computational costs that are competitive with classical methods, at the expense of
fully Bayesian estimation for the solution of the ODE. Such methods are of interest as nonBayesian PNMs, but their role in pipelines of PNMs is unclear. Our contribution serves to
make this explicit.
Computational Cost The present research focused on the more fundamental cost of access to
the information A(x), rather than the additional CPU time required to obtain the PNM output.
Indeed, numerical disintegration constituted the predominant computational cost in the applications that were reported. However, we stress that in many challenging applications gated by
discretisation error, such as occur with climate models, the fundamental cost of the information
A(x) will be dominant. Furthermore, the Monte Carlo methods that were employed for numerical disintegration admit substantial improvements (e.g. in a similar vein to Botev and Kroese,
2012; Koskela et al., 2016). The objective of this paper was to establish statistical foundations
that will permit the development of more sophisticated and efficient Bayesian PNMs.
Prior Elicitation Throughout this work we assumed that a belief distribution µ was provided.
The question of whose belief is represented in µ has been discussed by several authors and
a chronology is included in the Electronic Supplement. Of these perspectives we mention in
particular Hennig et al. (2015), wherein µ is the belief of an agent that “we get to design”.
This offers a connection to frequentist statistics, in that an agent can be designed to ensure
favourable frequentist properties hold.
A robust statistics perspective is also relevant and one such approach would be to consider a
generalised Bayes risk (Eq. (3.1)) wherein the state variable X used for assessment is assumed
to be drawn from a distribution µ̃ 6= µ. This offers an opportunity to derive Bayesian PNMs
that are robust to certain forms of prior mis-specification. This direction was not considered in
the present paper, but has been pursued in the ACA literature for classical numerical methods
(see Chapter IV, Section 4 of Ritter, 2000).
In general, the specification of prior distributions for robust inference on an infinite-dimensional
state space can be difficult. The consistency and robustness of Bayesian inference procedures
— particularly with respect to perturbations of the prior such as those arising from numerical
approximations — in such settings is a subtle topic, with both positive (Castillo and Nickl,
2014; Doob, 1949; Kleijn and van der Vaart, 2012; Le Cam, 1953) and negative (Diaconis and
Freedman, 1986; Freedman, 1963; Owhadi et al., 2015) results depending upon fine topological
and geometric details.
In the context of computational pipelines, the challenge of eliciting a coherent prior is closely
connected to the challenge of eliciting a single unified prior based on the conflicting input of
multiple experts (French, 2011; Albert et al., 2012).
Consistent Estimation The present paper focused on foundations. Further methodological
work will be required to establish sufficient conditions for when B(µ, An (x† )) collapses to an
atom on a single element q † = Q(x† ) representing the data-generating QoI in the limit as the
amount of information, n, is increased. There are two questions here; (i) when is q † identifiable
from the given information, and (ii) at what rate does B(µ, An (x† )) concentrate on q † .

37

Generalisation and Extensions Two more directions are highlighted for extension of this work.
First, note that in this paper the information operator A : X → A was treated as a deterministic
object. However, in some applications there is auxiliary randomness in the acquisition of information. For our integration example, nodes ti might arise as random samples from a reference
distribution on [0, 1]. Or, observations x(ti ) themselves might occur with measurement error,
for example due to finite precision arithmetic. Then a more elaborate model A : X × Ω → A
would be required, where Ω is a probability space that injects randomness into the information
operator. This is the setting of, for instance, randomised quasi-Monte Carlo methods. Future
work will extend the framework of PNMs to include randomised information operators of this
kind.
As a second direction, recall that in an adaptive algorithm the choice of the information is
made in an iterative procedure that is informed by the information observed up to that point.
For the canonical illustration in Example 3.4 and its generalisations discussed there, it can be
proven that adaptive algorithms do not out-perform non-adaptive algorithms in average case
error (Lee and Wasilkowski, 1986). However, outside this setting adaptation can be beneficial
and should be investigated in the context of Bayesian PNM.
Connection with Probabilistic Programming The central goal of probabilistic programming
(PP) is to automate statistical computation, through symbolic representation of statistical objects and operations on those objects. The formalism of pipelines as graphical models presented
in this work can be compared to similar efforts to establish PP languages (Goodman et al.,
2012). For instance, a method node in a pipeline can be related to a monad aggregating several
distributions into a single output distribution (Ścibior et al., 2015). An important challenge in
PP is the automation of computing conditional distributions (Shan and Ramsey, 2017). Numerical disintegration and extensions thereof might be of independent interest to this field (e.g.
extending Wood et al., 2014).
Acknowledgements CJO was supported by the Australian Research Council (ARC) Centre
of Excellence for Mathematical and Statistical Frontiers. TJS was supported by the Excellence
Initiative of the German Research Foundation (DFG) through the Free University of Berlin.
MG was supported by the Engineering and Physical Sciences (EPSRC) grants EP/J016934/1,
EP/K034154/1, an EPSRC Mathematical Sciences Established Career Research Fellowship and
a Lloyds Register Foundation grant for Programme on Data-Centric Engineering. This material
was based upon work partially supported by the National Science Foundation under Grant DMS1127914 to the Statistical and Applied Mathematical Sciences Institute. Any opinions, findings,
and conclusions or recommendations expressed in this material are those of the author(s) and
do not necessarily reflect the views of the National Science Foundation.
The authors are grateful to Amazon for the provision of AWS credits and to the authors of
the Eigen and Eigency libraries in Python.

Appendices
A. Proofs
Proof of Theorem 3.3. The following observation will be required; the joint density of X and
A = A(X) can be expressed in two ways:
δ(A(x))(da)µ(dx) = µa (dx)A# µ(da)

38

(A.1)

which holds almost everywhere from the definition of a disintegration {µa }a∈A . Note that our
integrability assumption justifies the interchange of integrals from Fubini’s theorem.
The Bayes risk for a Bayesian PNM MBPNM = (A, BBPNM ), BBPNM (µ, a) = Q# µa , can be
expressed as:
Z
R(µ, MBPNM ) = r(x, B(µ, A(x)))µ(dx)
ZZ
=
L(Q(x), q)Q# µA(x) (dq)µ(dx) (since M Bayesian)
ZZZ
=
L(Q(x), q)Q# µa (dq)δ(A(x))(da)µ(dx)
ZZZ
=
L(Q(x), Q(x0 ))µa (dx0 )µa (dx)A# µ(da) (from Eq. (A.1))
On the other hand, let
b(a) ∈ arg min
q∈Q

Z

L(Q(x), q)µa (dx)

be a Bayes act. Then the Bayes risk associated with such a method MBR = (A, BBR ),
BBR (µ, a) = δ(b(a)), can be expressed as:
Z
R(µ, MBR ) = L(Q(x), b(A(x)))µ(dx)
ZZ
=
L(Q(x), b(a))δ(A(x))(da)µ(dx)
ZZ
=
L(Q(x), b(a))µa (dx)A# µ(da) (from Eq. (A.1))
Next we use the inner product structure on Q and the form of the loss function as L(q, q 0 ) =
kq − q 0 k2Q to argue that R(µ, MBPNM ) = 2R(µ, MBR ), which in turn implies that the optimal
information Aµ for Bayesian PNM and A∗µ for ACA are identical.
For this final step, fix a ∈ A and denote the random variables Qa (X) = Q(X) − b(a) that
are induced according to X ∼ µa . Denote by Q̃a an independent copy of Qa generated from
X̃ ∼ µa . The notation E will be used to refer to the expectation taken over X, X̃. Then we
have
Q(X) − Q(X̃) = (Q(X) − b(a)) − (Q(X̃) − b(a))
= Qa (X) − Q̃a (X̃)

and moreover, from Theorem 3.2 the posterior mean of Q(X) is b(a) and thus E[Qa ] = E[Q̃a ] = 0.
Then
Z
R(µ, MBPNM ) = E[kQa − Q̃a k2Q ]A# µ(da)
Z
= E[kQa k2Q − 2hQa , Q̃a iQ + kQ̃aA k2Q ]A# µ(da)
Z
= 2 E[kQa k2Q ]A# µ(da) (since E[Qa ] = 0 and Qa ⊥⊥ Q̃a )
= 2R(µ, MBR )

as required.

39

Proof of Theorem 4.3. Fix f ∈ F and a ∈ A. Then:


Z
kA(x) − akA
1
a
µ(dx)
µδ (f ) = a f (x)φ
Zδ
δ


ZZ
1
kã − akA
= a
f (x)φ
µã (dx)A# µ(dã)
Zδ
δ

Z 
1
kã − akA
= a φ
µã (f )A# µ(dã)
Zδ
δ
Z
= µã (f )A# µaδ (dã).

(from Eq. (A.1))

Thus
|µaδ (f )

a

Z

[µã (f ) − µa (f )]A# µaδ (dã)
Z
α
≤ Cµ kf kF kã − akαA A# µaδ (dã)

− µ (f )| =

(Assumption 4.2).

(A.2)

Now consider the random variable
R :=

kA(X) − akA
δ

(A.3)

induced from X ∼ µ. The existence of a continuous and positive density pA implies that R
also admits a density on [0, ∞), denoted pR,δ . The fact that pA is uniform on an infinitesimal
neighbourhood of a implies that pR,δ (r) is proportional to the surface area of a hypersphere of
radius δr centred on a ∈ A:
pR.δ (r) =

2π n/2
(δr)n−1 (pA (a) + o(1))
Γ( n2 )

(A.4)

This is valid since A is open and the hypersphere will be contained in A for r sufficiently small.
Eq. (A.2) can then be evaluated:


R
A
Z
kã − akαA φ kã−ak
A# µ(dã)
δ
kã − akαA A# µaδ (dã) =
R  kã−akA 
φ
A# µ(dã)
δ
R α
r φ(r)pR,δ (r)dr
= δα R
(change of variables; Eq. A.3).
(A.5)
φ(r)pR,δ (r)dr
R α+n−1
r
φ(r)dr
−−→ R n−1
(from Eq. A.4)
δ↓0
r
φ(r)dr
Cφα
= 0 (< ∞ from Assumption 4.1).
Cφ
Thus, for δ sufficiently small, Eq. (A.5) can be bounded above by δ α (1+ C̄φα ) where C̄φα := Cφα /Cφ0
and “1” is in this case an arbitrary positive constant. This establishes the upper bound
|µaδ (f ) − µa (f )| ≤ Cµα (1 + C̄φα )kf kF δ α
for δ sufficiently small and completes the proof.

40

Proof of Theorem 5.9. To reduce the notation, suppose that the random variables Y1 , . . . , YJ
admit a joint density p(y1 , . . . , yJ ), However, we emphasise that existence of a density is not
required for the proof to hold. To further reduce notation, denote ya:b = (ya , . . . , yb ).
The output of the computation P (M1 , . . . , Mn ) was defined algorithmically in Definition 5.5
and illustrated in Example 5.4. Our aim is to show that this algorithmic output coincides with
the distribution (Qn )# µa on Qn , which is identified in the present notation with p(yJ |y1:I ).
For j ∈ {I + 1, . . . , J}, the coherence condition on Y1 , . . . , YJ translates into the present
notation as p(yj |y1:j−1 ) = p(yj |yπ(j) ). This allows us to deduce that:
Z
Z
p(yJ |y1:I ) = · · · p(yI+1:J |y1:I )dyI+1:J−1
=

=

Z

Z

···
···

Z

Z

J
Y

j=I+1
J
Y

j=I+1

p(yj |y1:j−1 )dyI+1:J−1
p(yj |yπ(j) )dyI+1:J−1 .

The right hand side is recognised as the output of the computation P (M1 , . . . , Mn ), as defined
in Definition 5.5. This completes the proof.

References
N. L. Ackerman, C. E. Freer, and D. M. Roy. On computability and disintegration. Mathematical
Structures in Computer Science, 2017. To appear.
I. Albert, S. Donnet, C. Guihenneuc-Jouyaux, S. Low-Choy, K. Mengersen, and J. Rousseau.
Combining expert opinions in prior elicitation. Bayesian Anal., 7(3):503–531, 2012.
doi:10.1214/12-BA717.
T. V. Anderson. Efficient, accurate, and non-gaussian error propagation through nonlinear,
closed-form, analytical system models. Master’s thesis, Department of Mechanical Engineering, Brigham Young University, 2011.
I. Babuška and G. Söderlind. On round-off error growth in elliptic problems, 2016. In preparation.
S. Bartels and P. Hennig. Probabilistic approximate least-squares. In Proceedings of Artificial
Intelligence and Statistics (AISTATS), 2016.
J. O. Berger. Statistical Decision Theory and Bayesian Analysis. Springer Series in Statistics.
Springer-Verlag, New York, second edition, 1985. doi:10.1007/978-1-4757-4286-2.
A. Beskos, D. Crisan, and A. Jasra. On the stability of sequential Monte Carlo methods in high
dimensions. Ann. Appl. Probab., 24(4):1396–1445, 2014. doi:10.1214/13-AAP951.
A. Beskos, A. Jasra, E. A. Muzaffer, and A. M. Stuart. Sequential Monte Carlo methods for
Bayesian elliptic inverse problems. Stat. Comput., 25(4):727–737, 2015. doi:10.1007/s11222015-9556-7.
A. Beskos, M. Girolami, S. Lan, P. E. Farrell, and A. M. Stuart.
Geometric
MCMC for infinite-dimensional inverse problems. J. Comput. Phys., 335:327–351, 2017.
doi:10.1016/j.jcp.2016.12.041.
P. G. Bissiri, C. C. Holmes, and S. G. Walker. A general framework for updating belief distributions. J. R. Stat. Soc. Ser. B. Stat. Methodol., 78(5):1103–1130, 2016. doi:10.1111/rssb.12158.
V. I. Bogachev. Gaussian Measures, volume 62 of Mathematical Surveys and Monographs.
American Mathematical Society, Providence, RI, 1998. doi:10.1090/surv/062.

41

Z. I. Botev and D. P. Kroese. Efficient Monte Carlo simulation via the generalized splitting
method. Stat. Comput., 22(1):1–16, 2012. doi:10.1007/s11222-010-9201-4.
D. Bradley. The Hydrocyclone: International Series of Monographs in Chemical Engineering,
volume 4. Elsevier, 2013.
M. Briers, A. Doucet, and S. S. Singh. Sequential auxiliary particle belief propagation. In
International Conference on Information Fusion, 2005.
F.-X. Briol, C. J. Oates, M. Girolami, M. A. Osborne, and D. Sejdinovic. Probabilistic integration: A role for statisticians in numerical analysis?, 2016. arXiv:1512.00933v4.
A.-P. Calderón. On an inverse boundary value problem. In Seminar on Numerical Analysis
and its Applications to Continuum Physics (Rio de Janeiro, 1980), pages 65–73. Soc. Brasil.
Mat., Rio de Janeiro, 1980.
M. A. Capistrán, J. A. Christen, and S. Donnet. Bayesian analysis of ODEs: solver optimal accuracy and Bayes factors. SIAM/ASA J. Uncertain. Quantif., 4(1):829–849, 2016.
doi:10.1137/140976777.
I. Castillo and R. Nickl. On the Bernstein–von Mises phenomenon for nonparametric Bayes
procedures. Ann. Statist., 42(5):1941–1969, 2014. doi:10.1214/14-AOS1246.
F. Cérou, P. Del Moral, T. Furon, and A. Guyader. Sequential Monte Carlo for rare event
estimation. Stat. Comput., 22(3):795–808, 2012. doi:10.1007/s11222-011-9231-6.
J. T. Chang and D. Pollard. Conditioning as disintegration. Statist. Neerlandica, 51(3):287–317,
1997. doi:10.1111/1467-9574.00056.
O. A. Chkrebtii, D. A. Campbell, B. Calderhead, and M. A. Girolami. Bayesian solution
uncertainty quantification for differential equations. Bayesian Anal., 11(4):1239–1267, 2016.
doi:10.1214/16-BA1017.
J. Cockayne, C. Oates, T. J. Sullivan, and M. Girolami. Probabilistic meshless methods for
partial differential equations and Bayesian inverse problems, 2016. arXiv:1605.07811v1.
P. R. Conrad, M. Girolami, S. Särkkä, A. M. Stuart, and K. C. Zygalakis. Statistical analysis of differential equations: introducing probability measures on numerical solutions. Stat.
Comput., 2016. doi:10.1007/s11222-016-9671-0.
S. L. Cotter, M. Dashti, and A. M. Stuart. Approximation of Bayesian inverse problems for
PDEs. SIAM J. Numer. Anal., 48(1):322–345, 2010. doi:10.1137/090770734.
T. Cui, Y. Marzouk, and K. Willcox. Scalable posterior approximations for large-scale Bayesian
inverse problems via likelihood-informed parameter and state reduction. J. Comput. Phys.,
315:363–387, 2016. doi:10.1016/j.jcp.2016.03.055.
M. Dashti, S. Harris, and A. Stuart. Besov priors for Bayesian inverse problems. Inverse Probl.
Imaging, 6(2):183–200, 2012. doi:10.3934/ipi.2012.6.183.
M. de Carvalho, G. L. Page, and B. J. Barney. On the geometry of Bayesian inference, 2017.
arXiv:1701.08994.
P. Del Moral. Feynman–Kac Formulae: Genealogical and Interacting Particle Systems with
Applications. Probability and its Applications (New York). Springer-Verlag, New York, 2004.
doi:10.1007/978-1-4684-9393-1.
P. Del Moral, A. Doucet, and A. Jasra. Sequential Monte Carlo samplers. J. R. Stat. Soc. Ser.
B Stat. Methodol., 68(3):411–436, 2006. doi:10.1111/j.1467-9868.2006.00553.x.
P. Del Moral, A. Doucet, and A. Jasra. An adaptive sequential Monte Carlo method for approximate Bayesian computation. Stat. Comput., 22(5):1009–1020, 2012. doi:10.1007/s11222011-9271-y.
C. Dellacherie and P.-A. Meyer. Probabilities and Potential. North-Holland Publishing Co.,
Amsterdam-New York, 1978.

42

P. Diaconis. Bayesian numerical analysis. Statistical Decision Theory and Related Topics IV,
1:163–175, 1988.
P. Diaconis and D. Freedman. Frequency properties of Bayes rules. In Scientific inference, data
analysis, and robustness (Madison, Wis., 1981), volume 48 of Publ. Math. Res. Center Univ.
Wisconsin, pages 105–115. Academic Press, Orlando, FL, 1983.
P. Diaconis and D. A. Freedman. On the consistency of Bayes estimates. Ann. Statist., 14(1):
1–67, 1986. doi:10.1214/aos/1176349830. With a discussion and a rejoinder by the authors.
J. Dick and F. Pillichshammer. Digital Nets and Sequences: Discrepancy Theory and
Quasi–Monte Carlo Integration. Cambridge University Press, 2010.
J. L. Doob. Application of the theory of martingales. In Le Calcul des Probabilités et ses
Applications, Colloques Internationaux du Centre National de la Recherche Scientifique, no.
13, pages 23–27. Centre National de la Recherche Scientifique, Paris, 1949.
A. Doucet, N. de Freitas, and N. Gordon, editors. Sequential Monte Carlo Methods in Practice. Statistics for Engineering and Information Science. Springer-Verlag, New York, 2001.
doi:10.1007/978-1-4757-3437-9.
M. M. Dunlop and A. M. Stuart. The Bayesian formulation of EIT: analysis and algorithms.
Inverse Probl. Imaging, 10(4):1007–1036, 2016. doi:10.3934/ipi.2016030.
L. Ellam, N. Zabaras, and M. Girolami. A Bayesian approach to multiscale inverse
problems with on-the-fly scale determination. J. Comput. Phys., 326:115–140, 2016.
doi:10.1016/j.jcp.2016.08.031.
P. E. Farrell, A. Birkisson, and S. W. Funke. Deflation techniques for finding distinct solutions
of nonlinear partial differential equations. SIAM J. Sci. Comput., 37(4):A2026–A2045, 2015.
doi:10.1137/140984798.
G. E. Fasshauer. Solving differential equations with radial basis functions: multilevel methods
and smoothing. Adv. Comput. Math., 11(2-3):139–159, 1999. doi:10.1023/A:1018919824891.
Radial basis functions and their applications.
D. A. Freedman. On the asymptotic behavior of Bayes’ estimates in the discrete case. Ann.
Math. Statist., 34:1386–1403, 1963. doi:10.1214/aoms/1177703871.
S. French. Aggregating expert judgement. Rev. R. Acad. Cienc. Exactas Fı́s. Nat. Ser. A Math.
RACSAM, 105(1):181–206, 2011. doi:10.1007/s13398-011-0018-6.
N. Garcia Trillos and D. Sanz-Alonso. Gradient flows: Applications to classification, image
denoising, and Riemannian MCMC, 2017. arXiv:1705.07382.
A. Gelman and X.-L. Meng.
Simulating normalizing constants: from importance
sampling to bridge sampling to path sampling.
Statist. Sci., 13(2):163–185, 1998.
doi:10.1214/ss/1028905934.
C. J. Geyer. Markov chain Monte Carlo maximum likelihood. Computing Science and Statistics,
Proceedings of the 23rd Symposium on the Interface, 1991.
M. Girolami and B. Calderhead. Riemann manifold Langevin and Hamiltonian Monte Carlo
methods. J. R. Stat. Soc. Ser. B Stat. Methodol., 73(2):123–214, 2011. doi:10.1111/j.14679868.2010.00765.x. With discussion and a reply by the authors.
N. Goodman, V. Mansinghka, D. M. Roy, K. Bonawitz, and J. B. Tenenbaum. Church: a
language for generative models, 2012. arXiv:1206.3255.
T. Gunter, M. A. Osborne, R. Garnett, P. Hennig, and S. J. Roberts. Sampling for inference
in probabilistic models with fast Bayesian quadrature. In Proceedings of Advances in Neural
Information Processing Systems (NIPS), pages 2789–2797, 2014.
J. A. Gutierrez, T. Dyakowski, M. S. Beck, and R. A. Williams. Using electrical impedance
tomography for controlling hydrocyclone underflow discharge. Powder Technology, 108(2):
180–184, 2000.

43

R. Harvey and D. Verseghy. The reliability of single precision computations in the simulation
of deep soil heat diffusion in a land surface model. Clim. Dynam., 46(3865):3865–3882, 2015.
doi:10.1007/s00382-015-2809-5.
P. Hennig. Probabilistic interpretation of linear solvers. SIAM J. Optim., 25(1):234–260, 2015.
doi:10.1137/140955501.
P. Hennig and M. Kiefel. Quasi-Newton methods: a new direction. J. Mach. Learn. Res., 14:
843–865, 2013.
P. Hennig, M. A. Osborne, and M. Girolami. Probabilistic numerics and uncertainty in computations. Proceedings of the Royal Society A, 471(2179):20150142, 2015.
P. Henrici. Error Propagation for Difference Method. John Wiley and Sons, Inc., New YorkLondon, 1963.
N. J. Higham.
Accuracy and Stability of Numerical Algorithms.
Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2002.
doi:10.1137/1.9780898718027.
M. Horstein. Sequential transmission using noiseless feedback. IEEE Transactions on Information Theory, 9(3):136–143, 1963.
T. E. Hull and J. R. Swenson. Tests of probabilistic models for propagation of roundoff errors.
Communications of the ACM, 9(2):108–113, 1966.
A. T. Ihler and D. A. McAllester. Particle belief propagation. In Proceedings of Artificial
Intelligence and Statistics (AISTATS), 2009.
M. John and Y. Wu. Confidence intervals for finite difference solutions, 2017. arXiv:1701.05609.
J. B. Kadane. Principles of Uncertainty. Texts in Statistical Science Series. CRC Press, Boca
Raton, FL, 2011. doi:10.1201/b11322.
J. B. Kadane and G. W. Wasilkowski. Average case -complexity in computer science: A
Bayesian view. Technical report, Columbia University, 1983.
J. B. Kadane and G. W. Wasilkowski. Bayesian Statistics, chapter Average Case -Complexity
in Computer Science: A Bayesian View, pages 361–374. Elsevier, North-Holland, 1985.
W. Kahan. The Improbability of Probabilistic Error Analyses for Numerical Computations. In
UCB Statistics Colloquium, 1996.
M. Kanagawa, B. K. Sriperumbudur, and K. Fukumizu. Convergence guarantees for kernelbased quadrature rules in misspecified settings. In Advances in Neural Information Processing
Systems, pages 3288–3296, 2016.
T. Karvonen and S. Särkkä. Fully symmetric kernel quadrature, 2017. arXiv:1703.06359.
H. Kersting and P. Hennig. Active uncertainty calibration in Bayesian ODE solvers, 2016.
arXiv:1605.03364.
G. S. Kimeldorf and G. Wahba.
A correspondence between Bayesian estimation on
stochastic processes and smoothing by splines. Ann. Math. Statist., 41:495–502, 1970a.
doi:10.1214/aoms/1177697089.
G. S. Kimeldorf and G. Wahba. Spline functions and stochastic processes. Sankhyā Ser. A, 32:
173–180, 1970b.
D. Kinderlehrer and G. Stampacchia. An Introduction to Variational Inequalities and their
Applications, 2000. doi:10.1137/1.9780898719451. Reprint of the 1980 original.
B. J. K. Kleijn and A. W. van der Vaart. The Bernstein–Von-Mises theorem under misspecification. Electron. J. Stat., 6:354–381, 2012. doi:10.1214/12-EJS675.
A. N. Kolmogorov. Foundations of Probability. Ergebnisse Der Mathematik, 1933.
A. Kong, P. McCullagh, X.-L. Meng, D. Nicolae, and Z. Tan. A theory of statistical models
for Monte Carlo integration. J. R. Stat. Soc. Ser. B Stat. Methodol., 65(3):585–618, 2003.
doi:10.1111/1467-9868.00404. With discussion and a reply by the authors.

44

A. Kong, P. McCullagh, X.-L. Meng, and D. L. Nicolae. Further explorations of likelihood theory for Monte Carlo integration. In Advances in statistical modeling and inference, volume 3 of Ser. Biostat., pages 563–592. World Sci. Publ., Hackensack, NJ, 2007.
doi:10.1142/9789812708298 0028.
J. Koskela, D. Spano, and P. A. Jenkins. Inference and rare event simulation for stopped Markov
processes via reverse-time sequential Monte Carlo, 2016. arXiv:1603.02834.
J. T. N. Krebs. Consistency and asymptotic normality of stochastic Euler schemes for ordinary
differential equations, 2016. arXiv:1609.06880.
J. Kuelbs, F. M. Larkin, and J. A. Williamson. Weak probability distributions on reproducing
kernel Hilbert spaces. Rocky Mt. J. Math., 2(3):369–378, 1972. doi:10.1216/RMJ-1972-2-3369.
F. M. Larkin. Estimation of a non-negative function. BIT Numerical Mathematics, 9(1):30–52,
1969.
F. M. Larkin. Optimal approximation in Hilbert spaces with reproducing kernel functions.
Mathematics of Computation, 24(112):911–921, 1970.
F. M. Larkin. Gaussian measure in Hilbert space and applications in numerical analysis. Rocky
Mt. J. Math., 2(3):379–421, 1972. doi:10.1216/RMJ-1972-2-3-379.
F. M. Larkin. Probabilistic error estimates in spline interpolation and quadrature. In Information processing 74 (Proc. IFIP Congress, Stockholm, 1974), pages 605–609. North-Holland,
Amsterdam, 1974.
F. M. Larkin. A modification of the secant rule derived from a maximum likelihood principle.
BIT, 19(2):214–222, 1979a. doi:10.1007/BF01930851.
F. M. Larkin. Bayesian Estimation of Zeros of Analytic Functions. Queen’s University of
Kingston. Department of Computing and Information Science, 1979b.
S. Lauritzen. Graphical Models. Oxford University Press, 1991.
L. Le Cam. On some asymptotic properties of maximum likelihood estimates and related Bayes’
estimates. Univ. California Publ. Statist., 1:277–329, 1953.
D. Lee and G. W. Wasilkowski. Approximation of linear functionals on a Banach space with a
Gaussian measure. J. Complexity, 2(1):12–43, 1986. doi:10.1016/0885-064X(86)90021-X.
T. Lienart, Y. W. Teh, and A. Doucet. Expectation particle belief propagation. In Proceedings
of Advances in Neural Information Processing Systems (NIPS), 2015.
D. V. Lindley. Understanding Uncertainty. Wiley Series in Probability and Statistics. John
Wiley & Sons, Inc., Hoboken, NJ, revised edition, 2014. doi:10.1002/9781118650158.indsp2.
F. Lindsten, A. M. Johansen, C. A. Naesseth, B. Kirkpatrick, T. B. Schön, J. A. D. Aston,
and A. Bouchard-Côté. Divide-and-conquer with sequential Monte Carlo. J. Comput. Graph.
Statist., 26(2):445–458, 2017. doi:10.1080/10618600.2016.1237363.
D. J. C. MacKay. Bayesian interpolation. Neural Computation, 4(3):415–447, 1992.
M. Mahsereci and P. Hennig. Probabilistic line searches for stochastic optimization. In Proceedings of Advances In Neural Information Processing Systems (NIPS), 2015.
J. Mockus. Bayesian Approach to Global Optimization: Theory and Applications. Springer
Science & Business Media, 1989.
S. Mosbach and A. G. Turner. A quantitative probabilistic investigation into the accumulation
of rounding errors in numerical ODE solution. Computers & Mathematics with Applications,
57(7):1157–1167, 2009.
A. Müller. Integral probability metrics and their generating classes of functions. Adv. in Appl.
Probab., 29(2):429–443, 1997. doi:10.2307/1428011.
S. Niederer, L. Mitchell, N. Smith, and G. Plank. Simulating human cardiac electrophysiology
on clinical time-scales. Front. in Physiol., 2:14, 2011. doi:10.3389/fphys.2011.00014.

45

E. Novak and H. Woźniakowski. Tractability of Multivariate Problems: Standard Information
for Functionals. European Mathematical Society, 2010.
C. Oates, F.-X. Briol, and M. Girolami. Probabilistic integration and intractable distributions,
2016a. arXiv:1606.06841.
C. J. Oates, T. Papamarkou, and M. Girolami. The controlled thermodynamic integral for
Bayesian model evidence evaluation. J. Amer. Statist. Assoc., 111(514):634–645, 2016b.
doi:10.1080/01621459.2015.1021006.
C. J. Oates, J. Cockayne, and R. G. Aykroyd. Bayesian probabilistic numerical methods for
industrial process monitoring. In preparation., 2017.
W. L. Oberkampf and C. J. Roy. Verification and Validation in Scientific Computing. Cambridge University Press, Cambridge, 2013.
A. O’Hagan. Bayes–Hermite quadrature. J. Statist. Plann. Inference, 29(3):245–260, 1991.
doi:10.1016/0378-3758(91)90002-V.
M. Osborne, R. Garnett, Z. Ghahramani, D. K. Duvenaud, S. J. Roberts, and C. E. Rasmussen.
Active learning of model evidence using Bayesian quadrature. In Proceedings of Advances in
Neural Information Processing Systems (NIPS), 2012a.
M. A. Osborne, R. Garnett, S. J. Roberts, C. Hart, S. Aigrain, N. Gibson, and S. Aigrain.
Bayesian quadrature for ratios. In Proceedings of Artificial Intelligence and Statistics (AISTATS), 2012b.
H. Owhadi. Bayesian numerical homogenization. Multiscale Model. Simul., 13(3):812–828, 2015.
doi:10.1137/140974596.
H. Owhadi. Multigrid with rough coefficients and multiresolution operator decomposition from
hierarchical information games. SIAM Rev., 59(1):99–149, 2017. doi:10.1137/15M1013894.
H. Owhadi, C. Scovel, and T. J. Sullivan. On the brittleness of Bayesian inference. SIAM Rev.,
57(4):566–582, 2015. doi:10.1137/130938633.
B. Paige and F. Wood. Inference networks for sequential Monte Carlo in graphical models. In
Proceedings of NIPS, 2015. arXiv:1602.06701.
J. Pfanzagl. Conditional distributions as derivatives. Ann. Probab., 7(6):1046–1050, 1979.
H. Poincaré. Calcul des Probabilités. Gauthier-Villars, 1912.
W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes: The
Art of Scientific Computing. Cambridge University Press, Cambridge, third edition, 2007.
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Inferring solutions of differential equations
using noisy multi-fidelity data. arXiv, 2016. arXiv:1607.04805.
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Numerical Gaussian processes for timedependent and non-linear partial differential equations, 2017. arXiv:1703.10230.
K. Ritter. Average-Case Analysis of Numerical Problems, volume 1733 of Lecture Notes in
Mathematics. Springer-Verlag, Berlin, 2000. doi:10.1007/BFb0103934.
C. Robert and G. Casella. Monte Carlo Statistical Methods. Springer Science & Business Media.,
2013.
G. O. Roberts and R. L. Tweedie. Exponential convergence of Langevin distributions and their
discrete approximations. Bernoulli, 2(4):341–363, 1996. doi:10.2307/3318418.
C. Roy. Review of discretization error estimators in scientific computing. In Proceedings of AIAA
Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition,
2010.
J. Sacks and D. Ylvisaker. Statistical designs and integral approximation. In Proc. Twelfth
Biennial Sem. Canad. Math. Congr. on Time Series and Stochastic Processes; Convexity and
Combinatorics (Vancouver, B.C., 1969), pages 115–136. Canad. Math. Congr., Montreal,
Que., 1970.

46

S. Särkkä, J. Hartikainen, L. Svensson, and F. Sandblom. On the relation between Gaussian
process quadratures and sigma-point methods. Journal of Advances in Information Fusion,
11(1):31–46, 2016. arXiv:1504.05994.
M. Schober, D. K. Duvenaud, and P. Hennig. Probabilistic ODE solvers with Runge–Kutta
means. In Proceedings of Advances in Neural Information Processing Systems (NIPS), 2014.
M. Schober, S. Särkkä, and P. Hennig. A probabilistic model for the numerical solution of initial
value problems, 2016. arXiv:1610.05261v1.
A. Ścibior, Z. Ghahramani, and A. D. Gordon. Practical probabilistic programming with monads. SIGPLAN Notices, 50(12):165–176, 2015.
G. Shafer. A Mathematical Theory of Evidence. Princeton University Press, Princeton, N.J.,
1976.
C. Shan and N. Ramsey. Exact Bayesian inference by symbolic disintegration. In Proceedings
of the 44th ACM SIGPLAN Symposium on Principles of Programming Languages, pages
130–144. ACM, 2017.
J. Skilling. Bayesian solution of ordinary differential equations. In C. R. Smith, G. J. Erickson,
and P. O. Neudorfer, editors, Maximum Entropy and Bayesian Methods, volume 50 of Fundamental Theories of Physics, pages 23–37. Springer, 1992. doi:10.1007/978-94-017-2219-3.
R. Sripriya, M. Kaulaskar, S. Chakraborty, and B. Meikap. Studies on the performance of
a hydrocyclone and modeling for flow characterization in presence and absence of air core.
Chemical Engineering Science, 62(22):6391–6402, 2007.
G. Strang and G. Fix. An Analysis of the Finite Element Method. Englewood Cliffs, NJ:
Prentice-Hall., 1973.
S. H. Strogatz. Nonlinear Dynamics and Chaos. Westview Press, 2014.
A. M. Stuart. Inverse problems: a Bayesian perspective. Acta Numer., 19:451–559, 2010.
doi:10.1017/S0962492910000061.
A. V. Sul0 din. Wiener measure and its applications to approximation methods. I. Izv. Vysš.
Učebn. Zaved. Matematika, 6(13):145–158, 1959.
A. V. Sul0 din. Wiener measure and its applications to approximation methods. II. Izv. Vysš.
Učebn. Zaved. Matematika, 5(18):165–179, 1960.
T. J. Sullivan. Well-posed Bayesian inverse problems and heavy-tailed stable quasi-Banach
space priors, 2016. arXiv:1605.05898.
Z. Tan. On a likelihood approach for Monte Carlo integration. J. Amer. Statist. Assoc., 99
(468):1027–1036, 2004. doi:10.1198/016214504000001664.
O. Teymur, K. Zygalakis, and B. Calderhead. Probabilistic linear multistep methods. In
Proceedings of Advances in Neural Information Processing Systems (NIPS), 2016.
A. Törn and A. Žilinskas. Global Optimization, volume 350 of Lecture Notes in Computer
Science. Springer-Verlag, Berlin, 1989. doi:10.1007/3-540-50871-6.
J. F. Traub, G. W. Wasilkowski, and H. Woźniakowski. Information-Based Complexity. Computer Science and Scientific Computing. Academic Press, Inc., Boston, MA, 1988. With
contributions by A. G. Werschulz and T. Boult.
R. Waeber, P. I. Frazier, and S. G. Henderson. Bisection search with noisy responses. SIAM J.
Control Optim., 51(3):2261–2279, 2013. doi:10.1137/120861898.
R. M. West, S. Meng, R. G. Aykroyd, and R. A. Williams. Spatial-temporal modeling for
electrical impedance imaging of a mixing process. Review of Scientific Instruments, 76(7):
073703, 2005.
F. Wood, J.-W. van de Meent, and V. Mansinghka. A new approach to probabilistic programming inference. In Proceedings of Artificial Intelligence and Statistics (AISTATS), 2014.

47

∗

University of Warwick, j.cockayne@warwick.ac.uk
Newcastle University and Alan Turing Institute, chris.oates@ncl.ac.uk
‡
Free University of Berlin and Zuse Institute Berlin, sullivan@zib.de
§
Imperial College London and Alan Turing Institute, m.girolami@imperial.ac.uk
†

48

Electronic Supplement
to the paper Bayesian Probabilistic Numerical Methods

S1. Philosophical Status of the Belief Distribution
The aim of this section is to discuss in detail the semantic status of the belief distribution µ
in a probabilistic numerical method (PNM). In Section S1.1 we survey historical work on this
topic, while in Section S1.2 more recent literature is covered. Then in Section S1.3 we highlight
some philosophical objections and their counter-arguments.

S1.1. Historical Precedent
The use of probabilistic and statistical methods to model a deterministic mathematical object
can be traced back to Poincaré (1912), who used a stochastic model to construct interpolation
formulae. In brief, Poincaré formulated a polynomial
f (x) = a0 + a1 x + · · · + am xm
whose coefficients ai were modelled as independent Gaussian random variables. Thus Poincaré
in effect constructed a Gaussian measure over the Hilbert space with basis {1, x, . . . , xm }. This
pre-empted Kimeldorf and Wahba (1970a,b) and others, which associated spline interpolation
formulae to the means of Gaussian measures over Hilbert spaces.
The first explicit statistical model for numerical error (of which we are aware) was in the
literature on rounding error in the numerical solution of ordinary differential equations (ODE),
as summarised in Hull and Swenson (1966). Therein it was supposed that rounding, by which
we mean representation of a real number
x = 0.a1 a2 a3 a4 . . .

∈ [0, 1]

in a truncated form
x̂ = 0.a1 a2 a3 a4 . . . an ,
is such that the error e = x − x̂ can be reasonably modelled by a uniform random variable on
[−5 × 10−(n+1) , 5 × 10−(n+1) ]. This implies a distribution µ over the unknown value of x given
x̂. The contribution of Hull and Swenson (1966) and others was to replace the last digit an , in
each stored number that arises in the numerical solution of an ODE, with a uniformly chosen
element of {0, . . . , 9}. This performs approximate propagation of the numerical uncertainty due
to rounding error through further computation and, in their case, induces a distribution over
the solution space of the ODE. Note that this work focused on rounding error, rather than
the (time) discretisation error that is intrinsic to numerical ODE solvers; this could reflect the
limited precision arithmetic that was available from the computer hardware of the period.
Larkin (1972) was an important historical paper for PNMs, being the first to set out the
modern statistical agenda for PNMs:

49

In any particular problem situation we are given certain specific properties of the
solution, e.g. a finite number of ordinate or derivative values at fixed abscissae. If
we can assume no more than this basic information we can conclude only that our
required solution is a member of that class of functions which possesses the given
properties - a tautology which is unlikely to appeal to an experimental scientist!
Clearly, we need to be given, or to assume, extra information in order to make more
definite statements about the required function.
Typically, we shall assume general properties, such as continuity or non-negativity of
the solution and/or its derivatives, and use the given specific properties in order to
assist in making a selection from the class K of all functions possessing the assumed
general properties. We shall choose K either to be a Hilbert space or to be simply
related to one.
This description defines a set K of permissible functions, rather than an explicit distribution
over K, but it is clear that Larkin envisaged numerical analysis as an instance of statistical
estimation:
In the present approach, an a priori localisation is achieved effectively by making
an assumption about the relative likelihoods of elements of the Hilbert space of
possible candidates for the solution to the original problem. Among other things,
this permits, at least in principle, the derivation of joint probability density functions
for functionals on the space and also allows us to evaluate confidence limits on the
estimate of a required functional (in terms of given values of other functionals)
without any extra information about the norm of the function in question.
Later, Diaconis (1988) re-iterated this argument for the construction of K more explicitly,
considering numerical integration of the function



x + x2 + cos(x)
f (x) = exp cosh
.
3 + sin(x3 )
over the unit interval. In particular, Diaconis asked:
“What does it mean to ‘know’ a function?” The formula says some things (e.g. f is
smooth, positive and bounded by 20 on [0, 1]) but there are many other facts about
f that we don’t know (e.g. is f monotone, unimodal or convex?)
This argument was provided as justification for belief distributions that encode certain basic
features, such as the smoothness of the integrand. The belief distributions that were then
considered in Diaconis’ paper were Gaussian distributions on K. Diaconis, as well as Larkin
(1972); Kadane and Wasilkowski (1983), observed that some classical numerical methods are
Bayes rules in this context.
The arguments of these papers are intrinsic to modern PNMs. However, the associated
theoretical analysis of computation under finite information has proceeded outside of statistics,
in the applied mathematical literature, where it is usually presented without a statistical context.
That research is reviewed next.

S1.2. Contemporary Outlook
The mathematical foundations of computation based on finite information are established in
the field of information-based complexity (IBC). The monograph of Traub et al. (1988) presents
the foundations of IBC. In brief, the starting point for IBC is the mantra that

50

To compute fast you need to compute with partial information (∼ Houman Owhadi,
SIAM UQ 2016)
This motivates the search for optimal approximations based on finite information, in either
the worst-case or average-case sense of optimal. The particular development of PNMs that we
presented in the main text is somewhat aligned to average-case analysis (ACA) and we focus
on that literature in what follows.
Among the earliest work on ACA, Sul0 din (1959, 1960) studied numerical integration and L2
function approximation in the setting where µ was induced from the Weiner process, with a
focus on optimal linear methods. Later, Sacks and Ylvisaker (1970) moved from analysis with
fixed µ to analysis over a class of µ defined by the smoothness properties of their covariance
kernels. At the same time Kimeldorf and Wahba (1970a,b) established optimality properties
of splines in reproducing kernel Hilbert spaces in the ACA context. Kadane and Wasilkowski
(1985); Diaconis (1988) discussed the connection between ACA and Bayesian statistics. A
general framework for ACA was formalised in the IBC monograph of Traub et al. (1988), while
Ritter (2000) provides a more recent account.
Game theoretic arguments have recently been explored in Owhadi (2015), who argued that the
optimal prior for probabilistic meshless methods (Cockayne et al., 2016) is a particular Gaussian
measure under a game theoretic framework where the energy norm is the loss function. This
provides one route to the specification of default or objective priors for PNMs which deserves
further exploration in general.
The question of “whose” belief is captured in µ was addressed in Hennig et al. (2015), where
it was argued that the prior information in µ represents that of a hypothetical agent (numerical
analyst) which
[. . . ] we are allowed to design (∼ Michael Osborne, personal correspondence, 2016).
This represents a more pragmatic approach to the design of PNM.

S1.3. Paradise Lost?
Typical numerical algorithms contain several different sources of discretisation error. Consider
the solution of the wave equation: A standard finite element method involves both spatial and
temporal discretisations, a series of numerical quadrature problems, as well as the use of finite
precision arithmetic for all numerical calculations. Yet, decades of numerical analysis have led
to highly optimised computer codes such that these methods can be routinely used. To develop
PNM for solution of the wave equation, which accounts for each separate source of discretisation
error, is it required to unpick and reconstruct such established numerical algorithms? This would
be an unattractive prospect that would detract from further research into PNMs.
Our view is that there is a choice for which discretisation errors to model. In practice the
PNMs implemented in this work were run on floating point precision machines, yet we did not
model rounding error in their output. This was because, in our examples, floating point error
is insignificant compared to discretisation error and so we chose not to model it. This is in line
with the view that a model is a useful simplification of the real world.

S2. Existence of Non-Randomised Bayes Rule
In this section we recall an argument for the general existence of non-randomised Bayes rules,
that was stated without proof in the main text. Sufficient conditions for Fubini’s theorem to
hold are assumed.

51

Proposition S2.1. Let B(A) be non-empty. Then B(A) contains a classical numerical method
of the form B(µ, a) = δ ◦ b(a) where b(a) is a Bayes act for each a ∈ A.
Proof. Let C be the set of belief update operators of the classical form B(µ, a) = δ ◦ b(a).
Suppose there exists a belief update operator B ∗ ∈ B(A) \ C. Then B ∗ can be characterised as
a non-atomic distribution π over the elements of C. Its risk can be computed as:
Z
∗
R(µ, (A, B )) = r(Q(x), B ∗ (µ, A(x)))µ(dx)
ZZ
=
L(Q(x), b(A(x)))π(db)µ(dx)
Z
= R(µ, (A, δ ◦ b))π(db).
If we had R(µ, (A, B ∗ )) < R(µ, (A, δ ◦ b)) for all δ ◦ b ∈ C we would have a contradiction, so it
follows that B(A) ∩ C is non-empty. This completes the proof.

S3. Optimal Information: A Counterexample
In this section we demonstrate that the optimal information Aµ for Bayesian PNM and the
optimal information A∗µ from average case analysis are different in general.
Let X = {♠, ♦, ♥, ♣} be a discrete set, with quantity of interest Q(x) = 1[x = ♠] and
information operator A(x) = 1[x ∈ S] so that Q = A = {0, 1}. In particular, Q is not a vector
space and hence not an inner product space as specified in Theorem 3.3.
Consider two possible choices, S = {♠, ♦} and S = {♠, ♦, ♥}. Assume a uniform prior over
X . Consider the 0-1 loss function L(q, q 0 ) = 1[q 6= q 0 ]. It will be shown that ACA optimal
information for this example can be based on either S = {♠, ♦} or S = {♠, ♦, ♥} whereas
PNM optimal information must be based on S = {♠, ♦, ♥}. Thus Bayesian PNM optimal
information Aµ and ACA optimal information A∗µ need not coincide in general.
The classical case considers a method of the form MBR = (A, BBR ), BBR = δ ◦ b, where
b(a) = 1[a = 0]c0 + 1[a = 1]c1
for some c0 , c1 ∈ {0, 1}. The Bayes risk is
R(µ, MBR ) =

Case of S = {♠, ♦}:

1
4

X

x∈{♠,♦,♥,♣}

1[x ∈
/ S] L(c0 , 1[x = ♠]) + 1[x ∈ S] L(c1 , 1[x = ♠]).

We have

4 R(µ, MBR ) = L(c1 , 1) + L(c1 , 0) + L(c0 , 0) + L(c0 , 0)
= 1[c1 = 0] + 1[c1 = 1] + 2 × 1[c0 = 1]
which is minimised by c1 ∈ {0, 1} and c0 = 0 to obtain a minimum Bayes risk of 41 .
Case of S = {♠, ♦, ♥}: We have
4 R(µ, MBR ) = L(c1 , 1) + L(c1 , 0) + L(c1 , 0) + L(c0 , 0)
= 1[c1 = 0] + 2 × 1[c1 = 1] + 1[c0 = 1]

52

which is minimised by c0 = 0 and c1 = 0 to again obtain a minimum Bayes risk of 14 . Thus the
ACA optimal information can be based on either S = {♠, ♦} or S = {♠, ♦, ♥}.
On the other hand, for the Bayesian PNM we have that MBPNM = (A, BBPNM ), BBPNM =
Q# µA and
R(µ, MBPNM ) =

1
4

X

x∈{♠,♦,♥,♣}

1[x ∈
/ S]L(0, 0)



1
1
+ 1[x ∈ S] (1 −
)L(0, 1[x = ♠]) +
L(1, 1[x = ♠]) .
|S|
|S|

Case of S = {♠, ♦}:

We have
4 R(µ, MBPNM ) =

1 1
+ +0+0
2 2

=

1.

2 1 1
+ + +0
3 3 3

=

4
.
3

Case of S = {♠, ♦, ♥}: We have
4 R(µ, MBPNM ) =

Thus the PNM optimal information is S = {♠, ♦} and not S = {♠, ♦, ♥}. Hence, PNM and
ACA optimal information differ in general.

S4. Monte Carlo Methods for Numerical Disintegration
In this section, Monte Carlo methods for sampling from the distribution µaδ (or µaδ,N ; the N
subscript will be suppressed to reduce notation in the sequel) are considered. The Monte Carlo
approximation of µaδ is, in effect, a problem in rare event simulation as most of the mass of
µaδ will be confined to a set S such that µ(S) is small. Rare events pose some difficulties for
classical Monte Carlo, as an enormous number of draws can be required to study the rare event
of interest.
In the literature there are two major solutions proposed. Importance sampling (Robert and
Casella, 2013) samples from a modified process, under which the event of interest is more likely,
then re-weights these samples to compensate for the adjustment. Conversely, in splitting (Botev
and Kroese, 2012) trajectories of the process are constructed in a genetic fashion, by retaining
and duplicating those which approach the events of interest and discarding others. Splitting is
closely related to SMC (Cérou et al., 2012) and Feynman–Kac models (Del Moral, 2004).
The splitting approach is described in the following section, while in Section S4.3 a parallel
tempering (PT) algorithm is described. In spirit these approaches are similar in that they
employ a tempering approach to ease sampling the relaxed posterior distribution for a small
value of δ. The SMC method employs a particle approximation to accomplish this, while the
PT algorithm uses coupled Markov chains.

S4.1. Sequential Monte Carlo Algorithms for Numerical Disintegration
Let {δi }m
i=0 be such that δ0 = ∞, δm = δ and δi > δi+1 > 0 for all i < m − 1. Furthermore let
a
{Ki }m
i=1 be some set of Markov transition kernels that leave µδi invariant, for which Ki (·, S) is
measurable for all S ∈ ΣX and Ki (x, ·) is an element of PX for all x ∈ X . Then our SMC for
numerical disintegration (SMC-ND) algorithm, based on P particles, is given in Algorithm 1.

53

Sample x0j ∼ µ for j = 1, . . . , P [Initialise]
for i = 1, . . . , m do
Sample xi−1
∼ Ki (xi−1
j
j , ·) for j = 1, . . . , P [Move]
−1
i−1
φ
δ
kA(x
)−akA )
(
i
j
for j = 1, . . . , P [Re-weight]
Set wji ←
−1
φ(δi−1
kA(xi−1
)−akA )
j
P
i P
Sample xij ∼ Discrete({xi−1
j }j=1 ; {wj }j=1 ) for j = 1, . . . , P [Re-sample]
end
Algorithm 1: Sequential Monte Carlo for Numerical Disintegration (SMC-ND).
Here we have used Discrete({xj }Pj=1 ; {wj }Pj=1 ) to denote the discrete distribution which puts
mass proportional to wj on the state xj ∈ X .
The output of the SMC-ND algorithm is an empirical approximation1
µaδm ,P

P
1 X
δ(xm
=
j )
P
j=1

P
to µaδm based on a population of P particles {xm
j }j=1 . There is substantial room to extend and
improve the SMC-ND algorithm based on the wide body of literature available on this subject
(e.g. Doucet et al., 2001; Del Moral et al., 2006; Beskos et al., 2017; Ellam et al., 2016), but
we defer all such improvements for future work. Our aim in the remainder is to establish the
approximation properties of the SMC-ND output. This will be based on theoretical results in
Del Moral et al. (2006).

Assumption S4.1. φ > 0 on R+ .
Assumption S4.2. For all i = 0, . . . , m − 1 and all x, y ∈ X , it holds that Ki+1 (x, ·) 
Ki+1 (y, ·). Furthermore there exist constants i > 0 such that the Radon–Nikodým derivative
dKi+1 (x, ·)
≥ i .
dKi+1 (y, ·)
Assumption S4.1 ensures that Algorithm 1 is well-defined, else it can happen that all particles
are assigned zero weight and re-sampling will fail. However, the result that we obtain in Theorem
S4.3 below can also be established in the special case of an indicator function φ(r) = 1[r < 1].
The details for this variation of the results are also included in the sequel.
The interpretation of Assumption S4.2 is that, for fixed i, transition kernels do not allocate
arbitrarily large or small amounts of mass to different areas of the state space, as a function of
their first argument. This poses a constraint on the choice of Markov kernels for the SMC-ND
algorithm.
Theorem S4.3. For all δ ∈ {δi }m
i=0 and fixed p ≥ 1 it holds that
E



µaδ,P (f ) − µaδ (f )

p  1

p

≤

Cp kf kF
√
P

m−1
for some constant Cp independent of P but dependent on {δi }m
i=0 , p and {i }i=0 .

1

The bandwidth parameter δ and the use of δ to denote an atomic distribution should not be confused.

54

The proof of Theorem S4.3 is presented next. Note that the established bound is independent of
δ ∈ {δi }m
i=0 ; this is therefore a uniform convergence result. The assumptions and the conclusion
of Theorem S4.3 can be weakened in several directions, as discussed in detail in (Del Moral
et al., 2006). Development of SMC methods in the context of high-dimensional and infinitedimensional state spaces has also been considered in Beskos et al. (2014, 2015).

S4.2. Proof of Theorem S4.3
In this section we establish the uniform convergence of the SMC-ND algorithm as claimed in
Theorem S4.3. This relies on a powerful technical result from Del Moral (2004), whose context
is now established.
S4.2.1. Feynman–Kac Models
Let (Ei , Ei ) for i = 0, . . . , m be a collection of measurable spaces. Let η0 be a measure on E0 and
let Γi index a collection of Markov transition kernels from Ei−1 to Ei . Let Gi : Ei → (0, 1] be a
collection of functions, which are referred to as potentials. The triplets (η0 , Gi , Γi ) are associated
with Feynman–Kac measures ηi on Ei defined as, for bounded and measurable functions fi on
Ei ;
ηi (fi ) =

γi (fi )
γi (1)


γi (fi ) = Eη0 fi (X i )

i−1
Y

j=0



Gj (X j )

where the expectation is taken with respect to the Markov process X i defined by X 0 ∼ η0 and
X i |X i−1 ∼ Γi (X i−1 , ·).
The Feynman–Kac measures can be associated with a (non-unique) McKean interpretation
of the form ηi+1 = ηi Λi+1,ηi where the Λi+1,η are a collection of Markov transitions for which
the following compatibility condition holds:
ηΛi+1,η =

Gi
ηΓi+1
η(Gi )

Then the ηi can be interpreted as the ith step marginal distribution of the non-homogeneous
Markov chain defined by X 0 ∼ η0 and X i+1 |X i ∼ Λi+1,ηi (X i , ·). The corresponding P -particle
model is defined on EiP = Ei × · · · × Ei and has
X0 ∼ η0P
P(Xi ∈ dxi |Xi ) =

P
Y

Λi,ηP (Xji−1 , dxij )
i−1

j=1

P
where ηiP = P1 Pj=1 δ(Xji ) is an empirical (random) measure on Ei . The SMC-ND algorithm
can be cast as an instance of such a P -particle model, as is made clear later.
The result that we require from Del Moral (2004) is given next. Denote by Osc1 (Ei ) the set
of measurable functions fi on Ei for which sup{|fi (xi ) − fi (y i )| : xi , y i ∈ Ei } ≤ 1.
Theorem (Theorem 7.4.4 in Del Moral (2004)). Suppose that:

55

i
G
i
i i
(G) There exist G
i ∈ (0, 1] such that Gi (x ) ≥ i Gi (y ) > 0 for all x , y ∈ Ei .

(M1 ) There exist Γi ∈ (0, 1) such that Γi+1 (xi , ·) ≥ Γi Γi+1 (y i , ·) for all xi , y i ∈ Ei .
Then for p ≥ 1 and any valid McKean interpretation Λi,η , the associated P -particle model ηiP
satisfies the uniform (in i) bound
√
sup
sup
P E[|ηiP (fi ) − ηi (fi )|p ]1/p ≤ Cp
0≤i≤m fi ∈Osc1 (Ei )

m
Γ m−1
for some constant Cp independent of P but dependent on {G
i }i=0 and {i }i=0 .

The actual statement in Del Moral (2004) contains a more general version of (M1 ) and a
more explicit decomposition of the constant Cp ; however the simpler version presented here is
sufficient for the purposes of the present paper.
S4.2.2. Case A: Positive Function φ(r) > 0
First we prove Theorem S4.3 as it is stated. Later the assumption of φ > 0 will be relaxed.
SMC-ND as a Feynman–Kac Model The aim here is to demonstrate that the SMC-ND
algorithm fits into the framework of Section S4.2.1 for a specific McKean interpretation. This
connection will then be used to establish uniform convergence for the SMC-ND algorithm as a
consequence of Theorem 7.4.4 in Del Moral (2004).
For the state spaces we associate each Ei = X and Ei = ΣX . For the potentials we associate


1
kA(xi ) − akA
φ δi+1


Gi (xi ) =
φ δ1i kA(xi ) − akA

which clearly does not vanish and takes values in (0, 1] since δi > δi+1 and φ is decreasing. For
the Markov transitions we associate Γi+1 with Ki+1 .
The Feynman–Kac measures associated with the SMC-ND algorithm can be cast as a nonhomogeneous Markov chain with transitions Λi+1,η . Here Λi+1,ηi acts on the current measure
ηi on X by first propagating as ηi Ki+1 and then “warping” this measure with the potential Gi ;
i.e.
Gi
ηΛi+1,η =
ηΓi+1 .
η(Gi )
This demonstrates that the SMC-ND algorithm is the P -particle model corresponding to the
McKean interpretation Λi+1,η of the Feynman–Kac triplet (η0 , Gi , Γi ). Thus the SMC-ND
algorithm can be studied in the context of Section S4.2.1, which we report next.
Note that it is common in applications of SMC to perform the “Re-sample” step before the
“Move” step - our choice of order was required for the McKean framework that is the basis of
the theoretical results in Del Moral et al. (2006). It is known in the SMC “folk lore” that the
order of these steps can be interchanged.
Proof of Uniform Convergence Result for SMC-ND It remains to verify the hypotheses of
Theorem 7.4.4 in Del Moral (2004). Condition (G) is satisfied if and only if


1
i
φ
kA(x ) − akA
δi+1
56

is bounded below, since
φ



1
kA(xi ) − akA
δi



is bounded above by 1. Since φ is continuous, decreasing and satisfies φ > 0 (Assumption S4.1),
1
it suffices to show that its argument δi+1
kA(x) − akA is upper-bounded. This is the content of
Assumption 4.6 in the main text, which shows that
1
1
kA(x) − akA ≤ sup kA(x)kA + kakA
δi
δi x∈X
1
=: G < ∞.
i
Condition (M1 ) requires that
Γi+1 (xi , S) ≥ Γi Γi+1 (y i , S)

for all xi , y i ∈ Ei and S ∈ Ei+1 . From construction this is equivalent to
Ki+1 (xi , S) ≥ Γi Ki+1 (y i , S)
for all xi , y i ∈ X and S ∈ ΣX . This is the content of Assumption S4.2.
Thus we have established the hypotheses of Theorem 7.4.4 in Del Moral (2004) for the SMCND algorithm. Theorem S4.3 is a re-statement of this result. For the statement of the result
we used the kf kF norm, based on the fact that (from Assumption 4.7) kfi kOsc(Ei ) ≤ 2kf k∞ ≤
2CF kf kF .
S4.2.3. Case B: Indicator Function φ(r) = 1[r < 1]
The previous analysis required that φ > 0 on R+ . However, the most basic choice for φ is
the indicator function φ(r) = 1[r < 1] which can take the value 0. The case of an indicator
function demands special attention, since Algorithm 1 can fail in this case if all particles are
assigned zero weight. If this occurs, then we just define µαδ,P (f ) = 0. To be specific, the SMCND algorithm associated to the indicator function φ for approximation of the integral µaδ (f ) is
stated as Algorithm 2 next.
Let Xδa = {x ∈ X : kA(x) − akA < δ}. If there is some iteration i at which, after applying
the kernel Ki to each particle, no particle lies within Xδai , the algorithm fails. As a result it
is critical to ensure that the distance between successive δi is small so that the probability of
failure is controlled. This requirement is made formal next. To establish the approximation
properties of the random measure µaδm ,P , two assumptions are required. These are intended to
replace Assumptions S4.1, S4.2 and Assumption 4.6 from the main text:
Assumption S4.4. For all i = 0, . . . , m − 1 and all xi ∈ Xδai , it holds that Ki+1 (xi , Xδai+1 ) > 0.
Assumption S4.5. For all i = 0, . . . , m − 1 and all xi , y i ∈ Xδai , Ki+1 (xi , ·)  Ki+1 (y i , ·).
Furthermore there exist constants i > 0 such that the Radon–Nikodým derivative
dKi+1 (xi , ·)
≥ i .
dKi+1 (y i , ·)
Assumption S4.4 requires that the probability of reaching Xδai+1 when starting in Xδai and applying the transition kernel Ki+1 , is bounded away from zero. Assumption S4.5 ensures that, for
fixed i, transition kernels do not allocate arbitrarily large or small amounts of mass to different
areas of the state space, as a function of their first argument.

57

Sample x0j ∼ µ for j = 1, . . . , P [Initialise]
for i = 1, . . . , n do
Sample xij ∼ Ki (xi−1
j , ·) for j = 1, . . . , P [Sample]
i
i
Ei ← {xj : xj ∈ Xδai }
if Ei = ∅ then
Return µaδ,P (f ) ← 0
end
for j = 1, . . . , P do
if xij ∈
/ Ei then
xij ∼ Uniform(Ei ) [Re-sample]
end
end
end
P
Return µaδ,P (f ) ← P1 Pj=1 f (xnj ).
Algorithm 2: Sequential Monte Carlo for Numerical Disintegration (SMC-ND), for the case
where φ(r) = 1[r < 1].
Theorem S4.6. For the alternative situation of an indicator function, it holds that for all
δ ∈ {δi }m
i=0 and fixed p ≥ 1,
E



µaδ,P (f ) − µaδ (f )

p  1

p

≤

Cp kf kF
√
P

for some constant Cp independent of P but dependent on p and {i }m−1
i=0 .
Cérou et al. (2012) proposed an algorithm similar to the one herein but focussed on approximation of the probability of a rare event rather than sampling from the rare event itself. In
particular the theoretical results provided are in terms of these probabilities rather than how
well the measure restricted to the rare event is approximated. Furthermore, many of the results
therein focused upon an idealised version of the problem, in which it was assumed that the
intermediate restricted measures can be sampled directly; this avoids the issues with vanishing
potentials indicated in Del Moral (2004). A similar algorithm was discussed in Ścibior et al.
(2015) but was not shown to be theoretically sound.
The remainder of this Section establishes Theorem S4.6.
SMC-ND as a Feynman–Kac Model The aim here is to demonstrate that Algorithm 2 fits
into the framework of Section S4.2.1 for a specific McKean interpretation. This is analogous to
the proof of Theorem S4.3.
A technical complication is that the potentials Gi must take values in (0, 1], which precludes
the “obvious” choice of Ei = X and Gi (xi ) as indicator functions for the sets Xδai . Instead, we
associate Ei = Xδai and Ei with the corresponding restriction of ΣX . For the potentials we then
take Gi (xi ) = 1 for all xi ∈ Ei , which clearly does not vanish and takes values in (0, 1]. For the
Markov transitions Γi+1 from Ei to Ei+1 we consider
Γi+1 (xi , dxi+1 ) ∝ Ki+1 (xi , xi+1 )
which is the restriction of Ki+1 to Ei+1 . For the latter to be well-defined it is required that the

58

normalisation constant

Z

Ki+1 (xi , xi+1 )dxi+1 > 0

Ei+1

xi

for all
∈ Ei , so that there is a positive probability of reaching Ei+1 from Ei . This is the
content of Assumption S4.4.
The Feynman–Kac measures associated with Algorithm 2 can be cast as a non-homogeneous
Markov chain with transitions Λi+1,η . Here Λi+1,ηi acts on the current measure ηi on Ei by
first propagating as ηi Ki+1 and then restricting this measure to Ei+1 . This procedure is seen
to be identical to the Markov transition Γi+1 defined above and, since the potentials Gi ≡ 1, it
follows that
ηΛi+1,η = ηΓi+1
Gi
=
ηΓi+1 .
η(Gi )
This demonstrates that Algorithm 2 is the P -particle model corresponding to the McKean
interpretation Λi+1,η of the Feynman–Kac triplet (η0 , Gi , Γi ). Thus the SMC-ND algorithm can
be studied in the context of Section S4.2.1, which we report next.
Proof of Uniform Convergence Result for SMC-ND It remains to verify the hypotheses of
Theorem 7.4.4 in Del Moral (2004). Condition (G) is satisfied with no further assumption, since
Gi ≡ 1 and we can take G
i = 1. Condition (M1 ) requires that
Γi+1 (xi , S) ≥ Γi Γi+1 (y i , S)
for all xi , y i ∈ Ei and S ∈ Ei+1 . From construction this is equivalent to
Ki+1 (xi , S) ≥ Γi Ki+1 (y i , S)
for all xi , y i ∈ Ei and S ∈ Ei+1 . This is the content of Assumption S4.5.
Thus we have established the hypotheses of Theorem 7.4.4 in Del Moral (2004) for Algorithm
2 and in doing so have established Theorem S4.6.

S4.3. Parallel Tempering for Numerical Disintegration
a
Let Ki , {δi }m
i=1 be as in Section S4.1. The PT algorithm (Geyer, 1991) for sampling from µδm
runs m Markov chains in parallel, one for each temperature, by alternately applying Ki , then
randomly proposing to “swap” the current state of two of the chains. Commonly only swaps of
adjacent chains are considered; to this end suppose at iteration j an index q ∈ {0, . . . , m − 1}
has been selected. Denote by xq the state of the chain with µaδq as its invariant measure. Then
to ensure the correct invariant distribution of all chains is maintained, the swap of state xq and
xq+1 is accepted with probability

α(xq , xq+1 ) =

πq (xq+1 )πq+1 (xq )
πq (xq )πq+1 (xq+1 )

(S4.1)

where πq denotes the density of the target distribution µaδq with respect to a suitable reference
measure. The density notation can be justified since in our experiments the sampler was applied
to the finite-dimensional distributions µaδq ,N and so the reference measure can be taken to be
the Lebesgue measure on RN .

59

Given some initial xi0 for i = 1, . . . , m [Initialise]
for j = 1, . . . , P do
Sample x̂ij ∼ Ki (xij−1 , ·) for i = 1, . . . , m [Move]
Sample q ∼ Uniform(0, m − 1)
if U (0, 1) < α(xqj , xq+1
) then
j
q
q+1
Set xj = x̂j and xq+1
= x̂qj [Accept Swap]
j
else
Set xqj = x̂qj and xq+1
= x̂q+1
[Reject Swap]
j
j
i
i
For i 6= q, q + 1, set xj = x̂j [Update]
end
Algorithm 3: Parallel Tempering for Numerical Disintegration
The PT algorithm for numerical disintegration is described in Algorithm 3. The samples
P
a
{xm
j }j=1 are approximate draws from the distribution µδm .
Algorithms 1 and 3 are each valid for sampling from a target measure µaδ . The choice of which
algorithm to use is problem dependent, and each algorithm has been applied in the experiments
in Section 6.

S4.4. Estimation of Model Evidence
The model evidence pA (a) was estimated as a by-product of the numerical disintegration algorithm developed. Attention is restricted to the specific relaxation function φ(r) = exp(−r2 ).
Then the thermodynamic integral identity (Gelman and Meng, 1998) can be exploited to calculate the model evidence:
Z Z
1 1
log pA (a) = − lim 2
kA(x) − ak2A dµaδ/√t dt
δ↓0 δ
0
√
where the parameterisation δ 7→ δ/ t is such that t = 0 corresponds to the prior, while t = 1
corresponds to the distribution µaδ .
To approximate this integral, the outer integral is first discretised. To this end, fix a sequence
∞ = δ0 < δ1 < · · · < δm of relaxation parameters. For convenience this may be the same
√
sequence as used to apply numerical disintegration. Then for δm small, and letting ti = δm /δi :
Z
m
1 X
log pA (a) ≈ − 2
(ti − ti−1 ) kA(x) − ak2A dµaδm /√ti
δm
i=1

Thus we obtain a consistent approximation
Z
m 
X
1
1
log pA (a) ≈ −
kA(x) − ak2A dµaδi
2 − δ2
δ
i
i−1 |
i=1
{z
}
(∗)

The terms (∗) were estimated via Monte Carlo, based on samples from the distributions µaδi obtained through numerical disintegration. Higher-order quadrature rules and variance reduction
techniques can be used, but were not implemented for this work (Oates et al., 2016b).

60

S4.5. Monte Carlo Details for Painlevé Transcendental
Sampling of the posterior was performed for a temperature schedule of m = 1600 steps, equally
spaced on a logarithmic scale from 10 to 10−4 , for an ensemble of P = 200 particles.
Specification of appropriate transition kernels Ki for this problem was challenging due both
to the high dimension and the empirical observation that, for small δ, mixing of the chains
tends to be poor. This is likely due to the nonlinearity of the information operator which leads
to highly a complex posterior structure. For this reason a gradient-based sampler was used to
construct the transition kernel; the Metropolis-adjusted Langevin algorithm (MALA) (Roberts
and Tweedie, 1996).
Denote by uk the coefficients [ukj ]N
j=1 at iteration k of MALA. Then, recall that MALA has
proposals given by
p
uk+1 = uk + τi Γ∇ log πi (uk ) + 2τi ΓW

where W is a standard Gaussian distribution and Γ ∈ RN ×N is a positive definite preconditioning matrix. The τi were taken to be fixed for each kernel Ki to a value found empirically to
provide a reasonable acceptance rate. πi denotes the unnormalised target distribution for Ki ,
here given by
!
N −a
Ax
πi (uk ) = φ
q N (uk )
δi
P
N
N
where xN = N
i=0 ui φi and q (·) denotes the prior density of the coefficients [uj ]j=1 .
To ensure proposals were scaled to match the decay of the prior for the coefficients, we took
Γ = diag(γ), the diagonal matrix which has the coefficients γi on its diagonal. Even with
such a transition kernel, mixing is generally poor. To compensate k was taken to be large; for
n = 12, 17 we took k = 10, 000, while for n = 22 we took k = 40, 000. We note that such a
large number of temperature levels and transitions makes computation expensive, highlighting
the importance of future work toward methods for approximating the Bayesian posterior in a
more computationally efficient manner.

S4.6. Monte Carlo Details for Poisson Equation
The posterior distribution was obtained by use of the PT algorithm, for m = 20 temperatures
equally spaced on a logarithmic scale between 10−2 and 10−4 . The transition kernels Ki were
given by 10 iterations of a MALA sampler, with preconditioner as described earlier and parameter τ chosen to achieve a good acceptance rate. The number of iterations P was taken to be
106 when n = 25 and 107 when n = 25 or n = 36.

S5. Truncation of the Prior Distribution (Proof of Theorem 4.8)
In this section we present the proof of Theorem 4.8 in the main text. We use a general result
on the well-posedness of Bayesian inverse problems:
Theorem S5.1 (Theorem 4.6 in Sullivan (2016)). Let X and A be separable quasi-Banach
spaces over R. Suppose that
dµaδ
exp(−Φδ (x; a))
=
(S5.1)
dµ
Zδa
where the potential function Φδ satisfies:

61

S0 Φδ (x; ·) is continuous for each x ∈ X , Φδ (·; a) is measurable for each a ∈ A, and for
every r > 0, there exists M0,r,δ ∈ R such that, for all (x, a) ∈ X × A with kxkX < r and
kakA < r,
|Φδ (x; a)| ≤ M0,r,δ .
S1 For every r > 0, there exists a measurable M1,r,δ : R+ → R such that, for all (x, a) ∈ X ×A
with kakA < r,

Φδ (x; a) ≥ M1,r,δ kxkX .
S2 For every r > 0, there exists a measurable M2,r,δ : R+ → R+ such that, for all (x, a, ã) ∈
X × A × A with kakA < r, kãkA < r,

|Φδ (x; a) − Φδ (x; ã)| ≤ exp M2,r,δ kxkX ka − ãkA .

Let Φδ,N be an approximation to Φδ that satisfies (S1-S3) with Mi,r,δ independent of N , and
such that
S3 Ψ : N → R+ is such that, for every r > 0, there exists a measurable M3,r,δ : R+ → R+ ,
such that, for all (x, a) ∈ X × A with kakA < r,

|Φδ,N (x; a) − Φδ (x; a)| ≤ exp M3,r,δ kxkX Ψ(N ).

S4 For some r > 0,



EX∼µ exp(2M3,r,δ (kXkX ) − M1,r,δ (kXkX )) < ∞.

(S5.2)

Let dH denote the Hellinger distance on PX . Then there exists a constant Cδ , independent of
N , such that

dH µaδ,N , µaδ ≤ Cδ Ψ(N )

where µaδ,N is the posterior distribution based on the potential function Φδ,N instead of Φδ .

This allows us to establish conditions on A and µ that guarantee stability under truncation
of the prior:
Proof of Theorem 4.8. Let ϕ be as in Section 4.1, and let


kA(x) − akA
Φδ (x; a) = ϕ
δ


kA ◦ PN (x) − akA
Φδ,N (x; a) = ϕ
.
δ
Our task is to check the conditions of Theorem S5.1 hold for Φδ and Φδ,N .
S0 First, note that Φδ (x; ·) is continuous (since ϕ is continuous from Assumption 4.1 and
Φδ (x; ·) is a composition of continuous functions) and that Φδ (·; a) is measurable (since φ
is measurable and Φδ (·; a) is a composition of measurable functions). Second, note that
ϕ is a continuous bijection from (0, ∞) to itself with ϕ(0) = 0. Thus ϕ−1 exists and we
can consider
δϕ−1 sup{|Φδ (x; a)| : kxkX , kakA < r} = sup{kA(x) − akA : kxkX , kakA < r}
≤ sup kA(x)kA + r
x∈X

≤ ∞ (Assumption 4.6).

Thus we can take M0,r,δ = ϕ( 1δ supx∈X kA(x)kA + rδ ).

62

S1 Since Φδ (x; a) ≥ 0 we can take M1,r,δ = 0.
S2 Given r > 0 let R =
bound

1
δ

supx∈X kA(x)kA + rδ , which is finite by Assumption 4.6. The upper




kA(x) − ãkA
kA(x) − akA
−ϕ
|Φδ (x; a) − Φδ (x; ã)| = ϕ
δ
δ
kA(x) − akA kA(x) − ãkA
≤ CR
−
(Assumption 4.4)
δ
δ
CR
≤
ka − ãkA (reverse triangle inequality)
δ


demonstrates that we can take M2,r,δ = max{0, log( CδR )}.
Minor variation on the above arguments show that S1-3 also hold for Φδ,N with the same
constants Mi,r,δ .
S3 Let CR be defined as in S2. The upper bound




kA ◦ PN (x) − akA
kA(x) − akA
|Φδ,N (x; a) − Φδ (x; a)| = ϕ
−ϕ
δ
δ
kA ◦ PN (x) − akA kA(x) − akA
(Assumption 4.4)
≤ CR
−
δ
δ
CR
≤
kA ◦ PN (x) − A(x)kA (reverse triangle inequality)
δ
CR
≤
exp(m(kxkX ))Ψ(N ) (Assumption 4.5)
δ
demonstrates that we can take M3,r,δ (kxkX ) = max{0, log( CδR ) + m(kxkX )}.
S4 Let CR be defined as in S2. The upper bound
EX∼µ [exp(2M3,r,δ (kXkX ) − M1,r,δ (kXkX ))] = EX∼µ [exp(2 max{0, log(CR /δ) + m(kXkX )})]
CR
≤1+
EX∼µ [exp(2m(kXkX ))]
δ
< ∞ (Assumption 4.5)
establishes the last of the conditions for Theorem S5.1 to hold.

Thus from Theorem S5.1, dH µaδ,N , µaδ ≤ Cδ Ψ(N ). The proof is completed since Assumption
4.7 implies that dF ≤ CF−1 dTV where dTV is the total
√ variation distance based on F = {f :
kf k∞ ≤ 1}; in turn it is a standard fact that dTV ≤ 2dH .

63

