tablesection algorithmsection

# _Data Assimilation and Inverse Problems_

Daniel Sanz-Alonso ${}^{\mbox{\musFontLarge\kern
0.70004pt\raisebox{1.50694pt}{\char 90\relax}\kern-0.1pt}{}}$ Andrew Stuart
${}^{\mbox{\musFontLarge\kern 0.70004pt\raisebox{1.50694pt}{\char
92\relax}\kern-0.1pt}{}}$ and Armeen Taeb ${}^{\mbox{\musFontLarge\kern
0.70004pt\raisebox{1.50694pt}{\char 94\relax}\kern-0.1pt}{}}$
${}^{\mbox{\musFontLarge\kern 0.70004pt\raisebox{1.50694pt}{\char
90\relax}\kern-0.1pt}{}}$ Department of Statistics  University of Chicago
${}^{\mbox{\musFontLarge\kern 0.70004pt\raisebox{1.50694pt}{\char
92\relax}\kern-0.1pt}{}}$ Department of Computational and Mathematical
Sciences  Caltech
${}^{\mbox{\musFontLarge\kern 0.70004pt\raisebox{1.50694pt}{\char
94\relax}\kern-0.1pt}{}}$ Department of Electrical Engineering  Caltech

## Introduction

Overview of the Notes

These notes were first developed out of Caltech course ACM 159 in Fall 2017,
and then substantially modified for the University of Chicago course STAT
31550 in Winter 2019. The aim of the notes is to provide a clear and concise
introduction to the subjects of Inverse Problems and Data Assimilation, and
their inter-relations, together with citations to some relevant literature in
this area

In its most basic form, inverse problem theory is the study of how to estimate
model parameters from data. Often the data provides indirect information about
these parameters, corrupted by noise. The theory of inverse problems, however,
is much richer than just parameter estimation. For example, the underlying
theory can be used to determine the effects of noisy data on the accuracy of
the solution; it can be used to determine what kind of observations are needed
to accurately determine a parameter; and it can be used to study the
uncertainty in a parameter estimate and, relatedly, is useful, for example, in
the design of strategies for control or optimization under uncertainty, and
for risk analysis. The theory thus has applications in many fields of science
and engineering.

To apply the ideas in these notes, the starting point is a mathematical model
mapping the unknown parameters to the observations: termed the “forward” or
“direct” problem, and often a subject of research in its own right. A good
forward model will not only identify how the data is dependent on parameters,
but also what sources of noise or model uncertainty are present in the
postulated relationship between unknown parameters and data. For example, if
the desired forward problem cannot be solved analytically, then the forward
model may be approximated by a simulation; in this case discretization may be
considered as a source of error. Once a relationship between model parameters,
sources of error, and data is clearly defined, the inverse problem of
estimating parameters from data can be addressed. The theory of inverse
problems can be separated into two cases: (1) the ideal case where data is not
corrupted by noise and is derived from a known perfect model; and (2) the
practical case where data is incomplete and imprecise. The first case is
useful for classifying inverse problems and determining if a given set of
observations can, in principle, provide exact solutions; this provides insight
into conditions needed for existence, uniqueness, and stability of a solution.
The second case is useful for the formulation of practical algorithms to learn
about parameters, and uncertainties in their estimates, and will be the focus
of these notes.

A model that has the properties: (a) a solution map exists, (b) is unique, and
(c) its behavior changes continuously with input (stability) is termed “well-
posed”. Conversely, a model lacking any of these properties is termed “ill-
posed”. Ill-posedness is present is many inverse problems, and mitigating it
is an extensive part of the subject. Out of the different approaches for
formulating an inverse problem, the Bayesian framework naturally offers the
ability to assess quality in parameter estimation, and also leads to a form of
well-posedness at the level of probability distributions describing the
solution. The goal of the Bayesian framework is to find a probability measure
that assigns a probability to each possible solution for a parameter $u$,
given the data $y$. Bayes formula states that

$\operatorname{\mathbb{P}}(u|y)=\frac{1}{\operatorname{\mathbb{P}}(y)}\operatorname{\mathbb{P}}(y|u)\operatorname{\mathbb{P}}(u).$

It enables calculation of the posterior probability on $u|y$,
$\operatorname{\mathbb{P}}(u|y)$, in terms of the product of the data
likelihood $\operatorname{\mathbb{P}}(y|u)$ and the prior information on the
parameter encoded in $\operatorname{\mathbb{P}}(u)$. The likelihood describes
the probability of the observed data $y$, if the input parameter were set to
be $u$; it is determined by the forward model, and the structure of the noise.
The normalization constant $\operatorname{\mathbb{P}}(y)$ ensures that
$\operatorname{\mathbb{P}}(u|y)$ is a probability measure. There are five
primary benefits to this framework: (a) it provides a clear theoretical
setting in which the forward model choice, noise model and a priori
information are explicit; (b) it provides information about the entire
solution space for possible input parameter choices; (c) it naturally leads to
quantification of uncertainty and risk in parameter estimates; (d) it is
generalizable to a wide class of inverse problems, in finite and infinite
dimension and comes with a well-posedness theory useful in these contexts.

The first half of the notes is dedicated to studying the Bayesian framework
for inverse problems. Techniques such as importance sampling and Markov Chain
Monte Carlo (MCMC) methods are introduced; these methods have the desirable
property that in the limit of an infinite number of samples they reproduce the
full posterior distribution. Since it is often computationally intensive to
implement these methods, especially in high dimensional problems, approximate
techniques such as approximating the posterior by a Dirac or a Gaussian
distribution are discussed.

The second half of the notes covers data assimilation. This refers to a
particular class of inverse problems in which the unknown parameter is the
initial condition of a dynamical system, and in the stochastic dynamics case
the subsequent states of the system, and the data comprises partial and noisy
observations of that (possibly stochastic) dynamical system. A primary use of
data assimilation is in forecasting, where the purpose is to provide better
future estimates than can be obtained using either the data or the model
alone. All the methods from the first half of the course may be applied
directly, but there are other new methods which exploit the Markovian
structure to update the state of the system sequentially, rather than to learn
about the initial condition. (But of course knowledge of the initial condition
may be used to inform the state of the system at later times.) We will also
demonstrate that methods developed in data assimilation may be employed to
study generic inverse problems, by introducing an artificial time to generate
a sequence of probability measures interpolating from the prior to the
posterior.

Notation

Throughout the notes we use $\mathbb{N}$ to denote the positive integers
$\\{1,2,3,\cdots\\}$ and $\mathbb{Z}^{+}$ to denote the non-negative integers
$\mathbb{N}\cup\\{0\\}=\\{0,1,2,3,\cdots\\}.$ The matrix $I_{d}$ denotes the
identity on $\mathbb{R}^{d}.$ We use $|\cdot|$ to denote the Euclidean norm
corresponding to the inner-product $\langle\cdot,\cdot\rangle.$ A square
matrix $A$ is positive definite (resp. positive semi-definite) if the
quadratic form $\langle u,Au\rangle$ is positive (resp. non-negative) for all
$u\neq 0$. By $|\cdot|_{A}$ we denote the weighted norm defined by
$|v|_{A}^{2}=v^{T}A^{-1}v$. The corresponding weighted Euclidean inner-product
is given by $\langle\cdot,\cdot\rangle_{A}$ and
$\langle\cdot,A^{-1}\cdot\rangle.$ We use $\otimes$ to denote the outer
product between two vectors: $(a\otimes b)c=\langle b,c\rangle a.$ We let
$B(u,\delta)$ denote the open ball of radius $\delta$ at $u$, in the Euclidean
norm.

Throughout we denote by
$\operatorname{\mathbb{P}}(\cdot),\operatorname{\mathbb{P}}(\cdot|\cdot)$ the
pdf of a random variable and its conditional pdf, respectively. We write

$\rho(f)=\operatorname{\mathbb{E}}^{\rho}[f]=\int_{\mathbb{R}^{d}}f(u)\rho(u)du$

to denote expectation of $f:\mathbb{R}^{d}\mapsto\mathbb{R}$ with respect to
probability measure with probability density function (pdf) $\rho$ on
$\mathbb{R}^{d}.$ The distribution of the random variables in this book will
often have density with respect to Lebesgue measure, but occasional use of
Dirac masses will be required; we will use the notational convention that
Dirac mass at point $v$ has “density” $\delta(\cdot-v)$ or
$\delta_{v}(\cdot).$ When a random variable $u$ has pdf $\rho$ we will write
$u\sim\rho.$ We use $\Rightarrow$ to denote weak convergence of probability
measures, that is, $\rho_{n}\Rightarrow\rho$ if $\rho_{n}(f)\to\rho(f)$ for
all bounded and continuous $f:\mathbb{R}^{d}\mapsto\mathbb{R}$.

Acknowledgements

These notes were first created out of Caltech course ACM 159 in Fall 2017 and
were further developed for the University of Chicago course STAT 31550 in
Winter 2019. The notes were created in LaTeX by the students in ACM 159, based
on lectures presented by the instructor Andrew Stuart, and on input from the
course TA Armeen Taeb. The individuals responsible for the notes listed in
alphabetic order are: Blancquart, Paul; Cai, Karena; Chen, Jiajie; Cheng,
Richard; Cheng, Rui; Feldstein, Jonathan; Huang, De; Idíni, Benjamin;
Kovachki, Nikola; Lee, Marcus; Levy, Gabriel; Li, Liuchi; Muir, Jack; Ren,
Cindy; Seylabi, Elnaz, Schäfer, Florian; Singhal, Vipul; Stephenson, Oliver;
Song, Yichuan; Su, Yu; Teke, Oguzhan; Williams, Ethan; Wray, Parker; Zhan,
Eric; Zhang, Shumao; Xiao, Fangzhou. Furthermore, the following students have
added content beyond the class materials: Parker Wray – the Overview, Jiajie
Chen – alternative proof of Theorem 1.10 and proof idea for Theorem 14.3,
Fangzhou Xiao – numerical simulation of prior, likelihood $\&$ posterior,
Elnaz Seylabi $\&$ Fangzhou Xiao – catching typographical errors in a draft of
these notes, Cindy Ren – numerical simulations to enhance understanding of
importance sampling in Examples 6.2 and 6.5, Cindy Ren $\&$ De Huang –
improving the constants in Theorem 6.3 regarding the approximation error of
importance sampling, Richard Cheng $\&$ Florian Schäfer – illustrations to
enhance understanding of the coupling argument used to study convergence of
MCMC algorithms by presenting the finite state-space case, and Ethan Williams
$\&$ Jack Muir – numerical simulations and illustrations of Ensemble Kalman
Filter and Extended Kalman Filter. We would also like to thank Tapio Helin
(Helsinki) who used the notes in his own course and provided very helpful
feedback on an early draft.

The work of AS has been funded by the EPSRC (UK), ERC (Europe) and by AFOSR,
ARL, NIH, NSF and ONR (USA). The work of DSA is funded by NSF. This funded
research has helped to shape the presentation of the material here and is
gratefully acknowledged.

Warning

These are rough notes, far from being perfected. They are likely to contain
mathematical errors, incomplete bibliographical information, inconsistencies
in notation and typographical errors. We hope that the notes are nonetheless
useful. Please contact the authors with any feedback from typos, through
mathematical errors and bibliographical omissions, to comments on the
structural organization of the material.

###### Contents

  1. 1 Bayesian Inverse Problems and Well-Posedness
    1. 1.1 Formulation of Bayesian Inverse Problems
    2. 1.2 Formula for Posterior pdf: Bayes Theorem
    3. 1.3 Well-posedness of Bayesian Inverse Problems
      1. 1.3.1 Metrics on Probability Densities
      2. 1.3.2 Approximation Theorem
      3. 1.3.3 Example
    4. 1.4 Discussion and Bibliography
  2. 2 The Linear-Gaussian Setting
    1. 2.1 Derivation of the Posterior Distribution
    2. 2.2 Small Noise Limit of the Posterior Distribution
      1. 2.2.1 Overdetermined Case
      2. 2.2.2 Determined Case
      3. 2.2.3 Underdetermined Case
    3. 2.3 Discussion and Bibliography
  3. 3 Optimization Perspective
    1. 3.1 The Setting
    2. 3.2 Theory
    3. 3.3 Examples
    4. 3.4 Discussion and Bibliography
  4. 4 The Gaussian Approximation
    1. 4.1 The Kullback-Leibler Divergence
    2. 4.2 Best Gaussian Fit By Minimizing $d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$
    3. 4.3 Best Gaussian Fit By Minimizing $d_{\mbox{\tiny{\rm KL}}}(\pi\|p)$
    4. 4.4 Comparison between $d_{\mbox{\tiny{\rm KL}}}(\pi\|p)$ and $d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$
    5. 4.5 Variational Formulation of Bayes Theorem
    6. 4.6 Discussion and Bibliography
  5. 5 Monte Carlo and Importance Sampling
    1. 5.1 Monte Carlo Sampling
    2. 5.2 Importance Sampling
    3. 5.3 Discussion and Bibliography
  6. 6 Markov Chain Monte Carlo
    1. 6.1 The Idea Behind MCMC
    2. 6.2 The Metropolis-Hastings Algorithm
    3. 6.3 Invariance of the Target Distribution $\pi$
      1. 6.3.1 Detailed Balance and its Implication
      2. 6.3.2 Detailed Balance and the Metropolis-Hastings Algorithm
    4. 6.4 Convergence to the Target Distribution
      1. 6.4.1 Finite State Space
      2. 6.4.2 The pCN Method
    5. 6.5 Discussion and Bibliography
  7. 7 Filtering and Smoothing Problems and Well-Posedness
    1. 7.1 Formulation of Filtering and Smoothing Problems
    2. 7.2 The Smoothing Problem
      1. 7.2.1 Formula for pdf of the Smoothing Problem
      2. 7.2.2 Well-Posedness of the Smoothing Problem
    3. 7.3 The Filtering Problem
      1. 7.3.1 Formula for pdf of the Filtering Problem
      2. 7.3.2 Well-Posedness of the Filtering Problem
    4. 7.4 Discussion and Bibliography
  8. 8 The Linear-Gaussian Setting
    1. 8.1 Kalman Filter
      1. 8.1.1 Kalman Filter: Alternative Formulation
      2. 8.1.2 Optimization Perspective: Mean of Kalman Filter
      3. 8.1.3 Optimality of Kalman Filter
    2. 8.2 Kalman Smoother
      1. 8.2.1 Defining Linear System
      2. 8.2.2 Solution of the Linear System
    3. 8.3 Discussion and Bibliography
  9. 9 Optimization for Filtering and Smoothing: 3DVAR and 4DVAR
    1. 9.1 The Setting
    2. 9.2 3DVAR
    3. 9.3 4DVAR
    4. 9.4 Discussion and Bibliography
  10. 10 The Extended and Ensemble Kalman Filters
    1. 10.1 The Setting
    2. 10.2 The Extended Kalman Filter
    3. 10.3 Ensemble Kalman Filter
    4. 10.4 Derivation of EnKF
    5. 10.5 Subspace Property of EnKF
    6. 10.6 Discussion and Bibliography
  11. 11 Particle Filter
    1. 11.1 The Setting
    2. 11.2 The Bootstrap Particle Filter
    3. 11.3 Bootstrap Particle Filter Convergence
    4. 11.4 The Bootstrap Particle Filter as a Random Dynamical System
    5. 11.5 Discussion and Bibliography
  12. 12 Optimal Particle Filter
    1. 12.1 Introduction
    2. 12.2 The Bootstrap and Optimal Particle Filters Compared
    3. 12.3 Implementation of Optimal Particle Filter: Linear Observation Setting
    4. 12.4 “Optimality” of the Optimal Particle Filter
    5. 12.5 Particle Filters for High Dimensions
    6. 12.6 Filtering Overview
      1. 12.6.1 Dynamical Model
      2. 12.6.2 Data Model
      3. 12.6.3 Interaction of the Dynamics and Data
      4. 12.6.4 Summary of Applicability of Discrete Filtering Methods
    7. 12.7 Discussion and Bibliography
  13. 13 Filtering Approach to the Inverse Problem
    1. 13.1 General Formulation
    2. 13.2 Ensemble Kalman Inversion
    3. 13.3 Linking Ensemble Kalman Inversion and SMC
      1. 13.3.1 Continuous Time Limit
      2. 13.3.2 Linear Setting
    4. 13.4 Discussion and Bibliography

###  1 Bayesian Inverse Problems and Well-Posedness

In this chapter we introduce the Bayesian approach to inverse problems. We
will show that the Bayesian formulation leads to a form of well-posedness:
small perturbations of the forward model or the observed data translate into
small perturbations of the Bayesian solution. Well-posedness will be
established in Hellinger distance. The total variation and Hellinger distances
between probability densities play an important role in the analysis of
Bayesian methodology. We introduce both distances in this chapter, as well as
some characterizations and bounds between them that will be used throughout
this book.

####  1.1 Formulation of Bayesian Inverse Problems

We consider the following setting. We let $G:\mathbb{R}^{d}\to\mathbb{R}^{k}$
be a forward model and aim to recover an unknown parameter
$u\in\mathbb{R}^{d}$ from data $y\in\mathbb{R}^{k}$ given by

$y=G(u)+\eta,$ (1.1)

where $\eta\in\mathbb{R}^{k}$ represents observation noise. We view
$(u,y)\in\mathbb{R}^{d}\times\mathbb{R}^{k}$ as a random variable, whose
distribution is specified by means of the following assumption.

######  Assumption 1.1.

We make the following probabilistic assumptions:

  * •

$u\sim\rho(u),u\in\mathbb{R}^{d}$.

  * •

$\eta\sim\nu(\eta),\eta\in\mathbb{R}^{k}$.

  * •

$u$ and $\eta$ are independent, written $u\perp\eta.$

Here $\rho$ and $\nu$ describe the (Lebesgue) probability density functions
(pdfs) of the random variables $u$ and $\eta$ respectively. Then $\rho(u)$ is
called the _prior_ pdf, $y|u\sim\nu\bigl{(}y-G(u)\bigr{)}$ for each fixed
$u\in\mathbb{R}^{d}$ determines the _likelihood_ function. In this
probabilistic perspective, the solution to the inverse problem is the
conditional distribution of $u$ given $y$, which is called the _posterior_
distribution and will be denoted by $u|y\sim\pi^{y}(u).$ From the posterior
pdf one can infer parameter values that are consistent with both the data and
the prior pdf. The posterior also contains information on the uncertainty
remaining in the parameter recovery: for instance, large posterior covariance
may indicate that the data contains insufficient information to recover the
input parameter.

####  1.2 Formula for Posterior pdf: Bayes Theorem

Bayes theorem is a bridge connecting the prior, the likelihood and the
posterior.

######  Theorem 1.2 (Bayes theorem).

Let Assumption 1.1 hold and assume that

$Z=Z(y):=\int_{\mathbb{R}^{d}}\nu\bigl{(}y-G(u)\bigr{)}\rho(u)du>0.$

Then $u|y\sim\pi^{y}(u),$ where

$\pi^{y}(u)=\frac{1}{Z}\nu\bigl{(}y-G(u)\bigr{)}\rho(u).$ (1.2)

###### Proof.

Denote by $\operatorname{\mathbb{P}}(\cdot)$ the pdf of a random variable and
by $\operatorname{\mathbb{P}}(\cdot|\cdot)$ its conditional conditional pdf.
We have

$\displaystyle\operatorname{\mathbb{P}}(u,y)$
$\displaystyle=\operatorname{\mathbb{P}}(u|y)\operatorname{\mathbb{P}}(y),\textrm{
if }\operatorname{\mathbb{P}}(y)>0,$
$\displaystyle\operatorname{\mathbb{P}}(u,y)$
$\displaystyle=\operatorname{\mathbb{P}}(y|u)\operatorname{\mathbb{P}}(u),\textrm{
if }\operatorname{\mathbb{P}}(u)>0.$

Note that the marginal pdf on $y$ is given by

$\operatorname{\mathbb{P}}(y)=\int_{\mathbb{R}^{d}}\operatorname{\mathbb{P}}(u,y)du,$

and similarly for $\operatorname{\mathbb{P}}(u).$ Assume
$\operatorname{\mathbb{P}}(y)>0$. Then

$\operatorname{\mathbb{P}}(u|y)=\frac{1}{\operatorname{\mathbb{P}}(y)}\operatorname{\mathbb{P}}(y|u)\operatorname{\mathbb{P}}(u)=\frac{1}{Z}\nu\bigl{(}y-G(u)\bigr{)}\rho(u)$
(1.3)

for both $\operatorname{\mathbb{P}}(u)>0$ and
$\operatorname{\mathbb{P}}(u)=0$. Here we remark that

$\operatorname{\mathbb{P}}(y)=Z=\int_{\mathbb{R}^{d}}Z\operatorname{\mathbb{P}}(u|y)du=\int_{\mathbb{R}^{d}}\nu\bigl{(}y-G(u)\bigr{)}\rho(u)du>0.$

Thus the assumption that $\operatorname{\mathbb{P}}(y)>0$ is justified and the
desired result follows from equation (1.3). ∎

We will often denote the likelihood function by
$g(u):=\nu\bigl{(}y-G(u)\bigr{)}$; we then write

$\pi(u)=\frac{1}{Z}g(u)\rho(u),$

omitting the data $y$ in both the likelihood function and the posterior pdf.

######  Remark 1.3.

The proof of Theorem 1.2 shows that in order to apply Bayes formula (1.2) one
needs to guarantee that the normalizing constant
$\operatorname{\mathbb{P}}(y)=Z$ is positive. In other words, the marginal
density of the observed data $y$ needs to be positive, i.e. the observed data
needs to be consistent with the probabilistic assumptions (1.1). From now on
it will be assumed that $\operatorname{\mathbb{P}}(y)=Z>0$ without further
notice.

The posterior distribution $\pi^{y}(u)$ contains all the knowledge on the
parameter $u$ available in the prior and the data. In applications it is often
useful, however, to summarize the posterior distribution through a few
numerical values or through parameter regions of prescribed posterior
probability (known as credible intervals). Summarizing the posterior is
particularly important if the parameter is high dimensional, since then
visualizing the posterior or detecting regions of high posterior probability
is nontrivial. Two natural numerical summaries are the posterior mean and the
posterior mode.

######  Definition 1.4.

The posterior mean estimator of $u$ given data $y$ is the mean of the
posterior distribution:

$u_{\mbox{\tiny{\rm PM}}}=\int_{\mathbb{R}^{d}}u\pi^{y}(u)\,du.$

The maximum a posteriori (MAP) estimator of $u$ given data $y$ is the mode of
the posterior distribution $\pi^{{y}}({u})$, defined as

$u_{\mbox{\tiny{\rm MAP}}}=\arg\max_{{u}\in\mathbb{R}^{d}}\pi^{{y}}({u}).$

The posterior mean and the MAP estimators already suggest the importance of
computing maxima and integrals in the practical implementation of Bayesian
formulations of inverse problems and data assimilation. For this reason,
optimization and sampling will play an important role in this book. An
alternative way to make Bayesian formulations tractable is to approximate the
posterior by a simple distribution, often a Gaussian or a combination of Dirac
masses. An optimization perspective for inverse problems and data assimilation
will be studied in Chapters 3 and 7, respectively, and the use of Gaussian
approximations will be discussed in Chapters 4 and 10; Dirac approximations
constructed via sampling will be studied in Chapters 5 and 6 (inverse
problems) and in Chapters 11 and 12 (data assimilation).

We next consider two simple examples of a direct application of Bayes formula.

######  Example 1.5.

Let $d=k=1$, $\eta\sim\nu=\mathcal{N}(0,\gamma^{2})$ and

$\rho(u)=\begin{cases}\frac{1}{2},&u\in(-1,1),\\\ 0,&u\in(-1,1)^{c}.\\\
\end{cases}$

Suppose that the observation is generated by $y=u+\eta$. Using Bayes Theorem
1.2, we derive the posterior

$\pi^{y}(u)=\begin{cases}\frac{1}{2Z}\exp(-\frac{1}{2\gamma^{2}}|y-u|^{2}),&u\in(-1,1),\\\
0,&u\in(-1,1)^{c},\\\ \end{cases}$

where $Z$ is a normalizing constant ensuring
$\int_{\mathbb{R}}\pi^{y}(u)du=1$. The support of $\pi^{y}$, i.e. $(-1,1)$,
agrees with the support of the prior $\rho$. Now we find the MAP estimator.
From the explicit formula for $\pi^{y}$, we have

$u_{\mbox{\tiny{\rm
MAP}}}=\arg\max_{u\in\mathbb{R}}\pi^{y}(u)=\begin{cases}y&\textrm{if
}y\in(-1,1),\\\ -1&\textrm{if }y\leq-1,\\\ 1&\textrm{if }y\geq 1.\\\
\end{cases}$

In this example, the prior on $u$ is supported on $(-1,1)$ and the posterior
on $u|y$ is supported on $(-1,1)$. If the data lies in $(-1,1)$ then the MAP
estimator is the data itself; otherwise it is the extremal point of the prior
support which matches the sign of the data. The posterior mean is

$u_{\mbox{\tiny{\rm
PM}}}=\frac{1}{2Z}\int_{-1}^{1}u\exp\Bigl{(}-\frac{1}{2\gamma^{2}}|y-u|^{2}\Bigr{)}\,du,$

which may be approximated for instance by using the sampling methods described
in Chapters 5 and 6.

The following example illustrates once again the application of Bayes formula,
and shows that the posterior may be concentrated near a low-dimensional
manifold of the input parameter space $\mathbb{R}^{d}$. In such a case it is
important to understand the geometry of the support of the posterior density,
which cannot be captured by point estimation or Gaussian approximations.

######  Example 1.6.

Let $d=2,k=1,\rho\in
C(\mathbb{R}^{2},\mathbb{R}^{+}),0<\rho(u)\leq\rho_{\max}<+\infty$, for all
$u\in\mathbb{R}^{2}$ and

$y=G(u)+\eta=u_{1}^{2}+u_{2}^{2}+\eta,\quad\eta\sim\nu=\mathcal{N}(0,\gamma^{2}),\quad
0<\gamma\ll 1.$

Assume that $y>0$. Using Bayes theorem, we obtain the posterior distribution

$\pi^{y}(u)=\frac{1}{Z}\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2\gamma^{2}}|u_{1}^{2}+u_{2}^{2}-y|^{2}}\right)\rho(u).$

We now show that the posterior concentrates on a manifold: the circumference
$\\{u\in\mathbb{R}^{2}:u_{1}^{2}+u_{2}^{2}=y\\}.$ Denote by
$A^{\pm}:=\\{u\in\mathbb{R}^{2}:|u_{1}^{2}+u_{2}^{2}-y|^{2}\leq\gamma^{2\pm\delta}\\}$,
for some $\delta\in(0,2)$, and let $\rho_{\min}=\inf_{u\in B}\rho(u)$, where
$B$ is the closed ball of radius $2\sqrt{y}$ centered at the origin. Since
$\rho(u)$ is positive and continuous and $B$ is compact, $\rho_{\min}>0$. Let
$u^{+}\in A^{+}\subset B,u^{-}\in(A^{-})^{c}$. Taking the small noise limit
yields

$\frac{\pi^{y}(u^{+})}{\pi^{y}(u^{-})}\geq\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}\gamma^{\delta}+\frac{1}{2}\gamma^{-\delta}}\right)\frac{\rho_{\min}}{\rho_{\max}}\to\infty,\
\textrm{as }\gamma\to 0^{+}.$

Therefore, conditional on $y>0$, the posterior $\pi^{y}$ concentrates, as
$\gamma\to 0,$ on the circumference with radius $\sqrt{y}$.

Figure 1: The posterior measure concentrates on a circumference with radius
$\sqrt{y}$. Here, the blue shadow area is $A^{+}$ and the green shadow area is
$(A^{-})^{c}$.

####  1.3 Well-posedness of Bayesian Inverse Problems

In this section we show that the Bayesian formulation of inverse problems
leads to a form of well-posedness. More precisely, we study the sensitivity of
the posterior pdf to perturbations of the forward model $G.$ In many inverse
problems the ideal forward model $G$ is not accessible but can be approximated
by some computable $G_{\delta};$ consequently $\pi^{y}$ is replaced by
$\pi_{\delta}^{y}$. An example that is often found in applications, to which
the theory contained herein may be generalized, is when $G$ is an operator
acting on an infinite-dimensional space which is approximated, for the
purposes of computation, by some finite-dimensional operator $G_{\delta}$. We
seek to prove that, under certain assumptions, the small difference between
$G$ and $G_{\delta}$ (forward error) leads to similarly small difference
between $\pi^{y}$ and $\pi^{y}_{\delta}$ (inverse error):

Meta Theorem: Well-posedness

$|G-G_{\delta}|=O(\delta)\ \Longrightarrow\
d(\pi^{y},\pi_{\delta}^{y})=O(\delta),$

for small enough $\delta>0$ and some metric $d(\cdot,\cdot)$ on probability
densities.

This result will be formalized in Theorem 1.14 below, which shows that the
$O(\delta)$-convergence of $\pi_{\delta}^{y}$ with respect to some distance
$d(\cdot,\cdot)$ can be guaranteed under certain assumptions on the
likelihood. We will conclude the chapter by showing an example where these
assumptions hold true. In order to discuss these issues we will need to
introduce metrics on probability densities.

#####  1.3.1 Metrics on Probability Densities

Here we introduce the total variation and the Hellinger distance, both of
which have been used to show well-posedness results. In this chapter we will
use the Hellinger distance to establish well-posedness of Bayesian inverse
problems, and in Chapter 7 we employ the total variation distance to establish
well-posedness of Bayesian formulations of filtering and smoothing in data
assimilation.

######  Definition 1.7.

The _total variation distance_ between two probability densities
$\pi,\pi^{\prime}$ is defined by

$d_{\mbox{\tiny{\rm
TV}}}(\pi,\pi^{\prime})\coloneqq\frac{1}{2}\int|\pi(u)-\pi^{\prime}(u)|du=\frac{1}{2}\|\pi-\pi^{\prime}\|_{L^{1}}.$

The _Hellinger distance_ between two probability densities $\pi,\pi^{\prime}$
is defined by

$d_{\mbox{\tiny{\rm
H}}}(\pi,\pi^{\prime})\coloneqq\Big{(}\frac{1}{2}\int|\sqrt{\pi(u)}-\sqrt{\pi^{\prime}(u)}|^{2}du\Big{)}^{1/2}=\frac{1}{\sqrt{2}}\|\sqrt{\pi}-\sqrt{\pi^{\prime}}\|_{L^{2}}.$

In the rest of this subsection we will establish bounds between the Hellinger
and total variation distance, and show how both distances can be used to bound
the difference of expected values computed with two different densities; these
results will be used in subsequent chapters. Before doing so, the next lemma
motivates our choice of normalization constant $1/2$ for total variation
distance and $1/\sqrt{2}$ for Hellinger distance: they are chosen so that the
maximum possible distance between two densities is one.

######  Lemma 1.8.

For any probability densities $\pi,\pi^{\prime}$,

$0\leq d_{\mbox{\tiny{\rm TV}}}(\pi,\pi^{\prime})\leq 1,\quad 0\leq
d_{\mbox{\tiny{\rm H}}}(\pi,\pi^{\prime})\leq 1.$

###### Proof.

The lower bounds follow immediately from the definitions. We only need to
prove the upper bounds:

$d_{\mbox{\tiny{\rm
TV}}}(\pi,\pi^{\prime})=\frac{1}{2}\int|\pi(u)-\pi^{\prime}(u)|du\leq\frac{1}{2}\int\pi(u)du+\frac{1}{2}\int\pi^{\prime}(u)du=1,$
$\displaystyle d_{\mbox{\tiny{\rm H}}}(\pi,\pi^{\prime})=$ $\displaystyle\
\biggl{(}\frac{1}{2}\int\Bigl{|}\sqrt{\pi(u)}-\sqrt{\pi^{\prime}(u)}\Bigr{|}^{2}du\biggr{)}^{1/2}$
$\displaystyle=$ $\displaystyle\
\bigg{(}\frac{1}{2}\int\Big{(}\pi(u)+\pi^{\prime}(u)-2\sqrt{\pi(u)\pi^{\prime}(u)}\,\,\Big{)}du\bigg{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle\
\Big{(}\frac{1}{2}\int\big{(}\pi(u)+\pi^{\prime}(u)\big{)}du\Big{)}^{1/2}$
$\displaystyle=$ $\displaystyle\ 1.$

∎

The proofs also show that $\pi$ and $\pi^{\prime}$ have total variation and
Hellinger distance equal to one if and only if they have disjoint supports,
that is, if $\int\pi(u)\pi^{\prime}(u)du=0.$

The following result gives bounds between total variation and Hellinger
distance.

######  Lemma 1.9.

For any probability densities $\pi,\pi^{\prime}$,

$\frac{1}{\sqrt{2}}d_{\mbox{\tiny{\rm TV}}}(\pi,\pi^{\prime})\leq
d_{\mbox{\tiny{\rm H}}}(\pi,\pi^{\prime})\leq\sqrt{d_{\mbox{\tiny{\rm
TV}}}(\pi,\pi^{\prime})}.$

###### Proof.

Using the Cauchy–Schwartz inequality

$\displaystyle d_{\mbox{\tiny{\rm TV}}}(\pi,\pi^{\prime})=$ $\displaystyle\
\frac{1}{2}\int\Big{|}\sqrt{\pi(u)}-\sqrt{\pi^{\prime}(u)}\Big{|}\Big{|}\sqrt{\pi(u)}+\sqrt{\pi^{\prime}(u)}\Big{|}du$
$\displaystyle\leq$ $\displaystyle\
\biggl{(}\frac{1}{2}\int\Big{|}\sqrt{\pi(u)}-\sqrt{\pi^{\prime}(u)}\Big{|}^{2}du\biggr{)}^{1/2}\biggl{(}\frac{1}{2}\int\Big{|}\sqrt{\pi(u)}+\sqrt{\pi^{\prime}(u)}\Big{|}^{2}du\biggr{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle\ d_{\mbox{\tiny{\rm
H}}}(\pi,\pi^{\prime})\biggl{(}\frac{1}{2}\int\big{(}2\pi(u)+2\pi^{\prime}(u)\big{)}du\biggr{)}^{1/2}$
$\displaystyle=$ $\displaystyle\ \sqrt{2}d_{\mbox{\tiny{\rm
H}}}(\pi,\pi^{\prime}).$

Notice that
$|\sqrt{\pi(u)}-\sqrt{\pi^{\prime}(u)}|\leq|\sqrt{\pi(u)}+\sqrt{\pi^{\prime}(u)}|$
since $\sqrt{\pi(u)},\sqrt{\pi^{\prime}(u)}\geq 0$. Thus we have

$\displaystyle d_{\mbox{\tiny{\rm H}}}(\pi,\pi^{\prime})=$ $\displaystyle\
\biggl{(}\frac{1}{2}\int\Big{|}\sqrt{\pi(u)}-\sqrt{\pi^{\prime}(u)}\Big{|}^{2}du\biggr{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle\
\biggl{(}\frac{1}{2}\int\Big{|}\sqrt{\pi(u)}-\sqrt{\pi^{\prime}(u)}\Big{|}\Big{|}\sqrt{\pi(u)}+\sqrt{\pi^{\prime}(u)}\Big{|}du\biggr{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle\
\Big{(}\frac{1}{2}\int\big{|}\pi(u)-\pi^{\prime}(u)\big{|}du\Big{)}^{1/2}$
$\displaystyle=$ $\displaystyle\ \sqrt{d_{\mbox{\tiny{\rm
TV}}}(\pi,\pi^{\prime})}.$

∎

The following two lemmas show that if two densities are close in total
variation or in Hellinger distance, expectations computed with respect to both
densities are also close. In addition, the following lemma also provides a
useful characterization of the total variation distance which will be used
repeatedly throughout this book.

######  Lemma 1.10.

Let $f$ be a function such that
$\sup_{u\in\mathbb{R}^{d}}|f(u)|=:|f|_{\infty}<\infty$, then

$\big{|}\operatorname{\mathbb{E}}^{\pi}[f]-\operatorname{\mathbb{E}}^{\pi^{\prime}}[f]\big{|}\leq
2|f|_{\infty}d_{\mbox{\tiny{\rm TV}}}(\pi,\pi^{\prime}).$

Moreover, the following variational characterization of the total variation
distance holds:

$d_{\mbox{\tiny{\rm
TV}}}(\pi,\pi^{\prime})=\frac{1}{2}\sup_{|f|_{\infty}<1}\big{|}\operatorname{\mathbb{E}}^{\pi}[f]-\operatorname{\mathbb{E}}^{\pi^{\prime}}[f]\big{|}.$
(1.4)

###### Proof.

For the first part of the lemma, note that

$\displaystyle\Big{|}\operatorname{\mathbb{E}}^{\pi}[f]-\operatorname{\mathbb{E}}^{\pi^{\prime}}[f]\Big{|}=$
$\displaystyle\
\Big{|}\int_{\mathbb{R}^{d}}f(u)\bigl{(}\pi(u)-\pi^{\prime}(u)\bigr{)}du\Big{|}$
$\displaystyle\leq$ $\displaystyle\
2|f|_{\infty}\cdot\frac{1}{2}\int_{\mathbb{R}^{d}}|\pi(u)-\pi^{\prime}(u)|du$
$\displaystyle=$ $\displaystyle\ 2|f|_{\infty}d_{\mbox{\tiny{\rm
TV}}}(\pi,\pi^{\prime}).$

This in particular shows that, for any $f$ with $|f|_{\infty}=1,$

$d_{\mbox{\tiny{\rm
TV}}}(\pi,\pi^{\prime})\geq\frac{1}{2}\big{|}\operatorname{\mathbb{E}}^{\pi}[f]-\operatorname{\mathbb{E}}^{\pi^{\prime}}[f]\big{|}.$

Our goal now is to show a choice of $f$ with $|f|_{\infty}=1$ that achieves
equality. Define $f(u):=\text{sign}\Bigl{(}\pi(u)-\pi^{\prime}(u)\Bigr{)}$, so
that $f(u)\Bigl{(}\pi(u)-\pi^{\prime}(u)\Bigr{)}=|\pi(u)-\pi^{\prime}(u)|.$
Then, $|f|_{\infty}=1$ and

$\displaystyle d_{\mbox{\tiny{\rm TV}}}(\pi,\pi^{\prime})$
$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{d}}|\pi(u)-\pi^{\prime}(u)|\,du$
$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{d}}f(u)\Bigl{(}\pi(u)-\pi^{\prime}(u)\Bigr{)}\,du$
$\displaystyle=\frac{1}{2}\Big{|}\operatorname{\mathbb{E}}^{\pi}[f]-\operatorname{\mathbb{E}}^{\pi^{\prime}}[f]\Big{|}.$

This completes the proof of the variational characterization. ∎

######  Lemma 1.11.

Let $f$ be a function such that
$\operatorname{\mathbb{E}}^{\pi}[|f|^{2}]+\operatorname{\mathbb{E}}^{\pi^{\prime}}[|f|^{2}]=:f_{2}^{2}<\infty$,
then

$\big{|}\operatorname{\mathbb{E}}^{\pi}[f]-\operatorname{\mathbb{E}}^{\pi^{\prime}}[f]\big{|}\leq
2f_{2}d_{\mbox{\tiny{\rm H}}}(\pi,\pi^{\prime}).$

###### Proof.

$\displaystyle\Big{|}\operatorname{\mathbb{E}}^{\pi}[f]-\operatorname{\mathbb{E}}^{\pi^{\prime}}[f]\Big{|}=$
$\displaystyle\
\biggl{|}\int_{\mathbb{R}^{d}}f(u)\Bigl{(}\sqrt{\pi(u)}-\sqrt{\pi^{\prime}(u)}\Bigr{)}\Bigl{(}\sqrt{\pi(u)}+\sqrt{\pi^{\prime}(u)}\,\Bigr{)}\,du\biggr{|}$
$\displaystyle\leq$ $\displaystyle\
\biggl{(}\frac{1}{2}\int\Big{|}\sqrt{\pi(u)}-\sqrt{\pi^{\prime}(u)}\Big{|}^{2}du\biggr{)}^{1/2}\biggl{(}2\int|f(u)|^{2}\Big{|}\sqrt{\pi(u)}+\sqrt{\pi^{\prime}(u)}\Big{|}^{2}\,du\biggr{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle\ d_{\mbox{\tiny{\rm
H}}}(\pi,\pi^{\prime})\Big{(}4\int|f(u)|^{2}\bigl{(}\pi(u)+\pi^{\prime}(u)\bigr{)}du\Big{)}^{1/2}$
$\displaystyle=$ $\displaystyle\ 2f_{2}\,d_{\mbox{\tiny{\rm
H}}}(\pi,\pi^{\prime}).$

∎

Note that the result for Hellinger only assumes that $f$ is square integrable
with respect to $\pi$ and $\pi^{\prime}$. In contrast, the result for total
variation distance assumes that $f$ is bounded, which is a stronger condition.

#####  1.3.2 Approximation Theorem

We denote by

$g(u)=\nu\bigl{(}y-G(u)\bigr{)}\quad\text{and}\quad
g_{\delta}(u)=\nu\bigl{(}y-G_{\delta}(u)\bigr{)}$

the likelihoods associated with $G(u)$ and $G_{\delta}(u),$ so that

$\pi^{y}(u)=\frac{1}{Z}g(u)\rho(u)\quad\text{and}\quad\pi_{\delta}^{y}(u)=\frac{1}{Z_{\delta}}g_{\delta}(u)\rho(u),$

where $Z,Z_{\delta}>0$ are the corresponding normalizing constants. Before we
proceed to our main result we make some assumptions:

######  Assumption 1.12.

There exist $\delta^{+}>0$ and $K_{1},K_{2}<\infty$ such that, for all
$\delta\in(0,\delta^{+}),$

  * (i)

$|\sqrt{g(u)}-\sqrt{g_{\delta}(u)}|\leq\varphi(u)\delta$, for some
$\varphi(u)$ such that $\operatorname{\mathbb{E}}^{\rho}[\varphi^{2}(u)]\leq
K_{1}$;

  * (ii)

$\sup_{u\in\mathbb{R}^{d}}(|\sqrt{g(u)}|+|\sqrt{g_{\delta}(u)}|)\leq K_{2}.$

######  Remark 1.13.

Assumption 1.12 only involves conditions on the likelihood $g$ and the
approximate likelihood $g_{\delta}.$ While our presentation in this chapter
emphasizes that this approximation may arise due to the need of approximating
the forward model $G,$ another important scenario that the theory covers is
approximation due to perturbations of the data $y.$ Well-posedness results
that guarantee the stability under data perturbations of Bayesian data
assimilation will be established in Chapter 7.

Now we state the main result of this section:

######  Theorem 1.14 (Well-posedness of Posterior).

Under Assumption 1.12 we have

$d_{\mbox{\tiny{\rm H}}}(\pi^{y},\pi_{\delta}^{y})\leq
c\delta,\quad\delta\in(0,\tilde{\delta}^{+}),$

for some $\tilde{\delta}^{+}>0$ and some $c\in(0,+\infty)$ independent of
$\delta$.

Notice that this theorem together with Lemma 1.11 guarantee that expectations
computed with respect to $\pi^{y}$ and $\pi_{\delta}^{y}$ are order $\delta$
apart. To prove Theorem 1.14, we first show a lemma which characterizes the
normalization factor $Z_{\delta}$ in the small $\delta$ limit.

######  Lemma 1.15.

Under Assumption 1.12 there exist $\tilde{\delta}^{+}>0$,
$c_{1},c_{2}\in(0,+\infty)$ such that

$|Z-Z_{\delta}|\leq c_{1}\delta\quad\text{and}\quad
Z,Z_{\delta}>c_{2},\quad\text{for}\ \delta\in(0,\tilde{\delta}^{+}).$

###### Proof.

Since $Z=\int g(u)\rho(u)du$ and $Z_{\delta}=\int g_{\delta}(u)\rho(u)du$ we
have

$\displaystyle|Z-Z_{\delta}|=$ $\displaystyle\
\int\bigl{(}g(u)-g_{\delta}(u)\bigr{)}\rho(u)du$ $\displaystyle\leq$
$\displaystyle\
\Big{(}\int\Bigl{|}\sqrt{g(u)}-\sqrt{g_{\delta}(u)}\Bigr{|}^{2}\rho(u)du\Big{)}^{1/2}\Big{(}\int\Bigl{|}\sqrt{g(u)}+\sqrt{g_{\delta}(u)}\Bigr{|}^{2}\rho(u)du\Big{)}^{1/2}$
$\displaystyle\leq$ $\displaystyle\
\Big{(}\int\delta^{2}\varphi(u)^{2}\rho(u)du\Big{)}^{1/2}\Big{(}\int
K_{2}^{2}\rho(u)du\Big{)}^{1/2}$ $\displaystyle\leq$ $\displaystyle\
\sqrt{K_{1}}K_{2}\delta,\quad\delta\in(0,\delta^{+}).$

And when
$\delta\leq\tilde{\delta}^{+}\coloneqq\min\\{\frac{Z}{2\sqrt{K_{1}}K_{2}},\delta^{+}\\}$,
we have

$Z_{\delta}\geq Z-|Z-Z_{\delta}|\geq\frac{1}{2}Z.$

The lemma follows by taking $c_{1}=\sqrt{K_{1}}K_{2}$ and
$c_{2}=\frac{1}{2}Z$.
∎

###### Proof of Theorem 1.14.

We break the distance into two error parts, one caused by the difference
between $Z$ and $Z_{\delta}$, the other caused by the difference between $g$
and $g_{\delta}$:

$\displaystyle d_{\mbox{\tiny{\rm H}}}(\pi^{y},\pi_{\delta}^{y})=$
$\displaystyle\
\frac{1}{\sqrt{2}}\Bigl{\|}\sqrt{\pi^{y}}-\sqrt{\pi_{\delta}^{y}}\Bigr{\|}_{L^{2}}$
$\displaystyle=$ $\displaystyle\
\frac{1}{\sqrt{2}}\Big{\|}\sqrt{\frac{g\rho}{Z}}-\sqrt{\frac{g\rho}{Z_{\delta}}}+\sqrt{\frac{g\rho}{Z_{\delta}}}-\sqrt{\frac{g_{\delta}\rho}{Z_{\delta}}}\Big{\|}_{L^{2}}$
$\displaystyle\leq$ $\displaystyle\
\frac{1}{\sqrt{2}}\Big{\|}\sqrt{\frac{g\rho}{Z}}-\sqrt{\frac{g\rho}{Z_{\delta}}}\Big{\|}_{L^{2}}+\frac{1}{\sqrt{2}}\Big{\|}\sqrt{\frac{g\rho}{Z_{\delta}}}-\sqrt{\frac{g_{\delta}\rho}{Z_{\delta}}}\Big{\|}_{L^{2}}.$

Using Lemma 1.15, for $\delta\in(0,\tilde{\delta}^{+})$, we have

$\displaystyle\Big{\|}\sqrt{\frac{g\rho}{Z}}-\sqrt{\frac{g\rho}{Z_{\delta}}}\Big{\|}_{L^{2}}=$
$\displaystyle\
\Big{|}\frac{1}{\sqrt{Z}}-\frac{1}{\sqrt{Z_{\delta}}}\Big{|}\Big{(}\int
g(u)\rho(u)du\Big{)}^{1/2}$ $\displaystyle=$ $\displaystyle\
\frac{|Z-Z_{\delta}|}{(\sqrt{Z}+\sqrt{Z_{\delta}})\sqrt{Z_{\delta}}}$
$\displaystyle\leq$ $\displaystyle\ \frac{c_{1}}{2c_{2}}\delta,$

and

$\Big{\|}\sqrt{\frac{g\rho}{Z_{\delta}}}-\sqrt{\frac{g\rho}{Z_{\delta}}}\Big{\|}_{L^{2}}=\frac{1}{\sqrt{Z_{\delta}}}\Big{(}\int\Bigl{|}\sqrt{g(u)}-\sqrt{g_{\delta}(u)}\Bigr{|}^{2}\rho(u)du\Big{)}^{1/2}\leq\sqrt{\frac{K_{1}}{c_{2}}}\delta.$

Therefore

$d_{\mbox{\tiny{\rm
H}}}(\pi^{y},\pi_{\delta}^{y})\leq\frac{1}{\sqrt{2}}\frac{c_{1}}{2c_{2}}\delta+\frac{1}{\sqrt{2}}\sqrt{\frac{K_{1}}{c_{2}}}\delta=c\delta,$

with
$c=\frac{1}{\sqrt{2}}\frac{c_{1}}{2c_{2}}+\frac{1}{\sqrt{2}}\sqrt{\frac{K_{1}}{c_{2}}}$
independent of $\delta$. ∎

#####  1.3.3 Example

Many inverse problems arise from differential equations with unknown input
parameters. Here we consider a simple but typical example where $G(u)$ comes
from the solution of an ODE, which needs to be solved numerically. Let $x(t)$
be the solution to the initial value problem

$\frac{dx}{dt}=F(x;u),\quad x(0)=0,$ (1.5)

where $F:\mathbb{R}^{k}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{k}$ is a
function such that $F(x;u)$ and the partial Jacobian $D_{x}F(x;u)$ are
uniformly bounded with respect to $(x,u)$, i.e.

$|F(x;u)|,|D_{x}F(x;u)|<F_{{\rm max}},\quad\text{ for all
}(x,u)\in\mathbb{R}^{k}\times\mathbb{R}^{d},$

for some constant $F_{{\rm max}}$, and thus $F(x;u)$ is Lipschitz in $x$ in
that

$|F(x;u)-F(x^{\prime};u)|\leq F_{{\rm max}}|x-x^{\prime}|,\quad\text{ for all
}x,x^{\prime}\in\mathbb{R}^{k}.$

Now consider the inverse problem setting

$y=G(u)+\eta,$

where

$G(u):=x(1)=x(t)|_{t=1},$

and $\eta\sim\mathcal{N}(0,\gamma^{2}I_{k})$. We assume that the exact mapping
$G(u)$ is replaced by some numerical approximation $G_{\delta}(u)$. In
particular, $G_{\delta}(u)$ is given by using the forward Euler method to
solve the ODE (1.5). Define $X_{0}=0$, and

$X_{\ell+1}=X_{\ell}+\delta F(X_{\ell};u),\quad\ell\geq 0,$

where $\delta=\frac{1}{L}$ for some large integer $L$. Finally define
$G_{\delta}(u):=X_{L}$.

In what follows, we will prove that $G_{\delta}(u)$ is uniformly bounded and
close to $G(u)$ when $\delta$ is small, and then we will use these results to
show that Assumption 1.12 is satisfied. Therefore, we can apply Theorem 1.14
to this example and claim that the approximate posterior $\pi_{\delta}^{y}$ is
close to the unperturbed one $\pi^{y}$.

Define $t_{\ell}=\ell\delta$, $x_{\ell}=x(t_{\ell})$. The following lemma
gives an estimate on the error generated from using the forward Euler method.

######  Lemma 1.16.

Let $E_{\ell}:=x_{\ell}-X_{\ell}$. Then there is $c<\infty$ independent of
$\delta$ such that

$|E_{\ell}|\leq c\delta,\quad 0\leq\ell\leq L.$

In particular

$|G(u)-G_{\delta}(u)|=|E_{L}|\leq c\delta.$

###### Proof.

For simplicity of exposition we consider the case $k=1$; the case $k>1$ is
almost identical, simply requiring the integral form for the remainder term in
the Taylor expansion. Using Taylor expansion in the case $k=1,$ there is
$\xi_{\ell}\in[t_{\ell},t_{\ell+1}]$ such that

$\displaystyle x_{\ell+1}=$ $\displaystyle\
x_{\ell}+\delta\frac{dx}{dt}(t_{\ell})+\frac{\delta^{2}}{2}\frac{d^{2}x}{dt^{2}}(\xi_{\ell})$
$\displaystyle=$ $\displaystyle\ x_{\ell}+\delta
F(x_{\ell};u)+\frac{\delta^{2}}{2}D_{x}F\bigl{(}x(\xi_{\ell});u\bigr{)}F\bigl{(}x(\xi_{\ell});u\bigr{)}.$

Thus we have

$\displaystyle|E_{\ell+1}|=$ $\displaystyle\ |x_{\ell+1}-X_{\ell+1}|$
$\displaystyle=$ $\displaystyle\
\Big{|}x_{\ell}-X_{\ell}+\delta\Bigl{(}F(x_{\ell};u)-F(X_{\ell};u)\Bigr{)}+\frac{\delta^{2}}{2}D_{x}F\bigl{(}x(\xi_{\ell});u\bigr{)}F\bigl{(}x(\xi_{\ell});u\bigr{)}\Big{|}$
$\displaystyle\leq$ $\displaystyle\
|x_{\ell}-X_{\ell}|+\delta\bigl{|}F(x_{\ell};u)-F(X_{\ell};u)\bigr{|}+\frac{\delta^{2}}{2}\bigl{|}D_{x}F\bigl{(}x(\xi_{\ell});u\bigr{)}\bigr{|}\bigl{|}F\bigl{(}x(\xi_{\ell});u\bigr{)}\bigr{|}$
$\displaystyle\leq$ $\displaystyle\ |E_{\ell}|+\delta F_{{\rm
max}}|E_{\ell}|+\frac{\delta^{2}}{2}F_{{\rm max}}^{2}.$

Noticing that $|E_{0}|=0$, the discrete Gronwall inequality gives

$\displaystyle|E_{\ell}|\leq$ $\displaystyle\ (1+\delta F_{{\rm
max}})^{\ell}|E_{0}|+\frac{(1+\delta F_{{\rm max}})^{\ell}-1}{\delta F_{{\rm
max}}}\cdot\frac{\delta^{2}}{2}F_{{\rm max}}^{2}$ $\displaystyle\leq$
$\displaystyle\ \bigg{(}\Big{(}1+\frac{F_{{\rm
max}}}{L}\Big{)}^{L}-1\bigg{)}\cdot\frac{F_{{\rm max}}\delta}{2}$
$\displaystyle\leq$ $\displaystyle\ \frac{(e^{F_{{\rm max}}}-1)F_{{\rm
max}}}{2}\delta.$

The lemma follows by taking $c=\frac{(e^{F_{{\rm max}}}-1)F_{{\rm max}}}{2}$.
∎

######  Lemma 1.17.

For any $u\in\mathbb{R}^{d},$

$|G(u)|,|G_{\delta}(u)|<F_{{\rm max}}.$

###### Proof.

For $G(u)$ we use that $F(x;u)$ is uniformly bounded, so that

$|G(u)|=|x(1)|=\Big{|}\int_{0}^{1}F(x(t);u)dt\Big{|}\leq\int_{0}^{1}\bigl{|}F(x(t);u)\bigr{|}dt\leq
F_{{\rm max}}.$

As for $G_{\delta}(u)$, we first notice that

$|X_{\ell+1}|=|X_{\ell}+\delta
F(X_{\ell};u)|\leq|X_{\ell}|+\delta|F(X_{\ell};u)|\leq|X_{\ell}|+\delta
F_{{\rm max}},$

and by induction

$|X_{\ell}|\leq|X_{0}|+\ell\delta F_{{\rm max}}=\ell\delta F_{{\rm max}}.$

In particular,

$|G_{\delta}(u)|=|X_{L}|\leq L\delta F_{{\rm max}}=F_{{\rm max}}.$

∎

To conclude this chapter we show that in this example Assumption 1.12 is
satisfied. Recall that $\eta\sim\mathcal{N}(0,\gamma^{2}I_{k})$, and thus

$\sqrt{g(u)}=\sqrt{\nu\bigl{(}y-G(u)\bigr{)}}=\frac{1}{(2\pi)^{k/4}\gamma^{k/2}}\exp\Big{(}-\frac{1}{4\gamma^{2}}|y-G(u)|^{2}\Big{)},$
$\sqrt{g_{\delta}(u)}=\sqrt{\nu\bigl{(}y-G_{\delta}(u)\bigr{)}}=\frac{1}{(2\pi)^{k/4}\gamma^{k/2}}\exp\Big{(}-\frac{1}{4\gamma^{2}}|y-G_{\delta}(u)|^{2}\Big{)}.$

  * •

For Assumption 1.12(i) notice that the function $e^{-w}$ is Lipschitz for
$w>0$, with Lipschitz constant $1$. Therefore we have

$\displaystyle\Big{|}\sqrt{g(u)}-\sqrt{g_{\delta}(u)}\Big{|}\leq$
$\displaystyle\
\frac{1}{(2\pi)^{k/4}\gamma^{k/2}}\cdot\frac{1}{4\gamma^{2}}\cdot\big{|}|y-G(u)|^{2}-|y-G_{\delta}(u)|^{2}\big{|}$
$\displaystyle=$ $\displaystyle\
\frac{1}{(2\pi)^{k/4}\gamma^{k/2}}\cdot\frac{1}{4\gamma^{2}}\cdot|2y-G(u)-G_{\delta}(u)||G(u)-G_{\delta}(u)|$
$\displaystyle\leq$ $\displaystyle\
\frac{1}{(2\pi)^{k/4}\gamma^{k/2}}\cdot\frac{1}{4\gamma^{2}}\cdot(2|y|+2F_{{\rm
max}})c\delta$ $\displaystyle=$ $\displaystyle\ \tilde{c}\delta.$

That is to say, Assumption 1.12(i) is satisfied with $\varphi(u)=\tilde{c}$
and $\int_{\mathbb{R}^{d}}\varphi^{2}(u)\rho(u)du=\tilde{c}^{2}<\infty$.

  * •

Assumption 1.12(ii) is satisfied since

$\sqrt{g(u)}=\frac{1}{(2\pi)^{k/4}\gamma^{k/2}}\exp\Big{(}-\frac{1}{4\gamma^{2}}|y-G(u)|^{2}\Big{)}\leq\frac{1}{(2\pi)^{k/4}\gamma^{k/2}},$
$\sqrt{g_{\delta}(u)}=\frac{1}{(2\pi)^{k/4}\gamma^{k/2}}\exp\Big{(}-\frac{1}{4\gamma^{2}}|y-G_{\delta}(u)|^{2}\Big{)}\leq\frac{1}{(2\pi)^{k/4}\gamma^{k/2}}.$

####  1.4 Discussion and Bibliography

The book by Kaipio and Somersalo [63] provides an introduction to the Bayesian
approach to inverse problems, especially in the context of differential
equations. An overview of the subject of Bayesian inverse problems in
differential equations, with a perspective informed by the geophysical
sciences, is the book by Tarantola [108] (see, especially, Chapter 5).

In the paper [105] the Bayesian approach to regularization is reviewed,
developing a function space viewpoint on the subject. A well-posedness theory
and some algorithmic approaches which are used when adopting the Bayesian
approach to inverse problems are introduced. The function space viewpoint on
the subject is developed in more detail in the chapter notes of Dashti and
Stuart [23]. An application of this function space methodology to a large-
scale geophysical inverse problem is considered in [83]. The paper [76]
demonstrates the potential for the use of dimension reduction techniques from
control theory within statistical inverse problems.

See [46] for more detail on the subject of metrics, and other distance-like
functions, on probability measures. See [105, 23] for more detailed
discussions on the well-posedness of Bayesian inverse problems with respect to
perturbations in the data; and see [20] for applications concerning numerical
approximation of partial differential equations appearing in the forward
model. Related results, but using divergences rather than the Hellinger
metric, may be found in [84]. The paper [59] contains an interesting set of
examples where the Meta Theorem stated in this chapter fails in the sense
that, whilst well-posedness holds, the posterior is Hölder with exponent less
than one, rather than Lipschitz, with respect to perturbations.

###  2 The Linear-Gaussian Setting

Recall the inverse problem of estimating an unknown parameter
${u}\in\mathbb{R}^{d}$ from data $y\in\mathbb{R}^{k}$ under the model
assumption

$y=G(u)+\eta.$ (2.1)

In this chapter we study the linear-Gaussian setting, where the forward map
$G(\cdot)$ is linear and both the prior on $u$ and the distribution of the
observation noise $\eta$ are Gaussian. This setting is highly amenable to
analysis and arises frequently in applications. Moreover, as we will see
throughout this book, many methods employed in nonlinear or non-Gaussian
settings build on ideas from the linear-Gaussian case by performing
linearization or invoking Gaussian approximations. Having established a
formula for the posterior, we investigate the effect that the choice of prior
has on our solution. We do this by quantifying the spread of the posterior
distribution in the small noise (approaching zero) limit. This provides
intuitive understanding concerning the impact of the prior for overdetermined,
determined, and underdetermined regimes, corresponding to $d<k,d=k,$ and
$d>k,$ respectively.

The following will be assumed throughout this chapter.

######  Assumption 2.1.

The relationship between unknown $u\in\mathbb{R}^{d}$ and data
$y\in\mathbb{R}^{k}$ defined by equation (2.1) holds. Moreover,

1

Linearity of the forward model: $G({u})={A}{u}$, for some
${A}\in\mathbb{R}^{k\times d}$;

2

Gaussian prior: ${u}\sim\rho({u})=\mathcal{N}(\widehat{m},\widehat{C})$, where
$\widehat{C}$ is positive definite;

3

Gaussian noise: ${\eta}\sim\nu({\eta})=\mathcal{N}({0},{\Gamma})$, where
${\Gamma}$ is positive definite;

4

$u$ and $\eta$ are independent.

####  2.1 Derivation of the Posterior Distribution

Under Assumption 2.1 the likelihood on ${y}$ given ${u}$ is a Gaussian,

$\displaystyle{y}|{u}\sim\mathcal{N}({A}{u},{\Gamma}).$ (2.2)

Therefore, using Bayes formula (1.2) we see that the posterior $\pi^{y}(u)$ is
given by

$\displaystyle\pi^{{y}}({u})$
$\displaystyle=\frac{1}{Z}\nu({y}-{A}{u})\rho({u})$
$\displaystyle=\frac{1}{Z}\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|{y}-{A}{u}|_{{\Gamma}}^{2}}\right)\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|{u}-\widehat{m}|_{\widehat{C}}^{2}}\right)$
$\displaystyle=\frac{1}{Z}\exp\mathopen{}\mathclose{{}\left(-\frac{1}{2}|{y}-{A}{u}|_{{\Gamma}}^{2}-\frac{1}{2}|{u}-\widehat{m}|_{\widehat{C}}^{2}}\right)$
$\displaystyle=\frac{1}{Z}\exp\bigl{(}-{\mathsf{J}}({u})\bigr{)},$

with

${\mathsf{J}}({u})=\frac{1}{2}|{y}-{A}{u}|_{{\Gamma}}^{2}+\frac{1}{2}|{u}-\widehat{m}|_{\widehat{C}}^{2}.$
(2.3)

Since the posterior distribution can be written as the exponential of a
quadratic in $u$ it follows that the posterior is Gaussian. Its mean and
covariance are given in the following result.

######  Theorem 2.2 (Posterior is Gaussian).

Under Assumption 2.1 the posterior distribution is Gaussian,

$\displaystyle{u}|{y}\sim\pi^{{y}}({u})=\mathcal{N}(m,C).$ (2.4)

The posterior mean $m$ and covariance $C$ are given by the following formulae:

$\displaystyle m$
$\displaystyle=({A}^{T}{\Gamma}^{-1}{A}+\widehat{C}^{-1})^{-1}({A}^{T}{\Gamma}^{-1}{y}+\widehat{C}^{-1}\widehat{m}),$
(2.5) $\displaystyle C$
$\displaystyle=({A}^{T}{\Gamma}^{-1}{A}+\widehat{C}^{-1})^{-1}.$ (2.6)

###### Proof.

Since $\pi^{{y}}({u})=\frac{1}{Z}\exp\bigl{(}-{\mathsf{J}}({u})\bigr{)}$ with
${\mathsf{J}}(u)$ given by (2.3), a quadratic function of ${u}$, it follows
that the posterior distribution $\pi^{{y}}({u})$ is Gaussian. Denoting the
mean and variance of $\pi^{{y}}({u})$ by $m$ and $C$, we can write
${\mathsf{J}}({u})$ in the following form

$\displaystyle{\mathsf{J}}({u})$ $\displaystyle=\frac{1}{2}|{u}-m|_{C}^{2}+q,$
(2.7)

where the term $q$ does not depend on ${u}$. Now matching the coefficients of
the quadratic and linear terms in equations (2.3) and (2.7), we get

$\displaystyle C^{-1}$
$\displaystyle={A}^{T}{\Gamma}^{-1}{A}+\widehat{C}^{-1},$ $\displaystyle
C^{-1}m$ $\displaystyle={A}^{T}{\Gamma}^{-1}{y}+\widehat{C}^{-1}\widehat{m}.$

Therefore equations (2.5) and (2.6) follow. ∎

Equation (2.7) shows that the posterior mean $m$ minimizes ${\mathsf{J}}(u)$
given by (2.3). This demonstrates that the posterior mean is found by
compromising between maximizing the likelihood (by making small the _loss_
term $\frac{1}{2}|{y}-{A}{u}|_{{\Gamma}}^{2}$) and minimizing deviations from
the prior mean (by making small the _regularization_ term
$\frac{1}{2}|{u}-\widehat{m}|_{\widehat{C}}^{2}$). The relative importance
given to both objectives is determined by the relative size of the prior
covariance $\widehat{C}$ and the noise covariance $\Gamma.$ An important
feature of the linear-Gaussian setting is that the posterior covariance $C$
does not depend on the data $y;$ this is not true in general.

As we saw in the previous chapter, the posterior mean estimator and the MAP
estimator are typically different. However, in the linear-Gaussian setting the
posterior is Gaussian, and therefore both estimators agree.

######  Corollary 2.3 (Characterization of Bayes Estimators).

The posterior mean and MAP estimators under Assumptions 2.1 agree, and are
given by $u_{\mbox{\tiny{\rm MAP}}}=u_{\mbox{\tiny{\rm PM}}}=m$ defined in
equation (2.5).

######  Example 2.4.

Let ${\Gamma}=\gamma^{2}{I}$, $\widehat{C}=\sigma^{2}{I},$ $\widehat{m}=0,$
and set $\lambda=\frac{\sigma^{2}}{\gamma^{2}}$. Then

${\mathsf{J}}_{\lambda}({u}):=\gamma^{2}{\mathsf{J}}({u})=\frac{1}{2}|{y}-{A}{u}|^{2}+\frac{\lambda}{2}|{u}|^{2}.$

Since $m$ minimizes ${\mathsf{J}}_{\lambda}(\cdot)$ it follows that

$\displaystyle({A}^{T}{A}+\lambda{I})m={A}^{T}{y}.$ (2.8)

Example 2.4 provides a link between Bayesian inversion and optimization
approaches to inversion: ${\mathsf{J}}_{\lambda}({u})$ can be seen as the
objective functional in a linear regression model with a regularizer
$\frac{\lambda}{2}|{u}|^{2}$, as used in ridge regression. Equation (2.8) for
$m$ is exactly the normal equation with regularizer in the least square
problem. In fact, in the general case, equation (2.5) can be also viewed as a
generalized normal equation. This point of view helps us understand the
structure of Bayesian regularization by linking it to the deep understanding
of optimization approaches to inverse problems. A more extensive account of
the optimization perspective and its interplay with Bayesian formulations will
be given in the following chapter.

####  2.2 Small Noise Limit of the Posterior Distribution

In this section we study the small observational noise limit of the posterior
distribution in the linear-Gaussian setting. While most of the ideas can be
extended beyond this setting, the explicit calculations that the linear-
Gaussian setting allows for provide helpful intuition. Throughout this section
we assume the following.

######  Assumption 2.5.

In addition to the linear-Gaussian setting described in Assumption 2.1,
${\eta}:=\gamma{\eta}_{0}$ where ${\eta}_{0}\sim\mathcal{N}({0},{\Gamma}_{0})$
so that $\Gamma=\gamma^{2}\Gamma_{0}.$

Note that substituting $\Gamma=\gamma^{2}\Gamma_{0}$ into (2.5) and (2.6) we
obtain that

$\displaystyle m=$
$\displaystyle(A^{T}\Gamma_{0}^{-1}A+\gamma^{2}\widehat{C}^{-1})^{-1}(A^{T}\Gamma_{0}^{-1}y+\gamma^{2}\widehat{C}^{-1}\widehat{m}),$
(2.9) $\displaystyle C=$
$\displaystyle\gamma^{2}(A^{T}\Gamma_{0}^{-1}A+\gamma^{2}\widehat{C}^{-1})^{-1}.$
(2.10)

In the next three subsections we study the small noise limiting behavior as
$\gamma\to 0$ of the posterior mean and covariance in the overdetermined,
determined, and underdetermined regimes. We recall that $\Rightarrow$ denotes
weak convergence of probability measures. We will use repeatedly that weak
convergence of Gaussian distributions is equivalent to the convergence of
their means and covariances. In particular, the weak limit of a sequence of
Gaussians with means converging to $m^{+}$ and covariance matrices converging
to zero is a Dirac mass $\delta_{m^{+}}.$

#####  2.2.1 Overdetermined Case

We start with the overdetermined case $d<k$.

######  Theorem 2.6 (Small Noise Limit of Posterior Distribution -
Overdetermined).

Suppose that Assumption 2.5 holds, that ${\rm Null}(A)={0}$ and that $d<k.$
Then in the limit $\gamma^{2}\to 0,$

$\pi^{y}\Rightarrow\delta_{m^{+}},$

where $m^{+}$ is the solution of the least-squares problem

$m^{+}=\arg\min_{u\in\mathbb{R}^{d}}|\Gamma_{0}^{-1/2}(y-Au)|^{2}.$ (2.11)

###### Proof.

Since Null$(A)={0}$ and $\Gamma_{0}$ is invertible we deduce that there is
$\alpha>0$ such that, for all $u\in\mathbb{R}^{d},$

$\langle
u,A^{T}\Gamma^{-1}_{0}Au\rangle=|\Gamma_{0}^{-1/2}Au|^{2}\geq\alpha|u|^{2}.$

Thus $A^{T}\Gamma_{0}^{-1}A$ is positive definite and invertible. It follows
that as $\gamma\to 0$, the posterior covariance converges to the zero matrix,
$C\to 0,$ and the posterior mean satisfies the limit

$m\to m^{*}=(A^{T}\Gamma_{0}^{-1}A)^{-1}A^{T}\Gamma_{0}^{-1}y.$

This proves the weak convergence of $\pi^{y}$ to $\delta_{m^{*}}$. It remains
to characterize $m^{*}$. Since Null$(A)={0}$, the minimizers of

$\mathsf{L}(u):=\frac{1}{2}|\Gamma_{0}^{-1/2}(y-Au)|^{2}$

are unique and satisfy the normal equations
$A^{T}\Gamma_{0}^{-1}Au=A^{T}\Gamma_{0}^{-1}y$. Hence $m^{*}$ solves the
desired least-squares problem and coincides with $m^{+}$ given in (2.11). ∎

######  Remark 2.7.

In the overdetermined case where $A^{T}\Gamma_{0}A$ is invertible, the small
observational noise limit leads to a posterior which is a Dirac, centered at
the solution of the least-square problem (2.11). Therefore, in this limit the
prior plays no role in the Bayesian inference.

######  Theorem 2.8 (Posterior Consistency).

Suppose that the assumptions of Theorem 2.18 hold and that the data satisfies

$y=Au^{\dagger}+\gamma\eta_{0}^{\dagger},\quad\quad{\rm for\,\,fixed}\quad
u^{\dagger}\in\mathbb{R}^{d},\eta_{0}^{\dagger}\in\mathbb{R}^{k}.$ (2.12)

Then, for any sequence $M(\gamma)\rightarrow\infty$ as $\gamma\rightarrow 0$,

$\displaystyle\operatorname{\mathbb{P}}^{\pi^{y}}\bigl{\\{}|{u}-{u}^{\dagger}|^{2}>M(\gamma)\gamma^{2}\bigr{\\}}\rightarrow
0,$ (2.13)

where $\operatorname{\mathbb{P}}^{\pi^{y}}$ denotes probability under the
posterior distribution.

######  Remark 2.9.

For any $\varepsilon>0$, set $M(\gamma)=\frac{\varepsilon^{2}}{\gamma^{2}}$ in
Theorem 2.11 to obtain

$\operatorname{\mathbb{P}}^{\pi_{y}}\bigl{\\{}|{u}-{u}^{\dagger}|>\varepsilon\bigr{\\}}\rightarrow
0.$

###### Proof.

(Theorem 2.6) Throughout this proof we let $c$ be a constant independent of
$\gamma$ that may change from line to line, and we denote by
$\operatorname{\mathbb{E}}$ expectation with respect to the posterior
distribution, which is Gaussian with mean $m$ and covariance $C$ given by
equations (2.9) and (2.10). Denote

$m^{*}=(A^{T}\Gamma_{0}^{-1}A)^{-1}A^{T}\Gamma_{0}^{-1}y$

as in the proof of the previous theorem. We have that

$\operatorname{\mathbb{E}}|u-u^{\dagger}|^{2}\leq
c\Bigl{(}\operatorname{\mathbb{E}}|u-m|^{2}+|m-m^{*}|^{2}+|m^{*}-u^{\dagger}|^{2}\Bigr{)}.$
(2.14)

We now bound each of the three terms in the right-hand side.

For the first one,

$\displaystyle\operatorname{\mathbb{E}}|u-m|^{2}$
$\displaystyle=\operatorname{\mathbb{E}}\bigl{[}(u-m)^{T}(u-m)\bigr{]}=\operatorname{\mathbb{E}}\bigl{[}{\rm
Tr}[(u-m)\otimes(u-m)]\bigr{]}$ $\displaystyle={\rm
Tr}\operatorname{\mathbb{E}}\bigl{[}(u-m)\otimes(u-m)\bigr{]}$
$\displaystyle={\rm Tr}(C)\leq\gamma^{2}{\rm Tr}(A^{T}\Gamma_{0}^{-1}A).$

For the second term, note that

$\displaystyle(A^{T}\Gamma_{0}^{-1}A)m^{*}$
$\displaystyle=A^{T}\Gamma_{0}^{-1}y,$
$\displaystyle(A^{T}\Gamma_{0}^{-1}A+\gamma^{2}\widehat{C}^{-1})m$
$\displaystyle=A^{T}\Gamma_{0}^{-1}y+\gamma^{2}\widehat{C}^{-1}\widehat{m}.$

Therefore

$m-m^{*}=\gamma^{2}(A^{T}\Gamma_{0}A)^{-1}(\widehat{C}^{-1}\widehat{m}-\widehat{C}^{-1}m).$

Since $m$ converges it is bounded, and so there is $c>0$ such that

$|m-m^{*}|^{2}\leq c\gamma^{4}.$

Finally, for the third term we write

$\displaystyle m^{*}$
$\displaystyle=(A^{T}\Gamma_{0}^{-1}A)^{-1}A^{T}\Gamma_{0}^{-1}Au^{\dagger}+\gamma(A^{T}\Gamma_{0}^{-1}A)^{-1}A^{T}\Gamma_{0}^{-1}\eta_{0}^{\dagger}$
$\displaystyle=u^{\dagger}+\gamma(A^{T}\Gamma_{0}^{-1}A)^{-1}A^{T}\Gamma_{0}^{-1}\eta_{0}^{\dagger},$

which gives

$|m^{*}-u^{\dagger}|^{2}\leq c\gamma^{2}.$

Using Markov’s inequality and the three bounds above,

$\operatorname{\mathbb{P}}^{\pi^{y}}\bigl{\\{}|{u}-{u}^{\dagger}|^{2}>M(\gamma)\gamma^{2}\bigr{\\}}\leq\frac{\operatorname{\mathbb{E}}(|{u}-{u}^{\dagger}|^{2})}{M(\gamma)\gamma^{2}}\leq\frac{c}{M(\gamma)}\rightarrow
0,~{}{\rm as}~{}\gamma\rightarrow 0.$

∎

#####  2.2.2 Determined Case

As a byproduct of the proof of Theorem 2.6, we can determine the limiting
behavior of $\pi^{y}$ in the boundary case $d=k$.

######  Theorem 2.10 (Small Noise Limit of Posterior Distribution -
Determined).

Suppose that Assumption 2.5 holds, ${\rm Null}(A)=0,$ and $d=k$. Then in the
small noise limit $\gamma^{2}\to 0$,

$\pi^{y}\Rightarrow\delta_{A^{-1}y}.$

###### Proof.

In the proof of Theorem 2.6, the assumption $d<k$ is used only in that $A$ is
not a square matrix and thus $A,A^{T}$ are not invertible. Denote by $(m,C)$
the mean and variance of the posterior $u|y$. Using the same argument, we have
$C\to 0$ and

$m\to m^{*}=(A^{T}\Gamma_{0}^{-1}A)^{-1}A^{T}\Gamma_{0}^{-1}y.$

Using that $A,A^{T}$ are square invertible matrices we obtain

$m^{*}=(A^{-1}\Gamma_{0}(A^{T})^{-1})A^{T}\Gamma_{0}^{-1}y=A^{-1}y.$

Therefore, $\pi^{y}(u)\Rightarrow\delta_{m^{*}}=\delta_{A^{-1}y}$. ∎

Note that here, as in the overdetermined case, the prior plays no role in the
small noise limit. Moreover, it can be shown as above that posterior
consistency holds. The proof is identical and therefore omitted.

######  Theorem 2.11 (Posterior Consistency).

Suppose that the assumptions of Theorem 2.10 hold, and that the data satisfies

$y=Au^{\dagger}+\gamma\eta_{0}^{\dagger},\quad\quad{\rm for\,\,fixed}\quad
u^{\dagger},\eta_{0}^{\dagger}\in\mathbb{R}^{d}.$ (2.15)

Then for any sequence $M(\gamma)\rightarrow\infty$ as $\gamma\rightarrow 0$,

$\displaystyle\operatorname{\mathbb{P}}^{\pi^{y}}\bigl{\\{}|{u}-{u}^{\dagger}|^{2}>M(\gamma)\gamma^{2}\bigr{\\}}\rightarrow
0.$ (2.16)

#####  2.2.3 Underdetermined Case

Finally we consider the underdetermined case $d>k$. We assume that
$A\in\mathbb{R}^{k\times d}$ with ${\rm Rank}(A)=k$ and write

$A=(A_{0}\ 0)Q^{T}=(A_{0}\ 0)(Q_{1}\ Q_{2})^{T}=A_{0}Q_{1}^{T},$ (2.17)

with $A_{0}\in\mathbb{R}^{k\times k}$ an invertible matrix, $Q=(Q_{1}\
Q_{2})\in\mathbb{R}^{d\times d}$ an orthogonal matrix so that $Q^{T}Q=I$,
$Q_{1}\in\mathbb{R}^{d\times k},Q_{2}\in\mathbb{R}^{d\times(d-k)}$. We have
the following result:

######  Theorem 2.12 (Small Noise Limit of Posterior Distribution -
Underdetermined).

Suppose that Assumption 2.5 holds, that ${\rm Rank}(A)=k,$ and $d>k$. In the
small noise limit $\gamma^{2}\to 0$,

$\pi^{y}\Rightarrow\mathcal{N}(m^{+},C^{+}),$

where

$\displaystyle m^{+}$
$\displaystyle=\widehat{C}Q_{1}(Q_{1}^{T}\widehat{C}Q_{1})^{-1}A_{0}^{-1}y+Q_{2}(Q_{2}^{T}\widehat{C}^{-1}Q_{2})^{-1}Q_{2}^{T}\widehat{C}^{-1}\widehat{m},$
$\displaystyle C^{+}$
$\displaystyle=Q_{2}(Q_{2}^{T}\widehat{C}^{-1}Q_{2})^{-1}Q_{2}^{T}.$

Since ${\rm Rank}(C^{+})={\rm Rank}(Q_{2})=d-k<d$ this theorem demonstrates
that, in the small observational noise limit, the posterior retains
uncertainty in a subspace of dimension $d-k$, and has no uncertainty in a
subspace of dimension $k$. As a consequence there is no posterior consistency
in the underdetermined case.

######  Example 2.13.

To help understand the result in Theorem 2.12, we consider a simple explicit
example. Assume that $A=(A_{0}\ 0)\in\mathbb{R}^{k\times
d},\Gamma=\gamma^{2}\Gamma_{0}=\gamma^{2}I_{k},\widehat{C}=I_{d}$,
$\widehat{m}=0$. Let
$u=\mathopen{}\mathclose{{}\left(\begin{array}[]{c}u_{1}\\\ u_{2}\\\
\end{array}}\right)\sim\mathcal{N}(0,I_{d}),u_{1}\in\mathbb{R}^{k},u_{2}\in\mathbb{R}^{d-k}$.
The data then satisfies

$y=Au+\eta=A_{0}u_{1}+\eta,\ \eta\sim\mathcal{N}(0,\gamma^{2}I_{k}).$

The posterior $u|y$ is
$\pi_{\gamma}^{y}(u)=\frac{1}{Z_{\gamma}}\exp(-{\mathsf{J}}_{\gamma}(u))$,
where

$\displaystyle{\mathsf{J}}_{\gamma}(u)$
$\displaystyle=\frac{1}{2\gamma^{2}}|y-A_{0}u_{1}|^{2}+\frac{1}{2}|u|^{2}$
$\displaystyle=\mathopen{}\mathclose{{}\left(\frac{1}{2\gamma^{2}}|y-A_{0}u_{1}|^{2}+\frac{1}{2}|u_{1}|^{2}}\right)+\frac{1}{2}|u_{2}|^{2}.$
(2.18)

It is clear that

$\pi_{\gamma}^{y}(u_{1})\Rightarrow\delta_{A_{0}^{-1}y}(u_{1}).$

Once $u_{1}$ is fixed as $A_{0}^{-1}y$, the first term in (2.18) is a constant
$\frac{1}{2}|A_{0}^{-1}y|^{2}$. Since $u_{1}$ and $u_{2}$ are independent we
can derive, formally, the limiting posterior as follows

$\pi_{\gamma}^{y}(u)\Rightarrow\delta_{A_{0}^{-1}y}(u_{1})\otimes\frac{1}{Z}\exp(-\frac{1}{2}|u_{2}|^{2})=\delta_{A_{0}^{-1}y}(u_{1})\otimes\mathcal{N}(0,I_{d-k}),$

where $Z=\int_{\mathbb{R}^{d-k}}\exp(-\frac{1}{2}|u_{2}|^{2})du_{2}$. In fact,
this is exactly the limiting posterior measure given in Theorem 2.12.

To prove Theorem 2.12, we use the following decomposition of the identity
$I_{d}.$

######  Lemma 2.14.

Let $\widehat{C}\in\mathbb{R}^{d\times d}$ be invertible and $Q=(Q_{1}\
Q_{2})$ be an orthonormal matrix with $Q_{1}\in\mathbb{R}^{d\times
k},Q_{2}\in\mathbb{R}^{d\times(d-k)}$. We have the following decomposition of
$I_{d}$

$I_{d}=\widehat{C}Q_{1}(Q_{1}^{T}\widehat{C}Q_{1})^{-1}Q_{1}^{T}+Q_{2}(Q_{2}^{T}\widehat{C}^{-1}Q_{2})^{-1}Q_{2}^{T}\widehat{C}^{-1}.$
(2.19)

###### Proof.

Denote by $R$ the right-hand side of (2.19). Since $Q$ is orthonormal, we have
$Q_{1}^{T}Q_{2}=0,Q_{2}^{T}Q_{1}=0$ and thus

$Q_{1}^{T}(R-I)=0,\quad Q^{T}_{2}\widehat{C}^{-1}(R-I)=0.$

If $B:=(Q_{1}\ \widehat{C}^{-1}Q_{2})$ is full rank, the above identities
imply that $B^{T}(R-I)=0$ and thus $R=I$. Note that

$Q^{T}B=\mathopen{}\mathclose{{}\left(\begin{array}[]{c}Q_{1}^{T}\\\
Q_{2}^{T}\\\ \end{array}}\right)(Q_{1}\
\widehat{C}^{-1}Q_{2})=\mathopen{}\mathclose{{}\left(\begin{array}[]{cc}I_{k}&Q_{1}^{T}\widehat{C}^{-1}Q_{2}\\\
0&Q_{2}^{T}\widehat{C}^{-1}Q_{2}\\\ \end{array}}\right).$

Since the last matrix is invertible, $B$ is invertible and the proof is
complete. ∎

###### Proof of Theorem 2.12.

Using (2.19), we can decompose $u$ as follows

$\displaystyle u$
$\displaystyle=\underbrace{\widehat{C}Q_{1}(Q_{1}^{T}\widehat{C}Q_{1})^{-1}}_{S}\underbrace{Q_{1}^{T}u}_{u_{1}}+\underbrace{Q_{2}(Q_{2}^{T}\widehat{C}^{-1}Q_{2})^{-1}}_{T}\underbrace{Q_{2}^{T}\widehat{C}^{-1}u}_{u_{2}}=Su_{1}+Tu_{2}.$

Here $u_{1}$ and $u_{2}$ are Gaussian with
$u_{2}\sim\mathcal{N}(Q_{2}^{T}\widehat{C}^{-1}\widehat{m},Q_{2}^{T}\widehat{C}^{-1}Q_{2})$.
The identity

${\rm Cov}(u_{1},u_{2})=Q_{1}^{T}{\rm
Cov}(u,u)\widehat{C}^{-1}Q_{2}=Q_{1}^{T}Q_{2}=0$

shows that $u_{1}$ and $u_{2}$ are independent, written $u_{1}\perp u_{2}$.
From (2.17), we have

$y=Au+\eta=A_{0}Q_{1}^{T}u+\eta=A_{0}u_{1}+\eta.$ (2.20)

Since $u\perp\eta$ and $u_{1}\perp u_{2}$, we have that $u_{2}\perp y,u_{1}$.
We apply conditional probability to yield

$\pi^{y}(u_{1},u_{2}):=\operatorname{\mathbb{P}}(u_{1},u_{2}|y)=\operatorname{\mathbb{P}}(u_{2})\operatorname{\mathbb{P}}(u_{1}|y).$

Equation (2.20) suggests that $(y,u_{1},\eta)$ with
$A_{0}\in\mathbb{R}^{k\times k}$ invertible is Gaussian distributed and its
posterior is exactly $\operatorname{\mathbb{P}}(u_{1}|y)$. Theorem 2.10 shows
that
$\operatorname{\mathbb{P}}(u_{1}|y)\Rightarrow\delta_{A_{0}^{-1}y}(u_{1})$ as
the noise vanishes, that is, as $\gamma^{2}\to 0$. Note that $u_{2}\perp
u_{1}$ and $u_{2}\perp y$. The limiting posterior measure $(u_{1},u_{2})|y$ is

$\pi^{y}(u_{1},u_{2})\Rightarrow\operatorname{\mathbb{P}}(u_{2})\otimes\delta_{A_{0}^{-1}y}(u_{1})$
(2.21)

as $\gamma^{2}\to 0$. Recall $u=Su_{1}+Tu_{2}$ and
$u_{2}\sim\mathcal{N}(Q_{2}^{T}\widehat{C}^{-1}\widehat{m},Q_{2}^{T}\widehat{C}^{-1}Q_{2})$.
The mean and variance of the limiting posterior measure $u|y$ is

$\displaystyle m^{+}$
$\displaystyle=\operatorname{\mathbb{E}}(Su_{1}+Tu_{2}|y)=SA_{0}^{-1}y+T\operatorname{\mathbb{E}}(u_{2})=SA_{0}^{-1}y+TQ_{2}^{T}\widehat{C}^{-1}\widehat{m},$
$\displaystyle C^{+}$ $\displaystyle={\rm Var}(Su_{1}+Tu_{2}|y)={\rm
Var}(Tu_{2})=TQ_{2}^{T}\widehat{C}^{-1}Q_{2}T^{T}=Q_{2}(Q_{2}^{T}\widehat{C}^{-1}Q_{2})^{-1}Q_{2}^{T}.$

We have thus completed the proof. ∎

######  Remark 2.15.

Equation (2.21) shows that in the limit of zero observational noise, the
uncertainty is only in the variable $u_{2}$. Since ${\rm Span}(T)={\rm
Span}(Q_{2})$ and $u=SA_{0}^{-1}y+Tu_{2}$, the uncertainty we observed is in
${\rm Span}(Q_{2})$. The prior plays a role in the posterior measure, in the
limit of zero observational noise, but only in the variables $u_{2}.$

####  2.3 Discussion and Bibliography

The linear Gaussian setting plays, for several reasons, a central role in the
study of inverse problems. One is that it allows explicit solutions which can
be used to give insight into the subject area more generally. The second is
that in the large data limit many Bayesian posteriors are approximately
Gaussian. The paper [39], which is in the linear-Gaussian setting, plays an
important role in the history of Bayesian inversion as it was arguably the
first to formulate Bayesian inversion in function space.

We have also employed the Gaussian setting to present a basic form of
posterior consistency. For a treatment in infinite dimensions see [71, 2, 89].
For the consistency problem in the classical statistical setting, see the
books [48, 27]. The book [27] also contains definition and properties of
convergence in probability as used here. For the non-statistical approach to
inverse problems, and consistency results, see [31] and the references
therein.

###  3 Optimization Perspective

In this chapter we explore the properties of Bayesian inversion from the
perspective of an optimization problem which corresponds to maximizing the
posterior probability, in a sense which we will make precise. We demonstrate
the properties of the point estimator resulting from this optimization
problem, showing its positive and negative attributes, the latter motivating
our work in the following chapter.

####  3.1 The Setting

Once again we work in the inverse problem setting of finding
$u\in\mathbb{R}^{d}$ from $y\in\mathbb{R}^{k}$ given by

$y=G(u)+\eta$

with noise $\eta\sim\nu(\cdot)$ and prior $u\sim\rho,$ as in Assumption 1.1.
The posterior $\pi^{y}(u)$ on $u|y$ is given by Theorem 1.2 and has the form

$\pi^{y}({u})=\frac{1}{Z}\nu({y}-{G}({u}))\rho({u}).$

We may define a loss:

$\mathsf{L}({u})=-\log\nu\bigl{(}{y}-{G}({u})\bigr{)},$

and a regularizer

$\mathsf{R}({u})=-\log\rho({u}).$

When added together these two functions of $u$ comprise an objective function
of the form

${\mathsf{J}}({u})=\mathsf{L}({u})+\mathsf{R}({u}).$

Furthermore

$\pi^{y}({u})=\frac{1}{Z}\nu\bigl{(}{y}-{G}({u})\bigr{)}\rho({u})\propto
e^{-{\mathsf{J}}({u})}.$

We see that minimizing the objective function ${\mathsf{J}}(\cdot)$ is
equivalent to maximizing the posterior $\pi^{y}(\cdot)$. Therefore, recalling
Definition 1.4, the MAP estimator can be rewritten in terms of ${\mathsf{J}}$
as follows:

$\displaystyle u_{\mbox{\tiny{\rm MAP}}}$
$\displaystyle=\arg\max_{{u}\in\mathbb{R}^{d}}\pi^{{y}}({u})$
$\displaystyle=\arg\min_{u\in\mathbb{R}^{d}}{\mathsf{J}}(u).$

We will provide conditions under which the MAP estimator is attained in
Theorem 3.5 and we will give an interpretation of MAP estimators in terms of
maximizing the probability of infinitesimal balls in Theorem 3.8. This
interpretation allows to generalize the definition of MAP estimators to
measures that do not possess a Lebesgue density.

######  Example 3.1.

Consider the Gaussian setting of Assumption 2.1. Then, since the posterior is
Gaussian, its mode agrees with its mean, which is given by $m$ as defined in
Theorem 2.2.

######  Example 3.2.

If $\eta=\mathcal{N}({0},{\Gamma})$, then
$\nu\bigl{(}{y}-{G}({u})\bigr{)}\propto\exp(-\frac{1}{2}|{y}-{G}({u})|_{{\Gamma}}^{2})$.
So the loss in this case is
$\mathsf{L}({u})=\frac{1}{2}|{y}-{G}({u})|_{{\Gamma}}^{2}$, a
${\Gamma}$-weighted $\mathsf{L}_{2}$ loss.

######  Example 3.3.

If we have prior $\rho({u})=\mathcal{N}({0},{\widehat{C}})$, then ignoring
$u-$independent normalization factors, which appear as constant shifts in
${\mathsf{J}}(\cdot)$, we may take the regularizer as
$\mathsf{R}({u})=\frac{1}{2}|{u}|_{{\widehat{C}}}^{2}$. In particular, if
${\widehat{C}}=\lambda^{-1}{I}$, then
$\mathsf{R}({u})=\frac{\lambda}{2}|{u}|^{2}$, an $\mathsf{L}_{2}$ regularizer.

If we combine Example 3.2 and Example 3.3, we obtain a canonical objective
function

${\mathsf{J}}({u})=\frac{1}{2}|{y}-{G}({u})|_{{\Gamma}}^{2}+\frac{\lambda}{2}|{u}|^{2}.$

To connect with future discussions, here $\lambda$ corresponds to prior
precision, and may be learned from data, such as in hierarchical methods.

######  Example 3.4.

As an alternative to the $\mathsf{L}_{2}$ regularizer, consider
${u}=(u_{1},\ldots,u_{d})$ with $u_{i}$ having prior distribution i.i.d.
Laplace. Then
$\rho({u})\propto\exp(-\lambda\sum_{i=1}^{d}|u_{i}|)=\exp(-\lambda|{u}|_{1})$.
In this case $\mathsf{R}({u})=\lambda|{u}|_{1}$, an $\mathsf{L}_{1}$
regularizer. If we combine this prior with the weighted $\mathsf{L}_{2}$ loss
above then we have objective function

${\mathsf{J}}({u})=\frac{1}{2}|{y}-{G}({u})|_{{\Gamma}}^{2}+\lambda|{u}|_{1}.$

Even though this objective function promotes sparse solutions, samples from
the underlying posterior distribution are not typically sparse.

####  3.2 Theory

For any optimization problem for an objective function with a finite infimum,
it is of interest to determine whether the infinimum is attained. We have the
following result.

######  Theorem 3.5 (Attainable MAP Estimator).

Assume that ${\mathsf{J}}$ is non-negative, continuous and that
${\mathsf{J}}(u)\to\infty$ as $|u|\to\infty.$ Then ${\mathsf{J}}$ attains it
infimum. Therefore, the MAP estimator of $u$ based on the posterior
$\pi^{y}(u)\propto\exp\bigl{(}-{\mathsf{J}}(u)\bigr{)}$ is attained.

###### Proof.

By the assumed growth and non-negativity of ${\mathsf{J}}$, there is $R$ such
that $\inf_{u\in\mathbb{R}^{d}}{\mathsf{J}}(u)=\inf_{u\in B(0,R)}$ where
$B(0,R)$ denotes the closed ball of radius $R$ around the origin. Since
${\mathsf{J}}$ is assumed to be continuous, its infimum over $B(0,R)$ is
attained and the proof is complete. ∎

######  Remark 3.6.

Suppose that

  1. 1.

${G}\in C(\mathbb{R}^{d},\mathbb{R}^{k})$, i.e. ${G}$ is a continuous
function;

  2. 2.

the objective function ${\mathsf{J}}({u})$ has $\mathsf{L}_{2}$ loss as
defined in Example 3.2 and $\mathsf{L}_{p}$ regularizer
$\mathsf{R}({u})=\frac{\lambda}{p}|{u}|_{p}^{p}$, $p\in(0,\infty)$.

Then the assumptions on ${\mathsf{J}}$ in Theorem 3.5 are satisfied. This
shows that if ${G}$ is continuous, the infinimum of ${\mathsf{J}}$ defined
with $\mathsf{L}_{2}$ loss and $\mathsf{L}_{p}$ regularizer is attained at the
MAP estimator of the corresponding Bayesian problem with posterior pdf
proportional to $\exp\bigl{(}-{\mathsf{J}}(u)\bigr{)}$.

######  Remark 3.7.

Notice that the assumption that ${\mathsf{J}}(u)\to\infty$ is not restrictive:
this condition needs to hold in order to be able to normalize
$\pi^{y}(u)\propto\exp\bigl{(}-{\mathsf{J}}(u)\bigr{)}$ into a pdf, which is
implicitly assumed in the second part of the theorem statement.

Intuitively the MAP estimator maximizes posterior probability. We make this
precise in the following theorem which links the objective function
${\mathsf{J}}(\cdot)$ to small ball probabilities.

######  Theorem 3.8 (Objective Function and Posterior Probability).

Making the same assumptions as in Theorem 3.5, let

$\alpha({u},\delta):=\int_{{v}\in
B({u},\delta)}\pi^{y}({v})d{v}=\operatorname{\mathbb{P}}^{\pi^{y}}\bigl{(}B({u},\delta)\bigr{)},$

be the posterior probability of a ball with radius $\delta$ centered at $u$.
Then, for all ${u}_{1},{u}_{2}\in\mathbb{R}^{d}$, we have

$\lim_{\delta\rightarrow
0}\frac{\alpha({u}_{1},\delta)}{\alpha({u}_{2},\delta)}=e^{{\mathsf{J}}({u}_{2})-{\mathsf{J}}({u}_{1})}.$

###### Proof.

Let $u_{1},u_{2}\in\mathbb{R}^{d}$ and let $\epsilon>0.$ By continuity of
${\mathsf{J}}$ we have that, for all $\delta$ sufficiently small,

$\displaystyle e^{-{\mathsf{J}}(u_{1})-\epsilon}\leq e^{-{\mathsf{J}}(v)}$
$\displaystyle\leq e^{-{\mathsf{J}}(u_{1})+\epsilon}\quad\text{for
all}\,\,v\in B(u_{1},\delta),$ $\displaystyle
e^{-{\mathsf{J}}(u_{2})-\epsilon}\leq e^{-{\mathsf{J}}(v)}$ $\displaystyle\leq
e^{-{\mathsf{J}}(u_{2})+\epsilon}\quad\text{for all}\,\,v\in B(u_{2},\delta).$

Therefore, for all $\delta$ sufficiently small,

$\displaystyle B_{\delta}e^{-{\mathsf{J}}({u}_{1})-\epsilon}\leq\int_{v\in
B(u_{1},\delta)}e^{-{\mathsf{J}}({v})}d{v}\leq
B_{\delta}e^{-{\mathsf{J}}({u}_{1})+\epsilon},$ $\displaystyle
B_{\delta}e^{-{\mathsf{J}}({u}_{2})-\epsilon}\leq\int_{v\in
B(u_{2},\delta)}e^{-{\mathsf{J}}({v})}d{v}\leq
B_{\delta}e^{-{\mathsf{J}}({u}_{2})+\epsilon},$

where $B_{\delta}$ is the Lebesgue measure of a ball with radius $\delta$.
Taking the ratio of $\alpha$’s and using the above bounds we obtain that, for
all $\delta$ sufficiently small,

$e^{{\mathsf{J}}({u}_{2})-{\mathsf{J}}({u}_{1})-2\epsilon}\leq\frac{\alpha({u}_{1},\delta)}{\alpha({u}_{2},\delta)}\leq
e^{{\mathsf{J}}({u}_{2})-{\mathsf{J}}({u}_{1})+2\epsilon}.$

Since $\epsilon$ was arbitrary the desired result follows.

∎

######  Remark 3.9.

This theorem shows that maximizing the probability of an infinitesimally small
ball is the same as minimizing the objective function ${\mathsf{J}}(\cdot).$
This is obvious in finite dimensions, but the proof above generalizes beyond
measures which possess a Lebesgue density, and may be used in infinite
dimensions.

####  3.3 Examples

By means of examples we now probe the question of whether or not the MAP
estimator captures useful information about the posterior distribution.

######  Example 3.10.

If the posterior is single-peaked, such as a Gaussian or a Laplace
distribution, as shown in Figure 2, the MAP estimator, i.e. minimizer of the
objective functions, reasonably summarizes the most likely value of the
unknown parameter.

Figure 2: Posterior (left) and objective function (right) for
$\mathcal{N}(0,1)$ posterior (orange) and $\mathrm{Laplace}(0,1)$ posterior
(blue).

We next consider several examples where a point estimator —or a
$\delta$-radius ball with $\delta\rightarrow 0$— fails to summarize adequately
the posterior distribution.

######  Example 3.11.

If the posterior is rather unevenly distributed, such as a slab-and-spike
distribution, as shown in Figure 3, then it is less clear that the MAP
estimator usefully summarizes the posterior. For example, for the case in
Figure 3 we may want the solution output of our Bayesian problem to be a
weighted average of two Gaussian distributions, or two point estimators each
with a separate mean located at one of the two minima of the objective
functions, and weight describing the probability mass associated with each of
those two points.

Figure 3: Posterior (left) and objective function (right) for a posterior that
is a sum of two Gaussian distributions, $\mathcal{N}(0,0.1^{2})$ with
probability $0.1$ and $\mathcal{N}(10,3^{2})$ with probability $0.9$.

######  Example 3.12.

In addition to a multiple-peak posterior, there are cases where the objective
function and the associated posterior density are simply very rough. In these
cases, the small-scale roughness should be ignored, while the large-scale
variation should be captured. For example in Figure 4, the objective function
is very rough, and has a unique minimizer at a point far from $0$. However, it
also has a larger pattern: it tends to be smaller around $0$, while larger
away from $0$. The MAP estimator cannot capture this large scale pattern, as
it is found by minimizing the objective function. It is arguably the case that
$u=0$ is a better point estimate. An alternative way to interpret this
phenomena is that there is a natural “temperature” to this problem, in the
sense that variations lower than this temperature could be viewed as random
noise that do not capture meaningful information.

Figure 4: Posterior (left) and objective function (right) from an objective
function that is very rough in the small scale, but contains a regular patter
on the larger scale. This specific example is generated by white noise summed
with a quadratic function for the objective function, and the posterior is
computed from the objective function.

The preceding examples suggest that multi-peak distributions, or multi-minimum
objective functions, can cause problems for MAP estimation as a point
estimator. Next we illustrate that if the dimension $d$ of the parameter
${u}\in\mathbb{R}^{d}$ is high, then a single point estimator, even if a MAP
estimator, is typically not a good summary of the posterior.

Figure 5: Empirical density of $\ell_{2}$ norm of $\mathcal{N}({0},{I})$
random vectors for various dimension: $d=1$ (blue), $d=5$ (orange), $d=10$
(green), $d=50$ (red), and $d=100$ (purple). The empirical density is obtained
from $10000$ samples for each distribution.

######  Example 3.13.

We consider what is the “typical size” of a vector ${u}$ drawn from the
standard Gaussian distribution $\mathcal{N}({0},{I})$, as the dimension
increases. In Figure 5 we display the empirical density of the norm of such
random vectors. We can see that at low dimensions, such as when $d=1$,
obtaining a value close to the mode ${u}=0$ is highly likely. In higher
dimensions, however, the probability for a vector from this distribution to
have a small $\ell_{2}-$norm becomes increasingly small as $d$ grows. For
example, let us consider the probability for the norm to be less than $5$.
Then $\operatorname{\mathbb{P}}(|{u}|<5)$ is $0.99999943$ when $d=1$,
$0.99986$ when $d=5$, $0.99465$ when $d=10$, $0.001192$ when $d=50$, and
$1.135\times 10^{-15}$ when $d=100$. So we see that, as the dimension
increases, with probability close to $1$ a sample from the posterior would
have a norm far from $0$. Indeed, for $d=1000$, the $5$th and $95$th
percentiles are respectively $30.3464$ and $32.7823$. This means when
$d=1000$, we most likely will find a vector with size around $31$, not $0$.
Another way to see this is that, since the components $u_{i}$ of $u$ are
i.i.d. standard unit Gaussians we have that, by the strong law of large
numbers,

$\frac{1}{d}\sum_{i=1}^{d}u_{i}^{2}\to 1$

as $d\to\infty$ almost surely. Thus, with high probability, the
$\ell_{2}-$norm is of size $\sqrt{d}.$ This example suggests that in high
dimension, a point estimatormay not capture enough information about the
density.

The preceding examples demonstrate that MAP estimators should be treated with
caution as they may not capture the desired posterior information in many
cases. This motivates the study of alternative ways —beyond MAP estimators— to
capture information from the posterior distribution. One such approach is to
fit one or several Gaussian distributions to the posterior by minimizing an
appropriate distance-like measure between distributions. This will be
discussed in the next chapter.

####  3.4 Discussion and Bibliography

The optimization perspective on inversion predates the development of the
Bayesian approach as a computational tool, because it is typically far cheaper
to implement. The subject of classical regularization techniques for inversion
is discussed in [31]. The concept of MAP estimators, which links probability
to optimization, is discussed in the books [63, 108] in the finite dimensional
setting. The paper [24] studies this connection precisely: it defines the MAP
estimator for infinite dimensional Bayesian inverse problems, and the
corresponding variational formulation, in a Gaussian setting. The paper [55]
studies related ideas, but in the non-Gaussian setting, and [3] generalizes
the variational formulation of MAP estimators to non-Gaussian priors that are
sparsity promoting. The paper [109] shows an application of optimization based
inversion to a large-scale geophysical application.

###  4 The Gaussian Approximation

Recall the inverse problem of finding $u$ from $y$ given by (1.1), and the
Bayesian formulation which follows from Assumption 1.1. In the previous
chapter we explored the idea of obtaining a point estimator using an
optimization perspective arising from maximizing the posterior pdf. We related
this idea to finding the center of a ball of radius $\delta$ with maximal
probability in the limit $\delta\to 0.$ Whilst the idea is intuitively
appealing, and reduces the complexity of Bayesian inference from determination
of an entire distribution to determination of a single point, the approach has
a number of limitations, in particular for noisy, multi-peaked or high
dimensional posterior distributions; the examples at the end of the previous
chapter illustrated this.

In this chapter we again adopt an optimization approach to the problem of
Bayesian inference, but instead seek a Gaussian distribution
$p=\mathcal{N}({\mu},{\Sigma})$ that minimizes the Kullback-Leibler divergence
from the posterior $\pi^{y}({u})$; since the Kullback-Leibler divergence is
not symmetric this leads to two distinct problems, both of which we will
study.

####  4.1 The Kullback-Leibler Divergence

######  Definition 4.1.

Let $\pi,\pi^{\prime}>0$ be two probability distributions on $\mathbb{R}^{d}$.
The Kullback-Leibler (K-L) divergence, or relative entropy, of $\pi$ with
respect to $\pi^{\prime}$ is defined by

$\begin{split}d_{\mbox{\tiny{\rm
KL}}}(\pi\|\pi^{\prime})&:=\int_{\mathbb{R}^{d}}\log{\mathopen{}\mathclose{{}\left(\frac{\pi(u)}{\pi^{\prime}(u)}}\right)\pi(u)du}\\\
&=\operatorname{\mathbb{E}}{{}^{\pi}\mathopen{}\mathclose{{}\left[\log{\mathopen{}\mathclose{{}\left(\frac{\pi}{\pi^{\prime}}}\right)}}\right]}\\\
&=\operatorname{\mathbb{E}}{{}^{\pi^{\prime}}\mathopen{}\mathclose{{}\left[\log{\mathopen{}\mathclose{{}\left(\frac{\pi}{\pi^{\prime}}}\right)\frac{\pi}{\pi^{\prime}}}}\right]}.\end{split}$

Kullback-Leibler is a divergence in that $d_{\mbox{\tiny{\rm
KL}}}(\pi,\pi^{\prime})\geq 0,$ with equality if and only if
$\pi=\pi^{\prime}.$ However, unlike Hellinger and total variation it is not a
distance. In particular, the K-L divergence is not symmetric: in general

$d_{\mbox{\tiny{\rm KL}}}(\pi\|\pi^{\prime})\neq d_{\mbox{\tiny{\rm
KL}}}(\pi^{\prime}\|\pi),$

a fact that will be important in this chapter. Nevertheless, it is useful for
at least four reasons: (1) it provides an upper bound for many distances; (2)
its logarithmic structure allows explicit computations that are difficult
using actual distances; (3) it satisfies many convenient analytical properties
such as being convex in both arguments and lower-semicontinuous in the
topology of weak convergence (4) it has an information theoretic and physical
interpretation.

######  Lemma 4.2.

The K-L divergence provides the following upper bounds for Hellinger and total
variation distance:

$d_{\mbox{\tiny{\rm
H}}}(\pi,\pi^{\prime})^{2}\leq\frac{1}{2}d_{\mbox{\tiny{\rm
KL}}}(\pi\|\pi^{\prime}),\quad d_{\mbox{\tiny{\rm
TV}}}(\pi,\pi^{\prime})^{2}\leq d_{\mbox{\tiny{\rm KL}}}(\pi\|\pi^{\prime}).$

###### Proof.

The second inequality follows from the first one by Lemma 1.9; thus we prove
only the first inequality. Consider the function
$\varphi:\mathbb{R}^{+}\mapsto\mathbb{R}$ defined by

$\varphi(x)=x-1-\log{x}.$

Note that

$\begin{split}\varphi^{\prime}(x)&=1-\frac{1}{x},\\\
\varphi^{\prime\prime}(x)&=\frac{1}{x^{2}},\\\
\varphi(\infty)&=\varphi(0)=\infty.\end{split}$

Thus the function is convex on its domain. As the minimum of $\varphi$ is
attained at $x=1$, and as $\varphi(1)=0$, we deduce that $\varphi(x)\geq 0$
for all $x\in(0,\infty).$ Hence,

$\begin{split}x-1&\geq\log{x}\hskip 28.45274pt\text{ for all }x\geq 0,\\\
\sqrt{x}-1&\geq\frac{1}{2}\log{x}\hskip 28.45274pt\text{for all }x\geq
0.\end{split}$

We can use this last inequality to bound the Hellinger distance:

$\begin{split}d_{\mbox{\tiny{\rm
H}}}(\pi,\pi^{\prime})^{2}&=\frac{1}{2}\int\mathopen{}\mathclose{{}\left(1-\sqrt{\frac{\pi^{\prime}}{\pi}}}\right)^{2}\pi
du\\\
&=\frac{1}{2}\int\mathopen{}\mathclose{{}\left(1+\frac{\pi^{\prime}}{\pi}-2\sqrt{\frac{\pi^{\prime}}{\pi}}}\right)\pi
du\\\
&=\int\mathopen{}\mathclose{{}\left(1-\sqrt{\frac{\pi^{\prime}}{\pi}}}\right)\pi
du\leq-\frac{1}{2}\int\log\Bigl{(}{\frac{\pi^{\prime}}{\pi}}\Bigr{)}\pi
du=\frac{1}{2}d_{\mbox{\tiny{\rm KL}}}(\pi\|\pi^{\prime}).\end{split}$

∎

####  4.2 Best Gaussian Fit By Minimizing $d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$

In this section we prove the existence of an approximating Gaussian. To
streamline the exposition we work under the following assumption, which can be
relaxed:

######  Assumption 4.3.

The posterior distribution is constructed under the following assumptions:

  * •

The loss function $\mathsf{L}(u):=-\log\nu\bigl{(}y-G(u)\bigr{)}$ is non-
negative and bounded above.

  * •

The prior is a centered isotropic Gaussian:
$\rho\sim\mathcal{N}(0,\lambda^{-1}{I}).$

Let $\mathcal{A}$ be the set of Gaussian distributions on $\mathbb{R}^{d}$
with positive definite covariance,

$\mathcal{A}=\\{\mathcal{N}({\mu},{\Sigma}):{\mu}\in\mathbb{R}^{d},{\Sigma}\in\mathbb{R}^{d\times
d}\text{ positive-definite symmetric}\\}.$

We have the following theorem:

######  Theorem 4.4 (Best Gaussian Approximation).

Under Assumption 4.3, there exists at least one probability distribution
$p\in\mathcal{A}$ at which the infimum

$\inf_{p\in\mathcal{A}}d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$

is attained.

###### Proof.

The K-L divergence can be computed explicitly as

$\begin{split}d_{\mbox{\tiny{\rm
KL}}}(p\|\pi)=&\operatorname{\mathbb{E}}^{p}{\log{p}}-\operatorname{\mathbb{E}}^{p}\log\pi\\\
=&\operatorname{\mathbb{E}}^{p}\mathopen{}\mathclose{{}\left(-\frac{1}{2}|{u}-{\mu}|^{2}_{{\Sigma}}-\frac{1}{2}\log\mathopen{}\mathclose{{}\left((2\pi)^{d}\text{det}{\Sigma}}\right)}\right.+\mathopen{}\mathclose{{}\left.\mathsf{L}({u})+\frac{\lambda}{2}|{u}|^{2}+\log
Z}\right).\end{split}$

Note that $Z$ is the normalization constant for $\pi$ and is independent of
$p$ and hence of ${\mu}$ and ${\Sigma}.$ We can represent a given random
variable ${u}\sim p$ by writing ${u}={\mu}+{\Sigma}^{1/2}\xi$, where
$\xi\sim\mathcal{N}({0},{{I}})$, and hence

$|{u}|^{2}=|{\mu}|^{2}+|\Sigma^{1/2}\xi|^{2}+2\langle{\mu},{\Sigma}^{1/2}\xi\rangle$

to obtain

$\begin{split}d_{\mbox{\tiny{\rm
KL}}}(p\|\pi)&=-\frac{d}{2}-\frac{d}{2}\log(2\pi)-\frac{1}{2}\log\text{det}{\Sigma}+\operatorname{\mathbb{E}}^{p}\mathsf{L}({u})+\frac{\lambda}{2}|{\mu}|^{2}+\frac{\lambda}{2}\text{tr}({\Sigma})+\log
Z.\\\ \end{split}$

Define
${\mathsf{I}}(\mu,\Sigma)=\operatorname{\mathbb{E}}^{p}\mathsf{L}({u})+\frac{\lambda}{2}|{\mu}|^{2}+\frac{\lambda}{2}\text{tr}({\Sigma})-\frac{1}{2}\log\text{det}{\Sigma}.$
Clearly, there is a correspondence between minimizing $d_{\mbox{\tiny{\rm
KL}}}(p\|\pi)$ over $p\in{\cal A}$ and minimizing ${\mathsf{I}}(\mu,\Sigma)$
over $\mu\in\mathbb{R}^{d}$ and positive definite $\Sigma.$ Note that:

  * •

${\mathsf{I}}(0,{I})<\infty.$

  * •

For any $\Sigma,$ ${\mathsf{I}}(\mu,\Sigma)\to\infty$ as $|\mu|\to\infty.$

  * •

For any $\mu,$ ${\mathsf{I}}(\mu,\Sigma)\to\infty$ as $\text{tr}(\Sigma)\to 0$
or $\text{tr}(\Sigma)\to\infty.$

Therefore, there are $M,r,R>0$ such that the infimum of
${\mathsf{I}}(\mu,\Sigma)$ over $\mu\in\mathbb{R}^{d}$ and positive definite
$\Sigma$ is equal to the infimum of ${\mathsf{I}}(\mu,\Sigma)$ over

$\tilde{\mathcal{A}}:=\\{({\mu},{\Sigma}):{\mu}\in\mathbb{R}^{d},{\Sigma}\in\mathbb{R}^{d\times
d}\text{ positive-definite symmetric},|\mu|\leq M,r\leq\text{tr}(\Sigma)\leq
R\\}.$

Since ${\mathsf{I}}$ is continuous in $\tilde{\mathcal{A}}$ it achieves its
infimum and the proof is complete. ∎

We remark that the theorem establishes the existence of a best Gaussian
approximation. However, minimizers need not be unique.

####  4.3 Best Gaussian Fit By Minimizing $d_{\mbox{\tiny{\rm KL}}}(\pi\|p)$

In this section we show that the best Gaussian approximation in Kullback-
Leibler with respect to its second argument is unique and given by moment
matching.

######  Theorem 4.5 (Best Gaussian by Moment Matching).

Assume that $\bar{\mu}:=\operatorname{\mathbb{E}}^{\pi}[u]$ is finite and that
$\bar{\Sigma}:=\operatorname{\mathbb{E}}^{\pi}[(u-\bar{\mu})\otimes(u-\bar{\mu})]$
is positive-definite. Then the infimum

$\inf_{p\in\mathcal{A}}d_{\mbox{\tiny{\rm KL}}}(\pi\|p)$

is attained at the element in $\mathcal{A}$ with mean $\bar{\mu}$ and
covariance $\bar{\Sigma}.$

###### Proof.

By definition

$\displaystyle d_{\mbox{\tiny{\rm
KL}}}(\pi\|p)=-\operatorname{\mathbb{E}}^{\pi}{[\log{p}]}+\operatorname{\mathbb{E}}^{\pi}{[\log{\pi}]}.$
(4.1)

Since the second term does not involve $p$, we study minimization of

$\displaystyle-\operatorname{\mathbb{E}}^{\pi}{[\log{p}]}$
$\displaystyle=-\operatorname{\mathbb{E}}^{\pi}{\mathopen{}\mathclose{{}\left[\log\mathopen{}\mathclose{{}\left(\frac{1}{\sqrt{(2\pi)^{d}\text{det}\Sigma}}\text{exp}\mathopen{}\mathclose{{}\left(-\frac{1}{2}\big{|}u-\mu\big{|}_{\Sigma}^{2}}\right)}\right)}\right]}$
$\displaystyle=\frac{1}{2}\operatorname{\mathbb{E}}^{\pi}{\mathopen{}\mathclose{{}\left[\big{|}u-\mu\big{|}_{\Sigma}^{2}}\right]}+\frac{1}{2}\log{\text{det}\Sigma}+\frac{d}{2}\log{2\pi}.$

Let $\Omega=\Sigma^{-1}$. Then our task is equivalent to minimizing the
following function of $\mu$ and $\Omega$:

$\displaystyle{\mathsf{I}}(\mu,\Omega)=\frac{1}{2}\operatorname{\mathbb{E}}^{\pi}{\mathopen{}\mathclose{{}\left[\langle
u-\mu,\Omega(u-\mu)\rangle}\right]-\frac{1}{2}\log{\text{det}\Omega}}.$

First we find the critical points of ${\mathsf{I}}$ by taking its first order
partial derivative with respect to $\mu$ and $\Omega$ and setting both to
zero:

$\displaystyle\partial_{\mu}{\mathsf{I}}$
$\displaystyle=-\operatorname{\mathbb{E}}^{\pi}\mathopen{}\mathclose{{}\left[\Omega(u-\mu)}\right]=0;$
$\displaystyle\partial_{\Omega}{\mathsf{I}}$
$\displaystyle=\frac{1}{2}\partial_{\Omega}\mathopen{}\mathclose{{}\left(\operatorname{\mathbb{E}}^{\pi}\mathopen{}\mathclose{{}\left[(u-\mu)\otimes(u-\mu):\Omega}\right]}\right)-\frac{1}{2\text{det}\Omega}\partial_{\Omega}\text{det}\Omega$
$\displaystyle=\frac{1}{2}\operatorname{\mathbb{E}}^{\pi}\mathopen{}\mathclose{{}\left[(u-\mu)\otimes(u-\mu)}\right]-\frac{1}{2}\Omega^{-1}=0;\quad\quad$

where we have used the relation
$\partial_{\Omega}\text{det}\Omega=\text{det}{\Omega}\cdot\Omega^{-1}.$
Solving the above two equations gives us the critical point, expressed in
terms of mean and covariance,

$(\bar{\mu},\bar{\Sigma})=(\operatorname{\mathbb{E}}^{\pi}[u],\operatorname{\mathbb{E}}^{\pi}[(u-\bar{\mu})\otimes(u-\bar{\mu})]).$

Next we will show $(\bar{\mu},\bar{\Sigma}^{-1})$ is a minimizer of
${\mathsf{I}}$ or, equivalently, that $(\bar{\mu},\bar{\Sigma})$ is a
minimizer of the expression given in (4.1). To this end, if we deal with a
vector $\theta$ parametrizing the distribution $p$, it will be sufficient to
show that the Hessian of equation (4.1) with respect to $\theta$ is positive
definite. Recalling the label $d$ of the dimension of the parameter $u$ of
unknowns we then have:

$\displaystyle p_{\theta}(u)$
$\displaystyle=\sqrt{\frac{\text{det}\Omega}{(2\pi)^{d}}}\text{exp}\mathopen{}\mathclose{{}\left(-\frac{1}{2}(u-\mu)^{T}\Omega(u-\mu)}\right)$
$\displaystyle=\sqrt{\frac{\text{det}\Omega}{(2\pi)^{d}}}\text{exp}\mathopen{}\mathclose{{}\left[-\mathopen{}\mathclose{{}\left(\frac{1}{2}u^{T}\Omega
u-\mu^{T}\Omega u+\frac{1}{2}\mu^{T}\Omega\mu}\right)}\right]$
$\displaystyle=\sqrt{\frac{\text{det}\Omega}{(2\pi)^{d}}}\text{exp}\mathopen{}\mathclose{{}\left(-\frac{1}{2}\mu^{T}\Omega\mu}\right)\text{exp}\mathopen{}\mathclose{{}\left(-\frac{1}{2}u^{T}\Omega
u+\mu^{T}\Omega u}\right).$

Utilizing vectorization of matrices we may write:

$\displaystyle u^{T}\Omega u$
$\displaystyle=\Omega:uu^{T}=[\text{vec}(\Omega)]^{T}\text{vec}(uu^{T})$
$\displaystyle\mu^{T}\Omega u$ $\displaystyle=(\Omega\mu)^{T}u$
$\displaystyle\Rightarrow-\frac{1}{2}u^{T}\Omega u+\mu^{T}\Omega u$
$\displaystyle=\begin{bmatrix}\Omega\mu\\\
-\frac{1}{2}\text{vec}(\Omega)\end{bmatrix}^{T}\begin{bmatrix}u\\\
\text{vec}(uu^{T})\end{bmatrix}.$

Then we let:

$\displaystyle\theta=\begin{bmatrix}\Omega\mu\\\
-\frac{1}{2}\text{vec}(\Omega)\end{bmatrix},\quad\quad
T(u)=\begin{bmatrix}u\\\ \text{vec}(uu^{T})\end{bmatrix}.$

The pdf $p_{\theta}(u)$ can then be written in the following form:

$\displaystyle p_{\theta}(u)$
$\displaystyle=\frac{1}{Z(\theta)}\text{exp}\mathopen{}\mathclose{{}\left(\theta^{T}T(u)}\right)$
(4.2) $\displaystyle\text{with}\quad Z(\theta)$
$\displaystyle=\int_{\mathbb{R}^{d}}\text{exp}\mathopen{}\mathclose{{}\left(\theta^{T}T(u)}\right)du.$
(4.3)

Using equation (4.2) we can rewrite equation (4.1) as:

$\displaystyle H(\theta)=d_{\mbox{\tiny{\rm
KL}}}(\pi\|p_{\theta})=-\theta^{T}\operatorname{\mathbb{E}}^{\pi}\mathopen{}\mathclose{{}\left[T(u)}\right]+\log\mathopen{}\mathclose{{}\left(Z(\theta)}\right)+\operatorname{\mathbb{E}}^{\pi}\mathopen{}\mathclose{{}\left[\pi(u)}\right].$
(4.4)

Notice that

$\displaystyle\nabla_{\theta}\log\mathopen{}\mathclose{{}\left(Z(\theta)}\right)$
$\displaystyle=\frac{1}{Z(\theta)}\int_{\mathbb{R}^{d}}\nabla_{\theta}\mathopen{}\mathclose{{}\left[\text{exp}\mathopen{}\mathclose{{}\left(\theta^{T}T(u)}\right)}\right]du$
$\displaystyle=\frac{1}{Z(\theta)}\int_{\mathbb{R}^{d}}T(u)\text{exp}\mathopen{}\mathclose{{}\left(\theta^{T}T(u)}\right)du$
$\displaystyle=\operatorname{\mathbb{E}}^{p_{\theta}}\mathopen{}\mathclose{{}\left[T(u)}\right].$

Therefore we can calculate the gradient and Hessian of $H(\theta)$ as follows:

$\displaystyle\nabla_{\theta}H(\theta)$
$\displaystyle=-\operatorname{\mathbb{E}}^{\pi}\mathopen{}\mathclose{{}\left[T(u)}\right]+\operatorname{\mathbb{E}}^{p_{\theta}}\mathopen{}\mathclose{{}\left[T(u)}\right],$
$\displaystyle\mathopen{}\mathclose{{}\left[\nabla^{2}_{\theta}H(\theta)}\right]_{ij}$
$\displaystyle=\frac{\partial\log\mathopen{}\mathclose{{}\left(Z(\theta)}\right)}{\partial\theta_{i}\theta_{j}}=\frac{\partial}{\partial\theta_{j}}\mathopen{}\mathclose{{}\left(\frac{1}{Z(\theta)}\int_{\mathbb{R}^{d}}T_{i}(u)e^{\theta^{T}T(u)}du}\right)$
$\displaystyle=\frac{1}{Z(\theta)}\frac{\partial}{\partial\theta_{j}}\int_{\mathbb{R}^{d}}T_{i}(u)e^{\theta^{T}T(u)}dx+\int_{\mathbb{R}^{d}}T_{i}(u)e^{\theta^{T}T(u)}du\cdot\frac{\partial}{\partial\theta_{j}}Z(\theta)^{-1}$
$\displaystyle=\frac{1}{Z(\theta)}\int_{\mathbb{R}^{d}}T_{i}(u)T_{j}(u)e^{\theta^{T}T(u)}du$
$\displaystyle\quad\quad\quad-\frac{1}{Z(\theta)^{2}}\int_{\mathbb{R}^{d}}T_{i}(u)e^{\theta^{T}T(u)}du\cdot\int_{\mathbb{R}^{d}}T_{j}(u)e^{\theta^{T}T(u)}du$
$\displaystyle=\operatorname{\mathbb{E}}^{p_{\theta}}\mathopen{}\mathclose{{}\left[T_{i}T_{j}}\right]-\operatorname{\mathbb{E}}^{p_{\theta}}\mathopen{}\mathclose{{}\left[T_{i}}\right]\operatorname{\mathbb{E}}^{p_{\theta}}\mathopen{}\mathclose{{}\left[T_{j}}\right]$
$\displaystyle=\mathopen{}\mathclose{{}\left[\text{Cov}^{p_{\theta}}(T)}\right]_{ij}.$

Therefore the Hessian of the objective function is the covariance matrix of
$T(u)$ under the multivariate normal distribution $p_{\theta}$ and by
construction, is positive semi-definite. Therefore $d_{\mbox{\tiny{\rm
KL}}}(\pi\|p_{\theta})$ is convex in $\theta$ and the critical point
$(\bar{\mu},\bar{\Sigma})$ is a corresponding minimizer. ∎

######  Remark 4.6.

Notice that the preceding proof of convexity holds for any distribution $p$
that can alternatively be parametrized by a vector $\theta$ with the general
form of equation (4.2). In particular for such problems the convexity of
$d_{\mbox{\tiny{\rm KL}}}(\pi\|p_{\theta})$ follows, ensuring the existence of
a minimizer. In fact, equation (4.2) is a reduced form of the following more
general expression:

$\displaystyle p_{\theta}(u)$
$\displaystyle=h(u)\text{exp}\mathopen{}\mathclose{{}\left(\theta^{T}T(u)-A(\theta)}\right)$
(4.5) $\displaystyle\text{with}\quad\quad A(\theta)$
$\displaystyle=\log\mathopen{}\mathclose{{}\left[\int_{\mathbb{R}^{d}}h(u)\text{exp}\mathopen{}\mathclose{{}\left(\theta^{T}T(u)}\right)du}\right].$
(4.6)

Since $h(u)$ is independent of $\theta$, the conclusion of the previous
theorem carries over to distributions with the form of (4.5). Such
distributions are said to termed the exponential family in the statistics
literature, and $\theta$ is called the natural parameter, $T(u)$ the
sufficient statistic, $h(u)$ the base measure, and $A(\theta)$ the log-
partition. The Gaussian distribution is a special case in which $h(u)$ is
constant with respect to $u$. In what follows we show some examples of other
exponential family members.

######  Example 4.7.

The Bernoulli distribution
Let $u$ be a random variable following the Bernoulli distribution with
$\operatorname{\mathbb{P}}(u=1)=q$ and $\operatorname{\mathbb{P}}(u=0)=1-q$.
Then the point mass function describing the distribution of $u$ is:

$\displaystyle p(u;q)$ $\displaystyle=q^{u}(1-q)^{1-u}$
$\displaystyle=\text{exp}\mathopen{}\mathclose{{}\left[u\log\mathopen{}\mathclose{{}\left(\frac{q}{1-q}}\right)+\log(1-q)}\right]$
$\displaystyle=(1-q)\text{exp}\mathopen{}\mathclose{{}\left[u\log\mathopen{}\mathclose{{}\left(\frac{q}{1-q}}\right)}\right].$

Therefore $\theta=\log\mathopen{}\mathclose{{}\left(\frac{q}{1-q}}\right)$,
$T(u)=u$, $A(\theta)=\log(1+e^{\theta})$, and $h(u)=1$.

######  Example 4.8.

The Poisson distribution
The point mass function describing the distribution of $u$ is, in this case,

$\displaystyle p(u;\lambda)$
$\displaystyle=\frac{\lambda^{u}e^{-\lambda}}{u!}$
$\displaystyle=\frac{1}{u!}e^{u\log\lambda-\lambda}.$

Therefore $\theta=\log{\lambda}$, $T(u)=u$, $A(\theta)=e^{\theta}$ and
$h(u)=\frac{1}{u!}$.

####  4.4 Comparison between $d_{\mbox{\tiny{\rm KL}}}(\pi\|p)$ and
$d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$

It is instructive to compare the two different minimization problems, both
leading to a “best Gaussian”, that we described in the preceding two sections.
We write the two relevant divergences as follows and then explain the
nomenclature:

$\displaystyle d_{\mbox{\tiny{\rm
KL}}}(p\|\pi)=\operatorname{\mathbb{E}}^{p}{\mathopen{}\mathclose{{}\left[\log\mathopen{}\mathclose{{}\left(\frac{p}{\pi}}\right)}\right]}=\operatorname{\mathbb{E}}^{p}{[\log{p}]}-\operatorname{\mathbb{E}}^{p}{[\log{{\pi}}]}\quad\text{"Mode-
seeking"}$ $\displaystyle d_{\mbox{\tiny{\rm
KL}}}(\pi\|p)=\operatorname{\mathbb{E}}^{\pi}{\mathopen{}\mathclose{{}\left[\log\mathopen{}\mathclose{{}\left(\frac{\pi}{p}}\right)}\right]}=\operatorname{\mathbb{E}}^{\pi}{[\log{\pi}]}-\operatorname{\mathbb{E}}^{\pi}{[\log{{p}}]}\quad\text{"Mean-
seeking"}$

(a) Minimizing $d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$

(b) Minimizing $d_{\mbox{\tiny{\rm KL}}}(\pi\|p)$

Figure 6: (a) Minimizing $d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$ can lead to
serious information loss while (b) minimizing $d_{\mbox{\tiny{\rm
KL}}}(\pi\|p)$ ensures a comprehensive consideration of all components of
$\pi$.

Note that when minimizing $d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$ we want
$\log{\frac{p}{\pi}}$ to be small, which can happen when $p\simeq\pi$ or
$p\ll\pi$. This illustrates the fact that minimizing $d_{\mbox{\tiny{\rm
KL}}}(p\|\pi)$ may miss out components of $\pi$. For example, in Figure 6(a)
$\pi$ is a bi-modal like distribution but minimizing $d_{\mbox{\tiny{\rm
KL}}}(p\|\pi)$ over Gaussians $p$ can only give a single mode approximation
which is achieved by matching one of the modes; we may think of this as “mode-
seeking”. In contrast, when minimizing $d_{\mbox{\tiny{\rm KL}}}(\pi\|p)$ over
Gaussians $p$ we want $\log{\frac{\pi}{p}}$ to be small where $p$ appears as
the denominator. This implies that wherever $\pi$ has some mass we must let
$p$ also have some mass there in order to keep $\frac{\pi}{p}$ as close as
possible to one. Therefore the minimization is carried out by allocating the
mass of $p$ in a way such that on average the divergence between $p$ and $\pi$
attains its minimum as shown in Figure 6(b); hence the label “mean-seeking.”
Different applications will favor different choices between the mean and mode
seeking approaches to Gaussian approximation.

####  4.5 Variational Formulation of Bayes Theorem

This chapter has been concerned with finding the best Gaussian approximation
to a measure with respect to KL divergences. Bayes Theorem 1.2 itself can be
formulated through a closely related minimization principle. Consider a
posterior $\pi^{y}(u)$ in the following form:

$\displaystyle\pi^{y}(u)=\frac{1}{Z}\text{exp}\bigl{(}-\mathsf{L}(u)\bigr{)}\rho(u),$

where $\rho(u)$ is the prior, $\mathsf{L}(u)$ is the negative log-likelihood,
and $Z$ the normalization constant. We assume here for exposition that all
densities are positive. Then we can express $d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$
in terms of the prior as follows:

$\displaystyle d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$
$\displaystyle=\int_{\mathbb{R}^{d}}\text{log}\mathopen{}\mathclose{{}\left(\frac{p}{\pi}}\right)pdu$
$\displaystyle=\int_{\mathbb{R}^{d}}\text{log}\mathopen{}\mathclose{{}\left(\frac{p}{\rho}\frac{\rho}{\pi}}\right)pdu$
$\displaystyle=\int_{\mathbb{R}^{d}}\text{log}\mathopen{}\mathclose{{}\left(\frac{p}{\rho}\text{exp}\bigl{(}\mathsf{L}(u)\bigr{)}Z}\right)pdu$
$\displaystyle=d_{\mbox{\tiny{\rm
KL}}}(p\|\rho)+\operatorname{\mathbb{E}}^{p}[\mathsf{L}(u)]+\log{Z}.$

If we define

${\mathcal{J}}(p)=d_{\mbox{\tiny{\rm
KL}}}(p\|\rho)+\operatorname{\mathbb{E}}^{p}[\mathsf{L}(u)]$

then we have the following:

######  Theorem 4.9 (Bayes Theorem as an Optimization Principle).

The posterior distribution $\pi$ is given by the following minimization
principle:

$\pi=\text{\rm argmin}_{p\in{\mathscr{P}}}{\mathcal{J}}(p),$

where ${\mathscr{P}}$ contains all probability densities on $\mathbb{R}^{d}$.

###### Proof.

Since $Z$ is the normalization constant for $\pi$ and is independent of $p$,
the minimizer of $d_{\mbox{\tiny{\rm KL}}}(p\|\pi)$ will also be the minimizer
of ${\mathcal{J}}(p).$ Since the global minimizer of $d_{\mbox{\tiny{\rm
KL}}}(p\|\pi)$ is attained at $p=\pi$ the result follows. ∎

Why is it useful to view the posterior as the minimizer of an energy? There
are at least three advantages of this viewpoint. First, the variational
formulation provides a natural way to approximate the posterior by restricting
the minimization problem to distributions satisfying some computationally
desirable property. For instance, variational Bayes methods often restrict the
minimization to densities with product structure and in this chapter we have
studied restriction to the class of Gaussian distributions. Second,
variational formulations allow to show convergence of posterior distributions
indexed by some parameter of interest by studying the Gamma-convergence of the
associated objective functionals. Third, variational formulations provide
natural paths, defined by a gradient flow, towards the posterior.
Understanding these flows and their rates of convergence is helpful in the
choice of sampling algorithms.

####  4.6 Discussion and Bibliography

The definition of the K-L divergence, and upper-bounds in terms of probability
metrics, can be found in [46]. For a basic introduction to variational
methods, and the moment-matching version of Gaussian approximation, see [13].
The problem of finding a Gaussian approximation of a general finite
dimensional probability distribution is studied in [80], and infinite
dimensional formulations are considered in [94] and the companion paper [93].
The approximation in Theorem 4.4 consists of a single Gaussian distribution.
If the posterior has more than one mode, a single Gaussian may not be
appropriate. For an approximation composed of Gaussian mixtures, the reader is
referred to [80]. The paper [44] highlights how minimization of K-L divergence
arises naturally in the optimization of local entropy and heat regularized
costs in deep learning. The formulation of Bayes Theorem as an optimization
principle is well-known and widely used in the machine learning community
where it goes under the name “variational Bayes”; see the book [81] and the
paper [8] for clear expositions of this subject. The variational formulation
has been used to establish convergence of Bayesian procedures in [43] and [41]
and to define gradient flows that converge to the posterior as a way to
develop and analyze Markov chain Monte Carlo proposals in [42]. For more
information about the properties and examples of exponential family
distributions, we refer to [90]. and for background on matrix calculations
that were used in this chapter we refer to [92].

###  5 Monte Carlo and Importance Sampling

In this chapter we introduce Monte Carlo sampling and importance sampling.
They are two general techniques for estimating expectations computed with
respect to a particular distribution by generating samples from (a possibly
different) distribution; they may also be viewed as approximations of a given
measure via sums of Dirac measures.

Throughout this chapter we focus interest on computing expectations with
respect to a probability distribution with pdf $\pi$ given in the form

$\pi(u)=\frac{1}{Z}g(u)\rho(u),$ (5.1)

where $Z$ is a normalizing constant. Monte Carlo sampling will use samples
from $\pi$ itself; importance sampling will use samples from $\rho$. A
particular application of this setting is where $\rho$ is the prior, $\pi$ the
posterior and $g$ the likelihood. Let
$f:\mathbb{R}^{d}\longrightarrow\mathbb{R}$ denote the function whose
expectation we are interested in computing. For any pdf $p$ on
$\mathbb{R}^{d}$ we write

$p(f)=\operatorname{\mathbb{E}}^{p}\mathopen{}\mathclose{{}\left[f(u)}\right]=\int_{\mathbb{R}^{d}}f(u)p(u)du.$
(5.2)

In this chapter we will generalize the concept of pdf to include Dirac mass
distributions. A Dirac mass at $v$ will be viewed as having pdf
$\delta(\cdot-v)$ where $\delta(\cdot)$ integrates to one and takes the value
zero everywhere except at the origin.

####  5.1 Monte Carlo Sampling

If we have $N$ random samples $u^{(1)},\ldots,u^{(N)}$, generated i.i.d.
according to $\pi$, then we can estimate $\pi$ by the Monte Carlo estimator
$\pi_{MC}^{N}$, which is defined as

$\pi^{N}_{MC}:=\frac{1}{N}\sum_{n=1}^{N}\delta(u-u^{(n)}).$ (5.3)

This gives rise to the following estimator of $\pi(f):$

$\pi^{N}_{MC}(f)=\frac{1}{N}\sum_{n=1}^{N}f(u^{(n)}),\,\,u^{(n)}\sim\pi\quad{\rm
i.i.d.}$

######  Theorem 5.1 (Monte Carlo Error).

For $f:\mathbb{R}^{d}\longrightarrow\mathbb{R}$ denote
$|f|_{\infty}:=\sup_{u\in\mathbb{R}^{d}}|f(u)|.$ We have

$\begin{split}\sup_{|f|_{\infty}\leq
1}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\pi^{N}_{MC}(f)-\pi(f)}\right]}}\right|=0,\\\
\sup_{|f|_{\infty}\leq
1}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{MC}(f)-\pi(f)}\right)^{2}}\right]}}\right|\leq\frac{1}{N}.\end{split}$

###### Proof.

Define

$\bar{f}(u)=f(u)-\pi(f).$

To prove the first result, namely that the estimator is unbiased, note that

$\begin{split}&\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\pi^{N}_{MC}(f)-\pi(f)}\right]}\\\
&=\frac{1}{N}\sum_{n=1}^{N}\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[f\mathopen{}\mathclose{{}\left(u^{(n)}}\right)-\pi(f)}\right]}\\\
&=\frac{1}{N}\sum_{n=1}^{N}\bigl{(}\pi\mathopen{}\mathclose{{}\left(f}\right)-\pi\mathopen{}\mathclose{{}\left(f}\right)\bigr{)}\\\
&=\frac{1}{N}\cdot 0=0.\end{split}$

Therefore the supremum of its absolute value is also zero. For the second
result, which bounds the variance of the estimator, we observe that
$\operatorname{\mathbb{E}}{[\bar{f}]}=0$ and, then,

$\begin{split}\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{MC}(f)-\pi(f)}\right)^{2}}\right]}&=\frac{1}{N^{2}}\sum_{n=1}^{N}\sum_{n=1}^{N}\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\bar{f}\mathopen{}\mathclose{{}\left(u^{(n)}}\right)\bar{f}\mathopen{}\mathclose{{}\left(u^{(n)}}\right)}\right]}\\\
&=\frac{1}{N^{2}}\sum_{n=1}^{N}\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\bar{f}\mathopen{}\mathclose{{}\left(u^{(n)}}\right)^{2}}\right]}\\\
&=\frac{1}{N}\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\bar{f}\mathopen{}\mathclose{{}\left(u^{(1)}}\right)^{2}}\right]}=\frac{1}{N}\text{Var}_{\pi}[f]\end{split}$

since $u^{(n)}$ are i.i.d. In particular we have

$\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{MC}(f)-\pi(f)}\right)^{2}}\right]}=\frac{1}{N}\text{Var}_{\pi}[f]\leq\frac{1}{N}\pi(f^{2})$
(5.4)

since

$\text{Var}_{\pi}[f]=\pi(f^{2})-\pi(f)^{2}\leq\pi(f^{2}).$

Therefore

$\sup_{|f|_{\infty}\leq
1}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{MC}(f)-\pi(f)}\right)^{2}}\right]}}\right|=\sup_{|f|_{\infty}\leq
1}\mathopen{}\mathclose{{}\left|\frac{1}{N}\text{Var}_{\pi}[f]}\right|\leq\frac{1}{N}.$

∎

The theorem shows that the Monte Carlo estimator $\pi^{N}_{MC}$ is an unbiased
approximation for the posterior $\pi$ and that, by choosing $N$ large enough,
expectation of any bounded function $f$ can in principle be approximated by
Monte Carlo sampling to arbitrary accuracy. Furthermore, although the
convergence is slow with respect to $N$ there is no dependence on the
dimension of the problem or on the properties of $f$, other than its supremum.

######  Example 5.2 (Approximation of an Integral).

Let $f:\mathbb{R}\longrightarrow\mathbb{R}$ be a sigmoid function defined on
$\mathbb{R}$ and shown in Figure 7(a) below as the blue solid curve. Let
$u\sim\pi=0.1~{}\mathcal{N}(-5,1)+0.9~{}\mathcal{N}(5,1)$ be a Gaussian
mixture consisting of two Gaussian distributions. We wish to approximate the
expected value of $f(u)\times\mathbb{I}_{[a,b]}(u)$ and
$\mathbb{I}_{[a,b]}(u)=\begin{cases}1&\text{ if }u\in[a,b]\\\ 0&\text{
otherwise}\end{cases}$. We use Monte Carlo sampling to generate $N$ random
samples $u^{(1)},\ldots,u^{(N)}$ and compute the error between the actual
integral and Monte Carlo estimator. The integral and estimator are in the
form:

$\begin{split}&\pi(f)=\int_{a}^{b}f(u)\pi(u)\mathop{}\\!\mathrm{d}u,\\\
&\pi^{N}_{MC}(f)=\frac{1}{N}\sum_{n=1}^{N}f(u^{(n)})\mathbb{I}_{[a,b]}(u^{(n)}).\end{split}$

The results of a set of numerical experiments with $a=-5,b=5$ and varying $N$
are shown in Figure 7(b). A randomly chosen subset of the samples used when
$N=100$ is displayed in Figure 7(a).

---

---

Figure 7: Large sampling number reduces the estimation error by Monte Carlo
method.

####  5.2 Importance Sampling

Standard Monte Carlo sampling can only be used when it is possible to sample
from the desired target distribution $\pi$. When it is not possible to sample
from $\pi$, we can draw samples from another proposal distribution $\rho$
instead. We then need to evaluate $g$ in (5.1) at each sample and use it as an
importance weight in the approximation of the desired distribution. This is
the idea of importance sampling.

Consider $\pi$ as in equation (5.1). Given a test function
$f:\mathbb{R}^{d}\to\mathbb{R}$ we can rewrite its expectation with respect to
$\pi$ in terms of expected values with respect to $\rho$ as follows:

$\pi(f)=\frac{\rho(fg)}{\rho(g)},$

where we used that $Z=\rho(g).$ Approximating with standard Monte Carlo the
numerator and the denominator gives

$\displaystyle\pi(f)$
$\displaystyle\approx\sum_{n=1}^{N}w^{(n)}f\mathopen{}\mathclose{{}\left(u^{(n)}}\right),\,\,u^{(n)}\sim\rho\,\,{\rm
i.i.d.}$ $\displaystyle=\pi^{N}_{IS}(f),$

where

$w^{(n)}:=\frac{g\mathopen{}\mathclose{{}\left(u^{(n)}}\right)}{\sum_{m=1}^{N}g\mathopen{}\mathclose{{}\left(u^{(m)}}\right)},\quad\quad\pi^{N}_{IS}:=\sum_{n=1}^{N}w^{(n)}\delta(u-u^{(n)}).$

Thus, given $N$ samples $u^{(1)},\ldots,u^{(N)}$ generated i.i.d. according to
$\rho$ we can estimate $\pi$ with the _particle approximation measure_
$\pi^{N}_{IS}.$ The quality of this estimator may be assessed by considering
the worst-case bias and mean squared error when using the particle
approximation measure to estimate expectations over the class of bounded test
functions $\\{f:\mathbb{R}^{d}\to\mathbb{R}:|f|_{\infty}=1\\}.$ We will show
that worst-case error upper-bounds can be obtained in terms of the $\chi^{2}$
divergence between the target and the proposal, suggesting that the
performance of importance sampling depends on the closeness between target and
proposal.

######  Definition 5.3.

Let $\pi>0$ and $\pi^{\prime}>0$ be two probability distributions on
$\mathbb{R}^{d}.$ The _$\chi^{2}$ divergence_ of $\rho$ with respect to
$\rho^{\prime}$ is

$d_{\mbox{\tiny{$\chi^{2}$}}}(\pi\|\pi^{\prime}):=\int_{\mathbb{R}^{d}}\Bigl{(}\frac{\pi(u)}{\pi^{\prime}(u)}-1\Bigr{)}^{2}\,\pi^{\prime}(u)du.$
(5.5)

The $\chi^{2}$ is a divergence but not a distance as it is, in general, not
symmetric. The next result shows that, similarly as for standard Monte Carlo,
the mean squared error of $\pi^{N}_{IS}(f)$ as an estimator of $\pi(f)$ is
order $N^{-1}.$ However, there are two main differences: the estimator is now
biased, and the constant in the mean squared error depends on the $\chi^{2}$
divergence between the target and the proposal.

######  Theorem 5.4 (Importance Sampling Error).

We have

$\begin{split}\sup_{|f|_{\infty}\leq
1}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\pi^{N}_{IS}(f)-\pi(f)}\right]}}\right|\leq
2\frac{1+d_{\mbox{\tiny{$\chi^{2}$}}}(\pi\|\rho)}{N},\\\
\sup_{|f|_{\infty}\leq
1}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{IS}(f)-\pi(f)}\right)^{2}}\right]}}\right|\leq
4\frac{1+d_{\mbox{\tiny{$\chi^{2}$}}}(\pi\|\rho)}{N}.\end{split}$

###### Proof.

Given

$\pi(u)=\frac{1}{Z}g(u)\rho(u)=\frac{1}{\rho(g)}g(u)\rho(u),$

we have that
$d_{\mbox{\tiny{$\chi^{2}$}}}(\pi\|\rho)=\frac{\rho(g^{2})}{\rho(g)^{2}}-1.$
To ease the notation we denote $\zeta:=\frac{\rho(g^{2})}{\rho(g)^{2}}.$ We
rewrite

$\begin{split}\pi(f)&=\frac{\rho(gf)}{\rho(g)}\simeq\frac{\rho^{N}_{MC}(gf)}{\rho^{N}_{MC}(g)}=\pi^{N}_{IS}(f).\\\
\end{split}$

Then we have

$\begin{split}\pi^{N}_{IS}(f)-\pi(f)&=\pi^{N}_{IS}(f)-\frac{\rho(gf)}{\rho(g)}\\\
&=\frac{\pi^{N}_{IS}(f)\mathopen{}\mathclose{{}\left(\rho(g)-\rho^{N}_{MC}(g)}\right)}{\rho(g)}-\frac{\mathopen{}\mathclose{{}\left(\rho(gf)-\rho^{N}_{MC}(gf)}\right)}{\rho(g)}.\end{split}$
(5.6)

The expectation of the second term is zero and hence

$\begin{split}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\pi^{N}_{IS}(f)-\pi(f)}\right]}}\right|&=\frac{1}{\rho(g)}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\pi^{N}_{IS}(f)\mathopen{}\mathclose{{}\left(\rho(g)-\rho^{N}_{MC}(g)}\right)}\right]}}\right|\\\
&\leq\frac{1}{\rho(g)}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{IS}(f)-\pi(f)}\right)\mathopen{}\mathclose{{}\left(\rho(g)-\rho^{N}_{MC}(g)}\right)}\right]}}\right|,\\\
\end{split}$

since
$\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\rho(g)-\rho^{N}_{MC}(g)}\right]}=0$.
Using the Cauchy-Schwarz inequality, the second result from this theorem
(whose proof follows) and (5.4) from the proof of Theorem 5.1 we have, for all
$|f|_{\infty}\leq 1$,

$\begin{split}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\pi^{N}_{IS}(f)-\pi(f)}\right]}}\right|&\leq\frac{1}{\rho(g)}\mathopen{}\mathclose{{}\left(\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{IS}(f)-\pi(f)}\right)^{2}}\right]}}\right)^{1/2}\mathopen{}\mathclose{{}\left(\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\rho(g)-\rho^{N}_{MC}(g)}\right)^{2}}\right]}}\right)^{1/2}\\\
&\leq\frac{1}{\rho(g)}\mathopen{}\mathclose{{}\left(\frac{4\zeta}{N}}\right)^{1/2}\mathopen{}\mathclose{{}\left(\frac{\rho(g^{2})}{N}}\right)^{1/2}=\frac{2\zeta}{N}.\end{split}$

We now prove the second result. We use the splitting of
$\pi^{N}_{IS}(f)-\pi(f)$ into the sum of two terms as derived in equation 5.6.
Using (5.4), the basic inequality $(a-b)^{2}\leq 2(a^{2}+b^{2})$ and that for
all $|f|_{\infty}\leq 1,|\pi^{N}_{IS}(f)|\leq 1$, we have, for all
$|f|_{\infty}\leq 1$,

$\begin{split}&\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{IS}(f)-\pi(f)}\right)^{2}}\right]}}\right|\\\
&\leq\frac{2}{\rho(g)^{2}}\mathopen{}\mathclose{{}\left(\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{IS}(f)}\right)^{2}\mathopen{}\mathclose{{}\left(\rho(g)-\rho^{N}_{MC}(g)}\right)^{2}}\right]}+\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\rho(gf)-\rho^{N}_{MC}(gf)}\right)^{2}}\right]}}\right)\\\
&\leq\frac{2}{\rho(g)^{2}}\mathopen{}\mathclose{{}\left(\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\rho(g)-\rho^{N}_{MC}(g)}\right)^{2}}\right]}+\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\rho(gf)-\rho^{N}_{MC}(gf)}\right)^{2}}\right]}}\right)\\\
&=\frac{2}{\rho(g)^{2}N}\bigl{(}\text{Var}_{\rho}\mathopen{}\mathclose{{}\left[g}\right]+\text{Var}_{\rho}\mathopen{}\mathclose{{}\left[gf}\right]\bigr{)}\\\
&\leq\frac{2}{\rho(g)^{2}N}\mathopen{}\mathclose{{}\left(\rho(g^{2})+\rho(g^{2}f^{2})}\right)\\\
&\leq\frac{4\rho(g^{2})}{\rho(g)^{2}N}=\frac{4\zeta}{N}.\end{split}$

Therefore,

$\begin{split}\sup_{|f|_{\infty}\leq
1}\mathopen{}\mathclose{{}\left|\operatorname{\mathbb{E}}{\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\pi^{N}_{IS}(f)-\pi(f)}\right)^{2}}\right]}}\right|\leq\frac{4\zeta}{N}.\end{split}$

∎

This theorem shows that, unlike Monte Carlo, the Importance Sampling estimator
$\pi^{N}_{IS}$ is biased for $\pi$. The rate of convergence of the bias is
twice that of the standard deviation, however. As for Monte Carlo the rate of
convergence of the variance of the error is governed by the inverse of $N$,
but is independent of the dimension $d$ of the state space. However for
Importance Sampling to be accurate (with a limited number of samples $N$) it
is important that the $\chi^{2}$ divergence between the target and the
proposal is not too large.

######  Example 5.5 (Change of Measurement).

We consider a similar set-up as in Example 5.2, integrating a sigmoid
function, shown in blue in Figure 8, with respect to a probability measure
$\pi$ which is bimodal, shown in red in Figure 8; we again restrict the
support of the desired integral. We estimate the integral using Importance
Sampling based on $N$ random samples $u^{(1)},\ldots,u^{(N)}$ from the measure
$\rho=\mathcal{N}(\mu,\,\sigma^{2})$, shown in green in Figure 8. The
estimator of the integral is given by

$\displaystyle\begin{split}\pi^{N}_{IS}(f)&=\sum_{n=1}^{N}w^{(n)}f(u^{(n)}){\mathbb{1}}_{[a,b]}(u^{(n)}),\\\
w^{(n)}&=\frac{g(u^{(n)})}{\sum_{m=1}^{N}g(u^{(m)})}.\end{split}$

Here $g$ is a function proportional to the ratio of the densities of $\pi$ and
$\rho.$ If $\pi(u^{(n)})>\rho(u^{(n)})$, the samples should have been denser,
so we raise the weight on $f(u^{(n)})$ in proportion to
$\frac{\pi(u^{(n)})}{\rho(u^{(n)})}>1$. If $\pi(u^{(n)})<\rho(u^{(n)})$, the
samples should have been less dense, so we lower the weight on $f(u^{(n)})$ in
proportion to $\frac{\pi(u^{(n)})}{\rho(u^{(n)})}<1$.

---

---

Figure 8: Importance Sampling is a change of measurement via the importance
weights. The red curve shows a bimodal distribution $\pi$ and the green curve
shows a Gaussian distribution $\rho$. The blue curve is the function to be
integrated, on its support $[-5,5].$ The upper figure shows samples from the
posterior $\pi$ itself; these would be used for Monte Carlo sampling; the
lower curve shows samples from the prior $\rho$, as used for Importance
Sampling. The importance weights capture and compensate for the difference of
sampling from these two distributions.

####  5.3 Discussion and Bibliography

A classic reference on the Monte Carlo method which includes discussion of
several variance-reduction techniques is [52]. The chapter notes [7] give a
comparison of Monte Carlo and Importance Sampling with examples. The paper
[68] further explores advanced Importance Sampling via adaptive algorithms.
When $N$ is large enough, $\pi^{N}_{MC}$ arising from a Monte Carlo simulation
should be close to $\pi$ [36]. In practice, all probabilities, integrals and
summations can be approximated by the Monte Carlo method [110]. A review of
importance sampling, from the perspective of filtering and sequential
importance resampling, may be found in [4]; the proofs in this chapter closely
follow the presentation in that paper. Necessary sample size results for
importance sampling in terms of several divergences between target and
proposal were established in [100]. The subject of multi level Monte Carlo
(MLMC) has made the use of Monte Carlo methods practical in new areas of
application; see [47] for an overview. The methodology applies when
approximating expectations over infinite dimensional spaces, and distributes
the computational budget over different levels of approximation, with the goal
of optimizing the cost per unit error, noting that the latter balances
sampling and approximation based sources.

###  6 Markov Chain Monte Carlo

In this chapter we study Markov chain Monte Carlo (MCMC), a methodology that
allows to obtain approximate samples from a given target distribution $\pi.$
As with Monte Carlo and importance sampling, MCMC may be viewed as
approximating the target distribution by a sum of Dirac measures, thus
allowing to approximate expectations with respect to the target.
Implementation of Monte Carlo pressuposes that independent samples from the
target can be obtained. Importance sampling and MCMC bypass this restrictive
assumption: importance sampling by appropriately weighting independent samples
from a proposal distribution, and MCMC by drawing correlated samples from a
Markov kernel that has the target as invariant distribution. After some
discussion of the general MCMC methodology, we will specify to the case where
$\pi$ is a posterior distribution given via Bayes theorem from the product of
the likelihood function and the prior distribution. In this context, we will
analyze the convergence of the pCN algorithm, a popular MCMC method for
computing posterior expectations in high dimensional inverse problems which
uses these two ingredients, the prior and the likelihood, as part of its
design.

####  6.1 The Idea Behind MCMC

The idea of MCMC is that, given a target distribution $\pi,$ it is possible to
construct a Markov kernel that can be sampled from and has $\pi$ as its
invariant distribution (a formal definition of the invariance property will be
given below). Samples $\\{u^{(n)}\\}_{n=1}^{N}$ drawn iteratively from the
kernel may be used to approximate posterior expectations

$\pi(f)\approx\frac{1}{N}\sum_{n=1}^{N}f(u^{(n)}).$

The samples are given uniform weights $1/N$ but, in contrast to standard Monte
Carlo, they are not independent and they are not drawn from the target $\pi.$
However, if the chain is guaranteed to satisfy certain sample path ergodicity
then the estimator in the right hand-side is asymptotically unbiased for
$\pi(f)$ and satisfies a central limit theorem for suitable test functions
$f.$

Addressing the design and analysis of MCMC methods in generality and depth is
beyond the scope of a single chapter; entire books are devoted to this
subject. We will restrict our discussion to a particular class of MCMC