Joint estimation of Robin coefficient and domain boundary for the Poisson problem

In the current paper, we denote by $\tilde{{\Omega}}$ the (known) uncorroded reference geometry, which is assumed to be a slab, i.e., $\tilde{{\Omega}}=(0,L)\times (0,H)$, with L ≫ H > 0. The Poisson equation, as outlined in (2.1), is known to be invariant under certain deformations of the domain Ω [36, 45, 46]. More specifically, we let ψ denote the diffeomorphism mapping Ω to $\tilde{{\Omega}}$ given by,

$$\begin{equation}\psi :({x}_{1},{x}_{2}){\mapsto}\left({x}_{1},{x}_{2}{(f({x}_{1}))}^{-1}\right),\quad ({x}_{1},{x}_{2})\in {\Omega}.\end{equation} \tag{ 2.3 }$$

We use x to denote points in Ω, while points in $\tilde{{\Omega}}$ will be denoted by $\tilde{\boldsymbol{x}}$, i.e., x ∈ Ω and $\tilde{\boldsymbol{x}}\in \tilde{{\Omega}}$. Then, taking into account that such a diffeomorphisms leaves Γ_A unchanged ⁴ , i.e., ${\left.\psi \right\vert }_{{{\Gamma}}_{\text{A}}}=\mathrm{I}\mathrm{d}$ (the identity mapping), the functions ${\tilde{u}}_{k}(\tilde{\boldsymbol{x}})={u}_{k}{\circ}{\psi }^{-1}(\tilde{\boldsymbol{x}})$ are known to satisfy a transformed Poisson equation [36, 45, 46]:

$$\begin{equation}\begin{aligned}\hfill -\nabla \cdot (\tilde{\sigma }\nabla {\tilde{u}}_{k})& =0\hfill & \hfill & \quad \text{in}\enspace \tilde{{\Omega}}\hfill \\ \hfill \tilde{\sigma }\nabla {\tilde{u}}_{k}\cdot \tilde{\boldsymbol{\eta }}& ={g}_{k}\hfill & \hfill & \quad \text{on}\enspace {\tilde{{\Gamma}}}_{\text{A}}\hfill \\ \hfill \tilde{\sigma }\nabla {\tilde{u}}_{k}\cdot \tilde{\boldsymbol{\eta }}+\mathrm{exp}(\tilde{\beta }){\tilde{u}}_{k}& =0\hfill & \hfill & \quad \text{on}\enspace {\tilde{{\Gamma}}}_{\text{I}}\hfill \\ \hfill {\tilde{u}}_{k}& =0\hfill & \hfill & \quad \text{on}\enspace {\tilde{{\Gamma}}}_{\text{D}},\hfill \end{aligned}\end{equation} \tag{ 2.4 }$$

for k = 1, 2, ..., ℓ, where $\tilde{\boldsymbol{\eta }}$ is the outward facing normal of $\tilde{{\Omega}}$, the push forward of σ (by ψ) is given by

$$\begin{equation}\tilde{\sigma }(\tilde{\boldsymbol{x}})=\sigma (\boldsymbol{x})\frac{\boldsymbol{J}(\boldsymbol{x})\boldsymbol{J}{(\boldsymbol{x})}^{T}}{\vert \mathrm{det}(\boldsymbol{J}(\boldsymbol{x}))\vert },\qquad \boldsymbol{x}={\psi }^{-1}(\tilde{\boldsymbol{x}}),\end{equation} \tag{ 2.5 }$$

with J the Jacobian matrix of ψ, and

$$\begin{equation}\mathrm{exp}\left(\tilde{\beta }(\tilde{s})\right)=\mathrm{exp}\left(\beta (s)\right)\sqrt{1+{\left(\frac{\mathrm{d}f}{\mathrm{d}s}\right)}^{2}{H}^{2}},\quad s={\varphi }^{-1}(\tilde{s}),\end{equation} \tag{ 2.6 }$$

where φ is the restriction of ψ to Γ_I. It is worth pointing out that even if the conductivity σ in Ω is isotropic (as is the case in the current paper), the pushed forward conductivity $\tilde{\sigma }$ in $\tilde{{\Omega}}$ will in general be anisotropic. We show a schematic of the invariance property of the Poisson problem in figure 1.

The invariance of the Poisson equation, as outlined above, means that for any domain, Ω, which is diffeomorphic to $\tilde{{\Omega}}$, the solution to (2.1) can be found by solving (2.4). If we only require the solution to (2.1) in a single domain, Ω, there is little use in solving (2.4) instead. However, when considering an estimation problem, such as that considered in the current paper, the forward problem will need to be solved repeatedly. Then, since we want to estimate the domain boundary, we will require the solution to (2.1) in a number of different domains. Computing each of these solutions would require re-meshing of the domain, or at least of the inaccessible part of the boundary Γ_I. This increases the computational overhead, which can become significant if the number of forward solutions required is large. Furthermore, as the Robin boundary coefficient is defined on Γ_I, special care would need to be taken to ensure that spatial prior information, such as smoothness of the coefficient, was kept consistent across each of the domains. This would further increase the computational requirements, and could also lead to convergence issues for the estimates, see e.g. [31, 47, 48]. To reduce these issues, in the current paper, we propose only ever solving the forward problem in the fixed reference domain, $\tilde{{\Omega}}$, by exploiting the discussed invariance property.

Key to further reducing the computational costs of solving the estimation problem is efficient calculation of generalised derivatives of the forward model with respect to both the shape of the inaccessible part of the domain boundary Γ_A and the Robin coefficient β, which we now outline. A related approach is provided in [34, 43], though the procedure used there requires, among other things, a dual procedure to cope with distributional boundary conditions.

In what follows, we will slightly abuse notation and set g = g_k and u = u_k for some fixed $k\in \left\{1,2,\dots ,\ell \right\}$. Furthermore, it will be useful to consider writing the weak form of (2.1) as a(u, v) = F(v), where

$$\begin{equation}a(u,v){:=}{\int }_{{\Omega}}\sigma \nabla u\cdot \nabla v\enspace \mathrm{d}\enspace \boldsymbol{x}+{\int }_{{{\Gamma}}_{\text{I}}}\enspace \mathrm{exp}(\beta )uv\enspace \mathrm{d}\enspace {\boldsymbol{s}}_{\text{I}},\quad \text{and}\quad F(v){:=}{\int }_{{{\Gamma}}_{\text{A}}}gv\enspace \mathrm{d}\enspace {\boldsymbol{s}}_{\text{A}},\end{equation} \tag{ 2.7 }$$

for $v\in {H}_{D}^{1}({\Omega})$, with ${H}_{D}^{1}({\Omega})=\left\{v\in {H}^{1}({\Omega}):{\left.v\right\vert }_{{{\Gamma}}_{\text{D}}}=0\right\}$.

To calculate the generalised derivative for the boundary, we first take $h\in {C}_{\text{per}}^{1}([0,L])$ and introduce the perturbed problem: a^h (u^h , v^h ) = F(v^h ), for which the weak solution, u^h ∈ H¹(Ω^h ), satisfies

$$\begin{equation}{\int }_{{{\Omega}}^{h}}\sigma \nabla {u}^{h}\cdot \nabla {v}^{h}\enspace \mathrm{d}\enspace \boldsymbol{x}+{\int }_{{{\Gamma}}_{\text{I}}^{h}}\enspace \mathrm{exp}(\beta ){u}^{h}{v}^{h}\enspace \mathrm{d}\enspace {\boldsymbol{s}}_{\text{I}}={\int }_{{{\Gamma}}_{\text{A}}}g{v}^{h}\enspace \mathrm{d}\enspace {\boldsymbol{s}}_{\text{A}}\quad \forall \enspace {v}^{h}\in {H}_{D}^{1}({{\Omega}}^{h}),\end{equation} \tag{ 2.8 }$$

where ${{\Omega}}^{h}=\left\{({x}_{1},{x}_{2}):{x}_{1}\in (0,L)\enspace \text{and}\enspace 0< {x}_{2}< f({x}_{1})+th({x}_{1})\right\}$ with $t\in {\mathbb{R}}^{+}$. We next introduce the diffeomorphism which maps Ω^h to $\tilde{{\Omega}}$,

$$\begin{equation*}{\psi }^{h}:({x}_{1},{x}_{2}+th({x}_{1})){\mapsto}\left({x}_{1},({x}_{2}+th({x}_{1})){(f({x}_{1})+th({x}_{1}))}^{-1}\right).\end{equation*}$$

Finally, the Gateaux derivative, see e.g. [49, chapter 8], of u ∈ H¹(Ω) evaluated at Ω in the direction h is then

$$\begin{equation*}\mathrm{d}u({\Omega};h)=\underset{t{\searrow}0}{\mathrm{lim}}\enspace \frac{{u}^{h}{\circ}{\psi }^{h}-u{\circ}\psi }{t}.\end{equation*}$$

The potential u can be written ⁵ as u = T⁻¹ F, see e.g. [52]. Furthermore, as a in (2.7) is the symmetric form associated with T, straightforward differentiation then yields

$$\begin{align}\hfill \mathrm{d}u({\Omega};h)& =\mathrm{d}\left({T}^{-1}F\right)({\Omega};h)=\mathrm{d}({T}^{-1})({\Omega};h)F+{T}^{-1}\enspace \mathrm{d}F({\Omega};h)\hfill \\ \hfill & =-{T}^{-1}\mathrm{d}a({\Omega};h)u+{T}^{-1}\enspace \mathrm{d}F({\Omega};h).\hfill \end{align} \tag{ 2.9 }$$

Due to the nature of the admissible set of boundary perturbations, it follows that dF(Ω; h) = 0. That is to say, finding da(Ω; h) renders calculating du(Ω; h) essentially trivial.

With the two diffeomorphisms, ψ (see (2.3)) and ψ^h , in hand, we can rewrite both a and a^h using the invariance property of the Poisson problem (see section 2.1) as

$$\begin{align*}\hfill \tilde{a}(\tilde{u},\tilde{v})& ={\int }_{\tilde{{\Omega}}}\tilde{\sigma }\nabla \tilde{u}\cdot \nabla \tilde{v}\enspace \mathrm{d}\tilde{\boldsymbol{x}}+{\int }_{{\tilde{{\Gamma}}}_{\text{I}}}\enspace \mathrm{exp}(\tilde{\beta })\tilde{u}\tilde{v}\enspace \mathrm{d}{\tilde{\boldsymbol{s}}}_{\text{I}},\hfill \\ \hfill {\tilde{a}}^{h}({\tilde{u}}^{h},{v}^{h})& ={\int }_{\tilde{{\Omega}}}{\tilde{\sigma }}^{h}\nabla {\tilde{u}}^{h}\cdot \nabla {\tilde{v}}^{h}\enspace \mathrm{d}\tilde{\boldsymbol{x}}+{\int }_{{\tilde{{\Gamma}}}_{\text{I}}}\enspace \mathrm{exp}({\tilde{\beta }}^{h}){\tilde{u}}^{h}{\tilde{v}}^{h}\enspace \mathrm{d}{\tilde{\boldsymbol{s}}}_{\text{I}},\hfill \end{align*}$$

respectively, where $\tilde{\sigma }$ and $\mathrm{exp}(\tilde{\beta })$ are defined in (2.5) and (2.6) respectively, $\tilde{u}=u{\circ}\psi$, while ${\tilde{u}}^{h}={u}^{h}{\circ}{\left({\psi }^{h}\right)}^{-1}$, and

$$\begin{equation*}\begin{aligned}\hfill {\tilde{\sigma }}^{h}(\tilde{\boldsymbol{x}})& =\sigma ({\boldsymbol{x}}^{h})\frac{{\boldsymbol{J}}^{h}({\boldsymbol{x}}^{h}){\boldsymbol{J}}^{h}{({\boldsymbol{x}}^{h})}^{T}}{\vert \mathrm{det}({\boldsymbol{J}}^{h}({\boldsymbol{x}}^{h}))\vert },\qquad \mathrm{exp}\left({\tilde{\beta }}^{h}(\tilde{s})\right)=\mathrm{exp}\left(\beta (s)\right)\sqrt{1+{\left(\frac{\mathrm{d}\left(f(s)+th(s)\right)}{\mathrm{d}s}\right)}^{2}{H}^{2}},\hfill \end{aligned}\end{equation*}$$

with J ^h the Jacobian matrix of ψ^h , and where ${\boldsymbol{x}}^{h}={\left({\psi }^{h}\right)}^{-1}(\tilde{\boldsymbol{x}})$ and $s={\varphi }^{-1}(\tilde{s})$. We can then write the derivative as,

$$\begin{align}\hfill \mathrm{d}\tilde{a}({\Omega};h)\tilde{u}& =\underset{t\to 0}{\mathrm{lim}}\left({\int }_{\tilde{{\Omega}}}{\tilde{\sigma }}^{h}\nabla \tilde{u}\cdot \nabla \tilde{v}\enspace \mathrm{d}\tilde{\boldsymbol{x}}+{\int }_{{\tilde{{\Gamma}}}_{\text{I}}}\enspace \mathrm{exp}({\tilde{\beta }}^{h})\tilde{u}\tilde{v}\enspace \mathrm{d}{\tilde{\boldsymbol{s}}}_{\text{I}}\right.\hfill \\ \hfill & \quad -\left.{\int }_{\tilde{{\Omega}}}\tilde{\sigma }\nabla \tilde{u}\cdot \nabla \tilde{v}\enspace \mathrm{d}\tilde{\boldsymbol{x}}-{\int }_{{\tilde{{\Gamma}}}_{\text{I}}}\enspace \mathrm{exp}(\tilde{\beta })\tilde{u}\tilde{v}\enspace \mathrm{d}{\tilde{\boldsymbol{s}}}_{\text{I}}\right){t}^{-1}\hfill \\ \hfill & =\underset{t\to 0}{\mathrm{lim}}\left({\int }_{\tilde{{\Omega}}}\left({\tilde{\sigma }}^{h}-\tilde{\sigma }\right)\nabla \tilde{u}\cdot \nabla \tilde{v}\enspace \mathrm{d}\tilde{\boldsymbol{x}}\right){t}^{-1}+\underset{t\to 0}{\mathrm{lim}}\left({\int }_{{\tilde{{\Gamma}}}_{\text{I}}}\left(\mathrm{exp}({\tilde{\beta }}^{h})-\mathrm{exp}(\tilde{\beta })\right)\tilde{u}\tilde{v}\enspace \mathrm{d}{\tilde{\boldsymbol{s}}}_{\text{I}}\right){t}^{-1}\hfill \\ \hfill & =\mathrm{d}\tilde{a}(\tilde{\sigma };\mathrm{d}\tilde{\sigma }({\Omega};h))+\mathrm{d}\tilde{a}(B;\mathrm{d}B({\Omega};h)),\hfill \end{align} \tag{ 2.10 }$$

where we have used the substitution $B=\mathrm{exp}(\tilde{\beta })$, and where $\mathrm{d}\tilde{\sigma }({\Omega};h)\in {S}^{2}\left(C\left([0,L]\right)\right)$, the space of symmetric matrices with elements in $C\left([0,L]\right)$ and $\mathrm{d}B({\Omega};h)\in C\left([0,L]\right)$.

In the current set up, due to the assumptions on the domain shape perturbations, each of the derivatives in (2.10) are straightforward to compute as follows. First, since the push forward of σ is given by

we have $\mathrm{d}\tilde{\sigma }({\Omega};h)=\sigma (\boldsymbol{x})\mathbf{\Sigma }(\boldsymbol{x})h({x}_{1})$, with

On the other hand, inline with e.g. [53], by setting $\hat{h}=\mathrm{d}\tilde{\sigma }({\Omega};h)$ we have

$$\begin{equation*}\mathrm{d}\tilde{a}(\tilde{\sigma };\hat{h})\tilde{u}={\int }_{\tilde{{\Omega}}}\hat{h}\nabla \tilde{u}\cdot \nabla \tilde{v}\enspace \mathrm{d}\tilde{\boldsymbol{x}}.\end{equation*}$$

Next,

$$\begin{equation*}\mathrm{d}B({\Omega};h)=\mathrm{exp}(\beta (s)){\left(1+{\left(\frac{\mathrm{d}f(s)}{\mathrm{d}s}\right)}^{2}{H}^{2}\right)}^{-\frac{1}{2}}{H}^{2}\frac{\mathrm{d}f(s)}{\mathrm{d}s}\frac{\mathrm{d}h(s)}{\mathrm{d}s},\end{equation*}$$

and finally, by setting $\check{h}=\mathrm{d}B({\Omega};h)$, we have

$$\begin{equation*}\mathrm{d}\tilde{a}(B;\check{h})\tilde{u}={\int }_{{\tilde{{\Gamma}}}_{\text{I}}}\check{h}\tilde{u}\tilde{v}\enspace \mathrm{d}{\tilde{\boldsymbol{s}}}_{\text{I}}.\end{equation*}$$

Notice that in the process of determining the derivatives of the forward problem with respect to perturbations of the inaccessible part of the boundary, we have also determined the derivatives of the forward problem with respect to perturbations of the Robin coefficient.

As stated previously, the inaccessible part of the domain boundary is treated as a function of x₁ only. Similarly to [34, 43], we take ${{\Gamma}}_{\text{I}}=\left\{(s,f(s)):s\in [0,L]\right\}$ with $f\in {C}_{\text{per}}^{\infty }([0,L])$ and take

$$\begin{align}\hfill {{\Gamma}}_{\text{I}}& =H\left(1+{\alpha }_{0}+\sum\limits _{n=1}^{p}{\alpha }_{2n}\mathrm{cos}\left(n2\pi \frac{{x}_{1}}{L}\right)+{\alpha }_{2n-1}\mathrm{sin}\left(n2\pi \frac{{x}_{1}}{L}\right)\right),\quad {\alpha }_{0},{\alpha }_{1},\dots .{\alpha }_{2p}\in \mathbb{R}.\hfill \end{align} \tag{ 3.2 }$$

As such, the problem of inferring the shape of the inaccessible part of the domain is posed as inferring the Fourier coefficients ⁶ , $\boldsymbol{\alpha }={[{\alpha }_{0},{\alpha }_{1},\dots ,{\alpha }_{2p}]}^{T}\in {\mathbb{R}}^{2p+1}$. As in [34, 43], we encode our prior beliefs about the geometry of the domain through a Gaussian prior on α , with mean α _* and covariance matrix Γ_α . As such, the prior density for α can be written as

$$\begin{equation*}\pi (\boldsymbol{\alpha })\propto \mathrm{exp}\left\{-\frac{1}{2}{{\Vert}\boldsymbol{\alpha }-{\boldsymbol{\alpha }}_{\ast }{\Vert}}_{{\mathbf{\Gamma }}_{\alpha }^{-1}}^{2}\right\}.\end{equation*}$$

Specification of α _* and Γ_α is left to section 4.

We postulate a Gaussian prior measure on the Robin coefficient, i.e., ${\mu }_{\beta }=\mathcal{N}({\beta }_{\ast },{\mathcal{C}}_{\beta })$. Inline with e.g. [48, 55, 56], to ensure the inverse problem is well-posed in infinite dimensions [31], while simultaneously promoting smoothness, we use a squared inverse elliptic operator to define the prior covariance operator. Specifically, we take ${\mathcal{C}}_{\beta }={\mathcal{A}}^{-1}$, where $\mathcal{A}$ is the second order elliptic differential operator defined by

$$\begin{equation}\mathcal{A}\beta {:=}-\nabla \cdot ({\delta }_{\beta }^{-2}\nabla \beta )+{l}^{2}{\delta }_{\beta }^{-2}\beta \quad \text{on}\enspace {{\Gamma}}_{\text{I}},\end{equation} \tag{ 3.3 }$$

where δ_β controls the marginal variance and l is the inverse of the correlation length. As discussed in [57–59], suitable boundary conditions should also be applied to reduce so-called boundary artefacts. To this end, letting η _I denote the outward facing unit vector normal to ∂Γ_I, we equip $\mathcal{A}$ with a Robin boundary condition,

$$\begin{equation*}{\delta }_{\beta }^{-2}\nabla \beta \cdot {\boldsymbol{\eta }}_{\text{I}}+l{\delta }_{\beta }^{-2}\beta =0\quad \text{at}\;\partial {{\Gamma}}_{\text{I}},\end{equation*}$$

which leads to a constant (i.e., homogeneous) marginal variance for β over Γ_I [58, 59].

The Robin coefficient is discretised using continuous piece-wise linear Lagrange basis functions ${\left\{{\phi }_{i}\right\}}_{i=1}^{q}$, leading to a discrete approximation for β of the form ${\beta }_{h}={\sum }_{i=1}^{q}\;{\beta }_{i}{\phi }_{i}$. Hence, the parameters of interest for the Robin coefficient are $\boldsymbol{\beta }={[{\beta }_{1},{\beta }_{2},\dots ,{\beta }_{q}]}^{T}\in {\mathbb{R}}^{q}$, and the resulting discrete representation of the prior covariance operator, denoted Γ_β , is

$$\begin{equation*}{\left[{\mathbf{\Gamma }}_{\beta }^{-1}\right]}_{ij}={\int }_{{{\Gamma}}_{\text{I}}}{\phi }_{i}(s)\mathcal{A}{\phi }_{j}(s)\mathrm{d}s\quad i,j\in \left\{1,2,\dots ,q\right\},\end{equation*}$$

see for example [48, 56, 60]. That is, the covariance matrix can be written as

$$\begin{equation}{\mathbf{\Gamma }}_{\beta }={\delta }_{\beta }^{2}{(\boldsymbol{K}+\boldsymbol{M}+\boldsymbol{R})}^{-1},\end{equation} \tag{ 3.4 }$$

where

$$\begin{equation*}{\boldsymbol{K}}_{ij}={\int }_{{{\Gamma}}_{\text{I}}}\nabla {\phi }_{i}\cdot \nabla {\phi }_{j}\enspace \mathrm{d}s,\qquad {\boldsymbol{M}}_{ij}={\int }_{{{\Gamma}}_{\text{I}}}{l}^{2}{\phi }_{i}{\phi }_{j}\enspace \mathrm{d}s,\qquad {\boldsymbol{R}}_{ij}={\int }_{\partial {{\Gamma}}_{\text{I}}}l{\phi }_{i}{\phi }_{j}\enspace \mathrm{d}t,\end{equation*}$$

for i, j = 1, 2, ..., q. The prior distribution for β can then be written as

$$\begin{equation*}\pi (\boldsymbol{\beta })\propto \mathrm{exp}\left\{-\frac{1}{2}{{\Vert}\boldsymbol{\beta }-{\boldsymbol{\beta }}_{\ast }{\Vert}}_{{\mathbf{\Gamma }}_{\beta }^{-1}}^{2}\right\}.\end{equation*}$$

In the current paper, we take the boundary shape and Robin coefficient to be (statistically) mutually independent a priori, allowing us to decompose the prior as π_prior( α , β ) = π( α )π( β ). The joint prior density is thus given by

$$\begin{equation*}{\pi }_{\text{prior}}(\boldsymbol{\alpha },\boldsymbol{\beta })=\pi (\boldsymbol{\alpha })\pi (\boldsymbol{\beta })\propto \mathrm{exp}\left\{-\frac{1}{2}{{\Vert}\boldsymbol{\alpha }-{\boldsymbol{\alpha }}_{\ast }{\Vert}}_{{\mathbf{\Gamma }}_{\alpha }^{-1}}^{2}-\frac{1}{2}{{\Vert}\boldsymbol{\beta }-{\boldsymbol{\beta }}_{\ast }{\Vert}}_{{\mathbf{\Gamma }}_{\beta }^{-1}}^{2}\right\}.\end{equation*}$$

As is fairly standard, see e.g. [7, 24, 25], we assume that the data, $\boldsymbol{y}\in {\mathbb{R}}^{m}$, is corrupted by additive noise e and is related to the parameters through

$$\begin{equation*}\boldsymbol{y}=\mathcal{G}(\boldsymbol{\alpha },\boldsymbol{\beta })+\boldsymbol{e},\end{equation*}$$

where $\mathcal{G}:{\mathbb{R}}^{2p+1}\times {\mathbb{R}}^{q}\to {\mathbb{R}}^{m}$ is termed the parameters-to-observable mapping ⁷ , and $\boldsymbol{e}\in {\mathbb{R}}^{m}$ denotes the additive noise. As discussed in section 2, we take the data to be comprised of point-wise noisy measurements of u, satisfying (2.1), at points along the accessible part of the boundary, Γ_A. As such, we can rewrite the parameters-to-observable mapping as $\mathcal{G}(\boldsymbol{\alpha },\boldsymbol{\beta })=\boldsymbol{\mathcal{B}}\boldsymbol{u}\in {\mathbb{R}}^{m}$, where u = vec(u₁, u₂, ..., u_l ), and $\boldsymbol{\mathcal{B}}$ denotes the linear observation (interpolation) operator which gives the values of each u_i at the measurement locations. In practice, each of the u_i are discretised and thus with a slight abuse of notation, we take the parameter-to-observable mapping to be of $\boldsymbol{\mathcal{B}}:{\mathbb{R}}^{r}\to {\mathbb{R}}^{m}$.

We assume that the additive noise is normally distributed, i.e., $\boldsymbol{e}\sim \mathcal{N}({\boldsymbol{e}}_{\ast },{\mathbf{\Gamma }}_{e})$, and mutually independent of the parameters. These assumptions lead to the likelihood taking the form

$$\begin{equation*}\pi (\boldsymbol{y}\vert \boldsymbol{\alpha },\boldsymbol{\beta })\propto \mathrm{exp}\left\{-\frac{1}{2}{{\Vert}\boldsymbol{y}-\boldsymbol{\mathcal{B}}\boldsymbol{u}{\Vert}}_{{\mathbf{\Gamma }}_{e}^{-1}}^{2}\right\},\end{equation*}$$

see for example [29–31].

By employing the above parameterisations for the inaccessible part of the boundary, Γ_I, and the Robin coefficient, β, the solution to the inference problem is the joint posterior density of α and β . This can be expressed using the finite dimensional version of Bayes' formula, as

$$\begin{align*}\hfill {\pi }_{\text{post}}(\boldsymbol{\alpha },\boldsymbol{\beta }\vert \boldsymbol{y})& \propto {\pi }_{\text{like}}(\boldsymbol{y}\vert \boldsymbol{\alpha },\boldsymbol{\beta }){\pi }_{\text{prior}}(\boldsymbol{\alpha },\boldsymbol{\beta })\hfill \\ \hfill & \propto \mathrm{exp}\left\{-\frac{1}{2}\left({{\Vert}\boldsymbol{y}-\boldsymbol{\mathcal{B}}\boldsymbol{u}{\Vert}}_{{\mathbf{\Gamma }}_{e}^{-1}}^{2}+{{\Vert}\boldsymbol{\alpha }-{\boldsymbol{\alpha }}_{\ast }{\Vert}}_{{\mathbf{\Gamma }}_{\alpha }^{-1}}^{2}+{{\Vert}\boldsymbol{\beta }-{\boldsymbol{\beta }}_{\ast }{\Vert}}_{{\mathbf{\Gamma }}_{\beta }^{-1}}^{2}\right)\right\}.\hfill \end{align*}$$

We define the posterior potential (of the posterior density),

$$\begin{equation}\mathcal{J}{:=}\frac{1}{2}\left({{\Vert}\boldsymbol{y}-\boldsymbol{\mathcal{B}}\boldsymbol{u}{\Vert}}_{{\mathbf{\Gamma }}_{e}^{-1}}^{2}+{{\Vert}\boldsymbol{\alpha }-{\boldsymbol{\alpha }}_{\ast }{\Vert}}_{{\mathbf{\Gamma }}_{\alpha }^{-1}}^{2}+{{\Vert}\boldsymbol{\beta }-{\boldsymbol{\beta }}_{\ast }{\Vert}}_{{\mathbf{\Gamma }}_{\beta }^{-1}}^{2}\right),\end{equation} \tag{ 3.5 }$$

so that the posterior can be written as ${\pi }_{\text{post}}(\boldsymbol{\alpha },\boldsymbol{\beta }\vert \boldsymbol{y})\propto \mathrm{exp}\left\{-\mathcal{J}\right\}$.

3.3.1. The Laplace approximation to the posterior

In the context of Bayesian inference for problems governed by partial differential equations, the most common Gaussian approximation to the posterior is the Laplace approximation, see e.g. [7, 48, 61]. Assuming the parameter-to-observable mapping is Fréchet differentiable with respect to α and β, the Laplace approximation is ${\hat{\pi }}_{\text{post}}(\boldsymbol{\alpha },\boldsymbol{\beta }\vert \boldsymbol{y})=\mathcal{N}(({\boldsymbol{\alpha }}_{\text{MAP}},{\boldsymbol{\beta }}_{\text{MAP}}),{\hat{\mathbf{\Gamma }}}_{\text{post}})$, where ( α _MAP, β _MAP) denotes the MAP estimate of ( α , β ),

$$\begin{align*}\hfill ({\boldsymbol{\alpha }}_{\text{MAP}},{\boldsymbol{\beta }}_{\text{MAP}})& {:=}\mathrm{arg}\underset{(\boldsymbol{\alpha },\boldsymbol{\beta })\in {\mathbb{R}}^{q+2p+1}}{\mathrm{max}}\enspace {\pi }_{\text{post}}(\boldsymbol{\alpha },\boldsymbol{\beta }\vert \boldsymbol{y})\hfill \\ \hfill & \propto \mathrm{arg}\underset{(\boldsymbol{\alpha },\boldsymbol{\beta })\in {\mathbb{R}}^{q+2p+1}}{\mathrm{min}}\enspace \mathcal{J},\hfill \end{align*}$$

and the approximate posterior covariance matrix is given by

$$\begin{equation}{\hat{\mathbf{\Gamma }}}_{\text{post}}={\left(\mathsf{G}{({\boldsymbol{\alpha }}_{\text{MAP}},{\boldsymbol{\beta }}_{\text{MAP}})}^{\ast }{\mathbf{\Gamma }}_{e}^{-1}\mathsf{G}({\boldsymbol{\alpha }}_{\text{MAP}},{\boldsymbol{\beta }}_{\text{MAP}})+{\mathbf{\Gamma }}_{\boldsymbol{m}}^{-1}\right)}^{-1},\end{equation} \tag{ 3.6 }$$

where Γ _m = diag(Γ_α , Γ_β ), and $\mathsf{G}({\boldsymbol{\alpha }}_{\text{MAP}},{\boldsymbol{\beta }}_{\text{MAP}})$ denotes the generalised derivative of $\mathcal{G}$ with respect to α and β evaluated at the MAP estimate, i.e.,

$$\begin{equation}\mathsf{G}({\boldsymbol{\alpha }}_{\text{MAP}},{\boldsymbol{\beta }}_{\text{MAP}})=[\mathrm{d}\mathcal{G}(({\boldsymbol{\alpha }}_{\text{MAP}},{\boldsymbol{\beta }}_{\text{MAP}});\boldsymbol{\alpha })\quad \mathrm{d}\mathcal{G}(({\boldsymbol{\alpha }}_{\text{MAP}},{\boldsymbol{\beta }}_{\text{MAP}});\beta )].\end{equation} \tag{ 3.7 }$$

It is worth noting that the approximate posterior covariance matrix coincides with the inverse of the Gauss–Newton approximation to the Hessian of the negative log-posterior (referred to simply as the Hessian in what follows), denoted H , i.e.,

$$\begin{equation*}{\hat{\mathbf{\Gamma }}}_{\text{post}}={\boldsymbol{H}}^{-1}.\end{equation*}$$

3.3.2. Full characterisation of the posterior

In some cases, the MAP estimate and a localised Gaussian approximation, such as the Laplace approximation, may be sufficient for the quantification of uncertainty in the parameters. However, without fully characterising the posterior it is essentially impossible to determine the feasibility, and thus relevance, of the local Gaussian approximation. As stated previously, fully characterising the posterior for nontrivial problems requires sampling based methods such as MCMC.

Sampling in high dimensional problems can often become infeasible in practice. However, this curse of dimensionality can often be offset to some extent by implementing sampling schemes which exploit the geometry of the posterior, see for example [62, 63]. One popular MCMC sampling scheme which takes into account the geometry of the posterior is MALA [64–66], which can be motivated using Langevin dynamics. The basic idea in MALA is to guide the proposal towards an area of higher probability, thus speeding up the convergence. Use of MALA requires the gradient of the negative log-posterior with respect to the unknown parameters.

We now briefly recall the key steps to MALA, for more details see for example [64, 65]. Langevin diffusion is defined by the following stochastic differential equation

$$\begin{equation}\mathrm{d}{\boldsymbol{m}}_{t}=-\boldsymbol{A}{\nabla }_{\boldsymbol{m}}\mathcal{J}(\boldsymbol{m})\mathrm{d}t+\sqrt{2}{\boldsymbol{A}}^{1/2}\enspace \mathrm{d}{\boldsymbol{B}}_{t},\end{equation} \tag{ 3.8 }$$

where we have denoted all the unknowns as m = ( α , β), A is a positive definite (preconditioning) matrix, and B _t is standard Brownian motion.

A key feature of the above process is that in the continuous-time case its stationary and limiting distribution is π_post. However, discrete approximations to (3.8) can have vastly different asymptotic behaviours from the diffusion process they try to approximate [64]. As such, it is necessary to introduce a Metropolis–Hastings accept/reject step that ensures the convergence to π_post. In MALA, the proposal for the next step, ${\boldsymbol{m}}_{k+1}^{\ast }$, is generated as

$$\begin{equation*}{\boldsymbol{m}}_{k+1}^{\ast }={\boldsymbol{m}}_{k}-\tau \boldsymbol{A}{\nabla }_{\boldsymbol{m}}\mathcal{J}({\boldsymbol{m}}_{k})+\sqrt{2\tau }{\boldsymbol{A}}^{1/2}\enspace \boldsymbol{\xi },\end{equation*}$$

where now the preconditioning matrix A represents the covariance of the proposal density $q({\boldsymbol{m}}_{k+1}^{\ast }\vert {\boldsymbol{m}}_{k})$, τ > 0 can be seen as a step size parameter that globally scales the proposal density, and ξ is a standard normal (q + 2p + 1-dimensional) probability density. The proposal is then accepted or rejected according to the Metropolis–Hastings algorithm, where

$$\begin{equation}\alpha =\mathrm{min}\left\{1,\frac{\pi ({\boldsymbol{m}}_{k+1}^{\ast })q({\boldsymbol{m}}_{k}\vert {\boldsymbol{m}}_{k+1}^{\ast })}{\pi ({\boldsymbol{m}}_{k})q({\boldsymbol{m}}_{k+1}^{\ast }\vert {\boldsymbol{m}}_{k})}\right\},\end{equation} \tag{ 3.9 }$$

is the probability of accepting the proposal. In practice, we work with the logarithm of α to avoid numerical underflow.

In MALA, the drift towards the posterior gradient makes the proposal density q(⋅|⋅) nonsymmetric, and therefore the ratio $q({\boldsymbol{m}}_{k}\vert {\boldsymbol{m}}_{k+1}^{\ast })/q({\boldsymbol{m}}_{k+1}^{\ast }\vert {\boldsymbol{m}}_{k})$ does not cancel out. The log-ratio of the proposal densities is given by

$$\begin{align*}\hfill \mathrm{log}\left(\frac{q({\boldsymbol{m}}_{k}\vert {\boldsymbol{m}}_{k+1}^{\ast })}{q({\boldsymbol{m}}_{k+1}^{\ast }\vert {\boldsymbol{m}}_{k})}\right)& =-\frac{1}{4\tau }\left({{\Vert}\boldsymbol{L}(-\delta \boldsymbol{m}+\tau \boldsymbol{A}{\nabla }_{\boldsymbol{m}}\mathcal{J}({\boldsymbol{m}}_{k+1}^{\ast })){\Vert}}^{2}-{{\Vert}\boldsymbol{L}(\delta \boldsymbol{m}+\tau \boldsymbol{A}{\nabla }_{\boldsymbol{m}}\mathcal{J}({\boldsymbol{m}}_{k})){\Vert}}^{2}\right),\hfill \end{align*}$$

where $\delta \boldsymbol{m}{:=}{\boldsymbol{m}}_{k+1}^{\ast }-{\boldsymbol{m}}_{k}$ is the taken step, and L ^T L = A ⁻¹.

Let us now discuss the selection of the proposal covariance A and the proposal scaling parameter τ. To achieve optimal sampling efficiency, we should have A = Σ_π , where Σ_π denotes the covariance of the true posterior [67]. In addition, it has been shown that, under various assumptions, the asymptotically optimal acceptance rate for MALA is 0.574 [65]. The acceptance rate can be adjusted using the scaling parameter τ. Since Σ_π and the optimal value of τ are unknown prior to carrying out the inversion, we adapt them continuously during the sampling, see [67–69]. A reasonable starting location for MCMC is the MAP estimate given by the Gauss–Newton iterations, and as the initial proposal covariance before the start of the adaptation we use the Laplace approximation (3.6).

Before any conclusions can be drawn from a MCMC run, we need to be reasonably sure that the chain has converged to its limiting distribution. Naturally, the convergence can never be guaranteed since any method of diagnosing convergence can by definition only diagnose the non-convergence of a chain. In practice, however, we can arrive at a reasonable certainty by testing the convergence with multiple complementary methods. These methods include visual inspection of the chains and measuring the sampling quality by autocorrelation, comparing within-chain variances between multiple runs (the Gelman–Rubin diagnostic) [70, 71], and estimating sampling uncertainty by computing the Monte Carlo standard error (MCSE) [72, 73].

In this paper, we implement a stopping rule based on the estimated MCSE and its relation to the posterior variance of the parameters, as in [74]. First, we run the sampler for enough steps so that the MCSE can be reliably calculated using the non-overlapping batch means method, see [72]. Then, at regular intervals we check, for each unknown individually, if the MCSE estimate at 98% confidence level is smaller than 10% of the posterior standard deviation. When this is true for each unknown, we terminate the run.

Calculation of the MAP estimate, ( α _MAP, β _MAP), and construction of the Laplace approximation of the posterior covariance matrix, ${\hat{\mathbf{\Gamma }}}_{\text{post}}$, requires the generalised derivative $\mathsf{G}$ (see (3.7)), details of which were provided in section 2.2. On the other hand, the use of MALA necessitates computation of only the gradient of the negative log-posterior with respect to both the inaccessible part of the boundary and the Robin coefficient. Although the gradient is essentially trivial to calculate with the Jacobian in hand, computing the Jacobian can be costly, requiring as many forward solves as the number of measurements. To avoid this potential bottleneck, we employ the so-called adjoint approach [39–41], which computes the gradient at a cost of one additional (adjoint) forward solve. In what follows, we will set g = g_k and u = u_k for some fixed $k\in \left\{1,2,\dots ,\ell \right\}$. We begin by introducing a Lagrangian, $\mathcal{L}:{H}_{D}^{1}({\Omega})\times {H}_{D}^{1}({\Omega})\times {\mathbb{R}}^{2p+1}\times {H}^{1}([0,L])\to \mathbb{R}$, as

$$\begin{align*}\hfill \mathcal{L}(u,v,\boldsymbol{\alpha },\beta )& =\mathcal{J}+{\int }_{{\Omega}(\boldsymbol{\alpha })}\sigma \nabla u\cdot \nabla v\enspace \mathrm{d}\enspace \boldsymbol{x}+{\int }_{{{\Gamma}}_{\text{I}}(\boldsymbol{\alpha })}\enspace \mathrm{exp}(\beta )uv\enspace \mathrm{d}\enspace {\boldsymbol{s}}_{\text{I}}-{\int }_{{{\Gamma}}_{\text{A}}}gv\enspace \mathrm{d}\enspace {\boldsymbol{s}}_{\text{A}}\hfill \\ \hfill & =\mathcal{J}+{\int }_{\tilde{{\Omega}}}{\tilde{\sigma }}_{\boldsymbol{\alpha }}\nabla \tilde{u}\cdot \nabla \tilde{v}\enspace \mathrm{d}\tilde{\boldsymbol{x}}+{\int }_{{\tilde{{\Gamma}}}_{\text{I}}}\enspace \mathrm{exp}({\tilde{\beta }}_{\boldsymbol{\alpha }})\tilde{u}\tilde{v}\enspace \mathrm{d}{\tilde{\boldsymbol{s}}}_{\text{I}}-{\int }_{{{\Gamma}}_{\text{A}}}g\tilde{v}\enspace \mathrm{d}\enspace {\boldsymbol{s}}_{\text{A}},\hfill \end{align*}$$

with $\mathcal{J}$ defined in (3.5), and where we have used the invariance property of the Poisson equation. The subscript α is used to denote the dependence on α .

Determining the gradient of $\mathcal{J}$ (with respect to the unknowns) is achieved by requiring that variations of the Lagrangian with respect to the forward potential, u, and the so-called adjoint potential, v, are 0, see e.g. [41]. The variation of the Lagrangian with respect to the adjoint potential gives

$$\begin{align*}\hfill \mathrm{d}\mathcal{L}((u,v,\boldsymbol{\alpha },\beta );\tilde{w})& =\mathcal{J}+{\int }_{\tilde{{\Omega}}}{\tilde{\sigma }}_{\boldsymbol{\alpha }}\nabla \tilde{u}\cdot \nabla \tilde{w}\enspace \mathrm{d}\tilde{\boldsymbol{x}}+{\int }_{{\tilde{{\Gamma}}}_{\text{I}}}\enspace \mathrm{exp}({\tilde{\beta }}_{\boldsymbol{\alpha }})\tilde{u}\tilde{w}\enspace \mathrm{d}{\tilde{\boldsymbol{s}}}_{\text{I}}-{\int }_{{{\Gamma}}_{\text{A}}}g\tilde{w}\enspace \mathrm{d}\enspace {\boldsymbol{s}}_{\text{A}}=0,\hfill \end{align*}$$

for $\tilde{w}\in {H}_{D}^{1}(\tilde{{\Omega}})$, i.e., u must satisfy the forward problem. The variation with respect to the forward potential gives the adjoint equation,

$$\begin{align*}\hfill \mathrm{d}\mathcal{L}((u,v,\boldsymbol{\alpha },\beta );\tilde{z})& =-{\boldsymbol{\mathcal{B}}}^{\ast }{\mathbf{\Gamma }}_{e}^{-1}(\boldsymbol{y}-\boldsymbol{\mathcal{B}}\tilde{z})+{\int }_{\tilde{{\Omega}}}{\tilde{\sigma }}_{\boldsymbol{\alpha }}\nabla \tilde{z}\cdot \nabla \tilde{v}\enspace \mathrm{d}\tilde{\boldsymbol{x}}+{\int }_{{\tilde{{\Gamma}}}_{\text{I}}}\enspace \mathrm{exp}({\tilde{\beta }}_{\boldsymbol{\alpha }})\tilde{z}\tilde{v}\enspace \mathrm{d}\tilde{{\boldsymbol{s}}_{\text{I}}}=0\hfill \end{align*}$$

for $\tilde{z}\in {H}_{D}^{1}(\tilde{{\Omega}})$, which the adjoint variable v must satisfy (this is the one additional forward solve required). Finally, the gradients with respect to the Fourier coefficients and the Robin coefficient are then given (in strong form) as

$$\begin{equation*}{\nabla }_{\boldsymbol{\alpha }}\mathcal{L}(u,v,\boldsymbol{\alpha },\beta )={\mathbf{\Gamma }}_{\alpha }^{-1}(\boldsymbol{\alpha }-{\boldsymbol{\alpha }}_{\ast })+{\int }_{\tilde{{\Omega}}}{\boldsymbol{\kappa }}_{\boldsymbol{\alpha }}\nabla \tilde{u}\cdot \nabla \tilde{v}\enspace \mathrm{d}\tilde{\boldsymbol{x}}+{\int }_{{\tilde{{\Gamma}}}_{\text{I}}}{\boldsymbol{b}}_{\boldsymbol{\alpha }}\tilde{u}\tilde{v}\enspace \mathrm{d}{\tilde{\boldsymbol{s}}}_{\text{I}},\end{equation*}$$

and

$$\begin{equation*}{\nabla }_{\beta }\mathcal{L}(u,v,\boldsymbol{\alpha },\beta )={\mathcal{A}}_{\beta }({\beta }_{\boldsymbol{\alpha }}-{\beta }_{\ast })+\mathrm{exp}({\beta }_{\boldsymbol{\alpha }})uv,\end{equation*}$$

respectively, where

where $\boldsymbol{x}={\psi }^{-1}(\tilde{\boldsymbol{x}})$ and $s={\varphi }^{-1}(\tilde{s})$, and i = 0, 1, ..., 2p.

Finally, the gradient of the cost function evaluated at (u, v, α , β) is then given by

$$\begin{equation}\nabla \mathcal{J}=\left(\begin{matrix}\hfill {\nabla }_{\boldsymbol{\alpha }}\mathcal{L}(u,v,\boldsymbol{\alpha },\beta )\hfill \\ \hfill {\nabla }_{\beta }\mathcal{L}(u,v,\boldsymbol{\alpha },\beta )\hfill \end{matrix}\right).\end{equation} \tag{ 3.10 }$$

For each example considered, we take the diffeomorphism (from the true domain Ω to $\tilde{{\Omega}}$) to be of the form

$$\begin{equation*}\psi :({x}_{1},{x}_{2}){\mapsto}\left({x}_{1},{x}_{2}{(f({x}_{1}))}^{-1}\right),\quad ({x}_{1},{x}_{2})\in {\Omega},\end{equation*}$$

with

$$\begin{equation*}f({x}_{1})=1+{\alpha }_{0}+\sum\limits _{n=1}^{p}{\alpha }_{2n}\mathrm{cos}(2\pi n{x}_{1})+{\alpha }_{2n-1}\mathrm{sin}(2\pi n{x}_{1}).\end{equation*}$$

The choice of p can be viewed as making a trade-off between computational speed and the maximum possible accuracy of the solution. That is, with a smaller p we have fewer parameters to infer but smaller values of p cannot represent as complex shapes. In the current paper we take p = 7, meaning we have 15 Fourier coefficients to estimate.

For examples 2 and 3, representing the true boundary shapes accurately requires substantially more basis functions than the 15 used. Thus, to better assess the estimates, we can compare them to the true boundary shape projected onto the truncated Fourier basis. This is illustrated in figure 3, where we show the reference inversion mesh (used in all examples) as well as the inversion mesh corresponding to the true domain for example 3 using the first 15 Fourier basis functions. It is clear that steeper changes in the boundary shape cannot be represented by the truncated Fourier series.

Zoom In Zoom Out Reset image size

We use the same prior distribution for all examples. Specifically, for the boundary shape we take α _* = 0 and

$$\begin{equation}{\mathbf{\Gamma }}_{\alpha }=\mathrm{diag}({\sigma }_{0}^{2},{\sigma }_{1}^{2},{\sigma }_{1}^{2},\dots ,{\sigma }_{p}^{2},{\sigma }_{p}^{2})={\sigma }_{\alpha }^{2}\mathrm{diag}(1,{2}^{{s}_{\alpha }},{2}^{{s}_{\alpha }},\dots ,{(p+1)}^{{s}_{\alpha }},{(p+1)}^{{s}_{\alpha }}),\end{equation} \tag{ 4.1 }$$

where ${\sigma }_{\alpha }^{2}$ controls the overall variance of the parameters, and s_α < 0 controls how quickly the variance diminishes towards the higher frequency components. The faster the Fourier coefficients converge towards zero, the smoother the resulting function is [76], and as such, s_α can be interpreted as a smoothness parameter. We choose ${\sigma }_{\alpha }^{2}=0.01$ and s_α = − 1 for all considered examples. For the Robin coefficient we set β _* = 0 and Γ_β as in (3.4), where we set ${\delta }_{\beta }^{2}=50$ and l = 10. Draws from both priors are shown in figure 4.

Zoom In Zoom Out Reset image size

The meshes used for the generation of the synthetic data are different in each case (see figure 2), and each mesh has between 2000 and 3000 nodes. On the other hand, for all examples, in the data-generation mesh the Robin coefficient is represented using 230 finite element basis functions. The discretisation of the mesh used for inference has a total of 640 nodes, and results in 78 degrees of freedom for the Robin coefficient, i.e., $\boldsymbol{\beta }\in {\mathbb{R}}^{78}$. Thus, the total number of the parameters to be inferred is 78 + 15 = 93.

Finally, we remark that, inline with e.g. [77], the Gauss–Newton method used to compute the MAP estimate is terminated when the norm of the gradient has decreased by a factor of 10⁵, while line search is carried out using the Armijo condition with c₁ = 10⁻⁴ and c₂ = 0.9, see [33]. All computations were carried out using MATLAB 2018b on a laptop with an Intel i7-8850H CPU and 16 GB RAM.

Here we discuss and compare the Laplace approximation of the posterior and the accurate posterior found using MCMC for each of the examples outlined above.

Example 1. For the first example, the MAP estimate was found after 24 Gauss–Newton iterations, while the MCMC convergence criterion was reached in approximately 70 000 steps after burn-in (10 000 steps in all examples), which took approximately two hours. The MAP estimate, Laplace approximation of the (marginal) posterior, and the conditional mean and the (marginal) posterior found using MCMC are all shown in figure 5.

Both the Laplace approximation of the posterior and accurate posterior support the true boundary shape and Robin coefficient well. Comparing the Laplace approximation and accurate posterior for the boundary shape, it is clear that the accurate posterior has significantly narrower confidence intervals, particularly towards the lateral boundaries. Moreover, for the boundary shape, the CM and MAP estimates are in fairly good agreement, while the accurate posterior seems to be symmetric, indicating that the posterior is well approximated by a normal distribution. On the other hand, when comparing the Laplace approximation and accurate posterior for the Robin coefficient, we see that the accurate posterior is negatively skewed (evident from the confidence intervals and the fact that the CM estimate is generally smaller than the MAP estimate).

Zoom In Zoom Out Reset image size

Example 2. For the second example, calculating the MAP estimate took 30 Gauss–Newton iterations, while the MCMC convergence criterion was reached in approximately 140 000 steps after burn-in, which took approximately three and a half hours. In figure 6 we show the MAP estimate, Laplace approximation of the posterior and the posterior and conditional mean found using MCMC.

As in example 1, the Laplace approximation of the posterior and accurate posterior support the true boundary shape and Robin coefficient well. It also clear that for both the boundary shape and Robin coefficient, estimates are fairly insensitive to the small, highly oscillatory, variations in corresponding truths. Similarly to example 1, the Laplace approximation has slightly wider confidence intervals compared to the accurate posterior, particularly towards the lateral boundaries, and the CM and MAP estimates are in good agreement for the boundary shape. Again, the accurate posterior for the Robin coefficient is negatively skewed.

Zoom In Zoom Out Reset image size

Example 3. For the final example 31 Gauss–Newton iterations were needed to compute the MAP estimate, while the MCMC convergence criterion was reached in approximately 80 000 steps after burn-in, taking approximately two hours. Figure 7 shows the MAP estimate, Laplace approximation of the posterior and the posterior and CM found using MCMC. Since the true boundary shape lies well outside the prior distribution, we also show the projection of the true boundary shape onto the first 15 Fourier basis functions. This provides a reference for the best possible reconstruction.

Zoom In Zoom Out Reset image size

As in the previous examples, the Laplace approximation to the posterior and the accurate posterior for the Robin coefficient support the true coefficient well. Conversely, though as expected, the sharp jumps in the boundary shape are poorly supported by both the approximate and accurate posterior densities as these are impossible to estimate using the given prior on the boundary shape. However, we do see that both posterior densities support the projection of the true boundary shape onto the first 15 Fourier basis functions fairly well. It is also evident that the uncertainty in the boundary shape is significantly lower in regions where the domain is narrowest. This is inline with our intuition as the data (collected on the bottom of the slab) is much more informative on the shape of the boundary which is closer to where the data is measured. Finally, the accurate posterior does show some improvement over the approximate posterior for the boundary shape, particularly at the sharp jumps, as well as being slightly positively skewed.