High-dimensional nonlinear Bayesian inference of poroelastic fields from pressure data

Abstract

We investigate solution methods for large-scale inverse problems governed by partial differential equations (PDEs) via Bayesian inference. The Bayesian framework provides a statistical setting to infer uncertain parameters from noisy measurements. To quantify posterior uncertainty, we adopt Markov Chain Monte Carlo (MCMC) approaches for generating samples. To increase the efficiency of these approaches in high-dimension, we make use of local information about gradient and Hessian of the target potential, also via Hamiltonian Monte Carlo (HMC). Our target application is inferring the field of soil permeability processing observations of pore pressure, using a nonlinear PDE poromechanics model for predicting pressure from permeability. We compare the performance of different sampling approaches in this and other settings. We also investigate the effect of dimensionality and non-gaussianity of distributions on the performance of different sampling methods.

Keywords

Hamiltonian Monte Carlo high-dimensional inference Markov Chain Monte Carlo poroelastic model

1. Introduction and background

Many problems in science and engineering problems can be modeled by partial differential equations (PDEs). Inverse problems constrained by PDEs are a challenging class of these problems, which play an essential role in investigating many physical systems, including geomechanical engineering, medical engineering, and astrophysics [1 –11]. Inverse problems aim at inferring the unknown model parameters of a physical system from observations, which may be limited, indirectly related to the parameters, and affected by noise. When the parameter domain is high-dimensional, and the relation between parameters and observations is defined by a complex mechanics model, solving inverse problems is computationally expensive. One approach for addressing these problems is to identify a point-based estimator via the minimizing of a function quantifying the discrepancies between model predictions and observations, e.g., the negative log-likelihood function. However, in this formulation the solution may not be unique, and the point-based estimate may be not informative on the posterior knowledge. One classical method to overcome the ill-posedness and to decrease the sensitivity of the solution to measurements’ noise is to add a regularization term, enforcing continuity and smoothness of the solution [1,12,13]. In the Bayesian framework, the unknown parameters are treated as random variables, with an assigned prior probability, while the likelihood function relies on the forward model and the noisy measurements [2,6]. Prior and likelihood are integrated in the posterior distribution. A point-based estimator for the posterior evaluation is the Maximum a Posteriori (MAP), but the posterior uncertainty can also be represented, describing the posterior distribution via samples or approximations, when it cannot be derived in exact form.

In this paper, we adopt a Bayesian inference framework to infer the uncertain parameters of a PDE-based forward poromechanics model. For solving the Bayesian inverse problem numerically, the infinite-dimensional parameter space has to be discretized to a finite-dimensional domain. Different discretization methods have been studied in the context of infinite-dimensional inverse problems. Finite element discretization methods use a finite number of continuous Lagrange basis functions to approximate the infinite-dimensional parameter space [3,6]. The Karhunen–Loeve (K-L) expansion provides an alternative approach, representing the infinite-dimensional parameter set in terms of its eigenvalues, and allowing of truncating after a finite number of terms and approximate dimensionality reduction [14,15].

We discuss the effectiveness of sampling methods based on Markov Chain Monte Carlo (MCMC) for representing the posterior distribution in the high-dimensional parameter space. The standard Metropolis-Hastings algorithm is computationally too expensive for our high-dimensional application [16]. There have been many efforts to accelerate the speed of sampling for high-dimensional inverse problems, such as developing reduced-order models [17 –20], using local gradient information as Hamiltonian Monte Carlo (HMC) method [21], Hessian information as Riemannian manifold HMC [22,23], and Hessian-based MCMC methods [2 –4,24]. In this paper, we investigate the efficiency of Hessian-based MCMC and Hessian-based HMC methods in high-dimensional problems.

Using Hessian information can significantly improve the performance of MCMC methods in infinite-dimensional inverse problems. Qi and Minka [24] proposed a Hessian-based Metropolis-Hasting (HMH) algorithm. They formulated an adaptive proposal density by approximating a Gaussian distribution using the local gradient and Hessian information. Martin et al. [2] applied the HMH method for high-dimensional inverse problems by approximating the low-rank local Hessian to ensure the positive definiteness of the matrix, obtaining what they call the Stochastic Newton MCMC (SNMCMC) method. Petra et al. [4] modified the SNMCMC method by using the Hessian at the MAP point to decrease its computational cost.

The HMC method can explore the posterior distribution faster than regular MCMC alternatives, as it uses local gradient information to make long-distance moves through the parameter space [21]. HMC algorithm can be used for efficient sampling in high-dimensional parameter spaces by selecting appropriate step size, total number of samples, and mass matrix [16]. Several works have been proposed to modify the HMC algorithm by tuning the mass matrix using Hessian information [16,22,25 –27]. Lee and Vempala [22] proved that the Riemannian Manifold HMC (RMHMC) accelerates the convergence rate of HMC by using the local Hessian information as the mass matrix, and Zhang and Sutton [27] developed a modified HMC method by using the approximated BFGS Hessian as the mass matrix. Bui-Thanh and Girolami [23] investigated a RMHMC method using the exact and low-rank Fisher information matrix.

We particularly study the performance of Hessian-informed MCMC and HMC methods that have been proposed recently for high-dimensional inverse problems. We compare the SN-MAP method presented by Petra et al. [4], the Metropolis- Adjusted Langevin algorithm (MALA) [28] using the Hessian information at the MAP point, and the H-HMC method presented by Bui-Thanh and Girolami [23]. We investigate several numerical examples to show how the dimension and non-Gaussian nature of the problem can affect the performance of different Hessian-informed methods. We avoid investigating Hessian-informed methods using the local Hessian information since those are computationally expensive and require calculating the second derivative information in each iteration.

There are many applications for high-dimensional inverse problems [2,4,5,23]. Our analysis targets the inference of soil permeability in regions close to injection sites related to a broad range of energy activities such as waste- water injection, CO2 sequestration, geothermal energy activities, and hydraulic fracking. Inferring these parameters is crucial for estimating the capacity of reservoirs, predicting the pore pressure and stress around the injection centers, and assessing the hazards of sliding and injection-induced seismic events [29 –39]. We infer the unknown permeability field from the collected pressure measurements. In our work, we adopt the Bayesian framework to infer the unknown poroelastic properties of a nonlinear poromechanics model from the noisy and sparse measurements. The MAP is identified using an inexact Newton solver and then samples are generating via MCMC approaches, starting from that point. We apply the developed model to solve a large-scale nonlinear inverse problem and determine the unknown properties of the deep underground layers.

1.1. Organization

Section 2 provides the general Bayesian formulation for high-dimensional inverse problems and discusses the choice of prior. Section 3 reviews and discusses the standard and accelerated sampling methods for exploring the posterior distribution. Section 4 provides the governing equations of the forward poromechanics model. Section 5 illustrates the numerical results of identifying the MAP point and generating samples from prior and posterior distributions.

2. Bayesian framework for high-dimensional inverse problems

We infer field $θ (x)$ , where $x$ is a spatial coordinate on reference domain $Ω$ , following the Bayesian framework. The prior knowledge about the field is modeled by a function $π (θ) = \log (p (θ))$ , and $p (θ)$ is assumed to be Gaussian, so that:

π (θ) = - \frac{1}{2} | | A (θ - m_{π}) | |^{2} + const .,

(1)

where $m_{π} (x)$ is the prior mean function and $A$ is a Laplacian-like operator, which are defined in Section 2.1. For any function $f \in L^{2} (Ω)$ , $L^{2} (Ω)$ being the space of square integrable functions, the norm is defined as:

| | f (x) | |^{2} = \int_{Ω} f^{2} (x) d x .

(2)

A dataset of $n$ noisy observations is available, listed in a vector $y = [y_{1}, y_{2}, \dots, y_{n}]^{⊤}$ . The observation $i$ is related to the field via equation $y_{i} = Q_{θ} (x_{i}, t_{i}) + ϵ_{i}$ , where the function $Q_{θ}$ predicts the value to be measured, as a function of the parameter field $θ$ , time $t$ , and location $x$ , and it is computed by solving the forward predictive model. Noise terms are listed in a vector $ε = [ϵ_{1}, ϵ_{2}, \dots, ϵ_{n}]^{⊤}$ , modeled as a Gaussian random variable, independent of $θ$ , with a zero mean and covariance matrix $Γ_{ϵ}$ , i.e., $ε ~ N (0, Γ_{ϵ})$ .

Probabilistically, the agreement between predictions and measurements is modeled by the log-likelihood function $l (θ) = \log (p (y | θ))$ :

l (θ) = - \frac{1}{2} | | Γ_{ϵ}^{- 1 / 2} (Q_{θ} - y) | |^{2} + const .,

(3)

where $Q_{θ} = [Q_{θ} (x_{1}, t_{1}), \dots, Q_{θ} (x_{n}, t_{n})]^{⊤}$ is the vector of predictions. In this paper, for each positive definite and symmetric matrix $G$ , the matrix $G^{- 1 / 2}$ is a decomposition of inverse of that matrix, so that $G^{- 1} = G^{- 1 / 2} G^{- ⊤ / 2}$ (e.g. the decomposition can be based on eigenvalue analysis or on Cholesky’s method). For column vectors, the norm is defined as $| | v | |^{2} = v^{⊤} v$ . If $Γ_{ϵ} = σ_{ϵ}^{2} I$ (where $I$ is the identity matrix), then $Γ_{ϵ}^{- 1 / 2}$ can be replaced by $σ_{ϵ}^{- 1}$ in equation (3).

Following the Bayes’ rule, the log-posterior density $ω (θ) = \log (p (θ | y))$ is as follows:

ω (θ) = l (θ) + π (θ) + const .,

(4)

Motivated by the above equation, we define the objective function as $J (θ) = - ω (θ) + const .$ :

J (θ) = \frac{1}{2} | | Γ_{ϵ}^{- 1 / 2} (Q_{θ} - y) | |^{2} + \frac{1}{2} | | A (θ - m_{π}) | |^{2} .

(5)

The MAP point is $θ_{MAP} = argmin (J (θ))$ , and it is a parameter field that maximizes the posterior density. To identify $θ_{MAP}$ , we minimize $J$ : while least squares approaches minimize the first term (based on the likelihood), the second one (based on prior knowledge) usually acts as a regularization.

2.1. Choice of prior model

The choice of $p (θ)$ is crucial as it allows to include the prior knowledge on the parameter field $θ$ in a Bayesian inference setting. In this paper, we use a continuous Gaussian measure $N (m_{π}, C_{π})$ instead of adopting a discretized Gaussian distribution and covariance matrix [3, 23]. Here, $C_{π}$ is the covariance operator which is defined as below:

(C_{π} ϕ) (x) = \int_{Ω} c (x, y) ϕ (y) d y,

(6)

where $c (x, z)$ is covariance function:

c (x, z) = E [(θ (x) - m_{π} (x)) (θ (z) - m_{π} (z))] .

(7)

Given measure function $μ_{π}$ , the operator is defined as:

C_{π} = \int_{Ω} (θ - m_{π}) \otimes (θ - m_{π}) μ_{π} (d θ) = E [(θ - m_{π}) \otimes (θ - m_{π})] .

(8)

The operator $C_{π}$ is fast-to-apply and allows a simple discretization in high dimension [3]. As it is described in Appendix 1, the prior function $π$ (equation (1)) is written as follows:

π (θ) = - \frac{1}{2} | | C_{π}^{- 1 / 2} (θ - m_{π}) | |^{2} + const .,

(9)

Commonly, the covariance operator is considered as the fractional power of Laplacian-like operator $A$ , i.e., $C_{π} = A^{- α}, α > 0$ , where $A$ is positive definite, self adjoint and invertible [40 –42]. Here, we choose the covariance operator as $C_{π} = A^{- 2}$ , so that $C_{π}^{- 1 / 2} = A$ , and the action of operator $A$ on $θ$ , i.e., $A θ$ is defined by the following elliptic Boundary Value Problem $(A = - γ Ψ Δ + δ I)$ :

\begin{matrix} - γ div (Ψ \nabla θ) + δ θ = s in Ω \\ γ (\nabla θ) \cdot n = 0 on \partial Ω, \end{matrix}

(10)

where $Δ$ is the Laplacian operator, $γ$ and $δ$ are positive (hyper)parameters, $s$ is a white noise field, and the random $θ$ field can be expressed as linear function of the random noise, as $A^{- 1} s$ . Here, $n$ is the unit normal vector on $\partial Ω$ . The symmetric and positive-definite second-order tensor $Ψ$ defines the anisotropy and indicates the direction $\hat{m}$ in which the parameter space $θ$ has a larger variation. $Ψ$ can then be defined in terms of the unit vector $\hat{m}$ as follows:

Ψ = a (I - (1 - b / a) \hat{m} \otimes {\hat{m}}^{⊤}) .

(11)

In a two-dimensional setting, $Ψ$ is can be represented as:

Ψ = [\begin{matrix} a \sin {(β)}^{2} + b \cos {(β)}^{2} & (b - a) \sin (β) \cos (β) \\ (b - a) \sin (β) \cos (β) & b \sin {(β)}^{2} + a \cos {(β)}^{2} \end{matrix}],

(12)

where $a$ , $b$ , and $β$ are constants. The larger $γ$ and $δ$ are, the smaller becomes the variance (equation (10)). Based on the above formulation, the continuous form of prior term can be written as:

A θ = {\begin{matrix} - γ div (Ψ \nabla θ) + δ θ & in Ω \\ γ (\nabla θ) \cdot n & on \partial Ω \end{matrix} .

(13)

2.2. Discretization of the field

We use a finite element discretization with continuous Lagrangian basis functions ${ϕ_{j}}_{j = 1}^{m}$ , and the corresponding nodal values ${x_{j}}_{j = 1}^{m}$ , to obtain the finite-dimensional approximation of the problem. Thus, the discretized field can be described as a vector of coefficients $ψ = [ψ_{1}, \dots, ψ_{m}]^{⊤}$ . The inner product between nodal coefficient vectors can be weighted by a positive-definite and symmetric matrix $N$ , with entries given by:

N_{ij} = \int_{Ω} ϕ_{i} (x) ϕ_{j} (x) d x .

(14)

Then, for every pair of functions $θ_{1}, θ_{2} \in L^{2} (Ω)$ , we can define the weighted inner product in the discretized form as $(θ_{1}, θ_{2})_{L^{2} (Ω)} \approx (ψ_{1}, ψ_{2})_{N} = ψ_{1}^{⊤} N ψ_{2}$ . The matrix representation of the Laplacian operator $A$ can be defined as $A = N^{- 1} K$ [34], with entry ${i, j}$ of $K$ as:

K_{ij} = \int_{Ω} [γ Ψ \nabla ϕ_{i} (x) \cdot \nabla ϕ_{j} (x) + δ ϕ_{i} (x) ϕ_{j} (x)] d x .

(15)

Therefore, the log prior term can be approximated as:

π (θ) ≅ - \frac{1}{2} | | A^{- 1 / 2} ψ | |^{2} + const .,

(16)

and the objective function can be expressed in terms of $ψ$ as:

J (ψ) = \frac{1}{2} | | Γ_{ϵ}^{- 1 / 2} (Q_{ψ - y}) | |^{2} + \frac{1}{2} | | A^{- 1 / 2} ψ | |^{2},

(17)

where $Q_{ψ}$ is the predictive vector, reconstructing field $θ = \sum_{j = 1}^{m} ψ_{j} ϕ_{j}$ . To minimize this nonlinear function with respect to $θ$ , quasi-Newton methods can be generally used.

3. Sampling methods for high-dimension inverse problems

In this section, we review some key methods for drawing samples from probability density functions. We discuss the Metropolis–Hastings MCMC and Hessian- based MCMC methods, and compare their computational costs and efficiency for sampling high-dimensional parameter space.

3.1. Metropolis Hastings method

The Metropolis Hastings MCMC (MH-MCMC) method (Algorithm 1) relies on proposal density $q$ , which samples next value of random variables $z$ , along a chain, centered at current location $c$ . One of the most popular choices of proposal (log) density is the isotropic Gaussian distribution:

q (c, z) = - \frac{1}{2} | | Δ t^{- 1} (c - z) | |^{2} + const .,

(18)

where $Δ t$ is a fixed step-size.

Algorithm 1 Metropolis-Hastings MCMC.
Input: Initial parameter $ψ_{0}$
1: Compute $J_{0} = J (ψ_{0})$ ,
2: for $k = 0, \dots, N$ do
3: Draw sample $z$ from proposal density $q (ψ_{k}, \cdot)$
4: Compute $J_{z} = J (z)$ , $Δ q = q (z, ψ_{k}) - q (ψ_{k}, z)$ ,
5: Compute $α = - J_{z} + J_{k} + Δ q$
6: Compute $a = min [1, \exp (α)]$
7: Draw $u$ from $U ([0, 1))$
8: if $u < a$ then
9: Accept: Set $ψ_{k + 1} = z$ , $J_{k + 1} = J_{z}$
10: else
11: Reject: Set $ψ_{k + 1} = ψ_{k}$ , $J_{k + 1} = J_{k}$
12: end if
13: end for
Output: Random walk: ${ψ_{k}}_{k = 1}^{N}$

In high-dimensional problems governed by PDEs, the posterior standard deviation of a linear combination of the field variables varies significantly depending of the direction of the combination, and this poses some challenges in choosing the step size. Selecting a large step size will induce a low acceptance rate, due to the small standard deviation of some directions, and a step size as small as the smallest standard deviation will induce a slow mixing (and, since we need to solve the forward model in each iteration, the simulation campaign will be computationally expensive).

In the following sections, we discuss some accelerated sampling methods using the local gradient and Hessian information of the objective function $J$ .

3.2. Hessian-based MCMC methods

In this section, we review three Hessian informed MCMC methods proposed by Ghattas and colleagues [2 –4]:

The local quadratic form of equation (5), based on the Newton method at the vicinity of a current sample $ψ_{k}$ can be approximated using the Laplace approximation (based on a second order Taylor series) as follows [43]:

J (ψ) \approx \tilde{J} (ψ) = \frac{1}{2} | | H_{k}^{1 / 2} (ψ - ψ_{k}) | |^{2} + g_{k}^{⊤} (ψ - ψ_{k}) + J (ψ_{k}),

(19)

where $g_{k} = \nabla J (ψ_{k})$ and $H_{k} = \nabla^{2} J (ψ_{k})$ are the gradient and the Hessian matrix of the objective function, respectively, and $\tilde{J}$ is the approximate function. We note that the Hessian matrix is positive definite at the MAP point. However, we cannot guarantee a positive-definite matrix at any arbitrary point. Therefore, to guarantee a positive-definite Hessian matrix, we can use a low-rank approximation ${\tilde{H}}_{k}$ (Appendix 2), which is calculated based on the eigenvalue decomposition, truncating negative and zero eigenvalues. By expanding equation (19), and adding and subtracting the term $\frac{1}{2} g_{k}^{T} {\tilde{H}}_{k}^{- 1} g_{k}$ , we obtain:

\tilde{J} (ψ) = \frac{1}{2} (ψ - ψ_{k})^{⊤} {\tilde{H}}_{k} (ψ - ψ_{k}) + g_{k}^{⊤} (ψ - ψ_{k}) + J (ψ_{k}) + \frac{1}{2} g_{k}^{T} {\tilde{H}}_{k}^{- 1} g_{k} - \frac{1}{2} g_{k}^{T} {\tilde{H}}_{k}^{- 1} g_{k} .

(20)

Therefore, equation (20) can be re-written as [2]:

\tilde{J} (ψ) = \frac{1}{2} | | {\tilde{H}}_{k}^{1 / 2} (ψ - ψ_{k} + {\tilde{H}}_{k}^{- 1} g_{k}) | |^{2} + const .,

(21)

where term ${\tilde{H}}_{k}^{- 1} g_{k}$ can be interpreted as the Newton step. Inspired by this approximation, the Stochastic Newton MCMC (SN-MCMC) method adopts the following proposal log density:

q (ψ_{k}, z) = - \frac{1}{2} | | {\tilde{H}}_{k}^{1 / 2} (z - ψ_{k} + {\tilde{H}}_{k}^{- 1} g_{k}) | |^{2} + const .

(22)

When implementing the method, we need to solve the forward problem at each iteration of the sampling process, and calculate the local Hessian, local gradient, and approximated Hessian at each iteration.

Since recalculation of Hessian and approximated Hessian at each sample point is expensive, Bui-Thanh et al. and Petra et al. [3,4] proposed a modified method which is based on always using the Hessian at the MAP point. Therefore, in the Stochastic Newton MCMC with MAP-based Hessian (SN-MAP) method, we rewrite the proposal log density function as follows:

q (ψ_{k}, z) = - \frac{1}{2} | | H_{MAP}^{1 / 2} (z - ψ_{k} + H_{MAP}^{- 1} g_{k}) | |^{2} + const .

(23)

Thus, implementing this approach at each iteration we need to run the forward problem, and compute the local gradient at each sampled point, but we calculate the Hessian matrix just once, at the MAP point.

Using the locally approximated proposal density can increase the convergence speed with respect to MH-MCMC. However, when the target distribution is far from a Gaussian distribution, changing the mean value of the proposal density by adding term $(H_{MAP}^{- 1} g_{k})$ , see equation (19), can decrease the speed of convergence or the sampling process can be stuck in a low probability region. In order to overcome these problems, we can define a small random learning rate $γ \leq 1$ , substituting term $(H_{MAP}^{- 1} g_{k})$ with term $(γ H_{MAP}^{- 1} g_{k})$ , or use a trust region method [24].

Another Hessian-informed MCMC method proposed by Bui et al. [3,4], called Independence Sampling with MAP point-based Gaussian proposal (IS-MAP), neglects recalculating the local gradient at each sample by adopting a proposal density centered at the MAP point:

q (ψ_{MAP}, z) = - \frac{1}{2} | | H_{MAP}^{1 / 2} (z - ψ_{MAP}) | |^{2} + const .,

(24)

where $ψ_{MAP}$ is the parameter value at the MAP point. Therefore, by avoiding recalculation of local Hessian and gradient at each sampling point, we need to calculate the Hessian at the MAP point only once and solve forward problem at each sampling point. Note that, in this method, samples are generated independently at each step, so the sample generation is simpler than that of any other MCMC method. However, the acceptance criteria is still the traditional one, for MH-MCMC. This method can be related to the classical rejection sampling approach.

In the Metropolis-Adjusted Langevin Algorithm (MALA) method [28,44,45], the proposal density is defined as:

q (ψ_{k}, z) = - \frac{1}{2} | | (2 τ B)^{- 1 / 2} (z - ψ_{k} + τ B g_{k}) | |^{2} + const .,

(25)

where $B$ is a positive-definite preconditioning matrix, and $τ$ is a fixed time step, which has to be taken as $0 < τ \leq 1$ . If we use the inverse of local Hessian information at the MAP as the preconditioning matrix, the proposal density becomes:

q (ψ_{k}, z) = - \frac{1}{2} | | (H_{MAP} / (2 τ))^{1 / 2} (z - ψ_{k} + τ H_{MAP}^{- 1} g_{k}) | |^{2} + const .

(26)

The scalar parameter $τ$ in equations (25) and (26) can play a similar role as the small learning rate $γ$ , by preventing the algorithm from being stuck in the low probability regions and yield better results for non-Gaussian distributions.

3.3. Hamiltonian Monte Carlo method

HMC is one of the efficient sampling methods relying on local gradient information. In the HMC algorithm (Algorithm 2), we define an auxiliary momentum value and then we update the position and momentum according to the Hamiltonian’s system of differential equations. For a general system, the position is the variable of interest (e.g. the parameters to be inferred). The potential energy is defined as the objective function of Section 2, as the negative logarithm of the target distribution [46, 47]. The Hamiltonian function can be defined as:

H (ψ, p) = K (ψ, p) + J (ψ),

(27)

where $K (ψ, p)$ is the kinetic energy and $p$ is the momentum vector. In order to update variables $p$ and $ψ$ , we need to calculate the derivatives of these variables with respect to the “pseudo time” variable $t$ . Therefore, Hamiltonian equations are as follows:

\frac{d ψ}{dt} = \frac{\partial H}{\partial p} = \frac{\partial K}{\partial p} + \frac{\partial J}{\partial p} = \frac{\partial K}{\partial p}

(28)

\frac{d p}{dt} = - \frac{\partial H}{\partial ψ} = - \frac{\partial K}{\partial ψ} - \frac{\partial J}{\partial ψ} = - \frac{\partial J}{\partial ψ} .

(29)

Algorithm 2 Hamiltonian Monte Carlo.
Input: Initial parameter $ψ_{0}$
1: for $k = 0, \dots, N$ do
2: Draw momentum sample $p_{k}$ from proposal density $N (0, M)$
3: Compute $H (ψ_{k}, p_{k})$ , $\frac{\partial J}{\partial ψ_{k}}$ , $\frac{\partial J}{\partial ψ_{k + 1}}$
4: Update $ψ$ (based on leapfrog algorithm)
5: $p_{k + \frac{1}{2}} = p_{k} - \frac{Δ t}{2} \frac{\partial J}{\partial ψ_{k}}$
6: $ψ_{k + 1} = ψ_{k} + Δ t M^{- 1} p_{k + \frac{1}{2}}$
7: $p_{k + 1} = p_{k + \frac{1}{2}} - \frac{Δ t}{2} \frac{\partial J}{\partial ψ_{k + 1}}$
8: Compute $H (ψ_{k + 1}, p_{k + 1})$
9: Compute $α = H (ψ_{k}, p_{k}) - H (ψ_{k + 1}, p_{k + 1})$
10: Compute $a = \min [1, \exp (α)]$
11: Draw $u$ from $U ([0, 1))$
12: if $u < a$ then
13: Accept: $ψ_{k + 1}$
14: else
15: Reject: Set $ψ_{k + 1} = ψ_{k}$
16: end if
17: end for

Different choices of $K$ function have been proposed in the literature. One of the most common choices is the Euclidean-Gaussian kinetic energy [21]:

K (ψ, p) = \frac{1}{2} | | M^{- 1 / 2} p | |^{2} + \frac{1}{2} \log | M | + const .,

(30)

where $M$ is a positive-definite matrix known as the “mass matrix.” The choice of the mass matrix influences the convergence speed and the performance of the sampling process.

3.3.1. Hessian-based HMC methods

It has been proved that we can reach a better convergence speed in linear problems by setting the mass matrix as the inverse of posterior covariance [25]. In a more general setting, Riemannian manifold HMC method uses a Riemannian-Gaussian kinetic energy [48] to choose the mass matrix as $M = C (ψ)$ , where $C (ψ)$ is a function of the parameter of interest which has to be updated with the local information [21].

When we consider $C (ψ) = \nabla^{2} J (ψ)$ , the Hamiltonian updating equations will induce an updating step similar to that of Newton method [22], which can significantly improve the sampling performance. Using the Riemannian manifold HMC for high-dimensional inverse problems is computationally expensive since we need to calculate the local Hessian of the posterior at each iteration.

To decrease the computational cost, we can use the Hessian of the posterior at the MAP point as a mass matrix [23], which prevent calculating Hessian in each iteration.

3.4. Comparison of SN-MCMC, MALA, and H-HMC sampling methods

This section compares the formulation of MCMC sampling methods relying on Hessian information for high-dimension parameter spaces. We discuss their similarities and differences.

3.4.1. SN-MAP

Following equation (23), the updating step of SN-MAP method is (here, we directly indicate the next candidate point as $ψ_{k + 1}$ , instead of as $z$ ):

ψ_{k + 1} = ψ_{k} - H_{MAP}^{- 1} g_{k} + H_{MAP}^{- 1 / 2} u,

(31)

where $u$ is a random sample from a multivariate standard Gaussian distribution (i.e., with identity covariance matrix and zero mean vector). For MCMC methods, as seen above, the acceptance criterion is based on log-ratio $α$ . For MH-MCMC (described in Algorithm 1), it is:

α = J (ψ_{k}) - J (ψ_{k + 1}) + Δ q .

(32)

The extended form of equation (32), using equation (23), becomes:

\begin{matrix} α = J (ψ_{k}) - J (ψ_{k + 1}) - \frac{1}{2} | | H_{MAP}^{1 / 2} (ψ_{k} - ψ_{k + 1} + H_{MAP}^{- 1} g_{k + 1}) | |^{2} + \frac{1}{2} | | H_{MAP}^{1 / 2} (ψ_{k + 1} - ψ_{k} + H_{MAP}^{- 1} g_{k}) | |^{2} \\ = J (ψ_{k}) - J (ψ_{k + 1}) + (ψ_{k + 1} - ψ_{k})^{⊤} (g_{k} + g_{k + 1}) + \frac{1}{2} | | H_{MAP}^{- 1 / 2} g_{k} | |^{2} - \frac{1}{2} | | H_{MAP}^{- 1 / 2} g_{k + 1} | |^{2} . \end{matrix}

3.4.2. MALA

Following equation (26), the updating step of MALA can be derived as:

ψ_{k + 1} = ψ_{k} - τ H_{MAP}^{- 1} g_{k} + \sqrt{(2 τ)} H_{MAP}^{- 1 / 2} u,

(33)

where $u$ is defined as above. Considering $\sqrt{(2 τ)} = Δ t$ , the above equation can be rewritten as:

ψ_{k + 1} = ψ_{k} - \frac{Δ t^{2}}{2} H_{MAP}^{- 1} g_{k} + Δ t H_{MAP}^{- 1 / 2} u .

(34)

The extended form of the log ratio for the acceptance criterion for MALA method, considering equation (26), can be derived as:

\begin{matrix} α = J (ψ_{k}) - J (ψ_{k + 1}) - \frac{1}{4 τ} | | H_{MAP}^{1 / 2} (ψ_{k} - ψ_{k + 1} + τ H_{MAP}^{- 1} g_{k + 1}) | |^{2} + \frac{1}{4 τ} | | H_{MAP}^{1 / 2} (ψ_{k + 1} - ψ_{k} + τ H_{MAP}^{- 1} g_{k}) | |^{2} \\ = J (ψ_{k}) - J (ψ_{k + 1}) + \frac{1}{2} (ψ_{k + 1} - ψ_{k})^{⊤} (g_{k} + g_{k + 1}) + \frac{τ}{4} | | H_{MAP}^{- 1 / 2} g_{k} | |^{2} - \frac{τ}{4} | | H_{MAP}^{- 1 / 2} g_{k + 1} | |^{2} . \end{matrix}

(35)

3.4.3. H-HMC

Similarly, the updating step of H-HMC method using the Hessian information at the MAP point and leapfrog algorithm 2 can be written as:

\begin{matrix} p_{k + \frac{1}{2}} = p_{k} - \frac{Δ t}{2} g_{k} \\ ψ_{k + 1} = ψ_{k} + Δ t H_{MAP}^{- 1} p_{k + \frac{1}{2}} \\ p_{k + 1} = p_{k + \frac{1}{2}} - \frac{Δ t}{2} g_{k + 1}, \end{matrix}

(36)

where $p_{k}$ is a random sample from a Gaussian distribution with zero mean and covariance $H_{MAP}$ . Therefore, the updating step is as equation (34). The extended form of the log ratio, for H-HMC method, considering equation (27), is as follows:

α = J (ψ_{k}) + \frac{1}{2} | | H_{MAP}^{- 1 / 2} p_{k} | |^{2} - J (ψ_{k + 1}) - \frac{1}{2} | | H_{MAP}^{- 1 / 2} p_{k + 1} | |^{2},

(37)

where from the leapfrog algorithm (equation (36)) we can define:

p_{k} = \frac{{(ψ_{k + 1} - ψ_{k})}^{⊤}}{Δ t} H_{MAP} + \frac{Δ t}{2} g_{k}

(38)

p_{k + 1} = \frac{{(ψ_{k + 1} - ψ_{k})}^{T}}{Δ t} H_{MAP} - \frac{Δ t}{2} g_{k + 1}

and by substituting equation (38) in equation (37) we derive:

α = J (ψ_{k}) - J (ψ_{k + 1}) + \frac{1}{2} (ψ_{k + 1} - ψ_{k})^{⊤} (g_{k} + g_{k + 1}) + \frac{Δ t^{2}}{8} | | H_{MAP}^{- 1 / 2} g_{k} | |^{2} - \frac{Δ t^{2}}{8} | | H_{MAP}^{- 1 / 2} g_{k + 1} | |^{2},

which is similar to the acceptance criterion of MALA method (equation (35)).

We summarize the important points from the comparison of Hessian-based sampling methods. The above analysis proves that, with a proper choice of positive-definite matrix $B$ in the MALA method and Mass matrix $M$ in the HMC method, the updating step and acceptance coefficient of these two methods are the same. Choice of Mass matrix $M$ in HMC method and of positive-definite matrix $B$ in MALA method plays a crucial role in tuning and performance of HMC and MALA methods, respectively. In the above analysis, we propose using the Hessian information at the MAP point to tune HMC and MALA methods, which can significantly improve the performance of these methods compared to standard sampling methods. Comparison of these three Hessian-based methods shows that, although the SN-MAP method performs perfectly for high-dimensional Gaussian distributions, H-HMC and MALA methods with a variable step size can perform better for exploring non-Gaussian parameter spaces. These methods will be compared numerically in Section 5.

4. Poroelastic forward model

In this section, we present the forward model, based on a continuum energetic formulation of poroelasticity [49], to predict the water pore pressure $p$ in the soil as a function of the permeability $κ$ and accounting for fluid transport in a deformable medium. Relating this to the notation of Section 2, the observation $i$ is a noisy pressure measure $y_{i} = p_{obs}^{i}$ , and the uncertain field is $θ (x) = - \underset{e}{\log} κ (x)$ , where $κ$ has the units $m^{2}$ . Below, we describe the governing PDEs of this model.

The governing equations of the porous medium, wherein the solid matrix is incompressible, can be written as:

\frac{dtr (ε)}{d t} + \frac{d p}{d t} \frac{ϕ_{f}}{K_{f}} = - \frac{1}{ρ_{f}} div Q

(39)

div σ + b = 0,

(40)

where $ϵ$ , $ϕ_{f}$ , $K_{f}$ , and $γ$ are the strain tensor, fluid volume fraction, fluid bulk modulus, and dynamic viscosity of fluid, respectively.

The flux vector $Q$ is defined as:

Q = - \frac{κ}{γ} ρ_{f} (\nabla p + P_{f} g),

(41)

where $ρ_{f}$ and $P_{f}$ are the true fluid density and fluid density, respectively, with the relation $P_{f} = ϕ_{f} ρ_{f}$ . In equation (40) $σ$ is the Cauchy stress tensor, which is defined as:

σ = (1 - ϕ_{f}) (2 μ ε + λ tr (ε) I) - ϕ_{f} p I,

(42)

where $μ$ and $λ$ are Lame constants, and the body force is defined as $b : = (P_{s} + P_{f}) g$ . where $P_{s}$ is the density of solid phase.

5. Numerical results and discussion

In this section, we illustrate the results of applying the Bayesian formulation presented in Sections 2-3, and we compare the performance of Hessian-based sampling methods, to different examples related to the processing of pore pressure data to infer permeability. We compare the sampling methods on a Gaussian target distribution in Section 5.1, on a non-Gaussian target distribution in Section 5.2, on the posterior distribution defined by a likelihood function embedding the coupled poroelastic model of Section 4 in Section 5.3. This latter analysis is implemented in FEniCS, and the quasi-Newton solver from the dolfin-adjoint package is used to solve the nonlinear minimization problem.

5.1. Gaussian posterior distribution

In this example, we consider an analytical Gaussian posterior distribution with the covariance matrix $Σ$ and mean $m$ . We applied four sampling methods: (1) MH-MCMC, (2) HMC, (3) SN-MCMC, and (4) H-HMC. The acceptance rate for SN-MCMC method is one as, considering equations (22) and (32), it can be shown that, when the distribution is Gaussian, $q (ψ_{k}, z)$ is identical to $ω (z)$ . Consequently, the log ratio is as follows:

α = ω (z) - ω (ψ_{k}) + q (z, ψ_{k}) - q (ψ_{k}, z) = 0

(43)

and the acceptance coefficient becomes:

a = min {1, \exp (α)} = 1 .

(44)

Figure 1 compares the autocorrelation vs. lag of these four sampling methods, for a one-dimensional and for a high-dimensional Gaussian distribution. In the one-dimensional distribution, we assume $m = 0.5$ and $σ^{2} = 2.0$ . We generate $2000$ samples, all the sampling processes start at $θ = 0.5$ , and the step size for MH-MCMC, HMC, and H-HMC methods is $Δ t = 1$ .

For the high-dimensional problem, we assumed a domain $Ω$ discretized into $936$ meshpoints. The mean value and the covariance matrix derived from the problem solved in Section 5.3 (the covariance matrix is as the inverse of Hessian at the MAP point). We generate $10, 000$ samples, and all the sampling processes start at the MAP (i.e. $ψ = m$ ). The step size for MH-MCMC method is $Δ t = 0.01$ , for HMC methods is $Δ t = 0.1$ , and for H-HMC method it is assumed $Δ t = 0.3$ , and considering these step sizes the acceptance rates for MH-MCMC, HMC, and H-HMC methods are approximately $85 %$ , $94 %$ , and $92 %$ , respectively. Moreover, the acceptance rate for SN-MAP method is unitary (as shown in equation (44)).

We note that, in the standard HMC and MH-MCMC sampling methods, we need to take a small step size to keep the acceptance rate high and prevent sampling from trapping in the low probable regions, which is inefficient and computationally expensive. However, using accelerated sampling methods such as H-HMC, we can choose a bigger step size while keeping the acceptance rate high, which can significantly decrease the number of samples we need to explore the parameter space.

As it can be seen in Figure 1, for Hessian-based sampling methods (SN-MCMC and H-HMC) the autocorrelation in both one-dimensional and high-dimensional settings quickly goes to zero, which demonstrates a faster convergence in comparison to MH-MCMC and HMC methods. *Also, comparison of H-HMC and HMC methods in high-dimension shows that we can use a bigger step size for the H-HMC method without decreasing the acceptance rate, and this increases the convergence speed. Moreover, comparison of the autocorrelation functions for the one-dimensional and high-dimensional cases shows that the performance of MH-MCMC method decreases significantly with the dimension of the problem. For the high-dimensional setting, Figure 2 shows the $95 %$ credible interval and the sample average along the section indicated by the dashed line in the domain, based on $10, 000$ samples. The results show that the credible interval and sample average for SN-MAP and H-HMC methods match well with the analytical solution, which are indistinguishable from the analytical solutions. However, standard sampling methods as MH-MCMC and HMC need more samples to reach the same level of accuracy (Figure 2).

Table 1 represents a summary of the analysis, via correlation time $τ$ defined as $τ = 1 + 2 \sum_{t = 1}^{\infty} ρ_{t}$ . where $ρ_{t}$ is the correlation between two chains with lag equal to $t$ . $N_{eff}$ is the effective sample size which is $(N / τ)$ , and SE is the standard error which is defined as $\tilde{σ} / N_{eff}$ where $\tilde{σ}$ is the sample chain standard deviation.

Table 1.

Summary analysis of high-dimensional normal example: summary of convergence analysis of different sampling methods for high-dimensional normal distribution.

Method	$τ$	$N_{eff}$	SE
HMC	1149.3	9	0.0036
H-HMC	62.5	160	0.00027
SN-MAP	1	10,000	1.4 × $10^{- 5}$

SE: standard error; HMC: Hamiltonian Monte Carlo; H-HMC: Hessian-based Hamiltonian Monte Carlo; SN-MAP: Stochastic Newton MCMC with MAP-based Hessian.

The table shows how the SN-MAP sampling method has the lowest standard error and it is the most efficient method for the high-dimensional normal distribution.

We do not consider the MH-MCMC method in this analysis, since the sampling chains of this method are yet far from convergence after $10, 000$ iterations.

Table 2 represents five different cases, with different levels of uncertainty, to investigate the effect of measurement number and noise level on the credible interval amplitude. We increase the noise level by decreasing the number of measurements and increasing the noise level. For each case, the MAP point and Hessian at the MAP point are computed. Then, considering the normal distribution approximating the posterior distribution at the MAP, we take $10, 000$ samples using the H-HMC method and compute the $95 %$ credible interval at four random points in the domain, for each case. The results are reported in Figure 3. As we expected, the results show that by increasing the uncertainty level, the amplitude of the credible interval increases (Figure 3).

Table 2.

Summary analysis of high-dimensional normal example: cases with different uncertainty levels, level of noise ( $σ_{ϵ}$ ), and number of measurements.

Case #	$σ_{ϵ}$	Meas.
1	0.5 MPa	60
2	2 MPa	60
3	4 MPa	60
4	4 MPa	20
5	10 MPa	20

5.2. Nearly Gaussian posterior distribution

In this section, we apply the sampling methods to a log-normal target distribution. The objective function is defined as:

J (ψ) = \frac{1}{2} | | Λ^{1 / 2} (\log (ψ) - m_{l}) | |^{2} - \log (Π_{i = 1}^{N} ψ_{i}^{- 1}),

(45)

where $Λ = Σ_{l}^{- 1}$ . $Σ_{l}$ and $m_{l}$ are covariance matrix and mean vector inputs for the log-normal distributions, respectively. The local gradient and Hessian of the objective function can be derived as follows:

\frac{\partial J (ψ)}{\partial ψ_{j}} = ψ_{j}^{- 1} [n_{j}^{T} Λ (\log (ψ) - m_{l}) + 1],

(46)

where $n_{j} = {n_{j}^{i}}_{i = 1}^{N}$ is a vector with components $n_{j}^{i} = 0, i \neq j$ and $n_{j}^{i} = 1, i = j$ :

\frac{\partial^{2} J (ψ)}{\partial ψ_{j} \partial ψ_{i}} = {\begin{matrix} - ψ_{j}^{- 2} [n_{j}^{T} Λ (\log (ψ) - m_{l}) + 1 + Λ_{j j}] & if i = j \\ ψ_{j}^{- 1} [Λ_{j i} ψ_{i}^{- 1}] & if i \neq j \end{matrix}

(47)

Based on equation (46), the MAP point is $ψ_{MAP} = \exp (m_{l} - Σ_{l} 1)$ . We discuss the performance of standard and accelerated sampling methods for one-dimensional and high-dimensional log-normal distributions. In the one-dimensional case, we assume $m_{l} = 0.5$ and $σ_{l}^{2} = 2.0$ . We generate $4000$ samples, and the step size for MH-MCMC, HMC, and H-HMC methods is $Δ t = 1$ . All the sampling processes start from the mode of log-normal distribution, $\exp (m_{l} - σ_{l}^{2})$ .

As a second setting, we consider a high-dimensional log-normal distribution. The domain and distribution are assumed as Section 5.1, and the MAP point $ψ_{MAP} = \exp (m_{l} - Σ_{l} 1)$ is assumed similar to the MAP point of previous example (Figure 2). Similarly, the step size for MH-MCMC method is $Δ t = 0.01$ , for HMC method is $Δ t = 0.1$ , and for H-HMC method it is $Δ t = 0.3$ . The correspondence acceptance rates after taking $10, 000$ samples for MH-MCMC, HMC, and H-HMC methods are $82 %$ , $94 %$ , and $92 %$ , respectively. In addition, the acceptance rate for SN-MAP method is approximately $8.3 %$ , which is significantly lower than that of the SN-MAP sampling method for high-dimensional Gaussian distribution.

Figure 4 represents the autocorrelation function vs. lag, for one-dimensional and high-dimensional settings. The results show that the performance of the four methods in one-dimensional distribution is quite similar. In addition, as it is shown in Figure 4 the convergence speed of SN-MAP method and H-HMC method with the Hessian information at the MAP point in high-dimensional distribution are significantly faster than standard MH-MCMC method. We also considered another example, where we take sample from a log-normal distribution and the MAP point is assumed $ψ_{MAP} = \exp (m)$ . where $m$ is assumed similar to the MAP point of the normal example. We have reported the results in Appendix 5.

Figure 5 presents the $95 %$ credibility interval and sample average along the dashed line for MH-MCMC, HMC, SN-MAP, and H-HMC methods. As the results show, the HMC method converges and explores the parameter space faster than the MH-MCMC method. H-HMC and SN-MAP methods convergence speeds are significantly higher than standard sampling methods. After $10, 000$ samples, the sample credible interval matches well with the exact credible interval for both H-HMC and SN-MAP methods.

Table 3 shows the summary of convergence analysis for the high-dimensional log-normal distribution, after taking $10, 000$ samples using MH-MCMC, HMC, and H-HMC sampling methods. Results of Table 3 show a decrease in the value of $N_{eff}$ compared to applying similar sampling methods to high-dimensional normal distribution (Table 1), which proves that taking samples from the log-normal distribution is more difficult.

Table 3.

Summary analysis of high-dimensional log-normal example: summary of convergence analysis of different sampling methods for high-dimensional log-normal distribution.

Method	$τ$	$N_{eff}$	SE
HMC	1458.1	7	0.0062
H-HMC	70.7	141	0.0003
SN-MAP	87.7	114	0.0004

SE: standard error; HMC: Hamiltonian Monte Carlo; H-HMC: Hessian-based Hamiltonian Monte Carlo; SN-MAP: Stochastic Newton MCMC with MAP-based Hessian.

Moreover, comparison of SE values in Tables 1 and 3 shows that although the SN-MAP method works well for normal distributions and it has the lowest estimation error value, for log-normal distribution using the H-HMC method is more efficient, and the estimation error is lower.

Table 4 presents the Kullback–Leibler divergence (KLD) values between normal and log-normal distributions and the acceptance rates for different log-normal distributions. This analysis considers different log-normal distributions and calculates the KLD value, which measures the difference between each log-normal distribution and its corresponding normal approximation at the MAP. For each log-normal distribution, we run the SN-MAP sampling method and take $10, 000$ samples. As the results show, by increasing the KLD value, the acceptance rate decreases, which shows that, for distributions that are not close to Gaussian, generating samples is more expensive as the acceptance rate is lower. The details related to computing the KLD between normal and log-normal distributions is presented in Appendix 4.

Table 4.

Summary analysis of high-dimensional log-normal example: results of SN-MAP sampling method for different log-normal distributions.

Case #	KLD val.	acceptance $%$
1	2.47	8.9
2	3.15	8.3
3	5.04	7.3
4	13.55	4.3
5	69.03	1.7
6	195.02	0.3

KLD: Kullback-Leibler divergence; SN-MAP: Stochastic Newton MCMC with MAP-based Hessian.

5.3. Non-Gaussian posterior distribution: inferring permeability from pressure data

In this section, we infer the permeability field $κ$ from the sparse pore pressure data. The PDE governing equations of forward model are described in section 4. The properties of assumed underground layer are listed in Table 5. The boundary conditions are assumed as below:

top BC : p = const ., bottom BC : δ = 0 .

Table 5.

Properties of solid and fluid phases in the poroelastic model, see Section 4.

Property	Value
Solid phase Lame constant, $λ$	40 MPa
Solid phase Lame constant, $μ$	40 MPa
Fluid bulk modulus, $K_{f}$	2270 MPa
Water density, $ρ_{f}$	1000 kg/m³
Water viscosity, $γ$	0.001 Pa · s
Solid density, $ρ_{s}$	2000 kg/m³

The initial porosity is assumed $ϕ_{0 f} = 0.2$ . We consider a $8000 \times 4000$ deep underground layer with the constant fluid pressure of $p = 500 M$ on the top boundary. Figure 6 shows the assumed (hidden) permeability field $κ$ , used to generate artificial measures by running the forward model. In this example, we generated $52$ synthetic observations of pore pressure, affected by random noise with standard deviation $σ_{ϵ} = 1.0 M$ . The observations are assumed in a single time horizon 10 days after start of the injection.

As discussed in Section 2, the prior covariance is defined based on the Laplacian operator $A^{- 1} = - γ Ψ Δ + δ I$ , with $γ = 0.5$ , $δ = 5 \times 10^{- 3}$ (1/km²), $a = 0.018$ , $b = 0.97$ , and $β = 1.017 π$ (the prior parameters are defined in Section 2.1).

5.3.1. Computing the MAP point

To find the MAP point, we apply a quasi Newton algorithm (BFGS) to solve the nonlinear minimization problem. The detailed formulation of Lagrangian method is presented in Appendix 3. The prior mean is assumed constant $(m_{π} = 33 m^{2})$ . We start with an initial guess for the field and by applying the Newton method, we iteratively update the field value. Figure 7 shows the MAP point (which is similar to actual permeability shown in Figure 6).

Figure 8 shows the identified MAP points by considering different noise levels. By increasing the level of noise $σ_{ϵ}$ , the prior term becomes more dominant and consequently the MAP field is smoother.

Similarly, Figure 9 shows the effect of number of measurements on computing the MAP point. By decreasing the number of observations, the MAP point becomes less informed, more smooth and far from the target distribution.

5.3.2. Exploring the posterior distribution

We use the HMC and H-HMC methods to generate samples from the posterior distribution. In the H-HMC method, we use the Hessian information at the MAP point to improve the sampling performance. To achieve similar acceptance rates, the step size $Δ t$ for HMC and H-HMC methods are set as $0.1$ and $0.3$ , respectively. After taking $20, 000$ samples for HMC, the correspondence acceptance rate and H-HMC methods are approximately $85 %$ and $83 %$ , respectively. Figure 10 shows the autocorrelation function vs. lag, computed for the permeability value at three points in the domain using HMC and H-HMC methods.

In this analysis, we don’t use the MH-MCMC method due to the inefficiency and computational costs. Also, since the posterior distribution is highly nonlinear and non-Gaussian, the acceptance rate for the SN-MAP method is very low, which makes this method inefficient.

Furthermore, Figure 11 presents the credible interval of the inverse problem’s solution along the vertical dashed line, after taking $50, 000$ samples using H-HMC method. The graph also shows the prior samples’ mean and $95 %$ credible interval. The comparison of prior and posterior intervals shows that the posterior variance has decreased significantly, ensuring that the posterior distribution is highly influenced by observations. In addition, posterior samples’ mean, MAP point, and the target value are shown in Figure 11. As the results show, the posterior samples’ mean and MAP point match well.

6. Conclusion

We have presented a Bayesian inference framework for high-dimension inverse problems governed by PDE equations. We used a continuous Gaussian prior distribution (based on a Laplacian-like operator) to ensure the well-posedness of the infinite-dimensional inverse problem. The main advantage of this prior function is that it can be applied to continuous domains, and it allows simple discretization. We implemented the inverse problem in the FEniCS library and used a quasi-Newton solver (BFGS) of dolfin-adjoint package to solve the minimization problem.

We have investigated several sampling methods to describe the posterior distribution. We discussed the complexities of sampling methods in high-dimensional parameter spaces and compared the performance of MH-MCMC and HMC sampling methods with the accelerated methods using the Hessian and gradient information, such as SN-MCMC, MALA, and H-HMC. By considering several one-dimensional and high-dimensional problems with Gaussian and non-Gaussian posterior distributions, we showed that using the modified Hessian-based sampling methods can significantly increase the speed of convergence and exploring in high-dimensional inverse problems. However, using the local Hessian information in nonlinear high-dimensional problems is computationally expensive. We considered Hessian sampling methods using the Hessian information calculated at the MAP point to overcome this problem. The results revealed that the MAP-based Hessian sampling methods are both computationally efficient and fast in exploring high-dimensional distributions.

We also applied the developed framework to a high-dimensional inverse problem governed by poroelastic PDE equations to infer the unknown permeability distribution from the point-wise pore pressure observations. We calculated the MAP point and the posterior credible intervals by applying the HMC and H-HMC methods, using the Hessian information at the MAP point. Our results indicate that the H-HMC method using the Hessian information at the MAP point has a better performance in exploring this non-Gaussian high-dimensional distribution.

Footnotes

Appendix 1

Appendix 2

Appendix 3

Appendix 4

Appendix 5 Software availability

A version of the code developed for this work is available at .

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: We thank the National Science Foundation for support through XSEDE resources provided by Pittsburgh Supercomputing Center. M.K. acknowledges financial support from the Scott Institute. K.D. acknowledges financial support from NSF (CMMI MOMS 1635407, DMS 2108784), ARO (MURI W911NF-19-1-0245), ONR (N00014-18-1-2528), BSF (2018183), and an appointment to the National Energy Technology Laboratory sponsored by the U.S. Department of Energy. M.P. acknowledges financial support from NSF (CMMI 1638327). This work was funded (in part) by the Dowd Fellowship from the College of Engineering at Carnegie Mellon University.

ORCID iD

Kaushik Dayal

References

Petra

Stadler

Model variational inverse problems governed by partial differential equations. Technical report, Institute for Computational Engineering and Sciences, Texas University at Austin, Austin, TX, 2011.

Martin

Wilcox

Burstedde

, et al. A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion. SIAM J Sci Comput 2012; 34(3): A1460–A1487.

Bui-Thanh

Ghattas

Martin

, et al. A computational framework for infinite-dimensional Bayesian inverse problems, part I: the linearized case, with application to global seismic inversion. SIAM J Sci Comput 2013; 35(6): A2494–A2523.

Petra

Martin

Stadler

, et al. A computational framework for infinite-dimensional Bayesian inverse problems, part II: stochastic Newton MCMC with application to ice sheet flow inverse problems. SIAM J Sci Comput 2014; 36(4): A1525–A1555.

Alghamdi

Hesse

Chen

, et al. Bayesian poroelastic aquifer characterization from InSAR surface deformation data. Part I: maximum A posteriori estimate. arXiv preprint arXiv:2002.10706, 2020.

Villa

Petra

Ghattas

hIPPYlib: An extensible software framework for large-scale inverse problems governed by PDEs; part I: deterministic inversion and linearized Bayesian inference. arXiv preprint arXiv:1909.03948, 2019.

Funke

Farrell

PE.

A framework for automated PDE-constrained optimisation. arXiv preprint arXiv:1302.3894, 2013.

Cui

Martin

Marzouk

, et al. Likelihood-informed dimension reduction for nonlinear inverse problems. Inverse Probl 2014; 30(11): 114015.

Jha

Walker

Mitchell

, et al. Mutual-information based optimal experimental design for hyperpolarized

^{13}

C-Pyruvate MRI. arXiv preprint arXiv:2206.12509, 2022.

10.

Jha

Cao

Oden

JT.

Bayesian-based predictions of COVID-19 evolution in Texas using multispecies mixture-theoretic continuum models. Comput Mech 2020; 66(5): 1055–1068.

11.

Jha

Oden

JT.

Goal-oriented A-posteriori estimation of model error as an aid to parameter estimation. J Comput Phys 2022; 470: 111575.

12.

Engl

Hanke

Neubauer

Regularization of inverse problems, vol. 375. Dordrecht: Springer Science & Business Media, 1996.

13.

Tarantola

Inverse problem theory and methods for model parameter estimation. Philadelphia, PA: SIAM, 2005.

14.

Uribe

Papaioannou

Betz

, et al. Bayesian inference of random fields represented with the Karhunen–Loève expansion. Comput Method Appl M 2020; 358: 112632.

15.

Uribe

Papaioannou

Latz

, et al. Bayesian inference with subset simulation in varying dimension applied to the Karhunen–Loève expansion. Int J Numer Meth Eng 2021; 122: 5100–5127.

16.

Gebraad

Boehm

Fichtner

Bayesian elastic full-waveform inversion using Hamiltonian Monte Carlo. J Geophys Res: Sol Ea 2020; 125(3): e2019JB018428.

17.

Arridge

Kaipio

Kolehmainen

, et al. Approximation errors and model reduction with an application in optical diffusion tomography. Inverse probl 2006; 22(1): 175.

18.

Bui-Thanh

Willcox

Ghattas

Parametric reduced-order models for probabilistic analysis of unsteady aerodynamic applications. AIAA J 2008; 46(10): 2520–2529.

19.

Galbally

Fidkowski

Willcox

, et al. Non-linear model reduction for uncertainty quantification in large-scale inverse problems. Int J Numer Meth Eng 2010; 81(12): 1581–1608.

20.

Patsialis

Kyprioti

Taflanidis

AA.

Bayesian calibration of hysteretic reduced order structural models for earthquake engineering applications. Eng Struct 2020; 224: 111204.

21.

Betancourt

. A conceptual introduction to Hamiltonian Monte Carlo. arXiv preprint arXiv:1701.02434, 2017.

22.

Lee

Vempala

. Convergence rate of Riemannian Hamiltonian Monte Carlo and faster polytope volume computation. In: Proceedings of the 50th annual ACM SIGACT symposium on theory of computing, Los Angeles, CA, 25–29 June 2018, pp. 1115–1121 New York, NY: Association for Computing Machinery.

23.

Bui-Thanh

Girolami

Solving large-scale PDE-constrained Bayesian inverse problems with Riemann manifold Hamiltonian Monte Carlo. Inverse Probl 2014; 30(11): 114014.

24.

Minka

TP.

Hessian-based Markov chain Monte-Carlo algorithms. In: 1st cape cod workshop on Monte Carlo methods Citeseer, 2002 http://www.cs.purdue.edu/homes/alanqi/papers/qiminka-HMH-AMIT-02.pdf.

25.

Fichtner

Zunino

Gebraad

A tutorial introduction to the Hamiltonian Monte Carlo solution of weakly nonlinear inverse problems. EarthArXiv, 2018.

26.

Wang

Mohamed

Freitas

Adaptive Hamiltonian and Riemann manifold Monte Carlo. In: International conference on machine learning, Atlanta, GA, 16–21 June. 2013, pp. 1462–1470. JMLR.org.

27.

Zhang

Sutton

CA.

Quasi-Newton methods for Markov chain Monte Carlo. In: Jordan

LeCun

Solla

(eds) Advances in neural information processing systems. Cambridge, MA: MIT Press, 2011, pp. 2393–2401.

28.

Roberts

Rosenthal

JS.

Optimal scaling of discrete approximations to Langevin diffusions. J Roy Stat Soc B 1998; 60(1): 255–268.

29.

Sumy

Neighbors

Cochran

, et al. Low stress drops observed for aftershocks of the 2011 mw 5.7 Prague, Oklahoma, earthquake. J Geophys Res: Sol Ea 2017; 122(5): 3813–3834.

30.

Schoenball

Ellsworth

WL.

A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and southern Kansas. J Geophys Res: Sol Ea 2017; 122(12): 10–189.

31.

Keranen

Weingarten

Induced seismicity. Annu Rev Earth Pl Sc 2018; 46: 149–174.

32.

Segall

Injection-induced seismicity: poroelastic and earthquake nucleation effects. J Geophys Res: Sol Ea 2015; 120(7): 5082–5103.

33.

Wang

Pozzi

Small

, et al. Statistical method for early detection of changes in seismic rate associated with wastewater injections. B Seismol Soc Am 2015; 105(6): 2852–2862.

34.

Wang

Small

Harbert

, et al. A Bayesian approach for assessing seismic transitions associated with wastewater injections. B Seismol Soc Am 2016; 106(3): 832–845.

35.

Langenbruch

Weingarten

Zoback

MD.

Physics-based forecasting of man-made earthquake hazards in Oklahoma and Kansas. Nat Commun 2018; 9(1): 1–10.

36.

Keranen

Savage

Abers

, et al. Potentially induced earthquakes in Oklahoma, USA: links between wastewater injection and the 2011 mw 5.7 earthquake sequence. Geology 2013; 41(6): 699–702.

37.

Dong

, et al. Co2 sequestration in depleted oil and gas reservoirs – caprock characterization and storage capacity. Energ Convers Manage 2006; 47(11-12): 1372–1382.

38.

Shapiro

SA.

Fluid-induced seismicity. Cambridge: Cambridge University Press, 2015.

39.

Broccardo

Mignan

Grigoli

, et al. Induced seismicity risk analysis of the hydraulic stimulation of a geothermal well on Geldinganes, Iceland. Nat Hazard Earth Sys 2020; 20(6): 1573–1593.

40.

Dashti

Stuart

AM.

The Bayesian approach to inverse problems. arXiv preprint arXiv:1302.6989, 2013.

41.

Pandit

Infield

Comparative analysis of gaussian process power curve models based on different stationary covariance functions for the purpose of improving model accuracy. Renew Energ 2019; 140: 190–202.

42.

Cressie

Wikle

CK.

Statistics for spatio-temporal data. Hoboken, NJ: John Wiley & Sons, 2015.

43.

Murphy

KP.

Machine learning: a probabilistic perspective. Cambridge, MA: MIT Press, 2012.

44.

Girolami

Calderhead

Riemann manifold Langevin and Hamiltonian Monte Carlo methods. J Roy Stat Soc B 2011; 73(2): 123–214.

45.

Roberts

Tweedie

RL.

Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli 1996; 2: 341–363.

46.

Gelman

Carlin

Stern

, et al. Bayesian data analysis. Boca Raton, FL: CRC Press, 2013.

47.

Neal

RM.

MCMC using Hamiltonian dynamics. In: Brooks

Gelman

Jones

, et al. (eds) Handbook of Markov chain Monte Carlo. Boca Raton, FL: CRC Press, 2011, pp. 113–162.

48.

Betancourt

. A general metric for Riemannian manifold Hamiltonian Monte Carlo. In: Nielsen

Barbaresco

(eds) International conference on geometric science of information. Berlin: Springer, 2013, pp. 327–334.

49.

Karimi

Massoudi

Walkington

, et al. Energetic formulation of large-deformation poroelasticity. Int J Numer Anal Met 2022; 46(5): 910–932.