Reconciling Bayesian and perimeter regularization for binary inversion

TL;DR

This paper explores the connection between Bayesian inversion methods and perimeter regularization, demonstrating how Bayesian approaches can incorporate perimeter concepts for reconstructing discontinuous functions with uncertainty quantification.

Contribution

It establishes a theoretical link between Bayesian phase-field and level set methods and perimeter regularization, enhancing understanding of their regularization properties.

Findings

MAP objective function relates to least squares plus perimeter regularization

Bayesian level set samples have finite perimeter

Bayesian methods can learn about the true perimeter

Abstract

A central theme in classical algorithms for the reconstruction of discontinuous functions from observational data is perimeter regularization via the use of the total variation. On the other hand, sparse or noisy data often demands a probabilistic approach to the reconstruction of images, to enable uncertainty quantification; the Bayesian approach to inversion, which itself introduces a form of regularization, is a natural framework in which to carry this out. In this paper the link between Bayesian inversion methods and perimeter regularization is explored. In this paper two links are studied: (i) the maximum a posteriori (MAP) objective function of a suitably chosen Bayesian phase-field approach is shown to be closely related to a least squares plus perimeter regularization objective; (ii) sample paths of a suitably chosen Bayesian level set formulation are shown to possess finite perimeter and to have the ability to learn about the true perimeter.

Figures

Equations

$$B{V}_{\mathrm{binary}}(D)=\{\psi \in BV(D):\psi (x)\subset \{\pm 1\},x\in D\}$$

$$y=Ku+{\epsilon }^c\eta .$$

$$J(u):=\frac{1}{2\varepsilon^{2c}}\Bigl{|}\Sigma^{-\frac{1}{2}}(y-Ku)\Bigr{|}^{2}$$

$$\inf_{u\in\mathcal{P}}\Big{(}J(u)+\mathcal{R}(u)\Big{)}.$$

$$\inf_{u\in BV_{\rm binary}(D)}\Big{(}J(u)+\sigma\int_{D}|\nabla u|\Big{)}$$

$$\inf_{W^{1,p}(D)}\Big{(}J(u)+\sigma_{W}\int_{D}\Big{(}\tilde{\varepsilon}|\nabla u|^{2}+\frac{1}{\tilde{\varepsilon}}W(u)\Big{)}dx\Big{)}$$

$$\delta {\epsilon }^{-2a_1}{\Delta }^2\eta -q\delta {\epsilon }^{-2a_2}\Delta \eta +{\tau }^2\delta {\epsilon }^{-2a_3}\eta =f.$$

$$\|\xi\|_{\mathcal{E}_{\varepsilon,\tilde{\delta},\tilde{a}}}^{2}:=\langle(\mathcal{C}_{\varepsilon,\tilde{\delta},\tilde{a}})^{-\frac{1}{2}}\xi,(\mathcal{C}_{\varepsilon,\tilde{\delta},\tilde{a}})^{-\frac{1}{2}}\xi\rangle=\delta\int_{D}\bigl{(}\varepsilon^{-2a_{1}}|\triangle\xi|^{2}+q\varepsilon^{-2a_{2}}|\nabla\xi|^{2}+\tau^{2}\varepsilon^{-2a_{3}}\xi^{2}\bigr{)}\mathrm{d}x$$

$$u=\text{mathds}{1}_{D_{+}}-\text{mathds}{1}_{D_{+}^c}$$

$$S(v)=1,v>0,S(0)=0\mathrm{and}S(v)=-1,v<0.$$

$${\mathcal{E}}^S:=\{u=S(v)|v\in \mathcal{E}\}.$$

$$y=KS(v)+{\epsilon }^c\eta$$

$$\exp(-J(S(v)))=\exp\left(-\frac{1}{2\varepsilon^{2c}}\Bigl{|}\Sigma^{-\frac{1}{2}}\bigl{(}y-KS(v)\bigr{)}\Bigr{|}^{2}\right).$$

$$\frac{\mathrm{d}\mu^{y}}{\mathrm{d}\mu_{0,\alpha}}=\frac{1}{Z}\exp\left(-\frac{1}{2\varepsilon^{2c}}\Bigl{|}\Sigma^{-\frac{1}{2}}\bigl{(}y-KS(v)\bigr{)}\Bigr{|}^{2}\right);$$

$$d_{\mathrm{hell}}(\mu ,{\mu }^{\prime })=\sqrt{\left({\int }_X{\left(\sqrt{\frac{\mathrm{d}\mu }{\mathrm{d}{\mu }_0}}-\sqrt{\frac{\mathrm{d}{\mu }^{\prime }}{\mathrm{d}{\mu }_0}}\right)}^2\mathrm{d}{\mu }_0\right)}.$$

$$d_{\mathrm{hell}}({\mu }^y,{\mu }^{y^{\prime }})\le C(\rho )|y-y^{\prime }|;$$

$$\Psi(u)=\frac{r}{\varepsilon^{b}}\int_{D}\frac{1}{4}\bigl{(}1-u(x)^{2}\bigr{)}^{2}\,\mathrm{d}x.$$

$$\frac{\mathrm{d}{\nu }_0}{\mathrm{d}{\mu }_0}=\frac{1}{Z_0}\exp \left(-\Psi (u)\right).$$

$$\exp(-J(u))=\exp\left(-\frac{1}{2\varepsilon^{2c}}\Bigl{|}\Sigma^{-\frac{1}{2}}\bigl{(}y-Ku)\bigr{)}\Bigr{|}^{2}\right).$$

$$\frac{\mathrm{d}\nu^{y}}{\mathrm{d}\nu_{0}}=\frac{1}{Z}\exp\left(-\frac{1}{2\varepsilon^{2c}}\Bigl{|}\Sigma^{-\frac{1}{2}}(y-Ku)\Bigr{|}^{2}\right),$$

$$d_{\mathrm{hell}}({\nu }^y,{\nu }^{y^{\prime }})\le C(\rho )|y-y^{\prime }|.$$

$$\frac{\mathrm{d}\nu^{y}}{\mathrm{d}\mu_{0}}=\frac{1}{Z_{0}Z}\exp\left(-\frac{r}{\varepsilon^{b}}\int_{D}\frac{1}{4}\bigl{(}1-u(x)^{2}\bigr{)}^{2}\,\mathrm{d}x-\frac{1}{2\varepsilon^{2c}}\Bigl{|}\Sigma^{-\frac{1}{2}}\bigl{(}y-Ku\bigr{)}\Bigr{|}^{2}\right),$$

$$m_y=C{K}^{\ast }({\epsilon }^{2c}\Sigma +KC{K}^{\ast }{)}^{-1}y,C_y=C-C{K}^{\ast }({\epsilon }^{2c}\Sigma +KC{K}^{\ast }{)}^{-1}KC.$$

$$J(u)+\mathsf{R}(u)=\frac{1}{2\varepsilon^{2c}}\Big{|}\Sigma^{-\frac{1}{2}}(y-Ku)\Big{|}^{2}+\frac{1}{2}\|u\|_{\mathcal{E}}^{2}.$$

$$\underset{\rho \to 0}{\lim }\frac{{\nu }^y(B^{\rho }(\bar{z}))}{{\nu }^y(B^{\rho }(z^{\rho }))}=1$$

$$z^{\rho }=\underset{z\in X}{\mathrm{argmax}}{\nu }^y(B^{\rho }(z))$$

$$\Phi (u;y)=\Psi (u)+\frac{1}{2{\epsilon }^{2c}}|{\Sigma }^{-\frac{1}{2}}(y-Ku){|}^2,$$

$$J^{\varepsilon}(u)=\left\{\begin{array}[]{rl}\frac{1}{2}\|u\|_{\mathcal{E}}^{2}+\Phi(u;y)&\text{if }u\in\mathcal{E},\\ \infty&\text{if }u\notin\mathcal{E}.\end{array}\right.$$

$$J^{\epsilon }(u)={\epsilon }^{-2a_1-3}{I}^{\epsilon }(u),$$

$$I^{\epsilon }(u)=$$

Full text

{AMS}

35J35, 62G08, 62M40, 94A08.

1 Introduction

1.1 Problem Statement

Let $D$ be the unit cube $(0,1)^{d}\subset\mathbb{R}^{d},d=2,3$. Let K:L^{1}\bigl{(}D\bigr{)}\to\mathbb{R}^{J} be a bounded linear operator. We consider the problem of recovering a binary-valued function $u\in BV_{\rm binary}(D)$ where

[TABLE]

from finite dimensional data $y\in\mathbb{R}^{J}$ satisfying

[TABLE]

Here the finite number of observations are corrupted by noise $\varepsilon^{c}\eta$ of size $\varepsilon^{c}$ for which we assume that $\eta$ is a centred Gaussian $\mathcal{N}(0,\Sigma)$ with positive definite covariance $\Sigma\in\mathbb{R}^{J\times J}$. Here $\varepsilon$ and $c$ denote constants, with $\varepsilon\ll 1$ and $c>0$ (small noise) or $c=0$ (noise on the order of the observations). Observe that from an application perspective the space $BV_{\rm binary}(D)$ is a natural model for binary images. The problem of determining $u\in BV_{\rm binary}(D)$ from $y$ given by (1.1) thus constitutes a canonical binary inverse problem for which the issue is to recover the interface between the different domains in $D$ in which function $u$ takes its two values. For a given $u$ a measure of the discrepancy with the data is the following scaled misfit functional

[TABLE]

A typical deterministic approach to recover of $u$, based on this misfit, would be to specify a model prior space $\mathcal{P}$ for $u$ and to regularize the misfit functional by addition of a functional $\mathcal{R}(u)$ defined on $\mathcal{P}$, and then to seek a solution to the optimization problem

[TABLE]

An intuitive and common method of regularization for binary problems is to penalize the perimeter of the interface. In the case of $\mathcal{P}=BV_{\rm binary}(D)$ this leads to

[TABLE]

where $\int_{D}|\nabla u|$ denotes the total variation of $u$ and $\sigma>0$ is a parameter. In this way, the minimization over $BV_{\rm binary}(D)$ identifies perimeter regularization with total variation [57] and the Mumford-Shah [50] approaches, since these regularizations coincide on binary functions.

Often perimeter regularization is relaxed to a convex regularization by allowing values of $u\in[-1,1]$. An alternative, non-convex, relaxation is to use a Cahn-Hilliard functional to approximate the perimeter functional. For example one might consider, for a small parameter $\tilde{\varepsilon}>0$,

[TABLE]

where $W(\cdot)$ is a double well potential. The minimizers of this Cahn-Hilliard functional are known as phase-fields, and this relaxation is often referred to as a phase-field regularization. In appropriate circumstances this $\Gamma-$converges to (1.4) in the limit $\tilde{\varepsilon}\to 0$.

However since the unknown observational noise $\eta$ has an assumed Gaussian distribution, it is natural to take a probabilistic approach to the recovery of $u$ and model uncertainty about $u$, and hence the interface between different domains, through a probability distribution. This leads to Bayesian formulations of the problem in which a prior probability distribution is specified on the unknown function, and the likelihood of the data is used to compute a posterior probability distribution on the unknown function, given the data. The prior probability distribution imposes a prior space ${\mathcal{Q}}$ where, almost surely, samples from the posterior distribution live; the mean or mode of the posterior distribution will typically live in a smoother space ${\mathcal{P}}\subset{\mathcal{Q}}$; this space ${\mathcal{P}}$ will be analogous to the prior space ${\mathcal{P}}$ described above.

We adopt two approaches. In the first, the level set method, we reformulate the inverse problem as determining smooth functions $v$ whose zero level set defines the interface in the unknown function $u$. Specifically the sign of $v$ defines $u$ and the prior probability distribution on $v$ yields a space ${\mathcal{Q}}$ containing $C^{1}$ functions. The pushforward measure on $u$, defined by the sign function, has support in $BV_{\rm binary}(D)$. In the second approach, motivated by phase-field regularization, we relax the prior measure to allow for smooth functions $u$ with sharp interfaces near zero; the implied space ${\mathcal{Q}}$ for functions $u$ contains $H^{s}$ functions for any $s<2-d/2$, whilst ${\mathcal{P}}$ contains $H^{2}$ functions.

We define a prior distribution over the smooth function $v$ or over $u$, and formulate an associated likelihood determined by $J(u)$. Bayes’ theorem is then employed in a form which implies that the posterior probability measure is absolutely continuous with respect to the prior probability measure. Taking the sign of such distributions on $v$ yields a distribution for $u$. By sampling the distribution one can approximate the mean. It is then interesting to make a connection between this mean and the solution to the perimeter regularization problem (1.4).

The main goals of the paper are to investigate how the connection to perimeter regularization appears for different (Bayesian) formulations of the inverse problem, and to demonstrate the performance and applicability of these formulations for both the linear inverse problem (1.1), and nonlinear generalizations.

1.2 Background

There are many problems in the physical sciences where piecewise constant reconstruction is of interest, for example in subsurface inversion and imaging [12, 14, 29, 39, 48, 10] and other problems in the physical sciences [28]; the problem of image deblurring [34] (in particular, for barcodes and QR codes [16, 43, 65, 63, 42, 55]) is also of interest in the context of piecewise constant reconstruction. We draw our motivation from these problems and our numerical experiments are based on imaging problems possessing a variety of geometric interfaces, smooth and including edges.

A seminal paper linking probabilistic approaches to classical numerical methods is [27], and a review describing developments since then can be found in [18]. In the context of inverse problems the link between Bayesian and classical approaches is well-understood in the setting of Gaussian random field priors: the Bayesian maximum a posteriori (MAP) estimator [44, 21] is then the solution of a Tikhonov-Phillips regularized least squares problem [32]. When more complex priors are used the connection between classical and Bayesian perspectives is more subtle, even for linear inverse problems [6, 36, 7, 8, 46]; see [1, 17, 35] for recent work generalizing [21] beyond the Gaussian prior setting. Two interesting approaches to Bayesian inversion, both using thresholding as we do in this paper, may be found in [51] and [38]. In the one dimensional setting an interesting construction of random functions with finite TV norm may be found in [19]; a Poisson process is used to define points of discontinuity, with smooth processes between these discontinuity points.

In interface reconstruction, classical methods have been dominated by inversion techniques which penalize the length of the perimeter between different subdomains. Two contrasting approaches for describing interfaces are the use of a level set of a continuous function or a characteristic function which takes just two values. Total variation (TV) regularization has played a central role [57] and has been shown to lead to empirically effective methods which are computationally efficient. The phase-field representation of interfaces is described in [23]. The method approximates the perimeter using a scaled gradient energy and double well potential for which minimizers have diffuse interfaces of width a small length scale and which encloses a zero level set of the minimizer. In contrast, the level set approach of [60] represents interfaces as level sets of continuous fields. See [5, 26, 58] for applications of phase-field and level-set ideas in classical, non-Bayesian, inversion for interfaces. For simplicity this paper focuses primarily on recovery of a binary function, taking two known values, with unknown interface. In the more general setting of recovering unknown piecewise continuous function, in which the interface and the values of the function off the interface, are unknown, the classical TV and Mumford-Shah perimeter regularization methods become distinct. For an elliptic problem, [13] use TV regularization on a level set function to penalize both perimeter length and jump discontinuities; the Mumford-Shah minimization can be written over a suitable space to jointly minimize the function and its set of discontinuity [54].

Computer power has started to render Bayesian inversion techniques tractable in some applications [44, 64, 22], enabling uncertainty quantification. In this paper we address the question of how perimeter regularization appears within Bayesian inversion techniques for the reconstruction of binary function $u$. This is a notoriously difficult problem, as made transparent in the paper [47] which showed that use of discrete total variation regularization, in a Bayesian setting, does not lead to a meaningful problem in the continuum limit; this work led to the development of new Besov priors in [46], and the combination of TV and Gaussian priors [67]. Other approaches, with demonstrable numerical success, make use of the introduction of hyperparameters [48, 9]. Instead of approximating the TV regularization, one can derive Bayesian approaches based on the Mumford-Shah functional [36]; these methods can be extended to higher order functionals such as Blake-Zisserman [11]. Our probabilistic level set based method generalizes to a hierarchical Bayesian approach that learns the unknown continuous function off the unknown interfaces, as well as the interface itself; the interface may be viewed as a nonparametric hyperparameter.

1.3 Our Contribution

In detail our contributions are as follows:

We formulate a Bayesian level set approach and establish conditions under which this leads to posterior samples with almost surely finite perimeter, and, hence, almost surely finite total variation (TV) norm. This demonstrates that TV regularization arises naturally out of appropriately chosen Bayesian formulations of inversion. 2. 2.

We formulate a Bayesian based phase-field approach and establish a link with perimeter regularization through its MAP estimator. We prove, for appropriate choice of prior distribution and parameters carefully scaled with respect to $\varepsilon$ (and so not strictly Bayesian), that the maximum a posteriori (MAP) estimator for this phase-field approach has a $\Gamma-$limit as $\varepsilon\to 0$. This limit is exactly the perimeter (TV) regularization of the least squares fidelity objective function. 3. 3.

For a linear inverse problem we provide numerical investigations of these approaches; we also compare to (widely used) Gaussian process regression which is a natural method in this linear setting. These investigations demonstrate that the level set approach may be implemented quite cheaply in comparison with the phase-field approach, for similar levels of reconstruction accuracy. Also it is demonstrated that the level set approach can learn the true perimeter. Gaussian process regression also performs well at a low computational cost for the linear problem, but is not readily extended to nonlinear problems. 4. 4.

We provide numerical evidence for the flexibility of the Bayesian level set approach, by showing an application to a nonlinear inverse problem arising from the eikonal equation. Within this context, we also show that hyperparameters contained in the statistical model may also be efficiently learnt.

1.4 Some Notation

We use $|\cdot|$ to denote the Euclidean norm on $\mathbb{R}^{J}$. Let $C_{\#}^{k,\gamma}(\bar{D}),k\geq 0$ denote the space of real valued continuous periodic functions on $\bar{D}$ whose derivatives up to order $k^{\rm th}$ derivative are Hölder continuous with exponent $\gamma$. By virtue of continuous embedding $K$ is a bounded linear operator on $C_{\#}^{k,\gamma}(\bar{D})$ for any integer $k\geq 0.$ Also let $H^{k}_{\#}(D),k\geq 0$ denote the restriction to periodic functions of the Sobolev space of $H^{k}(D)$ of $k$-times weakly differentiable real-valued functions on $D$. These Sobolev spaces are readily characterized as weighted $\ell^{2}$ spaces on Fourier coefficients [56]. Let $X$ denote the space $C(\bar{D})$, restricted to periodic functions, and let $H$ denote $L^{2}(D).$

Fix constants $\delta>0,\tau>0,q\geq 0$ and $a_{i}\geq 0,i=1,2,3$. Denote $\tilde{\delta}=(\delta,\tau,q)$ and $\tilde{a}=(a_{1},a_{2},a_{3})$. We define a covariance operator $\mathcal{C}_{\varepsilon,\tilde{\delta},\tilde{a}}$ implicitly as the solution operator corresponding to the weak formulation of the following elliptic boundary value problem: given $f\in H$ find $\eta\in H^{2}_{\#}(D),$ so that

[TABLE]

Elliptic regularity gives $\eta\in H^{4}_{\#}(D)$ and so we may define $(\mathcal{C}_{\varepsilon,\tilde{\delta},\tilde{a}})^{-1}:H^{4}_{\#}(D)\rightarrow H$ by $(\mathcal{C}_{\varepsilon,\tilde{\delta},\tilde{a}})^{-1}\eta=f$ for $f\in H.$ The Hilbert space $\mathcal{E}_{\varepsilon,\tilde{\delta},\tilde{a}}$, is defined to be $H^{2}_{\#}(D)$ endowed with the norm

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the standard $L^{2}(D)$ inner-product. By polarization an inner-product is defined on $\mathcal{E}_{\varepsilon,\tilde{\delta},\tilde{a}}$. The three parameters $\delta,q$ and $\tau$ weight the contributions of the $H^{2}_{\#}(D),H^{1}_{\#}(D)$ and $L^{2}(D)$ terms appearing in the Hilbert space $\mathcal{E}_{\varepsilon,\tilde{\delta},\tilde{a}}$ norm whereas the parameters $a_{1},a_{2}$ and $a_{3}$ scale these terms with respect to powers of $\varepsilon$. The Hilbert space $\mathcal{E}_{\varepsilon,\tilde{\delta},\tilde{a}}$ is exactly the Cameron-Martin space for the Gaussian $\mathcal{N}(0,\mathcal{C}_{\varepsilon,\tilde{\delta},\tilde{a}})$, [22, Definition 6.26]. In the following we write $C$ and $\mathcal{E}$, with the dependence on the parameters being understood. In the computations it is made clear which values of the parameters are chosen.

Remark 1**.**

We note that including an $H^{2}_{\#}(D)$ contribution in the Cameron-Martin norm is required in dimensions $d=2,3$ in order to ensure that the underlying Gaussian is supported on continuous functions. In dimension $d=1$ it is possible to remove this requirement. [61].

Remark 2**.**

*The requirement that $\tau>0$ is made to ensure that $(\mathcal{C}_{\varepsilon,\tilde{\delta},\tilde{a}})^{-1}$ is invertible on $L^{2}_{\#}(D)$. This could also be addressed when $\tau=0$ by working on spaces of functions where the mean-value is zero. *

1.5 Outline Of The Paper

In Section 2.1 we formulate three inversion approaches for the linear inverse problem (1.1), a Bayesian level set based approach, a Bayesian based phase-field regularization and Gaussian process regression. We introduce level set and phase-field priors, both of which are non-Gaussian and state a well-posedness result for the resulting posterior distributions, as well as some relevant properties. Section 3 characterizes the MAP estimator for the phase-field prior. Under appropriate parameter scalings, we obtain perimeter regularization as a $\Gamma-$limit for the MAP estimator in the small noise regime, using the analysis in [37]. This $\Gamma-$limit links the MAP estimator to classical deterministic perimeter regularization. In Section 4 we describe testing the approaches with numerical experiments based on MCMC. We also discuss the properties of the length of level sets of Gaussian random fields, and of random fields whose law has density with respect to a Gaussian random field, and use this to demonstrate that the level set approach learns the perimeter. These numerical results establish interesting links between Bayesian level set inversion and perimeter regularization. Within Section 5 we go beyond linear inverse problems, showing that the Bayesian level set approach readily extends both to a nonlinear inverse problem arising from the eikonal equation, and to the learning of unknown hyperparameters from the prior; for hyperparameter learning we use algorithms introduced in [31] and further developed in [15]. Appendix A contains proofs of the main results relating to MAP estimators for the phase-field approach.

2 Inversion Approaches

2.1 A Bayesian Level Set Based Approach To The Inverse Problem

2.1.1 Prior And Likelihood

Let $\mathds{1}_{\cdot}$ denotes the characteristic function of a set. Our model prior for $u$ is that

[TABLE]

where $D_{+}\subset D$ is a random set defined in such a way that that ${\rm leb}({\overline{D}_{+}}\backslash D_{+})=0$, almost surely; this ensures that $u$ is almost surely in $BV_{\rm binary}$. This is achieved by working with an auxiliary variable $v$ and recovering $u$ by an application of a thresholding (sign) function $S:\mathbb{R}\mapsto\{-1,0,+1\}$ defined by

[TABLE]

The prior on $v$ is chosen to be $\mu_{0,\alpha}=\mathcal{N}(0,C^{\frac{\alpha}{2}})$, with the power satisfying $\alpha>d/2$. This implies almost sure continuity of $v$. The prior for $u$ is then defined by the pushforward of $\mu_{0,\alpha}$ under $S$

[TABLE]

This is justified since the level sets of the Gaussian random field $v$ have Lebesgue measure zero [41, Proposition 2.5]. Also by the lemma which follows, it holds almost surely that if $\alpha>1+d/2$ then $u\in\{\pm 1\}$ $D$ a.e. and has bounded total variation.

Lemma 2.1**.**

*If function $v$ is drawn from measure $\mu_{0,\alpha}$ with $\alpha>1+d/2$ then almost surely function $u$ defined by (2.1) has finite total variation norm. *

Proof 2.2**.**

*If $\alpha>1+d/2$ then almost surely $v\sim\mu_{0,\alpha}$ will be a $C^{1}$ function. The paper [45] establishes that the level set $v=0$ will then have finite length, almost surely. Since $u$ is a binary function given by (2.1) this establishes that $u$ will have finite total variation, almost surely. *

If we set $u=S(v)$ it follows that model Eq. 1.1 becomes

[TABLE]

and hence that $y|v$ is distributed as the Gaussian $\mathcal{N}(KS(v),\varepsilon^{2c}\Sigma)$. The likelihood is given by the Gaussian density proportional to

[TABLE]

$\Phi(v;y)=J\circ S(v)$ is the negative log-likelihood function.

2.1.2 Posterior

In the following proposition we to establish a relationship between the (posterior) distribution $\mu^{y}$ of the random variable $v|y$, the prior on $v$ and the likelihood on $y|v$, by means of an infinite dimensional Bayes’ Theorem, [22]:

Proposition 2.3**.**

Let $\alpha>d/2.$ Then the posterior probability $\mu^{y}$ on random variable $v|y$ is a probability measure supported on $C_{\#}^{k,\gamma}(\bar{D})$ for all $\gamma\in[0,\gamma^{\prime})$, where $k\in\{0,1,2,\cdots\}$ is chosen so that $\gamma^{\prime}:=\alpha-\frac{d}{2}-k\in(0,1],$ and is determined by

[TABLE]

*here $Z\in(0,\infty)$ is the normalization constant that makes $\mu^{y}$ a probability measure. *

Proof 2.4**.**

*This follows from an application of the theory in [41]. *

We state some beneficial properties of the the posterior, making it fit for purpose in this application. Specifically, the posterior $\mu^{y}$ has a continuous dependence on $y$, and the pushforward under $S$ defines an implied posterior $\nu^{y}$, whose samples have finite total variation.

Recall the Hellinger distance between measures $\mu$ and $\mu^{\prime}$, defined with respect to any common reference measure $\mu_{0}$ (but independent of it) and given by

[TABLE]

Proposition 2.5**.**

Given the setting of Proposition 2.3 (i) the posterior measure $\mu^{y}$ is locally Lipschitz continuous with respect to $y~{}\in~{}\mathbb{R}^{J}$; more precisely: if $|y|<\rho$ and $|y^{\prime}|<\rho$ for a constant $\rho>0$ then there is a constant $C=C(\rho)$ such that

[TABLE]

*(ii) if $\alpha>1+d/2$ then $u=S(v)$, with $v\sim\mu^{y}$, has finite total variation norm, almost surely. *

Proof 2.6**.**

*(i) follows from application of the theory in [41]. (ii) follows by noting that, since $\mu^{y}$ has density with respect to $\mu_{0,\alpha}$, anything which holds almost surely under $\mu_{0,\alpha}$, will also hold almost surely under $\mu^{y}$. Application of Lemma 2.1 gives the desired result. *

2.2 A Phase-Field Regularization Based Bayesian Approach To The Inverse Problem

The Bayesian level set approach of Section 2.1 is formulated in terms of a prior on a smooth variable $v$ whose push forward under the thresholding map gives a function $u$ taking values in $\{-1,1\}$, respecting the fact that the data takes values of the form $\{-1,1\}~{}+~{}\mbox{noise}$. Here we describe a different approach, one in which the $\varepsilon-$dependent prior on $u$ may take values anywhere in $\mathbb{R}$, but concentrates close to $\{-1,1\}$ when $\varepsilon\ll 1.$ This leads to a connection with phase-field regularization.

2.2.1 Prior And Likelihood

Fixing constants $r,b>0$ define $\Psi:X\to\mathbb{R}^{+}$ by

[TABLE]

We define the prior probability measure $\nu_{0}$ on $X$ via the Radon-Nikodym derivative

[TABLE]

where $\mu_{0}$ is the Gaussian measure $\mu_{0}=\mathcal{N}(0,C)$ on the Hilbert space $H$.The normalization $Z_{0}$ is chosen so that $\nu_{0}$ is a probability measure. Since the Gaussian measure $\mu_{0}$ is supported on continuous functions in dimensions $2$ and $3$, so is the non-Gaussian measure $\nu_{0}$. Furthermore, since $r,b>0$ and $\varepsilon\ll 1,$ this measure will concentrate on functions taking values close to $\pm 1.$ In what follows the choice of parameter $b$ will be crucial, and will be explained below; the precise value of the positive parameter $r$ is less significant.

The random variable $y|u$, given by Eq. 1.1, has a Gaussian distribution $\mathcal{N}(Ku,\varepsilon^{2c}\Sigma)$ and the likelihood is the Gaussian density proportional to

[TABLE]

$\Phi(u;y)=J(u)$ is the negative log-likelihood function.

2.2.2 Posterior

We let $\nu^{y}(\mathrm{d}u)$ denote the probability of the conditioned random variable $u|y$.The following propositions are the analogue of Proposition 2.3 and Proposition 2.5. They are proved by a straightforward application of the theory in [22, 41]:

Proposition 2.7**.**

The posterior probability $\nu^{y}$ on random variable $u|y$ is a probability measure supported on $C_{\#}^{k,\gamma}(\bar{D})$ for any $\gamma<2-d/2$ and determined by

[TABLE]

*where $Z\in(0,\infty)$ is the normalization constant that makes $\nu^{y}$ a probability measure. *

Proposition 2.8**.**

In the setting of Proposition 2.7, the posterior measure $\nu^{y}$ is locally Lipschitz continuous with respect to $y~{}\in~{}\mathbb{R}^{J}$; more precisely: if $|y|<\rho$ and $|y^{\prime}|<\rho$ for a constant $\rho>0$ then there is a constant $C=C(\rho)$ such that

[TABLE]

Remark 2.9**.**

For computations, it is convenient to draw samples from the Gaussian prior $\mu_{0}$, not the prior $\nu_{0}$; to this end we note that the posterior $\nu^{y}$ may be written as

[TABLE]

*for normalization constants $Z_{0},Z\in(0,\infty)$. *

2.3 A Gaussian Process Regression On The Inverse Problem

A popular approach for linear inverse problems is to use a Gaussian process (GP) regression to find posterior parameter distributions [52, 66]; this maybe combined with thresholding to perform classification [66], an approach we adapt here to learning a binary function.

A Gaussian process is a collection of random variables, with all finite subsets being described by a joint Gaussian distribution. Adapted to our specific inverse problem, Gaussian process regression proceeds by imposing a Gaussian prior on the unknown function, and then conditioning this on $y$ given by (1.1). We take as prior $\mu_{0}$ the Gaussian $\mathcal{N}(0,C).$ The posterior is completely described as $\nu^{y}=\mathcal{N}(m_{y},C_{y})$, with mean $m_{y}$ and covariance function $C_{y}$. This is attractive as with a closed form for $m_{y}$ and $C_{y}$, we may sample the posterior directly without the need for MCMC, thus it is very computationally efficient. Closed forms are found to be

[TABLE]

In the linear setting $m_{y}$ is the MAP estimator of the posterior, thus is the unique minimizer of the the functional

[TABLE]

This approach is unnatural from a modelling point of view, as the difference $y-Ku$ appearing in (2.4) contains comparison between data produced from a binary field and the forward map of a non-binary field. It is therefore difficult to interpret the results from this approach. Nonetheless we proceed in a fashion standard in machine learning, namely to threshold the solution of the regression to obtain a classifier; in our particular setting this corresponds to application of $S$ to $m_{y}$, or to samples from the Gaussian posterior distribution. We also note that the Gaussian process regression methodology is very specific to the linear inverse problems and does not extend directly to nonlinear forward mappings.

3 MAP Estimators, Phase Field Regularization And $\Gamma$-convergence

A MAP estimator of a Bayesian posterior distribution maximizes the posterior probability. Intuitively, the MAP estimator locates points in $X$ at which arbitrarily small balls will have maximal probability. It is defined as follows, [21].

Definition 3.1**.**

A point $\bar{z}\in X$ is a MAP estimator for the posterior measure $\nu^{y}$ if

[TABLE]

where

[TABLE]

*and $B^{\rho}(z),\rho>0$ is the ball centred at $z\in X$ with radius $\rho$. *

We explore this concept in the context of phase-field regularization. Set $\Phi:X\times\mathbb{R}^{J}\to\mathbb{R}^{+}$ to be the sum

[TABLE]

where $\Psi$ is defined in (2.2). We define the Onsager-Machlup functional $J^{\varepsilon}$, associated with the measure $\nu^{y}$ by

[TABLE]

Recall $(\mathcal{E},\|\cdot\|_{\mathcal{E}})$ is the Hilbert space and corresponding norm defined in Section 1.4. From a probabilistic perspective, $(\mathcal{E},\|\cdot\|_{\mathcal{E}})$ is the Cameron-Martin space associated with the Gaussian measure $\mathcal{N}(0,C)$ [22, Definition 6.26]. We have the following result demonstrating the role of the Onsager-Machlup functional defined on the Cameron-Martin space from [21, Theorem 3.5].

Proposition 3.2**.**

*There exists a MAP estimator for the posterior measure $\nu^{y}$ which is a minimizer of the functional $J^{\varepsilon}$. *

The functional $J^{\varepsilon}(u)$ can be written as

[TABLE]

where

[TABLE]

We consider the case where

[TABLE]

With these parameter constraints the functional $I^{\varepsilon}(u)$ becomes, for $u\in H^{2}_{\#}(D)$,

[TABLE]

and $I^{\varepsilon}(u)=+\infty$ when $u\in H\setminus H^{2}_{\#}(D)$.

Definition 3.3**.**

Define the following two functionals and constant

[TABLE]

[TABLE]

[TABLE]

Based on the work of Hilhorst et al [37]. we have the following theorem for $\Gamma-$ convergence of the functional $I^{\varepsilon}$:

Theorem 3.4**.**

Then

[TABLE]

*in the sense of $\Gamma-$convergence in the strong $L^{1}(D)$ topology. *

Proof 3.5**.**

*See Appendix A. *

This theorem shows that the MAP estimator is, for small observational noise $\varepsilon^{c}\eta$, close to a perimeter regularization. Furthermore, since $2a_{1}+3>0$, Eq. 3.2 together with the preceding $\Gamma-$limit theorem suggest that, when $\varepsilon\ll 1$, the measure will approximately concentrate on a single point close to a minimizer of $I_{0}^{\delta}.$ Our numerical results, presented in the next section, support this conjecture.

Remark 3.6**.**

*This establishes a link with perimeter regularization for the Bayesian phase-field approach at the level of the MAP estimator. This is not available for the Bayesian level set method because MAP estimators do not exist [41]. Conversely, the Bayesian level set approach has a link to perimeter length at the level of samples from the posterior (see Lemma 2.1), which is not true for the phase-field approach. We explain why this is the case. Recall that almost sure properties of the prior are inherited in the posterior. Prior samples are drawn from the centred Gaussian with covariance $C$, which corresponds to choosing $\alpha=2$. In dimension $d\geq 2$ we thus do not have $\alpha>1+d/2$ and we cannot deduce that samples have finite perimeter almost surely. (Numerical results, reported in subsection 4.5, demonstrate that the lower bound $\alpha=1+d/2$ is sharp and that random draws beneath this value do not have finite perimeter almost surely). MAP estimators, on the other hand, will live in the Cameron-Martin space of the underlying Gaussian reference measure and are necessarily smoother than draws from the measure itself [21]. Thus there is no contradiction between the fact that the MAP estimator, for small $\varepsilon$, approximately penalizes the perimeter, whilst samples from the posterior may have infinite perimeter. *

Remark 3.7**.**

*The phase-field approach has the (undesirable) property that, for $c>0$, the prior construction depends on the noise-level $\varepsilon^{c}$ through (3.3) meaning that it is not strictly Bayesian. *

4 Numerical Simulations

4.1 Test Problems

We test the three inversion techniques on three images referred to as Truth A, Truth B and Truth C. These are three fields $u^{\dagger}$ lying in the set of $BV_{\rm binary}$ images and are illustrated in the obvious way in Fig. 1. Truth A and Truth B are observed on a uniform grid of $15\times 15$ points, Truth C is observed at 50 uniformly distributed points, and all of these observations are corrupted by additive Gaussian noise as in equation Eq. 1.1. Pointwise observation does not fit our theory as we assume $K$ is linear on $L^{1}(D)$; however mollification can be used to address this and leads to results which are not different in any substantive way, as noted in [40].

In order to avoid an inverse crime [44], Truth A and Truth B are generated on a square mesh of $2^{16}$ points, and Truth C is generated on a square mesh of $320^{2}$ points, but the Gaussian random fields are constructed over a mesh of $2^{14}$ points ($N=2^{7}$ in Section 4.2). We perform numerical experiments in both the small noise and order one noise regimes.

4.1.1 Small Observational Noise Set-Up

We set $c=3/2$ and $\varepsilon=0.01.$ The implied standard deviation of the observational noise is thus $0.001.$ We make the choices of parameters in the prior covariance operator $C(=C_{\varepsilon,\tilde{\delta},\tilde{a}})$, $a_{1}=0$, $a_{2}=1$, $a_{3}=0$, $b=4$ and $r=1$ for both the Bayesian level set approach and the phase-field approach. For the Bayesian level set approach we set $\delta=1$, $q=0$, $\tau=50$ and $\alpha=3$ whereas for the phase-field approach $\delta=0.01$, $q=0.1$, $\tau=1$ and $\alpha=2$. Thus we ensure that the relations Eq. 3.3 hold so that the phase-field MAP estimator for $\nu^{y}$ approximates the minimizer of $I_{0}^{\delta}$ as given in Theorem 3.4 and we expect the posterior mass to concentrate fairly close to this MAP estimator. Note that in general we need not insist on the parameters being related via Eq. 3.3 for the level set formulation; this is because, unlike the phase-field formulation, there is no MAP estimator whose properties we are seeking to control via parameter choices. For these small noise experiments the GP regression used the same parameters as for the Bayesian level set method.

4.1.2 Order One Observational Noise Set-Up

We set $c=0.$ Note that now $\varepsilon$ does not enter the observational noise; it is simply a parameter that enters the prior and so the phase-field formulation is truly Bayesian. With this choice of $c$ we require, for the phase-field formulation, $a_{1}=-3/2,a_{2}=-1/2,a_{3}=-1,b=1$. We also set $\delta=100$, $q=0.1,\tau=1$, $\alpha=2$ and $r=1$. For the level set formulation we retain the same choice of parameters as for the small noise case above. For the GP regression we use the same parameters as for the phase-field approach for these order one noise experiments.

4.2 Sampling From

The Gaussian Prior Space

We describe how to sample numerically from Gaussian priors $\zeta_{0}=\mathcal{N}(0,\mathsf{C})$ that are key to the techniques outlined in the preceding two subsections. Here $\mathsf{C}$ is either $C^{\frac{\alpha}{2}}$ or $C$. We consider the case that $D$ is the unit square $(0,1)^{2}$. Let $\{\lambda_{k}\}$ denote the eigenvalues of $\mathsf{C}$ in increasing order with corresponding $L^{2}(D)$-normalized eigenfunctions (which are Fourier modes) $\{\varphi_{k}\}$. Then samples $z$ from $\zeta_{0}$ may be expressed through the Karhunen-Loève expansion as

[TABLE]

We implement an approximation to this by jointly approximating the field via spectral truncation and evaluation on a discrete grid of points. Such an approximation may be efficiently implemented using the Fast Fourier Transform. We work on a uniformly spaced grid $\{x_{i}\}$ of $N^{d}$ points in $D$. An approximate sample on this grid is then given by

[TABLE]

All of our numerical results are performed on the two dimensional grid which arises from this approach to generating Gaussian random fields in dimension $d=2.$ In practice we choose $N=2^{7}$, and so the discrete grid for our inversion is $2^{14}$ points.

4.3 MCMC Simulation

Markov Chain Monte Carlo (MCMC) simulations may be used to sample measures $\nu^{y}$.

In all MCMC runs we generate $10^{6}$ samples and, when computing means, discard the first $5\times 10^{5}$ samples as burn-in .

We employ the preconditioned Crank-Nicolson (pCN) algorithm [2, 20] which may be used to sample any measure $\sigma$ of the form

[TABLE]

without computing derivatives of $A(\cdot).$ Both of our posterior measures can be written in this form. For the Bayesian level set approach we have

[TABLE]

whereas for the phase-field approach we have

[TABLE]

where $\mu_{0},\mu_{0,\alpha}$ are Gaussian measures, and with the appropriate normalization constants $Z$. The pCN method has the advantage that, unlike the standard Random Walk Metropolis MCMC algorithm, its rate of convergence to equilibrium can be bounded below independently of the number of terms used in the truncated Karhunen-Loève expansion described in Section 4.2 [33].

In the notation of [20] for the pCN algorithm, $\beta\in(0,1]$ denotes the proposal variance parameter. Note that larger $\beta$ tends to lead to smaller acceptance probability, but to greater exploration of state space when steps are accepted; the optimal $\beta$ is a trade-off between these two competing effects. Depending on the noise model and number of observations, we take the proposal standard deviation parameter $\beta$ between $0.02$ and $0.1$ for level set simulations and $\beta$ between $0.002$ and $0.02$ for phase-field simulations. These choices are made in order to balance acceptance rate and size of proposed move with a view to optimizing the convergence rate of the Markov chain.

We simply assume that the resulting Markov chains $\{\eta^{(m)}\}$ are ergodic and make the approximation for an associated measure $\nu^{y}$ that

[TABLE]

where the error $e_{M}$ is Gaussian with variance $c_{g}/M$, and $c_{g}$ is the integrated auto-correlation of $g(\eta^{(m)}).$ We do not impose specific stopping criteria on the Markov chains, rather we will examine the approximation qualities derived from the chains, as a function of $M$, and study the convergence to equilibrium of quantities of interest; in particular in Section 4.3.1 we compare the acceptance probability of the chain, as a function of $M$, for the level set and phase-field approaches. The samples up to step $M$ can then be used to produce point estimates for the fields, by calculating, for example, their mean or the sign of their mean. We compare the cost of sampling versus the quality of reconstruction with these point estimates, for differing formulations.

Remark 4.1**.**

*The theory in [33] demonstrates ergodicity for problems similar to those arising in the phase-field formulation. Developing an analogous theory for the level set formulation is an open and interesting research direction; however our numerics do suggest that ergodicity holds in this case too. *

Remark 4.2**.**

*Preliminary numerical results for the one dimensional analogue of the problem studied here may be found in [61]. They are consistent with what we report here in dimension two. *

4.3.1 Computational Cost

For MCMC sampling, which we use for both the phase-field and level set approaches, every set of the Markov chain requires an evaluation of $A(u).$ Due to the presence of an extra integral term, this evaluation will typically be more expensive for the phase-field model than the level set model; for the simulations performed here, evaluation of $A(u)$ is approximately twice as expensive for the phase-field model than for the level set model. For the GP regression simulations no sampling is required and so the computational cost is significantly cheaper; the means were calculated from the expression in Section 2.3, with the cost arising from the matrix multiplications and inversion involved.

For the phase-field and level set approaches, the most significant discrepancy in computational cost arises from the statistical properties of the Markov chain used to sample the posterior approximately. In Fig. 2 we show the evolution of the local acceptance rates of proposed states for Truth A with small observational noise, as a function of $M$. The evolutions are similar for the other datasets and so are not presented for brevity.The parameter $\beta$, the proposal variance, is chosen so that the acceptance probability is neither close to one nor zero; recall that this results in an order of magnitude smaller value of $\beta$ for the phase-field method in comparison with level set, meaning that the former method makes a much slower exploration of the posterior distribution. Fig. 2 suggests that the phase-field chains have not reached equilibrium until after at least $5\times 10^{5}$ samples, whereas the level set chains converge much earlier. This is illustrated in Fig. 3, which shows a selection of samples for $M=\mathcal{O}(10^{4})$ for Truth B with small observational noise. With $M=10^{4}$ samples, the three inclusions have already been identified by the level set chain, however after $M=3\times 10^{4}$ samples the phase-field chain has only started to identify a second inclusion. Thus, even though for both models we produced the same number of samples, it would have sufficed to terminate the level set chains much earlier, significantly reducing the computational cost.

The fact that the acceptance rates for the phase-field chains are lower than those for the level set chains, despite the proposal standard deviation parameter $\beta$ being one tenth of the size can be understood as follows. Note that the measure $\nu^{y}$ can informally be thought of as having Lebesgue density proportional to $\exp(-J^{\varepsilon}(u))=\exp(-\varepsilon^{-2a_{1}-3}I^{\varepsilon}(u))$. Thus for small $\varepsilon$ the probability mass is concentrated in a small neighborhood of critical points of $I^{\varepsilon}\approx I^{0}_{\delta}$. The MCMC simulations for $\nu^{y}$ could hence be viewed as a form of derivative-free optimization for the functional $J^{\varepsilon}$.

4.4 Reconstruction Of The Means

4.4.1 Small Observational Noise

In Fig. 4 we present sample means associated with small-noise observations for the phase-field, level set and GP regression models, both with and without thresholding by $S$. Note that the phase-field and GP regression models attempt to fit the un-thresholded field to the data points, whereas the level set method attempts to fit the thresholded field; the un-thresholded field for the level set method is hence on a different scale to the other two models.

For Truth A and Truth B, the general quality of the reconstruction is similar for all three models after thresholding, though the level set method does not overfit to the datapoints as significantly as the other two methods; this overfitting for the phase-field and GP regression is manifest in a boundary for the largest inclusion in Truth B which has variations on the scale of the observational noise. Another noticeable effect in the quality of the phase-field and GP regression, manifest in Truth A, is that the edges of the circular inclusion are rendered approximately piecewise linear; this might be ameliorated by using a small mesh increment to $\varepsilon$ ratio. The level set method has no small length scale to resolve, and hence does not suffer from this effect.

For all three models reconstruction of Truth C is fairly inaccurate as the sparse observation network does not resolve the length scale on which the true field varies. The level set and GP regression models perform similarly, whereas the phase-field model places much more mass into the positive class; it is likely that this reflects a lack of convergence of the Markov chain for the phase-field model.

4.4.2 Order One Noise Reconstructions

In Fig. 5 the sample means associated with order one observational noise are shown. As would be expected, reconstruction quality is generally poorer than for the small-noise observations, though overfitting to the observational noise is no longer an issue for the phase-field and GP approaches. The three models perform similarly, though there seems to be an increased amount of penalization on the length of the interface from left-to-right. Without thresholding, the GP regression means provide poor estimates of the truth in terms of scale, due to the far weaker influence of the likelihood and lack of prior information enforcing values close to $\pm 1$.

4.5 Perimeter Learning

Here we study perimeter learning for the Bayesian level set method. The length of the zero level set of $\eta^{N}$ may then be approximated by using the discrete variation of $w^{N}:=S(\eta^{N})$,

[TABLE]

where the operator $\mathsf{D}^{N}$ approximates the gradient on the grid $\{x_{i},y_{j}\}$ via central differences. Using this we investigate numerically the length of the level sets of these samples. We have shown in Lemma 2.1 that choosing $\alpha>1+d/2$ is sufficient to ensure almost sure finite length of level sets. Numerical experiments indicate that this is a sharp result. In Fig. 6 interface lengths for prior distributions approximated as described above for a single realization of $\{\xi_{k}\}$ in Eq. 4.1 are displayed as a function of $N$ for varying values of $\alpha$. We use $d=2$ and observe that for $\alpha<2$ the length of the interface diverges algebraically with $N$ (left hand panel), for $\alpha=2$ it diverges logarithmically (right hand panel shows this best), and for $\alpha>2$ it converges to a constant (both left and right hand panels show this). The results, then, suggest that level sets have finite length almost surely if and only if $\alpha>1+d/2.$

In order to compare the perimeter distribution between the prior and posterior, the choice $\alpha>1+d/2$ ensures that the length of the zero level set is well-defined so in two dimensions we take $\alpha=3$. The results for recovery of Truth B are shown in Fig. 7. Whilst the perimeter still retains some variation under the posterior, the variation is much lower and, in contrast to the prior, it is concentrated around the true value of the perimeter. We see that though the Bayesian level set approach does not explicitly penalize the perimeter, it has the ability to estimate the perimeter, and quantify uncertainty in the estimation.

5 Eikonal Equation

In this section we build on what we have learned so far for the linear inverse problem defined by (1.1) and use it to study a nonlinear inverse problem from the eikonal equation. Gaussian process regression is fast to implement, and appears to give qualitatively comparable accuracy to the Bayesian level set method; but it does not generalize to nonlinear problems and so we do not consider it further. The experiments in the previous section, set-up in Section 4.1, suggest that the Bayesian level set approach to binary recovery has two advantages over the Bayesian phase-field formulation, for the linear inverse problem considered: the level set method is faster and draws from the posterior contain information about the true perimeter. Thus we focus attention purely on the Bayesian level set method. Within this context we also demonstrate the benefits of hierarchical Bayesian inversion.

5.1 The Forward Equation

Let $x_{0}\in D$ be a wave emitting source, and define the first arrival time of the wave at $x\in D$ as $T(x)$. The wave passes through a medium which adjusts the wave speed according to a scalar function $u\colon D\to\mathbb{R}_{+}$ known as the slowness. It is shown in [49, 62] that $T(x)$ may be viewed as solution of a stationary Hamilton-Jacobi equation, namely the following eikonal equation:

[TABLE]

where $n$ is the outward pointing unit normal. The Soner boundary condition (5.3) ensures ray paths terminate at $\partial D$ [62]. The recovery of the slowness function $u$ from observations of arrival times is known as first arrival traveltime tomography. Extensive discussion of the well-posedness of the forward problem can be found in [24, 25, 30]. For our application, we consider a binary slowness function $u\colon D\to\{u_{\min},u_{\max}\}$ where $0<u_{\min}\leq u_{\max}$.

We define the solution map $\mathsf{G}$ mapping the slowness $u$ to the travel times $T$ via solution of the eikonal equation with source $x_{0}$. Because we are interested in binary slowness functions we also introduce

[TABLE]

here $S(v)$, is the sign function defined in Section 2.1.

To solve the forward problem we first discretize using an upwind finite difference scheme. We then use a fast marching procedure (see [59]) to solve the discrete eikonal equation. A formulation of the discretization and marching algorithm, along with a proof of numerical convergence is found in [25].

5.2 The Bayesian Inverse Problem

Let $\eta\sim\mathcal{N}(0,\Sigma)$ be a normal random variable in $\mathbb{R}^{J}$. Fix $x_{0}\in\bar{D}$. Defining the observation map $K$ taking traveltimes to $\mathbb{R}^{J}$ we define the inverse problem of finding $v$, given observed data $y$ satisfying

[TABLE]

We will assume that $K$ is defined so that the data $y$ is a set of observed first hitting times at fixed known receiver locations $\{z_{j}\}_{j=1}^{J}\in\bar{D}$. The random variable $y|v\sim\mathcal{N}(K\circ\mathsf{G}\circ\mathsf{S}(v),\varepsilon^{2c}\Sigma)$, leading to negative log likelihood defined, up to an additive constant, by

[TABLE]

We will treat problems of multiple sources $\{x_{0}^{m}\}_{m=0}^{M}$ as multiple experiments, with solution maps $\mathsf{G}^{m}$ each producing data $y^{m}\in\mathbb{R}^{J}$. The natural extension is to consider the following negative log likelihood,

[TABLE]

for the $m$ observations $y=(y^{1},\dots,y^{m}).$ We assume that we are given prior measure \zeta_{0}=\mathcal{N}\bigl{(}0,\mathsf{C}\bigr{)} and let $\zeta^{y}(\mathrm{d}v)$ denote the probability distribution of the conditioned random variable $v|y$. Then using Bayes’ theorem (see Proposition 2.3), we deduce that $\zeta^{y}$ is a probability measure supported on continuous functions, determined by

[TABLE]

with normalization constant $Z$. This problem has been formulated for piecewise constant slowness function with Ginzburg-Landau type regularization in the deterministic setting [30], where simulations and proofs of convergence of numerical schemes can be found. We choose here to use a level set formulation for the reasons discussed at the start of the section.

5.2.1 Hierarchical Inference For Inverse Lengthscale

In hierarchicalBayesian inference we employ a Gaussian prior measure $\zeta_{0}=\mathcal{N}(0,\mathsf{C}(\tau))$ in which the covariance $\mathsf{C}$ depends on parameter $\tau>0$ which we interpret as an additional unknown to be learned during the inversion process. In particular we will work with settings in which $\tau$ has interpretation as an inverse length-scale. To this end, define the hierarchical prior $\zeta_{0}$ by decomposing as follows:

[TABLE]

We call $\pi_{0}(\tau)$ the hyperprior. Generalizing the derivation of the posterior in the preceding subsection, we now find that the distribution of $v,\tau|y$ is determined by probability measure $\zeta^{y}(\mathrm{d}v,\mathrm{d}\tau)$ defined by

[TABLE]

for normalization constant $Z$. For reasons discussed in [53, 68] it can be advantageous to reparameterise the hierarchical inverse problem in terms of $(\xi,\tau)$, rather than $(v,\tau)$, where $\xi$ is a Gaussian white noise distributed as $\mathcal{N}(0,I),$ so that $v=\sqrt{\mathsf{C}(\tau)}\xi;$ the underlying latent Gaussian white noise $\xi$ may be identified with the collection of i.i.d. unit Gaussians $\{\xi_{k}\}$ used to construct prior samples in Section 4.2. Abusing notation we may write the posterior distribution $\zeta^{y}$, now for the variable $(\xi,\tau)$, as

[TABLE]

Working with variables $(v,\tau)$ as in (5.7) is refered to as the centred problem formulation; using variables $(\xi,\tau)$ as in (5.8) is refered to as the non-centred problem formulation.

Numerical evidence described in [31, 15] demonstrates that for level-set based inverse problems use of the non-centered formulation in (5.8) confers considerable advantages in terms of speed of convergence of MCMC. We thus employ the non-centered formulation. To sample from this distribution we use [15, Algorithm 6.1], known as the non-centred pCN-within-Gibbs. This algorithm updates $\xi$ and $\tau$ in separate substeps of a pCN sampling method, and we perform this in a random order each iteration.

5.2.2 Hierarchical Inference For Contrast

We also consider a model hyperparameter that originates in the application (whereas $\tau$ appears when regularizing through the prior). We investigate the contrast between the binary materials, and rewriting the relationship between $u$ and $v$, we obtain

[TABLE]

where $\kappa=(u_{\max}-u_{\min})>0$ is the contrast. Knowing the parameter $u_{\min}$, we consider the parameter $\kappa$ as a positive random variable to be learnt from the data. This increases flexibility of techniques in application as we can apply our methods to scenarios where slowness contrast is uncertain. We may include this in a non-centred pCN-within-Gibbs algorithm as above, by exploiting the form

[TABLE]

substituting this into the likelihood (5.6). We update each of $\xi,\tau,\kappa$ separately and in random order in each iteration.

5.2.3 Numerical Results: Lengthscale Hyperparameter Only

For this first test we wish to demonstrate recovery of hyperparameters with the nonlinear eikonal forward model as described in Section 5.2.1. Our domain is given by $D=[0,1]\times[0,1]$. We choose Truth D as seen in Fig. 8; here the truth has been produced by applying the slowness function $\mathsf{S}(v^{*})$ for an instance $v^{*}\sim\mu_{0}(\tau^{*})$ with hyperparameter $\tau^{*}=\exp(6.5)$. We take $u_{\min}=1.0$, and a known contrast $\kappa=0.2$ (thus $u_{\max}=1.2$).

The discretization uses an equidistributed mesh with grid spacing $h=128^{-1}$. To avoid an inverse crime we produce the data on a numerical mesh with spacing $h/2$. We take 4 sources $\{x^{m}_{0}\}_{m=1}^{4}$ and 100 receivers uniformly distributed in the (coarser) discrete domain.

We work with the parameters same as taken in Section 4.1.1 with the exception of setting $c=1$ and, of course, viewing $\tau$ as an unknown. For consistency across different source–receiver combinations we additionally scale the noise by the range of the traveltime observations, so the effective observational noise is $10^{-2}$. We take $\alpha=3$ to ensure finite perimeter.

For the hyperparameter $\tau$, we choose a lognormal prior to ensure positivity, we take prior $\pi(\log\tau)=\mathcal{N}(5,2.5)$ and initialize the Markov chain at $\log(3)$. We use the non-centred pCN-within-Gibbs algorithm [15, Algorithm 6.1] described in Section 5.2.1 with a random walk Metropolis proposal for the hyperparameter step. The step sizes were chosen to achieve $20\%-25\%$ acceptance rates for $\xi$ and $\tau$.

The result of the recovery of truth D after $2\times 10^{5}$ iterations run is displayed in Fig. 9. The recovery of the field is fairly faithful; we note that information can only be learnt between source and receiver pairs, see Truth D in Fig. 8 for their random distribution. We additionally provide the distribution of the perimeter in Fig. 10 with true perimeter marked, and one can see that this falls well within the posterior distribution with high probability mass given to a small neighborhood of the truth. In Fig. 11 we see the results of learning the hyperparameter distribution. We see the marked truth and the prior and posterior distribution, and again we see the posterior peaks near the true value and gives large mass to a small neighborhood of the truth. We note that the parameter $\log\tau$ is sampling close to its posterior within $10^{4}$ iterations in this example.

5.3 Numerical Results: Lengthscale And Contrast Hyperparameters

For this test we demonstrate recovery of the contrast in the medium. We assume that we are in the situation of performing the inverse eikonal problem where we do not have exact information on the contrast between the binary phases. We choose Truth E as seen in Fig. 8, comprising four circles; three of diameter $0.1$ and one of diameter $0.15$. In the figure we see the choice of four sources, and take $15^{2}$ equally spaced receivers over the domain. We take $u_{\min}=1.0$, and again we work with the parameters as detailed in Section 4.1.1 with the exception of setting $c=1$, and viewing $\tau$ as unknown and, as in the previous subsection, we set $\alpha=3$. We assume the contrast is lognormal with prior $\pi(\log\kappa)=\mathcal{N}(\log(0.2),0.3)$, as we require a positive prior, we also treat $\tau$ as a random unknown quantity, and take a lognormal $\pi(\log\tau)=\mathcal{N}(5,2.5)$. We initialize the MCMC method to sample these variables at $\log\tau=8$ and $\log\kappa=\log(0.1)$.

Our recovery of Truth E after $2\times 10^{5}$ iterations is recorded in Fig. 12. We see an excellent recovery of the simple geometry, in particular the colorscale shows the approximation of the mean recovered contrast to the true contrast, where $u_{\min}=1,\ u_{\max}=1.2$. This recovery is demonstrated in Fig. 13 where we see the prior and posterior densities for the interfacial length. The posterior places most weight in a small neighbourhood of the truth, and peaks nearby.

The hyperparameter recovery is displayed in Fig. 14 and Fig. 15. We see the profiles for prior and posterior for $\log\tau$ in Fig. 14, demonstrating that the posterior concentrates within a sensible range; note that there is no true $\tau$ for this example. In Fig. 15 we see the contrast recovery, and observe that the posterior places large weight close to the true value. It again takes less than $10^{4}$ iterations for both hyperparameters to draw approximately from their posterior distributions.

6 Conclusions

The paper investigates the reconciliation of perimeter and Bayesian regularization for the reconstruction of functions with interfaces, from direct or indirect noisy measurements. Three approaches are studied: Formulation 1 is based on Bayesian level method; Formulation 2 is based on Bayesian phase-field regularization; and Formulation 3 is based on Gaussian process regression and classification.

By studying a class of linear inverse problems we show that Formulation 2 exhibits perimeter regularization in the context of its MAP estimator, but not at the level of samples from the posterior distribution. Formulation 1 exhibits perimeter regularization at the level of individual samples from the posterior; there is no MAP estimator in this context. Both Formulations 1 and 2 require careful choices of constants in construction of the prior, but Formulation 2 is far more constrained in this regard. Furthermore, as a consequence of these constraints, Formulation 2 exhibits a measure concentration phenomenon meaning that MCMC based algorithms using Formulation 1 are considerably faster than those based on Formulation 2. Formulation 3 is competitive with Formulation 1 in terms of both sample properties and speed, but does not generalize beyond linear problems. We study Formulation 1 for a nonlinear inverse problem, demonstrating that it is effective in this context and, in addition, showing how hierarchical methods may be used to learn model hyper-parameters appearing in the prior.

The ideas in this paper can be combined in different ways: the methodology may be extended beyond binary-valued functions to a variety of piecewise continuous problems; or other limiting functionals could be contemplated, such as Mumford-Shah [36]; and other smoothed thresholding functions could be contemplated within the level set method, such as the double-obstacle approximation to the signum function [3, 4]. The success of the Bayesian level set method suggests that further analysis of it, as well as its deployment in new application domains, would be very valuable.

Appendix A Proofs Of Main Results

Proof of Proposition 3.2. Throughout this proof $C$ is a universal constant whose value may change between occurrences. To apply Theorem 4.12 from [21], we need to show that the function $\Phi(\cdot,y)$ is bounded from below, is locally bounded from above and is locally Lipschitz. We note that $\Phi(\cdot,y)$ is always non-negative so is bounded from below. If $\|u\|_{X}=\max_{x\in\bar{D}}|u(x)|\leq\rho$ then we may bound $|\Phi(u,y)|$ by a constant depending on $\rho$ i.e. $\Phi(\cdot,y)$ is locally bounded. For the local Lischiptz continuity, we have

[TABLE]

Assume that $\|u\|_{X}\leq\rho$ and $\|v\|_{X}\leq\rho$. Then, since $K$ is a bounded linear operator on $L^{1}(D),$

[TABLE]

The desired result follows. $\Box$

Proof of Theorem 3.4. We adapt the proof of Hilhorst et al. to allow for periodic boundary conditions and the additional $L^{2}$ norm appearing in the functional to be infimized. From Hilhorst et al., we have that if $u^{\varepsilon}\to u$ in $L^{1}(D)$ then

[TABLE]

Now we show that for each $u\in L^{1}(D)$, there is a sequence $\{u^{\varepsilon}\}\subset H^{2}_{\#}(D)$ which converges strongly to $u$ in $L^{1}(D)$ such that $\limsup_{\varepsilon\to 0}I^{\varepsilon}(u^{\varepsilon})\leq I_{0}^{\delta}(u)$. We first review the main points in the proof of Hilhorst et al. for functions $u\in H^{2}(D)$. Considering the case $I^{\delta}(u)<\infty$, without loss of generality, we assume that

[TABLE]

where $Q$ is a bounded domain, with $\partial Q\in C^{\infty}$ and $Q\subset\subset D$. The sign distance function is defined as

[TABLE]

Let $N_{h}$ be an $h$ neighbourhood of $\partial Q$ (we choose $h$ so that $h$ is less than the distance between $\partial Q$ and $\partial D$.) We choose a function $\eta\in C^{2}(\bar{D})$ such that $\eta(x)=d(x)$ for $x\in N_{h}$, $\eta(x)\geq h$ when $x\in Q\setminus N_{h}$ and $\eta(x)\leq-h$ when $x\in D\setminus(Q\bigcup N_{h})$. Let $U$ be an odd minimizer of the functional $e^{\delta}(U)$ with $\lim_{t\to\infty}U(t)=1$ and $\lim_{t\to-\infty}U(t)=-1$. We let

[TABLE]

We note that $u^{\varepsilon}(x)$ is uniformly bounded pointwise and $u^{\varepsilon}(x)\to u(x)$ for all $x\in D$. From the Lebesgue dominated convergence theorem, $u^{\varepsilon}\to u$ in $L^{1}(D)$ and in $L^{2}(D)$. Thus

[TABLE]

and, since $a>0$,

[TABLE]

To show that $\lim_{\varepsilon\to 0}I^{\varepsilon}(u)=I^{\delta}_{0}(u)$, we follow the approach of Hilhorst et al.. The integral

[TABLE]

is written as

[TABLE]

Using the exponential decay of $U,U^{\prime}$ and $U^{\prime\prime}$ at $\infty$ and $-\infty$, we deduce that the integral over $D\setminus N_{h}$ goes to [math] when $\varepsilon\to 0$ (note that $|\eta(x)|/\varepsilon>h/\varepsilon$ which goes to $\infty$ when $\varepsilon\to 0$ for $x\in D\setminus N_{h}$). The integral over $N_{h}$ is shown to converge to $I^{\delta}(u)$ when $\varepsilon\to 0$.

To adapt this proof of Hilhorst et al. to functions with periodic boundary condition on $D$, we only need to choose the function $\eta$ so that $\eta$ is periodic and $\eta(x)\geq h$ when $x\in Q\setminus N_{h}$ and $\eta(x)\leq-h$ when $x\in D\setminus(Q\bigcup N_{h})$. Such a function can be constructed as follows. Let $\psi(x)\in C^{\infty}_{0}(D)$ be such that $\psi(x)=1$ when $x$ is in a neighbourhood of $Q\bigcup N_{h}$, and $0\leq\psi(x)\leq 1$ for all $x\in D$. Let $\eta_{1}(x)$ be a smooth periodic function with $\eta_{1}(x)\leq-h$ for all $x\in D$. Using the function $\eta$ of Hilhorst et al., we define a new function

[TABLE]

The function $\bar{\eta}(x)$ satisfies the requirement. $\Box$