POSTERIOR CONSISTENCY AND CONVERGENCE
RATES FOR BAYESIAN INVERSION WITH
HYPOELLIPTIC OPERATORS

arXiv:1507.01772v2 [] 24 May 2016

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

Abstract. Bayesian approach to inverse problems is studied in the
case where the forward map is a linear hypoelliptic pseudodifferential
operator and measurement error is additive white Gaussian noise. The
measurement model for an unknown Gaussian random variable U (x, ω)
is
Mδ (y, ω) = A(U (x, ω)) + δ E(y, ω),
where A is a finitely many orders smoothing linear hypoelliptic operator
and δ > 0 is the noise magnitude. The covariance operator CU of U is
smoothing of order 2r, self-adjoint, injective and elliptic pseudodifferential operator.
If E was taking values in L2 then in Gaussian case solving the conditional mean (and maximum a posteriori) estimate is linked to solving
the minimisation problem

−1/2
Tδ (mδ ) = arg min kAu − mδ k2L2 + δ 2 kCU uk2L2 .
u∈H r

However, Gaussian white noise does not take values in L2 but in H −s
where s > 0 is big enough. A modification of the above approach to
solve the inverse problem is presented, covering the case of white Gaussian measurement noise. Furthermore, the convergence of conditional
mean estimate to the correct solution as δ → 0 is proven in appropriate
function spaces using microlocal analysis. Also the frequentist posterior
contractions rates are studied.

Keywords: Posterior consistency, convergence rate, Bayesian inverse problem, white noise, pseudodifferential operator
Contents
1. Introduction
Notations
2. Convergence results
2.1. Main result
3. Generalised random variables
3.1. White noise paradox
3.2. The smoothness of the prior
4. Proof of the main Theorems
5. Posterior distribution and confidence regions
5.1. Contraction of the posterior distribution

2
4
4
6
10
11
11
13
19
21

Department of Mathematics and Statistics P.O. Box 68, 00014 University
of Helsinki FINLAND, e-mails: hanne.kekkonen@helsinki.fi, matti.lassas@helsinki.fi,
samuli.siltanen@helsinki.fi.
1

2

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

5.2. Convergence of the posterior distribution in Bayesian settings
5.3. Discussion
Appendix A. Some examples of hypoelliptic operators
Appendix B. Computational example
References

24
27
27
29
31

1. Introduction
Practical inverse problems arise from the need to extract information from
indirect data. For example, consider a device designed for measuring point
values of a physical quantity u(x). Technological imperfections cause the
values of u at nearby points to merge together in the measurement. Mathematically this corresponds to convolution Φ ∗ u by a point spread function
Φ. The inverse problem is to recover the function u approximately from a
finite number of point values of Φ ∗ u corrupted by random white noise.
Computational inversion requires a finite representation of the quantity
u(x). In this paper we promote the view that it is a good idea to design a
continuous model for u which then can be discretized with a desired number
of degrees of freedom. This paves the way for the analysis of convergence
as the discretization becomes finer. Such convergence enables switching
between different discretizations consistently; this is crucial for multigrid
methods and for certain parameter-choice strategies.
In ill-posed inverse problems the measurement data alone is not sufficient
for noise-robust recovery of the quantity of interest. For instance, Fourier
b u, so frequency-domain information is lost in areas
transforming Φ∗u gives Φb
b
where Φ is very close to zero. Therefore, successful computational inversion
requires some a priori information in addition to the measurement data.
Practical inversion is all about combining measurement data and a priori
information in a noise-robust way. The classical approach to do this is regularization that assumes that the noise is deterministic and small in norm.
Regularization involves defining a family of continuous maps, parametrized
by the norm knoisekL2 , from the data space to the space of unknown quantities. This must be done so that as knoisekL2 → 0, the reconstruction
approaches the true solution along a stable path. This methodology was
originated by Tikhonov [53, 54]. Both continuous and discrete cases have
been studied in depth in [9, 16, 24, 40, 38, 55, 47].
There is a serious drawback in the above noise model in the continuous
limit. Namely, continuous white noise on Rd is not square integrable. We
discuss this in detail below in Section 3.1. The goal of this paper is to use
Bayesian inversion to construct a consistent continuous-discrete framework
covering the case of white noise.
Bayesian inversion is a flexible framework for combining measurement
data and a priori information in the form of a posterior distribution [6, 10,
22, 33, 50]. Computational exploration of the finite-dimensional posterior
distribution yields useful estimates of the quantity of interest and enables

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

3

uncertainty quantification. Furthermore, analytic results about the continuous model can then be restricted to a given resolution in a discretizationinvariant way.
Our approach to Bayesian inversion follows a general strategy of computational mathematics: we consider a continuous model which can be discretized for any practical setting.
We study the following continuous model for indirect measurements:
(1.1)

Mδ = AU + Eδ,

where the random variables Mδ (data) and U (quantity of interest) take
values in the Sobolev spaces H −s (N ) and H τ (N ), respectively. Here N
is a d-dimensional compact manifold e.g. a torus corresponding to a ddimensional cube with opposite sides glued together. The real parameter τ
is related to our a priori information about the smoothness of the unknown
quantity of interest.
The measurement operator A in our model (1.1) is quite general: we
assume it to be a finitely smoothing, injective hypoelliptic pseudodifferential
operator (ΨDO). See section 4 for precise definition. This class includes
convolution operators with finitely smooth kernel Φ. One example of an
operator that is hypoelliptic but not elliptic is the heat operator. For more
examples of hypoelliptic operators see Appendix A. The measurement noise
E is assumed to be normalised white Gaussian noise with mean zero and
unit variance, and δ > 0 models the noise amplitude.
We model practical measurement data by
(1.2)

Mk = Pk (AU ) + Pk (E)δ.

Here Pk is a linear operator related to measurement device; we assume that
Pk is an orthogonal projection with k-dimensional range. We discretize the
unknown U using some computationally feasible approximation of the form
Un = Sn U . Now we can study an inverse problem
(1.3)

given a realisation of Mk , estimate Un .

We are interested to know what happens to the approximated solutions
of (1.1) when δ → 0. The analysis of small noise limit, also known as
the theory of posterior consistency, has attracted a lot of interest in the
last decade. Posterior convergence rates were first studied in [14, 48] and
further in papers [1, 2, 7, 18, 21, 25, 26, 28, 42, 48, 51, 58]. However, much
remains to be done. Developing a comprehensive theory is important since
posterior consistency justifies the use of the Bayesian approach the same
way as convergence results justify the use of regularisation techniques.
In the above mentioned papers the problem is studied from the frequentist
point of view, that is, the data is thought to be generated by a fixed ’true’
solution u† instead of random draw U (ω) from the prior distribution. This
means that all the randomness in M comes from the randomness of the
noise E. The interest is then on the contraction of the posterior distribution
around the ’true’ solution u† as the noise goes to zero, see Subsection 2.1.2
and Theorem 3. The main emphasis of this paper, however, is in the purely
Bayesian approach where it is assumed that also U is random. Since U and

4

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

E are assumed independent we can write
(1.4)

Mδ (ω) = AU (ω1 ) + E(ω2 )δ,

where ω = (ω1 , ω2 ) ∈ Ω1 × Ω2 . In Bayesian case the posterior distribution
is a function of (ω1 , ω2 ). Also the probability measure dP can be written in
the following form
dP = dP1 (ω1 )dP2 (ω2 ).
We will denote the expected value over the joint distribution of U and E by
E. The expected value over the noise is defined by
Z
(1.5)
F (ω1 , ω2 )dP2 (ω2 ).
EU (F (ω1 , ω2 )) =
Ω2

Our paper provides a conceptual advantage over much of the existing literature. In many earlier studies A and CU are perturbations of negative
powers of operator (I − ∆). Our assumption, formulated in terms of hypoelliptic operators, means roughly speaking that the measurement operator A
and the covariance operator CU do not need to have a common basis in their
singular value decomposition.
The rest of this paper is organised as follows. In section 2 we will introduce the Bayesian setting we are using and present our main result about
convergence rates. In section 3 we take a closer look to the generalised
Gaussian random variables in Sobolev spaces and introduce so called white
noise paradox. We will also show that the distribution U takes values in
H τ , where τ is related to the smoothness of the solution and depends on the
dimension of the space and the covariance of the prior. In section 4 we will
introduce hypoelliptic operators and prove Theorems 1 and 2. In section 5
we characterise credible sets and frequentist confidence regions and present
and prove two theorems about the contraction of them. In Appendix A we
will give examples of some hypoelliptic operators and in Appendix B we give
a computational example.
Notations.
Sm
Ψm
HΨm,m0
Ψm
p
0
HΨm,m
p

T rH q →H q C

Class of pseudodifferential symbols of order m. See Definition 1.
Space of pseudodifferential operators (ΨDO) of order m.
See Definition 2.
Space of hypoelliptic ΨDOs of type (m, m0 ). See Definition
3.
Space of ΨDOs of order m depending on spectral variable
with order p. See Definition 4.
Space of hypoelliptic ΨDOs of type (m, m0 ) depending on
spectral variable with order p. See Definition 5.
The trace of the operator C : H q → H q . See (3.8).
2. Convergence results

Let us return to our indirect measurement problem
Mδ = AU + Eδ,

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

5

where we model U = U (x, ω), Mδ = Mδ (y, ω) and E = E(y, ω) as random
functions. Here ω ∈ Ω is an element of a complete probability space (Ω, Σ, P)
and x and y denote the variables in domains of Euclidean spaces. The reason
why we model also U as a random variable is that even though the unknown
quantity is assumed to be deterministic we have only incomplete data of it.
All information available about U before performing the measurements is
included in a prior distribution that is independent of the measurement.
The Bayesian inversion theory is based on the Bayes formula. To solve
the inverse problem (1.3) we have to express available a priori information
of Un in the form of a prior distribution πpr in an n-dimensional subspace.
Let Mk and Ek be random vectors taking values in Rk , and denote their distributions by πMk and πEk , respectively. The solution of the inverse problem
after performing the measurements is the posterior distribution of the unknown random variable. Given a realisation of the discrete measurement
the posterior density for Un taking values in the n-dimensional subspace is
given by the Bayes formula
πpr (u)πEk (mδ | u)
π(u | mδ ) =
πMk (mδ )
(2.1)


1
= cπpr (u) exp − 2 kmδ − Auk22 , u ∈ Rn , mδ ∈ Rk
2δ
where A = Pk ASn is a k × n matrix approximation to the operator A.
An approximated solution for the inverse problem is often given as a point
estimate for (2.1). Let us assume that also Un has Gaussian distribution.
The maximum a posteriori (MAP) estimate TδM AP : Rk → Rn is defined by
(2.2)

TδM AP (Mk (ω)) := arg maxn π(u | Mk (ω)).
u∈R

Note that the MAP estimate depends on ω through the realisation of the
noise Ek (ω) and unknown Un (ω). When Un is Gaussian distributed the MAP
estimate coincides almost surely with the conditional mean estimate (CM)
(2.3)

TδCM (Mk (ω)) = E(Un |Mk )(ω)

a.s.

where Mk is the σ-algebra generated by Mk .
Since in our case CM=MAP a.s. we will consider below the MAP estimate. Let us denote the covariance matrix of Un by CUn . Solving the
maximisation problem (2.2) with a fixed realisation of noise and unknown
corresponds to solving the minimisation problem
 1
1 −1/2
(2.4)
Tδ (mδ ) = arg minn
kAu − mδ k22 + kCUn uk22 .
2
u∈R
2δ
2
Constructing Sn and πpr is the core difficulty in Bayesian inversion. In
many inverse problems there is no natural discretisation for the continuum
quantity U , so n can be freely chosen. Consequently, Sn and πpr should in
principle be described for all n > 0. This raises the following questions: do
the chosen Sn and πpr represent the same a priori knowledge consistently at
all resolutions n > 0? Does the estimate Tδ (mδ ) converge as n → ∞? See
e.g. [29, 34] Also, the number of data points may change, for example due to
an updated measurement device. The aim of this paper is to build a rigorous

6

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

theory that allows us to connect discrete models to their infinite-dimensional
limit models in a consistent way.
We achieve consistent representation of a priori knowledge by constructing the prior distribution for U in the infinite-dimensional space X. Then
the random variable Un = Sn U takes values in the finite-dimensional subspace Xn ⊂ X and represents approximately the same a priori knowledge
as U . The same way we construct distributions for M and E in the infinitedimensional space Y in which case the random variables Mk and Ek take
values in the finite-dimensional subspace Yk ⊂ Y .
The finite-dimensional problem (2.4) Γ-converges as n, k → ∞, under
certain assumptions (including that m should be an L2 -function), to the
following infinite-dimensional minimisation problem in a Sobolev space H r :
 1
1 −1/2
kmδ − Auk2L2 (N ) + kCU uk2L2 (N ) .
(2.5)
arg min
2
r
2δ
2
u∈H
−1/2

−1/2

Above CU
∈ Ψr , that is, CU
is −r orders smoothing pseudodifferential
operator. See [23] for a proof. If we are thinking the above as a MAP
estimate to a Bayesian problem we have to assume that U has formally the
following distribution


1 −1/2 2
πpr (u) = c exp − kCU ukL2 (N ) .
formally
2
Formula (2.5) only makes sense if the noise is square integrable. Even though
kεk k`2 < ∞
with any k ∈ N the limit, when k → ∞, is infinity. We will return to this so
called ‘white noise paradox’ in section 3.1.
2.1. Main result. Let us now modify formula (2.5) to arrive at something
useful for white Gaussian noise. When ε ∈ L2 we can write
(2.6)

kmδ − Auk2L2 (N ) = kAuk2L2 (N ) − 2(mδ , Au)L2 (N ) + kmδ k2L2 (N ) .

Now omitting the infinite ‘constant term’ kmδ k2L2 (N ) we get a new minimisation problem which is well-defined also when m is not an L2 function

−1/2
(2.7)
Tδ (mδ ) := arg min kAuk2L2 (N ) − 2hmδ , Aui + δ 2 kCU uk2L ,
u∈H r (N )

where hmδ , Aui is interpreted as a suitable duality pairing instead of L2 (N )
inner product. When A ∈ Ψ−t , t ≥ −τ + s, we can define hmδ , Aui =
hmδ , AuiH −s (N )×H s (N ) . Note that the forward operator A, the prior distribution and the noise depend on on each other only through assumption
t ≥ −τ + s.
It is well-known that the solution of the finite-dimensional problem (2.4)
can be calculated using the following formula:
(2.8)

−1 T
Tδ (mδ ) = (ATA + δ 2 C−1
Un ) A m δ .

We can write the approximated solution uδ := Tδ (mδ ) of the continuous
problem (2.7) by
(2.9)

Tδ (mδ ) = (A∗ A + δ 2 CU−1 )−1 A∗ mδ .

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

7

Before the main result of the paper we will study a simple example to
give a reader an insight to Bayesian settings.
Example 1. Let N be a 1-dimensional torus T1 . We are interested of the
inverse problem
Mδ = AU + δE
where we assume that E ∼ N (0, I) and U ∼ N (0, I), that is both the noise
and the unknown are assumed to be normalised white Gaussian noise, see
section 3 for rigorous definition. The white noise takes values in H τ with
some τ < −1/2. On the other hand white noise has formally the following
distribution


1
2
πpr (u) = c exp − kukL2 (T1 ) .
formally
2
Hence we want to solve

arg min kAuk2L2 − 2hmδ , Aui + δ 2 kuk2L2 (T1 ) .
u∈L2

Note that we are looking for a solution in L2 (T1 ) even though the realisations
of U are in L2 (T1 ) with probability zero. In general if we are interested in
finding a solution in H r (N ) then we can show that the prior should take
values in H τ (N ) where τ = r − s, see section 3.
2.1.1. Convergence results in Bayesian setting. We will now formulate the
main theorem of this paper about the convergence of the continuous solution
(2.9). The precise definitions are discussed in more detail in section 3 and
Theorem 1 is proved in section 4.
Theorem 1. Let r, s ∈ [0, ∞) and N be a d-dimensional closed manifold.
Let U (x, ω) be a generalised Gaussian random function taking values in
H τ (N ), τ = r − s, with zero mean and covariance operator CU . Assume
that the operator CU is a self-adjoint, injective and elliptic pseudodifferential
operator (ΨDO) of order −2r. Let E(y, ω) be white Gaussian noise on N .
Consider the measurement
Mδ (y, ω) = A(U (·, ω)) + δE(y, ω),

ω ∈ Ω,

where A ∈ HΨ−t,−t0 , t > max{0, −τ +s} and t ≤ t0 < 2t+r, is a hypoelliptic
pseudodifferential operator on the manifold N and A : L2 (N ) → L2 (N ) is
injective. Above δ ∈ R+ is the noise level and E takes values in H −s (N )
with some s > d/2.
Take ζ < τ − 3(t0 − t). Then we have the following convergence
(2.10)

EkUδ (ω) − U (ω)kH ζ → 0,

as δ → 0.

The expectation is taken with respect to the joint distribution of (U, E). We
have the following estimates for the speed of convergence:
(i) If ζ ≤ t − s − 2t0 then there is such C > 0 independent of δ that
(2.11)

EkUδ (ω) − U (ω)kH ζ ≤ Cδ

2t−t0 +r
t0 +r

.

8

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

(ii) If t − s − 2t0 ≤ ζ < τ − 3(t0 − t) then there is such C > 0 independent
of δ that
(2.12)

EkUδ (ω) − U (ω)kH ζ ≤ Cδ

−

ζ−τ +3(t0 −t)
t0 +r

.

The different convergence speeds (i) and (ii) show the trade-off between
the smoothness of the space and the speed of convergence. In case (i) we
get better convergence rates but in case (ii) we can use a stronger norm. We
see also that the smoother the forward operator A is the worse convergence
rates we get.
We note that instead of the estimates (2.11) or (2.12), we could alternatively take the expected value only with respect to the noise in which case
the constant C would depend on the realisation of U (ω). That is, proof of
Theorem 1 also shows that, we have almost surely
(2.13)

lim sup
δ→0

EU kUδ (ω) − U (ω)kH ζ
<∞
δν

ζ−τ +3(t0 −t)
0 +r
where ν = 2t−t
when t − s −
t0 +r when ζ ≤ t − s − 2t0 and ν = −
t0 +r
2t0 ≤ ζ < τ − 3(t0 − t).

Remark 1. The MAP-estimate Uδ takes values in the Cameron-Martin
space of the prior. The Cameron-Martin space is the intersection of all
linear subspaces where the random variable U belongs with probability one,
and since there may be uncountable many such linear spaces, the CameronMartin space may be a zero measurable subset of the space where U takes
values, see [3]. In the above settings where U ∼ N (0, CU ) and CU is of order
−2r the random variable U takes values in H τ , τ = r − s where s > d/2,
and the Cameron-Martin space containing the MAP-estimate is H r . In the
Bayesian setting it is natural that the MAP-estimate can not converge in a
smaller space than the one U takes values. However, the same behaviour
can be seen also in the deterministic setting when the unknown is in H r and
the MAP-estimate is thought as a Tikhonov regularised solution, see [23].
Example 2. Let us study a simple example in two dimensional torus T2 .
We consider a problem
Mδ = (I − ∆)−1 U + δE
where E is normalised white Gaussian noise that takes values in H −s (T2 ),
s > 1 and δ > 0 is the noise amplitude. The model operator A = (I − ∆)−1
is elliptic operator, smoothing of order 2. Let us consider the case when U
has a priori distribution N (0, CU ) where CU = (I − ∆)−2 , that is, r = 2.
Then U takes values in H τ (T2 ), where τ = r − s < 1, almost surely and
Uδ ∈ H 2 . Theorem 1 guarantees us convergence rate Cδ when ζ ≤ −3 and
(2.12) when −3 ≤ ζ ≤ 1. For example we get the following convergence in
L2 (T2 )
1

EkUδ (ω) − U (ω)kL2 (T2 ) ≤ Cδ 4 −
with  > 0 arbitrarily small.

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

9

2.1.2. Convergence results in frequentists setting. In the frequentist case one
is often interested in the model
(2.14)

Mδ† (ω) = A(u† ) + δE(ω)

where the data Mδ† is generated by a ‘true’ solution u† ∈ H τ (N ). Above E
is normalised white Gaussian noise and δ ∈ R+ is the noise level. In (2.14)
all the randomness of the Mδ† comes from the randomness of the noise E.
We denote
(2.15)

Uδ† (ω) = Tδ (Mδ† (ω)).

In the frequentist setting we consider the case where u† is an arbitrary element of the space H τ (N ) where the random variable U takes values, instead
of considering almost every element. Note that even though a set is measuretheoretically large, that is, has probability 1, it can sill be topologically small
(meager, or a set of first Baire category). For a discussion on these issues
see [13].
We can then study the convergence Eu† kUδ† (ω) − u† kH ζ or in more frequentist spirit Eu† kUδ† (ω) − u† k2H ζ where we use notation
Z
†
†
(2.16)
F (Tδ (Au† + δE(ω)), u† , ω)dP(ω),
Eu† F (Uδ (ω), u , ω) =
Ω

that is, the expectation Eu† is taken with respect to the noise E(ω) and the
other terms depending on ω, and u† is considered as a fixed parameter, c.f.
(1.5). This means that after computing the estimator Uδ† using Bayesian
methods we will consider the convergence of the estimator Uδ† to the ‘true’
solution u† which is not thought to be a random draw from the prior any
more.
Remark 2. Next we consider the frequentist case when in addition it is
assumed that u† ∈ H τ (N ) where τ ≥ 0. Note that then r ≥ d2 . The mean
b is defined
integrated squared error (MISE) of an estimator U
b , u† ) = Eu† kU
b − u† k2 2 .
(2.17)
R(U
L

The minimax risk rδ (H τ (N )) on the Sobolev space H τ (N ) is then given by
b , u† )
rδ (H τ (N )) = inf sup R(U
b u† ∈H τ (N )
U

b = g(M † ) where
where the infimum is taken over all estimators of the form U
δ
g ∈ B(H −s , H τ ). Here B(H −s , H τ ) is the set of Borel measurable functions
from H −s to H τ .
Theorem 2. Let r > s > d/2 and N be a d-dimensional closed manifold.
Let u† ∈ H τ (N ) where τ = r − s > 0. Assume that CU , the covariance
operator of the Gaussian prior, is a self-adjoint, injective and elliptic pseudodifferential operator of order −2r. Let E(y, ω) be white Gaussian noise on
N . Consider the measurement
Mδ† (y, ω) = A(u† (·)) + δE(y, ω),

ω ∈ Ω,

where A ∈ HΨ−t,−t0 , t > max{0, −τ + s} and t ≤ t0 ≤ t + τ /3, is a
hypoelliptic pseudodifferential operator on the manifold N and A : L2 (N ) →

10

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

L2 (N ) is injective. Above δ ∈ R+ is the noise level and E takes values in
H −s (N ) with some s > d/2.
Then there is C > 0, independent of δ and u† , such that
(2.18)

Eu† kUδ† (ω) − u† k2L2 ≤ C(1 + ku† k2H τ )δ

2(τ −3(t0 −t))
t0 +r

.

Note that the assumptions on CU in Theorem 2 imply that CU is a covariance operator of a random variable taking values in H τ (N ), see Remark
3 below.
In the elliptic case t = t0 we can write (2.18) as
(2.19)

2τ

Eu† kUδ† (ω) − u† k2L2 ≤ C(1 + ku† k2H τ )δ s+τ +t .

Since s = d2 +  the above convergence rate (2.19) agrees, up to  > 0 arbitrarily small, with the minimax convergence rate, see [6]. The convergence
of confidence regions is considered in Section 5.
3. Generalised random variables
This section is largely based on the work of Lasanen [31, 32]; see also
Piiroinen [41].
For any s ∈ R, let H s (N ) be the L2 -based Sobolev space equipped with
Hilbert space inner product
Z
(3.1)
((I − ∆)s/2 φ)(x) ((I − ∆)s/2 ψ)(x) dx.
(φ, ψ)H s (N ) =
N

We also define a dual pairing between H −s (N ) and H s (N )
Z
(3.2)
hφ, ψiH −s (N )×H s (N ) =
φ(x) ψ(x) dx
N

when φ, ψ ∈ C0∞ (N ). Note
H s (N ) and L2 = L2 (N ).

that

H 0 (N )

= L2 (N ). We often denote H s =

A generalised Gaussian random variable V takes values in the space of
generalised functions, and the pairing hV, φi with any test function φ ∈ D =
C ∞ (N ) is a Gaussian random variable taking values in Rd , see [45]. The
generalised Gaussian random variables we will consider below are assumed
to take values in some Hilbert space, typically in a Sobolev space H s (N ),
where the smoothness index s ∈ R may also be negative. Now, if V takes
values in H s we say that V has the covariance operator BV : H s → H s if

(3.3)
E (V − E V, φ)H s (V − E V, ψ)H s = (BV φ, ψ)H s ,
with any φ, ψ ∈ H s , see [45]. We can also define covariance operator to be
a mapping CV : H −s → H s


(3.4) E hV − E V, φiH s ×H −s hV − E V, ψiH s ×H −s = hCV φ, ψiH s ×H −s ,
with any φ, ψ ∈ H −s , see [3]. The connection between BV and CV is
BV = CV (I − ∆)s : H s → H s .
Next we will take a closer look to the generalised white Gaussian noise and
introduce the ’white noise paradox’ by a simple example.

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

11

3.1. White noise paradox. White noise E can be considered as a measurable map E : Ω → D0 (N ) where Ω is the probability space. Then normalised white noise is a random generalised function E(y, ω) on N for which
the pairings hE, φiD0 ×D are Gaussian random variables for all test functions
φ ∈ D = C ∞ (N ), EE = 0, and


(3.5)
E hE, φiD0 ×D hE, ψiD0 ×D = hIφ, ψiD0 ×D for φ, ψ ∈ D.
We will denote this by E ∼ N (0, I). A realisation of E is the generalised
function ε = E( · , ω0 ) on N with a fixed ω0 ∈ Ω.
The probability density function of white noise E is often formally written
in the form


1
(3.6)
πE (ε) = c exp − kεk2L2 (N ) .
formally
2
However, the realisations of the white Gaussian noise are almost surely not
in L2 (N ). This brings us back to the problem in formula (2.5).
Example 3. Let E be normalised white Gaussian noise defined on the ddimensional flat torus Td = (R/(2πZ))d . Let e~` ∈ L2 (Td ), ~` = (`1 , `2 , . . . , `d )
∈ Zd be an orthonormal basis of L2 (Td ) consisting of eigenfunctions of Laplacian, numbered so that −∆e~` = |~`|2 e~`. Such functions e~`(x) can be chosen
to be normalised products of the sine and cosine functions sin(`j xj ) and
cos(`j xj ) that form the standard Fourier basis of L2 (Td ). The Fourier coefficients of E with respect to this basis are independent, normally distributed
R-valued random variables with variance one, that is, hE, e~`i ∼ N (0, 1).
Then
X
X
EkEk2L2 (Td ) =
E|hE, e~`i|2 =
1 = ∞.
~
`∈Zd

~
`∈Zd

This implies that realisations of E are in L2 (Td ) with probability zero. However, when s > d/2
X
(3.7)
EkEk2H −s (Td ) =
(1 + |~`|2 )−s E|hE, e` i|2 < ∞
~k∈Zd

and hence E takes values in H −s (Td ) almost surely (that is, with probability
one)
On the other hand [44, Theorem 2] implies that if kEk2H −s (Td ) < ∞ almost
surely then EkEk2H −s (Td ) < ∞ which yields s > d/2. This concludes that
the realisations of white noise E are almost surely in the space H −s (Td ) if
and only if s > d/2. In particular for s ≤ d/2 the function x 7→ E(x, ω) is in
H −s (Td ) only when ω ∈ Ω0 ⊂ Ω where P(Ω0 ) = 0.
3.2. The smoothness of the prior. Consider the continuum measurement model Mδ = AU + Eδ where the operator is now viewed as a smoothing map A : H τ (N ) → H τ +t (N ) for all τ ∈ R. We construct the prior
by choosing U to be a generalised Gaussian random variable taking values
in H τ (N ) and having expectation E U = 0. First we will, however, give a
definition for pseudodifferential operators.

12

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

Definition 1. Let m ∈ R. We define the symbol class S m (Rd , Rd ) to consist
of a(x, ξ) ∈ C ∞ (Rd , Rd ) such that for all multi-indices α and β and any
compact set K ⊂ Rd there is such constant Cα,β,K > 0 that
|∂ξα ∂xβ a(x, ξ)| ≤ Cα,β,K (1 + |ξ|)m−|α| ,

ξ ∈ Rd , x ∈ K.

Definition 2. Let Y : U → Rd be local coordinates of the manifold N . A
bounded linear operator A : D0 (N ) → D0 (N ) is called a pseudodifferential
operator if for any local coordinates Y : U → Rd , U ⊂ N , there is a symbol
a ∈ S m (Rd × Rd ) such that for u ∈ C0∞ (U ) we have
Z
kA (y1 , y2 )u(y2 )dVg (y2 )
Au(y1 ) =
N

C ∞ (N ×N \diag(N ))

where kA |N ×N \diag(N ) ∈
and diag(N ) = {(y, y) ∈ N ×
N | y ∈ N }. Also when Y : U → V ⊂ Rd are local C ∞ -smooth coordinates
kA (y1 , y2 ) is given on U × U by
Z
−1
−1
kA (Y (x1 ), Y (x2 )) =
ei(x1 −x2 )·ξ a(x1 , ξ)dξ
Rd

where x1 , x2 ∈ V ⊂ Rd and a = aV ∈ S m (V, Rd ). In this case we will write
A ∈ Ψm (N ),
and say that in local coordinates Y : U → V ⊂ Rd the operator A has the
symbol a(x, ξ) ∈ S m (V × Rd ).
We assume that the covariance operator CU ∈ Ψ−2r , that is CU is smoothing of order 2r, self-adjoint and elliptic. With given r ∈ R we have to choose
τ ∈ R so that BU = CU (I − ∆)τ ∈ Ψ−2(r−τ ) is a trace class operator. An
operator BU is in the trace class S1 (H τ ) if
(3.8)

T rH τ →H τ (BU ) =

∞
X

|λj | < ∞

j=1

where λj are the eigenvalues of the operator BU . Condition BU ∈ S1 (H τ )
guarantees that E(U, U )H τ < ∞.
Let νj be the eigenvalues of BU−1 ∈ Ψ2(r−τ ) . Counting the geometric multiplicity of the eigenvalues, we arrange the eigenvalues of BU−1 in ascending
order as
ν1 ≤ ν2 ≤ · · · ≤ νk ≤ · · ·
Since BU−1 is a self-adjoint elliptic operator with smooth coefficients Weyl’s
law for elliptic operators tells us that the number N (ν) = #{νj | νj ≤ ν}
of the eigenvalues of BU−1 in a closed manifold less than or equal to ν has
asymptotics
d

− 1
N (ν) ∼ cν 2(r−τ ) 1 + O(ν 2(r−τ ) )
when ν → ∞.
Hence for the eigenvalues λ of the operator BU ∈ Ψ−2(r−τ )
λj ∼ c0 j −

2(r−τ )
d

(1 + o(1))

when j → ∞.

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

13

To satisfy condition (3.8) we require
∞
X
j=1

|λj | ≤ C

∞
X

j−

2(r−τ )
d

<∞

j=1

which gives us the condition τ < −(d/2 − r). From here on we will assume
that τ = r − s < r − d/2.
Remark 3. Any elliptic operator CU ∈ Ψ−2r that defines a non-negative
symmetric operator CU : D(N ) = C ∞ (N ) → D0 (N ) has the property that
BU = CU (I − ∆)τ is in S1 (H τ (N )). By [3], these yield that CU is a covariance operator of a random variable taking values in H τ (N ).
The operator BU corresponds formally to the smoothness prior




1 −1/2 2
1 −1
πpr (u) = c exp − (BU u, u)H τ = c exp − kCU ukL2
formally
2
2
Notice that the realisations of U are almost surely not r times differentiable.
In a case r < d/2 the realisations of U are almost surely not even in L2 let
alone differentiable. This is why we need to consider U as taking values in
some space H τ (N ) with possible negative smoothness index τ .
4. Proof of the main Theorems
Before we move to prove theorem 1 we will give a short introduction to
hypoelliptic pseudodifferential operators.
Definition 3. Let t, t0 ∈ R. We define symbol class HS −t,−t0 (Rd , Rd ) to
consist of a(x, ξ) ∈ C ∞ (Rd , Rd ) for which
(1) For an arbitrary compact set K ⊂ Rd we can find such positive
constants R, c1 and c2 that
c1 (1 + |ξ|)−t0 ≤ |a(x, ξ)| ≤ c2 (1 + |ξ|)−t ,

|ξ| ≥ R, x ∈ K.

(2) For any compact set K ⊂ Rd there exist constants R and Cα,β,K
such that for all multi-indices α and β
|∂ξα ∂xβ a(x, ξ)| ≤ Cα,β,K |a(x, ξ)|(1 + |ξ|)−|α| ,

|ξ| ≥ R, x ∈ K.

We will denote by HΨ−t,−t0 (N ) the class of ΨDO with local symbol a(x, ξ) ∈
HS −t,−t0 (V × Rd ), see Definition 2.
We denote H r (N ) = H r and L2 (N ) = L2 where N is a closed manifold
and dim N = d.
The proof of Theorem 1 is rather long and technical so we will start by
going through the main steps of it in a nutshell. The approximated solution
we are studying is of the form

−1/2
Tδ (mδ ) := arg min kAuk2L2 (N ) − 2hmδ , Aui + δ 2 kCU uk2L2 (N ) .
u∈H r (N )

As mentioned before the solution to this is
(4.1)

Tδ (mδ ) = (A∗ A + δ 2 CU−1 )−1 A∗ mδ .

14

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

We can rewrite the above as
(4.2)

Tδ (mδ ) = Zδ−1 A∗ Au + Zδ−1 A∗ (εδ)
= u − δ 2 Zδ−1 CU−1 u + Zδ−1 A∗ (εδ)

where Zδ = A∗ A + δ 2 CU−1 .
To study the convergence of the last term on the right hand side of (4.2)
we would like to write it in the form
(4.3)

Zδ−1 A∗ (εδ) = δ −1 Fδ−1 (A∗ A)−1 A∗ ε

where Fδ = (A∗ A)−1 CU−1 + δ −2 . In order to show that (4.3) is well-defined
we first need to prove that A∗ A and Fδ are invertible.
Lastly we study the converge of δ 2 Zδ−1 CU−1 u and δ −1 Fδ−1 (A∗ A)−1 A∗ ε to
zero in appropriate Sobolev spaces and show that the latter term is always
dominating.
We will start by showing that A∗ A ∈ HΨ−2t,−2t0 is invertible and (A∗ A)−1 ∈
HΨ2t0 ,2t . Define A∗ : L2 (N ) → L2 (N ) as the adjoint of an operator
A : L2 → L2 . We assumed in Theorem 1 that A : L2 → L2 is one-toone. Since A∗ A : H r → H r+2t , r ∈ R, is hypoelliptic [49, Propositions 5.2
and 5.3] A∗ Au = 0 ∈ C ∞ implies that u ∈ C ∞ and hence Au ∈ L2 . Now we
see that if A∗ Au = 0 then
0 = (A∗ Au, u)L2 = (Au, Au)L2 = kAuk2L2
which implies Au = 0 and furthermore u = 0. Thus the operator A∗ A :
H r (N ) → H r+2t (N ) is one-to-one.
To study the mapping A∗ A ∈ HΨ−2t,−2t0 from some Sobolev space H r ,
r ∈ R, we define Lr = A∗ A : H r → H r+2t . The adjoint of Lr is denoted by
L0r = (A∗ A)0 : H −(r+2t) → H −r . Let φ, ψ ∈ C ∞ . Then
hA∗ Aφ, ψiH r+2t ×H −(r+2t) = hφ, A∗ AψiH r ×H −r ,
that is, L0r = L−(r+2t) and hence the adjoint is one-to-one. Now we can
conclude that Lr (H r ) ⊂ H r+2t is a dense subset. Next we will prove two
lemmas that show that the operator A∗ A : D0 → D0 is also surjective.
Since A∗ A ∈ HΨ−2t,−2t0 is a hypoelliptic pseudodifferential operator it
has a parametrix B1 ∈ HΨ2t0 ,2t [49, Theorem 5.1]. Hence for any r0 > 0 we
get norm estimates

kA∗ AukH r+2t ≤ C1 kukH r
(4.4)
kukH r ≤ C2 kA∗ AukH r+2t0 + C3 kukH r−r0
for all u ∈ C ∞ . Next we will show that C3 is zero.
Lemma 1. Let A∗ A ∈ HΨ−2t,−2t0 be an injective hypoelliptic pseudodifferential operator. Then we have the following estimates
C1 kA∗ AukH r+2t ≤ kukH r ≤ C2 kA∗ AukH r+2t0
Proof. We get the first inequality since A∗ A is continuous linear operator.
If the second estimate in (4.4) is not valid with C3 = 0 for any C2 > 0
then we can choose a sequence uj such that kuj kH r−r0 = 1 and kuj kH r ≥

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

15

jkA∗ Auj kH r+2t0 . When j > 2C2 then
kuj kH r ≤ C2 kA∗ Auj kH r+2t0 + C3 kuj kH r−r0
1
≤ kuj kH r + C3 kuj kH r−r0 .
2
This gives us
kuj kH r ≤ 2C3 kuj kH r−r0 = 2C3 .
Since r0 > 0, the embedding H r ,→ H r−r0 is compact. Now there exists a
subsequence uj` and such a w ∈ H r−r0 that lim`→∞ uj` = w in H r−r0 . We
assumed that kuj` kr−r0 = 1 which implies kwkr−r0 = 1. On the other hand
j` kA∗ Auj` kr+2t0 ≤ kuj` kr ≤ 2C3 ,
that is, lim`→∞ kA∗ Auj` kr+2t0 = 0 and because −r0 + 2t ≤ 2t0
lim kA∗ Auj` kr−r0 +2t = 0.

`→∞

Since uj` → w in H r−r0 we also have A∗ Auj` → A∗ Aw in H r−r0 +2t . Combining the above results we see that kA∗ AwkH r−r0 +2t = 0. Operator A∗ A is
one to one and hence w = 0. This is a contradiction since kwkH r−r0 = 1. 
Lemma 2. Let A∗ A : H r → H r+2t be an injective hypoelliptic operator.
Then the image of H r in the map A∗ A satisfies
H r+2t0 ⊂ A∗ A(H r ) ⊂ H r+2t .
Proof. The second inclusion is a direct consequence of the mapping properties of A∗ A. Let f ∈ H r+2t0 . Since C ∞ ⊂ H r+2t0 is a dense subset we
can find such a sequence fj ∈ C ∞ that limj→∞ fj = f in H r+2t0 . Since
A∗ A(H r ) ⊂ H r+2t is dense we can also choose a sequence hj,` = A∗ Agj,` ∈
A∗ A(H r ) such that lim`→∞ hj,` = fj in H r+2t0 . Denote gj = gj,`j ∈ H r for
which limj→∞ A∗ Agj = f in H r+2t0 . Using Lemma 1 we see
lim kgj − gk kH r ≤ C2 lim kA∗ Agj − A∗ Agk kH r+2t0 = 0.

j,k→∞

j,k→∞

Hr

Hence also gj ∈
is a Cauchy sequence. Thus there exists such g ∈ H r
that limj→∞ gj = g in H r . On the other hand,
lim kA∗ Agj − A∗ AgkH r+2t ≤ C1 lim kgj − gkH r = 0.

j→∞

k→∞

Combining the above we get

A∗ Ag

= f.



Using Lemma 2 we see that A∗ A(D0 ) = D0 , that is the operator A∗ A
is also onto. Now we can conclude that there exists an inverse operator
(A∗ A)−1 : D0 → D0 . It remains to show that the inverse operator is a
hypoelliptic pseudodifferential operator.
Lemma 3. A self-adjoint, smoothing, one-to-one hypoelliptic operator A∗ A ∈
HΨ−2t,−2t0 has an inverse operator (A∗ A)−1 ∈ HΨ2t0 ,2t .
Proof. Denote B = (A∗ A)−1 : D0 → D0 . For an operator A∗ A : H r → H r+2t
we define B0 ⊂ B with domain
D(B0 ) = {f ∈ H r+2t | Bf ∈ H r } = A∗ A(H r ).

16

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

Using the hypoellipticity of A∗ A we see that A∗ Au = f ∈ C ∞ implies
u ∈ C ∞ . This gives us B0 : C ∞ → C ∞ . Since C ∞ is a Frechét space and
A∗ A is continuous and linear A∗ A : C ∞ → C ∞ is an open mapping [46,
Theorem 2.11]. Hence the operator B0 : C ∞ → C ∞ is continuous.
Since A∗ A is hypoelliptic it has a parametrix B1 ∈ HΨ2t0 ,2t [49, Theorem
5.1]

B1 (A∗ A) = I + K1 , K1 ∈ Ψ−∞
(A∗ A)B1 = I + K2 , K2 ∈ Ψ−∞
and we can write
B0 = B0 ((A∗ A)B1 − K2 ) = B1 − B0 K2 .
The operator B0 K2 : D0 → C ∞ is continuous and thus we have shown that
B0 = B1

mod Ψ−∞ .

That is B0 ∈ HΨ2t0 ,2t .



Next we will examine ΨDOs that depend on spectral variable λ = δ −2 .
For the general theory see [49].
Definition 4. The symbol class Spm (Rd × Rd , R+ ) consist of the functions
a(x, ξ, λ) such that
(1) a(x, ξ, λ0 ) ∈ C ∞ (Rd × Rd ) for every fixed λ0 ≥ 0 and
(2) for arbitrary multi-indices α and β and any compact set K ⊂ Rd
there exists a constant Cα,β,K such that
|∂ξα ∂xβ a(x, ξ, λ)| ≤ Cα,β,K (1 + |ξ| + |λ|1/p )m−|α|
for x ∈ K, ξ ∈ Rd and λ ≥ 0.
We denote by Ψm
p (N, R+ ) the class of pseudodifferential operators Aλ for
which the local symbol a(x, ξ, λ) ∈ Spm (V × Rd , R+ ), see Definition 2.
Definition 5. If there are constants C1 , C2 , R > 0 such that the symbol
a(x, ξ, λ) ∈ Spm (Rd × Rd , R+ ) satisfies
C1 (|ξ| + |λ|1/p )m0 ≤ |a(x, ξ, λ)| ≤ C2 (|ξ| + |λ|1/p )m ,
for |ξ| + |λ| ≥ R, we say that a is hypoelliptic with parameter λ and denote
0
a(x, ξ, λ) ∈ HSpm,m0 (Rd × Rd , R+ ). We will denote by HΨm,m
(N, R+ ) the
p
class of ΨDOs depending on the parameter λ whose local symbol belongs to
HSpm,m0 (V × Rd , R+ ), see Definition 2.
Next we will prove that
Fλ = (A∗ A)−1 CU−1 + λ
is invertible. Operator Fλ ∈ HΨ2(t0 +r),2(t+r) (N ) is hypoelliptic since (A∗ A)−1 CU−1 ∈
HΨ2(t0 +r),2(t+r) (N ) is hypoelliptic and λ I ∈ Ψ0 (N ). Denote Q = (A∗ A)−1 CU−1
and its symbol q(x, ξ) ∈ HS 2(t0 +r),2(t+r) (N ). Then for the symbol σ(Fλ )(x, ξ) =
q(x, ξ) + λ of the operator Fλ
|∂ξα ∂xβ (q(x, ξ) + λ)| ≤ Cα,β (1 + |ξ| + |λ|1/(2(t0 +r)) )2(t0 +r)−|α| .

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

17

By [49, Theorem 9.2.] there exist R > 0 such that for |λ| ∈ [R, ∞) the
2(t +r),2(t+r)
operator Fλ ∈ HΨ2(t00 +r)
(N, R+ ) is invertible with
−2(t+r),−2(t0 +r)

Fλ−1 ∈ HΨ2(t0 +r)

(N, [R, ∞)).

We have now shown that the operator Zδ can be written

−1
−1
∗
−1 −1
(4.5)
Zδ = λ (A A) CU + λI
(A∗ A)−1
where λ = δ −2 . Hence we can rewrite (4.2)
√
(4.6)
Tλ (m) = u − λ−1 Zδ−1 CU−1 u + λFλ−1 (A∗ A)−1 A∗ ε.
Now we will proceed to study the convergence of the second and third
term on the right hand side of (4.6). For the third term of (4.6) we have
(A∗ A)−1 A∗ : H −s → H k , k = −s + t − 2t0 and Fλ−1 : H k → H k+2(t+r) .
Hence when ζ ≤ k + 2(t + r) we have
kFλ−1 (A∗ A)−1 A∗ εkH ζ

≤ kFλ−1 kk,ζ k(A∗ A)−1 A∗ εkH k

where kFλ−1 kk,ζ is the norm of Fλ−1 : H k (N ) → H ζ (N ) and k, ζ ∈ R. Next
we want to study what happens to the norm when λ → ∞.
We have the following norm estimates for Fλ−1 ∈ Ψm
p (N, R+ ) when ` ≥ m
and λ large enough [49, Theorem 9.1.]
(4.7)

kFλ−1 kk,k−` ≤ Ck,` (1 + |λ|1/p )m ,

if ` ≥ 0

(4.8)

kFλ−1 kk,k−` ≤ Ck,` (1 + |λ|1/p )−(`−m) ,

if ` ≤ 0.

In our case Fλ−1 ∈ Ψm
p where m = −2(t + r) and p = 2(t0 + r). We will write
−1
−1
kFλ kk,ζ = kFλ kk,k−` where ` = k − ζ ≥ m.
First we study the case when ` ≥ 0 that is ζ ≤ k. Inequality (4.7) gives
us the norm estimate
kFλ−1 kk,k−` ≤ C(1 + |λ|1/p )m .

√
Because we want λkFλ−1 (A∗ A)−1 A∗ εkH ζ to converge when λ → ∞ we
have to require that
m
−2(t + r)
1
=
<− .
p
2(t0 + r)
2
This is true when t0 < 2t + r.
When ` ≤ 0 we have k ≤ ζ ≤ k + 2(t + r) and can use (4.8)
kFλ−1 kk,k−` ≤ C(1 + |λ|1/p )−(`−m) .
For convergence we need
m−`
−2(t + r) − k + ζ
1
=
<−
p
2(t0 + r)
2
that is k ≤ ζ < k + 2t + r − t0 = r − s − 3(t0 − t) which can be true only if
t0 < 2t + r.
Next we will prove the convergence of the term δ 2 Zδ−1 CU−1 u in H ζ . Since
we got above that ζ < τ − 3(t0 − t) we can write ζ = τ − θ where θ ≥
3(t0 − t) ≥ 0. We need to find such η ≥ 0 and γ ≥ 0 that γ + η = 1 and

18

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

t0 γ − rη + r − θ/2 = 0. Define γ = θ/2(t0 + r) and η = 1 − θ/2(t0 + r). Now
η ≥ 0 only if θ ≤ 2(t0 + r). Hence we will choose θ = min{τ − ζ, 2(t0 + r)}.
Since Zδ = A∗ A + δ 2 CU−1 where A∗ A ≥ c1 (I − ∆)−t0 and c2 (I − ∆)r ≤
−1
CU ≤ c3 (I − ∆)r we get
ζ

kδ 2 Zδ−1 CU−1 ukH ζ ≤ δ 2 k(A∗ A)−γ (δ 2 CU−1 )−η c3 (I − ∆)r+ 2 ukL2
θ

τ

≤ δ 2 k(c1 (I − ∆)−t0 )−γ (c2 δ 2 (I − ∆)r )−η (I − ∆)r− 2 + 2 ukL2

(4.9)

θ

= Cδ t0 +r kukH τ
where θ = min{τ − ζ, 2(t0 + r)}.
Adding the above results together we can prove Theorem 1.
Proof of Theorem 1. To get the speed of convergence we use the fact that
U and E are independent. Similarly to (4.6) we get
(4.10)

1
∗
−1 ∗
Uδ (ω) − U (ω) = −δ 2 Zδ−1 CU−1 U (ω) + Fδ−1
A E(ω),
−2 (A A)
δ

where by (4.9),
θ

kδ 2 Zδ−1 CU−1 U kH ζ ≤ Cδ t0 +r kU kH τ .

(4.11)

For the second part on the right hand side of (4.10) we can write
1
1
∗
−1 ∗
∗
−1 ∗ p
(4.12) Ek Fδ−1
A E(ω)kpH ζ ≤ k Fδ−1
A kH −s →H ζ EkE(ω)kpH −s ,
−2 (A A)
−2 (A A)
δ
δ
with p ∈ {1, 2} and θ = min{τ − ζ, 2(t0 + r)}.
When ζ ≤ t − s − 2t0 we get
EkUδ (ω) − U (ω)kH ζ

≤ δ 2 EkZδ−1 CU−1 U (ω)kH ζ + δ −1 EkFδ−1 (A∗ A)−1 A∗ E(ω)kH ζ
θ

≤ C1 δ t0 +r EkU (ω)kH τ + C2 δ
n θ
2t+r−t0 o
≤ C max δ t0 +r , δ t0 +r

−1+

2(t+r)
t0 +r

EkEkH −s

where θ = min{τ − ζ, 2(t0 + r)}. Next we will study which of the terms is
dominating. The noise term δ −1 Fδ−1 (A∗ A)−1 A∗ E(ω) is dominating if
2t + r − t0 ≤ θ.
Assume first that θ = τ − ζ. Then
θ = τ − ζ ≥ 2t0 + r − t ≥ 2t + r − t0 .
If θ = 2(t0 + r) we get
θ = 2(t0 + r) ≥ 2t0 + t0 − 2t + r ≥ 2t + r − t0
since t ≤ t0 < 2t + r. Hence the noise term is dominating in both cases and
we have proven
EkUδ (ω) − U (ω)kH ζ ≤ Cδ

2t−t0 +r
t0 +r

.

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

19

If t − s − 2t0 ≤ ζ < τ − 3(t0 − t) we get
EkUδ (ω) − U (ω)kH ζ

≤ δ 2 EkZδ−1 CU−1 U (ω)kH ζ + δ −1 EkFδ−1 (A∗ A)−1 A∗ E(ω)kH ζ
τ −ζ

≤ C1 δ t0 +r EkU (ω)kH τ + C2 δ
n τ −ζ 2t+r−t0 +` o
≤ C max δ t0 +r , δ t0 +r

−1+

2(t+r)+`
t0 +r

EkE(ω)kH −s

Above ` = t−s−2t0 −ζ. Note that when ζ ≥ t−s−2t0 then τ −ζ ≤ 2(r+t0 ).
The noise term is dominating if
2t + r − t0 + ` ≤ τ − ζ.
This is always true since τ = r − s and
2t + r − t0 + ` = r − s − ζ − 3(t0 − t) ≤ τ − ζ.
Hence we can conclude
EkUδ (ω) − U (ω)kH ζ ≤ Cδ

2t+r−t0 +`
t0 +r

.


Proof of Theorem 2. Similarly to (4.10) we get
1
∗
−1 ∗
Uδ† (ω) − u† = −δ 2 Zδ−1 CU−1 u† + Fδ−1
(4.13)
A E(ω),
−2 (A A)
δ
where the first term on the right side satisfies, by (4.9),
θ

(4.14)

kδ 2 Zδ−1 CU−1 u† kH ζ ≤ Cδ t0 +r ku† kH τ

with θ = min{τ − ζ, 2(t0 + r)}. The expectation of the second term in the
right side is estimated in (4.12). Analysing the obtained terms as in the
proof of Theorem 1, we obtain
EkUδ† (ω) − u† k2H ζ ≤ C(1 + ku† k2H τ )δ

2(τ −3(t0 −t))
t0 +r

.


5. Posterior distribution and confidence regions
One advantage Bayesian inversion offers over deterministic regularization
is uncertainty quantification. Since the solution to the Bayesian inverse
problem is the posterior distribution of the unknown we can study its credible sets and their contraction in some Sobolev space H ζ when δ → 0. A
Bayesian credible set is a region in the posterior distribution that contains a
large fraction of the posterior mass, for instance, 95%. We are dealing with
Gaussian distributions so we define the credible sets to be central regions.
This means these sets are defined as central balls with uδ as a centre.
The above mentioned credible sets are often used to visualise the remaining Bayesian uncertainty in the estimate. Frequentists use another kind of
uncertainty quantification called confidence region. A confidence region is
a range of values that frequently includes the unknown of interest if the experiment is repeated. We can define confidence regions as central balls with
u†δ as the centre. Here u†δ is the frequentist approximated solution generated

20

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

by a true solution u† . How frequently the ball around the approximated
solution, with different realisation of the noise, contains the true solution is
determined by the confidence level. See for example [15, 56].
In the finite-dimensional parametric case and under mild conditions on
the prior Bernstein–von Mises theorem provides that the credible sets of
smooth models are asymptotically equivalent with the frequentist confidence
regions based on the maximum likelihood estimator, see [56]. In infinitedimensional case there is no corresponding theorem and Bayesian credible
sets are not automatically frequentist confidence sets. This means that
if we assume that the data is generated by a ‘true parameter’, it is not
automatically true that credible sets contain that truth with probability
at least the credible level. However the correspondence of Bayesian and
frequentist uncertainty has been studied in many recent papers see e.g. [4,
5, 27, 35, 43, 52]. These results are important since they show that some
credible sets can give a good idea of the uncertainty of the estimate in
the classical sense. In this section we show that the posterior distribution
converges and we give some convergence rates. We also prove that in the
elliptic case the frequentist posterior contractions rates agrees, up to ε > 0
arbitrarily small, with the minimax convergence rate. We do not address
the question about the frequentist coverage of the credible sets.
We will start by studying the convergence of the posterior covariance
Cδ which, with the convergence of the posterior mean Uδ , guarantees the
convergence of the posterior distribution.
When U ∼ N (0, CU ), E ∼ N (0, I) and
(5.1)

Mδ = AU + δE

the conditional probability distribution of U with respect to the measurement Mδ is a Gaussian measure with mean Uδ and covariance [36, 37]
(5.2)

Cδ = CU − CU A∗ (ACU A∗ + δ 2 I)−1 ACU .

If A∗ : D0 (N ) → D0 (N ) is invertible we can rewrite the above
−1
Cδ = δ 2 A∗ A + δ 2 CU−1
−1

(5.3)
= (A∗ A)−1 CU−1 + δ −2 I
(A∗ A)−1 .
Note that the covariance operator is deterministic and thus independent of
Mδ .
We define Fλ = (A∗ A)−1 CU−1 + λI, where λ = δ −2 , as in section 4. Then
Fλ−1 ∈ Ψm
p (N, [R, ∞))
where m = −2(t + r) and p = 2(t0 + r). Using the norm estimate (4.8) we
get
−1
kδ 2 Zδ−1 k−τ,τ = k (A∗ A)−1 CU−1 + δ −2 I
k−τ −2t,τ k(A∗ A)−1 k−τ,−τ −2t
≤ ckFλ−1 kk,k−`
1

≤ c(1 + λ p )−(`−m) .

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

21

Above k = −τ − 2t and ` = −2(τ + t) < 0. Since τ = r − s we can write
`−m
−τ + r
=
p
t0 + r
s
=
t0 + r
and hence we get the following convergence rate for the posterior covariance
2s

kCδ k−τ,τ

(5.4)

≤ cδ t0 +r
d

≤ cδ t0 +r .
We see that the more smoothing the forward operator A is the worse
convergence we get. Note that r and s do not only affect the convergence
speed but also the spaces between which the norm is taken.
Remark 4. Observe that the random variable U takes values in H τ and the
estimate (5.4) concerns the mapping properties of the posterior covariance
operator Cδ from the dual space H −τ = (H τ )0 to the space H τ . For strictly
positive δ > 0 the MAP estimator Uδ belongs to the space H r , r = τ +s ≥ τ ,
but as δ → 0, the MAP estimators Uδ converge in a less regular space H ζ ,
ζ < τ − 3(t0 − t) ≤ τ , see (2.10).
5.1. Contraction of the posterior distribution. Next we consider the
inverse problem using the frequentist setting described in Subsection 2.1.2
with the additional assumption that τ > 0. We assume below that CU
satisfies the assumptions in Theorem 2 that in particular imply that CU is
the covariance operator of a random variable U taking values in H τ (N ), see
Remark 3. We recall that we consider a fixed ‘true’ solution u† ∈ H τ (N )
and the noise model Mδ† (ω) = Au† + δE(ω) as in (2.14). Also, note that the
MAP-estimate is then Uδ† = Tδ (Mδ† ).
In the frequentist case one is often interested in the the limiting behaviour
of the posterior measure PM † when δ → 0. Here, PM † is a random measure in
δ

δ

H τ (N ), depending on δ and the MAP estimator Uδ† = Uδ† (ω) = Tδ (Mδ† (ω))
(that further depends on the deterministic variable u† and the realisation
E(ω) of the random noise). Let Wδ be a Gaussian random variable, taking
values in H τ (N ), that is independent of the noise E, has zero mean and
the covariance operator Cδ , see (3.4). For a measurable set B ⊂ H τ (N ) we
define
n
o
(5.5)
PM † (B) = µδ wδ ∈ H τ (N ) | wδ = b − Uδ† (ω) with b ∈ B
δ

where µδ = N (0, Cδ ). Roughly speaking, PM † is a Gaussian measure in
δ

H τ (N ) with the mean Uδ† and the covariance operator Cδ .
Recall that we consider the probability space (Ω, Σ, P) and denote by χS
the indicator function of S. Let S ∈ Σ. We use the notations
(5.6) EM † (F (ω, u† )) := E(F (ω, u† )|F),
δ

PM † (S) := P(S|F) = E(χS |F),
δ

for the conditional expectation and conditional probability. Above F ⊂ Σ is
the σ-algebra generated by random variable Mδ† (ω) or equivalently, the noise
E(ω). Roughly speaking, in the notation Eu† F the subindex u† reminds that

22

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

u† is a fixed parameter and the expectation is taken only with respect the
noise. The notation PM † (B) indicates that the measure of B is computed
δ

using the posterior probability measure which mean Uδ† = Tδ (Mδ† ) depends
on the measurement Mδ† . Since the random variable Wδ has distribution µδ ,
we have by (5.5)
Z
(5.7) PM † (B) =
χB (wδ + Uδ† ) dµδ (wδ ) = PM † ({Wδ + Uδ† ∈ B}).
δ

δ

H τ (N )

Following the approach in [11, 14, 27, 57] we next show that the posterior
measure contracts to a Dirac measure centred on the fixed true solution u† .
Theorem 3. Let r > s > d/2 and N be a d-dimensional closed manifold.
Let u† ∈ H τ (N ) where τ = r − s > 0. Assume that CU , the covariance
operator of the Gaussian prior, is a self-adjoint, injective and elliptic pseudodifferential operator of order −2r. Let E(y, ω) be white Gaussian noise on
N . Consider the measurement
Mδ† (y, ω) = A(u† (·)) + δE(y, ω),

ω ∈ Ω,

where A ∈ HΨ−t,−t0 , t > max{0, −τ + s} and t ≤ t0 ≤ t + τ /3, is a
hypoelliptic pseudodifferential operator on the manifold N and A : L2 (N ) →
L2 (N ) is injective. We assume also that A∗ : D0 (N ) → D0 (N ) is invertible.
Above δ ∈ R+ is the noise level and E takes values in H −s (N ) with some
s > d/2. Let Uδ† be the MAP estimated given by (2.15).
0 −t))
Let κ < κ0 = 2(τ −3(t
, c0 > 0, and R > 0. Then there is c1 > 0 such
t0 +r
that
(5.8)
o
n
sup Eu† PM † u ∈ H τ (N ) ku − u† kL2 (N ) ≥ c0 δ κ ≤ c1 δ 2(κ0 −κ) → 0,
ku† kH τ ≤R

δ

as δ → 0.
Proof. Let u† ∈ H τ (N ) and Wδ be the Gaussian variable defined above.
Using the Markov inequality and (5.7), we get
o
o
n
n
PM † u ∈ H τ (N ) ku − u† kL2 ≥ c0 δ κ = PM † kWδ + Uδ† − u† kL2 (N ) ≥ c0 δ κ
δ

δ

≤

(5.9)


1
†
† 2
E
† kWδ + U − u kL2 (N ) .
δ
(c0 δ κ )2 Mδ

Since Wδ and Uδ† are independent and Wδ has the covariance operator Cδ ,
we obtain using notations (5.6)
Eu† EM † (kWδ + Uδ† − u† k2L2 (N ) ) = Eu† EM † kWδ k2L2 + Eu† EM † kUδ† − u† k2L2
δ

(5.10)

δ

δ

= T rL2 (N )→L2 (N ) (Cδ ) + Eu† kUδ† − u† k2L2 (N ) .

We have shown in Theorem 2 that the second term on the right side of
(5.10) can be estimated by c2 (1 + R2 )δ 2κ0 with some c2 > 0. Hence it is
enough to show that
T rL2 →L2 (Cδ ) ≤ cδ 2κ0

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

23

2(τ −3(t0 −t))
.
t0 +r

We can estimate the trace by writing

T rL2 →L2 (Cδ ) = T rL2 →L2 (I − ∆)−s (I − ∆)s Cδ

≤ T rL2 →L2 (I − ∆)−s k(I − ∆)s Cδ kL2 →L2

with κ0 =

Above (I − ∆)−s is trace class operator in L2 since
∞
X
(j 2/d )−s < ∞ when s > d/2.
j=1

As before we get
k(I − ∆)s Cδ k0,0 = k(I − ∆)s k2s,0 kFλ−1 k−2t,2s k(A∗ A)−1 k0,−2t
≤ ckFλ−1 k−2t,−2t−`
1

≤ c(1 + λ p )−(`−m) .
Above ` = −2(s + t) < 0. We can write
`−m
τ
=
p
t0 + r
and hence
2τ

T rL2 →L2 (Cδ ) ≤ cδ t0 +r ≤ cδ 2κ0 .

Note that in the elliptic case t0 = t we get contraction
o
n
Eu† PM † u ∈ H τ (N ) ku − u† kL2 (N ) ≥ c0 δ κ → 0
δ

2τ
d
when δ → 0 for all c0 > 0 and κ < s+τ
+t . Since s = 2 +  the above
convergence rate agrees, up to  > 0 arbitrarily small, with the minimax
convergence rate.

Remark 5. Above we have assumed that u† is in H τ . This correspond
to the fact that the random variable U , having the covariance operator
CU ∈ Ψ−2r , takes values in H τ . The L2 (N ) norm in the contraction formula
(5.8) can be considered as a loss function on H τ (N ). Note that the loss
function d(v1 , v2 ) = kv1 − v2 kL2 defines a distance function in the H τ (N ),
but the obtained metric space is not complete. When the direct map A is
the identity map, similar estimates with different loss functions have been
studied in a general setting in [18]. However, from the point of view of
inverse problems [18] corresponds to the case when the direct operator and
the covariance operator of the prior commute. This differs from the problem
analysed in our paper, where covariance operator CU and the operator A
may not commute, and are of quite different type in the sense that CU is an
elliptic operator but A is hypoelliptic operator. The phenomenon that the
solution u† is assumed to be in a smoother space, in our case in H τ , and the
convergence of the posterior distribution is analysed using a loss function
given by a less strict norm, in our case L2 -norm, appears in many frequentist

24

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

studies, see e.g. Theorems 2.2 and 2.3 and Remark 3.6 in [1]. Conditions
similar to the smoothness requirement u† ∈ H τ are also encountered in
classical regularisation theory [9] where this type of conditions are called
source conditions.
5.2. Convergence of the posterior distribution in Bayesian settings.
Next we will proceed to study the contraction of the posterior distribution
using Bayesian techniques, the measurement model Mδ = AU + δE and the
MAP estimator Uδ = Tδ (Mδ ). Let us write
Vδ = Uδ + Wδ
where Wδ ∼ N (0, Cδ ) is a Gaussian variable having the covariance operator
Cδ given in (5.3) and the zero mean. Random variables Wδ and Uδ are
assumed to be independent. Let Mδ = σ(Mδ ) be the σ-algebra generated
by the random variable Mδ . Then the distribution of the random variable Vδ
is the same as the posterior distribution of U with respect to the σ-algebra
Mδ .
Let νδ be the posterior distribution of U with respect to the σ-algebra
Mδ . Equivalently νδ is the distribution of Vδ in the Sobolev space H ζ1 (N )
where ζ1 ≤ τ . Let µδ be the distribution of the random variable Wδ which
is independent of Uδ = Tδ (Mδ ). Then the conditional expectation of the
indicator function χBζ1 (Uδ (ω),R(δ)) (U (ω)) with respect to Mδ is


E χBζ1 (Uδ ,R(δ)) | Mδ (ω) = P {U ∈ Bζ1 (Uδ , R(δ)), R(δ))} | Mδ (ω)

(5.11)
= P {Wδ ∈ Bζ1 (0, R(δ))} (ω)

= µδ Bζ1 (0, R(δ)) (ω).
Above Bζ1 (0, R(δ)) denotes a ball in H ζ1 of radius R(δ).
Let U be the σ-algebra generated by U . By [8, Theorem 10.2.2] there
are regular conditional probabilities P(K | Mδ )(ω) for all K ∈ B(H ζ1 ) and
ω ∈ Ω such that
P(K | Mδ )(ω) = E(χK (U ) | Mδ )(ω)

a.s.

Moreover, by applying [8, Theorem 10.2.1] to the joint distribution of (U, Mδ )
we see that there are such functions

(m, K) 7→ Pm (K) =: P {U ∈ K} | Mδ = m ,
defined for m ∈ H −s and K ∈ B(H ζ1 ), that
(5.12)

PMδ (ω) (K) = P(K | Mδ )(ω)

a.s.

Using (5.11) and (5.12) we see that



P U ∈ Bζ1 (Tδ (m), R(δ)) | Mδ = m = µδ Bζ1 (0, R(δ)) .
Note that the right hand side is in fact independent of m and depends only
on δ. Next we will give a theorem for the credible sets




P Vδ ∈ Bζ1 (Uδ , R(δ)) = P U ∈ Bζ1 (Tδ (mδ ), R(δ)) | Mδ = mδ .

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

25

Theorem 4. Let U , E, Mδ and A be defined as in Theorem 1 and assume
that A∗ : D0 (N ) → D0 (N ) is invertible.
Let Mδ = σ(Mδ ) and U = σ(U ) be the σ-algebras generated by the random variable Mδ and U respectively. Then the posterior distribution of the
random variable U with respect to the σ-algebra Mδ can be given in terms
of function

(m, K) 7→ P {U ∈ K} | Mδ = m ,
where m ∈ H −s (N ) and K ∈ B(H ζ1 ), cf. (5.12).
Take ζ1 < τ + t − t0 and α < γ/2. Then if R(δ) = C1 δ α we have the
following contraction:


P U ∈ Bζ1 (Tδ (mδ ), R(δ)) | Mδ = mδ ≥ 1 − Cδ γ−2α → 1
when δ → 0. The speed of contraction depends on ζ1 :
(i) If ζ1 ≤ −s − t0 then γ = 2(t + r)/(t0 + r).
(ii) If −s − t0 ≤ ζ1 < τ + t − t0 then γ = 2(τ + t − t0 − ζ1 )/(t0 + r).
Proof. We use below R = R(δ) = C1 δ α with some α > 0 and denote
pδ = 1 − µδ (Bζ1 (0, R(δ))).
To study what happens to pδ we first notice that
Z

2
kwk2H ζ1 dµδ (w)
E kWδ kH ζ1 =
ζ
1
ZH
≥
kwk2H ζ1 dµδ (w)
2
2
kwk ζ >(R(δ))
(5.13)
H 1
Z
≥ R(δ)2

kwk2

dµδ (w)

>(R(δ))2
H ζ1

= (R(δ))2 pδ .
Next we will prove that

E kWδ k2H ζ1 = T rH ζ1 →H ζ1 BWδ ≤ Cδ γ
with some γ that depends on ζ1 . Above we use the definition BWδ = CWδ (I −
∆)ζ1 where CWδ : H −ζ1 → H ζ1 and


hCWδ φ, ψiH ζ1 ×H −ζ1 = E hWδ , φiH ζ1 ×H −ζ1 hWδ , ψiH ζ1 ×H −ζ1 .
As noted before when A∗ is invertible we can write
CWδ = Fλ−1 (A∗ A)−1

(5.14)

where Fλ = (A∗ A)−1 CU−1 + λI and λ = δ −2 .
We want to estimate


E kWδ k2H ζ1 = E k(I − ∆)ζ1 /2 Wδ k2L2 .
Let us define
Wδ,ζ1 = (I − ∆)ζ1 /2 Wδ : Ω → L2 (N ).
We can write the covariance operator of Wδ,ζ1
CWδ,ζ1 = (I − ∆)ζ1 /2 CWδ (I − ∆)ζ1 /2 : L2 → L2 .

26

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

Note that in L2 we have BWδ,ζ1 = CWδ,ζ1 . Now we get for the trace
T rL2 →L2 ((I − ∆)ζ1 /2 CWδ (I − ∆)ζ1 /2 ) = T rL2 →L2 ((I − ∆)ζ1 CWδ )
= T rL2 →L2 ((I − ∆)−s (I − ∆)ζ1 +s CWδ )
≤ T rL2 →L2 ((I − ∆)−s )k(I − ∆)ζ1 +s CWδ kL2 →L2 .
Using (5.14) we get
k(I − ∆)ζ1 +s CWδ kL2 →L2

= k(I − ∆)ζ1 +s Fλ−1 (A∗ A)−1 k0,0
≤ k(I − ∆)ζ1 +s k2(ζ1 +s),0 kFλ−1 k−2t0 ,2(ζ1 +s) k(A∗ A)−1 k0,−2t0
≤ CkFλ−1 k−2t0 ,2(ζ1 +s) .

Above Fλ−1 = ((A∗ A)−1 CU−1 + λ)−1 ∈ Ψm
p where m = −2(t + r) and p =
2(t0 + r). We want to use the norm estimates (4.7) and (4.8) so we write
kFλ−1 k−2t0 ,2(ζ1 +s) = kFλ−1 k−2t0 ,−2t0 −`
where ` = −2(s + ζ1 + t0 ). To use the norm estimates we need to assume
` ≥ m, that is, ζ1 ≤ r − s + t − t0 = τ + t − t0 .
First we assume that ` ≥ 0 which is true when ζ ≤ −s − t0 . Then
1

kFλ−1 k−2t0 ,−2t0 −` ≤ C(1 + λ 2(t0 +r) )−2(t+r)
and kFλ−1 k−2t0 ,−2t0 −` → 0 when λ → ∞ with all t, t0 > 0 and r ≥ 0.
Next we assume ` ≤ 0. Then for −s − t0 ≤ ζ1 ≤ τ + t − t0 we get
1

kFλ−1 k−2t0 ,−2t0 −` ≤ C(1 + λ 2(t0 +r) )2(s+ζ1 +t0 )−2(t+r) .
Now kFλ−1 k−2t0 ,−2t0 −` → 0 when λ → ∞ if ζ1 < τ + t − t0 .
We have proven that

E kWδ k2H ζ1 ≤ Cδ γ
where γ = 2(t + r)/(t0 + r) if ζ1 ≤ −s − t0 and γ = 2(τ + t − to − ζ1 )/(t0 + r)
if −s − t0 ≤ ζ < τ + t − t0 . Hence using the above estimate and (5.13) we
see that
pδ ≤

Cδ γ
= Cδ γ−2α .
(C1 δ α )2

Above we have to assume α < γ/2 to have convergence pδ → 0 when δ → 0.
Finally, since we have denoted uδ = Tδ (mδ ), we can conclude that with
above choices for R(δ), γ and α




P Vδ ∈ Bζ1 (Uδ , R(δ)) = P U ∈ Bζ1 (Tδ (mδ ), R(δ)) | Mδ = mδ
≥ 1 − Cδ γ−2α → 1
when δ → 0.



POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

27

5.3. Discussion. Above, we have considered in the frequentist setting the
case when the solution u† is an element of H τ (N ) with τ > 0 and studied in Theorems 2 and 3 the convergence of the MAP estimators and the
contraction of the posterior distribution in L2 (N ).
In the Bayesian setting we have examined the case when the solution
is a realisation of the random variable U . In Theorems 1 and 4 we have
studied the the convergence of the MAP estimators and the contraction of
the posterior distribution in H ζ (N ) with various values of ζ.
In classical regularisation theory for linear inverse problems, one is usually
interested in the convergence of the optimisers of the minimisation problem
(2.5) to the true solution in the space H r (N ), r > τ + d/2, as the noise level
δ goes to zero. This gives restrictions to the measurement noise that can be
considered. Summarising, our above statistical considerations concern the
case where unknown and the noise are significantly less smooth than in the
standard setting of the regularisation theory. In the recent regularisation
theory inverse problems where the direct map A and the regularisation term
are non-linear have been studied extensively. It is interesting to ask how our
analysis on the contraction of the posterior distribution could be generalised
for such non-linear inverse problems that corresponds to non-Gaussian statistical problems.
Appendix A. Some examples of hypoelliptic operators
A linear partial differential operator P is hypoelliptic if for every distribution u such that P (u) is C ∞ smooth also u is C ∞ . Every elliptic operator
with smooths coefficients is hypoelliptic. The heat operator
P u(x, t) = ∂t u − k∆x u,

(x, t) ∈ Rd × R

and Kolmogorov operator [30, 19]
(A.1)

P u(x, y, t) = ∂xx u + x∂y u − ∂t u,

(x, y, t) ∈ R3

are examples of operators that are hypoelliptic but not elliptic. General
Kolmogorov type hypoelliptic diffusion operators are used e.g. in the theory
of kinetic equations, statistical physics and mathematical finance [17, 39].
The fact that (A.1) is hypoelliptic follows from Hörmander’s theorem on
hypoelliptic PDEs. Let (X0 , X1 , . . . , Xp ) be real C ∞ vector fields in the d
dimensional manifold N . If X and Y are two vector fields we define the
bracket of X and Y by
[X, Y ]f = X(Y f ) − Y (Xf ).
Note that [X, Y ] is a new vector field.
Definition 6 (Hörmander Condition). We say that the Hörmander condition is satisfied if the real C ∞ vector fields (X0 , X1 , . . . , Xp ), p ≤ n, in the
manifold N generate a Lie algebra of rank n = dim N at every point x ∈ N .
This means that the vector fields
Xj , [Xj1 , Xj2 ], [Xj1 , [Xj2 , Xj3 ]], . . .
span a space that has the same dimension n as the manifold N at every
point x ∈ N . Now we can formulate Hörmander’s classical theorem [20].

28

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

Theorem 5 (Hörmander’s theorem). The operator
P =

p
X

Xj2 + X0

j=1

defined on n dimensional manifold N is hypoelliptic if the vector fields
(X0 , X1 , . . . , Xp ) satisfy the Hörmander condition.
By writing
X1 = ∂x ,

X0 = x∂y − ∂t ,

[X1 , X0 ] = ∂y

we see that (A.1) is indeed hypoelliptic in R3 . Next we will give another,
important example of vector fields satisfying Hörmander’s condition
Example 4. Let us study Heisenberg group H. Let Γ be a discrete subgroup of Heisenberg group H such that H/Γ is compact,see e.g. [12] The
orthonormal frame on R3 is given by the Lie vector fields
1
X = ∂x + y∂t
2
1
Y = ∂y − x∂t
2
Z = ∂t
We can easily see that
[X, Y ] = Z.
Using Theorem 5 we get that the sub-Laplacian
1
P = (X 2 + Y 2 )
2
on H/Γ is hypoelliptic.
Example 5. One example of hypoelliptic inverse problem is the heat equation on a compact manifold N = M ×R where M is a closed two-dimensional
manifold. Note that in this paper we have considered the problem on a
compact manifold, that is, our results are applicable in the case when the
equation is periodic in time. We are interested in solving the heat sources
U (x, t) from the noisy measurements Mδ (x, t) of temperature T (x, t), that
is we want to solve U from
(A.2)
(A.3)

(∂t − ∆x )T (x, t) = U (x, t),
Mδ (x, t) = (∂t − ∆x )−1 U (x, t) + δE.

The operator A = (∂t − ∆x )−1 is not elliptic but it is hypoelliptic of type
(1,2).
Such situations arise in non-invasive monitoring. Consider, for example,
using a thermal camera to record video footage of a car with engine running. Let us model the metal surface of the car as a compact and closed
two-dimensional manifold M . The running engine produces heat which we
observe in the video data. The temperature on the car surface is modelled
as the solution T (x, t) defined on N = M × R. Equation (A.2) describes
the conduction of heat along the car surface. The effect of the engine is

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

29

simply modelled as the heat source term U (x, t); recovering U will provide
information about the state of the engine.
Appendix B. Computational example
Since the operator A does not have a continuous inverse operator L2 →
the condition number of the matrix approximation A of the operator
A grows when the discretisation is refined. This is the very reason why
regularisation is need in the (numerical) solutions of the inverse problems.
Next we demonstrate the above results numerically and consider twodimensional deblurring problem on T2 ,
L2 ,

M = AU + δE,
H −s ,

where E(ω) ∈
s > 1 a.s. is normalised white noise and A is elliptic
operator, smoothing of order 2,

(Au)(x) = F −1 (1 + |n|2 )−1 (Fu)(n) (x).
The true solution u† ∈ H 1 (T2 ), see Subsection 2.1.2 for the frequentist
interpretation, is a piecewise linear function presented in Figure 1. We
choose a priori distribution N (0, CU ) where CU = (I −∆)−1 . Then a random
draw U (ω) from the prior distribution belongs to H τ , where τ < 0, with
probability one. The Cameron-Martin space ZU of a measurable mapping
U : Ω → X is defined by

ZU = φ ∈ X kφk2ZU = hCU−1 φ, φiX ∗ ×X < ∞ .
Cameron-Martin space can also be defined as
\
ZU =
Y Y ⊂ H τ linear subspace, P({U ∈ Y }) = 1 .
The approximated solutions Uδ belongs to ZU and with the chosen a priori
distribution we have ZU = H 1 (T2 ).

Figure 1. On the left the original piecewise linear function
u† . On the right side the noiseless data m† = Au† .
Solving u from Au(x) = m(x) corresponds to the solution of ordinary
differential equation (1 − ∂x2 )m(x) = u(x) so A can be thought e.g. as a
blurring operator.

30

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

The approximated solution to the problem is
u†δ = (A∗ A + δ 2 (I − ∆))−1 A∗ m† .
We get from Subsection 2.1.2 that
lim Eu ku† − Uδ† kH ζ = 0

δ→0

when ζ < τ < 0. This behaviour can be seen even in numerical simulations
when the discretisation is fine enough, see Figure 2. In Figure 3 we have
compared the expected convergence rates given in formula (2.12) in Theorem 1 to the computational convergence rates. In the numerical simulations
in Figures 2 and 3 we see that for the test case presented in Figure 1 the
convergence u†δ → u† in different Sobolev spaces follows well the mean convergence predicted by Theorem 1.

Figure 2. Normalised errors c(ζ)ku† − u†δ kH ζ (T2 ) in logarithmic scale with different values of ζ. We use normalisation constants c(ζ) = 1/ku† − u†5·10−6 kH ζ . We observe that
u†δ ∈ H 1 does not converge to u† ∈ H 1 in H ζ when ζ > 0.
Acknowledgements. We would like to thank Petteri Piiroinen for helpful discussions. This work was supported by the Finnish Centre of Excellence
in Inverse Problems Research 2012-2017 (Academy of Finland CoE-project
284715). In addition, H.K. was supported by Emil Aaltonen Foundation
and EQUIP, grant EP/K034154/1, M.L. was supported by Academy of Finland, grant 273979, and S.S. was supported by Academy of Finland, project
141094.

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

31

Figure 3. Normalised errors c(ζ)ku† − u†δ kH ζ (T2 ) in logarithmic scale with different values of ζ. The normalized bounds (2.12) given in Theorem 1 for the expectations
c1 (ζ)EkU −Uδ kH ζ (T2 ) are plotted with dashed lines. The normalized errors c(ζ)ku† − u†δ kH ζ (T2 ) , for the example u† given
in Figure 1, are plotted with solid lines.
References
[1] S. Agapiou, S. Larsson, and A. M. Stuart, Posterior contraction rates for the
Bayesian approach to linear ill-posed inverse problems, Stochastic Processes and their
Applications, 123 (2013), pp. 3828–3860.
[2] S. Agapiou, A. M. Stuart, and Y.-X. Zhang, Bayesian posterior contraction
rates for linear severely ill-posed inverse problems, Journal of Inverse and Ill-posed
Problems, 22 (2014), pp. 297–321.
[3] V. Bogachev, Gaussian measures, vol. 62 of Mathematical Surveys and Monographs,
American Mathematical Society, Providence, RI, 1998.
[4] I. Castillo and R. Nickl, Nonparametric Bernstein–von Mises theorems in Gaussian white noise, The Annals of Statistics, 41 (2013), pp. 1999–2028.
[5]
, On the Bernstein–von Mises phenomenon for nonparametric Bayes procedures,
The Annals of Statistics, 42 (2014), pp. 1941–1969.
[6] L. Cavalier, Nonparametric statistical inverse problems, Inverse Problems, 24
(2008), p. 034004.
[7] M. Dashti, K. J. Law, A. M. Stuart, and J. Voss, MAP estimators and their
consistency in Bayesian nonparametric inverse problems, Inverse Problems, 29 (2013),
p. 095017.
[8] R. M. Dudley, Real analysis and probability, Mathematics Series, The Wadsworth
& Brooks/Cole, Pacific Grove, CA, 1989.
[9] H. Engl, M. Hanke, and A. Neubauer, Regularization of inverse problems,
Springer Netherlands, 1996.
[10] B. Fitzpatrick, Bayesian analysis in inverse problems, Inverse problems, 7 (1991),
pp. 675–702.
[11] J.-P. Florens and A. Simoni, Regularizing priors for linear inverse problems,
Econometric Theory, 32 (2016), pp. 71–121.
[12] G. Folland, Compact Heisenberg manifolds as CR manifolds, The Journal of Geometric Analysis, 14 (2004), pp. 521–532.

32

HANNE KEKKONEN, MATTI LASSAS AND SAMULI SILTANEN

[13] S. Ghosal, A review of consistency and convergence of posterior distribution, in
Varanashi Symposium in Bayesian Inference, Banaras Hindu University, 1997.
[14] S. Ghosal, J. K. Ghosh, and A. W. Van Der Vaart, Convergence rates of
posterior distributions, Annals of Statistics, 28 (2000), pp. 500–531.
[15] E. Giné and R. Nickl, Mathematical foundations of infinite-dimensional statistical
models, 2015.
[16] C. Groetsch, The Theory of Tikhonov Regularization for Fredholm Equations of the
First Kind, Pitman, London, 1984.
[17] B. Helffer and F. Nier, Hypoelliptic estimates and spectral theory for FokkerPlanck operators and Witten Laplacians, Springer, 2005.
[18] M. Hoffmann, J. Rousseau, and J. Schmidt-Hieber, On adaptive posterior concentration rates, The Annals of Statistics, 43 (2015), pp. 2259–2295.
[19] L. Hörmander, Hypoelliptic second order differential equations, Acta Mathematica,
119 (1967), pp. 147–171.
[20] L. Hörmander, The analysis of linear partial differential operators III: Pseudodifferential operators, Springer-Verlag, 1994.
[21] T.-M. Huang, Convergence rates for posterior distributions and adaptive estimation,
The Annals of Statistics, 32 (2004), pp. 1556–1593.
[22] J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems,
vol. 160 of Applied Mathematical Sciences, Springer Verlag, 2004.
[23] H. Kekkonen, M. Lassas, and S. Siltanen, Analysis of regularized inversion of
data corrupted by white Gaussian noise, Inverse Problems, 30 (2014), p. 045009.
[24] A. Kirsch, An introduction to the mathematical theory of inverse problems, SpringerVerlag, New York, 1996.
[25] B. Knapik and J.-B. Salomond, A general approach to posterior contraction in
nonparametric inverse problems, arXiv preprint arXiv:1407.0335, (2014).
[26] B. Knapik, B. Szabó, A. van der Vaart, and J. van Zanten, Bayes procedures
for adaptive inference in inverse problems for the white noise model, arXiv preprint
arXiv:1209.3628, (2012).
[27] B. Knapik, A. van Der Vaart, and J. Van Zanten, Bayesian inverse problems
with Gaussian priors, The Annals of Statistics, 39 (2011), pp. 2626–2657.
[28] B. Knapik, A. Van der Vaart, and J. Van Zanten, Bayesian recovery of the initial condition for the heat equation, Communications in Statistics-Theory and Methods, 42 (2013), pp. 1294–1313.
[29] V. Kolehmainen, M. Lassas, K. Niinimäki, and S. Siltanen, Sparsity-promoting
Bayesian inversion, Inverse Problems, 28 (2012), p. 025005.
[30] A. Kolmogoroff, Zufallige bewegungen (zur theorie der brownschen bewegung), Annals of Mathematics, (1934), pp. 116–117.
[31] S. Lasanen, Discretizations of generalized random variables with applications to inverse problems, PhD thesis, Ann. Acad. Sci. Fenn. Math. Diss., 2002.
[32] S. Lasanen, Non-Gaussian statistical inverse problems. part i: Posterior distributions., Inverse Problems & Imaging, 6 (2012).
[33] M. Lassas, E. Saksman, and S. Siltanen, Discretization-invariant Bayesian inversion and Besov space priors, Inverse Problems and Imaging, 3 (2009), pp. 87–122.
[34] M. Lassas and S. Siltanen, Can one use total variation prior for edge-preserving
Bayesian inversion?, Inverse Problems, 20 (2004), pp. 1537–1564.
[35] H. Leahu, On the Bernstein-von Mises phenomenon in the Gaussian white noise
model, Electronic Journal of Statistics, 5 (2011), pp. 373–404.
[36] M. Lehtinen, L. Päivärinta, and E. Somersalo, Linear inverse problems for
generalised random variables, Inverse Problems, 5 (1989), pp. 599–612.
[37] A. Mandelbaum, Linear estimators and measurable linear transformations on a
Hilbert space, Zeitung für Wahscheinlichkeitstheorie und verwandte Gebiete, 65
(1984), pp. 385–397.
[38] V. A. Morozov, Z. Nashed, and A. Aries, Methods for solving incorrectly posed
problems, Springer, 1984.
[39] A. Pascucci, PDE and martingale methods in option pricing, vol. 2, Springer, 2011.

POSTERIOR CONSISTENCY AND CONVERGENCE RATES FOR BAYESIAN INVERSION

33

[40] D. L. Phillips, A technique for the numerical solution of certain integral equations
of the first kind, Journal of the ACM (JACM), 9 (1962), pp. 84–97.
[41] P. Piiroinen, Statistical Measurements, Experiments and Applications, PhD thesis,
Ann. Acad. Sci. Fenn. Math. Diss, 2005.
[42] K. Ray, Bayesian inverse problems with non-conjugate priors, Electronic Journal of
Statistics, 7 (2013), pp. 2516–2549.
[43]
, Adaptive Bernstein-von Mises theorems in Gaussian white noise, (2014).
[44] I. A. Rozanov, Infinite-dimensional Gaussian distributions, no. 108, American
Mathematical Soc., 1971.
[45] Y. Rozanov, Markov random fields, Springer-Verlag, 1982.
[46] W. Rudin, Real and Complex Analysis, McGraw-Hill, third ed., 1987.
[47] T. Schuster, B. Kaltenbacher, B. Hofmann, and K. S. Kazimierski, Regularization methods in Banach spaces, vol. 10, Walter de Gruyter, 2012.
[48] X. Shen and L. Wasserman, Rates of convergence of posterior distributions, Annals
of Statistics, (2001), pp. 687–714.
[49] M. A. Shubin and S. I. Andersson, Pseudodifferential operators and spectral theory,
vol. 200, Springer, 1987.
[50] A. M. Stuart, Inverse problems: a Bayesian perspective, Acta Numerica, 19 (2010),
pp. 451–559.
[51] B. Szabo, A. van der Vaart, and H. van Zanten, Honest Bayesian confidence
sets for the l2-norm, Journal of Statistical Planning and Inference, 166 (2015), pp. 36–
51.
[52] B. Szabó, A. van der Vaart, J. van Zanten, et al., Frequentist coverage of
adaptive nonparametric Bayesian credible sets, The Annals of Statistics, 43 (2015),
pp. 1391–1428.
[53] A. N. Tikhonov, On the stability of inverse problems, in Dokl. Akad. Nauk SSSR,
vol. 39, 1943, pp. 195–198.
[54]
, Solution of incorrectly formulated problems and the regularization method, in
Soviet Mathematics Doklady, vol. 4, 1963, pp. 1035–1038.
[55] A. N. Tikhonov, Numerical methods for the solution of ill-posed problems, vol. 328,
Springer, 1995.
[56] A. W. van der Vaart, Asymptotic statistics, Cambridge series in statistical and
probabilistic mathematics, (2000).
[57] A. W. van der Vaart and J. H. van Zanten, Rates of contraction of posterior distributions based on Gaussian process priors, The Annals of Statistics, (2008),
pp. 1435–1463.
[58] S. J. Vollmer, Posterior consistency for Bayesian inverse problems through stability
and regression results, Inverse Problems, 29 (2013), p. 125011.
Department of Mathematics and Statistics P.O. Box 68, 00014 University
of Helsinki FINLAND

