[81]

S.L. Cotter, M. Dashti and A.M. Stuart,

Approximation of Bayesian Inverse Problems.

SIAM Journal of Numerical Analysis 48

(2010) 322­345.

© Society for Industrial and Applied Mathematics

SIAM J. NUMER. ANAL. c 2010 Society for Industrial and Applied Mathematics
!
Vol. 48, No. 1, pp. 322–345

APPROXIMATION OF BAYESIAN
INVERSE PROBLEMS FOR PDES∗
S. L. COTTER† , M. DASHTI† , AND A. M. STUART†

Abstract. Inverse problems are often ill posed, with solutions that depend sensitively on data.
In any numerical approach to the solution of such problems, regularization of some form is needed to
counteract the resulting instability. This paper is based on an approach to regularization, employing
a Bayesian formulation of the problem, which leads to a notion of well posedness for inverse problems,
at the level of probability measures. The stability which results from this well posedness may be
used as the basis for quantifying the approximation, in finite dimensional spaces, of inverse problems
for functions. This paper contains a theory which utilizes this stability property to estimate the
distance between the true and approximate posterior distributions, in the Hellinger metric, in terms
of error estimates for approximation of the underlying forward problem. This is potentially useful as
it allows for the transfer of estimates from the numerical analysis of forward problems into estimates
for the solution of the related inverse problem. It is noteworthy that, when the prior is a Gaussian
random field model, controlling differences in the Hellinger metric leads to control on the differences
between expected values of polynomially bounded functions and operators, including the mean and
covariance operator. The ideas are applied to some non-Gaussian inverse problems where the goal is
determination of the initial condition for the Stokes or Navier–Stokes equation from Lagrangian and
Eulerian observations, respectively.

Key words. inverse problem, Bayesian, Stokes flow, data assimilation, Markov chain–Monte
Carlo

AMS subject classifications. 35K99, 65C50, 65M32

DOI. 10.1137/090770734

1. Introduction. In applications it is frequently of interest to solve inverse prob-

lems [15, 28, 32]: to find u, an input to a mathematical model, given y an observation
of (some components of, or functions of) the solution of the model. We have an
equation of the form

(1.1) y = G(u)

to solve for u ∈ X, given y ∈ Y , where X and Y are Banach spaces. We refer to

evaluating G as solving the forward problem.1 We refer to y as data or observations. It
is typical of inverse problems that they are ill posed: there may be no solution, or the
solution may not be unique and may depend sensitively on y. For this reason some
form of regularization is often employed [8] to stabilize computational approximations.
We adopt a Bayesian approach to regularization [4, 12] which leads to the notion
of finding a probability measure µ on X, containing information about the relative
probability of different states u, given the data y. Adopting this approach is natural
in situations where an analysis of the source of data reveals that the observations y
∗ Received by the editors September 10, 2009; accepted for publication (in revised form) February

19, 2010; published electronically April 2, 2010. This research was supported by the EPSRC, ERC,
and ONR.
http://www.siam.org/journals/sinum/48-1/77073.html
† Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK (s.l.cotter@warwick.ac.

uk, M.Dashti@warwick.ac.uk, A.M.Stuart@warwick.ac.uk).

1 In the applications in this paper G is found from composition of the forward model with some

form of observation operator, such as pointwise evaluation at a finite set of points. The resulting
observation operator is often denoted with the letter H in the atmospheric sciences community [11];
because we need H for Hilbert space later on, we use the symbol G.
322

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 323

are subject to noise. A more appropriate model equation is then often of the form
(1.2) y = G(u) + η,
where η is a mean-zero random variable, whose statistical properties we might know,
or make a reasonable mathematical model for, but whose actual value is unknown to
us; we refer to η as the observational noise. We assume that it is possible to describe
our prior knowledge about u, before acquiring data, in terms of a prior probability
measure µ0 . It is then possible to use Bayes’s formula to calculate the posterior
probability measure µ for u given y.
In the infinite dimensional setting the most natural version of Bayes theorem is a
statement that the posterior measure is absolutely continuous with respect to the prior
[27] and that the Radon–Nikodým derivative (density) between them is determined
by the data likelihood. This gives rise to the formula
dµ 1 ! "
(1.3) (u) = exp −Φ(u; y) ,
dµ0 Z(y)
where the normalization constant Z(y) is chosen so that µ is a probability measure:
#
! "
(1.4) Z(y) = exp −Φ(u; y) dµ0 (u).
X

In the case where y is finite dimensional and η has Lebesgue density ρ this is simply
dµ
(1.5) (u) ∝ ρ(y − G(u)).
dµ0
More generally Φ is determined by the distribution of y given u. We call Φ(u; y) the
potential and sometimes, for brevity, refer to the evaluation of Φ(u; y) for a particular
u ∈ X as solving the forward problem as it is defined through G(·). Note that the
solution to the inverse problem is a probability measure µ which is defined through a
combination of solution of the forward problem Φ, the data y, and a prior probability
measure µ0 . Bayesian and classical regularization are linked via the fact that the
minimizer of the Tikhonov-regularized nonlinear least squares problem coincides with
MAP estimators, which maximize the posterior probability—see section 5.3 of [27].
In general it is hard to obtain information from a formula such as (1.3) for a prob-
ability measure. One useful approach to extracting information is to use sampling:
generate a set of points {u(k) }K k=1 distributed (perhaps only approximately) accord-
ing to µ. In this context it is noteworthy that the integral Z(y) appearing in formula
(1.3) is not needed to enable implementation of MCMC (Markov chain–Monte Carlo)
methods to sample from the desired measure. √ These methods incur an error which
is well understood and which decays as K [17]. However, for inverse problems on
function space there is a second source of error, arising from the need to approximate
the inverse problem in a finite dimensional subspace of dimension N . The purpose of
this paper is to quantify such approximation errors. The key analytical idea is that
we transfer approximation properties of the forward problem Φ into approximation
properties of the inverse problem defined by (1.3). Since the solution to the Bayesian
inverse problem is a probability measure, we will need to use metrics on probability
measures to quantify the effect of approximation. We will employ the Hellinger met-
ric because this leads directly to bounds on the approximation error incurred when
calculating the expectation of functions. The main general results concerning approx-
imation properties are Theorems 2.2 and 2.4, together with Corollary 2.3. The key

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

324 S. L. COTTER, M. DASHTI, AND A. M. STUART

practical implication of the theoretical developments is that it allows the practitioner

to apportion computational resources so as to balance the Monte Carlo error and the
error from finite dimensional approximation.
The main Theorem 2.2 is proved by application of the general “consistency plus
stability implies convergence” approach. By this, we mean that convergence of the
approximation of the forward problem (a form of consistency), along with well posed-
ness of the inverse problem, implies convergence of the posterior measure. Defining
well posedness is, of course, a serious obstacle in most inverse problem formulations.
However, in [5] it was demonstrated that, at the level of probability measures, there
is a useful notion of well posedness, and this was used to prove that the posterior
measure is Lipschitz in the data. Here we develop these ideas to encompass the more
complex perturbations arising from finite dimensional approximation. We also extend
the framework of [5], which allows only for finite data sets, to include potential func-
tions Φ which are not bounded below, hence allowing for data which is a function.
Examples of such inverse problems, and the Lipschitz continuity of the posterior with
respect to the data, are developed in [27], and we build on this framework here.
In section 2 we provide the general approximation theory, for measures µ given by
(1.3), upon which the following two sections build. In section 3 we study the inverse
problem of determining the initial condition for the Stokes equation, given a finite
set of observations of Lagrangian trajectories defined through the time-dependent
velocity field solving the Stokes equation. This section also includes numerical results
showing the convergence of the posterior distribution under refinement of the finite
dimensional approximation, as predicted by the theory. Section 4 is devoted to the
related inverse problem of determining the initial condition for the Navier–Stokes
equation, given direct observation of the time-dependent velocity field at a finite set
of points at positive times.
The underlying philosophy in this paper is to formulate the Bayesian inverse
problem on function space, and then approximate it. This philosophy forms the cor-
nerstone of the review article [27] and is also central to the book [28] where linear
Gaussian problems are studied from this viewpoint. Indeed, for such problems the
idea appeared four decades ago in [9]. An alternative philosophy is to discretize
the forward problem and then apply the Bayesian approach to the resulting inverse
problems in finite dimensions [14], often with some attempt to model the affect of
discretization error statistically. This approach has been very successful in practice
and is readily adapted to a variety of scenarios [2, 13, 21]. There has also been some
work on finite dimensional linear inverse problems, using the Bayesian approach to
regularization, and considering infinite dimensional limits [10, 18]. Another motiva-
tion for discretizing the forward problem before applying the Bayesian formulation is
that it easily allows for white noise priors (identity covariance operators) providing a
direct link to the basic form of Tikhonov regularization used in many inverse problems
[33]. In function space the Bayesian version of this setup leads to technical obstacles
because the identity is not trace class in the natural Hilbert space, and hence the re-
sulting probability measure is supported on (for example) spaces of distributions and
not on the Hilbert space itself [16]. Nonetheless there are significant advantages to
formulating inverse problems on function space, beyond the clarity of the mathemati-
cal formulation, including the possibility of designing sampling methods which do not
degenerate under mesh refinement, as demonstrated in [6], and the ability to transfer
approximation results from classical numerical analysis of forward problems into the
study of inverse problems, as demonstrated in this article. Considerable practical

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 325

experience will need to be gathered in order to fully evaluate the relative merits of
discretizing before or after formulation of the Bayesian inverse problem.
Overviews of inverse problems arising in fluid mechanics, such as those studied
in sections 3 and 4, may be found in [19]. The connection between the Tikhonov-
regularized least squares and Bayesian approaches to inverse problems in fluid me-
chanics is overviewed in [1].
2. General framework. In this section we establish three useful results which
concern the effect of approximation on the posterior probability measure µ given by
(1.3). These three results are Theorem 2.2, Corollary 2.3, and Theorem 2.4. The key
point to notice about these results is that they simply require the proof of various
bounds and approximation properties for the forward problem, and yet they yield
approximation results concerning the Bayesian inverse problem. The connection to
probability comes only through the choice of the space X, in which the bounds and
approximation properties must be proved, which must have full measure under the
prior µ0 .
The probability measure of interest (1.3) is defined through a density with respect
to a prior reference measure µ0 which, by shift of origin, we take to have mean
zero. Furthermore, we assume that this reference measure is Gaussian with covariance
operator C. We write µ0 = N (0, C). In fact, the only consequence of the Gaussian
prior assumption that we use is the Fernique Theorem A.3. Hence the results may be
trivially extended to all measures which satisfy the conclusion of this theorem. The
Fernique theorem holds for all Gaussian measures on a separable Banach space [3]
and also for other measures with tails which decay at least as fast as a Gaussian.
It is demonstrated in [27] that in many applications, including those considered
here, the potential Φ(·; y) satisfies certain natural bounds on a Banach space (X, %·%X ).
It is then natural to choose the prior µ0 so that µ0 (X) = 1. Such bounds on the
forward problem Φ are summarized in the following assumptions. We assume that the
data y lay in a Banach space (Y, % · %Y ). The key point about the form of Assumption
1(i) is that it allows the use of the Fernique theorem to control integrals against µ.
Assumption 1(ii) may be used to obtain lower bounds on the normalization constant
Z(y).
Assumption 1. For some Banach space X with µ0 (X) = 1, the function Φ :
X × Y → R satisfies the following:
(i) for every ε > 0 and r > 0 there is M = M (ε, r) ∈ R such that for all u ∈ X
and y ∈ Y with %y%Y < r
Φ(u; y) ! M − ε%u%2X ;
(ii) for every r > 0 there is a L = L(r) > 0 such that for all u ∈ X and y ∈ Y
with max{%u%X , %y%Y } < r
Φ(u; y) " L(r).
For Bayesian inverse problems in which a finite number of observations are made
and the observation error η is mean zero Gaussian with covariance matrix Γ, the
potential Φ has the form
1
(2.1) Φ(u; y) = |y − G(u)|2Γ ,
2
where y ∈ Rm is the data, G : X → Rm is the forward model, and | · |Γ is a covariance
1
weighted norm on Rm given by | · |Γ = |Γ− 2 · | and | · | denotes the standard Euclidean
norm. In this case it is natural to express conditions on the measure µ in terms of G.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

326 S. L. COTTER, M. DASHTI, AND A. M. STUART

Assumption 2. For some Banach space X with µ0 (X) = 1, the function G : X →

Rm satisfies the following: for every ε > 0 there is M = M (ε) ∈ R such that, for all
u ∈ X,
!
|G(u)|Γ " exp ε%u%2X + M ).

Lemma 2.1. Assume that Φ : X × Rm → R is given by (2.1) and let G satisfy

Assumption 2. Assume also that µ0 is a Gaussian measure satisfying µ0 (X) = 1.
Then Φ satisfies Assumption 1.
Proof. Assumption 1(i) is automatic since Φ is positive; Assumption 1(ii) follows
from the bound

Φ(u; y) " |y|2Γ + |G(u)|2Γ

and the use of the exponential bound on G.

Since the dependence on y is not relevant in this paper, we suppress it notationally
and study measures µ given by
dµ 1 ! "
(2.2) (u) = exp −Φ(u) ,
dµ0 Z
where the normalization constant Z is given by
#
! "
(2.3) Z= exp −Φ(u) dµ0 (u).
X

We approximate µ by approximating Φ over some N -dimensional subspace of X. In

particular, we define µN by

dµN 1 ! "
(2.4) (u) = N exp −ΦN (u) ,
dµ0 Z
where
#
N
! "
(2.5) Z = exp −ΦN (u) dµ0 (u).
X

The potential ΦN should be viewed as resulting from an approximation to the solution

of the forward problem. Our interest is in translating approximation results for Φ into
approximation results for µ. The following theorem proves such a result. Again the
particular exponential dependence of the error constant for the forward approximation
is required so that we may use the Fernique theorem to control certain expectations
arising in the analysis.
Theorem 2.2. Assume that Φ and ΦN satisfy Assumptions 1(i) and 1(ii) with
constants uniform in N . Assume also that, for any ε > 0, there is K = K(ε) > 0
such that
! "
(2.6) |Φ(u) − ΦN (u)| " K exp ε%u%2X ψ(N ),

where ψ(N ) → 0 as N → ∞. Then the measures µ and µN are close with respect to
the Hellinger distance: there is a constant C, independent of N , such that

(2.7) dHell (µ, µN ) " Cψ(N ).

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 327
! "
Consequently all moments of %u%X are O ψ(N ) close. In particular,
! " the mean and,
in the case X is a Hilbert space, the covariance operator are O ψ(N ) close.
Proof. Throughout the proof, all integrals are over X. The constant C may
depend upon r and changes from occurrence to occurrence. Using Assumption 1(ii)
gives
#
! " ! "
|Z| ! exp −L(r) dµ0 (u) ! exp −L(r) µ0 {%u%X " r}.
{$u$X !r}

This lower bound is positive because µ0 has a full measure on X and is Gaussian so
that all balls in X have positive probability. We have an analogous lower bound for
|Z N |.
From Assumptions 1(i) and (2.6), using the fact that µ0 is a Gaussian probability
measure so that the Fernique Theorem A.3 applies, we obtain
#
! " ! "
|Z − Z N | " Kψ(N ) exp ε%u%2X − M exp ε%u%2X dµ0 (u)

" Cψ(N ).

From the definition of Hellinger distance we have

# $ $ % $ %%2
N 2 − 12 1 N − 12 1 N
2dHell (µ, µ ) = Z exp − Φ(u) − (Z ) exp − Φ (u) dµ0 (u)
2 2
" I1 + I2 ,

where
# $ $ % $ %%2
2 1 1
I1 = exp − Φ(u) − exp − ΦN (u) dµ0 (u),
Z 2 2
#
1 1 "
I2 = 2|Z − 2 − (Z N )− 2 |2 exp(−ΦN (u) dµ0 (u).

Now, again using Assumption 1(i) and (2.6), together with the Fernique Theorem
A.3,
#
Z 1 2 ! "
I1 " K ψ(N )2 exp 3ε%u%2X − M dµ0 (u)
2 4
" Cψ(N )2 .

Note that the bounds on Z, Z N from below are independent of N . Furthermore,

# #
! " ! "
exp −Φ (u) dµ0 (u) " exp ε%u%2X − M dµ0 (u)
N

with bound independent of N , by the Fernique Theorem A.3. Thus

! "
I2 " C Z −3 ∨ (Z N )−3 |Z − Z N |2
" Cψ(N )2 .

Combining gives the desired continuity result in the Hellinger metric.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

328 S. L. COTTER, M. DASHTI, AND A. M. STUART

Finally all moments of u in X are finite under the Gaussian measure µ0 by the
Fernique Theorem A.3. It follows that all moments are finite under µ and µN because,
for f : X → Z polynomially bounded,
! "1 ! "1
Eµ %f % " Eµ0 %f %2 2 Eµ0 exp(−2Φ(u; y)) 2 ,

and the first term on the right-hand side is finite since all moments are finite under
µ0 , while the second term may be seen to be finite by the use of Assumption 1(i) and
the Fernique Theorem A.3.
For Bayesian inverse problems with finite data the potential Φ has the form given
in (2.1) where y ∈ Rm is the data, G : X → Rm is the forward model, and | · |Γ is a
covariance weighted norm on Rm . In this context the following corollary is useful.
Corollary 2.3. Assume that Φ is given by (2.1) and that G is approximated by
a function G N with the property that, for any ε > 0, there is K % = K % (ε) > 0 such
that
! "
(2.8) |G(u) − G N (u)| " K % exp ε%u%2X ψ(N ),

where ψ(N ) → 0 as N → ∞. If G and G N satisfy Assumption 2 uniformly in N ,

then Φ and ΦN := 12 |y − G N (u)|2Γ satisfy the conditions necessary for application of
Theorem 2.2 and all the conclusions of that theorem apply.
Proof. That (i) and (ii) of Assumption 1 hold follows as in the proof of Lemma
2.1. Also (2.6) holds since (for some K(·) defined in the course of the following chain
of inequalities)

1
|Φ(u) − ΦN (u)| " |2y − G(u) − G N (u)|Γ |G(u) − G N (u)|Γ
&2 ! "' ! "
" |y| + exp ε%u%2X + M × K % (ε) exp ε%u%2X ψ(N )
" K(2ε) exp(2ε%u%2X )ψ(N )

as required.
A notable fact concerning Theorem 2.2 is that the rate of convergence attained
in the solution of the forward problem, encapsulated in the approximation of the
function Φ by ΦN , is transferred into the rate of convergence of the related inverse
problem for measure µ given by (2.2) and its approximation by µN . Key to achieving
this transfer of rates of convergence is the dependence of the constant in the forward
error bound (2.6) on u. In particular, it is necessary that this constant is integrable
by use of the Fernique Theorem A.3. In some applications it is not possible to obtain
such dependence. Then convergence results can sometimes still be obtained, but at
weaker rates. We now describe a theory for this situation.
Theorem 2.4. Assume that Φ and ΦN satisfy Assumptions 1(i) and 1(ii) with
constants uniform in N . Assume also that for any R > 0 there is K = K(R) > 0
such that for all u with %u%X " R

(2.9) |Φ(u) − ΦN (u)| " Kψ(N ),

where ψ(N ) → 0 as N → ∞. Then the measures µ and µN are close with respect to
the Hellinger distance:

(2.10) dHell (µ, µN ) → 0

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 329

as N → ∞. Consequently all moments of %u%X under µN converge to corresponding

moments under µ as N → ∞. In particular the mean, and, in the case X is a Hilbert
space, the covariance operator converge.
Proof. Throughout the proof, all integrals are over X unless specified otherwise.
The constant C changes from occurrence to occurrence. The normalization constants
Z and Z N satisfy lower bounds which are identical to that proved for Z in the course
of establishing Theorem 2.2.
From Assumption 1(i) and (2.9),
#
! " ! "
|Z − Z N | " | exp −Φ(u) − exp −ΦN (u) |dµ0
#X
! "
" exp ε%u%2X − M |Φ(u) − ΦN (u)|dµ0 (u)
{$u$X !R}
#
! "
+ 2 exp ε%u%2X − M dµ0 (u)
{$u$X >R}
! "
" exp εR2 − M K(R)ψ(N ) + JR
:= K1 (R)ψ(N ) + JR .
Here
#
! "
JR = 2 exp ε%u%2X − M dµ0 (u).
{$u$X >R}

Now, again by the Fernique Theorem A.3, JR → 0 as R → ∞ so, for any δ > 0, we
may choose R > 0 such that JR < δ. Now choose N > 0 so that K1 (R)ψ(N ) < δ to
deduce that |Z − Z N | < 2δ. Since δ > 0 is arbitrary, this proves that Z N → Z as
N → ∞.
From the definition of Hellinger distance we have
# $ $ % $ %%2
N 2 − 12 1 N − 12 1 N
2dHell (µ, µ ) = Z exp − Φ(u) − (Z ) exp − Φ (u) dµ0 (u)
2 2
" I1 + I2 ,

where
# $ $ % $ %%2
2 1 1 N
I1 = exp − Φ(u) − exp − Φ (u) dµ0 (u),
Z 2 2
#
1 1 "
I2 = 2|Z − 2 − (Z N )− 2 |2 exp(−ΦN (u) dµ0 (u).

Now, again using Assumption 1(i) and (2.9),

#
1 ! "
I1 " K(R)2 ψ(N )2 exp ε%u%2X − M dµ0 (u)
2Z {$u$X !R}
#
4 ! "
+ 2 exp ε%u%2X − M dµ0 (u)
Z {$u$X >R}
1 4
" K2 (R)ψ(N )2 + JR
2Z Z
for suitably chosen K2 = K2 (R). An argument similar to the one above for |Z − Z N |
shows that I1 → 0 as N → ∞.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

330 S. L. COTTER, M. DASHTI, AND A. M. STUART

Note that the bounds on Z, Z N from below are independent of N . Furthermore,

# #
! " ! "
exp −ΦN (u) dµ0 (u) " exp ε%u%2X − M dµ0 (u)

with bound independent of N , by the Fernique Theorem A.3. Thus

1 1 ! "
|Z − 2 − (Z N )− 2 |2 " C Z −3 ∨ (Z N )−3 |Z − Z N |2 ,
and so I2 → 0 as N → ∞. Combining gives the desired continuity result in the
Hellinger metric.
The proof may be completed by the same arguments used in Theorem 2.2.
In section 4 of [27], Theorem 2.2 is applied to the linear Gaussian problem of de-
termining the initial condition of the heat equation from observation of the solution
at a positive time. In the next two sections we show how Corollary 2.3 and Theorem
2.4 may be applied to study finite dimensional approximation of non-Gaussian in-
verse problems arising in the determination of the initial condition for two nonlinear,
dissipative models arising in fluid mechanics.
3. Lagrangian data assimilation. In this section we study an idealized model
of Lagrangian data assimilation as practiced in oceanography. The setup is exactly
that in [5] except that we consider the Stokes equations in place of the Navier–Stokes
equations. The velocity field v given by the incompressible Stokes (ι = 0) or Navier–
Stokes (ι = 1) equations is given by
∂v
(3.1a) + ιv · ∇v = ν∆v − ∇p + f, (x, t) ∈ D × [0, ∞),
∂t

(3.1b) ∇ · v = 0, (x, t) ∈ D × [0, ∞),

(3.1c) v = u, (x, t) ∈ D × {0}.

Here D is the unit square. We impose periodic boundary conditions on the velocity
field v and the pressure p. We assume that f has zero average over D, noting that
this implies the same for v(x, t), provided that u(x) = v(x, 0) has zero initial average.
The PDE can be formulated as a linear dynamical system on the Hilbert space
( )# *
)
(3.2) 2 )
H = u ∈ Lper (D) ) udx = 0, ∇ · u = 0 ,
D

with the usual L2 (D) norm and inner product on this subspace of L2per (D). Through-
out this article A denotes the (self-adjoint, positive) Stokes operator on T2 and
P : L2per → H the Leray projector [29, 30]. The operator A is densely defined on
H and is the generator of an analytic semigroup. We denote by {(φk , λk )}k∈K a
complete orthonormal set of eigenfunctions/eigenvalues for A in H. We then define
fractional powers of A by
+
(3.3) Aα u = λα
k +u, φk ,φk .
k∈K

For any s ∈ R we define the Hilbert spaces Hs by

, -
+
s s 2
(3.4) H = u: λk |+u, φk ,| < ∞ .
k∈K

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 331

The norm in Hs is denoted by % · %s and is given by

+
%u%2s = λsk |+u, φk ,|2 .
k∈K

Of course, H0 = H.
If we let ψ = P f , then in the Stokes case ι = 0 we may write (3.1) as an ODE in
Hilbert space H:

dv
(3.5) + νAv = ψ, v(0) = u.
dt
Our aim is to determine the initial velocity field u from Lagrangian data. To be
precise we assume that we are given noisy observations of J tracers with positions zj
solving the integral equations
# t
(3.6) zj (t) = zj,0 + v(zj (s), s)ds.
0

For simplicity assume that we observe all the tracers z at the same set of positive
times {tk }K
k=1 and that the initial particle tracer positions zj,0 are known to us:

(3.7) yj,k = zj (tk ) + ηj,k , j = 1, . . . , J, k = 1, . . . , K,

where the ηj,k ’s are zero mean Gaussian random variables. Concatenating data we
may write

(3.8) y = G(u) + η
∗ ∗
with y = (y1,1 , . . . , yJ,K )∗ and η ∼ N (0, Γ) for some covariance matrix Γ capturing the
correlations present in the noise. Note that G is a complicated function of the initial
condition for the Stokes equations, describing the mapping from this initial condition
into the positions of Lagrangian trajectories at positive times. We will show that the
function G maps H into R2JK and is continuous on a dense subspace of H.
The objective of the inverse problem is thus to find the initial velocity field u,
given y. We adopt a Bayesian approach, place a prior µ0 (du) on u, and identify
the posterior µ(du) = P(u|y)du. We now spend some time developing the Bayesian
framework, culminating in Theorem 3.3 which shows that µ is well defined. The
reader interested purely in the approximation of µ can skip straight to Theorem 3.4.
The following result shows that the tracer equations (3.6) have a solution, under
mild regularity assumptions on the initial data. An analogous result is proved in [7]
for the case where the velocity field is governed by the Navier–Stokes equation, and
the proof may be easily extended to the case of the Stokes equations.
Theorem 3.1. Let ψ ∈ L2 (0, T ; H) and let v ∈ C([0, T ]; H) denote the solution of
(3.5) with initial data u ∈ H. Then the integral equation (3.6) has a unique solution
z ∈ C([0, T ], R2 ).
We assume throughout that ψ is sufficiently regular that this theorem applies. To
determine a formula for the probability of u given y, we apply the Bayesian approach
described in [5] for the Navier–Stokes equations and easily generalized to the Stokes
equations. For the prior measure we take µ0 = N (0, βA−α ) for some β > 0, α > 1,
with the condition on α chosen to ensure that draws from the prior are in H, by

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

332 S. L. COTTER, M. DASHTI, AND A. M. STUART

Lemma A.4. We condition the prior on the observations, to find the posterior measure
on u. The likelihood of y given u is
& 1 '
P (y | u) ∝ exp − |y − G(u)|2Γ .
2
This suggests the formula
dµ
(3.9) (u) ∝ exp(−Φ(u; y)),
dµ0
where
1
(3.10) Φ(u; y) := |y − G(u)|2Γ
2
and µ0 is the prior Gaussian measure. We now make this assertion rigorous. The first
step is to study the properties of the forward model G. Proof of the following lemma
is given after statement and proof of the main approximation result, Theorem 3.4.
Lemma 3.2. Assume that ψ ∈ C([0, T ]; Hγ ) for some γ ! 0. Consider the forward
model G : H → R2JK defined by (3.7) and (3.8).
• If γ ! 0, then for any . ! 0 there is C > 0 such that for all u ∈ H#
! "
|G(u)| " C 1 + %u%# .

• If γ > 0, then for any . > 0 and R > 0 and for all u1 , u2 with %u1 %# ∨%u2 %# <
R, there is L = L(R) > 0 such that

|G(u1 ) − G(u2 )| " L%u1 − u2 %# .

Furthermore, for any ε > 0, there is M > 0 such that L(R) " M exp(εR2 ).
Thus G satisfies Assumption 2 with X = Hs and any s ! 0.
Since G is continuous on H# for . > 0 and since, by Lemma A.4, draws from µ0
are almost surely in Hs for any s < α − 1, use of the techniques in [5], employing the
Stokes equation in place of the Navier–Stokes equation, shows the following.
Theorem 3.3. Assume that ψ ∈ C([0, T ]; Hγ ), for some γ > 0, and that the
prior measure µ0 = N (0, βA−α ) is chosen with β > 0 and α > 1. Then the measure
µ(du) = P(du|y) is absolutely continuous with respect to the prior µ0 (du), with the
Radon–Nikodým derivative given by (3.9).
In fact, the theory in [5] may be used to show that the measure µ is Lipschitz
in the data y, in the Hellinger metric. This well posedness underlies the following
study of the approximation of µ in a finite dimensional space. We define P N to
be the orthogonal projection in H onto the subspace spanned by {φk }|k|!N ; recall
that k ∈ K := Z2 \{0}. Since P N is an orthogonal projection in any Ha , we have
%P N u%X " %u%X . Define

G N (u) := G(P N u).

The approximate posterior measure µN is given by (3.9) with G replaced by G N . This

approximate measure coincides with the prior on the orthogonal complement of P N H.
On P N H itself the measure is finite dimensional and amenable to sampling techniques
as demonstrated in [6]. We now quantify the error arising from the approximation of
G in the finite dimensional subspace P N X.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 333

Theorem 3.4. Let the assumptions of Theorem 3.3 hold. Then, for any q <
α − 1, there is a constant c > 0, independent of N , such that dHell
! (µ, µ"N ) " cN −q .
N −q
Consequently the mean and covariance operator of µ and µ are O N close in the
H and H-operator norms, respectively.
Proof. We set X = Hs for any s ∈ (0, α − 1). We employ Corollary 2.3. Clearly,
since G satisfies Assumption 2 by Lemma 3.2, so too does G N , with constants uniform
in N. It remains to establish (2.8). Write u ∈ Hs as
+
u= uk φk
k∈K

and note that

+
|k|2s |uk |2 < ∞.
k∈K

We have, for any . ∈ (0, s),

+
%u − P N u%2# = |k|2# |uk |2
|k|>N
+
= |k|2(#−s) |k|2s |uk |2
|k|>N
+
" N −2(s−#) |k|2s |uk |2
|k|>N

" C%u%2s N −2(s−#) .

By the Lipschitz properties of G from Lemma 3.2 we deduce that, for any . ∈ (0, s),
! "
|G(u) − G(P N u)| " M exp ε%u%2# %u − P N u%#
1 ! "
" C 2 M exp ε%u%2s %u%s N −(s−#) .
N
This!establishes
" the desired error bound (2.8). It follows from Corollary 2.3 that µ
−(s−#)
is O N close to µ in the Hellinger distance. Choosing s arbitrarily close to its
upper bound, and . arbitrarily close to zero, yields the optimal exponent q as appears
in the theorem statement.
Proof of Lemma 3.2. Throughout the proof, the constant C may change from
instance to instance, but is always independent of the ui . It suffices to consider a
single observation so that J = K = 1. Let z (i) (t) solve
# t
(i)
z (i) (t) = z0 + v (i) (z (i) (τ ), τ )dτ,
0

where v (i) (x, t) solves (3.1) with u = ui .

Let . ∈ [0, 2 + γ). By (A.3),
$ %
1
(3.11) %v (t)%s " C (s−#)/2 %ui %# + %ψ%C([0,T ];Hγ )
(i)
t
for s ∈ [., 2 + γ). Also, by linearity and (A.2),
C
(3.12) %v (1) (t) − v (2) (t)%s " %u1 − u2 %# .
t(s−#)/2

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

334 S. L. COTTER, M. DASHTI, AND A. M. STUART

To prove the first part of the lemma note that, by the Sobolev embedding theorem,
for any s > 1,
# t
(i)
|z (i) (t)| " |z0 | + %v (i) (·, τ )%L∞ dτ
0
& # t '
"C 1+ %v (i) (·, τ )%s dτ
0
& # t '
1
"C 1+ (s−#)/2
%u %
i # dτ .
0 τ
.
For any γ ! 0 and . ∈ [0, 2 + γ) we may choose s such that s ∈ [., 2 + γ) (1, . + 2).
Thus the singularity is integrable and we have, for any t ! 0,
! "
|z (i) (t)| " C 1 + %ui %#

as required.
To prove the second part of the lemma, choose . ∈ (0, 2 + γ) and then choose
s ∈ [. − 1, 1 + γ) ∩ (1, . + 1); this requires γ > 0 to ensure a nonempty intersection.
Then
& 1 '
(3.13) %v (i) (t)%1+s " C (1+s−#)/2 %ui %# + %ψ%C([0,T ];Hγ ) .
t
Now we have
# t
|z (1)
(t) − z (2)
(t)| " |z (0) − z (0)| +
(1) (2)
|v (1) (z (1) (τ ), τ ) − v (2) (z (2) (τ ), τ )|dτ
0
# t
" %Dv (1) (·, τ )%L∞ |z (1) (τ ) − z (2) (τ )|dτ
0
# t
+ %v (1) (·, τ ) − v (2) (·, τ )%L∞ dτ
0
# t
" %v (1) (·, τ )%1+s |z (1) (τ ) − z (2) (τ )|dτ
0
# t
+ %v (1) (·, τ ) − v (2) (·, τ )%s dτ
0
# t & '
1
" C (1+s−#)/2 %u1 %# + %ψ%C([0,T ];Hγ ) |z (1) (τ ) − z (2) (τ )|dτ
0 τ
# t
C
+ (s−#)/2
%u1 − u2 %# dτ.
0 τ
Both time singularities are integrable and application of the Gronwall inequality from
Lemma A.1 gives, for some C depending on %u1 %# and %ψ%C([0,T ];Hγ ) ,

%z (1) − z (2) %L∞ ((0,T );R2 ) " C%u1 − u2 %# .

The desired Lipschitz bound on G follows. In particular, the desired dependence of

the Lipschitz constant follows from the fact that for any ε > 0 there is M > 0 with
the property that, for all θ ! 0,

1 + θ exp(θ) " M exp(εθ2 ).

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 335

We conclude this section with the results of numerical experiments illustrating

the theory. We consider the Stokes equations with viscosity ν = 0.05. We compute
the posterior distribution on the initial condition from observation of J Lagrangian
trajectories at one time t = 0.1. The prior measure is taken to be N (0, 400 × A−2 ).
The initial condition used to generate the data is found by making a single draw from
the prior measure and the observational noise on the Lagrangian data is independently
and identically distributed (i.i.d.) N (0, γ 2 ) with γ = 0.01.
Note that, in the periodic geometry assumed here, the Stokes equations can be
solved exactly by Fourier analysis [30]. Thus there are four sources of approximation
when attempting to sample the posterior measure on u. These are
(i) the effect of generating approximate samples from the posterior measure by
use of MCMC methods;
(ii) the effect of approximating u in a finite subspace found by truncating the
eigenbasis of the Stokes operator;
(iii) the effect of interpolating a velocity field on a grid, found from use of the
FFT, into values at the arbitrary locations of Lagrangian tracers;
(iv) the effect of the time-step in an Euler integration of the Lagrangian trajectory
equations.
The MCMC method that we use is a generalization of the random walk Metropo-
lis method and is detailed in [6]. The method is appropriate for sampling measures
absolutely continuous with respect to a Gaussian in the situation where it is straight-
forward to sample directly from the Gaussian itself. We control the error (i) simply
by running the MCMC method until time averages of various test statistics have con-
verged; the reader interested in the effect of this Monte Carlo error should consult [6].
The error in (ii) is precisely the error which we quantify in Theorem 3.4; for the partic-
ular experiments used here we predict an error of order N −q for any q ∈ (0, 1). In this
paper we have not analyzed the errors resulting from (iii) and (iv): these approxima-
tions are not included in the analysis leading to Theorem 3.4. However, we anticipate
that Theorem 2.2 or Theorem 2.4 could be used to study such approximations, and
the numerical evidence which follows below is consistent with this conjecture.
In the following three numerical experiments (each illustrated by a figure) we
study the effect of one or more of the approximations (ii), (iii), and (iv) on the
empirical distribution (“histogram”) found from marginalizing data from the MCMC
method onto the real part of the Fourier mode with wave vector k = (0, 1). Similar
results are found for other Fourier modes although it is important to note that at
high values of |k| the data are uninformative and the posterior is very close to the
prior (see [6] for details). The first two figures use J = 9 Lagrangian trajectories,
while the third uses J = 400. Figure 3.1 shows the effect of increasing the number
of Fourier modes2 used from 16, through 100 and 196, to a total of 400 modes and
illustrates Theorem 3.4 in that convergence to a limit is observed as the number of
Fourier modes increases. However, this experiment is conducted by using bilinear
interpolation of the velocity field on the grid, in order to obtain off-grid velocities
required for particle trajectories. At the cost of quadrupling the number of FFTs
it is possible to implement bicubic interpolation.3 Conducting the same refinement
of the number of Fourier modes then yields Figure 3.2. Comparison of Figures 3.1

2 Here by number of Fourier modes, we mean the dimension of the Fourier space approximation,

i.e., the number of grid points.

3 Bicubic interpolation with no added FFTs is also possible by using finite difference methods to

find the partial derivatives, but at a lower order of accuracy

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

336 S. L. COTTER, M. DASHTI, AND A. M. STUART

Fig. 3.1. Marginal distributions on Re(u0,1 (0)) with differing numbers of Fourier modes.

Fig. 3.2. Marginal distributions on Re(u0,1 (0)) with differing numbers of Fourier modes, bicu-
bic interpolation used.

and 3.2 shows that the approximation (iii) by increased order of interpolation leads
to improved approximation of the posterior distribution, and Figure 3.2 alone again
illustrates Theorem 3.4. Figure 3.3 shows the effect (iv) of reducing the time-step used
in the integration of the Lagrangian trajectories. Note that many more (400) particles
were used to generate the observations leading to this figure than were used in the
preceding two figures. This explains the quantitatively different posterior distribution;
in particular, the variance in the posterior distribution is considerably smaller because
more data are present. The result shows clearly that reducing the time-step leads to
convergence in the posterior distribution.
4. Eulerian data assimilation. In this section we consider a data assimilation
problem that is related to weather forecasting applications. In this problem, direct
observations are made of the velocity field of an incompressible viscous flow at some

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 337

Fig. 3.3. Marginal distributions on Re(u0,1 (0)) with differing time-step, Lagrangian data.

fixed points in space-time, the mathematical model is the two-dimensional Navier–

Stokes equations on a torus, and the objective is to obtain an estimate of the initial
velocity field. The spaces H and Hs are as defined in section 3, with % · %s the norm
in Hs and % · % = % · %0 . The definitions of A, the Stokes operator, and P , the Leray
projector, are also as in the previous section.
We write the Navier–Stokes equations (3.1) with ι = 1 as an ordinary differential
equation in H
dv
(4.1) + νAv + B(v, v) = ψ, v(0) = u,
dt
which is the same as the Stokes equations (3.5), up to the addition of the bilinear
form B(v, v). This term arises from projection of the nonlinear term under P . Our
aim is to determine u from noisy observations of the velocity field v at time t > 0 and
at points x1 , . . . , xK ∈ D:

yk = v(xk , t) + ηk , k = 1, . . . , K.

We assume that the noise is Gaussian and the ηk form an i.i.d. sequence with η1 ∼
N (0, γ 2 ). It is known (see Chapter 3 of [29], for example) that for u ∈ H and
f ∈ L2 (0, T ; Hs ) with s > 0 a unique solution to (4.1) exists which satisfies u ∈
L∞ (0, T ; H1+s ) ⊂ L∞ (0, T ; L∞(D)). Therefore for such initial condition and forcing
function the value of v at any x ∈ D and any t > 0 can be written as a function of u.
Hence, we can write

y = G(u) + η,

where y = (y1∗ , · · · , yK ) ∈ R2K and η = (η1∗ , . . . , ηk∗ )∗ ∈ R2K is distributed as

∗ ∗
2
N (0, γ I) and

(4.2) G(u) = (v(x1 , t)∗ , · · · , v(xK , t)∗ )∗ .

Now consider a Gaussian prior measure µ0 ∼ N (ub , βA−α ) with β > 0 and α > 1; the
second condition ensures that functions drawn from the prior are in H, by Lemma
A.4. In Theorem 3.4 of [5] it is shown that with such prior measure, the posterior
measure of the above inverse problem is well defined.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

338 S. L. COTTER, M. DASHTI, AND A. M. STUART

Theorem 4.1. Assume that f ∈ L2 (0, T, Hs ) with s > 0. Consider the Eulerian
data assimilation problem described above. Define a Gaussian measure µ0 on H, with
mean ub and covariance operator β A−α for any β > 0 and α > 1. If ub ∈ Hα , then
the probability measure µ(du) = P(du|y) is absolutely continuous with respect to µ0
with Radon–Nikodým derivative
$ %
dµ 1 2
(4.3) (u) ∝ exp − 2 |y − G(u)|Σ .
dµ0 2γ

We now define an approximation µN to µ given by (4.3). The approximation is

made by employing the Galerkin approximations of v to define an approximate G.
The Galerkin approximation of v, v N , is the solution of
dv N
(4.4) + νAv N + P N B(v N , v N ) = P N ψ, v N (0) = P N u,
dt
with P N as defined in the previous section. Let
! "T
G N (u) = v N (x1 , t), . . . , v N (xK , t) ,

and then consider the approximate prior measure µN defined via its Radon–Nikodým
derivative with respect to µ0 :
$ %
dµN 1 N 2
(4.5) ∝ exp − 2 |y − G (u)|Σ .
dµ0 2γ

Our aim is to show that µN converges to µ in the Hellinger metric. Unlike the examples
in the previous section, we are unable to obtain sufficient control on the dependence of
the error constant on u in the forward error bound to enable application of Theorem
2.2; hence we employ Theorem 2.4. In the following lemma we obtain a bound on
%v(t) − v N (t)%L∞ (D) and therefore on |G(u) − G N (u)|. Following the statement of the
lemma, we state and prove the basic approximation theorem for this section. The
proof of the lemma itself is given after the statement and proof of the approximation
theorem for the posterior probability measure.
Lemma 4.2. Let v N be the solution of the Galerkin system (4.4). For any t > t0

%v(t) − v N (t)%L∞ (D) " C(%u%, t0 ) ψ(N ),

where ψ(N ) → 0 as N → ∞.
The above lemma leads us to the following convergence result for µN .
Theorem 4.3. Let µN be defined according to (4.5) and let the assumptions of
Theorem 4.1 hold. Then

dHell (µ, µN ) → 0

as N → ∞.
Proof. We apply Theorem 2.4 with X = H. Assumption 2 (and hence Assumption
1) is established in Lemma 3.1 of [5]. By Lemma 4.2

|G(u) − G N (u)| " Kψ(N )

with K = K(%u%) and ψ(N ) → 0 as N → 0. Therefore the result follows by Theorem

2.4.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 339

Proof of Lemma 4.2. Let e1 = v − P N v and e2 = P N v − v N . Applying P N to

(4.1) yields

dP N v
+ νAP N v + P N B(v, v) = P N ψ.
dt
Therefore e2 = P N v − v N satisfies
de2
(4.6) + νAe2 = P N B(e1 + e2 , v) + P N B(v N , e1 + e2 ), e2 (0) = 0.
dt
Since for any l ! 0 and for m > l
1
(4.7) %e1 %2l " %v%2m ,
N 2(m−l)
we will obtain an upper bound for %e2 %1+l , l > 0, in terms of the Sobolev norms of
e1 and then use the embedding H1+l ⊂ L∞ to conclude the result of the lemma.
Taking the inner product of (4.6) with e2 , and noting that P N is self-adjoint,
N
P e2 = e2 , and (B(v, w), w) = 0, we obtain

1 d
%e2 %2 + ν%De2 %2 = (B(e1 + e2 , v), e2 ) + (B(v N , e1 ), e2 )
2 dt
1/2 1/2
" c%e1 %1/2 %e1 %1 %v%1 %e2 %1/2 %e2 %1 + c%e2 % %v%1 %e2 %1
1/2 1/2
+ c%v N %1/2 %v N %1 %e1 %1 %e2 %1/2 %e2 %1
" c%e1 %2 %e1 %21 + c%v%21 %e2 % + c%e2 %2 %v%21
ν
+ c%v N % %v N %1 %e1 %1 + c%e1 %1 %e2 % + %e2 %21 .
2
Therefore
d
(1 + %e2 %2 ) + ν %De2 %2
dt
" c (1 + %v%21 ) (1 + %e2 %2 ) + c(1 + %e1 %2 ) %e1 %21 + c %v N % %v N %1 %e1 %1 ,

which gives
# t $ # t %# t
%e2 (t)%2 + ν %De2 %2 " c β(t) 1 + %v N %2 %v N %21 dτ %e1 %21 dτ
0 0 0
# t
+ c β(t) (1 + %e1 %2 ) %e1 %21 dτ
0

with
$ # t %
2
β(t) = exp c 1 + %v%1 dτ .
0

Hence
# t # t
c+c$u$2
(4.8) 2
%e2 (t)% + ν %De2 % " c(1 + %u% ) e
2 4
(1 + %e1 %2 ) %e1 %21 dτ.
0 0

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

340 S. L. COTTER, M. DASHTI, AND A. M. STUART

To estimate %e2 (t)%s for s < 1, we take the inner product of (4.6) with As e2 ,
0 < s < 1, and write
1 d ! " ! "
%e2 %2s + ν%e2 %21+s " | ((e1 + e2 ) · ∇)v, As e2 | + | (v N · ∇)(e1 + e2 ), As e2 |.
2 dt
Using
! "
| (u · ∇)v, As w | " c%u%s %v%1 %w%1+s
and Young’s inequality we obtain
d
%e2 %2s + ν%e2 %21+s " c (%e1 %2s + %e2 %2s ) %v%21 + c %v N %2s (%e1 %21 + %e2 %21 ).
dt
Now integrating with respect to t over (t0 , t) with 0 < t0 < t we can write
# t # t
%e2 (t)%2s + ν %e2 %21+s dτ " %e2 (t0 )%2s + c sup %v(τ )%21 %e1 %2s + %e2 %2s dτ
t0 τ "t0 0
# t
+ c sup %v N (τ )%2s %e1 %21 + %e2 %21 dτ.
τ "t0 0

Therefore since for s " 1 and t ! t0

c(1 + %u%2 )
%v(t)%2s " ,
ts0
and noting that the same kind of decay bounds that hold for v can be shown similarly
for v N as well, we have
# t # t
c 2
%e2 (t)%2s + ν %e2 %21+s dτ " %e2 (t0 )%2s + (1 + %u%6 )ec+c$u$ (1 + %e1 %2 ) %e1 %21 dτ.
t0 t 0 0

Integrating the above inequality with respect to t0 in (0, t) we obtain

# t # t
c 2
(4.9) %e2 (t)%2s + ν %e2 %21+s dτ " (1 + %u%6 )ec+c$u$ (1 + %e1 %2 ) %e1 %21 dτ
t0 t0 0

for t > t0 .
Now we estimate %e2 (t)%s for s > 1. Taking the inner product of (4.6) with
A1+l e2 , 0 < l < 1, we obtain
1 d ! "
%e2 %21+l + ν%e2 %22+l " | ((e1 + e2 ) · ∇)v, A1+l e2 |
2 dt ! "
+ | (v N · ∇)(e1 + e2 ), A1+l e2 |.
Since (see [5])
! "
(u · ∇)v, A1+l w " c %u%1+l %v%1 %w%2+l + c %u%l %v%2 %w%2+l
and using Young’s inequality, we can write
d
%e2 %21+l + ν%e2 %22+l " c %e1 %21+l %v%21 + c %e1 %2l %v%22
dt
+ c %e2 %21+l %v%21 + c %e2 %2l %v%22
+ c %v N %21+l %e1 %21 + c %v N %2l %e1 %22
2/l
+ c %v N %21+l %e2 %21 + c %v N %l %e2 %21+l .

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 341

Now we integrate the above inequality with respect to t and over (t0 /2 + σ, t) with
0 < t0 < t and 0 < σ < t − t0 /2 and obtain (noting that %v N %s " %v%s for any s > 0)
# t
%e2 (t)%21+l " %e2 (t0 /2 + σ)%21+l + sup %v(τ )%21 %e1 %21+l + %e2 %21+l dτ
τ "t0 /2 t0 /2+σ
# t
+ sup (%e1 (τ )%2l + %e2 (τ )%2l ) %v%22 dτ
τ "t0 /2 t0 /2+σ
# t
+ sup (%e1 (τ )%21 + %e2 (τ )%21 ) %v N %21+l dτ
τ "t0 /2 t0 /2+σ
# t
N 2/l
+ sup (1 + %v (τ )%l ) %e1 %22 + %e2 %21+l dτ.
τ "t0 /2 t0 /2+σ

We have, for s > 1 and t > t0 [5],

c(1 + %u%4 )
%v(t)%2s " .
ts0
Therefore using (4.9) and (4.7) we conclude that
%e2 (t)%21+l " %e2 (t0 /2 + σ)%21+l
$ # t
1 1
+ Cp (%u%) + (1 + %e1 %2 ) %e1 %21 dτ
N 2(m−l) t1+m
0 t 1+l
0 0
%
1
+ 2(r−1) 1+r
N t0
with r > 1 and where Cp (%u%) is a constant depending on polynomials of %u%. Inte-
grating the above inequality with respect to σ over (0, t − t0 /2) we obtain
$ %# t
1 1
%e2 (t)%21+l " Cp (%u%) 1+l + 2+l (1 + %e1 %2 ) %e1 %21 dτ
t0 t0 0
$ %
1 1
+ Cp (%u%) + .
N 2(m−l) t2+m
0 N 2(r−1) t2+r
0
/t
Now to show that %e1 %2 + 0 %e1 %21 dτ → 0 as N → ∞, we note that e1 satisfies
1 d
%e1 %2 + ν%De1 %2 " %(I − PN )f % %e1% + %(B(v, v), e1 )%
2 dt
" %(I − PN )f % %e1% + %v%1/2 %Dv%3/2 %e1 %1/2 %De1 %1/2
ν
" %(I − PN )f % %e1% + c %v%2/3 %Dv%2 %e1 %2 + %De1 %2 .
2
Therefore
d
%e1 %2 + ν%De1 %2 " c %(I − PN )f %2 + c (1 + %v%2/3 %Dv%2 ) %e1 %2 ,
dt
and after integrating, we get
# T 0 # T 1
2
%e1 % + %e1 %1 dτ " exp(1 + Cp (%u%)) %e1 (0)% +
2 2 N
%(I − P )f % 2
dτ.
0 0

Since f ∈ L2 (0, T ; H), the above integral tends to zero as N → ∞ and the result
follows.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

342 S. L. COTTER, M. DASHTI, AND A. M. STUART

5. Conclusions. In this paper we have studied the approximation of inverse

problems which have been regularized by means of a Bayesian formulation. We have
developed a general approximation theory which allows for the transfer of approx-
imation results for the forward problem into approximation results for the inverse
problem. The theory clearly separates analysis of the forward problem, in which no
probabilistic methods are required, and the probabilistic framework for the inverse
problem itself: it is simply necessary that the requisite bounds and approximation
properties for the forward problem hold in a space with full measure under the prior.
Indeed, the spaces in which the forward approximation theory is developed may be
seen to place constraints on the prior, in order to ensure the desired robustness of the
inverse problem to finite dimensional approximation.
In applications there are two sources of error when calculating expectations of
functions of infinite dimensional random variables: the error for which we provide
an analysis in this paper, namely, the approximation of the measure itself in a finite
dimensional subspace, together with the error incurred through calculation of expec-
tations. The two sources of error must be balanced in order to optimize computational
cost and thus the analysis in this paper is of potential practical value to practitioners
in estimating one source of error. The expectations themselves can be computed by
MCMC methods [17], by quasi–Monte Carlo methods [24], or by the more recently
analyzed polynomial chaos methods [25, 26, 23, 31]. In all these cases there are
quantitative error estimates associated with the calculation of an expectation.
We have studied two specific applications, both concerned with determining the
initial condition of a dissipative PDE, from observations of various kinds, at positive
times. However, the general approach is applicable to a range of inverse problems for
functions when formulated in a Bayesian fashion. The article [27] overviews many
applications from this point of view. Furthermore we have limited our approximation
of the underlying forward problem to spectral methods. However, we anticipate that
the general approach will be useful for the analysis of other spatial approximations
based on finite element methods, for example, and to approximation errors resulting
from time discretization; indeed, it would be interesting to carry out analyses for such
approximations.
Appendix A. Here we collect together some basic facts in analysis and proba-
bility. First we state the well-known Gronwall inequality in the form in which we will
use it.4
Lemma A.1. Let I = [c, d) with d ∈ (c, ∞]. Assume / that α, u ∈ C(I; R+ ) and
that there is λ < ∞ such that, for all intervals J ⊆ I, J β(s)ds < λ. If
# t
u(t) " α(t) + β(s)u(s)ds, t ∈ I,
c

then
# t $# t %
u(t) " α(t) + α(s)β(s) exp β(r)dr ds, t ∈ I.
c s

In particular, if α(t) = u + 2at is positive in I and β(t) = 2b, then

a
u(t) " exp(2bt)u + (exp(2bt) − 1), t ∈ I.
b
4 See http://en.wikipedia.org/wiki/Gronwall’s inequality.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 343

Finally, if c = 0, and 0 < α(t) " K in I, then

u(t) " K + Kλ exp(λ), t ∈ I.
Now consider the Hilbert-space valued ODE
dv
(A.1) + Av = f, v(0) = u.
dt
We state some basic results in this area, provable by use of the techniques in [20], for
example, or by direct calculation using the eigenbasis for A. For f = 0 the solution
v ∈ C([0, ∞), H) ∩ C 1 ((0, ∞), D(A)) and
(A.2) %v%2s " Ct−(s−l) %u%2l ∀t ∈ (0, T ] s ! l ! 0.
If f ∈ C([0, T ], Hγ ) for some γ ! 0, then (A.1) has a unique mild solution u ∈
C([0, T ]; H) and, for 0 " . < γ + 2,
& %u% '
l
(A.3) %v(t)%s " C (s−l)/2 + %f %C([0,T ],Hγ )
t
for s ∈ [., 2 + γ).
It is central to this paper to estimate the distance between two probability mea-
sures. Assume that we have two probability measures µ and µ% , both absolutely
continuous with respect to the same reference measure ν. The Hellinger distance
between µ and µ% is
2
& 1 # &3 dµ 3 dµ% '2 '
%
dHell (µ, µ ) = − dν .
2 dν dν
The Hellinger distance is particularly useful for estimating the difference between
expectation values of functions of random variables under different measures. This is
illustrated in the following lemma with straightforward proof.
Lemma A.2. Assume that two measures µ and µ% on a Banach space (X, % · %X )
are both absolutely continuous with respect to a measure ν. Assume also that f : X →
Z, where (Z, % · %) is a Banach space, has second moments with respect to both µ and
µ% . Then
"
& "
' 12
%Eµ f − Eµ f % " 2 Eµ %f %2 + Eµ %f %2 dHell (µ, µ% ).

Furthermore, if (Z, +·, ·,) is a Hilbert space and f : X → Z has fourth moments, then
"
& "
' 12
%Eµ f ⊗ f − Eµ f ⊗ f % " 2 Eµ %f %4 + Eµ %f %4 dHell (µ, µ% ).

Note, in particular, that choosing X = Z, and with f chosen to be the identity

mapping, we deduce that the differences in mean and covariance operators under two
measures are bounded above by the Hellinger distance between the two measures.
The following Fernique theorem (see [22], Theorem 2.6) will be used repeatedly.
Theorem A.3. Let x ∼ µ = N (0, C) where µ is a Gaussian measure on Hilbert
space H. Assume that µ0 (X) = 1 for some Banach space (X, % · %X ) with X ⊆ H.
Then there exists α > 0 such that
#
! "
exp α%x%2X µ(dx) < ∞.
H

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

344 S. L. COTTER, M. DASHTI, AND A. M. STUART

The following regularity properties of Gaussian random fields will be useful to

us; the results may be proved by use of the Kolmogorov continuity criterion, together
with the Karhunen–Loeve expansion (see [22], section 3.2).
Lemma A.4. Consider a Gaussian measure µ = N (0, C) with C = βA−α where
A is as defined earlier in this Appendix A. Then u ∼ µ is almost surely s-Hölder
continuous for any exponent s < min{1, α − d2 } and u ∈ Hs , µ—almost surely, for
any s < α − d2 .

REFERENCES

[1] A. Apte, C. K. R. T. Jones, A. M. Stuart, and J. Voss, Data assimilation: Mathematical

and statistical perspectives, Internat. J. Numer. Methods Fluids, 56 (2008), pp. 1033–1046.
[2] S. R. Arridge, J. P. Kaipio, V. Kolehmainen, M. Schweiger, E. Somersalo, and T. Tar-
vainen, Approximation errors and model reduction with an application to optical diffusion
tomography, Inverse Problems, 22 (2006), pp. 175–195.
[3] V. I. Bogachev, Gaussian Measures, American Mathematical Society, Providence, RI, 1998.
[4] D. Calvetti and E. Somersalo, Introduction to bayesian scientific computing, Surveys and
Tutorials in the Applied Mathematical Sciences 2, Springer, New York, 2007.
[5] S. L. Cotter, M. Dashti, J. C. Robinson, and A. M. Stuart, Bayesian inverse problems
for functions and applications to fluid mechanics, Inverse Problems, 25 (2010), p. 115008.
[6] S. L. Cotter, M. Dashti, J. C. Robinson, and A. M. Stuart, MCMC Methods on Function
Space and Applications to Fluid Mechanics, in preparation, 2010.
[7] M. Dashti and J. C. Robinson, A simple proof of uniqueness of the particle trajectories for
solutions of the Navier-Stokes equations, Nonlinearity, 22 (2009), pp. 735–746.
[8] H. K. Engl, M. Hanke, and A. Neubauer, Regularization of inverse problems, Kluwer,
Dordrecht, the Netherlands, 1996.
[9] J. N. Franklin, Well-posed stochastic extensions of ill posed linear problems, J. Math. Anal.
Appl., 31 (1970), pp. 682–716.
[10] A. Hofinger and H. K. Pikkarainen, Convergence rates for the Bayesian approach to linear
inverse problems, Inverse Problems, 23 (2007), pp. 2469–2484.
[11] K. Ide, M. Courier, M. Ghil, and A. Lorenc, Unified notation for assimilation: Operational,
sequential and variational, J. Met. Soc. Japan, 75 (1997), pp. 181–189.
[12] J. Kaipio and E. Somersalo, Statistical and computational inverse problems, Applied Math-
ematical Sciences 160, Springer, New York, 2005.
[13] J. P. Kaipio and E. Somersalo, Approximation errors in nonstationary inverse problems,
Inverse Prob. Imag., 1 (2007), pp. 77–93.
[14] J. P. Kaipio and E. Somersalo, Statistical inverse problems: Discretization, model reduction
and inverse crimes, J. Comput. Appl. Math., 198 (2007), pp. 493–504.
[15] A. Kirsch, An Introduction to the Mathematical Theory of Inverse Problems, Springer, New
York, 1996.
[16] M. S. Lehtinen, L. Paivarinta, and E. Somersalo, Linear inverse problems for generalized
random variables, Inverse Problems, 5 (1989), pp. 599–612.
[17] J. Liu, Monte Carlo strategies in scientific computing, Springer Texts in Statistics, Springer,
New York, 2001.
[18] A. Neubauer and H. K. Pikkarainen, Convergence results for the Bayesian inversion theory,
J. Inverse Ill-Posed Probl., 16 (2008), pp. 601–613.
[19] N. K. Nichols, Data assimilation: Aims and basic concepts, Data Assimilation for the
Earth System, R. Swinbank, V. Shutyaev, and W. A. Lahoz, eds., Kluwer, Dordrecht, the
Netherlands, 2003, pp. 9–20.
[20] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations,
Springer, New York, 1983.
[21] H. K. Pikkarainen, State estimation approach to nonstationary inverse problems: Discretiza-
tion error and filtering problem, Inverse Problems, 22 (2006), pp. 365–379.
[22] G. Da Prato and J. Zabczyk, Stochastic equations in infinite dimensions, Encyclopedia of
Mathematics and its Applications 44, Cambridge University Press, Cambridge, 1992.
[23] C. Schwab and R. A. Todor, Karhunen-Loeve approximation of random fields in domains
by generalized fast multipole methods, J. Comput. Phys., 217 (2006), pp. 100–122.
[24] I. Sloan and H. Wozniakowski, When are quasi-Monte Carlo algorithms efficient for high-
dimensional integrals? J. Complexity, 14 (1998), pp. 1–33.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

 APPROX. OF BAYESIAN INVERSE PROBLEMS FOR PDES 345

[25] P. D. Spanos and R. Ghanem, Stochastic finite element expansion for random media, J. Eng.
Mech., 115 (1989), pp. 1035–1053.
[26] P. D. Spanos and R. Ghanem, Stochastic Finite Elements: A Spectral Approach, Dover, New
York, 2003.
[27] A. M. Stuart, Inverse problems: A Bayesian approach, Acta Numer., 19 (2010), to appear.
[28] A. Tarantola, Inverse Problem Theory. SIAM, Philadelphia, PA, 2005.
[29] R. Temam, Navier-Stokes equations and nonlinear functional analysis, CBMS-NSF Regional
Conference Series in Applied Mathematics 66, SIAM, Philadelphia, PA, 1995.
[30] R. Temam, Navier-Stokes equations. AMS Chelsea, Providence, RI, 2001.
[31] R. A. Todor and C. Schwab, Convergence rates for sparse chaos approximations of elliptic
problems with stochastic coefficients, IMA J. Numer. Anal., 27 (2007), pp. 232–261.
[32] C. R. Vogel, Computational Methods for Inverse Problems, SIAM, Philadelphia, PA, 2002.
[33] L. W. White, A study of uniqueness for the initialiazation problem for burgers’ equation, J.
Math. Anal. Appl., 172 (1993), pp. 412–431.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.