Past exam of the mathematics course of the University of Cambridge / 2019 / iii / Paper 350 (Bayesian Inverse Problems)

Let ${\pi }_0$ be the probability density function of the unknown and $\rho$ that of the noise, both relative to Lebesgue measure. By independence of random variables, the conditional distribution of the measurement given $u$ has density $\rho (m-Au)$. Thus its likelihood function and the joint density areThe Bayesian model evidence is the marginal data densityBayes' theorem gives the posterior density solving this Bayesian inverse problem:This formula applies when $0<Z(m)<\infty$. Since $\int Z(m)dm=1$ by Tonelli theorem, these conditions hold almost everywhere under the marginal data law. Values of a conditional distribution at exceptional data are not determined by the joint law. The posterior mean or a maximum a posteriori estimate can then give a point estimate, while the full posterior density describes the remaining uncertainty.

Write $a=(2,1{)}^T$, so the data map is $a^{T}u$. The Gaussian likelihood and prior distribution giveCompleting the square, or using Gaussian conjugacy for a normal linear model, identifies a multivariate normal distribution with precision matrix $C_{\delta }^{-1}=I_2+{\delta }^{-2}a{a}^T$. Since $a^{T}a=5$, its posterior mean and covariance matrix areThe parameters areFor example, $(I_2+{\delta }^{-2}a{a}^T)C_{\delta }=I_2$, which checks the square completion.

The orthogonal vectorsare eigenvectors of the covariance matrix. Their respective variances areThey are independent under the posterior distribution, by independence of uncorrelated jointly normal variables. Thus uncertainty is different in the two directions even for positive noise: the data constrain $e_{\parallel }$, whereas $e_{\perp }$ spans the null space of the data map.

As $\delta \to 0$, the posterior distribution converges to the law ofThis is supported on $2u_1+u_2=m$, and its posterior mean is the minimum-norm least-squares solution. The prior distribution still determines the distribution along the unobserved null space; it has not disappeared. This example of a Gaussian posterior in the zero-noise limit distinguishes exact recovery of an observed component from recovery of the whole unknown.

Use the Hellinger distance convention in the original PDF:The TeX transcription incorrectly places the integral outside the outer square root in the definition; the displayed formula for the squared distance and the original PDF both use the convention above. Although Lebesgue measure is not a probability measure, we may use it as the common dominating measure: the integral is unchanged by replacing it with any equivalent dominating probability measure.

For unit-variance normal distributions, put $\bar{\theta }=({\theta }_1+{\theta }_2)/2$ and $\Delta ={\theta }_1-{\theta }_2$. The product of the square roots of their probability density functions isThe final factor integrates to one, soFor the bound at arbitrary ${\sigma }_1,{\sigma }_2>0$, use the given formula for the Hellinger distance between normal distributions and writeSince $1-e^{-x}\le x$ for $x\ge 0$,Therefore one valid choice is $\boxed{C_{{\sigma }_1,{\sigma }_2}=1/({\sigma }_1^2+{\sigma }_2^2)}$. In particular this bound is uniform when the two standard deviations stay bounded away from zero.

The conditional mean estimator is the posterior meanprovided the first moment is finite. A maximum a posteriori estimate is any maximizer of the posterior density relative to Lebesgue measure:Existence of the posterior density alone does not guarantee a finite posterior mean or existence of a maximizing point.

For the stated multivariate normal distribution and independent standard Gaussian noise, the negative log of the posterior density, up to a constant, isIts precision matrix $A^{T}A+{\Sigma }^{-1}$ is a positive-definite matrix, so Gaussian conjugacy for a normal linear model givesThe two forms agree becauseThe covariance matrix is $(A^{T}A+{\Sigma }^{-1}{)}^{-1}$, and the maximum a posteriori estimate equals the posterior mean for this multivariate normal distribution.

A simple choice is an independent Laplace distribution at each pixel, giving a Laplace prior for sparse images:Its negative log density penalizes the L1 norm. Relative to a normal distribution, it combines a sharp peak at zero with heavier tails, allowing most pixels to lie near the background level while a few have appreciable intensity. Quantitatively, $\mathbb{P}(|u_j|>\tau )=e^{-\lambda \tau }$, so the expected number of appreciable pixels is $d{e}^{-\lambda \tau }$. Choosing $\lambda \tau$ large makes this small. This is a simple sparsity model: it does not enforce that nearby active pixels form connected objects, nor does its continuous density give exact zero pixels positive probability. Additional spatial modelling would be needed to insist on contiguous shapes.

To sample, draw independent $U_j\sim \mathcal{U}(0,1)$ and use inverse transform sampling:The endpoint events have probability zero and can be assigned arbitrary finite values. To justify the rule, the cumulative distribution function of a centered Laplace distribution isThe displayed rule is $F^{-1}(U_j)$, so $\mathbb{P}(u_j\le x)=\mathbb{P}(U_j\le F(x))=F(x)$. Independence of random variables then gives the stated product prior distribution.

Use the shape-rate convention for the Gamma distribution. DefineThe zero null space assumption makes $B$ a positive-definite matrix, and $u_{\ast }$ is the unique ordinary least squares solution. The Gaussian likelihood contributes ${\gamma }^{k/2}$, which must be retained because $\gamma$ is unknown. Multiplying it by the improper prior and the Gamma distribution gives the normal-gamma posterior with a flat prior:Here the constant density of the improper prior cancels on normalization. PutThe normal equation $A^{T}r=0$ makes the residual orthogonal to the range of $A$. The orthogonal decomposition by a closed subspace into that range and its orthogonal complement givesIntegration over $u$ by the Gaussian integral, followed by integration over $\gamma$ using the Gamma distribution, yieldsThus the posterior distribution is proper, despite the improper prior; here $k\ge d$ follows from the zero null space assumption.

Completing the square in $u$ and collecting the powers of $\gamma$ give the two conditional distributions:Explicitly their probability density functions areFor a joint maximum a posteriori estimate relative to $dud\gamma$, first fix $\gamma >0$. The unique maximizing $u$ is $u_{\ast }$. The remaining log density is $(a-1)\log \gamma -b_0\gamma$ plus a constant. Consequently, if $a>1$,The stated assumptions do not always imply $a>1$. If $a=1$, the supremum occurs as $\gamma \downarrow 0$ but is not attained on $\gamma >0$; if $a<1$, the density is unbounded there. In either case there is no joint maximum a posteriori estimate in the stated parameter space. Calling zero a boundary mode does not make it an admissible positive precision.

If “MAP estimators” instead means separate modes of the marginal distributions, integration givesThus the marginal maximum a posteriori estimate of $u$ is always $u_{\ast }$, and that of $\gamma$ is $(s-1)/b_0$ when $s>1$, with the same nonattainment at zero when $s\le 1$. Joint and marginal maximization are different operations; the usual joint interpretation gives the preceding boxed pair.

The Cameron-Martin space of a Gaussian measure $\mu =\mathcal{N}(0,\Sigma )$ isThe inverse is defined on the range of ${\Sigma }^{1/2}$. In this infinite-dimensional setting “positive definite” must mean $⟨\Sigma h,h⟩>0$ for every nonzero $h$, rather than a uniform lower bound: a trace-class operator cannot be uniformly positive on an infinite-dimensional Hilbert space.

By the spectral theorem for compact Hermitian operators, choose an orthonormal basis $(e_j)$ with $\Sigma e_j={\lambda }_j{e}_j$, ${\lambda }_j>0$, and ${\sum }_j{\lambda }_j<\infty$. (A strictly positive trace-class operator also forces the ambient Hilbert space to be separable.) In coordinates $h_j=⟨h,e_j⟩$,This is a Hilbert space with its indicated norm, even though its range need not be closed in the ambient norm.

The Cameron-Martin theorem for a Gaussian measure states that the translated law ${\mu }_h(B)=\mu (B-h)=\mathcal{N}(h,\Sigma )(B)$ is an equivalent probability measure to $\mu$ exactly when $h\in E$. For such $h$, its Radon-Nikodym derivative isThe series has mean-square convergence and converges almost surely under $\mu$, since its independent summands have total variance ${\sum }_j{h}_j^2/{\lambda }_j$. Its law is $\mathcal{N}(0,\Vert h{\Vert }_E^2)$, so the density has expected value one. It is sometimes formally written ${\ell }_h(x)=⟨h,x{⟩}_E$, but a typical infinite-dimensional Gaussian sample is not in $E$; the series interpretation is essential. If $h\notin E$, the two laws are mutually singular measures.

The notation $P_0=\mathcal{N}(0,I)$ denotes generalized Gaussian white noise over $H=L^2({\mathbb{R}}^2)$, not a Gaussian measure supported on $H$: the identity is not a trace-class operator there. For example, realize the observation on a suitable space of tempered distributions. Its Cameron-Martin space of a Gaussian measure is $H$. For deterministic $h\in H$, the stochastic pairing $W_m(h)$ has normal distribution $\mathcal{N}(0,\Vert h{\Vert }_H^2)$ under $P_0$ and is the isonormal Gaussian process indexed by $H$.

Interpret the stated maps between Sobolev spaces as bounded linear operators, as usual. In particular $K:H\to H^2({\mathbb{R}}^2)\subset H$ is bounded, so its adjoint operator $K^{\ast }$ is bounded on $H$. Hence $A=K^{\ast }K$ maps $H$ into $H$, which is the property needed here. The stronger smoothing follows by Sobolev duality: boundedness of $K:H^{-4}\to H^{-2}$ gives $K^{\ast }:H^2\to H^4$, and therefore $A:H\to H^4\subset H$. As $\Pi (H)=1$, the shift $Au$ belongs to the noise Cameron-Martin space of a Gaussian measure for $\Pi$-almost every $u$.

The white-noise likelihood for a square-integrable shift, obtained from the Cameron-Martin theorem for a Gaussian measure in its white-noise form, isHere $P_u$ is the law of $Au+\eta$. In an orthonormal basis $(e_j)$ of $H$, with white-noise coordinates $m_j$, the pairing under $P_0$ isFor fixed $u$ this has mean-square convergence and converges almost surely, since ${\sum }_j|⟨Au,e_j⟩{|}^2=\Vert Au{\Vert }_H^2$. The stipulated joint measurability allows the likelihood to be used under ${\nu }_0(du,dm)=\Pi (du)P_0(dm)$. One must not replace this expression by a finite $\Vert m-Au{\Vert }_H^2$, because white noise is not $H$-valued. On the unbounded domain ${\mathbb{R}}^2$, membership of white noise in a global unweighted negative Sobolev space must not be assumed either; the stochastic pairing avoids that issue.

Let $\nu (du,dm)=\Pi (du)P_u(dm)$ be the actual joint law. The shift formula gives $\nu \ll {\nu }_0$ with Radon-Nikodym derivative $L$. DefineFor each fixed admissible $u$, the moment-generating function of a normal distribution givesBy Tonelli theorem, $\int Z(m)P_0(dm)=1$, so $Z(m)<\infty$ for $P_0$-almost every $m$. The likelihood is finite and strictly positive for ${\nu }_0$-almost every $(u,m)$; Fubini's theorem then gives $Z(m)>0$ for $P_0$-almost every $m$. The marginal observation law is $Z(m)P_0(dm)$ and is an equivalent probability measure to $P_0$.

The Bayes formula for a dominated observation model therefore definesTo verify that it is the conditional distribution, for measurable sets $B$ of unknowns and $C$ of data,Thus the posterior distribution is well defined for almost every observation under the actual data law. This is the appropriate almost-everywhere assertion supplied by Bayes' theorem; the stated assumptions do not prescribe a canonical posterior at every exceptional generalized datum.

Choose the common dominating probability measure $\nu =(\mu +{\mu }^{\prime })/2$ and let $p=d\mu /d\nu$, $q=d{\mu }^{\prime }/d\nu$ be the Radon-Nikodym derivatives. The assumed second moments imply first-moment integrability by the Cauchy-Schwarz inequality. Since $Y$ is a separable Banach space, the measurable $f$ is a strongly measurable function, and both expected values exist as Bochner integrals.

The difference of these Bochner integrals satisfiesThe first inequality is the norm bound for a Bochner integral, and the second is the scalar Cauchy-Schwarz inequality. Since $(\sqrt{p}+\sqrt{q}{)}^2\le 2(p+q)$, the first factor is at mostWith the Hellinger distance convention of the original PDF, the second factor is $\sqrt{2}{d}_{\mathrm{Hell}}(\mu ,{\mu }^{\prime })$. Multiplication proves the Hellinger bound for differences of expectations:The proof uses only the two laws and their common dominating measure; no absolute continuity of measures assumption between $\mu$ and ${\mu }^{\prime }$ is needed.