hIPPYlib-MUQ: A Bayesian Inference Software Framework for Integration of Data with Complex Predictive Models under Uncertainty

2.1 Infinite-dimensional Bayesian Inverse Problems

The objective of the inverse problem is to determine an unknown input parameter field m that would give rise to given observational (or experimental) noisy data \(\boldsymbol {d}\) by means of a (physics-based) mathematical model. In other words, given \(\boldsymbol {d}\in \mathbb {R}^q\), we seek to infer \(m\in \mathcal {M}\) such that (1) \(\begin{equation} \boldsymbol {d}\approx \mathcal {F}(m), \end{equation}\) where \(\mathcal {F}: \mathcal {M}\rightarrow \mathbb {R}^q\) is the parameter-to-observable map that predicts observations from a given parameter m through a forward mathematical model, and \(\mathcal {M}\) is an infinite-dimensional Hilbert space of functions defined on a domain \(\mathcal {D} \subset \mathbb {R}^d\). Note that the evaluation of this map involves solving the forward PDE model given m, followed by extracting the observations from the solution of the forward problem. In what follows, we assume that the forward equation residual is continuously Fréchet-differentiable and its Jacobian a continuous linear operator with continuous inverse [Ghattas and Willcox 2021].

In the Bayesian approach, the inverse problem is framed as a statistical inference problem. The uncertain parameter m and the observational data \(\boldsymbol {d}\) are deemed as random variables and the solution is a conditional probability distribution \(\mu _{\mbox{post}}(m| \boldsymbol {d})\) that represents level of confidence in the estimation of the parameter given the data. The approach combines a prior model reflecting our prior knowledge or assumptions about the parameters before data are acquired, and a likelihood model representing the probability that a given set of parameters might give rise to the observed data.

Using the Radon-Nikodym derivative [Williams 1991] of the posterior measure \(\mu _{\mbox{post}}(m)\) with respect to the prior measure \(\mu _{\mbox{prior}}(m)\), Bayes’ theorem in infinite dimensions is stated as (2) \(\begin{equation} \frac{d\mu _{\mbox{post}}}{d\mu _{\mbox{prior}}} \propto \pi _{\mbox{like}}(\boldsymbol {d}| m), \end{equation}\) where \(\pi _{\mbox{like}}\) denotes the likelihood function. For detailed conditions under which the posterior measure is well-defined, we refer the reader to Stuart [2010].

Without loss of generality, the probability density function of the likelihood can be expressed as (3) \(\begin{equation} \pi _{\mbox{like}}(\boldsymbol {d}| m) \propto \exp \left\lbrace -\Phi (m; \boldsymbol {d})\right\rbrace . \end{equation}\) The negative log-likelihood function \(\Phi (m; \boldsymbol {d})\) has different forms depending on how one models the noise that stems from measurement uncertainties and/or modeling errors; for example, in the case of an additive Gaussian noise model \(\boldsymbol {d}= \mathcal {F}(m) + \boldsymbol {\eta }\) with a Gaussian noise random variable \(\boldsymbol {\eta } \in \mathbb {R}^q\) with zero mean and covariance matrix \(\boldsymbol {\Gamma }_{\!\mbox{noise}} \in \mathbb {R}^{q \times q}\), it has the form \(\Phi (m; \boldsymbol {d}) = \frac{1}{2}\Vert \mathcal {F}(m) - \boldsymbol {d} \Vert ^2_{\boldsymbol {\Gamma }_{\!\mbox{noise}} ^{-1}}\), where \(\Vert \cdot \Vert _{\boldsymbol {\Gamma }_{\!\mbox{noise}} ^{-1}}\) denotes the \(L^2\) norm weighted by the inverse noise covariance \(\boldsymbol {\Gamma }_{\!\mbox{noise}} ^{-1}\).

We take the prior to be a Gaussian measure,¹ i.e., \(\mu _{\mbox{prior}}= \mathcal {N}\!({m_{\text{pr}}}, {\mathcal {C}_{\text{{prior}}}})\), and assume that samples from the prior distribution are square-integrable functions in the domain \(\mathcal {D}\), i.e., belong to \(L^2(\mathcal {D})\). The prior covariance operator \(\mathcal {C}_{\text{{prior}}}\) is constructed to be a trace-class operator to guarantee bounded variance of samples from the prior distribution and well-posedness of the Bayesian inverse problem [Stuart 2010; Bui-Thanh et al. 2013; Villa et al. 2021] for detailed explanation. Specifically, we take the prior covariance operator as the inverse of the vth power of a Laplacian-like operator, namely \(\mathcal {C}_{\text{{prior}}}:= \mathcal {A}^{-v} = (-\gamma \Delta + \delta I)^{-v}; v \gt \frac{d}{2}\). Here \(\gamma\) and \(\delta \gt 0\) control the correlation length and the pointwise variance of the prior operator [Villa et al. 2021; Lindgren et al. 2011]. These choices of prior ensure that \(\mathcal {C}_{\text{{prior}}}\) is a trace-class operator, guaranteeing bounded pointwise variance and a well-posed infinite-dimensional Bayesian inverse problem [Stuart 2010; Bui-Thanh et al. 2013].

2.2 Discretization of the Bayesian Formulation

Here, we briefly present the finite-dimensional approximation of the Bayesian formulation described in the previous section. We consider a finite-dimensional subspace \(\mathcal {M}_h\) of \(\mathcal {M}\), defined by the span of a set of basis functions \(\lbrace \phi _j \rbrace _{j=1}^n\). For example, in hIPPYlib-MUQ, these basis functions are globally continuous piecewise polynomials on each element of a mesh discretization of the domain \(\mathcal {D}\) [Becker et al. 1981; Strang and Fix 1988]. These are the natural choice for the discretization of the elliptic operator \(\mathcal {A}\) used in the definition of the prior covariance. The parameter field m is then approximated as \(m\approx m_h = \sum _{j=1}^nm_j\phi _j\), and, in what follows, \(\boldsymbol {m} =\left(m_1,\dots ,m_n \right)^T\in \mathbb {R}^n\) denotes the vector of the finite element coefficients of \(m_h\).

The finite-dimensional approximation of the prior measure \(\mu _{\mbox{prior}}\) is now specified by the density (4) \(\begin{equation} \pi _{\mbox{prior}}(\boldsymbol {m}) \propto \exp \Big (-\frac{1}{2}\Vert \boldsymbol {m}-{\boldsymbol {m}}_{\mbox{pr}} \Vert ^2_{\boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1}} \Big), \end{equation}\) where \({\boldsymbol {m}}_{\mbox{pr}} \in \mathbb {R}^n\) and \(\boldsymbol {\Gamma }_{\!\mbox{prior}} \in \mathbb {R}^{n \times n}\) are the mean vector and the covariance matrix that arise upon discretization of \(m_{\text{pr}}\) and \(\mathcal {C}_{\text{{prior}}}\), respectively. We refer the reader to Bui-Thanh et al. [2013]; Villa et al. [2021] for the explicit expression of the prior covariance matrix \(\boldsymbol {\Gamma }_{\!\mbox{prior}}\).

Then the Bayes’ theorem for the density of the finite-dimensional approximation of the posterior measure \(\mu _{\mbox{post}}\) is given by (5) \(\begin{equation} \pi _{\mbox{post}}(\boldsymbol {m}) =\pi _{\mbox{post}}(\boldsymbol {m} | \boldsymbol {d}) \propto \pi _{\mbox{like}}(\boldsymbol {d}| \boldsymbol {m}) \pi _{\mbox{prior}}(\boldsymbol {m}). \end{equation}\) The finite-dimensional posterior probability density function can be expressed explicitly as (6) \(\begin{equation} \pi _{\mbox{post}}(\boldsymbol {m}) \propto \exp \Big (-\Phi (\boldsymbol {m}; \boldsymbol {d}) -\frac{1}{2}\Vert \boldsymbol {m}-{\boldsymbol {m}}_{\mbox{pr}} \Vert ^2_{\boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1}} \Big). \end{equation}\) Here, evaluating \(\Phi (\boldsymbol {m}; \boldsymbol {d})\) requires constructing a finite dimensional discretization of the parameter-to-observable map, \({\mathbf {{F}}}(\boldsymbol {m})\).

2.3 The Laplace Approximation of the Posterior Distribution

In general, the posterior probability distribution (6) is not Gaussian due to the nonlinearity of the parameter-to-observable map. In this section, we discuss the solution to the so-called linearized Bayesian inverse problem by use of the Laplace approximation. The Laplace approximation amounts to constructing a Gaussian distribution centered at the maximum a posteriori (MAP) point. The MAP point represents the most probable value of the parameter vector over the posterior, i.e., (7) \(\begin{equation} \boldsymbol {{m}_{\mbox{MAP}} } := \text{argmin}_{\boldsymbol {m}} (- \log \pi _{\mbox{post}}(\boldsymbol {m})) = \text{argmin}_{\boldsymbol {m}} \left[ \Phi (\boldsymbol {m}; \boldsymbol {d}) + \frac{1}{2}\Vert \boldsymbol {m}-{\boldsymbol {m}}_{\mbox{pr}} \Vert ^2_{\boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1}}\right]. \end{equation}\) The covariance matrix of the Laplace approximation is the inverse of the Hessian of the negative log-posterior evaluated at the MAP point. Then under the Laplace approximation, the solution of the linearized Bayesian inverse problem is given by (8) \(\begin{equation} \hat{\pi }_{\text{post}}(\boldsymbol {m}) \sim \mathcal {N}(\boldsymbol {{m}_{\mbox{MAP}} },\boldsymbol {\Gamma }_{\!\mbox{post}}), \end{equation}\) with (9) \(\begin{equation} \boldsymbol {\Gamma }_{\!\mbox{post}} = \mathbf {{H}}^{-1} (\boldsymbol {{m}_{\mbox{MAP}} }) = \left(\mathbf {{H}}_{\mbox{misfit}} (\boldsymbol {{m}_{\mbox{MAP}} }) + \boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1} \right)^{-1}, \end{equation}\) where \(\mathbf {{H}}(\boldsymbol {{m}_{\mbox{MAP}} })\) and \(\mathbf {{H}}_{\mbox{misfit}} (\boldsymbol {{m}_{\mbox{MAP}} })\) denote the Hessian matrices of, respectively, the negative log-posterior and the negative log-likelihood evaluated at the MAP point; see Villa et al. [2021] for a derivation of this Hessian using the adjoint-based method.

The quality of the Gaussian approximation of the posterior depends on the degree of nonlinearity in the parameter-to-observable map, the noise covariance matrix, and the number of observations [Gelman et al. 2004; Tarantola 2005; Bui-Thanh et al. 2013; Isaac et al. 2015; Evans and Swartz 2000; Press 2003; Stigler 1986; Tierney and Kadane 1986; Wong 2001]. When the parameter-to-observable map is linear and the additive noise and prior models are both Gaussian, the Laplace approximation is exact. Even if the parameter-to-observable map is significantly nonlinear, the Laplace approximation is a crucial ingredient to achieve scalable, efficient, and accurate posterior sampling with MCMC methods, as we will discuss in the following section.

Note that the Laplace approximation involves the Hessian of the negative log-likelihood (the data misfit part of the Hessian), which may not be explicitly constructed when the parameter dimension is large. However, the data typically provide only limited information about the parameter field, and thus the eigenspectrum of the Hessian matrix often decays very rapidly. We exploit this compact nature of the Hessian to overcome its prohibitive computational cost, and construct a low-rank approximation of the data misfit Hessian matrix using a matrix-free method (such as the randomized subspace iteration [Halko et al. 2011]).

Concretely, we construct a low-rank approximation of the data misfit Hessian, i.e., \(\mathbf {{H}}_{\mbox{misfit}} \approx \boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1} \mathbf {{V}}_r \boldsymbol {\Lambda }_r \mathbf {{V}}_r^T \boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1}\), where \(\boldsymbol {\Lambda }_r = \text{diag}(\lambda _1, \ldots , \lambda _r) \in \mathbb {R}^{r \times r}\) and \(\mathbf {{V}}_r = [\boldsymbol {v}_1, \ldots , \boldsymbol {v_r}] \in \mathbb {R}^{n \times r}\) contain the r largest eigenvalues and corresponding eigenvectors, respectively, of the generalized symmetric eigenvalue problems (10) \(\begin{equation} \mathbf {{H}}_{\mbox{misfit}} \boldsymbol {v}_i = \lambda _i \boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1} \boldsymbol {v}_i; \quad i = 1, \ldots , n. \end{equation}\) Note that the eigenvectors \(\boldsymbol {v}_i\) are orthonormal with respect to \(\boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1}\), that is \(\boldsymbol {v}_i^T \boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1} \boldsymbol {v}_j = \delta _{ij}\), where \(\delta _{ij}\) is the Kronecker delta. With this row-rank approximation and using the Sherman-Morrison-Woodbury formula [Golub and Van Loan 1996], we obtain, for the inverse of the Hessian in (9), (11) \(\begin{equation} \mathbf {{H}}^{-1} = \left(\mathbf {{H}}_{\mbox{misfit}} + \boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1}\right)^{-1} = \boldsymbol {\Gamma }_{\!\mbox{prior}} - \mathbf {{V}}_r \mathbf {{D}}_r \mathbf {{V}}_r^T + \mathcal {O}\left(\sum _{i=r+1}^n \frac{ \lambda _i }{1 + \lambda _i}\right), \end{equation}\) where \(\mathbf {{D}}_r = \operatorname{diag}(\lambda _1/(\lambda _1+1), \dots , \lambda _r/(\lambda _r+1)) \in \mathbb {R}^{r\times r}\). We refer the reader to Villa et al. [2021] for a detailed description of how the randomized algorithm [Halko et al. 2011] proceeds to construct the low-rank approximation of the Hessian, including the associated computational complexity.

We can see from the last remainder term in (11) that to obtain an accurate low-rank approximation of \(\mathbf {{H}}^{-1}\), we must keep eigenvectors corresponding to eigenvalues that are greater than 1. This approximation is used to efficiently perform various operations related to the Hessian, for example, applying the square-root inverse of the Hessian to a vector, which is needed to draw samples from the Gaussian approximation discussed in this section. We remark that the efficiency and scalability of our approach is based on the low rankness of the data misfit Hessian. We refer the reader to Ghattas and Willcox [2021], where this argument is made via model problems where low-rankness can be analytically shown, and for references to more complex problems where it can be shown empirically.

3.1 Markov Chain Monte Carlo

MCMC provides a flexible framework for exploring the posterior distribution. It generates samples from the posterior distribution that can be employed in Monte Carlo approximations of posterior expectations. For example, the posterior expectation of a quantity of interest \(\mathcal {G}(m)\) can be approximated by (12) \(\begin{equation} \int \mathcal {G}(m)\, d\mu _{\mbox{post}}\approx \frac{1}{N}\sum _{i=1}^N \mathcal {G}\left(m_i\right), \end{equation}\) where each \(m_i\sim \mu _{\mbox{post}}\) is a sample of the posterior distribution.

MCMC techniques construct ergodic Markov chains where the posterior distribution is the unique stationary distribution of the chain [Robert and Casella 2005]. Asymptotically, the states of the Markov chain are therefore exact samples of the posterior distribution and can be used in (12). Markov chains are defined in terms of a transition kernel, which is a position dependent probability distribution \(K(\cdot | \boldsymbol {m}_i)\) over state \(\boldsymbol {m}_{i+1}\) in the chain given the previous state \(\boldsymbol {m}_i\), i.e., \(\boldsymbol {m}_{i+1} \sim K(\cdot | \boldsymbol {m}_i)\). Note that chains of finite length must be employed in practice and the statistical accuracy of the Monte Carlo estimator is therefore highly dependent on the ability of the transition kernel to efficiently explore the parameter space.

There are several frameworks for constructing transition kernels that are appropriate for MCMC, including the well known Metropolis-Hastings (MH) rule [Metropolis et al. 1953; Hastings 1970], Gibbs sampler (e.g., [Casella and George 1992]), and delayed rejection (DR) [Mira et al. 2001]. MUQ provides implementations of these frameworks, as well as the generalized Metropolis-Hastings (gMH) kernel [Calderhead 2014] and multilevel MCMC framework of [Dodwell et al. 2019]. Most of these frameworks start by drawing samples from one or more proposal distributions \(q_1(\cdot | \boldsymbol {m}_i),\, \ldots ,\, q_K(\cdot | \boldsymbol {m}_i)\) that are easy to sample from (e.g., Gaussian) and then “correct” the proposed samples to obtain exact, but correlated, posterior samples. In the MH and DR kernels, corrections take the form of accepting or rejecting the proposed point. In the gMH kernel, the correction involves analytically sampling a finite state Markov chain over multiple proposed points. Intuitively, proposal distributions that capture the shape of the posterior, either locally around \(m_i\) or globally over the parameter space, tend to require fewer “corrections” and yield more efficient algorithms.

Proposal Distributions. Let \(q(\cdot | \boldsymbol {m}_i)\) denote a proposal distribution that is “parameterized” by the current state of the chain \(\boldsymbol {m}_i\). We require that the proposal distribution is easily sampled and that its density can be efficiently evaluated. The MH rule [Metropolis et al. 1953; Hastings 1970] defines a transition kernel \(K_{MH}(\cdot | \boldsymbol {m}_i)\) through a two step process: first draw a random sample \(\boldsymbol {m}^\prime \sim q(\cdot | \boldsymbol {m}_i)\) from the proposal distribution, and then accept the proposed sample \(\boldsymbol {m}^\prime\) as the next step in the chain \(\boldsymbol {m}_{i+1}\) with a probability \(\alpha = \min \lbrace 1, \gamma \rbrace\) where \(\gamma = \frac{ \pi _{\mbox{post}}(\boldsymbol {m}^\prime)}{\pi _{\mbox{post}}(\boldsymbol {m}_i)} \frac{q(\boldsymbol {m}_i | \boldsymbol {m}^\prime)}{q(\boldsymbol {m}^\prime | \boldsymbol {m}_i)}\). If rejected, set \(\boldsymbol {m}_{i+1}=\boldsymbol {m}_i\). Under mild technical conditions on the proposal distribution (see e.g., Roberts et al. [2004]), the MH rule defines a Markov chain that is ergodic and has \(\mu _{\mbox{post}}\) as a stationary distribution, thus enabling states in the chain to be used in Monte Carlo estimators. Note that the detailed balance condition (see e.g., Owen [2013]) is commonly employed to verify that a Markov chain has \(\mu _{\mbox{post}}\) as a stationary distribution, but this condition alone is not sufficient to guarantee that the chain will converge to the stationary distribution. See Roberts et al. [2004] for a detailed discussion of MH convergence and convergence rates.

While the MH rule will yield a valid MCMC kernel for a large class of proposal distributions, the dependence of the proposal on the previous state, combined with possible rejection of the proposed state, results in inter-sample correlations in the Markov chain. Because of these correlations, the error of the Monte Carlo approximation in (12) will be larger when using MCMC than in the classic Monte Carlo setting with independent samples. Markov chains with large correlations will result in larger estimator variance. To reduce correlation in the Markov chain, we seek proposal distributions that can take large steps with a high probability of acceptance. From the acceptance probability in the MH rule we see that this can occur when the proposal density \(q(\boldsymbol {m} | \boldsymbol {m}_i)\) is a good approximation to \(\pi _{\mbox{post}}(\boldsymbol {m})\), so that \(\gamma\) is close to 1.

We now turn to describing specific proposal distributions used in hIPPYlib-MUQ. First, we begin by describing common proposal mechanisms that exploit gradient and curvature information to accelerate sampling in finite-dimensional spaces. These algorithms comprise the left face of the cube in Figure 1. We then show how these ideas can be extended to construct proposals with performance that is independent of mesh-refinement (i.e., independent of dimension), thus “lifting” the derivative-accelerated proposals to an infinite-dimensional setting. This “lifting” operation transforms proposals on the left face in Figure 1 to their dimension-independent analogs on the right face of the proposal cube.

Fig. 1.

View Figure

Exploiting Gradient and Curvature Information. Perhaps the simplest and most common, but not generally efficient, proposal distribution takes the form of a Gaussian distribution centered at the current state in the chain, (13) \(\begin{equation} q_{\text{RW}}(\boldsymbol {m} | \boldsymbol {m}_i) = \mathcal {N}\left(\boldsymbol {m}_i, \boldsymbol {\Gamma }_{\text{prop}}\right), \end{equation}\) where \(\boldsymbol {\Gamma }_{\text{prop}} \in \mathbb {R}^{n \times n}\) is a user defined covariance matrix. When used with the MH rule, this random walk(RW) proposal yields an MCMC algorithm that is commonly called the random walk Metropolis algorithm. The adaptive Metropolis (AM) algorithm employs a variant of this proposal where the covariance \(\boldsymbol {\Gamma }_{\text{prop}}\) is adapted based on previous samples [Haario et al. 2001]. A proposal covariance \(\boldsymbol {\Gamma }_{\text{prop}}\) that matches posterior covariance increases efficiency, but the random walk proposal is still a poor approximation of the posterior density.

A slightly more efficient proposal can be obtained through a one-step Euler-Maruyama discretization of the Langevin stochastic differential equations [Roberts and Stramer 2003]. The resulting Langevin proposal takes the form (14) \(\begin{equation} q_{\text{MALA}}(\boldsymbol {m} | \boldsymbol {m}_i) = \mathcal {N}\left(\boldsymbol {m}_i + \tau \boldsymbol {\Gamma }_{\text{prop}} \nabla \log \pi _{\mbox{post}}(\boldsymbol {m}_i),\, 2\tau \boldsymbol {\Gamma }_{\text{prop}} \right), \end{equation}\) where \(\tau\) is the step size parameter. MH samplers with this proposal are called Metropolis-adjusted Langevin algorithms(MALA). Like the AM algorithm, adapting the covariance of the MALA proposal can also improve performance [Atchadé 2006; Marshall and Roberts 2012]. It is also common to use an approximation of the posterior covariance, such as the inverse of the log-posterior Hessian, to help the MALA proposal capture the posterior correlation. In this work for example, we employ a low rank-based approximation of the log-posterior Hessian at the MAP point (c.f. Equation (11)) (15) \(\begin{equation} q_{\text{H-MALA}}(\boldsymbol {m} | \boldsymbol {m}_i) = \mathcal {N}\left(\boldsymbol {m}_i + \tau \mathbf {{H}}^{-1} \nabla \log \pi _{\mbox{post}}(\boldsymbol {m}_i), 2\tau \mathbf {{H}}^{-1}\right). \end{equation}\) This metric is similar to the one used by Martin et al. [2012] and is equivalent to the preconditioned MALA proposal in (14) using the covariance of the Laplace approximation in (9).

Both (13) and (14) use a covariance that is constant across the parameter space. Allowing this covariance to adapt to the local correlation structure of the posterior density enables higher order approximations to be obtained, resulting in more efficient MCMC algorithms. In Girolami and Calderhead [2011], a differential geometric viewpoint was employed to define a family of proposal mechanisms on a Riemannian manifold. Adapting the MALA proposal in (14) to this manifold setting and ignoring the manifold’s curvature, results in (16) \(\begin{equation} q_{\text{sMMALA}}(\boldsymbol {m} | \boldsymbol {m}_i) = \mathcal {N}\left(\boldsymbol {m}_i + \tau \mathbf {G}^{-1}(\boldsymbol {m}_i) \nabla \log \pi _{\mbox{post}}(\boldsymbol {m}_i), 2\tau \mathbf {G}^{-1}(\boldsymbol {m}_i)\right), \end{equation}\) where \(\mathbf {G}(\boldsymbol {m})\) is a position-dependent metric tensor. This is known as the simplified Manifold MALA (sMMALA) proposal. Girolami and Calderhead [2011] defined the metric tensor \(\mathbf {G}(\boldsymbol {m})\) using the expected Fisher information metric, which provides a positive definite approximation of the posterior covariance at the point \(\boldsymbol {m}\).

Hamiltonian Monte Carlo techniques, including the No-U-Turn Sampler (NUTS) [Hoffman and Gelman 2014], are another important class of MCMC proposals. These techniques approximately solve a Hamiltonian system to take large jumps in parameter space. While efficient in many scenarios (see e.g., Neal [2010]), especially with purely statistical models, we have found that solving the Hamiltonian system involves an intractable number of posterior gradient evaluations on our PDE-based problems of interest. The transport map MCMC algorithms of Parno and Marzouk [2018] are also not considered here because of the challenge of building high-dimensional transformations.

Dimension-Independent Proposal Distributions. For finite-dimensional parameters, the random walk and MALA proposals defined above can be used with the MH rule for MCMC. However, their performance is not discretization invariant. As the discretization of the function m is refined, the performance of the samplers on the finite-dimensional posterior \(\pi _{\mbox{post}}(\boldsymbol {m})\) will worsen. As the dimension increases, the difference between the largest two eigenvalues of the MCMC transition kernel (i.e., the spectral gap), goes to zero and the mixing times of the chains grows indefinitely; see Hairer et al. [2014]; Cotter et al. [2012] for details. Some modifications to the proposals are necessary to obtain “dimension-independent” performance. The works of Cotter et al. [2012], Beskos et al. [2017], and Bardsley et al. [2020], for example, modify existing finite-dimensional proposals to ensure the algorithm performance is independent of mesh refinement.

The dimension-independent analog of the RW proposal is the preconditioned Crank-Nicolson (pCN) proposal introduced in Cotter et al. [2012]. It takes the form (17) \(\begin{equation} q_{\text{pCN}}(\boldsymbol {m} | \boldsymbol {m}_i) = \mathcal {N}\left({\boldsymbol {m}}_{\mbox{pr}} + \sqrt {1-\beta ^2}(\boldsymbol {m}_i-{\boldsymbol {m}}_{\mbox{pr}}),\, \beta ^2 \boldsymbol {\Gamma }_{\!\mbox{prior}} \right). \end{equation}\) Notice that when \(\beta =1\), the pCN proposal is equal to the prior distribution. The MALA proposal was also adapted in Cotter et al. [2012] to obtain the infinite-dimensional MALA (\(\infty\)-MALA) proposal (18) \(\begin{equation} q_{\text{MALA}}^{\infty }(\boldsymbol {m} | \boldsymbol {m}_i) = \mathcal {N}\left(\sqrt {1 - \beta ^2}\boldsymbol {m}_i + \beta \frac{\sqrt {h}}{2} \left({\boldsymbol {m}}_{\mbox{pr}} - \boldsymbol {\Gamma }_{\!\mbox{prior}} \nabla \Phi (\boldsymbol {m}_i)\right), \, \beta ^2 \boldsymbol {\Gamma }_{\!\mbox{prior}} \right), \end{equation}\) where \(\beta = 4\sqrt {h} / (4 + h)\) and h is a parameter that can be tuned. While the pCN and \(\infty\)-MALA proposals result in discretization-invariant Metropolis-Hastings algorithms, they suffer from the same deficiencies as their finite-dimensional RW and MALA analogs, i.e., they do not capture the posterior geometry.

Several efforts have worked to minimize this deficiency, see for example Beskos et al. [2017]; Rudolph and Sprungk [2018]; Pinski et al. [2015]; Petra et al. [2014]. We consider a generalization of the pCN proposal described in Pinski et al. [2015]. It incorporates the MAP point and the posterior curvature information at that point into the pCN proposal, which is denoted by H-pCN and takes the form (19) \(\begin{equation} q_{\text{H-pCN}}(\boldsymbol {m} | \boldsymbol {m}_i) = \mathcal {N}\left(\boldsymbol {{m}_{\mbox{MAP}} } + \sqrt {1-\beta ^2}(\boldsymbol {m}_i-\boldsymbol {{m}_{\mbox{MAP}} }),\, \beta ^2 \mathbf {{H}}^{-1} \right). \end{equation}\) Another method that can exploit the posterior geometry is an extension of the \(\infty\)-MALA proposal discussed in Beskos et al. [2017]: (20) \(\begin{equation} q_{\text{sMMALA}}^{\infty }(\boldsymbol {m} | \boldsymbol {m}_i) = \mathcal {N}\left(\mu ^\prime (\boldsymbol {m}_i),\, \Gamma ^\prime (\boldsymbol {m}_i)\right), \end{equation}\) where (21) \(\begin{eqnarray} \mu ^\prime (\boldsymbol {m}_i) &=& \sqrt {1 - \beta ^2}\boldsymbol {m}_i + \beta \frac{\sqrt {h}}{2} \left(\boldsymbol {m}_i - \mathbf {G}^{-1} \boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1} (\boldsymbol {m}_i-{\boldsymbol {m}}_{\mbox{pr}}) - \mathbf {G}^{-1} \nabla \Phi (\boldsymbol {m}_i) \right) \end{eqnarray}\) (22) \(\begin{eqnarray} \Gamma ^\prime (\boldsymbol {m}_i) &=& \beta ^2 \mathbf {G}^{-1}(\boldsymbol {m}_i). \end{eqnarray}\) This \(\infty\)-sMMALA proposal simplifies to \(\infty\)-MALA when \(\mathbf {G}^{-1}(\boldsymbol {m}_i)=\boldsymbol {\Gamma }_{\!\mbox{prior}}\). When \(\mathbf {G}(\boldsymbol {m}_i)\) is the Laplace approximation Hessian from (9), the \(\infty\)-sMMALA proposal simplifies to (23) \(\begin{equation} q_{\text{H-MALA}}^\infty (\boldsymbol {m} | \boldsymbol {m}_i) = \mathcal {N}\left(\sqrt {1 - \beta ^2}\boldsymbol {m}_i + \beta \frac{\sqrt {h}}{2} \left(\boldsymbol {m}_i - \mathbf {{H}}^{-1} \boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1} (\boldsymbol {m}_i-{\boldsymbol {m}}_{\mbox{pr}}) - \mathbf {{H}}^{-1} \nabla \Phi (\boldsymbol {m}_i) \right), \beta ^2 \mathbf {{H}}^{-1} \right), \end{equation}\) which we denote by H-\(\infty\)-MALA.

Alternative Transition Kernels. The proposal distributions above are classically considered in the context of a Metropolis-Hasting kernel. However, there are alternative transition kernels that also result in ergodic Markov chains. Here we consider transition kernels constructed from the delayed rejection approach of Mira et al. [2001] as well as Metropolis-within-Gibbs kernels, which repeatedly use the Metropolis-Hastings rule on different conditional slices of the posterior distribution to construct the Markov chain. In particular, we consider the family of dimension-independent likelihood-informed (DILI) approaches [Cui et al. 2016a; Cui and Zahm 2021], which define a Metropolis-within-Gibbs sampler that inherits dimension-independent properties from an appropriate dimension-independent proposal. By dimension-independence here we mean that the acceptance rate and mixing properties will not deteriorate when the dimension of the problem increases.

The delayed rejection kernel allows multiple proposals to be attempted in each step of the Markov chain. This can be advantageous when using multiple proposals with complementary properties. For example, it is possible to start with a proposal that attempts to make large ambitious jumps across the parameter space but may have low acceptance probability while falling back on a more conservative proposal that takes smaller steps with a larger probability of acceptance. Similarly, it is possible to start with a proposal that is more computationally efficient (e.g., does not require gradient information) but less likely to be accepted, while employing a more expensive proposal mechanism in a second stage to ensure the chain explores the space. In either case, if the first proposed move is rejected by the Metropolis-Hastings rule, another more expensive proposal that is more likely to be accepted can be tried with an adjusted acceptance probability. More than two stages can also be employed.

DILI divides the parameter space into a finite-dimensional subspace, which can be explored with standard proposal mechanisms, and a complementary infinite-dimensional space that can be explored with a dimension-independent approach, such as those described above. The resulting transition kernel is more complicated than the Metropolis-Hastings rule, but inherits the dimension-independent properties of the complementary space proposal. The likelihood-informed subspace is computed using the generalized eigenvalue problem (10). If the eigenvalue is larger than one, it indicates that the likelihood function dominates the prior density in that direction. The same low rank structure used to approximate the posterior Hessian can therefore be used to decompose the parameter space into a likelihood-informed subspace (LIS) spanned by the columns of \(\mathbf {{V}}_r\) and an orthogonal complementary space (CS). Within each subspace, a standard Metropolis-Hastings kernel is employed. As long as the kernel in the CS uses a dimension-independent proposal (typically pCN), then the DILI sampler will remain dimension-independent. Unlike the original implementation described in Cui et al. [2016a], the MUQ implementation does not use a whitening transform and thus avoids computing any symmetric decomposition of the prior covariance. In general, the Hessian used in (10) can be adapted to capture more correlation structure. However, we did not find this necessary in the numerical experiments below.

Assembling an MCMC Algorithm. It is possible to combine nearly any of the proposals and kernels described above, resulting in myriad possible MCMC algorithms. As suggested in Figure 2, there are three fundamental building blocks to an MCMC algorithm. The chain keeps track of previous points and allows computing Monte Carlo estimates. The kernel defines a mechanism for sampling the next state \(\boldsymbol {m}_{i+1}\) given the value of the current state \(\boldsymbol {m}_i\) and one or more proposal distributions. The proposal defines a position specific probability distribution that can be easily sampled and has a density that can be efficiently evaluated. We mimic these abstract interfaces in our software design to define and test a large number of kernel-proposal combinations.

Fig. 2.

View Figure

3.2 MCMC Diagnostics

Two questions naturally arise when analyzing a length I Markov chain \([\mathbf {m}_1,\ldots ,\mathbf {m}_{I}]\) produced by MCMC. First, has the chain converged to the stationary distribution? Second, what is the statistical efficiency of the chain, that is, how many independent samples does the chain have that actually contribute to the accuracy of Monte Carlo estimators? Most theoretical guarantees are asymptotic, and it is important to quantitatively answer these questions when employing finite-length MCMC chains. Based on these considerations, this section describes the diagnostics implemented in hIPPYlib-MUQ to check the convergence and statistical efficiency of high-dimensional MCMC chains.

5.1 Coefficient Field Inversion in a Two-dimensional Poisson Linear PDE

We first consider the coefficient field inversion in a Poisson partial differential equation given pointwise noisy state measurements. We begin by describing the forward model setup and quantity of interest, followed by the definition of the prior and the likelihood distributions. Next, we present the Laplace approximation of the posterior and apply several MCMC methods to characterize the posterior distribution, as well as the predictive posterior distribution of a scalar quantity of interest. Finally, the scalability of the proposed methods with respect to the parameter dimension is then assessed in a mesh refinement study.

5.1.2 Prior Model.

As discussed in Section 2, we choose the prior to be a Gaussian distribution \(\mathcal {N}\!({m_{\text{pr}}}, {\mathcal {C}_{\text{{prior}}}})\) with \(\mathcal {C}_{\text{{prior}}}= \mathcal {A}^{-2}\) where \(\mathcal {A}\) is a Laplacian-like operator given as (35) \(\begin{equation} \mathcal {A}m = {\left\lbrace \begin{array}{ll} - \gamma \nabla \cdot (\Theta \nabla m) + \delta m & \text{in } \Omega ,\\ \Theta \nabla m \cdot \mathbf {n} + \beta m & \text{on } \partial \Omega . \end{array}\right.} \end{equation}\)

Here, \(\beta \propto \sqrt {\gamma \delta }\) is the optimal Robin coefficient introduced to alleviate undesirable boundary effects [Daon and Stadler 2018], and an anisotropic tensor \(\Theta\) is of the form (36) \(\begin{equation} \Theta = \begin{bmatrix} \theta _1 \sin ^2\alpha + \theta _2 \cos ^2\alpha & (\theta _1 - \theta _2) \sin \alpha \cos \alpha \\ (\theta _1 - \theta _2) \sin \alpha \cos \alpha & \theta _1 \cos ^2\alpha + \theta _2\sin ^2\alpha \end{bmatrix}. \end{equation}\)

For this example we take \(\gamma = 0.1, \delta = 0.5, \beta = \sqrt {\gamma \delta } / 1.42, \theta _1 = 2.0, \theta _2 = 0.5\) and \(\alpha = \pi / 4\). Figure 6 shows the prior mean \(m_{\text{pr}}\) and three samples from the prior distribution.

Fig. 6.

View Figure

5.1.3 Observations with Noise and the Likelihood.

We generate state observations at l random locations uniformly distributed over \([0.05, 0.95]^2\) by solving the forward problem on the finest mesh with the true parameter field \(m_{\scriptscriptstyle \text{true}}\) (here a sample from the prior is used) and then adding a random Gaussian noise to the resulting state values; see Figure 7. The number of observations l is set to 300 for this example. The vector of synthetic observations is given by (37) \(\begin{equation} \boldsymbol {d}= \mathcal {B} u + \boldsymbol {\eta }, \end{equation}\) where \(\mathcal {B}\) is a linear observation operator, restricting the state solution to the l observation points. The additive noise vector \(\boldsymbol {\eta }\) has mutually independent components that are normally distributed with zero mean and standard deviation \(\sigma = 0.005\). The likelihood function is then given by (38) \(\begin{equation} \pi _{\rm like}(\boldsymbol {d}\,| \, m) \propto \exp \left(-\frac{1}{2}\Vert \mathcal {B}\,u(m) - \boldsymbol {d}\Vert ^2_{\boldsymbol {\Gamma }_{\!\mbox{noise}} ^{-1}}\right), \end{equation}\) where \(\boldsymbol {\Gamma }_{\!\mbox{noise}} = \sigma ^2 \mathbf {I}\).

Fig. 7.

View Figure

5.1.4 Laplace Approximation of the Posterior.

We next construct the Laplace approximation of the posterior, a Gaussian distribution \(\hat{\mu }_{\mbox{post}}\sim \mathcal {N}\!({m_{\scriptscriptstyle \text{MAP}}}, {\mathcal {H}(m_{\scriptscriptstyle \text{MAP}})^{-1}})\) with mean equal to the MAP point and covariance given by the Hessian of the negative log-posterior evaluated at the MAP point. The MAP point is obtained by minimizing the negative log-posterior, i.e., (39) \(\begin{equation} \underset{m \in \mathcal {M}}{\min } \; \mathcal {J}(m) \;:=\; \frac{1}{2} \Vert \mathcal {B}\,u(m)- \boldsymbol {d}\Vert ^2_{ {\bf \Gamma }_{\text{noise}}^{-1}} +\frac{1}{2} \Vert m -m_{\text{pr}}\Vert ^2_{\mathcal {C}_{\text{{prior}}}^{-1}}. \end{equation}\)

We employ the inexact Newton-CG algorithm implemented in hIPPYlib to solve the above optimization problem. We refer the reader to Villa et al. [2021] for a detailed description of the algorithm and the expressions for the gradient and Hessian actions of the negative log-posterior \(\mathcal {J}(m)\).

As pointed out in Section 2, explicitly computing the Hessian is prohibitive for large-scale problems, as this entails solving two forward-like PDEs as many times as the number of parameters. To make the operations with the Hessian scalable with respect to the parameter dimension, we invoke a low-rank approximation of the data misfit part of the Hessian, retaining only r eigenvectors that are the most significantly informed directions from the data [Villa et al. 2021].

Figure 8 shows the eigenspectrum of the prior-preconditioned data misfit Hessian. The double pass randomized algorithm provided by \(\texttt {hIPPYlib}\) with an oversampling factor of 20 is used to accurately compute the dominant eigenpairs. We see that eigenvalues are smaller than 1 after around the \(60th\) eigenvalue, indicating that keeping 60 eigenpairs is sufficient for the low-rank approximation. Figure 8 also shows four eigenvectors, which, as expected, illustrate that eigenvectors corresponding to smaller eigenvalues display more fluctuations.

Fig. 8.

View Figure

In Figure 9, we depict the MAP point and three samples drawn from the Laplace approximation of the posterior.

Fig. 9.

View Figure

5.1.5 Exploring the Posterior using MCMC Methods.

In this section, we implement the advanced MCMC algorithms discussed in the Section 3 to explore the posterior and compare their performance.

In particular, we consider pCN, MALA, \(\infty\)-MALA, DR, DILI, and their Hessian-informed counterparts. For each method, we simulate 20 independent MCMC chains, each with 25,000 samples (after discarding 2,500 samples as burn-in), and hence draw a total of 500,000 samples from the posterior. A sample from the Laplace approximation of the posterior is chosen as a starting point for each chain.

For checking the convergence and statistical efficiency of MCMC chains, we consider the subspace spanned by the r dominant eigenvectors of the generalized eigensystem in (10), instead of all components of the parameter vector \(\boldsymbol {m}\). Specifically, we compute the MPSRF, autocorrelation time, and ESS with respect to a coefficient vector \(\mathbf {c} \in \mathbb {R}^r\) defined by (40) \(\begin{equation} \mathbf {c} = \mathbf {V}_r^T \boldsymbol {\Gamma }_{\!\mbox{prior}} ^{-1} \mathbf {m}. \end{equation}\)

Table 1 shows the convergence diagnostics and computational efficiency of the MCMC samples. MPSRF and ESS are computed with respect to the projection of parameter samples along the first 25 dominant eigenvectors of the prior-preconditioned data misfit Hessian at the MAP point. Table 1 reports the mininum, maximum, and average ESS over all the 25 projections.

Table 1.

Method	AR (%)	MPSRF	Min. ESS (index)	Max. ESS (index)	Avg. ESS	NPS/ES
pCN (5.0E-3)	24	2.629	25 (24)	225 (8)	84	5,952
MALA (6.0E-6)	48	2.642	26 (22)	874 (5)	148	10,135
\(\infty\)-MALA (1.0E-5)	57	2.943	25 (23)	1,102 (5)	160	9,375
H-pCN (4.0E-1)	27	1.192	64 (1)	3,598 (15)	2,314	216
H-MALA (6.0E-2)	60	1.014	545 (1)	8,868 (19)	6,459	232
H-\(\infty\)-MALA (1.0E-1)	71	1.016	582 (1)	8,417 (18)	5,905	254
DR (H-pCN (1.0E0), H-MALA (6.0E-2))	(4, 61)	1.013	641 (1)	12,522 (17)	9,222	215
DR (H-pCN (1.0E0), H-\(\infty\)-MALA (2.0E-1))	(4, 48)	1.011	613 (1)	12,812 (17)	9,141	213
DILI-PRIOR (0.8, 0.1)	(60, 33)	1.064	314 (1)	4,667 (13)	3,216	548
DILI-LA (0.8, 0.1)	(83, 36)	1.017	562 (1)	10,882 (17)	7,192	245
DILI-MAP (0.8, 0.1)	(77, 22)	1.006	1,675 (1)	10,271 (20)	8,692	202

• Acceptance rate (AR), multivariate potential scale reduction factor (MPSRF), and effective sample sample size (ESS) are reported for convergence diagnostics. MPSRF and ESS are computed with respect to the projection of parameter samples along the first 25 dominant eigenvectors of the prior-preconditioned data misfit Hessian at the MAP point. Two values of AR are listed in DR and DILI-MAP, which are for the first and the second proposal moves, respectively. We also provide the number of forward and/or adjoint PDE solves per effective sample (NPS/ES) for sampling efficiency. We use 20 MCMC chains, each with 25,000 iterations (500,000 samples in total). The numbers in parentheses in each method name represent the parameter values used (\(\beta\) for pCN, \(\tau\) for MALA, \(h\) for \(\infty\)-MALA, and \(\beta\) and \(\tau\) for and DILI). The numbers in parentheses of the minimum ESS and the maximum ESS indicate the corresponding eigenvector index.

View Table

Table 1. Comparison of the Performance of Several MCMC Methods for the Poisson Problem: pCN, MALA, \(\infty\) -MALA, DR, DILI, and their Hessian-informed Versions

• Acceptance rate (AR), multivariate potential scale reduction factor (MPSRF), and effective sample sample size (ESS) are reported for convergence diagnostics. MPSRF and ESS are computed with respect to the projection of parameter samples along the first 25 dominant eigenvectors of the prior-preconditioned data misfit Hessian at the MAP point. Two values of AR are listed in DR and DILI-MAP, which are for the first and the second proposal moves, respectively. We also provide the number of forward and/or adjoint PDE solves per effective sample (NPS/ES) for sampling efficiency. We use 20 MCMC chains, each with 25,000 iterations (500,000 samples in total). The numbers in parentheses in each method name represent the parameter values used (\(\beta\) for pCN, \(\tau\) for MALA, \(h\) for \(\infty\)-MALA, and \(\beta\) and \(\tau\) for and DILI). The numbers in parentheses of the minimum ESS and the maximum ESS indicate the corresponding eigenvector index.

The last column in Table 1 represents the number of forward and/or adjoint PDE solves required to draw a single independent sample (average ESS is used). This quantity can be used to measure the sampling efficiency and rank the methods in terms of computational efficiency. Under this metric, DILI-MAP is the most efficient method and requires only 202 PDE solves for an effective sample. DR (213 NPS/ES for H-\(\infty\)-MALA and 215 NPS/ES for H-MALA) and H-pCN (216 NPS/ES) are close seconds.

We next assess the convergence of MCMC samples of the quantity of interest in (34) to the predictive posterior distribution of \(\mathcal {G}(m)\): the autocorrelation function estimates of the quantity of interest (34) are shown in Figure 10 (here, we use the formula (29) to account for the use of multiple chains), and trace plots from three independent MCMC chains and histograms of all the MCMC samples for pCN and H-pCN are depicted in Figure 11.

Fig. 10.

View Figure

Fig. 11.

View Figure

Lastly, we compare estimates of moments of the quantity of interest for the different sampling strategies. For each MCMC chain, the kth (\(k = 1,2,3\)) moment of the quantity of interest computed from parameter samples \(\mathbf {m}_i\) (\(i=1, 2, \ldots , N; N = 25,000\)) is computed as (41) \(\begin{equation} \hat{\mathcal {G}}_k = \frac{1}{N} \sum _{i=1}^{N} \mathcal {G}^k (\mathbf {m}_i). \end{equation}\) The results are reported in Figure 12 as box-and-whisker plots.

Fig. 12.

View Figure

From the results presented in this section, we draw the following conclusions:

• The Hessian information at the MAP point plays an important role in enhancing the sampling performance of the MCMC methods. In fact, MCMC chains without the Hessian information did not converge over the entire length of the chain and were localized around the starting point. The convergence was achieved or nearly achieved only when the MCMC proposal exploited the Laplace approximation of the posterior that incorporates the Hessian information.

• DILI-MAP shows the best sampling efficiency in terms of the number of forward and/or adjoint PDE solves per effective sample. Note that the parameter value used in the MCMC methods (e.g., \(\beta\) and/or \(\tau\)) was not the optimal and a different result may be obtained with different parameter values.

We further study the performance of MCMC methods under different problem settings to provide more insight into the practical use of the hIPPYlib-MUQ framework.

5.1.6 Scalability of Hessian-informed pCN.

Here we investigate the effect of mesh resolution on the sampling performance. A curvature aware MCMC method, the H-pCN is selected with \(\beta =0.4\) for the test. The dimensions of the parameter and the state variables from a coarse mesh (mesh 1) to the finest mesh (mesh 4) are (1,089, 4,225), (4,225, 16,641), (16,641, 66,049), and (66,049, 263,169), respectively.

We follow the same problem setting as before, and use the same synthetic observations (obtained from the true parameter field generated from the finest mesh) for all levels. The top-left part of Figure 13 shows the \(r=100\) dominant eigenvalues of the prior-preconditioned data misfit Hessian. One observes that the eigenspectrum is virtually independent of mesh refinement.

Fig. 13.

View Figure

To assess the convergence of the MCMC methods, in Table 2 we report the acceptance rate, MPSRF, and ESS of the posterior samples. The MPSRF and ESS are computed with respect to the projection of parameter samples along the first 25 dominant eigenvectors of the prior-preconditioned data misfit Hessian at the MAP point. We present the autocorrelation function estimates (29) in the bottom-left part of Figure 13, and show histograms for the quantity of interest \(\mathcal {G}\) (34) in the right part of Figure 13. The results show that the convergence of samples is almost independent with respect to the MPSRF and the autocorrelation function.

Table 2.

Dimension (state, parameter)	AR (%)	MPSRF	Min. ESS (index)	Max. ESS (index)	Avg. ESS
(4,225, 1,089)	27	1.192	64 (1)	3,598 (15)	2,314
(16,641, 4,225)	24	1.333	63 (1)	3,221 (18)	1,830
(66,049, 16,641)	23	1.075	209 (1)	3,073 (11)	1,940
(263,169, 66,049)	22	1.117	102 (2)	3,276 (15)	1,767

• We use \(\beta = 0.4\) for the H-pCN method and draw in total 500,000 samples (20 MCMC chains, each with 25,000 iterations). MPSRF and ESS are computed with respect to the projection of samples along the first 25 dominant eigenvectors of the prior-preconditioned data misfit Hessian at the MAP point. The numbers in parentheses of the minimum ESS and the maximum ESS indicate the corresponding eigenvector index.

View Table

Table 2. Acceptance Rate (AR), Multivariate Potential Scale Reduction Factor (MPSRF) and Effective Sample Size (ESS) of the Posterior Samples Generated by using the H-pCN Method with Different Dimensions

• We use \(\beta = 0.4\) for the H-pCN method and draw in total 500,000 samples (20 MCMC chains, each with 25,000 iterations). MPSRF and ESS are computed with respect to the projection of samples along the first 25 dominant eigenvectors of the prior-preconditioned data misfit Hessian at the MAP point. The numbers in parentheses of the minimum ESS and the maximum ESS indicate the corresponding eigenvector index.

5.2 Coefficient Field Inversion in a Robin Boundary Condition for a Three-dimensional p-Poisson Nonlinear PDE

So far, we restricted our attention to the additive Gaussian noise model. While this additive noise model is the most commonly used model, it would be inappropriate in some cases such as speckle noise found in synthetic aperture radar (SAR) images. In this example, we consider a different noise model where the noise is proportional to the value of the observations.

The example forward model is a nonlinear PDE in three space dimensions for which we seek to infer an unknown coefficient field in a Robin boundary condition. Specifically, the forward governing equations are given by (42) \(\begin{align} -\nabla \cdot \left(|\nabla u|^{p-2}_\epsilon \nabla u\right) &= f &\quad &\text{in } \Omega , \nonumber \nonumber\\ |\nabla u|^{p-2}_\epsilon \nabla u \cdot \mathbf {n} + e^m u &= 0 &&\text{on } \partial \Omega _R, \\ |\nabla u|^{p-2}_\epsilon \nabla u \cdot \mathbf {n} &= 0 &&\text{on } \partial \Omega \setminus \partial \Omega _R, \nonumber \nonumber \end{align}\) with \(1 \le p \le \infty\). Note that the p-Laplacian, \(\nabla \cdot (|\nabla u|^{p-2} \nabla u)\), is singular when \(p\lt 2\) and degenerates when \(p\gt 2\) at points \(\nabla u =0\) [Lindqvist 2017; Brown 2010], so a regularization term \(\epsilon\) (here we take \(\epsilon = 1.0 \times 10^{-8}\)) is introduced in the above equation as \(|\nabla u|_\epsilon = \sqrt {|\nabla u|^2 + \epsilon }\). The p-Laplacian is a nonlinear counterpart of the Laplacian operator, and appears in many nonlinear diffusion problems (e.g., non-Newtonian fluids), where a nonlinear diffusion is modeled as a power law type.

We assume \(p=3\) and consider a thin brick domain \(\Omega = [0, 1]^2 \times [0, 0.05]\) with a volume source term (\(f = 1\)) and a mixed boundary condition, e.g., we impose a Robin boundary condition on the bottom boundary surface and no normal flux on the remaining boundary surfaces.

The problem domain \(\Omega\) is discretized using a regular tetrahedral grid and using linear finite elements for all the state, adjoint, and parameter variables. After discretization, the dimension is 66,564 for the state and adjoint variables, and 16,641 for the parameter variable.

For the prior, we use the same Gaussian distribution as the one employed in the previous example. We create synthetic state observations at \(l=300\) random locations uniformly distributed on the top boundary surface by first solving the forward problem with the true parameter field \(m_{\textrm {true}}\) obtained from a sample of the prior, and then multiplying a Gamma-distributed noise. Specifically, the vector of synthetic observations is of the form (43) \(\begin{equation} \boldsymbol {d}= \boldsymbol {\eta } \odot \mathcal {B} u, \end{equation}\) where \(\odot\) denotes component-wise multiplication, and each component of \(\boldsymbol {\eta }\) is independently and identically Gamma-distributed with shape \(\kappa\) and scale \(\nu\), i.e., \(\eta _i = \text{Gamma}(\kappa , \nu), i=1, \ldots , q\). We take \(\kappa = 1 / \nu\) and in this case the negative log-likelihood function has the form (44) \(\begin{equation} \Phi (m; \boldsymbol {d}) = \kappa \sum _{i=1}^{q} \left(\log (\mathcal {B} u)_i + \frac{d_i}{(\mathcal {B} u)_i} \right), \end{equation}\) where the subscript i means ith component of the corresponding vector. In this example, we set \(\kappa = 10^{4}\).

Figure 14 illustrates the true parameter field on the bottom boundary, the locations of the observations on the top surface, and the MAP point obtained by solving the optimization problem of minimizing the negative log-posterior. The Laplace approximation of the posterior is then constructed based on the low-rank factorization of the data misfit Hessian at the MAP point. The spectrum of the prior-preconditioned data misfit Hessian indicates that the number of dominant eigenvalues (larger than 1) is about 55.

Fig. 14.

View Figure

5.2.1 MCMC Results for Characterizing the Posterior.

Here we present MCMC simulation results for the uncertain boundary coefficient vector. In this example, we consider the H-pCN method with \(\beta =0.2\) and run 40 independent MCMC chains, each with 25,000 iterations after discarding 2,500 samples as burn-in (1,000,000 samples are generated in total). For each MCMC run, a sample from the Laplace approximation of the posterior is taken as the starting point.

Table 3 lists convergence diagnostics of the MCMC simulation. The parameter samples are projected onto the first 25 eigenvectors of the prior-preconditioned data misfit Hessian at the MAP point, and the MPSRF and the ESS are evaluated based on this projection.

Table 3.

AR (%)	MPSRF	Min. ESS (index)	Max. ESS (index)	Avg. ESS
23	1.041	243 (0)	2,891 (22)	1,586

• We use H-pCN method (\(\beta = 0.2\)) with 40 chains, each with 25,000 iterations (1,000,000 samples in total). The numbers in parentheses of the minimum ESS and the maximum ESS indicate the corresponding eigenvector index.

View Table

Table 3. Convergence Diagnostics for the \(p\) -Poisson Problem: Acceptance Rate (AR), Multivariate Potential Scale Reduction Factor (MPSRF), and Effective Sample Sample Size (ESS) of the Projection of the Parameter Samples along the First 25 Eigenvectors of the Prior-preconditioned Data Misfit Hessian at the MAP Point

• We use H-pCN method (\(\beta = 0.2\)) with 40 chains, each with 25,000 iterations (1,000,000 samples in total). The numbers in parentheses of the minimum ESS and the maximum ESS indicate the corresponding eigenvector index.

Figure 15 shows marginal distributions of the posterior MCMC samples and of the Laplace approximation. As before, these marginals are computed with respect to the projection of the parameter samples onto the eigenvectors. We observe that there is a clear difference between the marginal distributions of the MCMC samples and those of the Laplace approximation, especially for the eigenvector corresponding to larger eigenvalues.

Fig. 15.

View Figure

The numerical studies have been carried out on the Multi-Environment Computer for Exploration and Discovery (MERCED) cluster at UC Merced. While there are some performance variations depending on the cluster CPU node that a job is assigned to, approximately, a single PDE solve required about 0.07 seconds for the first example, and about 4 seconds for the second three-dimensional example (both are the elapsed time).