Hessian-based adaptive sparse quadrature for infinite-dimensional Bayesian inverse problems

TL;DR

This paper introduces a Hessian-based adaptive sparse quadrature method for efficiently computing infinite-dimensional integrals in Bayesian inverse problems, improving accuracy by better capturing the posterior distribution.

Contribution

The paper develops a novel Hessian-based parametrization and adaptive sparse quadrature that achieve dimension-independent convergence rates for Bayesian inverse problems.

Findings

Achieves dimension-independent convergence rates.

Demonstrates faster convergence than O(N^{-1/2}) in experiments.

Effectively approximates posterior distributions in linear and nonlinear problems.

Abstract

In this work we propose and analyze a Hessian-based adaptive sparse quadrature to compute infinite-dimensional integrals with respect to the posterior distribution in the context of Bayesian inverse problems with Gaussian prior. Due to the concentration of the posterior distribution in the domain of the prior distribution, a prior-based parametrization and sparse quadrature may fail to capture the posterior distribution and lead to erroneous evaluation results. By using a parametrization based on the Hessian of the negative log-posterior, the adaptive sparse quadrature can effectively allocate the quadrature points according to the posterior distribution. A dimension-independent convergence rate of the proposed method is established under certain assumptions on the Gaussian prior and the integrands. Dimension-independent and faster convergence than $O(N^{-1/2})$ is demonstrated for a linear as well as a nonlinear inverse problem whose posterior distribution can be effectively approximated by a Gaussian distribution at the MAP point.

Equations

$$y=\mathcal{G}(m)+\eta .$$

$${\pi }_{\text{like}}^y:=\exp \left(-\Phi (m,y)\right)=\exp \left(-\frac{1}{2}\Vert y-\mathcal{G}(m){\Vert }_{{\Gamma }_{\text{noise}}}^2\right),$$

$$||y-\mathcal{G}(m)|{|}_{{\Gamma }_{\text{noise}}}^2={\left(y-\mathcal{G}(m)\right)}^T{\Gamma }_{\text{noise}}^{-1}\left(y-\mathcal{G}(m)\right).$$

$$\frac{d{\mu }^y}{d{\mu }_0}=\frac{1}{Z}{\pi }_{\text{like}}^y,$$

$$Z={\int }_X{\pi }_{\text{like}}^{y}d{\mu }_0(m).$$

$${\mathbb{E}}^{{\mu }^y}[Q]={\int }_{X}Q(m)d{\mu }^y(m).$$

$${\mathbb{E}}^{{\mu }_0}[Q{\pi }_{\text{like}}^y]=\frac{1}{Z}{\int }_{X}Q(m)\exp \left(-\Phi (m,y)\right)d{\mu }_0(m).$$

$$m:{\mathbb{R}}^{\infty }\to X,\text{ where }m(\mathbf{\xi })=m_0+\sum _{j\ge 1}\sqrt{{\lambda }_j}{\psi }_j{\xi }_j,{\xi }_j\stackrel{i.i.d.}{\sim }\mathcal{N}(0,1),$$

$${\mathbb{E}}^{{\mu }_0}[Q{\pi }_{\text{like}}^y]=\frac{1}{Z}{\int }_{\Xi }Q(m(\mathbf{\xi }))\exp \left(-\Phi (m(\mathbf{\xi }),y)\right)d\mu (\mathbf{\xi }),$$

$$Z={\int }_{\Xi }\exp \left(-\Phi (m(\mathbf{\xi }),y)\right)d\mu (\mathbf{\xi }),$$

$$\mu (\mathbf{\xi })=\underset{j\ge 1}{⨂}\rho ({\xi }_j)d{\xi }_j,\text{ where }\rho ({\xi }_j)=\frac{1}{\sqrt{2\pi }}{e}^{-{\xi }_j^2/2}.$$

$${\mathcal{C}}_0:={\mathbb{E}}^{{\mu }_0}[(m-m_0)\otimes (m-m_0)]\equiv {\mathbb{E}}^{\mu }[(m-m_0)\otimes (m-m_0)]=\sum _{j\ge 1}{\lambda }_j{\psi }_j\otimes {\psi }_j,$$

$${\mathcal{C}}_0{\psi }_k=\sum _{j\ge 1}{\lambda }_j⟨{\psi }_j,{\psi }_k⟩{\psi }_j={\lambda }_k{\psi }_k,k=1,2,\dots ,$$

$$\mathcal{A}u:=-\beta \nabla \cdot (\Theta \nabla u)+\gamma u,u\in H^1(D),$$

$$m_1:=\underset{ϵ\to 0}{\lim }\underset{m\in X}{\mathrm{argmax}}{\mu }^y(B(m,ϵ)).$$

$$\underset{m\in E}{\min }\mathcal{J}(m)$$

$$\mathcal{J}(m)=\frac{1}{2}||y-\mathcal{G}(m)|{|}_{{\Gamma }_{\text{noise}}}^2+\frac{1}{2}||m-m_0|{|}_{{\mathcal{C}}_0}^2,$$

$$H_{\text{MAP}}:=D_m^2\mathcal{J}(m){|}_{m=m_1}=H_{\text{misfit}}+{\mathcal{C}}_0^{-1},$$

$$\mathcal{G}(m)=Gm,$$

$$H_{\text{misfit}}=G^{\ast }{\Gamma }_{\text{noise}}^{-1}G,$$

$$\mathcal{G}(m)\approx \mathcal{G}(m_1)+{\mathcal{G}}^{\prime }(m_1)(m-m_1),$$

$${\mathcal{C}}_1=H_{\text{MAP}}^{-1}={\left(H_{\text{misfit}}+{\mathcal{C}}_0^{-1}\right)}^{-1},$$

$${\mathbb{E}}^{{\mu }_1}[Q]:={\int }_{X}Q(m)d{\mu }_1(m).$$

$$m^1(\mathbf{\xi })=m_1+\sum _{j\ge 1}\sqrt{{\lambda }_j^1}{\psi }_j^1{\xi }_j,{\xi }_j\stackrel{i.i.d.}{\sim }\mathcal{N}(0,1),$$

$${\mathbb{E}}^{{\mu }_1}[Q]={\int }_{\Xi }Q(m^1(\mathbf{\xi }))d\mu (\mathbf{\xi }).$$

$${\mathbb{E}}^{{\mu }^y}[Q]=\frac{1}{Z}{\int }_{X}Q(m){\pi }_{\text{like}}^{y,1}(m)d{\mu }_1(m),\text{ where }Z:={\int }_X{\pi }_{\text{like}}^{y,1}(m)d{\mu }_1(m),$$

$${\pi }_{\text{like}}^{y,1}(m)=\exp \left(\mathcal{J}(m_1)\right)\exp \left(-\Phi (m,y)\right)\frac{d{\mu }_0(m)}{d{\mu }_1(m)}=:\exp \left(-{\mathcal{J}}_1(m)\right),$$

$${\mathcal{J}}_1(m)=\mathcal{J}(m)-\mathcal{J}(m_1)-\frac{1}{2}||m-m_1|{|}_{{\mathcal{C}}_1}^2,$$

$${\mathbb{E}}^{{\mu }^y}[Q]=\frac{1}{Z}{\int }_{\Xi }Q(m^1(\mathbf{\xi })){\pi }_{\text{like}}^{y,1}(m^1(\mathbf{\xi }))d\mu (\mathbf{\xi })\text{ and }Z={\int }_{\Xi }{\pi }_{\text{like}}^{y,1}(m^1(\mathbf{\xi }))d\mu (\mathbf{\xi }).$$

$${\mathcal{Q}}_{\Lambda }(g)=\sum _{\mathbf{\nu }\in \Lambda }{△}_{\mathbf{\nu }}(g).$$

Full text

\slugger

siscxxxxxxxx–x

1 Introduction

In many practical applications, a quantity of interest (QoI) predicted by a mathematical model of a physical system may be stochastic rather than deterministic due to uncertainties arising from inadequate knowledge of input parameters that characterize material properties, computational geometries, initial and boundary conditions, source terms, etc. Solution of an inverse problem will reduce these uncertainties based on available observational data on some system outputs, leading to more reliable computational predictions of the QoI. In a Bayesian framework, the solution of the inverse problem entails the computation of the posterior probability distribution of the uncertain input parameter conditioned on (possibly) noisy observational data of the system output. A prior probability distribution, based on experts’ belief, is prescribed for the uncertain input parameter. Then, the posterior distribution of the uncertain parameter can be formally obtained by the Radon–Nikodym derivative by Bayes’ theorem [30]. Once the posterior distribution is known, one is often interested in evaluating the statistical moments of the QoI (e.g., expectation and variance) with respect to the posterior distribution for the assessment, design, control and optimization of the system.

We consider the case of a spatially heterogeneous uncertain parameter, i.e., a spatially correlated random field with a prescribed bounded covariance between any two points of the physical domain. This setting makes the uncertain parameter infinite-dimensional, which naturally leads to an infinite-dimensional integration problem for the evaluation of the statistical moments of the QoI with respect to the posterior distribution. Several computational challenges are commonly faced for such integration problems. Traditional deterministic integration methods face the curse of dimensionality, i.e., the computational complexity grows exponentially fast with respect to the parameter dimension so that only a limited number of dimensions can be resolved. On the other hand stochastic integration methods such as Monte Carlo (MC) or Markov-chain Monte Carlo (MCMC) do not depend on the parameter dimension but on the variance of the QoI; however the convergence of the quadrature errors, $O(N^{-1/2})$ with $N$ quadrature samples [18], is often very slow. In addition, the posterior distribution, unlike the prior distribution, is usually not explicitly available but rather is implicitly represented by means of the Bayes Theorem and becomes rather concentrated in the parameter space when the observational data are very informative and thus is difficult to sample efficiently. Finally, when the system is modeled by partial differential equations (PDE) in complex geometries, evaluation of the QoI at each quadrature sample require a full PDE solve that can involve expensive large-scale computations, so that only a small number of samples can be computed.

We tackle these computational challenges by exploiting the intrinsic sparsity of the integration problem with respect to the posterior distribution and the structure of the inverse problem. By sparsity, we mean that the QoI has markedly different sensitivity in different parameter dimensions under the posterior distribution; in other words, the dependence of the QoI on different parameter dimensions is rather anisotropic with suitable parametrization of the uncertain parameter, so that we can identify the most sensitive dimensions and allocate most of the computational effort in these dimensions [28, 26, 27, 7, 8, 11]. By structure, we mean that high-order derivatives of the parameter-to-observable map with respect to the parameter at the maximum a posteriori (MAP) point, in particular the Hessian, describe the inherent structure of the manifold of the posterior distribution about the MAP point, which is especially useful when the posterior distribution is relatively concentrated at the MAP point due to informative observational data [16, 21, 3, 24]. To exploit these opportunities and to make the solution of the infinite-dimensional Bayesian inverse problem computationally feasible, here we develop a novel Hessian-based adaptive sparse quadrature for integration of a QoI with respect to the posterior distribution.

More specifically, we consider a Gaussian distribution in a Hilbert space as the prior distribution for the uncertain parameter, where the covariance operator is given by a fractional power of a second order elliptic differential operator. We then compute the MAP point by solving an optimization problem using a Lagrangian variational approach and an inexact Newton–conjugate gradient solver. The covariance of the Gaussian approximation of the posterior distribution at the MAP point, to which we refer as the Gaussian posterior covariance, is given as the inverse of the Hessian of the negative log-posterior distribution. Like the prior covariance, this covariance is of trace class; indeed we show that the eigenvalues of this covariance are less than or equal to those of the prior covariance in every dimension as long as the Hessian of the data misfit is non-negative. However, the Gaussian posterior covariance operator is formally a large dense operator, and it is not available explicitly, but implicitly through its action on a given function. To apply this covariance operator to a function is very expensive, as it requires solution of the so-called incremental forward and adjoint problems. To make actions of the Gaussian posterior covariance computationally feasible, we compute its eigendecomposition by a Lanczos [4] or a randomized [17] method. The fast decay of the eigenvalues of the Gaussian posterior covariance allows us to express the infinite-dimensional uncertain parameter by a truncated Karhunen–Loève expansion on the eigenpairs, parametrized by a vector of independent and identically distributed Gaussian random variables. Our approach is in contrast with [27], in which the uncertain parameter is first represented by a truncated Karhunen–Loève expansion on the eigenpairs of the prior covariance, and then the coordinates corresponding to the Gaussian random variables are transformed according to the Hessian of the negative log-posterior with respect to these random variables. The advantage of our approach is that we only need to compute the eigenpairs that are important according to the Gaussian posterior covariance, which can save considerable computational cost if a large number of the eigenpairs of the prior covariance have to be kept in order to capture the main uncertainty in the Gaussian posterior distribution. For instance, an extremely large number of prior eigenpairs would have to be used to capture the localized fine scales of the parameter field for the Antarctic ice sheet flow inverse problem in [19], while many fewer eigenpairs of the Gaussian posterior covariance are needed.

The sensitivity of the QoI with respect to different parameter dimensions represented by the Gaussian random variables is characterized by the fast decay of the eigenvalues of the Gaussian posterior covariance. A dimension-adaptive sparse quadrature [15, 6, 26] is then constructed for the parametric integration problem. We establish under certain assumptions on the decay of the eigenvalues of the Gaussian prior and the regularity of the integrand that the convergence of the adaptive sparse quadrature is bounded by $\mathcal{O}(N^{-s})$, independent of the number of the parameter dimensions. The exponent $s$ depends only on the algebraic decay rate of the eigenvalues; the sparser the problem is, i.e., the faster the eigenvalues decay, the faster the convergence becomes. In particular, for $s$ larger than $1/2$, the deterministic sparse quadrature achieves faster convergence than the stochastic MC or MCMC quadrature. We demonstrate the concentration of the posterior distribution in the region of the prior distribution, the efficacy of the Hessian-based allocation of the quadrature points, and the convergence property of the sparse quadrature in both a linear and a nonlinear inverse problem.

The paper is organized as follows: Section 2 presents the infinite-dimensional setting of Bayesian inverse problems, computational tasks and challenges; the Hessian-based parametrization is developed in Section 3; Section 4 describes the adaptive sparse quadrature and provides a dimension-independent convergence analysis; Section 5 is devoted to solution of two inverse problems, one linear and the other nonlinear; and conclusions are drawn in Section 6.

2 Infinite-dimensional Bayesian Inversion

2.1 Overview

Inverse problems consist of finding an unknown system input parameter given observational (or experimental) data on some system output. These problems are most often ill posed in the sense of Hadamard, i.e., violation of either existence or uniqueness, or the solution depends sensitively on the data. A classical remedy in a deterministic approach is to solve a regularized least-squares minimization problem. Alternatively, a statistical approach, formulated in the Bayesian framework, seeks for the posterior probability distribution of the parameter that maximizes the likelihood of the parameter to generate the observational data as system model output conditioned on its prior distribution that incorporates available data or expert opinion/prior belief on the parameter. The latter approach is far richer than the former as it not only computes a maximum a posterior (MAP) estimator as the most likely estimation of the parameter, but also quantifies the confidence of this estimator, and provides the probability of a different estimator other than the MAP one. Moreover, different choices of the prior distribution make it versatile to model the available data or belief on the parameter. Furthermore, for a quantity of interest (QoI) that depends on the parameter, possibly through the system state variable, more information can be assessed in the form of its statistics, such as mean, variance, and failure probability, which can be computed by either deterministic or statistical integration of the QoI with respect to the posterior distribution of the parameter.

When the parameter is spatially correlated random field in the physical domain, the solution of the statistical inverse problems is the posterior measure of the parameter in a function space, or more precisely a separable Banach space, which is typically infinite-dimensional as outlined in [30]. In this framework, an appropriate definition of the posterior measure is realized by the Radon–Nikodym derivative with respect to the prior measure, where the derivative is proportional to the likelihood function based on conditional probability distribution. The prior measure of the parameter plays an important role to guarantee the well-posedness of the statistical inverse problems. Moreover, it is the sparsity of the parameter encapsulated in the prior measure and successively inherited by the posterior measure that contributes to the development of sparse algorithms in sampling the parameter according to its posterior measure and computing the statistics of a given QoI. Among many choices of the prior measures, such as uniform, Besov, and Gaussian priors (see e.g. in [13]), we consider the last one for its popularity in practical applications.

In what follows, we present the Bayes’ formula for the definition of the posterior measure in infinite dimensions in Section 2.2, then in Section 2.3 we address the central and computationally challenging task in Bayesian inverse problems, i.e. the computation of statistics of some QoI as an infinite-dimensional integration problem. Section 2.4 is devoted to a practical construction of Gaussian priors.

2.2 Bayes’ formula for infinite-dimensional inverse problems

By $X$ we denote a separable Banach space defined over the physical domain $D$, an open and bounded subset of ${\mathbb{R}}^{d}$ ($d=1,2,3$) with Lipschitz boundary $\partial D$. By $m\in X$ we denote the uncertain system input parameter that belongs to $X$. We assume there exists a prior measure on $X$ for the parameter, denoted as $\mu_{0}$ such that $\mu_{0}(X)=1$, i.e. functions randomly drawn from the prior measure lie in $X$ with probability 1. Moreover, we assume that the system output data $y$ live in a separable Banach space $Y$, which in most practical applications is a finite-dimensional Euclidean space $Y={\mathbb{R}}^{K}$, for some $K\in{\mathbb{N}}$. By $\mathcal{G}:X\to Y$ we denote the parameter-to-observable map that maps the parameter $m$ to the observation data $y$, i.e. $y=\mathcal{G}(m)$. Evaluations of $\mathcal{G}$ involve the solution of the forward model, typically described by a partial differential equation (PDE). Furthermore, we assume that the data $y$ is corrupted by some addictive observation noise represented by a random variable $\eta\in Y$, i.e.

[TABLE]

We assume the discrepancy $\eta$ to be Gaussian distributed with zero mean and covariance $\Gamma_{\text{noise}}$, i.e. $\eta\sim\mathcal{N}(0,\Gamma_{\text{noise}})$, and we write the likelihood (i.e. the probability density function of the observation data conditioned on the parameter, $y|m$) as

[TABLE]

where $\Phi(m,y):=-\frac{1}{2}\|y-\mathcal{G}(m)\|^{2}_{\Gamma_{\text{noise}}}$ denotes the negative log likelihood or potential. Here, we defined the $\Gamma_{\text{noise}}$ weighted norm as

[TABLE]

The Bayes’ formula establishes the relation between the posterior measure $\mu^{y}$ and the prior measure $\mu_{0}$ of the parameter through the Radon–Nikodym derivative

[TABLE]

where the normalization constant $Z$ is given by

[TABLE]

A sufficient condition for the well-posedness of (4) is the quadratic boundedness and the Lipschitz continuity of the map $\mathcal{G}$ with respect to the parameter $m$, see [30]. Such conditions, which need to be verified for the specific PDE model at the hand, guarantee the absolute continuity of the posterior measure $\mu^{y}$ with respect to the prior measure $\mu_{0}$ and the positivity of the constant $Z$.

2.3 Computation of statistics and challenges

A central task of Bayesian inversion is to compute some statistics of a quantity of interest (QoI) $Q$ that depends on the parameter $m$, possibly through the solution of the PDE forward model $u$. For instance, we may be interested in computing the $k$-th moment of a functional $f$ of the solution $u$, i.e., $Q(m)=f^{k}(u(m))$, $k=1,2,\dots$. More precisely, for a given integrand $Q$, one needs to evaluate integrals with respect to the posterior measure of the form

[TABLE]

By the Bayes’ formula (4), this is equivalent to evaluate the following integral with respect to the prior measure

[TABLE]

Evaluation of integrals (6) or (7) is a notoriously computational challenge. First, the parameter is a random function in $X$ thus an infinite-dimensional integration needs to be performed, which confronts the curse-of-dimensionality for most deterministic quadrature rules or slow convergence for statistical quadrature rules. Second, each evaluation of $Q$ at a quadrature point (sample) involves the solution of the PDE forward problem, and therefore, for complex forward problems, standard algorithms for infinite-dimensional integration are computationally unfeasible due to the large number of evaluations of $Q$ required. In order to harness the total computational cost, we exploit the sparsity of the prior measure to design a novel sparse quadrature algorithm that breaks the curse-of-dimensionality and achieves fast convergence of the quadrature error.

2.4 The Gaussian prior

We first consider an abstract construction of the infinite-dimensional Gaussian prior through the pushforward of the probability measure ${\mathbb{P}}$ on the i.i.d. standard Gaussian random sequence ${\boldsymbol{\xi}}=(\xi_{j})_{j\geq 1}$ under a map that takes the sequence into the random function, as outlined in [13]. In particular, we consider a real-valued separable Hilbert space $X$ defined over the physical domain $D$, which is equipped with the inner-product $\langle\cdot,\cdot\rangle$ and the induced norm $||\cdot||_{X}$. For such space $X$, we denote its orthonormal basis as $(\psi_{j})_{j\geq 1}$. Then we can definite the map

[TABLE]

where the mean $m_{0}\in X$, the sequence ${\boldsymbol{\lambda}}:=(\lambda_{j})_{j\geq 1}\in\ell^{1}({\mathbb{R}}^{\infty})$, and $\mathcal{N}(0,1)$ represents the standard Gaussian distribution. We say that the map defined above is a random parameter in $X$ equipped with Gaussian measure. For a Gaussian prior measure $\mu_{0}$ defined through the pushforward map (8), the integral (7) can be equivalently written after the change of variables as

[TABLE]

with the normalization constant

[TABLE]

where the integration domain becomes $\Xi={\mathbb{R}}^{\infty}$, and the Gaussian measure $\mu$ on $\Xi$ is defined by the tensor product of the univariate Gaussian measure as

[TABLE]

Thus the infinite-dimensional integration (7) in the function space $X$ with measure $\mu_{0}$ becomes an infinite-dimensional parametric integration in $\Xi$ with measure $\mu$.

Before proceeding to the numerical integration with respect to the parametrized random function $m$, we have to make it precise for the construction of the orthonormal basis $(\psi_{j})_{j\geq 1}$, which are functions in the physical domain $D$. For rectangular domain $D$, e.g. $(0,1)^{d}$, tensor product of trigonometric functions or orthonormal polynomials such as Legendre polynomials could be used. For more general domain $D$ with arbitrary shape, we turn to the covariance operator of Gaussian prior for its practical construction, which, in the Hilbert space $X$, is given by [30]

[TABLE]

where we have used the parametric expression of $m$ in (8) and the i.i.d. property of $\xi$. By the orthonormality of the basis $(\psi_{j})_{j\geq 1}$, it is ready to have

[TABLE]

which implies that $(\lambda_{j},\psi_{j})_{j\geq 1}$ are the eigenpairs of the covariance operator $\mathcal{C}_{0}$. Therefore, to obtain the parametric expression (8), which is known as Karhunen–Loève expansion of the random function $m$ with Gaussian measure $\mathcal{N}(m_{0},\mathcal{C}_{0})$, we only need to specify the mean $m_{0}$ and the covariance operator $\mathcal{C}_{0}$. A commonly used covariance operator is the fractional powers of a ‘Laplacian-like’ operator $\mathcal{A}$, i.e. $\mathcal{C}_{0}=\mathcal{A}^{-\alpha}$ for certain $\alpha>0$, where $\mathcal{A}$ is positive definite, self-adjoint, and invertible. A specific example of $\mathcal{A}$, as considered in [3], is the elliptic differential operator

[TABLE]

where $H^{1}(D)$ is the conventional Hilbert space; $\beta,\gamma>0$, which control the variance of the parameter, the larger $\beta$ and $\gamma$ are, the smaller the variance becomes; $\Theta\in{\mathbb{R}}_{s.p.d.}^{d\times d}$ is symmetric positive definite, which controls the anisotropic property of the parameter; the power $\alpha>d/2$ controls the asymptotic decay rate of the eigenvalues of $\mathcal{C}_{0}$, which determines the sparsity of the parametrization of the random function $m$ under the Gaussian prior measure $\mathcal{N}(m_{0},\mathcal{C}_{0})$. In practice, these coefficients can be calibrated based on direct available data on the parameter or expert opinions.

3 Hessian-based Parametrization

Upon the parameterization of the random function as presented in the last section, evaluation of statistics for any given QoI in (6) and (7) becomes the infinite-dimensional parametric integration in (9) and (10). However, direct numerical integration of these quantities brings great challenges not only because of the curse-of-dimensionality, but also due to the fact that the posterior measure might be rather different from the prior measure, for instance the mean of the posterior measure is driven far away from that of the prior by the observational data, or the posterior measure tends to be more concentrated in a small region than the prior measure. Therefore, the numerical quadrature points allocated either deterministically or statistically according to the prior measure might not efficiently capture the posterior measure, or in the parametric setting, the posterior ‘density’ $\pi_{\text{like}}^{y}({\boldsymbol{\xi}}){\boldsymbol{\rho}}({\boldsymbol{\xi}})$ (which is not a true density but zero in the infinite-dimensional setting, here we assume that $\xi$ is finitely supported) might be extremely small where the prior ‘density’ ${\boldsymbol{\rho}}({\boldsymbol{\xi}})$ is relatively big, and vice versa. In order to tackle this challenge, we take the observation data into account and push the prior measure to the posterior measure or the one at least close to it accordingly by developing a Hessian-based parametrization method. We describe the MAP point and the Hessian of the cost functional in subsection 3.1 and 3.2, followed by the Hessian-based parametrization of the random function in subsection 3.3.

3.1 The MAP point

The MAP point characterizes the point in the parameter space where the posterior measure might center or concentrate at. In the finite-dimensional setting, the MAP point is defined as the point where the probability density function (pdf) of the posterior distribution attains the maximum value. However, this definition does not directly generalize to the infinite-dimensional setting as the pdf of the posterior measure does not exist or it is zero at every point. Nevertheless, we can define the MAP point ${m_{1}}\in X$ as the one where the ball $B(m_{c},\epsilon)=\{m\in X:||m_{c}-m||_{X}\leq\epsilon\}$, centered at $m_{c}$ with radius $\epsilon$, attains the maximum probability given by $\mu^{y}(B(m_{c},\epsilon))$ when $\epsilon\to 0$, i.e.

[TABLE]

Under the assumption of Gaussian prior measure $\mathcal{N}(m_{0},\mathcal{C}_{0})$ with $\mathcal{C}_{0}=\mathcal{A}^{-\alpha}$ defined in the last section, when the mean $m_{0}$ lies in the range of $\mathcal{C}_{0}^{-1/2}$, or equivalently in the domain of the differential operator with fractional power $\mathcal{A}^{\alpha/2}$, we have that ${m_{1}}$ also lies in this range. In fact, it is the Cameron–Martin space denoted as $E$, which is a Hilbert space equipped with the inner-product $\langle\cdot,\cdot\rangle_{\mathcal{C}_{0}}=\langle\mathcal{C}_{0}^{-1/2}\cdot,\mathcal{C}_{0}^{-1/2}\cdot\rangle$ and the induced norm $||m||^{2}_{\mathcal{C}_{0}}=\langle m,m\rangle_{\mathcal{C}_{0}}$. By the definition of the posterior measure in (4), we have that ${m_{1}}$ is the solution of the following minimization problem

[TABLE]

where the cost functional (a term corresponding to the regularized least-squares minimization problem in the deterministic inverse problems) $\mathcal{J}$ is defined as

[TABLE]

where the first term is the potential defined in (2) and the second term arises from the Gaussian prior assumption. Under the quadratic boundedness and the Lipschitz continuity of the map $\mathcal{G}$, see in [30], the optimization problem (16) is well-defined and admit a unique solution. In the case where the parameter-to-observable map $\mathcal{G}$ depends on the parameter through a PDE solution, which is of particular interest in many applications, we can solve the optimization problem by a Lagrangian variational approach, which bypasses the direct computation of the derivative of the map $\mathcal{G}$ with respect to the parameter $m$. We will explain this approach in Section 5 for specific PDE models.

3.2 The Hessian

The local curvature of $\mathcal{J}$ at the MAP point indicates the shape of the distribution of the parameter, which can be assessed by the Hessian of $\mathcal{J}$ evaluated at the MAP point, i.e.

[TABLE]

where we denote the Hessian of the misfit term–the potential (2)–evaluated at the MAP point as $H_{\text{misfit}}$, i.e. $H_{\text{misfit}}=D_{m}^{2}\Phi(m)|_{m={m_{1}}}$. When the parameter-to-observable map $\mathcal{G}$ is linear with respect to $m$, which we write as

[TABLE]

where $G$ is a bounded linear map, i.e. $G\in\mathcal{L}(X,Y)$, then the Hessian misfit becomes

[TABLE]

where $G^{*}:Y\to X$ is the adjoint of $G$. When the map $\mathcal{G}$ is nonlinear with respect to $m$, we can perform the linearization as an approximation of the $\mathcal{G}$ as

[TABLE]

where we have assumed that $\mathcal{G}$ is Fréchet differentiable with respect to $m$ and by $\mathcal{G}^{\prime}$ we denote its Fréchet derivative. Under this linearization, the Hessian misfit can be obtained the same as in (20) with the linear map $G:=\mathcal{G}^{\prime}({m_{1}})$. The linear approximation (21) is reasonable when $\mathcal{G}$ is nearly linear in the parameter space or at least in the region that the posterior measure is non-negligible. For instance when the observation noise is very small, or size of observation data is very big, the posterior measure decays very fast to zero for the parameter away from the MAP point [12]. In the case that the nonlinear map $\mathcal{G}$ depends on the parameter through the PDE solution, we can also compute the Hessian misfit by a Lagrangian variational approach, which will be exemplified in Section 5.

3.3 Hessian-based parametrization

For a linear parameter-to-observable map $\mathcal{G}$, it is known (see [30]) that the covariance operator of the posterior measure $\mathcal{C}_{1}$ is given by the inverse of the Hessian at the MAP point, i.e.

[TABLE]

which is self-adjoint, positive definite and of a trace-class compact operator. Consequently, the posterior measure is exactly the Gaussian measure $\mu_{1}=\mathcal{N}({m_{1}},\mathcal{C}_{1})$, and the integration defined in (6) becomes

[TABLE]

Let $(\lambda^{1}_{j},\psi_{j}^{1})_{j\geq 1}$ denote the eigenpairs of the compact covariance operator $\mathcal{C}_{1}$, then the parameter $m$ can be expressed by the Karhunen–Loève expansion as

[TABLE]

Therefore, the integration defined in (23) now becomes a parametric integration

[TABLE]

We remark that the parametrization (24) is obtained by the Karhunen–Loève expansion of the posterior covariance operator at the MAP point. A re-parametric formulation based on a transformation of the parametrization w.r.t. the prior covariance operator is presented in Appendix A, which is equivalent to our development above in the abstract functional space. However, they are rather different from the computational point of view, see the comparison in Appendix A. We also present an efficient computation of the Hessian-based parametrization in Appendix B.

For a nearly linear map $\mathcal{G}$, we can also use this Gaussian measure as an approximation of the posterior measure and compute the integration (23) to approximate the true integration (6). For a nonlinear map $\mathcal{G}$, we can compute the true integration (6) by changing the measure as used in importance sampling [1]

[TABLE]

where the likelihood function $\pi^{y,1}_{\text{like}}(m)$ is given by

[TABLE]

where we have multiplied a rescaling constant $\exp\left(\mathcal{J}(m_{1})\right)$ and defined

[TABLE]

where the cost functional $\mathcal{J}$ is defined in (17). Note that, by construction, $\mathcal{J}_{1}(m)$ and its first and second derivative vanish for $m=m_{1}$, since $D_{m}\mathcal{J}|_{m=m_{1}}=0$ and $D_{mm}\mathcal{J}|_{m=m_{1}}={\mathcal{C}_{1}}^{-1}$. Under the parametrization (24), the integrals (26) become

[TABLE]

4 Adaptive Sparse Quadrature for Bayesian Inversion

In this section, we present a sparse quadrature method for the computation of the parametric integration (25) and (29) in the Bayesian inverse problems, following our work in [5]. We establish a dimension-independent convergence result of the sparse quadrature and present two algorithms for the construction of the sparse quadrature.

4.1 Sparse quadrature

For any integrable function $g:\Xi\to{\mathbb{R}}$, we define a sparse quadrature for the approximation of the integral ${\mathbb{E}}^{\mu}[g]=\int_{\Xi}g({\boldsymbol{\xi}})d\mu({\boldsymbol{\xi}})$ as

[TABLE]

Here, ${\boldsymbol{\nu}}=(\nu_{1},\nu_{2},\dots)\in\mathcal{F}$ is a multi-index in $\mathcal{F}=\{{\boldsymbol{\nu}}\in{\mathbb{N}}^{\infty}:|{\mathbb{J}}_{\boldsymbol{\nu}}|<\infty\}$, where ${\mathbb{J}}_{\boldsymbol{\nu}}$ is the support of ${\boldsymbol{\nu}}$, i.e., ${\mathbb{J}}_{\boldsymbol{\nu}}=\{j\in{\mathbb{N}}:\nu_{j}>0\}$. $\Lambda$ is an admissible index set satisfying

[TABLE]

where ${\boldsymbol{\mu}}\preceq{\boldsymbol{\nu}}$ means $\mu_{j}\leq\nu_{j}\text{ for all }j\geq 1$. $\triangle_{\boldsymbol{\nu}}$ is a tensor product difference quadrature operator defined by tensorization of a sequence of univariate difference operators

[TABLE]

where the univariate difference quadrature operator is defined as

[TABLE]

Here the univariate quadrature operator $\mathcal{Q}_{-1}(g)=0$ and $\mathcal{Q}_{\nu_{j}}(g)$ is given by

[TABLE]

where $(\xi_{j}^{(k_{j})},w_{j}^{(k_{j})})$, $k_{j}=0,\dots,m_{\nu_{j}}-1$, are the quadrature points and weights of certain quadrature rules, such as Gauss–Hermite and Genz–Keister, see [5] for details.

4.2 Dimension-independent convergence rates

In this section, we analyze the convergence of the quadrature error $|{\mathbb{E}}^{\mu}[g]-\mathcal{Q}_{\Lambda_{N}}(g)|$ with respect to $N$, the cardinality of the admissible index set $\Lambda_{N}$, i.e., $N=|\Lambda_{N}|$, for $g$ in the context of the Bayesian inverse problems.

We first present a regularity assumption of the parametric function $g:\Xi\to{\mathbb{R}}$ with respect to the parameter ${\boldsymbol{\xi}}$, which is [5, Assumption 2].

Assumption 1**.**

Let $0<q<1$, and $(\tau_{j})_{j\geq 1}$ be a positive sequence such that

[TABLE]

Let $r$ be the smallest integer such that $r>14/q$, we assume that the parametric function $g:\Xi\to{\mathbb{R}}$ satisfies

[TABLE]

where $|{\boldsymbol{\mu}}|_{\infty}=\max_{j\geq 1}\mu_{j}$, $\boldsymbol{\tau}^{2{\boldsymbol{\mu}}}=\prod_{j\geq 1}\tau_{j}^{2\mu_{j}}$, ${\boldsymbol{\mu}}!=\prod_{j\geq 1}\mu_{j}!$, and $\partial_{\boldsymbol{\xi}}^{\boldsymbol{\mu}}g=\prod_{j\geq 1}\partial_{\xi_{j}}^{\mu_{j}}g$.

The following theorem from [5] provides sufficient condition for $g$ such that the sparse quadrature error converges with dimension-independent convergence rate. \textcolorblackBefore stating the theorem, we introduce a sequence of useful constants

[TABLE]

which appear in the equivalence of (36) and a weighted sum of the coefficients of the Hermite expansion of $g$, see [5, Proposition 3.3], or [2, Theorem 3.1].

Theorem 1**.**

[5, Theorem 3.5]** Under Assumption 1, there exists an admissible index set $\Lambda_{N}\subset\mathcal{F}$, with the $N$ indexes ${\boldsymbol{\nu}}$ corresponding to the smallest $b_{\boldsymbol{\nu}}$ given by (37), and there exists a constant $C$ independent of $N$, such that the error of the sparse quadrature (30) based on the Gauss–Hermite quadrature rule is bounded by

[TABLE]

Remark 4.1**.**

The convergence rate in Theorem 1 does not depend on the parameter dimensions but only on the sparsity parameter $q$ that controls the regularity of $g$. We remark that this convergence rate is not necessarily optimal. In fact, $s=1/q+1/2$ is observed in numerics [5]. Note also that as the number of quadrature points \textcolorblack$N_{p}$ in $\Lambda_{N}$ grows as $O(N^{2})$, see [14, Proposition 18], so that the convergence with respect to the number of quadrature points is deteriorated to $N_{p}^{-s/2}$. However, similar convergence rate $N_{p}^{-s}$ is observed in practice because the Gauss–Hermite quadrature $\mathcal{Q}_{\nu}$ with $m_{\nu}=\nu+1$ points is exact for polynomials of degree $2\nu+1$, \textcolorblacksee [5].

The dimension-independent convergence result in Theorem 1 applies once Assumption 1 can be verified, see [5] for examples of parametric functions and parametric PDEs which satisfy this assumption. Here, we verify Assumption 1 for the parametric QoI $Q$ integrated with respect to the parameter $m^{1}$ in (25), which is parametrized according to the posterior measure. For the prior measure $\mathcal{N}(m_{0},\mathcal{C}_{0})$, where the covariance operator is constructed as $\mathcal{C}_{0}=\mathcal{A}^{-\alpha}$, we assume that $\mathcal{A}$ is a ‘Laplacian-like’ operator satisfying the following conditions as in [30].

Assumption 2**.**

$\mathcal{A}$* is defined on a Hilbert space $X$ over $D\subset{\mathbb{R}}^{d}$ and fulfills*

1. A.1

$\mathcal{A}$* is positive definite, self-adjoint, and invertible with eigenpairs $\{(\lambda_{j},\psi_{j})\}_{j\geq 1}$, where the eigenfunctions form an orthonormal basis of $X$.* 2. A.2

$\lambda=(\lambda_{j})_{j\geq 1}$* grows asymptotically as $j^{2/d}$, or $\lambda^{-1}\in\ell^{q}({\mathbb{R}}^{\infty})$, $\forall q>d/2$.*

As the prior covariance $\mathcal{C}_{0}=\mathcal{A}^{-\alpha}$, the larger $\alpha$ is, the more sparse or anisotropic the parameter $m$ becomes with respect to different dimensions of the random variables ${\boldsymbol{\xi}}=(\xi_{j})_{j\geq 1}$ in the Karhunen–Loève expansion (80). Here the sparsity is indicated by the decay of the eigenvalues $\lambda^{0}=(\lambda_{j}^{0})_{j\geq 1}\in\ell^{q}({\mathbb{R}}^{\infty})$ of the prior covariance $\mathcal{C}_{0}$, where we have $q>d/(2\alpha)$ according to assumption A.2. We expect that this sparsity is inherited in the Gaussian posterior measure, i.e., the eigenvalues $\lambda^{1}=(\lambda^{1}_{j})_{j\geq 1}$ of the Gaussian posterior covariance $\mathcal{C}_{1}$ is such that $\lambda^{1}\in\ell^{q}({\mathbb{R}}^{\infty})$. This is established in the following lemma.

Lemma 2**.**

Let $(\lambda_{j}^{0})_{j\geq 1}$ and $(\lambda^{1}_{j})_{j\geq 1}$ denote the non-increasing sequence of eigenvalues of the Gaussian prior covariance operator $\mathcal{C}_{0}=\mathcal{A}^{-\alpha}$ with $\alpha>d/2$ in Assumption 2 and the Gaussian posterior covariance operator $\mathcal{C}_{1}$ in (22), respectively, if the Hessian misfit $H_{\text{misfit}}$ in (18) is positive semi-definite, we have

[TABLE]

Moreover, under Assumption 2, we have $\lambda^{1}=(\lambda_{j}^{1})_{j\geq 1}\in\ell^{q}({\mathbb{R}}^{\infty})$, $\forall q>d/(2\alpha)$.

Proof.

For a positive semi-definite Hessian misfit $H_{\text{misfit}}$, we have

[TABLE]

where $\lambda_{1}$ is the smallest eigenvalue of $\mathcal{A}$ in Assumption 2. Therefore, we have

[TABLE]

which implies, by definition of $\mathcal{C}_{1}$ in (22), that $0\leq\lambda_{j}^{1}\leq\lambda_{1}^{-\alpha}$. Moreover, we can write

[TABLE]

As the Hessian misfit $H_{\text{misfit}}$ is self-adjoint (by definition) and positive semi-definite (by assumption), so is $\mathcal{C}_{0}^{1/2}H_{\text{misfit}}\mathcal{C}_{0}^{1/2}$ by Assumption 2. Hence, we can write

[TABLE]

where $(\tilde{\lambda}_{j},\tilde{\psi}_{j})_{j\geq 1}$ are the eigenpairs of $\mathcal{C}_{0}^{1/2}H_{\text{misfit}}\mathcal{C}_{0}^{1/2}$ with $\tilde{\lambda}_{j}\geq 0$, $\forall j\geq 1$. A combination of (42) and (43) yields

[TABLE]

Since $\mathcal{C}_{0}$ and $\mathcal{C}_{0}^{1/2}\mathcal{D}\mathcal{C}_{0}^{1/2}$ are both compact and self-adjoint, so is $\mathcal{C}_{1}$. By Courant–Fischer (min-max) theorem [25], the eigenvalues of the prior covariance operator $\mathcal{C}_{0}$ and the Gaussian posterior covariance operator $\mathcal{C}_{1}$ can be written as

[TABLE]

so that by (44) we have

[TABLE]

where the minimum eigenvalue of $\mathcal{C}_{0}^{1/2}\mathcal{D}\mathcal{C}_{0}^{1/2}$ satisfies $\lambda_{\text{min}}\left(\mathcal{C}_{0}^{1/2}\mathcal{D}\mathcal{C}_{0}^{1/2}\right)\geq 0$ since both $\mathcal{C}_{0}$ in (12) and $\mathcal{D}$ in (44) are positive. Therefore, (39) is established. Under Assumption 2, we have $\lambda^{0}=(\lambda_{j}^{0})_{j\geq 1}\in\ell^{q}({\mathbb{R}}^{\infty})$, $\forall q>d/(2\alpha)$, so that by (39) we obtain $\lambda^{1}=(\lambda_{j}^{1})_{j\geq 1}\in\ell^{q}({\mathbb{R}}^{\infty})$, $\forall q>d/(2\alpha)$. ∎

Corollary 3**.**

For linear (and linearized) inverse problems, where the Hessian misfit is given by (20), which is positive semi-definite, (39) holds.

Based on the above analysis for the decay of the eigenvalues of the Gaussian posteriori covariance, Theorem 5 is established for the dimension-independent convergence of the sparse quadrature error for the integration of certain quantities of interest with respect to the posterior measure (25). Before stating this theorem, let us present a summability result from [10] in the following lemma.

Lemma 4**.**

[10, Theorem 7.2]** For any $p\leq 1$, the sequence $(\frac{|{\boldsymbol{\mu}}|_{1}!}{{\boldsymbol{\mu}}!}\boldsymbol{\eta}^{\boldsymbol{\mu}})_{{\boldsymbol{\mu}}\in\mathcal{F}}\in\ell^{p}(\mathcal{F})$ for a sequence $\boldsymbol{\eta}\in{\mathbb{R}}^{\infty}$, if and only if $|\boldsymbol{\eta}|_{1}<1$ and $\boldsymbol{\eta}\in\ell^{p}({\mathbb{N}})$.

Theorem 5**.**

Under Assumption 2, suppose a functional $Q:X\to{\mathbb{R}}$ satisfies

[TABLE]

with a constant $C<\infty$, where $D^{k}_{m}Q(m^{1}({\boldsymbol{\xi}}))(\psi^{1}_{j_{1}},\cdots,\psi^{1}_{j_{k}})$ is the $k$-th derivative of $Q$ with respect to $m$ acting in the eigendirections $(\psi^{1}_{j_{1}},\cdots,\psi^{1}_{j_{k}})$ of $\mathcal{C}_{1}$ in (22), $j_{i}\in{\mathbb{N}}$, $i=1,\dots,k$, and evaluated at $m^{1}({\boldsymbol{\xi}})$ in (24). Then there exists an admissible index set $\Lambda_{N}$ with the $N$ indexes ${\boldsymbol{\nu}}\in\Lambda_{N}$ corresponding to the smallest $b_{\boldsymbol{\nu}}$ in (37), and there exists a constant $C$ independent of $N$, such that the convergence of the sparse quadrature error is bounded by

[TABLE]

Proof.

By Theorem 1, to obtain the bound of the convergence rate (48) we only need to find a sequence $(\tau_{j})_{j\geq 1}$ such that the regularity assumption (36) can be verified for the integrand $g({\boldsymbol{\xi}})=Q(m^{1}({\boldsymbol{\xi}}))$. In fact, we have

[TABLE]

where we have denoted $(\psi^{1})^{\boldsymbol{\mu}}=(\psi_{j_{1}},\cdots,\psi_{j_{k}})$ with $j_{i}\in{\mathbb{J}}_{\boldsymbol{\mu}}$ exhausting $|{\boldsymbol{\mu}}|_{1}$ indexes in the support of ${\boldsymbol{\mu}}$, $i=1,\dots,k$ with $k=|{\boldsymbol{\mu}}|_{1}$, and $\boldsymbol{\eta}=(\eta_{j})_{j\geq 1}$ with $\eta_{j}=\tau_{j}^{2}\lambda_{j}^{1}$. In the equality, we used the chain rule for the derivative $\partial_{\xi_{j}}Q(m^{1}({\boldsymbol{\xi}}))=D_{m}Q(m^{1}({\boldsymbol{\xi}}))\partial_{\xi_{j}}m^{1}({\boldsymbol{\xi}})$ and $\partial_{\xi_{j}}m^{1}({\boldsymbol{\xi}})=(\lambda_{j}^{1})^{1/2}\psi_{j}^{1}$. The first inequality is due to the assumption (47) and the second inequality is due to $\{|{\boldsymbol{\mu}}|_{\infty}\leq r\}\subset\mathcal{F}$, see [2, Theorem 4.1]. By Lemma 4 with $p=1$, we only need to find $\boldsymbol{\tau}$ such that $|\boldsymbol{\eta}|_{1}<1$. By Assumption 2 and Lemma 2, we set the sequence $(\tau_{j})_{j\geq 1}$ as $\tau_{j}\sim j^{\beta}$ with $2\beta-2\alpha/d<-1$, i.e., $\beta<\alpha/d-1/2$, such that $\sum_{j\geq 1}\tau_{j}^{2}\lambda_{j}^{1}<1$. Hence, we have $(\tau_{j}^{-1})_{j\geq 1}\in\ell^{q}({\mathbb{N}})$ for any $q>1/\beta$. Consequently, by Theorem 1 we obtain the convergence rate $N^{-s}$ for any $s=1/q-1/2<\alpha/d-1$. ∎

Remark 4.2**.**

In practice, the eigenfunctions $(\psi^{1}_{j})_{j\geq 1}$ of $\mathcal{C}_{1}$ are not given for the verification of the assumption (47). In such case, a stronger condition \textcolorblackI would keep stronger instead of using sufficient as (47) is already a sufficient condition for $Q$ is

[TABLE]

where $X^{\prime}$ is the dual space of $X$. (50) implies (47) by the orthonormality of $(\psi^{1}_{j})_{j\geq 1}$ in $X$. In the case $\max_{x\in D,j\geq 1}|\psi^{1}_{j}(x)|<C$ for some constant $C$, we can also assume

[TABLE]

then the convergence (48) is followed for $\tau_{j}\sim j^{\beta}$ such that $\sum_{j\geq 1}\tau_{j}^{2}\lambda_{j}^{1}<1/C^{2}$.

The result in the above theorem holds for general (nonlinear) Bayesian inverse problems as long as the integrand in (29) satisfies the assumption (47). It allows very quick growth, $O(k!)$, of the $k$-th derivative of $Q$ with respect to the parameter $m$. In particular, we obtain an explicit result for the integration (25) in the case of linear Bayesian inverse problems as follows.

Theorem 6**.**

For any bounded linear functional $f$ of the parameter $m^{1}$ (24), let $Q$ be its $k$-th power, $k=1,2,\dots$, i.e., $Q(m^{1})=f^{k}(m^{1})$, there exists an admissible index set $\Lambda_{N}$ with the $N$ indexes ${\boldsymbol{\nu}}\in\Lambda_{N}$ corresponding to the smallest $b_{\boldsymbol{\nu}}$ in (37), and there exists a constant $C$ independent of $N$, such that the convergence of the sparse quadrature error is bounded by

[TABLE]

Proof.

By Theorem 1, to obtain the bound of the convergence rate (52) we only need to find a sequence $(\tau_{j})_{j\geq 1}$ such that the regularity assumption (36) can be verified for the quantity of interest $Q$. In fact, for $Q(m^{1}({\boldsymbol{\xi}}))=f(m^{1}({\boldsymbol{\xi}}))$, we have

[TABLE]

where the first equality is due to the linearity of $f$ and the definition of $m^{1}$ in (24), the first inequality is due to the boundedness of $f$ and the orthonormality of $\psi_{j}^{1}$, i.e., $f^{2}(\psi_{j})\leq||f||^{2}_{X^{\prime}}||\psi_{j}^{1}||_{X}=||f||^{2}_{X^{\prime}}$. Therefore, the assumption (36) is verified if $\sum_{j\geq 1}\lambda_{j}^{1}<\infty$ and $\sum_{j\geq 1}\tau_{j}^{2}\lambda_{j}^{1}<\infty$. The former is ensured by Lemma 2. As for the latter, by Assumption 2 and Lemma 2, we find the sequence $(\tau_{j})_{j\geq 1}$ as $\tau_{j}\sim j^{\beta}$ with $2\beta-2\alpha/d<-1$, i.e., $\beta<\alpha/d-1/2$. Hence, we have $(\tau_{j}^{-1})_{j\geq 1}\in\ell^{q}({\mathbb{N}})$ for any $q>1/\beta$. Consequently, by Theorem 1 we obtain the convergence rate $N^{-s}$ for any $s=1/q-1/2<\alpha/d-1$.

For any $k>1$, with $g({\boldsymbol{\xi}})=f^{k}(m^{1}({\boldsymbol{\xi}}))$ we have

[TABLE]

which is bounded as long as ${\mathbb{E}}^{\mu}[f^{2i}]<\infty$, for all $i=0,\dots,k$, and $\beta<\alpha/d-1/2$ for $\tau_{j}\sim j^{\beta}$. ${\mathbb{E}}^{\mu}[f^{2i}]<\infty$ are verified for $i=0$ and $i=1$ in (53). Similarly, for any $i=2,\dots,k$, by setting $\xi_{0}=1$, $\lambda_{0}^{1}=1$, and $\psi_{0}^{1}=m_{1}$, we have

[TABLE]

where the equality is due to ${\mathbb{E}}^{\mu}[\xi_{j}^{n}]=0$ for any $j\geq 1$ and odd $n\in{\mathbb{N}}$, the first inequality is due to ${\mathbb{E}}^{\mu}[\xi^{2}_{l_{1}}\cdots\xi^{2}_{l_{i}}]\leq{\mathbb{E}}^{\mu}[\xi^{2i}]=(2i-1)!!$, $\forall l_{1},\dots,l_{i}\geq 0$, where $\xi\sim N(0,1)$, and the last inequality with $\sum_{l\geq 0}\lambda_{l}^{1}<\infty$ is due to Lemma 2. ∎

Remark 4.3**.**

Note that we need $\alpha>d$ in Theorem 5 and 6 to guarantee the convergence of the sparse quadrature and $\alpha>2d$ (so $s=\alpha/d-1\geq 1$) to ensure a convergence rate of $N_{p}^{-s/2}\leq N_{p}^{-1/2}$ with respect to the number of quadrature points $N_{p}$, see Remark 4.1. However, in practice, we observe faster convergence than predicted by (52) as it is only an upper bound, see computational evidence in [5], where the asymptotic convergence rate is observed to be $N_{p}^{-(s+1/2)}$ with respect to the number of quadrature points $N_{p}$, which implies a much weaker requirement of $\alpha>d$ to achieve the convergence as fast as $N_{p}^{-1/2}$, see the numerical tests in Sec 5.

4.3 Adaptive construction

We follow [5] to present two algorithms for the construction of the sparse quadrature, including a priori construction guarantees the dimension-independent convergence rate in Theorem 1, and a goal-oriented a posteriori construction based on a posteriori error indicator – the difference quadrature $\triangle_{\boldsymbol{\nu}}(g)$ in (32) that depends on each specific function $g$. Even though it can not guarantee the dimension-independent convergence rate in theory but turns out to be as accurate as the a priori construction in practice.

4.3.1 A priori construction

From Theorem 1 we observe that the dimension-independent convergence rate of the sparse quadrature can be achieved by choosing the admissible index set $\Lambda_{N}$ with indexes ${\boldsymbol{\nu}}\in\mathcal{F}$ corresponding to the largest value of $b_{{\boldsymbol{\nu}}}$. While we can compute $b_{\boldsymbol{\nu}}$ for all the indexes ${\boldsymbol{\nu}}\in\mathcal{F}_{r,J}$ where

[TABLE]

it is expensive/unfeasible if $r$ and $J$ are very large or infinite. Thanks to the monotonic increasing property of $b_{\boldsymbol{\nu}}$, see [5], we can adaptively construct the admissible index set $\Lambda_{N}$ by Algorithm 1 when $(\tau_{j})_{j\geq 1}$ is arranged in increasing order. We point out that this a priori construction does not depend on the function $g$ once the Assumption 1 can be verified for certain $q$, $\boldsymbol{\tau}$ and $r$. However, it is not always straightforward or possible to verify this assumption especially for nonlinear function with respect to the parameter, e.g. the integrand for the nonlinear inverse problems (29). Due to this difficulty, we adopt a goal-oriented a posteriori construction.

4.3.2 Goal-oriented a posteriori construction

We present a goal-oriented a posteriori construction of the sparse quadrature based on a dimension-adaptive tensor-product quadrature. It is initially developed in [15] by taking advantage of the different importance of different dimensions, or different regularity of $f$ with respect to different $y_{j}$, $j\in{\mathbb{J}}$, which we call adaptive sparse quadrature, whose associated grid $G_{\Lambda}$ is called adaptive sparse grid. The basic idea is based on the following adaptive process: given an admissible index set $\Lambda$, we search an index ${\boldsymbol{\nu}}\in\mathcal{F}$ among the forward neighbors of $\Lambda$ (${\boldsymbol{\nu}}\in\mathcal{F}$ is called a forward neighbor of $\Lambda$ if $\Lambda\cup{\boldsymbol{\nu}}$ is still admissible), at which $||\triangle_{\boldsymbol{\nu}}||_{\mathcal{S}}$ is maximized, and add this index to the index set $\Lambda=\Lambda\cup\{{\boldsymbol{\nu}}\}$. As the number of forward neighbors depends on the dimension $J$ (in fact, the forward neighbors of ${\boldsymbol{0}}$ are ${\boldsymbol{e}}_{j}$ for all $j\in{\mathbb{J}}$), in high or infinite dimensions, we can not search over all the forward neighbors. In such cases, it is usually reasonable to assume that the higher the dimensions, the less important they are, as determined e.g., by the fast decaying eigenvalues in Karhunen–Loève representation of the high/infinite dimensional random field. Therefore, we can explore the forward neighbors dimension by dimension in the set (see, e.g., [26, 9])

[TABLE]

where ${\mathbb{J}}_{\boldsymbol{\nu}}=\{j:\nu_{j}\neq 0\}$; $j(\Lambda)$ is the smallest $j$ such that $\nu_{j+1}=0$ for all ${\boldsymbol{\nu}}\in\Lambda$.

The adaptive sparse quadrature can be constructed following a basic greedy algorithm proposed in [15], which was improved on the data structure in [20] to copy with very high dimensions (e.g., up to $10^{4}$ dimensions in a personal laptop with $16$GB memory). We present the goal-oriented a posteriori construction in Algorithm 1. Note that for the a priori construction in this algorithm, we do not need to evaluate step 5 and 18 if the maximum number of indexes is imposed as the only stopping criterion. We can also replace the maximum number of indexes by the maximum number of points $|G_{\Lambda\cup\mathcal{N}(\Lambda)}|$. Moreover, for step 13 it is also a common practice to chose ${\boldsymbol{\nu}}$ as ${\boldsymbol{\nu}}=\operatornamewithlimits{argmax}_{{\boldsymbol{\mu}}\in\mathcal{N}(\Lambda_{N})}||\triangle_{{\boldsymbol{\mu}}}(f)||_{\mathcal{S}}/|G_{\boldsymbol{\mu}}|$ to balance the error and the work, e.g., [15, 23].

5 Numerical Experiments

In this section, we present both a linear and a nonlinear inverse problem to demonstrate the concentration of the posterior distribution and the convergence property of the sparse quadrature.

5.1 A linear inverse problem \textcolorblackwith analytic solution

5.1.1 Problem setup

We first consider an analytic linear inverse problem to explicitly demonstrate the variance reduction of the uncertain parameter from its prior distribution to its posterior distribution and to illustrate the dependence of the reduction on the prior distribution and the observation noise. In particular, we consider Poisson’s equation in a physical domain $D\subset{\mathbb{R}}^{d}\;(d=1,2,3)$ with homogeneous Dirichlet boundary condition: find $u\in H_{0}^{1}(D)$ such that

[TABLE]

where $m$ is a Gaussian random field with prior distribution $\mathcal{N}(m_{0},\mathcal{C}_{0})$, and $\triangle$ is the Laplace operator. We consider $\mathcal{C}_{0}=(-\beta\triangle)^{-\alpha}$ with the smoothness parameter $\alpha>d/2$, and the variance parameter $\beta>0$. By $\{(\lambda_{j},\psi_{j})\}_{j\geq 1}$ we denote the eigenpairs of $-\triangle$, so that the eigenpairs of $\mathcal{C}_{0}$ are given by $\{((\beta\lambda_{j})^{-\alpha},\psi_{j})\}_{j\geq 1}$, and we can write

[TABLE]

We assume that the observation data is given by

[TABLE]

where $\eta\sim\mathcal{N}(0,\sigma^{2}I)$, with $\sigma>0$ and $I$ being the identity operator in $L^{2}(D)$. Conditioned on the observation data, we seek the posterior distribution of $m$, which is an infinite-dimensional linear Bayesian inverse problem. Due to linearity, the posterior distribution is also Gaussian, i.e., $m\sim\mathcal{N}(m_{1},\mathcal{C}_{1})$ with

[TABLE]

and

[TABLE]

where the cost functional $\mathcal{J}(m)$ reads

[TABLE]

Note that $\mathcal{C}_{1}$ admits the eigendecomposition with the same eigenfunctions of $\mathcal{C}_{0}$

[TABLE]

Fig. 1 displays the decay of the square root of the eigenvalues of the prior and the posterior covariance operators with different noise level $\sigma=10^{-1},10^{-2}$, different prior variance parameter $\beta=10^{-1},10^{-2}$, and different smoothness parameter $\alpha=1,2$.

We can see that $\lambda_{j}^{1}<\lambda_{j}^{0}$ for all $j$, the uncertainty is reduced from the prior to the posterior in every dimension (eigendirection), and the reduction is more evident in the low dimensions, and less evident in high dimensions as $\lim_{j\to\infty}\lambda_{j}^{1}/\lambda_{j}^{0}=1.$ Moreover, we can see that the smaller $\sigma$ is (smaller noise), the larger the reduction becomes, which can be interpreted as smaller noise provides more informative observation data. Furthermore, the smaller $\beta$ is, i.e., the larger the prior variance, the larger the reduction becomes. Finally, as the smoothness parameter $\alpha$ becomes larger, $\lambda_{j}^{-\alpha-2}$ decays faster, so that the reduction effect is only evident in the first fewer dimensions. Note that the eigenvalues of the posterior covariance operator do not decay monotonically. To use the sparse quadrature that is adaptively constructed from low dimension to high dimension, we rearrange these eigenvalues in decreasing order.

5.1.2 Numerical test

We set the physical domain $D=(0,1)$ and use piecewise linear finite element in a uniform mesh with mesh width $h=2^{-10}$ for the discretization in physical space, which leads to 1023 parameter dimensions (note that $m(0)=m(1)=0$). The eigenpairs of the Laplace operator with Dirichlet boundary conditions in $D$ read $\lambda_{j}=\pi^{2}j^{2}$ and $\psi_{j}=\sin(j\pi x)$, $j\geq 1$. We write the Karhunen–Loève expansion of the parameter $m$ as follows

[TABLE]

where $i=0$ and $i=1$ represent the expansion w.r.t. the prior and the posterior, respectively. To generate the observational data $y$, we first take a random sample $m_{s}$ from the prior distribution $\mathcal{N}(m_{0},\mathcal{C}_{0})$ (with $m_{0}=0$, $\beta=5\times 10^{-2}$, and $\alpha=1$) and we compute the solution $u_{s}$ of the forward problem (58) using $m_{s}$ as forcing term. Then we set $y=u_{s}+\eta$, where $\eta\sim\mathcal{N}(0,\sigma^{2}I)$ and $\sigma=10^{-2}$.

In Fig. 2 we show one and two-dimensional anchored marginals from the prior and the posterior distribution (recall that due to linearity the posterior is also Gaussian). The anchored marginals are computed with respect to the first (i.e., dominant) 6 dimensions of the prior distribution. More explicitly, the one dimensional plots are obtained (up to a rescaling factor) by setting the $\xi_{j}$ in the expansion (65) to [math] for all indexes $j$ except index $i$ ($i=1,\ldots,6$) which is allowed to vary in ${\mathbb{R}}$; the two dimensional plots are obtained by setting all the $\xi_{j}$ in the expansion (65) to [math] expect the two indexes $k,l\in\{1,\dots,6\}$. Note how the posterior density is concentrated in a small region in the support of the prior density, especially for the first few dimensions. This concentration of the posterior may cause prior-based sparse quadrature to yield an inaccurate integration w.r.t. the posterior distribution; as shown in the middle-top part of the figure, Gauss–Hermite quadrature points w.r.t. the prior distribution may completely fail in detecting the posterior distribution. In contrast, quadrature points computed using the Hessian-based parametrization well capture the posterior distribution, as shown in the top-right part of the figure.

To demonstrate the convergence properties of the Hessian-based adaptive sparse quadrature, we consider two quantities of interest. The first is

[TABLE]

which is a nonlinear function of parameter $m$ evaluated at $x=0.5$, with all the derivatives bounded by

[TABLE]

which satisfies the assumption in Theorem 5. The second is

[TABLE]

which is the square of the first derivative of the solution $u$ at $x=0.5$. Note that $Q_{2}(m)$ \textcolorblackdepends quadratically on $m$ (as the derivative of the solution of the forward problem $u^{\prime}$ is bounded and linear with respect to $m$), align with Theorem 6. Note that the delta functional $\delta_{0.5}(m)=m(0.5)$ acting the eigenfunctions $\sin(\pi jx)$ is uniformly bounded by $1$. We use the approximation $u^{\prime}(0.5)\approx(u(0.5+h)-u(0.5-h))/(2h)$ for mesh size $h$. Both integrals ${\mathbb{E}}^{\mu^{y}}[Q_{1}]$ and ${\mathbb{E}}^{\mu^{y}}[Q_{2}]$ can be explicitly computed. More explicitly,

[TABLE]

where ${\boldsymbol{m}}_{1}$ is the coefficient vector of the MAP point $m_{1}$ in the finite element space, ${\boldsymbol{e}}_{0.5}=(0,\dots,0,1,0,\dots,0)^{\top}$ is a vector of the same size as ${\boldsymbol{m}}_{1}$ with one at $x=0.5$ and zero everywhere else, $\boldsymbol{\psi}_{j}$ is the coefficient vector of the eigenfunction. For $Q_{2}$,

[TABLE]

where ${\mathbb{C}}_{1}$ is the covariance matrix of ${\boldsymbol{m}}^{1}$, $\boldsymbol{d}_{0.5}=10\times(0,\dots,0,-1,0,1,0,\dots,0)^{\top}/(2h)$ corresponding to the derivative vector at $x=0.5$, ${\mathbb{A}}$ and ${\mathbb{M}}$ are the stiffness and mass matrix corresponding to the Laplace operator $-\triangle$ and the identity operator $I$.

We run the adaptive sparse quadrature by Algorithm 1 for both the prior-based (with parametrization $m^{0}\sim\mathcal{N}(m_{0},\mathcal{C}_{0})$) and the Hessian-based (with parametrization $m^{1}\sim\mathcal{N}(m_{1},\mathcal{C}_{1})$) integration formulae (9) and (25), with both a priori and a posteriori construction methods. In this example, the concentration of the posterior in a small region of the support of the prior distribution (see Fig. 2) caused the prior-based adaptive sparse quadrature algorithm to terminate prematurely without adequate sampling of the posterior distribution and providing an erroneous integration result. To compute the quadrature errors for the two QoI we use the values computed in (66) and (67) as their references. Fig. 3 displays the maximum sparse grid level in each dimension constructed by the prior-based and Hessian-based sparse quadrature for the integration of $Q_{1}$ with $\alpha=1$. We can observe that the prior-based sparse quadrature allocates most of the quadrature points in the first dimension due to the concentration of the posterior distribution in the domain of the prior distribution along the first dimension, while the Hessian-based sparse quadrature does not show such concentration. By increasing the number of quadrature points, more and more dimensions are activated. We can see that with $10^{5}$ points the first $617$ dimensions are activated by the Hessian-based sparse quadrature construction. We remark that \textcolorblack$Q_{1}$ does not depend on $j$ if $j$ is even since $\psi_{j}^{1}(0.5)=\sin(0.5\pi j)=0$, so only the first grid level (with one quadrature point) is activated in the even dimensions, see the right of Fig. 3.

In Fig. 4 we plot the convergence of the sparse quadrature errors as a function of the number of quadrature points (or PDE solves). We can see that the prior-based quadrature errors do not converge and remain large for both the two quantities of interest. In contrast, the Hessian-based quadrature errors by both the a priori and the a posteriori constructions converge with rate independent of the (active) parameter dimension. For $Q_{1}$, the a posteriori construction leads to smaller quadrature errors than the a priori construction, which is observed in most other cases, see [5]. However, the opposite is observed for $Q_{2}$ as the a priori construction activates more dimensions than the a posteriori construction for $Q_{2}$. Note that $Q_{2}$ is only quadratic with respect to $\xi_{j}$, $j\geq 1$, thus one level of Gauss–Hermite quadrature is exact for each dimension, so that more activated dimensions result in more accurate quadrature. The asymptotic convergence rate for $Q_{1}$ is $N_{p}^{-s}$ with $s=1/2$ for $\alpha=1$ and $s=3/2$ for $\alpha=2$, which is $s=\alpha/d-1/2$, larger than the bound $\alpha/d-1$ in Theorem 5 and 6. For $Q_{2}$ we observe $N_{p}^{-s}$ with $s=3/2$ for $\alpha=1$ and $s=5/2$ for $\alpha=2$. The faster convergence for $Q_{2}$ is because of the higher regularity of the derivative of the solution $u^{\prime}$ than that of the parameter $m$ itself in the sense of Assumption 1. This also indicates that the convergence property of the sparse quadrature does not only depend on the parametrization but also on the regularity of the quantity of interest. We also plot the decay of the averaged Monte Carlo (MC) quadrature errors of 100 trials for the Hessian-based integration, from which we can observe the convergence rate $N_{p}^{-1/2}$ in all cases, even for the integrand with high regularity, e.g., $\alpha=2$.

5.2 A nonlinear inverse problem

5.2.1 Problem setup

We consider a nonlinear inverse problem where the forward model is an elliptic equation describing Darcy flow in a porous medium $D\subset{\mathbb{R}}^{d}\;(d=1,2,3)$: find $u\in V:=H^{1}_{\Gamma_{D}}(D)$ such that

[TABLE]

with suitable Dirichlet boundary condition on the boundary $\Gamma_{D}$ and zero Neumann boundary condition on the rest of the boundary. We assume that the diffusion coefficient is a lognormal random field, i.e., the parameter $m$ obeys a Gaussian distribution with mean $m_{0}$ and covariance $\mathcal{C}_{0}$, i.e., $m\sim\mathcal{N}(m_{0},\mathcal{C}_{0})$. Moreover, we assume that measurements of $m$ are available at a few locations $x^{l}\in D$, $l=1,\dots,L$. The mean of the parameter $m_{0}$ is obtained as a solution of the optimization problem

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the inner product in $L^{2}(D)$, $\mathcal{A}=-\beta\triangle+\gamma I$, $m_{\rm true}$ is the true parameter field, $\kappa$ is a penalization parameter, and $\mathcal{M}$ is a measurement operator defined as

[TABLE]

with suitable radius $r>0$. The covariance of $m$ is defined as

[TABLE]

We specify the parameter-to-observable map $\mathcal{G}$ in (1) as

[TABLE]

where $\mathcal{B}=(b_{1},\dots,b_{K})^{T}$ is a vector of linear observation functionals $b_{k}\in V^{*}$, $k=1,\ldots,K$, where

[TABLE]

Therefore, $\mathcal{G}$ is nonlinear w.r.t. $m$ as $u$ is so, leading to a nonlinear inverse problem. We specify the covariance matrix $\Gamma_{\text{noise}}=\sigma^{2}I$ for the noise $\eta$ in (1), where $I$ is the identity matrix of size $K\times K$. To compute the MAP point, we apply the inexact Newton–conjugate gradient method [3], i.e., using Newton iteration method for the nonlinear optimization w.r.t. the parameter $m$ and conjugate gradient method for the solution of the incremental parameter $\hat{m}$ at each Newton step, for which we need to compute the gradient and Hessian of the cost functional $\mathcal{J}$ (defined in (17)) constrained by problem (68). We introduce the Lagrangian functional

[TABLE]

where $p$ is the Lagrangian multiplier. Then the gradient of $\mathcal{J}$ can be computed as

[TABLE]

where $u$ and $p$ are the solutions of the state and adjoint equations obtained by the first order variation of the Lagrangian functional

[TABLE]

and the Hessian of $\mathcal{J}$ evaluated at $\hat{m}\in X$ can be computed as

[TABLE]

where $\hat{u}$ and $\hat{p}$ are the solutions of the incremental state and adjoint equations obtained by the second order variation of the Lagrangian functional

[TABLE]

5.2.2 Numerical test

We consider the physical domain $D=(0,1)$ with Dirichlet boundary condition $u(0)=1$ and $u(1)=0$, and specify the parameters for the prior distribution as $\alpha=1,\beta=2,\gamma=1,\kappa=10^{3}$, and the measurement locations $x^{l}=(l-1)2^{-2}$ with $l=1,\dots,5$. For the observation noise we set $\sigma=5\times 10^{-2}$ at the observation locations $x_{obs}^{k}=(k-1)2^{-6}$ with $k=1,\dots 65$. We use piecewise linear finite element for the discretization at a uniform mesh of width $h=2^{-10}$, thus resulting in a $1025$-dimensional problem. We set the radius for the measurement and the observation operators $r=r_{obs}=h$.

To generate the prior mean $m_{0}$ and the observation vector $y$, we first sample the true parameter $m_{\rm true}$ from the Gaussian distribution $\mathcal{N}(0,\mathcal{A}^{-1})$. Then we let $m_{0}$ be the solution of the quadratic optimization problem (69), and we set $y=\mathcal{B}\,u_{\rm true}+\eta$, where $u_{\rm true}$ is the solution of the forward problem (68) with $m=m_{\rm true}$ and $\eta\sim\mathcal{N}({\bf 0},\Gamma_{\rm noise})$ is a random noise vector.

In Fig. 5, we plot the square root of the eigenvalues $(\sqrt{\lambda_{j}})_{j\geq 1}$ of the prior covariance $\mathcal{C}_{0}$ in (93), the Gaussian posterior covariance $\mathcal{C}_{1}$ in (96), as well as the prior-preconditioned Hessian misfit term $H_{\text{misfit}}$ in (94) for completeness. The former two converge with an asymptotic rate of $O(j^{-1})$ for the sparsity parameter $\alpha=1$. The latter decay very fast for the first few dimensions which are the most informed by the observation data. We can observe a reduction of the variance from the prior to the Gaussian approximate of the posterior in these dimensions. Note that, for this particular problem with $f=0$ in (68), the data do not inform the constant component of $m$; for this reason we do not observe reduction of the variance for the first dimension of the prior distribution whose corresponding eigenvector is nearly constant in space.

Fig. 6 displays the anchored marginal densities of the prior $\rho^{0}(m)=\exp(-\frac{1}{2}||m-m_{0}||^{2}_{\mathcal{C}_{0}})$, the posterior $\rho^{post}=\exp(-\frac{1}{2}||y-\mathcal{G}(m)||_{\Gamma_{noise}}^{2})\rho^{0}(m)$, and the Gaussian approximate of the posterior $\rho^{1}=\exp(-\frac{1}{2}||m-m_{1}||^{2}_{\mathcal{C}_{1}})$, as well as their mutual joint density in the first 6 dimension of the KL expansion of the parameter w.r.t. its prior distribution. Note that all the densities are divided (normalized) by their maximal value. The data are particularly informative in the second dimension, leading to a much smaller variance of the posterior distribution with respect to the prior distribution, as shown by how the posterior density concentrates in a small region. The prior-based quadrature points fail to capture the local shape of the posterior density, thus leading to inefficient sparse quadrature. On the contrary, the posterior and the Gaussian approximate of the posterior densities are close to each other in each of the 6 leading dimensions in the parametrization of the prior distribution. Moreover, as shown in the top-right part of Fig. 6, the Hessian-based quadrature points effectively capture the shape of the posterior distribution.

We run Algorithm 1 to compute the integration of the QoI

[TABLE]

using both the prior-based quadrature (7) and the Hessian-based quadrature (26). We need to compute the integrals with the two integrands $Q_{1}$ and $Q_{2}$ given by

[TABLE]

for the prior-based quadrature (7) and by

[TABLE]

for the Hessian-based quadrature (26). We remark that due to the complexity of the posterior density and the nonlinear dependence of the solution $u$ on the parameter $m$, we can not rigorously verify the assumption for the two integrands in Theorem 5. The decay of the adaptive sparse quadrature errors is shown in Fig. 7, where the reference values to compute the quadrature errors are taken as the adaptive sparse grid quadrature results with $10^{4}$ points for the prior-based quadrature and $10^{5}$ points for the Hessian-based quadrature, respectively. The error of prior-based quadrature rule does not decay and remains oscillating for both $Q_{1}$ and $Q_{2}$; on the contrary, the error of the Hessian-based quadrature decays with an asymptotic rate of $O(N_{p}^{-1})$ for $Q_{1}$ and $Q_{2}$. We remark that in Fig. 7 we only show the quadrature errors of the sparse quadrature obtained by the a posteriori construction in Algorithm 1, since the a priori construction based only on the coefficient $b_{\boldsymbol{\nu}}$ in (37) failed to effectively detect the complicated shape of two QoIs in the parameter space and resulted in too much oscillation of the quadrature errors. This numerical test suggests that a posteriori construction should be used for nonlinear inverse problems as the posterior distribution/density function could have rather complicated shape. Moreover, we point out that the Gaussian approximate of the posterior distribution is quite close to the true posterior distribution as shown in Fig. 6. If the Gaussian approximate is very poor, which occurs for informative observation data with very small noise or high non-linearity of the forward model, we caution that the Hessian-based allocation of the quadrature points may not capture the posterior distribution so that the sparse quadrature could fail to approximate the required integrals. How to tackle this difficulty is an ongoing research.

6 Conclusions

In this work, we proposed and analyzed a Hessian-based adaptive sparse quadrature to compute infinite-dimensional integrals with respect to the posterior distribution in the context of infinite-dimensional Bayesian inverse problems. Dimension-independent convergence of the proposed method was theoretically established for linear inverse problems and demonstrated numerically for both linear and nonlinear inverse problems. The convergence of the adaptive sparse quadrature is faster than that of the Monte Carlo or Markov-chain Monte Carlo method, especially for sparser problems, i.e., for problems in which the eigenvalues of the covariance operator decay faster and the parameter-to-quantity-of-interest map is smoother. We illustrated the concentration of the posterior distribution in a small region of the support of the prior distribution for both linear and nonlinear inverse problems. We showed that a prior-based parametrization quadrature may fail to capture the shape of the posterior distribution and result in inaccurate values of integrals with respect to the posterior distribution, while the Hessian-based parametrization effectively allocates quadrature points in regions of high density of the posterior distribution and leads to correct values of the integrals. Further investigation and application of our proposed method is ongoing for Bayesian inverse problems with very small observation noise and highly nonlinear forward problems. These lead to more complicated and localized posterior distributions that cannot be effectively captured by the Hessian-based transformation, but may be with a measure transport [22].

Appendix A A re-parametric formulation

In this section, we briefly review the re-parametric formulation presented in [27], where $m$ is first pre-parametrized according to the prior measure $\mu_{0}=\mathcal{N}(m_{0},\mathcal{C}_{0})$ and then re-parametrized using Hessian-based directions. We then compare our development in Section 3.3 with such re-parametric formulation both in terms of computational efficiency and accuracy for the solution of Bayesian inverse problems constrained by large-scale PDE models. To start, we denote the eigenpairs of the covariance operator $\mathcal{C}_{0}$ as $(\lambda_{j}^{0},\psi_{j}^{0})_{j\geq 1}$ and denote

[TABLE]

so that the parameter $m$ can be expressed by the Karhunen–Loève expansion as

[TABLE]

The prior measure of $m$ can be fully represented by that of ${\boldsymbol{\xi}}$, which is Gaussian $\mathcal{N}(0,I)$, where $0=(0,0,\dots,)$ and $I=\text{diag}(1,1,\dots,)$. The MAP point of the posterior measure is a parameter ${\boldsymbol{\xi}}_{1}\in\Xi$ that solves the following minimization problem

[TABLE]

where the ${\boldsymbol{\xi}}$-parametric cost functional is given explicitly by (17) and (80)

[TABLE]

Note that we have used the eigendecomposition $\mathcal{C}_{0}^{-1}=\Psi_{0}\Lambda_{0}^{-1}\Psi_{0}^{\top}$ in the second term. We remark that the solution ${\boldsymbol{\xi}}_{1}\in\Xi$ implies that $m_{1}\in E$, the Cameron–Martin space associated with the prior measure $\mathcal{N}(m_{0},\mathcal{C}_{0})$. The Hessian of the ${\boldsymbol{\xi}}$-parametric cost functional $\mathcal{J}$ evaluated at the MAP point is given by

[TABLE]

where Hessian of the misfit term is obtained by applying the chain rule for derivative

[TABLE]

Consequently, the covariance of the Gaussian measure $\mathcal{N}({{\boldsymbol{\xi}}_{1}},\mathcal{C}_{1}^{\boldsymbol{\xi}})$ is given by

[TABLE]

Since the covariance operator $\mathcal{C}_{1}$ can be decomposed as $\mathcal{C}_{1}=\Psi_{1}\Lambda_{1}\Psi_{1}^{\top}$, we have

[TABLE]

where by the compact and self-adjoint property of $\mathcal{C}_{1}^{\boldsymbol{\xi}}$, i.e. there exist the eigenpairs of $\mathcal{C}_{1}^{\boldsymbol{\xi}}$, denoted as $(\Lambda_{1}^{\boldsymbol{\xi}},\Psi_{1}^{\boldsymbol{\xi}})$, and there holds $\mathcal{C}_{1}^{\boldsymbol{\xi}}=\Psi_{1}^{\boldsymbol{\xi}}\Lambda_{1}^{\boldsymbol{\xi}}(\Psi_{1}^{\boldsymbol{\xi}})^{\top}$, so that we can identify $L=\Psi_{1}^{\boldsymbol{\xi}}(\Lambda_{1}^{\boldsymbol{\xi}})^{1/2}$, so that $L^{-1}=(\Lambda_{1}^{{\boldsymbol{\xi}}})^{-1/2}(\Psi_{1}^{\boldsymbol{\xi}})^{\top}$. We make the change of variables by shifting ${{\boldsymbol{\xi}}_{1}}$, rotating $(\Psi_{1}^{\boldsymbol{\xi}})^{\top}$ and rescaling $(\Lambda_{1}^{\boldsymbol{\xi}})^{-1/2}$ for ${\boldsymbol{\xi}}$, i.e.

[TABLE]

for which we have ${\mathbb{E}}^{\boldsymbol{\rho}}[\boldsymbol{\zeta}]=0$ and ${\mathbb{E}}^{\boldsymbol{\rho}}[\boldsymbol{\zeta}\boldsymbol{\zeta}^{\top}]$ = I, so that $\boldsymbol{\zeta}\in\mathcal{N}(0,I)$; moreover, we have

[TABLE]

Replacing this ${\boldsymbol{\xi}}$ as a function of $\boldsymbol{\zeta}$ in (80), we obtain

[TABLE]

where $m_{1}=m({\boldsymbol{\xi}}_{1})$. Denoting $m^{1}(\boldsymbol{\zeta})=m({\boldsymbol{\xi}}(\boldsymbol{\zeta}))$, we obtain the same parametrization of the parameter $m$ as in (24) derived in Section 3.3.

Our Hessian-based parametrization and the above re-parametric formulation are therefore equivalent in the idealized case of no truncation in the Karhunen–Loève expansion of $m$ with respect to the prior. However, in practice, only a finite number of eigenmodes in such expansion can be computed. In [27], the authors first perform a finite-dimensional truncation of the pre-parametric parameter (80), then solve the minimization problem for the MAP point (81) in the subspace spanned by the dominant eigenmodes of the prior, and finally compute the Hessian and its eigenmodes projected onto such subspace. This may lead to discretization bias with respect to the stochastic parameter, since directions that are well informed by the data may not be correctly captured by the subspace of the dominant eigenmodes of the prior. In contrast, our formulation does not involve truncation of the parameter according to the prior measure and, therefore, allows to capture all the directions — resolved by the spatial discretization of the parameter field — that are informed by the data.

Appendix B Computation of the parametrization

Let $X_{h}$ denote a finite-dimensional subspace of $X$ with $N_{h}$ basis functions $\{\phi_{h}^{n}\}_{n=1}^{N_{h}}$, for instance finite element space or spectral space, where $h$ represents the mesh size or the reciprocal of the spectral order, then we can approximate the parameter $m$ by

[TABLE]

where ${\boldsymbol{m}}=(m_{1},\dots,m_{N_{h}})^{T}$ is the vector of the coefficients in the basis function expansion. Let ${\mathbb{M}}\in{\mathbb{R}}^{N_{h}\times N_{h}}$ denote the mass matrix with ${\mathbb{M}}_{ij}=\langle\phi_{h}^{j},\phi_{h}^{i}\rangle,\;,i,j=1,\dots,N_{h}$. Let ${\mathbb{A}}\in{\mathbb{R}}^{N_{h}\times N_{h}}$ denote the matrix corresponding to the elliptic differential operator (14), given by

[TABLE]

Then we can compute the following inner products as

[TABLE]

where ${\mathbb{A}}_{\alpha}=({\mathbb{A}}{\mathbb{M}}^{-1})^{\alpha-1}{\mathbb{A}}$, $\alpha=1,2,\dots$. The discrete eigenpairs of the prior covariance operator $\mathcal{C}_{0}$ are obtained by solving the generalized eigenvalue problem

[TABLE]

for which we can employ a Lanczos method [4] or a randomized method [17] that involves $O(\alpha J)$ elliptic PDE solves for computing $J$ eigenpairs. Both methods get rid of direct assembling of the matrix ${\mathbb{M}}{\mathbb{A}}_{\alpha}^{-1}{\mathbb{M}}$, which could be very expensive for large-scale problems, i.e., $N_{h}$ large.

To compute the discrete eigenpairs of the posterior covariance operator $\mathcal{C}_{1}$ in (22), we first need to compute the MAP point as the solution of the optimization problem (16), for which we apply a Newton-conjugate gradient method [3, 31], where for the solution of the increment at each Newton step, we use conjugate gradient method. At the MAP point, evaluation of the Hessian misfit term $H_{\text{misfit}}$ in a given direction can be realized by solving one forward PDE and three additional linear PDEs corresponding to the adjoint, the incremental state and adjoint problems as illustrated in Section 5 for a specific model. Then the eigenpairs of the covariance operator $\mathcal{C}_{1}$ in (22) can be obtained by solving two generalized eigenvalue problems. By ${\mathbb{H}}_{\text{misfit}}$ we (formally) denote the matrix version of the Hessian misfit term $H_{\text{misfit}}$. We solve the first generalized eigenvalue problem with $\lambda_{1}\geq\lambda_{2}\geq\cdots$

[TABLE]

by Lanczos or randomized method, which involves $O(J)$ PDE solves for $J$ eigenpairs. Then by algebraic manipulation we can write the posterior covariance matrix as [29]

[TABLE]

where $D_{J}=\text{diag}(\lambda_{1}/(1+\lambda_{1}),\dots,\lambda_{J}/(1+\lambda_{J}))$, with $J$ determined according to certain criteria, e.g., $\lambda_{J}$ is small enough, and $\Psi_{J}=(\boldsymbol{\psi}_{1},\dots,\boldsymbol{\psi}_{J})$. The second generalized eigenvalue problem is also solved by Lanczos or randomized method

[TABLE]

where $(\lambda_{j}^{1},\boldsymbol{\psi}_{j}^{1})_{j\geq 1}$ are the eigenpairs corresponding to the covariance operator $\mathcal{C}_{1}$. We remark that by solving the two step generalized eigenvalue problems, we don’t need to assemble the Hessian matrix, which involves $O(N_{h})$ PDE solves, but only need to perform $O(J)$ PDE solves. Considerable computational saving can thus be achieved for small $J$ when the eigenvalues of the Hessian of the misfit term and the posterior covariance decay fast.