Random forward models and log-likelihoods in Bayesian inverse problems

Throughout, $(\Omega,\mathcal{F},\mathbb{P})$ is a fixed probability space that is rich enough to serve as a common domain for all random variables of interest.

The space of probability measures on the Borel $\sigma$-algebra of a topological space $\mathcal{U}$ will be denoted by $\mathcal{M}_{1}(\mathcal{U})$; in practice, $\mathcal{U}$ will be a separable Banach space.

When $\mu\in\mathcal{M}_{1}(\mathcal{U})$, integration of a measurable function (random variable) $f\colon\mathcal{U}\to\mathbb{R}$ will also be denoted by expectation, i.e. $\mathbb{E}_{\mu}[f]\coloneqq\int_{\mathcal{U}}f(u)\,\mathrm{d}\mu(u)$.

The space $\mathcal{M}_{1}(\mathcal{U})$ will be endowed with the Hellinger metric $d_{\textup{H}}\colon\mathcal{M}_{1}(\mathcal{U})^{2}\to\mathbb{R}_{\geq 0}$: for $\mu,\nu\in\mathcal{M}_{1}(\mathcal{U})$ that are both absolutely continuous with respect to a reference measure $\pi$,

The Hellinger distance is in fact independent of the choice of reference measure $\pi$ and defines a metric on $\mathcal{M}_{1}(\mathcal{U})$ (Bogachev, 2007, Lemma 4.7.35–36) with respect to which $\mathcal{M}_{1}(\mathcal{U})$ evidently has diameter at most $1$. The Hellinger topology coincides with the total variation topology (Kraft, 1955) and is strictly weaker than the Kullback–Leibler (relative entropy) topology (Pinsker, 1964); all these topologies are strictly stronger than the topology of weak convergence of measures.

As used in Sections 3–5, the Hellinger metric is useful for uncertainty quantification when assessing the similarity of Bayesian posterior probability distributions, since expected values of square-integrable functions are Lipschitz continuous with respect to the Hellinger metric:

when $f\in L^{2}_{\mu}(\mathcal{U})\cap L^{2}_{\nu}(\mathcal{U})$. In particular, for bounded $f$, $|\mathbb{E}_{\mu}[f]-\mathbb{E}_{\nu}[f]|\leq 2\sqrt{2}\|f\|_{\infty}d_{\textup{H}}(\mu,\nu)$.

By an inverse problem we mean the recovery of $u\in\mathcal{U}$ from an imperfect observation $y\in\mathcal{Y}$ of $G(u)$, for a known forward operator $G\colon\mathcal{U}\to\mathcal{Y}$. In practice, the operator $G$ may arise as the composition $G=O\circ S$ of the solution operator $S\colon\mathcal{U}\to\mathcal{V}$ of a system of ordinary or partial differential equations with an observation operator $O\colon\mathcal{V}\to\mathcal{Y}$, and it is typically the case that $\mathcal{Y}=\mathbb{R}^{J}$ for some $J\in\mathbb{N}$, whereas $\mathcal{U}$ and $\mathcal{V}$ can have infinite dimension. For simplicity, we assume an additive noise model

where the statistics but not the realisation of $\eta$ are known. In the strict sense, this inverse problem is ill-posed in the sense that there may be no element $u\in\mathcal{U}$ for which $G(u)=y$, or there may be multiple such $u$ that are highly sensitive to the observed data $y$.

The Bayesian perspective eases these problems by interpreting $u$, $y$, and $\eta$ all as random variables or fields. Through knowledge of the distribution of $\eta$, (2.3) defines the conditional distribution of $y|u$. After positing a prior probability distribution $\mu_{0}\in\mathcal{M}_{1}(\mathcal{U})$ for $u$, the Bayesian solution to the inverse problem is nothing other than the posterior distribution for the conditioned random variable $u|y$. This posterior measure, which we denote $\mu^{y}\in\mathcal{M}_{1}(\mathcal{U})$, is from the Bayesian point of view the proper synthesis of the prior information in $\mu_{0}$ with the observed data $y$. The same posterior $\mu^{y}$ can also be arrived at via the minimisation of penalised Kullback–Leibler, $\chi^{2}$, or Dirichlet energies (Dupuis and Ellis, 1997; Jordan and Kinderlehrer, 1996; Ohta and Takatsu, 2011), where the penalisation again expresses compromise between fidelity to the prior and fidelity to the data.

The rigorous formulation of Bayes’ formula for this context requires careful treatment and some further notation (Stuart, 2010). The pair $(u,y)$ is assumed to be a well-defined random variable with values in $\mathcal{U}\times\mathcal{Y}$. The marginal distribution of $u$ is the Bayesian prior $\mu_{0}\in\mathcal{M}_{1}(\mathcal{U})$. The observational noise $\eta$ is distributed according to $\mathbb{Q}_{0}\in\mathcal{M}_{1}(\mathcal{Y})$, independently of $u$. The random variable $y|u$ is distributed according to $\mathbb{Q}_{u}$, the translate of $\mathbb{Q}_{0}$ by $G(u)$, which is assumed to be absolutely continuous with respect to $\mathbb{Q}_{0}$, with

The function $\Phi\colon\mathcal{U}\times\mathcal{Y}\to\mathbb{R}$ is called the negative log-likelihood or simply potential. In the elementary setting of centred Gaussian noise, $\eta\sim\mathcal{N}(0,\Gamma)$ on $\mathcal{Y}=\mathbb{R}^{J}$, the potential is the non-negative quadratic misfit⁶⁶6Hereafter, to reduce notational clutter, we write both $\|\hbox to5.71527pt{\hss$\cdot$\hss}\|_{\mathcal{U}}$ and $\|\hbox to5.71527pt{\hss$\cdot$\hss}\|_{\mathcal{Y}}$ as $\|\hbox to5.71527pt{\hss$\cdot$\hss}\|$. $\Phi(u;y)=\tfrac{1}{2}\bigl{\|}\Gamma^{-1/2}(y-G(u))\bigr{\|}_{\mathcal{Y}}^{2}$. However, particularly for cases in which $\dim\mathcal{Y}=\infty$, it may be necessary to allow $\Phi$ to take negative values and even to be unbounded below (Stuart, 2010, Remark 3.8).

With this notation, Bayes’ theorem is then as follows (Dashti and Stuart, 2016, Theorem 3.4):

Note that, for (2.4) to make sense, it is essential to check that $0<Z(y)<\infty$. Hereafter, to save space, we regard the data $y$ as fixed, and hence write $\Phi(u)$ in place of $\Phi(u;y)$, $Z$ in place of $Z(y)$, and $\mu$ in place of $\mu^{y}$. In particular, we shall redefine the negative log-likelihood as a function $\Phi\colon\mathcal{U}\to\mathbb{R}$, instead of a function $\Phi\colon\mathcal{U}\times\mathcal{Y}\to\mathbb{R}$ as in Theorem 2.1 above.

From the perspective of numerical analysis, it is natural to ask about the well-posedness of the Bayesian posterior $\mu$: is it stable when the prior $\mu_{0}$, the potential $\Phi$, or the observed data $y$ are slightly perturbed, e.g. due to discretisation, truncation, or other numerical errors? For example, what is the impact of using an approximate numerical forward operator $G_{N}$ in place of $G$, and hence an approximate $\Phi_{N}\colon\mathcal{U}\to\mathbb{R}$ in place of $\Phi$? Here, we quantify stability in the Hellinger metric $d_{\textup{H}}$ from (2.1).

Stability of the posterior with respect to the observed data $y$ and the log-likelihood $\Phi$ was established for Gaussian priors by Stuart (2010) and for more general priors by many later contributions (Dashti et al., 2012; Hosseini, 2017; Hosseini and Nigam, 2017; Sullivan, 2017). (We note in passing that the stability of BIPs with respect to perturbation of the prior is possible but much harder to establish, particularly when the data $y$ are highly informative and the normalisation constant $Z(y)$ is close to zero; see e.g. the “brittleness” phenomenon of (Owhadi and Scovel, 2017; Owhadi et al., 2015).) Typical approximation theorems for the replacement of the potential $\Phi$ by a deterministic approximate potential $\Phi_{N}$, leading to an approximate posterior $\mu_{N}$, aim to transfer the convergence rate of the forward problem to the inverse problem, i.e. to prove an implication of the form

where $M\colon\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ is suitably well behaved, $\psi\colon\mathbb{N}\to\mathbb{R}_{\geq 0}$ quantifies the convergence rate of the forward problem, and $C$ is a constant. Following Stuart and Teckentrup (2017), the purpose of this article is to extend this paradigm and these approximation results to the case in which the approximation $\Phi_{N}$ is a random object.

This section considers the general setting in which the deterministic potential $\Phi$ is approximated by a random potential $\Phi_{N}\sim\nu_{N}$. Recall from (2.4) that $Z$ is the normalisation constant of $\mu$, and that for $\mu$ to be well-defined, we must have that $0<Z<\infty$. The following two results, Theorems 3.1 and 3.2, extend Theorems 4.9 and 4.11 respectively of Stuart and Teckentrup (2017), in which the approximation is a GP model:

In the proof of Theorem 3.1, we show that hypothesis (a) arises as an upper bound on the quantity $\|(e^{-\Phi}+\mathbb{E}_{\nu_{N}}[e^{-\Phi_{N}}])^{-1}\|_{L^{p_{1}}_{\mu_{0}}(\mathcal{U})}$. In order for the conclusion of Theorem 3.1 to hold, we need the latter to be finite. Thus, hypothesis (a) is an exponential decay condition on the positive tails of either $\Phi$ or $\Phi_{N}$, with respect to the appropriate measures. Alternatively, by applying Jensen’s inequality to $\mathbb{E}_{\nu_{N}}[e^{-\Phi_{N}}]^{-1}$, one can strengthen hypothesis (a) into the hypothesis of exponential integrability of either $\Phi$ with respect to $\mu_{0}$ or $\Phi_{N}$ with respect to $\nu_{N}\otimes\mu_{0}$; this yields the same interpretation. Thus, the parameter $p_{1}$ quantifies the exponential decay of the positive tail of either $\Phi$ or $\Phi_{N}$.

By comparing the quantity $\|(e^{-\Phi}+\mathbb{E}_{\nu_{N}}[e^{-\Phi_{N}}])^{-1}\|_{L^{p_{1}}_{\mu_{0}}(\mathcal{U})}$ from hypothesis (a) with the quantity in hypothesis (b), it follows that hypothesis (b) is an exponential decay condition on the negative tails of both $\Phi$ and $\Phi_{N}$. The two new parameters in this decay condition arise because we apply Hölder’s inequality twice in order to develop the desired bound (3.3) on $d_{\textup{H}}(\mu,\mu_{N}^{\textup{M}})$. The key desideratum here is that the bound is multiplicative in some $L^{p^{\prime}}_{\mu_{0}}(\mathcal{U})$-norm of $\mathbb{E}_{\nu_{N}}[|\Phi-\Phi_{N}|^{p_{2}^{\prime}}]^{1/p_{2}^{\prime}}$. The two new parameters $p_{2}$ and $p_{1}^{\prime}p_{3}$ quantify the decay with respect to $\nu_{N}$ and $\mu_{0}$ respectively. Note that the interaction between the hypotheses (a) and (b) as described by the conjugate exponent pair $(p_{1},p_{1}^{\prime})$ implies that one can trade off faster exponential decay of one tail with slower exponential decay of the other.

The two-sided condition on $\mathbb{E}_{\nu_{N}}[Z_{N}^{\textup{S}}]$ in hypothesis (c) ensures that both tails of $\Phi_{N}$ with respect to $\nu_{N}\otimes\mu_{0}$ decay sufficiently quickly. This hypothesis ensures that the Radon–Nikodym derivative in (3.2) is well-defined.

Finally, we note that the quantity on the right hand side of (3.3a) depends directly on the conjugate exponents of $p_{1},p_{2}$ and $p_{3}$ appearing in hypotheses (a) and (b). The more well behaved the quantities in these hypotheses are, the weaker the norm we can choose on the right hand side of (3.3a).

Hypothesis (a) of Theorem 3.2 arises during the proof as a result of developing an upper bound on $\|\mathbb{E}_{\nu_{N}}[(e^{-\Phi/2}-e^{-\Phi_{N}/2})^{2}]\|$ that is multiplicative in some $L^{p^{\prime}}_{\mu_{0}}(\mathcal{U})$-norm of $\|\mathbb{E}_{\nu_{N}}[|\Phi-\Phi_{N}|^{2q_{1}^{\prime}}]^{1/2q_{1}^{\prime}}\|$. Thus, it describes an exponential decay condition of the negative tails of both $\Phi$ or $\Phi_{N}$; in particular, hypothesis (a) is always satisfied when the potentials $\Phi$ or $\Phi_{N}$ are non-negative, as is usually the case for finite-dimensional data. The appearance of $q_{1}$ and $q_{2}$ arises due to one application of Hölder’s inequality for fulfulling the desideratum of multiplicativity, and $q_{1}$ and $q_{2}$ quantify the decay with respect to $\nu_{N}$ and $\mu_{0}$ respectively.

Hypothesis (b) of Theorem 3.2 arises as a result of developing an upper bound on the quantity $\mathbb{E}_{\nu_{N}}[Z_{N}^{\textup{S}}(Z^{-1/2}-(Z_{N}^{\textup{S}})^{-1/2})^{2}]$ that fulfills the desideratum of multiplicativity mentioned above. The presence of both $Z_{N}^{\textup{S}}$ and its reciprocal indicates that hypothesis (b) is analogous to hypothesis (c) of Theorem 3.1, in that hypothesis (b) is a condition on the tails of $\Phi_{N}$ with respect to $\mu_{0}$. The difference between hypothesis (b) of Theorem 3.2 and hypothesis (c) of Theorem 3.1 arises due to the fact that the Radon–Nikodym derivative in (3.1) features $Z_{N}^{\textup{S}}$ instead of $\mathbb{E}_{\nu_{N}}[Z_{N}^{\textup{S}}]$.

We now show that the assumptions of Theorems 3.1 and 3.2 are satisfied when the exact potential $\Phi$ and the approximation quality $\Phi_{N}\approx\Phi$ are suitably well behaved. Since $0<Z<\infty$, it follows that $C_{3}^{-1}<Z<C_{3}$ for some $0<C_{3}<\infty$.

The lower bound conditions in (3.5) ensure that the hypothesised exponential decay conditions on the negative tails of the true likelihood and the random likelihoods from Theorems 3.1 and 3.2 are satisfied. The uniform lower bound on $\Phi$ translates into a uniform upper bound of the Radon–Nikodym derivative of the posterior with respect to the prior, and is a very mild condition that is satisfied in many, if not most, BIPs. Given this fact, it is reasonable to demand that the $\Phi_{N}$ satisfy the same uniform lower bound, $\nu_{N}$-almost surely and for all $N\in\mathbb{N}$; this is the content of the second condition in (3.5). Condition (3.6) expresses the condition that, by choosing $N$ sufficiently large, one can approximate $\Phi$ arbitrarily well using the random $\Phi_{N}$, with respect to the $L^{1}_{\mu_{0}\otimes\nu_{N}}$ topology. This assumption ensures that the stated aims of this work are reasonable.

The uniform upper bound condition on $\Phi_{N}$ with respect to $\mu_{0}$ in (3.7) is rather strong; we use it to ensure that $Z_{N}^{\textup{S}}$ is bounded away from zero, uniformly with respect to $\Phi_{N}$ and $N\in\mathbb{N}$. Together with the condition on $\Phi_{N}$ in (3.5), this translates to uniform lower and upper bounds on $Z_{N}^{\textup{S}}$; the latter implies that hypothesis (b) in Theorem 3.2 holds with the stated values of $q_{1}$ and $q_{2}$. A sufficient condition for (3.7) is that the $\Phi_{N}$ are themselves uniformly bounded. This condition is of interest when the misfit $\Phi$ is associated to a bounded forward model and the data take values in a bounded subset.

By comparing the hypotheses and conclusions of Lemma 3.4 and Lemma 3.5, we observe that, by reducing the exponent of integrability from $q_{1}=+\infty$ to $q_{1}=\rho^{\ast}/2$, we can replace the strong uniform upper bound condition (3.7) on $\Phi_{N}$ from Lemma 3.4 with the weaker condition that $\exp(\Phi_{N})\in L^{\rho^{\ast}}_{\mu_{0}}(\mathcal{U})$ in Lemma 3.5, and thus increase the scope of applicability of the conclusion.

In Lemmas 3.4 and 3.5 above, we have specified the largest possible values of the exponents that are compatible with the hypotheses. This is because later, in Theorem 3.9, we will want to use the smallest possible values of the corresponding conjugate exponents in the resulting inequalities (3.3a) and (3.4).

In many settings, the potentials $\Phi$ and $\Phi_{N}$ have a common form and differ only in the parameter-to-observable map. In this section we shall assume that $\Phi$ and $\Phi_{N}$ are quadratic misfits of the form

corresponding to centred Gaussian observational noise with symmetric positive-definite covariance $\Gamma$. Again, we assume that $G$ is deterministic while $G_{N}$ is random. In this section, for this setting, we show how the quality of the approximation $G_{N}\approx G$ transfers to the approximation $\Phi_{N}\approx\Phi$, and hence to the approximation $\mu_{N}\approx\mu$ (for either the sample or marginal approximate posterior).

Pointwise in $u$ and $\omega$, the errors in the misfit and the forward model are related according to the following proposition.

By assuming that $\mathcal{Y}=\mathbb{R}^{J}$, we assume that the data live in a finite-dimensional space. This is a standard assumption in the area, and implies that the operator $\Gamma$ is simply a matrix. The assumption of the eigenvalues of $\Gamma$ being bounded away from zero is equivalent to assuming that $\Gamma$ is invertible, which follows immediately from the assumption stated earlier that $\Gamma$ is a symmetric and positive-definite covariance matrix.

The hypotheses ensure that the integrability of the misfit $\Phi$ determines the highest degree of integrability of the forward operators $G_{N}$ and $G$, and that for sufficiently large $N$, we may make the norm of the difference of $G-G_{N}$ in an appropriate topology small enough. The constraint (3.7) is used to combine the $\|G(u)-G_{N}(u)\|$ and $\|G(u)-G_{N}(u)\|^{2}$ terms in (3.9). The resulting simplification ensures that we may apply Lemma 3.8.

The lemma states that if the random forward model converges to the true forward model in the appropriate topology, then the conditions in Assumption 3.3 are satisfied by the corresponding random misfits. Since the misfits were assumed to be quadratic in (3.8), the key contribution of Lemma 3.8 is to ensure that the approximation quality condition (3.6) is satisfied.

We shall use the preceding results to obtain bounds on the Hellinger distance in terms of errors in the forward model, of the following form: for $C,D>0$ and $r_{1},r_{2},s_{1},s_{2}\geq 1$ that do not depend on $N$,

For brevity and simplicity, the following result uses one pair $q,s\geq 1$ in (3.11) in order to obtain convergence statements for both $\mu_{N}^{\textup{M}}$ and $\mu^{\textup{S}}_{N}$. If one is interested in only one of these measures, then one may optimise $q$ and $s$ accordingly.

The proof of Theorem 3.9 consists of tracking the dependence of the parameters over the sequential application of the preceding results, all of which are used.

Case (a) applies in the situation where the random approximations $\Phi_{N}$ are uniformly bounded from above; as discussed earlier, this condition is satisfied in the case that the misfit $\Phi$ is associated to a bounded forward model and the data take values in a bounded subset of $\mathcal{Y}=\mathbb{R}^{J}$. Note that the topology of the convergence of $G_{N}$ to $G$ is quantified by $s$ and $q$, and that $s$ depends on the parameter $p^{\ast}$ that quantifies the exponential $\mu_{0}$-integrability of the misfit $\Phi$. In particular, the faster the exponential decay of the positive tail of $\Phi$ (i.e. the larger the value of $p^{\ast}$), the stronger the topology of convergence of $G_{N}$ to $G$.

In contrast to case (a), case (b) does not assume that the misfit $\Phi$ is exponentially integrable or that the random approximations $\Phi_{N}$ are uniformly bounded from above $\nu_{N}$-almost surely. Instead, exponential integrability of the random misfit $\Phi_{N}$ is required. Another difference is that the exponential integrability parameter $\rho^{\ast}$ determines the strength of the topology of convergence of the random forward models, not only with respect to the $\mu_{0}$-topology, but also to the $\nu_{N}$-topology as well.