Solving large-scale PDE-constrained Bayesian inverse problems with Riemann manifold Hamiltonian Monte Carlo

Inverse problems are ubiquitous in science and engineering. Perhaps the most popular family of inverse problems is to determine a set of parameters (or a function) given a set of indirect observations, which are in turn provided by a parameter-to-observable map plus observation uncertainties. For example, if one considers the problem of determining the heat conductivity of a thermal fin given measured temperature at a few locations on the thermal fin, then: (i) the desired unknown parameter is the distributed heat conductivity, (ii) the observations are the measured temperatures, (iii) the parameter-to-observable map is the mathematical model that describes the temperature on the thermal fin as a function of the heat conductivity; indeed the temperature distribution is a solution of an elliptic partial differential equation (PDE) whose coefficient is the heat conductivity, and (iv) the observation uncertainty is due to the imperfection of the measurement device and/or model inadequacy.

The Bayesian inversion framework refers to a mathematical method that allows one to solve statistical inverse problems taking into account all uncertainties in a systematic and coherent manner. The Bayesian approach does this by reformulating the inverse problem as a problem in statistical inference, incorporating uncertainties in the observations, the parameter-to-observable map, and prior information on the parameter. In particular, we seek a statistical description of all possible (set of) parameters that conform to the available prior knowledge and at the same time are consistent with the observations. The solution of the Bayesian framework is the so-called posterior measure that encodes the degree of confidence on each set of parameters as the solution to the inverse problem under consideration.

Mathematically the posterior is a surface in high dimensional parameter space. The task at hand is therefore to explore the posterior by, for example, characterizing the mean, the covariance, and/or higher moments. The nature of this task is to compute high dimensional integrals for which most contemporary methods are intractable. Perhaps the most general method to attack these problems is the Markov chain Monte Carlo (MCMC) method which shall be introduced in subsequent sections.

Let us now summarize the content of the paper. We start with the description of the statistical inverse problem under consideration in section 2. It is an inverse steady state heat conduction governed by elliptic PDEs. We postulate a Gaussian measure prior on the parameter space to ensure that the inverse problem is well-defined. The prior itself is a well-defined object whose covariance operator is the inverse of an elliptic differential operator and with the mean function living in the Cameron–Martin space of the covariance. The posterior is given by its Radon–Nikodym derivative with respect to the prior measure, which is proportional to the likelihood. Since the Riemann manifold Hamiltonian Monte Carlo (RMHMC) simulation method requires the gradient, Hessian, and the derivative of the Fisher information operator, we discuss, in some depth, how to compute the derivatives of the potential function (the misfit functional) with PDE constraints efficiently using the adjoint technique in section 3. In particular, we define a Fisher information operator and show that it coincides with the well-known Gauss–Newton Hessian of the misfit. We next present a discretization scheme for the infinite Bayesian inverse problem in section 4. Specifically, we employ a standard continuous H¹-conforming finite element method (FEM) to discretize both the likelihood and the Gaussian prior. We choose to numerically compute the truncated Karhunen–Loève (KL) expansion which requires one to solve an eigenvalue problem with fractional Laplacian. In order to accomplish this task, we use a matrix transfer technique (MTT) which leads to a natural discretization of the Gaussian prior measure. In section 5, we describe the RMHMC and its variants at length, and its application to our Bayesian inverse problem. Section 6 presents a low rank approach to approximate the Fisher information matrix and its inverse efficiently. This is possibly due to the fact that the Gauss–Newton Hessian, and hence the Fisher information operator, is a compact operator. Various numerical results supporting our proposed approach are presented in section 7. We begin this section with an extensive study and comparison of Riemannian manifold MCMC methods for problems with two parameters, and end the section with 1025-parameter problem. Finally, we conclude the paper in section 8 with a discussion on future work.

In order to clearly illustrate the challenges arising in PDE-constrained inverse problems for MCMC based Bayesian inference, we consider the following heat conduction problem governed by an elliptic PDE in the open and bounded domain $\Omega \subset {{\mathbb{R}}^{n}}$

$$\begin{eqnarray*}&&\begin{array}{l} -\nabla \cdot \left( {{{\rm e}}^{u}}\nabla w \right)=\;0\qquad \quad {\rm in}\ \Omega , \\ -{{{\rm e}}^{u}}\nabla w\cdot {\bf n}=\;Bi\;u\qquad \ \ {\rm on}\ \partial \Omega \backslash {{\Gamma }_{R}}, \\ -{{{\rm e}}^{u}}\nabla w\cdot {\bf n}=-1\qquad \quad \ {\mkern 1mu} {\rm on}\ {{\Gamma }_{R}}, \\ \end{array}\end{eqnarray*}$$

where $w$ is the forward state, $u$ the logarithm of distributed thermal conductivity on Ω, ${\bf n}$ the unit outward normal on ∂Ω, and Bi the Biot number.

In the forward problem, the task is to solve for the temperature distribution $w$ given a description of distributed parameter $u$. In the inverse problem, the task is to reconstruct $u$ given some available observations, e.g, temperature observed at some parts/locations of the domain Ω. We initially choose to cast the inverse problem in the framework of PDE-constrained optimization. To begin, let us consider the following additive noise-corrupted pointwise observation model³

$$\begin{eqnarray}&&{{d}_{j}}:=\;w\left( {{{\bf x}}_{j}} \right)+{{\eta }_{j}},\quad j=1,\ldots ,K,\end{eqnarray} \tag{ 1 }$$

where K is the total number of observation locations, $\left\{ {{{\bf x}}_{j}} \right\}_{j=1}^{K}$ the set of points at which $w$ is observed, η_j the additive noise, and ${{d}_{j}}$ the actual noise-corrupted observations. In this paper we work with synthetic observations and hence there is no model inadequacy in (1). Concatenating all the observations, one can rewrite (1) as

$$\begin{eqnarray}&&{\bf d}:=\mathcal{G}\left( u \right)+\eta ,\end{eqnarray} \tag{ 2 }$$

with $\mathcal{G}:={{[w\left( {{{\bf x}}_{1}} \right),\ldots ,w\left( {{{\bf x}}_{K}} \right)]}^{T}}$ denoting the map from the distributed parameter $u$ to the noise-free observables, $\eta$ being random numbers normally distributed by $\mathcal{N}\left( 0,{\bf L} \right)$ with bounded covariance matrix ${\bf L},$ and ${\bf d}={{[{{d}_{1}},\ldots ,{{d}_{K}}]}^{T}}$. For simplicity, we take ${\bf L}={{\sigma }^{2}}{\bf I}$, where ${\bf I}$ is the identity matrix.

Our inverse problem can be now formulated as

$$\begin{eqnarray}&&\mathop{{\rm min} }\limits_{u}\mathcal{J}\left( u,{\bf d} \right):=\frac{1}{2}\left| {\bf d}-\mathcal{G}\left( u \right) \right|_{{\bf L}}^{2}=\frac{1}{2{{\sigma }^{2}}}\mathop{\sum }\limits_{j=1}^{K}{{\left( w\left( {{{\bf x}}_{j}} \right)-{{d}_{j}} \right)}^{2}},\end{eqnarray} \tag{ 3 }$$

subject to

$$\begin{eqnarray}&&-\nabla \cdot \left( {{{\rm e}}^{u}}\nabla w \right)=0\qquad \ {\rm in}\ \Omega ,\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&-{{{\rm e}}^{u}}\nabla w\cdot {\bf n}=Bi\;u\qquad \ {\rm on}\ \partial \Omega \backslash{{\Gamma }_{R}},\end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray}&&-{{{\rm e}}^{u}}\nabla w\cdot {\bf n}=-1\qquad \ {\rm on}\ {{\Gamma }_{R}},\end{eqnarray} \tag{ 4c }$$

where ${{\left| \cdot \right|}_{{\bf L}}}:=\mid {{{\bf L}}^{-\frac{1}{2}}}\cdot \mid$ denotes the weighted Euclidean norm induced by the canonical inner product $\left( \cdot ,\cdot \right)$ in ${{\mathbb{R}}^{K}}$. This optimization problem is however ill-posed. An intuitive reason is that the dimension of observations ${\bf d}$ is much smaller than that of the parameter $u$, and hence they provide limited information about the distributed parameter $u$. As a result, the null space of the Jacobian of the parameter-to-observation map $\mathcal{G}$ is non-empty. Indeed, we have shown that the Gauss–Newton approximation of the Hessian (which is the square of this Jacobian, and is also equal to the full Hessian of the data misfit $\mathcal{J}$ evaluated at the optimal parameter) is a compact operator [1–3], and hence its range space is effectively finite-dimensional.

One way to overcome the ill-posedness is to use Tikhonov regularization (see, e.g., [4]), which proposes to augment the cost functional (3) with a quadratic term, i.e.

$$\begin{eqnarray}&&\tilde{\mathcal{J}}:=\frac{1}{2}\left| {\bf d}-\mathcal{G}\left( u \right) \right|_{{\bf L}}^{2}+\frac{\kappa }{2}{\left\|{{R}^{1/2}}u\right\|}^{2},\end{eqnarray} \tag{ 5 }$$

where κ is a regularization parameter, R some regularization operator, and $\left\|\cdot \right\|$ some appropriate norm. This method is a representative of deterministic inverse solution techniques that typically do not take into account the randomness due to measurements and other sources, though one can equip the deterministic solution with a confidence region by post-processing (see, e.g., [5] and references therein). It should be pointed out that if the regularization term is replaced by the Cameron–Martin norm of $u$ (the second term in (7)), the Tikhonov solution is in fact identical to the maximum a posteriori (MAP) point in (7). However, such a point estimate is insufficient for the purpose of fully taking the randomness into account.

In this paper, we choose to tackle the ill-posedness using a Bayesian framework [6–10]. We seek a statistical description of all possible $u$ that conform to some prior knowledge and at the same time are consistent with the observations. The Bayesian approach does this by reformulating the inverse problem as a problem in statistical inference, incorporating uncertainties in the observations, the forward map $\mathcal{G}$, and prior information. This approach is appealing since it can incorporate most, if not all, kinds of randomness in a systematic manner. To begin, we postulate a Gaussian measure $\mu :=\mathcal{N}\left( {{u}_{0}},{{\alpha }^{-1}}\mathcal{C} \right)$ on $u$ in ${{L}^{2}}\left( \Omega \right)$ where

$$\begin{eqnarray*}&&\mathcal{C}:=\;{{\left( I-\Delta \right)}^{-s}}=:{{\mathcal{A}}^{-s}},\end{eqnarray*}$$

with the domain of definition

$$\begin{eqnarray*}&&D\left( \mathcal{A} \right):=\left\{ u\in {{H}^{2}}\left( \Omega \right):\frac{\partial u}{\partial {\bf n}}=0\ {\rm on}\ \partial \Omega \right\},\end{eqnarray*}$$

where ${{H}^{2}}\left( \Omega \right)$ is the usual Sobolev space. Assume that the mean function ${{u}_{0}}$ lives in the Cameron–Martin space of $\mathcal{C}$, then one can show (see [9]) that the measure μ is well-defined when $s\gt n/2$ (d is the spatial dimension), and in that case, any realization from the prior distribution μ is almost surely in the Hölder space $X:={{C}^{0,\;\beta }}\left( \Omega \right)$ with $0\lt \beta \lt s/2$. That is, $\mu \left( X \right)=1$, and the Bayesian posterior measure ν satisfies the Radon–Nikodym derivative

$$\begin{eqnarray}&&\frac{\partial \nu }{\partial \mu }\left( u|{\bf d} \right)\sim {\rm exp} \left( -\mathcal{J}\left( u,{\bf d} \right) \right)={\rm exp} \left( -\frac{1}{2}\left| {\bf d}-\mathcal{G}\left( u \right) \right|_{{\bf L}}^{2} \right),\end{eqnarray} \tag{ 6 }$$

if $\mathcal{G}$ is a continuous map from X to ${{\mathbb{R}}^{K}}$. Note that the Radon–Nikodym derivative is proportional to the the likelihood defined by

$$\begin{eqnarray*}&&{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)\sim {\rm exp} \left( -\mathcal{J}\left( u,{\bf d} \right) \right).\end{eqnarray*}$$

The MAP point is defined as

$$\begin{eqnarray}&&{{u}^{{\rm MAP}}}:={\rm arg} \mathop{{\rm min} }\limits_{u}\mathcal{J}\left( u,{\bf d} \right):=\frac{1}{2}\left| {\bf d}-\mathcal{G}\left( u \right) \right|_{{\bf L}}^{2}+\frac{\alpha }{2}\left\|u\right\|_{\mathcal{C}}^{2},\end{eqnarray} \tag{ 7 }$$

where ${\left\|\cdot \right\|}_{\mathcal{C}}:=\parallel {{\mathcal{C}}^{-\frac{1}{2}}}\cdot \parallel$ denotes the weighted ${{L}^{2}}\left( \Omega \right)$ norm induced by the ${{L}^{2}}\left( \Omega \right)$ inner product $\left\langle \cdot ,\cdot \right\rangle$.

In this section, we briefly present the adjoint method to efficiently compute the gradient, Hessian, and the third derivative of the cost functional (3). We start by considering the weak form of the (first order) forward equation (4):

$$\begin{eqnarray}&&{{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla {\rm w}\cdot \nabla \hat{\lambda }\;{\rm d}\Omega +{{\int }_{\partial \Omega\backslash {{\Gamma }_{R}}}}Bi\;{\rm w}\hat{\lambda }\;{\rm ds}={{\int }_{{{\Gamma }_{R}}}}\hat{\lambda }\;{\rm ds},\end{eqnarray} \tag{ 8 }$$

with $\hat{\lambda }$ as the test function. Using the standard reduced space approach (see, e.g., a general discussion in [11] and a detailed derivation in [12]) one can show that the gradient $\nabla \mathcal{J}\left( u \right)$, namely the Fréchet derivative of the cost functional $\mathcal{J}$, acting in any direction $\tilde{u}$ is given by

$$\begin{eqnarray}&&\left\langle \nabla \mathcal{J}\left( u \right),\tilde{u} \right\rangle ={{\int }_{\Omega }}\tilde{u}{{{\rm e}}^{u}}\nabla w\cdot \nabla \lambda \;{\rm d}\Omega ,\end{eqnarray} \tag{ 9 }$$

where the (first order) adjoint state λ satisfies the adjoint equation

$$\begin{eqnarray}&&{{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla \lambda \cdot \nabla \hat{w}\;{\rm d}\Omega +{{\int }_{\partial \Omega \backslash{{\Gamma }_{R}}}}Bi\;\lambda \hat{w}\;{\rm d}s=-\frac{1}{{{\sigma }^{2}}}\mathop{\sum }\limits_{j=1}^{K}\left( w\left( {{{\bf x}}_{j}} \right)-{{d}_{j}} \right)\hat{w}\left( {{{\bf x}}_{j}} \right),\end{eqnarray} \tag{ 10 }$$

with $\hat{w}$ as the test function. On the other hand, the Hessian, the Fréchet derivative of the gradient, acting in directions $\tilde{u}$ and ${{u}^{2}}$ (superscript '2' means the second variation direction) reads

$$\begin{eqnarray}&&\begin{array}{r} \left\langle \left\langle {{\nabla }^{2}}\mathcal{J}\left( u \right),\tilde{u} \right\rangle ,{{u}^{2}} \right\rangle ={{\int }_{\Omega }}\tilde{u}{{{\rm e}}^{u}}\nabla w\cdot \nabla {{\lambda }^{2}}\;{\rm d}\Omega +{{\int }_{\Omega }}\tilde{u}{{{\rm e}}^{u}}\nabla {{w}^{2}}\cdot \nabla \lambda \;{\rm d}\Omega \qquad \qquad \qquad \quad \, \\ +{{\int }_{\Omega }}\tilde{u}{{u}^{2}}{{{\rm e}}^{u}}\nabla w\cdot \nabla \lambda \;{\rm d}\Omega ,\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \\ \end{array}\end{eqnarray} \tag{ 11 }$$

where the second order forward state ${{w}^{2}}$ obeys the second order forward equation

$$\begin{eqnarray}&&{{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla {{w}^{2}}\cdot \nabla \hat{\lambda }\;{\rm d}\Omega +{{\int }_{\partial \Omega \backslash {{\Gamma }_{R}}}}Bi\;{{w}^{2}}\hat{\lambda }\;{\rm d}s=-{{\int }_{\Omega }}{{u}^{2}}{{{\rm e}}^{u}}\nabla w\cdot \nabla \hat{\lambda }\;{\rm d}\Omega ,\end{eqnarray} \tag{ 12 }$$

and the second order adjoint state ${{\lambda }^{2}}$ is governed by the second order adjoint equation

$$\begin{eqnarray}\begin{array}{rcl} {{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla {{\lambda }^{2}}\cdot \nabla \hat{w}\;{\rm d}\Omega +{{\int }_{\partial \Omega\backslash {{\Gamma }_{R}}}}Bi\;{{\lambda }^{2}}\hat{w}\;{\rm d}s & = & -\frac{1}{{{\sigma }^{2}}}\mathop{\sum }\limits_{j=1}^{K}{{w}^{2}}\left( {{{\bf x}}_{j}} \right)\hat{w}\left( {{{\bf x}}_{j}} \right)-{{\int }_{\Omega }}{{u}^{2}}{{{\rm e}}^{u}}\nabla \lambda \cdot \nabla \hat{w}\;{\rm d}\Omega . \\ \end{array}\end{eqnarray} \tag{ 13 }$$

We define the generalized Fisher information operator⁴ acting in directions $\tilde{u}$ and ${{u}^{2}}$ as

$$\begin{eqnarray}&&\left\langle \left\langle G\left( u \right),\tilde{u} \right\rangle ,{{u}^{2}} \right\rangle :={{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[\left\langle \left\langle {{\nabla }^{2}}\mathcal{J}\left( u \right),\tilde{u} \right\rangle ,{{u}^{2}} \right\rangle ],\end{eqnarray} \tag{ 14 }$$

where the expectation is taken with respect to the likelihood—the distribution of the observation ${\bf d}$. Now, substituting (11) into (14) and assuming that the integrals/derivatives can be interchanged we obtain

$$\begin{eqnarray*}\begin{array}{rcl} \left\langle \left\langle G\left( u \right),\tilde{u} \right\rangle ,{{u}^{2}} \right\rangle & = & {{\int }_{\Omega }}\tilde{u}{{{\rm e}}^{u}}\nabla w\cdot \nabla {{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[{{\lambda }^{2}}]\;{\rm d}\Omega +{{\int }_{\Omega }}\tilde{u}{{{\rm e}}^{u}}\nabla {{w}^{2}}\cdot \nabla {{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[\lambda ]\;{\rm d}\Omega \\ {} & {} & +{{\int }_{\Omega }}\tilde{u}{{u}^{2}}{{{\rm e}}^{u}}\nabla w\cdot \nabla {{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[\lambda ]\;{\rm d}\Omega , \\ \end{array}\end{eqnarray*}$$

where we have used the assumption that the parameter $u$ is independent of observation ${\bf d}$ and the fact that $w$ and ${{w}^{2}}$ do not depend on ${\bf d}$. The next step is to compute $\nabla {{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[\lambda ]$ and ${{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[{{\lambda }^{2}}]$. To begin, let us take the expectation the first order adjoint equation (10) with respect to ${{\pi }_{{\rm like}}}\left( {\bf d}|u \right)$ to arrive at

$$\begin{eqnarray*}\begin{array}{ccccccccccccccc} {{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla {{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[\lambda ]\cdot \nabla \hat{w}\;{\rm d}\Omega & + & {{\int }_{\partial \Omega \backslash{{\Gamma }_{R}}}}Bi\;{{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[\lambda ]\hat{w}\;{\rm d}s\;= \\ {} & - & \frac{1}{{{\sigma }^{2}}}\mathop{\sum }\limits_{j=1}^{K}{{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[w\left( {{{\bf x}}_{j}} \right)-{{d}_{j}}]\hat{w}\left( {{{\bf x}}_{j}} \right)=0, \\ \end{array}\qquad \qquad \ \ \ {\mkern 1mu} \end{eqnarray*}$$

where the second equality is obtained from (1) and the assumption ${{\eta }_{j}}\sim \mathcal{N}\left( 0,{{\sigma }^{2}} \right)$. We conclude that

$$\begin{eqnarray}&&\nabla {{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[\lambda ]=0.\end{eqnarray} \tag{ 15 }$$

On the other hand, if we take the expectation of the second order adjoint equation (13) and use (15) we have

$$\begin{eqnarray}&&\begin{array}{r} {{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla {{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[{{\lambda }^{2}}]\cdot \nabla \hat{w}\;{\rm d}\Omega +{{\int }_{\partial \Omega\backslash {{\Gamma }_{R}}}}Bi\;{{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[{{\lambda }^{2}}]\hat{w}\;{\rm d}s\qquad \qquad \qquad \qquad \quad \ \ \ {\mkern 1mu} \\ =-\frac{1}{{{\sigma }^{2}}}\mathop{\sum }\limits_{j=1}^{K}{{w}^{2}}\left( {{{\bf x}}_{j}} \right)\hat{w}\left( {{{\bf x}}_{j}} \right).\qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \ \ \ {\mkern 1mu} \\ \end{array}\end{eqnarray} \tag{ 16 }$$

Let us define

$$\begin{eqnarray*}&&\tilde{{{\lambda }^{2}}}:={{\mathbb{E}}_{{{\pi }_{{\rm like}}}\left( {\bf d}|u \right)}}[{{\lambda }^{2}}],\end{eqnarray*}$$

then (16) becomes

$$\begin{eqnarray}&&{{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla \tilde{{{\lambda }^{2}}}\cdot \nabla \hat{w}\;{\rm d}\Omega +{{\int }_{\partial \Omega \backslash{{\Gamma }_{R}}}}Bi\;\tilde{{{\lambda }^{2}}}\hat{w}\;{\rm d}s=-\frac{1}{{{\sigma }^{2}}}\mathop{\sum }\limits_{j=1}^{K}{{w}^{2}}\left( {{{\bf x}}_{j}} \right)\hat{w}\left( {{{\bf x}}_{j}} \right).\end{eqnarray} \tag{ 17 }$$

As a result, the Fisher information operator acting along directions $\tilde{u}$ and ${{u}^{2}}$ reads

$$\begin{eqnarray}&&\left\langle \left\langle G\left( u \right),\tilde{u} \right\rangle ,{{u}^{2}} \right\rangle ={{\int }_{\Omega }}\tilde{u}{{{\rm e}}^{u}}\nabla w\cdot \nabla \tilde{{{\lambda }^{2}}}\;{\rm d}\Omega ,\end{eqnarray} \tag{ 18 }$$

where $\tilde{{{\lambda }^{2}}}$ is the solution of (17), a variant of the second order adjoint equation (13). The Fisher information operator therefore coincides with the Gauss–Newton Hessian of the cost functional (3).

The procedure for computing the gradient acting on an arbitrary direction is clear. One first solves the first order forward equation (8) for $w$, then the first order adjoint (10) for λ, and finally evaluate (9). Similarly, one can compute the Hessian (or the Fisher information operators) acting on two arbitrary directions by first solving the second order forward equation (12) for ${{w}^{2}}$, then the second order adjoint equation (13) (or its variant (17)) for ${{\lambda }^{2}}$ (or $\tilde{{{\lambda }^{2}}}$), and finally evaluating (11) (or (18)).

One of the main goals of the paper is to study the RMHMC method in the context of Bayesian inverse problems governed by PDEs. It is therefore essential to compute the derivative of the Fisher information operator. This task is obvious for problems with available closed form expressions of the likelihood and the prior, but it is not so for those governed by PDEs. Nevertheless, using the adjoint technique we can compute the third order derivative tensor acting on three arbitrary directions with three extra PDE solves, as we now show. To that end, recall that the Fisher information operator acting on directions $\tilde{u}$ and ${{u}^{2}}$ is given by (18). The Fréchet derivative of the Fisher information operator along the additional direction ${{u}^{3}}$ (superscript '3' means the third variation direction) is given by

$$\begin{eqnarray}&&\begin{array}{l} \left\langle \left\langle \left\langle T\left( u \right),\tilde{u} \right\rangle ,{{u}^{2}} \right\rangle ,{{u}^{3}} \right\rangle :=\;\left\langle \nabla \left\langle \left\langle G\left( u \right),\tilde{u} \right\rangle ,{{u}^{2}} \right\rangle ,{{u}^{3}} \right\rangle \\ \quad ={{\int }_{\Omega }}\tilde{u}{{u}^{3}}{{{\rm e}}^{u}}\nabla w\cdot \nabla \tilde{{{\lambda }^{2}}}\;{\rm d}\Omega +{{\int }_{\Omega }}\tilde{u}{{{\rm e}}^{u}}\nabla {{w}^{3}}\cdot \nabla \tilde{{{\lambda }^{2}}}\;{\rm d}\Omega +{{\int }_{\Omega }}\tilde{u}{{{\rm e}}^{u}}\nabla w\cdot \nabla {{\lambda }^{2,3}}\;{\rm d}\Omega , \\ \end{array}\end{eqnarray} \tag{ 19 }$$

where ${{w}^{3}}$, ${{\lambda }^{2,3}}$ are the variation of $w$ and $\tilde{{{\lambda }^{2}}}$ in the direction ${{u}^{3}}$, respectively. One can show that ${{u}^{3}}$ satisfies another second order forward equation

$$\begin{eqnarray}&&{{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla {{w}^{3}}\cdot \nabla \hat{\lambda }\;{\rm d}\Omega +{{\int }_{\partial \Omega\backslash {{\Gamma }_{R}}}}Bi\;{{w}^{3}}\hat{\lambda }\;{\rm d}s=-{{\int }_{\Omega }}{{u}^{3}}{{{\rm e}}^{u}}\nabla w\cdot \nabla \hat{\lambda }\;{\rm d}\Omega .\end{eqnarray} \tag{ 20 }$$

Similarly, ${{\lambda }^{2,3}}$ is the solution of the third order adjoint equation

$$\begin{eqnarray}&&\begin{array}{r} {{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla {{\lambda }^{2,3}}\cdot \nabla \hat{w}\;{\rm d}\Omega +{{\int }_{\partial \Omega \backslash{{\Gamma }_{R}}}}Bi\;{{\lambda }^{2,3}}\hat{w}\;{\rm d}s=-\frac{1}{{{\sigma }^{2}}}\mathop{\sum }\limits_{j=1}^{K}{{w}^{2,3}}\left( {{{\bf x}}_{j}} \right)\hat{w}\left( {{{\bf x}}_{j}} \right)-{{\int }_{\Omega }}{{u}^{3}}{{{\rm e}}^{u}}\nabla \tilde{{{\lambda }^{2}}}\cdot \nabla \hat{w}\;{\rm d}\Omega , \\ \end{array}\end{eqnarray} \tag{ 21 }$$

and ${{w}^{2,3}}$, the variation of ${{u}^{2}}$ in direction ${{u}^{3}}$, satisfies the following third order forward equation

$$\begin{eqnarray}&&\begin{array}{l} {{\int }_{\Omega }}{{{\rm e}}^{u}}\nabla {{w}^{2,3}}\cdot \nabla \hat{\lambda }\;{\rm d}\Omega +{{\int }_{\partial \Omega\backslash {{\Gamma }_{R}}}}Bi\;{{w}^{2,3}}\hat{\lambda }\;{\rm d}s\;= \\ \quad -{{\int }_{\Omega }}{{u}^{3}}{{{\rm e}}^{u}}\nabla {{w}^{2}}\cdot \nabla \hat{\lambda }\;{\rm d}\Omega -{{\int }_{\Omega }}{{u}^{3}}{{u}^{2}}{{{\rm e}}^{u}}\nabla w\cdot \nabla \hat{\lambda }\;{\rm d}\Omega -{{\int }_{\Omega }}{{u}^{2}}{{{\rm e}}^{u}}\nabla {{w}^{3}}\cdot \nabla \hat{\lambda }\;{\rm d}\Omega . \\ \end{array}\end{eqnarray} \tag{ 22 }$$

Note that it would have required four extra PDE solves if one computes the third derivative of the full Hessian (11).

It is important to point out that the operator $T$ is only symmetric with respect to $\tilde{u}$ and ${{u}^{2}}$ since the Fisher information is symmetric, but not with respect to $\tilde{u}$ and ${{u}^{3}}$ or ${{u}^{2}}$ and ${{u}^{3}}$. The full symmetry only holds for the derivative of the full Hessian, that is, the true third derivative of the cost functional.

As presented in section 2, we view our inverse problem from an infinite dimensional point of view. As such, to implement our approach on computers, we need to discretize the prior, the likelihood and hence the posterior. We choose to use the FEM. In particular, we employ the standard ${{H}^{1}}\left( \Omega \right)$ FEM to discretize the forwards and adjoints (the likelihood), and the operator $\mathcal{A}$ (the prior). It should be pointed out that the Cameron–Martin space can be shown (see, e.g., [9]) to be a subspace of the usual fractional Sobolev space ${{H}^{s}}\left( \Omega \right)$, which is in turn a subspace of ${{H}^{1}}\left( \Omega \right)$. Thus, we are using a non-conforming FEM approach (outer approximation). For convenience, we further assume that the discretized state and parameter live on the same finite element mesh. Since FEM approximation of elliptic operators is standard (see, e.g., [13]), we will not discuss it here. Instead, we describe the MTT (see, e.g, [14] and the references therein) to discretize the prior.

Define $\mathcal{Q}:=\;{{\mathcal{C}}^{1/2}}={{\mathcal{A}}^{-s/2}}$, then the eigenpairs $\left( {{\lambda }_{i}},{{v}_{i}} \right)$ of $\mathcal{Q}$ define the KL expansion of the prior distribution as

$$\begin{eqnarray*}&&u={{u}_{0}}+\frac{1}{\sqrt{\alpha }}\mathop{\sum }\limits_{i=1}^{\infty }{{a}_{i}}{{\lambda }_{i}}{{v}_{i}},\end{eqnarray*}$$

where ${{a}_{i}}\sim \mathcal{N}\left( 0,1 \right)$. We need to solve

$$\begin{eqnarray*}&&\mathcal{Q}{{v}_{i}}={{\lambda }_{i}}{{v}_{i}},\end{eqnarray*}$$

or equivalently

$$\begin{eqnarray}&&{{\mathcal{A}}^{s/2}}{{v}_{i}}=\frac{1}{{{\lambda }_{i}}}{{v}_{i}}.\end{eqnarray} \tag{ 23 }$$

To solve (23) using the MTT, let us denote by ${\bf M}$ the mass matrix, and ${\bf K}$ the stiffness matrix resulting from the discretization of the Laplacian Δ. The representation of $\mathcal{A}$ in the finite element space (see, e.g., [15] and the references therein) is given by

$$\begin{eqnarray*}&&{\bf A}:=\;{{{\bf M}}^{-1}}{\bf K}+{\bf I}.\end{eqnarray*}$$

Let bold symbols denote the corresponding vector of FEM nodal values, e.g., ${\bf u}$ is the vector containing all FEM nodal values of $u$. If we define $\left( {{\sigma }_{i}},{{{\bf v}}_{i}} \right)$ as eigenpairs for ${\bf A},$ i.e.

$$\begin{eqnarray}&&{\bf A}{{{\bf v}}_{i}}={{\sigma }_{i}}{{{\bf v}}_{i}},\ {\rm or}\ {\bf AV}=\Sigma {\bf V},\end{eqnarray} \tag{ 24 }$$

where ${\bf v}_{i}^{T}{\bf M}{{{\bf v}}_{j}}={{\delta }_{ij}}$, and hence ${{{\bf V}}^{-1}}={{{\bf V}}^{T}}{\bf M}$, ${{\delta }_{ij}}$ is the Kronecker delta function, and $\Sigma$ is the diagonal matrix with entries σ_i. Since ${\bf A}$ is similar to ${{{\bf M}}^{-\frac{1}{2}}}\left( {\bf K}+{\bf M} \right){{{\bf M}}^{-\frac{1}{2}}}$, a symmetric positive definite matrix, ${\bf A}$ has positive eigenvalues. Using MTT method, the matrix representation of (23) reads

$$\begin{eqnarray*}&&{{{\bf A}}^{s/2}}{{{\bf v}}_{i}}=\frac{1}{{{\lambda }_{i}}}{{{\bf v}}_{i}},\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{{{\bf A}}^{s/2}}:=\;{\bf V}{{\Sigma }^{s/2}}{{{\bf V}}^{-1}}.\end{eqnarray*}$$

It follows that

$$\begin{eqnarray*}&&{{\lambda }_{i}}=\sigma _{i}^{-s/2}.\end{eqnarray*}$$

The Galerkin FEM approximation of the prior via truncated KL expansion reads

$$\begin{eqnarray}&&{\bf u}={{{\bf u}}_{0}}+\frac{1}{\sqrt{\alpha }}\mathop{\sum }\limits_{i=1}^{N}{{a}_{i}}{{\lambda }_{i}}{{{\bf v}}_{i}},\end{eqnarray} \tag{ 25 }$$

with ${\bf u}$ as the FEM nodal value of the approximate prior sample $u$ and N as the number of FEM nodal points. Note that for ease in writing, we have used the same notation $u$ for both infinite dimensional prior sample and its FEM approximation. Since $u\in {{L}^{2}}\left( \Omega \right)$, ${\bf u}$ naturally lives in $\mathbb{R}_{{\bf M}}^{N}$, the Euclidean space with weighted inner product ${{\left( \cdot ,\cdot \right)}_{{\bf M}}}:=\;\left( \cdot ,{\bf M}\cdot \right)$.

A question arises: What is the distribution of ${\bf u}$? Clearly ${\bf u}$ is a Gaussian with mean ${{{\bf u}}_{0}}$ since a_i are. The covariance matrix ${\bf C}$ for ${\bf u}$ is defined by

$$\begin{eqnarray*}&&{{\left( {\bf z},{\bf Cy} \right)}_{{\bf M}}}:=\;\mathbb{E}[\left( {\bf u}-{{{\bf u}}_{0}},{\bf Mz} \right)\left( {\bf u}-{{{\bf u}}_{0}},{\bf My} \right)]=\frac{1}{\alpha }{{{\bf z}}^{T}}{\bf MV}{{\Lambda }^{2}}{{{\bf V}}^{T}}{\bf My},\end{eqnarray*}$$

where we have used (25) to obtain the second equality and $\Lambda$ is the diagonal matrix with entries ${{\Lambda }_{ii}}=\lambda _{i}^{-1}$. It follows that

$$\begin{eqnarray}&&{\bf C}=\frac{1}{\alpha }{\bf V}{{\Lambda }^{2}}{{{\bf V}}^{T}}{\bf M},\end{eqnarray} \tag{ 26 }$$

as a map from $\mathbb{R}_{{\bf M}}^{N}$ to $\mathbb{R}_{{\bf M}}^{N}$, and its inverse can be shown to be

$$\begin{eqnarray*}&&{{{\bf C}}^{-1}}=\alpha {\bf V}{{\Lambda }^{-2}}{{{\bf V}}^{T}}{\bf M},\end{eqnarray*}$$

whence the distribution of ${\bf u}$ is

$$\begin{eqnarray}&&{\bf u}\sim \mathcal{N}\left( {{{\bf u}}_{0}},\alpha {\bf V}{{\Lambda }^{-2}}{{{\bf V}}^{T}}{\bf M} \right)\sim {\rm exp} [-\frac{\alpha }{2}{{\left( {\bf u}-{{{\bf u}}_{0}} \right)}^{T}}{\bf MV}{{\Lambda }^{-2}}{{{\bf V}}^{T}}{\bf M}\left( {\bf u}-{{{\bf u}}_{0}} \right)].\end{eqnarray} \tag{ 27 }$$

As a result, the FEM discretization of the prior can be written as

$$\begin{eqnarray*}&&\frac{\alpha }{2}\left\|u-{{u}_{0}}\right\|_{\mathcal{C}}^{2}:=\;\frac{\alpha }{2}{\left\|{{\mathcal{A}}^{s/2}}\left( u-{{u}_{0}} \right)\right\|}^{2}\mathop{\approx }\limits^{{\rm MTT}}\frac{\alpha }{2}{{\left( {\bf u}-{{{\bf u}}_{0}} \right)}^{T}}{\bf MV}{{\Lambda }^{-2}}{{{\bf V}}^{T}}{\bf M}\left( {\bf u}-{{{\bf u}}_{0}} \right).\end{eqnarray*}$$

Thus, the FEM approximation of the posterior is given by

$$\begin{eqnarray*}&&\pi \left( {\bf u}|{\bf d} \right)\sim {\rm exp} \left( -\frac{1}{2}\left| {\bf d}-\mathcal{G}\left( u \right) \right|_{{\bf L}}^{2} \right)\times {\rm exp} \left( -\frac{\alpha }{2}{{\left( {\bf u}-{{{\bf u}}_{0}} \right)}^{T}}{\bf MV}{{\Lambda }^{-2}}{{{\bf V}}^{T}}{\bf M}\left( {\bf u}-{{{\bf u}}_{0}} \right) \right).\end{eqnarray*}$$

The detailed derivation of the FEM approximation of infinite Bayesian inverse problems in general and the prior in particular will be presented elsewhere [16].

In this section we give a brief overview of the MCMC algorithms that we consider in this work. Some familiarity with the concepts of MCMC is required by the reader since an introduction to the subject is out of the scope of this paper.

5.1. Metropolis–Hastings

5.2. Riemann manifold Metropolis adjusted Langevin algorithm (RMMALA)

Denoting the log of the target density as $\mathcal{L}({\bf u})={\rm log} \pi ({\bf u})$, the manifold Metropolis Adjusted Langevin Algorithm method [19], defines a Langevin diffusion with stationary distribution $\pi ({\bf u})$ on the Riemann manifold of density functions with metric tensor ${\bf G}({\bf u})$. By employing a first order Euler integrator for discretising the stochastic differential equation a proposal mechanism with density $q({{{\bf u}}^{*}}|{{{\bf u}}^{t-1}})=\mathcal{N}({{{\bf u}}^{*}}|\mu ({{{\bf u}}^{t-1}},\epsilon ),{{\epsilon }^{2}}{{{\bf G}}^{-1}}({{{\bf u}}^{t-1}}))$ is defined, where is the integration step size, a parameter which needs to be tuned, and the kth component of the mean function $\mu {{({\bf u},\epsilon )}_{k}}$ is

$$\begin{eqnarray}&&\mu {{({\bf u},\epsilon )}_{k}}={{{\bf u}}_{k}}+\frac{{{\epsilon }^{2}}}{2}{{\left( {{{\bf G}}^{-1}}({\bf u}){{\nabla }_{{\bf u}}}\mathcal{L}({\bf u}) \right)}_{k}}-{{\epsilon }^{2}}\mathop{\sum }\limits_{i=1}^{N}\mathop{\sum }\limits_{j=1}^{N}{\bf G}({\bf u})_{i,j}^{-1}\Gamma _{i,j}^{k},\end{eqnarray} \tag{ 28 }$$

where $\Gamma _{i,j}^{k}$ are the Christoffel symbols of the metric in local coordinates. Note that we have used the Christoffel symbols to express the derivatives of the metric tensor, and they are computed using the adjoint method presented in section 3.

Due to the discretisation error introduced by the first order approximation convergence to the stationary distribution is not guaranteed anymore and thus the Metropolis–Hastings ratio is employed to correct for this bias. In [19] a number of examples are provided illustrating the potential of such a scheme for challenging inference problems.

One can interpret the proposal mechanism of RMMALA as a local Gaussian approximation to the target density where the effective covariance matrix in RMMALA is the inverse of the metric tensor evaluated at the current position. Furthermore a simplified version of the RMMALA algorithm, termed sRMMALA, can also be derived by assuming a manifold with constant curvature thus cancelling the last term in equation (28) which depends on the Christoffel symbols. Whilst this is a step forward in that much information about the target density is now embedded in the proposal mechanism it is still driven by a random walk. The next approach to be taken goes beyond the direct and scaled random walk by defining proposals which follow the geodesic flows on the manifold of densities and thus presents a potentially really powerful scheme to explore posterior distributions.

5.3. RMHMC

The RMHMC method defines a Hamiltonian on the Riemann manifold of probability density functions by introducing the auxiliary variables ${\bf p}\sim \mathcal{N}({\bf 0},{\bf G}({\bf u}))$ which are interpreted as the momentum at a particular position ${\bf u}$ and by considering the negative log of the target density as a potential function. More formally the Hamiltonian defined on the Riemann manifold is

$$\begin{eqnarray}&&\mathcal{H}({\bf u},{\bf p})=-\mathcal{L}({\bf u})+\frac{1}{2}{\rm log} \left( 2\pi |{\bf G}({\bf u})| \right)+\frac{1}{2}{{{\bf p}}^{T}}{\bf G}{{\left( {\bf u} \right)}^{-1}}{\bf p},\end{eqnarray} \tag{ 29 }$$

where the terms $-\mathcal{L}({\bf u})+\frac{1}{2}{\rm log} \left( 2\pi |{\bf G}({\bf u})| \right)$ and $\frac{1}{2}{{{\bf p}}^{T}}{\bf G}{{({\bf u})}^{-1}}{\bf p}$ are the potential energy and kinetic energy terms respectively. and the dynamics given by Hamiltons equations are

$$\begin{eqnarray}&&\frac{{\rm d}{{{\bf u}}_{k}}}{{\rm d}t}=\frac{\partial \mathcal{H}}{\partial {{{\bf p}}_{k}}}={{\left( {\bf G}{{\left( {\bf u} \right)}^{-1}}{\bf p} \right)}_{k}},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\begin{array}{l} \frac{{\rm d}{{{\bf p}}_{k}}}{{\rm d}t}=-\frac{\partial \mathcal{H}}{\partial {{{\bf u}}_{k}}}=\frac{\partial \mathcal{L}({\bf u})}{\partial {{{\bf u}}_{k}}}-\frac{1}{2}{\rm Tr}\left[ {\bf G}{{\left( {\bf u} \right)}^{-1}}\frac{\partial {\bf G}({\bf u})}{\partial {{{\bf u}}_{k}}} \right]+\frac{1}{2}{{{\bf p}}^{T}}{\bf G}{{\left( {\bf u} \right)}^{-1}}\frac{\partial {\bf G}({\bf u})}{\partial {{{\bf u}}_{k}}}{\bf G}{{\left( {\bf u} \right)}^{-1}}{\bf p}. \\ \end{array}\end{eqnarray} \tag{ 31 }$$

These dynamics define geodesic flows at a particular energy level and as such make proposals which follow deterministically the most efficient path across the manifold from the current density to the proposed one. Simulating the Hamiltonian requires a time-reversible and volume preserving numerical integrator. For this purpose the generalized leapfrog algorithm can be employed and provides a deterministic proposal mechanism for simulating from the conditional distribution, i.e. ${{{\bf u}}^{*}}|{\bf p}\sim \pi ({{{\bf u}}^{*}}|{\bf p})$. More details about the generalized lseapfrog integrator can be found in [19]. To simulate a path (which turns out to be a local geodesic) across the manifold, the Leapfrog integrator is iterated L times which along with the integration step size are parameters requiring tuning. Again due to the discrete integration errors on simulating the Hamiltonian in order to ensure convergence to the stationary distribution the Metropolis–Hastings acceptance ratio is applied.

The RMHMC method has been shown to be highly effective in sampling from posteriors induced by complex statistical models and offers the means to efficiently explore the hugely complex and high dimensional posteriors associated with PDE-constrained inverse problems.

As presented in section 5, we use the Fisher information matrix at the MAP point augmented with the Hessian of the prior as the metric tensor in our HMC simulations. It is therefore necessary to compute the augmented Fisher matrix and its inverse. In [1–3], we have shown that the Gauss–Newton Hessian of the cost functional (3), also known as the data misfit, is a compact operator, and that for smooth $u$ its eigenvalues decay exponentially to zero. Thus, the range space of the Gauss–Newton Hessian is effectively finite-dimensional even before discretization, i.e., it is independent of the mesh. In other words, the Fisher information matrix admits accurate low rank approximations and the accuracy can be improved as desired by simply increasing the rank of the approximation. We shall exploit this fact to compute the augmented Fisher information matrix and its inverse efficiently. We start with the augmented Fisher information matrix in $\mathbb{R}_{{\bf M}}^{N}$

$$\begin{eqnarray*}&&{\bf G}:=\;{{{\bf M}}^{-1}}{\bf H}+\alpha {\bf V}{{\Lambda }^{-2}}{{{\bf V}}^{T}}{\bf M}=\alpha {\bf V}{{\Lambda }^{-1}}\left( \frac{1}{\alpha }\Lambda {{{\bf V}}^{T}}{\bf HV}\Lambda +{\bf I} \right){{\Lambda }^{-1}}{{{\bf V}}^{-1}},\end{eqnarray*}$$

where ${\bf H}$ is the Fisher information matrix obtained from (18) by taking $\tilde{u}$ and ${{u}^{2}}$ as FEM basis functions.

Assume that ${\bf H}$ is compact (see, e.g., [1, 2]), together with the fact that ${{\Lambda }_{ii}}$ decays to zero, we conclude that the prior-preconditioned Fisher information matrix

$$\begin{eqnarray*}&&\tilde{{\bf H}}:=\;\frac{1}{\alpha }\Lambda {{{\bf V}}^{T}}{\bf HV}\Lambda ,\end{eqnarray*}$$

also has eigenvalues decaying to zero. Therefore it is expected that the eigenvalues of the prior-preconditioned matrix decays faster than those of the original matrix ${\bf H}$. Indeed, the numerical results in section 7 will confirm this observation. It follows that (see, e.g., [20, 21] for similar decomposition) $\tilde{{\bf H}}$ admits a r-rank approximation of the form

$$\begin{eqnarray*}&&\tilde{{\bf H}}=\frac{1}{\alpha }\Lambda {{{\bf V}}^{T}}{\bf HV}\Lambda \approx {{{\bf V}}_{r}}{\bf SV}_{r}^{T},\end{eqnarray*}$$

where ${{{\bf V}}_{r}}$ and ${\bf S}$ (diagonal matrix) contain the first r dominant eigenvectors and eigenvalues of $\tilde{{\bf H}}$, respectively. In this work, similar to [21], we use the one-pass randomized algorithm in [22] to compute the low rank approximation. Consequently, the augmented Fisher information matrix becomes

$$\begin{eqnarray*}&&{\bf G}\approx \alpha {\bf V}{{\Lambda }^{-1}}\left( {{{\bf V}}_{r}}{\bf SV}_{r}^{T}+{\bf I} \right){{\Lambda }^{-1}}{{{\bf V}}^{-1}},\end{eqnarray*}$$

from which we obtain the inverse, by using the Woodbury formula [23]

$$\begin{eqnarray*}&&{{{\bf G}}^{-1}}\approx \frac{1}{\alpha }{\bf V}\Lambda \left( {\bf I}-{{{\bf V}}_{r}}{\bf DV}_{r}^{T} \right)\Lambda {{{\bf V}}^{-1}},\end{eqnarray*}$$

where ${\bf D}$ is a diagonal matrix with ${{{\bf D}}_{ii}}={{{\bf S}}_{ii}}/\left( {{{\bf S}}_{ii}}+1 \right)$.

In the RMHMC method, we need to randomly draw the momentum variable as ${\bf p}\sim \mathcal{N}\left( {\bf 0},{\bf G} \right)$. If one considers

$$\begin{eqnarray*}&&{\bf p}=\sqrt{\alpha }{\bf V}{{\Lambda }^{-1}}{\bf b}+\sqrt{\alpha }{\bf V}{{\Lambda }^{-1}}{{{\bf V}}_{r}}{{{\bf S}}^{1/2}}{\bf c},\end{eqnarray*}$$

where ${{{\bf b}}_{i}},{{{\bf c}}_{i}}\sim \mathcal{N}\left( 0,1 \right)$, then one can show, by inspection, that ${\bf p}$ is distributed by $\mathcal{N}\left( {\bf 0},{\bf G} \right)$.

For convenience, let us recall that the FEM approximation of the posterior is given as

$$\begin{eqnarray}&&\pi \left( {\bf u}|{\bf d} \right)\sim {\rm exp} \left( -\frac{1}{2}\left| {\bf d}-\mathcal{G}\left( u \right) \right|_{{\bf L}}^{2} \right)\times {\rm exp} \left( -\frac{\alpha }{2}{{\left( {\bf u}-{{{\bf u}}_{0}} \right)}^{T}}{\bf MV}{{\Lambda }^{2}}{{{\bf V}}^{T}}{\bf M}\left( {\bf u}-{{{\bf u}}_{0}} \right) \right),\end{eqnarray} \tag{ 32 }$$

${{{\bf u}}_{0}}$ is the FEM nodal value of the prior mean function ${{u}_{0}}$, ${\bf M}$ is the mass matrix, ${\bf V}$ the matrix of eigenvectors defined in (24), $\Lambda$ the diagonal matrix introduced in (26), ${\bf L}={{\sigma }^{2}}{\bf I}$, ${\bf d}$ vector of observation data, and $\mathcal{G}\left( u \right)$ the forward map given by the forward equation

$$\begin{eqnarray*}&&\begin{array}{l} -\nabla \cdot \left( {{{\rm e}}^{u}}\nabla w \right)=0\qquad \quad \,{\rm in}\ \Omega , \\ -{{{\rm e}}^{u}}\nabla w\cdot {\bf n}=Bi\;u\qquad \ \ {\rm on}\ \partial \Omega\backslash {{\Gamma }_{R}}, \\ -{{{\rm e}}^{u}}\nabla w\cdot {\bf n}=-1\qquad \quad {\mkern 1mu} {\rm on}\ {{\Gamma }_{R}}, \\ \end{array}\end{eqnarray*}$$

that is discretized by the H¹-conforming FEM method.

In this section, we study Riemann manifold Monte Carlo methods and their variations to explore the posterior (32). In particular, we compare the performance of four methods: (i) sRMMALA obtained by ignoring the third derivative in RMMALA, (ii) RMMALA, (iii) sRMHMC obtained by first computing the augmented Fisher metric tensor at the MAP point and then using it as the constant metric tensor, (iv) RMHMC. For all methods, we form the augmented Fisher information matrix exactly using (18) with $\tilde{u}$ and ${{u}^{2}}$ as finite element basis vectors. For RMMALA and RMHMC we also need the derivative of the metric tensor which is a third order tensor. It can be constructed exactly using (19) with $\tilde{u},{{u}^{2}}$ and ${{u}^{3}}$ as finite element basis vectors. We also need extra work for the RMHMC method since each Stormer–Verlet step requires an implicit solve for both the first half of momentum and full position. For inverse problems such as those considered in this paper, the fixed point approach proposed in [19] does not seem to converge. We therefore have to resort to a full Newton method. Since we explicitly construct the metric tensor and its derivative, it is straightforward for us to develop the Newton scheme. For all problems considered in this section, we have observed that it takes at most five Newton iterations to converge.

Note that we limit ourselves in comparing these four methods in the Riemannian manifold MCMC sampling family. Clearly, other methods are available, we avoid 'unmatched comparison' in terms of cost and the level of exploiting the structure of the problem. Even in this limited family, RMHMC is most expensive since it requires not only third derivatives but also implicit solves, but the ability in generating almost independent samples is attractive and appealing as we shall show.

Though our proposed approach described in previous sections are valid for any spatial dimension d, we restrict ourselves to a one dimensional problem, i.e. d = 1, to clearly illustrate our points and findings. In particular, we take $\Omega =[0,1]$, ${{\Gamma }_{R}}=\left\{ 1 \right\}$. We set Bi = 0.1 for all examples. As discussed in section 2, for the Gaussian prior to be well-defined, we take $s=0.6\gt n/2=1/2$.

7.1. Two-parameter examples

We start our numerical experiments with two parameters. This will help demonstrate various aspects of RMHMC which are otherwise too computationally expensive for high dimensional problems. In particular, two-parameter example allows us to compute the complete third derivative tensor and perform the Newton method for each Stormer–Verlet step. This in turn allows us to show the capability of the full RMHMC over its simplified variants in tackling challenging posterior densities in which the metric tensor changes rapidly.

In order to construct the case with two parameters we consider FEM with one finite element. We assume that there is one observation point, i.e. K = 1, and it is placed at the left boundary x = 0. In the first example, we first take s = 0.6, $\sigma =0.1$, and $\alpha =0.1$. The posterior in this case is shown in figure 1(a). We start by taking a time step of 0.02 with 100 Stormer–Verlet steps for both sRMHMC and RMHMC. The acceptance rate for both methods is 1. One would take a time step of 2 for both sRMMALA and RMMALA to be comparable with sRMHMC and RMHMC, but the acceptance would be zero. Instead we take time step of 1 so that the acceptance rate is about 0.5 for sRMMALA and 0.3 for RMMALA. The MAP point is chosen as the initial state for all the chains with 5000 sample excluding the first 100 burn-ins. The result is shown in figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As can be seen, the RMHMC chain is the best in terms of mixing by comparing the second column (the trace plot) and the third column (the autocorrelation function (ACF)). Each RMHMC sample is almost uncorrelated to the previous ones. The sRMHMC is the second best, but the samples are strongly correlated compared to those of RMHMC, e.g. one uncorrelated sample for every 40. It is interesting to observe that the full RMMALA and sRMMALA have performance in terms of auto-correlation length that is qualitatively similar at least in the first 5000 samples. This is due to the RMMALA schemes being driven by a single step random walk that cannot exploit fully the curvature information available to the geodesic flows of RMHMC, see rejoinder of [19].

Note that it is not our goal to compare the behavior of the chains when they converge. Rather we would like to qualitatively study how fast the chains are well-mixed (mixing time). This is important for large-scale problems governed by PDEs since 'unpredicted' mixing time implies a lot of costly waste in PDE solves which one must avoid. Though RMHMC is expensive in generating a sample, the cost of generating an uncorrelated/independent sample seems to be comparable to sRMHMC for this example. In fact, if we measure the cost in terms of the number of PDE solves, the total number of PDE solves for RMHMC is 42 476 480 while it is 1 020 002 for sRMHMC, a factor of 40 more expensive. However, the cost in generating an almost uncorrelated/independent sample is the same since sRMHMC generates one uncorrelated sample out of 40 while it is one out of one for RMHMC.

To see how each method distributes the samples we plot one for every five samples in figure 3. All methods seem to explore the high probability density region very well. This explains why the sample mean and the 95% credibility region are similar for all methods in the first column of figure 2.

Zoom In Zoom Out Reset image size

In the second example we consider the combination s = 0.6, $\sigma =0.01$, and $\alpha =1$ which leads to the posterior shown in figure 1(b). For sRMHMC and RMHMC, we take time step $\varepsilon =0.04$ with L = 100 time steps, while it is 1 for both sRMMALA and RMMALA. Again, the acceptance rate is unity for both sRMHMC and RMHMC while it is 0.65 for sRMMALA and 0.55 for RMMALA, respectively. The result for four methods is shown in figure 4. As can be seen, this example seems to be easier than the first one since even though the time step is larger, the trace plot and the ACF looks better. It is interesting to observe that sRMHMC is comparable with RMHMC (in fact the ACF seems to be a bit better) for this example. As a result, RMHMC is more expensive than sRMHMC for less challenging posterior in figure 1(b). Here, by less challenging we mean that the posterior is quite well approximated by a Gaussian at the MAP point, e.g. the metric tensor is almost constant. This is true for the posterior in figure 1(b) in which the Gaussian prior contribution is significant, i.e., $\alpha =1$ instead of $\alpha =0.1$. Conversely, the posterior is challenging if the metric tensor changes rapidly. Similar to the first example, one also see that the sample mean and the 95% credibility region are almost the same for all methods.

Zoom In Zoom Out Reset image size

In the third example we consider the combination s = 0.6, $\sigma =0.01$, and $\alpha =0.1$ which leads to a skinny posterior with a long ridge as shown in figure 1(c). For sRMHMC and RMHMC, we take time step $\varepsilon =0.02$ with L = 100 time steps, while it is 1 for both sRMMALA and RMMALA. Again, the acceptance rate is unity for both sRMHMC and RMHMC while it is 0.45 for sRMMALA and RMMALA. The result for four methods is shown in figure 5. For this example, the RMHMC is more desirable than sRMHMC since the cost to generate an uncorrelated/independent sample is smaller for the former than the latter. The reason is that the total number of PDEs solves for the former is 40 times more than the latter, but one out of very sixty samples is uncorrelated/independent.

Zoom In Zoom Out Reset image size

7.2. Multi-parameter examples

In this section we choose to discretize $\Omega =[0,1]$ with ${{2}^{10}}=1024$ elements, and hence the number of parameters is 1025. For all simulations in this section, we choose s = 0.6, $\alpha =10$, and $\sigma =0.01$. For synthetic observations, we take K = 64 observations at ${{x}_{j}}=(j-1)/{{2}^{6}}$, $j=1,\ldots ,K$. Clearly, using the full blown RMHMC is out of the question since it is too expensive to construct the third derivative tensor and Newton method for each Stomer–Verlet step. For that reason, the sRMHMC becomes the viable choice. As studied in section 7.1, though sRMHMC loses the ability to efficiently sample from highly nonlinear posterior surfaces compared to the full RMHMC it is much less expensive to generate a sample since it does not require the derivative of the Fisher information matrix. In fact sRMHMC requires to (approximately) compute the Fisher information at the MAP point and then uses it as the fixed constant metric tensor throughout all leap-frog steps for all samples. Clearly, the gradient (9) has to be evaluated at each leap-frog step, but it can be computed efficiently using the adjoint method presented in section 3.

Nevertheless, constructing the exact Fisher information matrix requires 2 × 1025 PDEs solves. This is impractical if the dimension of the finite element space increases, e.g. by refining the mesh. Alternatively, due to the compactness of the Hessian of the prior-preconditioned misfit as discussed in section 6, we can use the randomized singular value decomposition (RSVD) technique [22] to compute its low rank approximations. Shown in figure 6 are the first 35 dominant eigenvalues of the Fisher information matrix and its prior-preconditioned counterpart. We also plot 20 approximate eigenvalues of the prior-preconditioned Fisher information matrix obtained from the RSVD method. As can be seen, the eigen spectrum of the prior-preconditioned Fisher information matrix decays faster than that of the original one. This is not surprising since the prior-preconditioned Fisher operator is a composition of the prior covariance, a compact operator, and the Fisher information operator, also a compact operator. The power of the RSVD is clearly demonstrated as the RSVD result for the first 20 eigenvalues is very accurate.

Zoom In Zoom Out Reset image size

Next, we perform the sRMHMC method using three different constant metric tensors: (i) the low rank Gauss–Newton Hessian, (ii) the exact Gauss–Newton Hessian, and (iii) the full Hessian. For each case, we start the Markov chain at the MAP point and compute 5100 samples, the first 100 of which is then discarded as burn-ins. The empirical mean (red line), the exact distributed parameter used to generate the observation (black line), and 95% credibility region are shown in figure 7. As can be seen, the results from the three methods are indistinguishable. The first sRMHMC is the most appealing since it requires $2\times 20=40$ PDE solves to construct the low rank Fisher information while the others need 2 × 1025 PDE solves. For large-scale problems with computationally expensive PDE solves, the first approach is the method of choice.

Zoom In Zoom Out Reset image size

To further compare the three methods we record the trace plot of the first two (1 and 2) and the last two (1024 and 1025) parameters in figure 8. As can be observed, the chains from the three methods seem to be well-mixed and it is hard to see the difference among them. We also plot the ACF for these four parameters in figure 9. Again, results for the three sRMHMC methods are almost identical, namely, they generate almost uncorrelated samples. We therefore conclude that low rank approach is the least computational extensive, yet it maintains the attractive features of the original RMHMC. As such, it is the most suitable method for large-scale Bayesian inverse problems with costly PDE solves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have proposed the adoption of a computationally inexpensive Riemann manifold Hamiltonian Monte method to explore the posterior of large-scale Bayesian inverse problems governed by PDEs in a highly efficient manner. We first adopt an infinite dimensional Bayesian framework to guarantee that the inverse formulation is well-defined. In particular, we postulate a Gaussian prior measure on the parameter space and assume regularity for the likelihood. This leads to a well-defined posterior distribution. Then, we discretize the posterior using the standard FEM and a MTT, and apply the RMHMC method on the resulting discretized posterior in finite dimensional parameter space. We present an adjoint technique to efficiently compute the gradient, the Hessian, and the third derivative of the potential function that are required in the RMHMC context. This is at the expense of solving a few extra PDEs: one for the gradient, two for a Hessian-vector product, and four for the product of third order derivative with a matrix. For large-scale problems, repeatedly computing the action of the Hessian and third order derivative is too computationally expensive and this motivates us to design a simplified RMHMC in which the Fisher information matrix is computed once at the MAP point. We further reduce the effort by constructing low rank approximation of the Fisher information using a RSVD technique. The effectiveness of the proposed approach is demonstrated on a number of numerical results up to 1025 parameters in which the computational gain is about two orders of magnitude while maintaining the quality of the original RMHMC method in generating (almost) uncorrelated/independent samples.

For more challenging inverse problems with significant change of metric tensor across the parameter space, we expect that sRMHMC with constant metric tensor is inefficient. In that case, RMHMC seems to be a better option, but it is too computational extensive for large-scale problems. Ongoing work is to explore approximation of the RMHMC methods in which we approximate the trace and the third derivative in (31) using adjoint and randomized techniques.

The first author is grateful for the support from DOE grants DE-SC0010518 and DE-SC0011118.