Probabilistic Numerical Methods for
PDE-constrained Bayesian Inverse Problems
Jon Cockayne∗

Chris Oates†
Mark Girolami§

Tim Sullivan

‡

arXiv:1701.04006v1 [stat.ME] 15 Jan 2017

January 17, 2017

This paper develops meshless methods for probabilistically describing discretisation error in the numerical solution of partial differential equations.
This construction enables the solution of Bayesian inverse problems while
accounting for the impact of the discretisation of the forward problem. In
particular, this drives statistical inferences to be more conservative in the
presence of significant solver error. Theoretical results are presented describing rates of convergence for the posteriors in both the forward and inverse
problems. This method is tested on a challenging inverse problem with a
nonlinear forward model.

1 Introduction
Partial differential equations (PDEs) are challenging problems which often have no analytical solution and must be solved numerically. In the style of Probabilistic Numerics
(PN) [7], in this work we describe methods for probabilistically modelling the uncertainty in the true solution arising from the numerical approximation. This uncertainty
can be thought of as arising from finite computation, as formalised in the Information
Complexity literature [11]; in solving a problem numerically, we are forced to discretise
some aspect of it. In the present work we model this uncertainty as arising from taking
a finite number of evaluations of the forcing terms of the system of PDEs.
One of the core principles of probabilistic numerics is that, in complex procedures in
which multiple numerical approximations must be composed to produce a final result,
the uncertainty from each procedure can combine in a nontrivial way which can lead
to incorrect inferences. The example we take here is that of PDE constrained Bayesian
∗

University of Warwick, j.cockayne@warwick.ac.uk
University of Technology Sydney, 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

inverse problems, in which we wish to estimate parameters of a PDE model in a Bayesian
framework, based on observations of a system which is believed to be described by the
underlying PDE. In such problems it has been shown that employing an inaccurate PDE
solver in the sampling can lead to incorrect inferences in the inverse problem [3].
There has been recent interest in construction of probabilistic solvers for PDEs. Work
by [3] constructs a nonparametric posterior distribution for ODEs and PDEs by injecting
noise into standard numerical solvers in such a way as to maintain the convergence
properties of these solvers. In [8], the authors discuss a meshless method which is
similar to the method discussed herein by modelling the forcing of the PDE. This is
developed in [9], which discusses a methodology for probabilistic solution of PDEs by
an hierarchical game-theoretic argument. These latter two approaches do not examine
application to inverse problems, however.
Work from [1] discusses the interpretation of symmetric collocation as the mean function of a Gaussian process prior after conditioning on observed values of the forcing, but
applies this methodology predominantly to stochastic differential equations.

1.1 Structure of the Paper
We begin by introducing the concept of a probabilistic meshless method and giving some
theoretical results related to it. We then show how the posterior measure over the forward
solution of the PDE can be propagated to the posterior measure over parameters in a
Bayesian inverse problem. Finally we present some numerical results for a challenging
nonlinear inverse problem given by the steady-state Allen–Cahn equations.
Proofs for the presented theorems are omitted, and can be found in [2].

2 The Probabilistic Meshless Method
We now introduce the concept of a probabilistic meshless method (PMM). Consider
an open, bounded subset D of Rd with Lipschitz boundary ∂D. We seek a solution
u ∈ H(D), some Hilbert space of functions defined over D, of the following system of
operator equations
Au(x) = g(x)

x∈D

Bu(x) = b(x)

x ∈ ∂D.

(1)

Here A : H(D) → HA (D) and B : H(D) → HB (D) with g ∈ HA (D) and b ∈ HB (D). A
is associated with a partial differential operator and B is associated with the boundary
conditions of the system. For notational simplicity we restrict attention to systems of
two operators, however the methods discussed can be generalised to an arbitrary number
of operator equations.
We proceed in a Bayesian setting by placing a prior measure Πu on u, and determining
its posterior distribution based on a finite number of observations of the system given in
Eq. 1. In this work we focus on the most direct observations of said system; namely, we
choose sets of design points {xi,A } = X0A ⊂ D, {xj,B } = X0B ⊂ ∂D for i = 1, . . . , mA ,

2

j = 1, . . . , mB . We then evaluate the right-hand-side corresponding to each of the
operators in the system at these points; g = [g(xi,A )], b = [b(xi,B )].
It remains to specify our prior distribution. Here we choose a Gaussian process prior
Πu = GP(m, k). Recall that a Gaussian Process is characterised by its mean function m
and its covariance function k, and the property that, if u ∼ GP(m, k) then for any set
of points {xi } ⊂ Rd , i = 1, . . . , n
u(X) ∼ N (µ, Σ)
[µ]i = m(xi )
[Σ]ij = k(xi , xj )
As is common in the literature we will use a centred Gaussian process prior; Πu =
GP(0, k). Define
 


A
L=
L̄ = Ā B̄
B
and furthermore for sets X = {xi }, i = 1, . . . , N , Y = {yj }, j = 1, . . . , M let K(X, Y )
denote the Gram matrix of K applied to X and Y ; [K(X, Y )]ij = k(xi , yj ). Similarly
[AK(X, Y )]ij = Ak(xi , yj ), etc. Then


AĀK(X0A , X0A ) AB̄K(X0A , X0B )
LL̄K(X0 , X0 ) =
B ĀK(X0B , X0A ) B B̄K(X0B , X0B )


AK(X0A , X)
LK(X0 , X) =
BK(X0B , X)


L̄K(X, X0 ) = ĀK(X, X0A ) B̄K(X, X0B )
Here X is to be interpreted as a set of points at which we evaluate those functions drawn
from the posterior distribution, in contrast with X0 = X0A ∪X0B which is the set of points
at which evaluations of the forcing terms are taken.
Proposition 1 (Probabilistic Meshless Method). Assume A and B are linear operators.
Then the posterior distribution Πg,b
u over the solution of the PDE, conditional on the data
g,b
g, b is such that, for u ∼ Πu we have
u(X) ∼ N (µ, Σ)

−1  > > >
g
b
µ = L̄K(X, X0 ) LL̄K(X0 , X0 )

−1
Σ = K(X, X) − L̄K(X, X0 ) LL̄K(X0 , X0 )
LK(X0 , X)

(2)

Note that the mean function in Eq. 2 is the same as the numerical solution to the PDE
that would be obtained using the method of symmetric collocation [5].
Thus far we not discussed the choice of prior covariance k. There are several interesting choices in the literature. Work in [8] proposes use of a covariance which encodes
information about the system through its Green’s function; [2] examined the properties
of this choice in more detail. However, reliance on the Green’s function, which is not

3

in general available in closed-form for complex systems, is a significant drawback. In
practice we will generally posit a prior covariance directly by examining the system in
question and selecting a prior which encodes a suitable level of differentiability.
We now present a theoretical result describing the rate of convergence of the posterior
measure Πg,b
u . Denote by ρ the differential order of the PDE; that is, the maximum
number of derivatives of u required. Furthermore denote by β the smoothness of the
prior; the number of weak derivatives that almost surely exist under the prior measure.
Lastly, define h to be the “fill distance” of the design points X0 :
h = sup min
0

x∈D x ∈X0

x − x0

2

Theorem 2 (Rate of Convergence). For a ball B (u0 ) of radius  centred on the true
solution u0 of the system (1):
 2β−2ρ−d 
h
g,b
c
Πu (B (u0 ) ) = O

where c denotes the set complement.

2.1 Illustrative Example: The Forward Problem
We conclude this section by examining the performance of the probabilistic meshless
method for a simple 1-dimensional PDE. Consider the system
−∇2 u(x) = sin(2πx) x ∈ (0, 1)
u(x) = 0

x = 0, 1

the solution to which can be computed by direct integration to be u(x) = −(2π 2 )−2 sin(2πx).
We compute the PMM solution to this PDE with varying number of design points. In
this setting the Green’s function for the system is available explicitly, and so we used its
associated prior covariance as suggested in [8]; full details are available in [2].
Samples from the posterior distribution can be seen in Fig. 1; note how, even with 20
design points, there is still significant posterior uncertainty. Convergence plots as the
number of design points is increased are shown in Fig. 2

3 Application to Bayesian Inverse Problems
We now turn to an examination of how the PMM, constructed in the previous section,
can be applied in Bayesian inverse problems. We now have a system in which we assume
the operator A depends upon some parameter θ, which we emphasise in the below
system:
Aθ u(x) = g(x)

x∈D

Bu(x) = b(x)

x ∈ ∂D.

4

0.06

0.04

0.02

0.02

0.00

0.00

u

u

0.04

0.06

Mean
Truth
Samples

0.02

0.02

0.04

0.04

0.06
0.0

0.2

0.4

x

0.6

0.8

0.06
0.0

1.0

Mean
Truth
Samples

0.2

0.4

x

0.6

0.8

1.0

Figure 1: Samples from the posterior distribution over the unkown solution to a onedimensional PDE, with mA = 10 (left) and mA = 40 (right).

10 -1

10 -1

10 -2

Tr(Σ)

kµ − uk

10 -2

10 -3

10 -3

10 -4
10 -5

0

20

40

mA

60

80

100

10 -4

0

20

40

mA

60

80

100

Figure 2: Convergence of mean function (left) and posterior covariance trace (right) as
the number of design points mA is increased.

5

In a Bayesian inverse problem we place a prior distribution over θ, θ ∼ Πθ , and seek to
determine its posterior distribution Πyθ based on data y collected at locations {xi } ⊂ D,
i = 1, . . . , n. Further details on Bayesian inverse problems can be found in [10].
Such a posterior distribution is usually intractable and must be investigated by sampling, which involves solution of the underlying system of PDEs as the sampler visits
different values of θ. We assume that the data is obtained by direct observation of the
solution u at these locations, corrupted with Gaussian noise
yi = u(xi ) + ξi
where ξ ∼ N (0, Γ). Our likelihood is thus given by
p(y|θ, u) = N (y; u, Γ)

(3)

where u, y are each vectors in Rn , with [u]i = u(xi ; θ) and [y]i = yi .
Since the solution u to the PDE system is inaccessible it is common to replace u with
an approximation û obtained by some numerical scheme. We instead use the PMM as
the forward solver, obtaining a measure Πg,b
u describing our uncertainty. We may them
marginalise u in Eq. 3 over this measure to obtain
Z
pPN (y|θ) = p(y|θ, u) Πg,b
u (du)
= N (y; µ(θ), Γ + Σ(θ))

(4)

where µ(θ), Σ(θ) are as in Prop. 1, and we have emphasised the dependence on θ. This is
thus similar to the standard approach of replacing u with û in Eq. 3, but we compensate
for the inaccuracy of the forward solver with an additive covariance term Σ incorporating
the uncertainty in the posterior distribution for the forward problem.
We now present a result which guarantees consistency in the inverse problem when
we replace the likelihood in Eq. 3 with that in Eq. 4.
Proposition 3. (Inverse Problem Consistency) Let Πyθ,PN be the posterior distribution
which uses the PN likelihood given in Eq. 4. Assume that the posterior distribution Πyθ
contracts such that Πyθ → δ(θ0 ) as n → ∞, a Dirac measure centred on the true value of
θ, θ0 . Then Πyθ,PN contracts such that Πyθ,PN → δ(θ0 ) provided
h = o(n−1/(β−ρ−d/2) )

3.1 Illustrative Example: The Inverse Problem
We now return to the previous illustrative example to demonstrate the use of a probabilistic solver in the inverse problem. Consideronsider the system
−∇ · θ∇u(x) = sin(2πx) x ∈ (0, 1)
u(x) = 0

x = 0, 1

6

Probabilistic

16

Standard

mA = 4
mA = 8
mA = 16
mA = 32
mA = 64

14
12
10
8
6
4
2
0
0.0

0.5

1.0
θ

1.5

2.0 0.0

0.5

1.0
θ

1.5

2.0

Figure 3: Posterior distributions over θ with varying numbers of design points, on the
left using the PMM, and on the right the standard approach of using a plug-in
estimate for the PDE solution, here given by symmetric collocation.
with the goal of inferring the parameter θ. Data yi was generated from the explicit
solution to this problem with θ = 1 at locations x = 0.25, 0.75, and corrupted with
Gaussian noise with distribution N (0, 0.012 ).
In Fig. 3 we compare posterior distributions for θ generated with the PMM versus the
standard approach of plugging a numerical solution of the PDE into the likelihood and
ignoring discretisation error. The numerical method used in the standard approach was
symmetric collocation, the most natural comparison. Note that when using collocation
the posterior distributions are peaked and biased for small mA , and that the posterior
uncertainty does not appear to depend on the number of design points. Conversely when
using the probabilistic method we see that for small mA the posterior distributions are
wide and flat, while as mA increases the distributions peak and centre on the true value
of θ. Thus, with a standard numerical method the posteriors over θ do not take into
account the quality of the numerical solver used; for poor forward solvers based on
coarse discretisations, the posteriors produced are as confident as those produced with
a fine, accurate numerical solver. With a probabilistic forward solver the variance in
the forward solver is propagated into the inverse problem, resulting in robust inferences
even when the discretisation is coarse.

4 A Nonlinear Example
We now present an application of the methods discussed herein to a nonlinear partial
differential equation known as the steady-state Allen–Cahn system, a model from mathematical physics describing the motion of boundaries between phases in iron alloys. This

7

Negative Stable

1.0

Unstable

1.0

Positive Stable

1.0

1.00
0.8

0.8

0.8

0.6

0.6

0.6

0.75

0.25
0.00

x2

x2

x2

0.50

0.4

0.4

0.4

0.2

0.2

0.2

0.25
0.50

0.0
0.0

0.2

0.4

x1

0.6

0.8

1.0

0.0
0.0

0.2

0.4

x1

0.6

0.8

1.0

0.0
0.0

0.75
1.00
0.2

0.4

x1

0.6

0.8

1.0

Figure 4: Solutions to the Allen–Cahn system for δ = 0.04
is given by
−δ∇2 u + δ −1 (u3 − u) = 0

x ∈ (0, 1)2

u = +1 x1 ∈ {0, 1} , x2 ∈ (0, 1)
u = −1

x2 ∈ {0, 1} , x1 ∈ (0, 1)

(5)

We phrase this as an inverse problem for determining δ. This system is noteworthy
for the fact that it does not admit a unique solution; the three solutions to this system
for δ = 0.04 are shown in Fig. 4. These were generated using the deflation technique
described in [4].
Since this is a nonlinear system the posterior distribution will not be Gaussian, and
we must resort to sampling techniques to explore the posterior distribution. In brief, we
introduce a latent function z and rearrange the system as follows:
−δ∇2 u − δ −1 u = z
δ

−1 3

u = −z

(6)
(7)

Note that by adding Eq. 6 and Eq. 7 we return to the original equation describing the
interior dynamics given in Eq. 5. However Eq. 7 is monotonic and thus invertible; by
inverting this we arrive at a new system:
−δ∇2 u − δ −1 u = z
u = (−δz)1/3
This system is equivalent to the original system but, importantly, is linear. Thus by the
introduction of z we are able to arrive at a new system which can be solved using the
PMM.
It remains to describe z, a latent function whose value is unknown. We seek to
marginalise z in the likelihood
Z
Z
p(y|δ) = p(z|δ) p(y|u) Πg,b,z
(du) dz
(8)
u

8

where Πg,b,z
is now additionally conditioned on a known value for z. This integral is
u
intractable. However by sampling from the posterior distribution over δ by pseudomarginal MCMC it is sufficient to produce an unbiased estimate of this quantity. This
is accomplished by importance sampling; we assume an improper prior p(z|δ) ∝ 1 and
approximate Eq. 8 by the Monte-Carlo estimate
M
1 X
p(y|δ) ≈
M
i=1

R

i
p(y|u) Πg,b,z
(du)
u
r(zi |y, δ)

for zi ∼ r(z|y, δ).
The importance distribution r(z|y, δ) is chosen by solving the original system in Eq. 5
using the techniques described in [4], with a coarse finite-element solver. This gives
estimates {û1 , û2 , û3 } for the solution given a value of δ. By applying Eq. 7 to these
estimates we obtain estimates of three values of z; {ẑ1 , ẑ2 , ẑ3 }.
To handle the multimodality in the solutions we extend the state-space of the inverse
problem to include the solution index j. The importance distribution is constructed as
a Gaussian distribution
z ∼ GP(ẑj , k)
with r(z|y, δ, j) thus the appropriate multivariate Gaussian density after the field for
z has been discretised. Discretisation points are necessarily chosen to match X0A , the
design points for u in the interior of the domain.
For application of the PMM we choose a squared-exponential prior covariance
!
kx − x0 k22
0
k(x, x ) = exp −
2`2
which is known to describe infinitely-differentiable functions. This choice is motivated
by the high differential order required by the PDE; since we must be able to apply both
the operator and the adjoint to the kernel, in this case we require that the covariance
be twice differentiable in each argument, which amounts to a four-times differentiable
covariance if the covariance chosen is isotropic.
The length-scale hyper-parameter ` was incorporated into the MCMC procedure, endowed with a half-Cauchy hyper-prior as recommended in [6]. The parameter of interest
δ was endowed with a uniform prior over the interval (0.02, 0.15), in which the PDE was
empirically found to consistently have three solutions.
Posterior distributions for δ generated using this methodology are shown in Fig. 5;
these are compared with posterior distributions generated using a finite-element forward
solver. In the finite-element case we see a more extreme version of the bias shown in
Fig. 3 for coarse grids, whereas when using a probabilistic forward solver the posteriors
are once again wider to account for an inaccurate forward solver.
We should also comment on the comparison to the finite-element method here; in
the previous example comparison was to the symmetric collocation method for solving
PDEs; in this case the comparison is more direct as the solution for the PDE in symmetric collocation is simply the posterior mean from the PMM. In this case we use a

9

1000
800
600

Probabilistic

1000

PMM ` = 5
PMM ` = 10
PMM ` = 20
PMM ` = 40
PMM ` = 80

800
600

400

400

200

200

0
0.02

0.03

0.04

0.05

δ

0.06

0.07

0.08

0.09

0
0.02

Standard
FEM 5x5
FEM 10x10
FEM 25x25
FEM 50x50

0.03

0.04

0.05

δ

0.06

0.07

0.08

0.09

Figure 5: Posterior distributions for δ obtained by use of the technique described herein
(left) versus a standard Finite Element solver that does not model discretisation error (right).
finite-element solver both to highlight the fact that the behaviour witnessed when using
symmetric collocation is not unique to that solver, and because in existing methods for
finding the multiple solutions to the Allen-Cahn equation the base numerical method
applied is the finite-element method. Furthermore we note that as the underlying numerical method becomes arbitrarily accurate, the posterior inferences made in the inverse
problem should be invariant to the forward solver used.

5 Discussion
We have shown how to construct probabilistic models for the solution of partial differential equations, which quantify the uncertainty arising from numerical discretisation of
the system. We have further shown how the uncertainty in the forward problem can be
propagated into posteriors over parameters in inverse problems. This allows robust inferences to be made in inverse problems, even when the numerical scheme used to solve the
forward problem is inaccurate, which is useful in cases where obtaining highly accurate
solutions is computationally expensive, or where we are willing to tolerate less certain
inferences in exchange for fast computation. In particular we have illustrated how this
might be used to make inferences in nonlinear systems where a variety of phenomena,
such as a non-unique solution could cause a numerical solver to fail.
Immediate extensions to this work lie in examining evolutionary systems in which the
solution is additionally a function of time; the added complexity from the additional
dimension demands more focussed attention. We also seek to examine a more generic
approach for sampling from posterior distributions for nonlinear PDEs. Furthermore we
note that the observations we have chosen for the forward problem are only one possible
choice; another attractive option is given by Galerkin schemes for approximating PDEs,
by choosing our observations to be Galerkin projections.
Lastly we seek to explore other choices of prior. The Gaussian measure is an unrealistic
option in general, as it penalises extreme values and prevents encoding such simple
properties as positivity of solutions.

10

6 Acknowledgements
TJS was supported by the Free University of Berlin within the Excellence Initiative of the
German Research Foundation (DFG). MG was supported by EPSRC [EP/J016934/1,
EP/K034154/1], an EPSRC Established Career Fellowship, the EU grant [EU/259348]
and a Royal Society Wolfson Research Merit Award.
The authors would like to thank John Skilling for useful discussion, Patrick Farrell
for providing code used in generating these results and François-Xavier Briol for helpful
feedback. In addition they express gratitude to the developers of the Python libraries
Autograd and GPyOpt.

References
[1] Igor Cialenco, Gregory E Fasshauer, and Qi Ye. Approximation of stochastic partial
differential equations by a kernel-based collocation method. International Journal
of Computer Mathematics, 89(18):2543–2561, December 2012.
[2] Jon Cockayne, Chris Oates, Tim Sullivan, and Mark Girolami. Probabilistic Meshless Methods for Partial Differential Equations and Bayesian Inverse Problems.
arXiv:1605.07811v1, May 2016.
[3] Patrick R Conrad, Mark Girolami, Simo Särkkä, Andrew Stuart, and Konstantinos Zygalakis. Statistical analysis of differential equations: introducing probability
measures on numerical solutions. Statistics and Computing, 2016.
[4] Patrick E Farrell, Asgeir Birkisson, and Simon W Funke. Deflation techniques for
finding distinct solutions of nonlinear partial differential equations. SIAM Journal
on Scientific Computing, 37(4):A2026–A2045, 2015.
[5] Gregory E Fasshauer. Solving differential equations with radial basis functions:
multilevel methods and smoothing. Advances in Computational Mathematics, 11
(2-3):139–159, 1999.
[6] A Gelman. Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper). Bayesian analysis, 1(3):515–534, 2006.
[7] Philipp Hennig, Michael A Osborne, and Mark Girolami. Probabilistic numerics
and uncertainty in computations. Proc R Soc A, 471(2179):20150142, July 2015.
[8] Houman Owhadi. Bayesian numerical homogenization. Multiscale Modeling & Simulation, 13(3):812–828, 2015.
[9] Houman Owhadi. Multigrid with rough coefficients and multiresolution operator
decomposition from Hierarchical Information Games. arXiv:1503.03467v4, March
2015.

11

[10] Andrew M. Stuart. Inverse problems: a Bayesian perspective. Acta Numer., 19:
451–559, 2010. ISSN 0962-4929.
[11] Henryk Woźniakowski. What is information-based complexity?
complexity of continuous problems, pages 89–95, 2009.

12

Essays on the

