Goal oriented optimal design of infinite-dimensional Bayesian inverse problems using quadratic approximations

The key components of a Bayesian inverse problem are the prior distribution, the data-likelihood, and the posterior distribution. The prior encodes our prior knowledge about the inversion parameter, which we denote by $m$. The likelihood, which incorporates the parameter-to-observable map, describes the conditional distribution of data for a given inversion parameter. Finally, the posterior is a distribution law for $m$ that is conditioned on the observed data and is consistent with the prior. These components are related via the Bayes formula Stuart10 . Here, we summarize the process for the case of linear Bayesian inverse problem.

The data likelihood. We consider a bounded linear parameter-to-observable map, $\mathcal{F}:\mathscr{M}\to\mathbb{R}^{d}$. In linear inverse problems governed by PDEs, we define $\mathcal{F}$ as the composition of a linear PDE solution operator $\mathcal{S}$ and a linear observation operator $\mathcal{B}$, which extracts solution values at a prespecified set of measurement points. Hence, $\mathcal{F}=\mathcal{B}\mathcal{S}$. In the present, work, we consider observation models of the form

We assume $m$ and $\boldsymbol{\eta}$ are independent, which implies, $\boldsymbol{y}|m\sim\mathsf{N}(\mathcal{F}m,\sigma^{2}\mathbf{I})$. This defines the data-likelihood.

Prior. Herein, $\mathscr{M}=L^{2}(\Omega)$, where $\Omega$ is a bounded domain in two- or three-space dimensions. This space is equipped with the $L^{2}(\Omega)$ inner product $\left\langle\cdot,\cdot\right\rangle$ and norm $\|\cdot\|=\left\langle\cdot,\cdot\right\rangle^{1/2}$. We consider a Gaussian prior law $\mu_{\text{pr}}:=\mathsf{N}(m_{\text{pr}},\mathcal{C}_{\text{pr}})$. To define the prior, we follow the approach in Stuart10 ; Bui-ThanhGhattasMartinEtAl13 . The prior mean is assumed to be a sufficiently regular element of $\mathscr{M}$ and the prior covariance operator $\mathcal{C}_{\text{pr}}$ is defined as the inverse of a differential operator. Specifically, let $\mathcal{E}$ be the mapping $s\mapsto m$, defined by the solution operator of

where $a_{1}$ and $a_{2}$ are positive constants. Then, the prior covariance is defined as $\mathcal{C}_{\text{pr}}:=\mathcal{E}^{2}$.

Posterior. For a Bayesian linear inverse problem with a Gaussian prior and a Gaussian noise model given by (2.1), it is well-known Stuart10 that the posterior is the Gaussian measure $\mu_{\text{post}}^{\boldsymbol{y}}\vcentcolon=\mathsf{N}\left(m_{\text{MAP}}^{\boldsymbol{y}},\mathcal{C}_{\text{post}}\right)$ with

where $\mathcal{F}^{*}$ denotes the adjoint of $\mathcal{F}$. Here, the posterior mean is the maximum a posteriori probability (MAP) point. Also, recall the variational characterization of this MAP point as the unique global minimizer of

in the Cameron–Martin space, $\mathrm{range}(\mathcal{C}_{\text{pr}}^{1/2})$; see DashtiStuart17 . The Cameron–Martin space plays a key role in the study of Gaussian measures on Hilbert spaces. In particular, this space is important in the theory of Bayesian inverse problems with Gaussian priors. Here, $\|\cdot\|_{\mathcal{C}_{\text{pr}}^{-1}}^{2}$ is the Cameron–Martin norm, $\|m\|_{\mathcal{C}_{\text{pr}}^{-1}}^{2}=\|\mathcal{C}_{\text{pr}}^{-1/2}m\|^{2}$.

It can be shown that the Hessian of $J$, denoted by $\mathcal{H}$, satisfies $\mathcal{H}=\mathcal{C}_{\text{post}}^{-1}$. In what follows, the Hessian of data-misfit term in (2.4) will be important. We denote this Hessian by $\mathcal{H}_{\text{mis}}\vcentcolon=\sigma^{-2}\mathcal{F}^{*}\mathcal{F}$. A closely related operator is the prior-preconditioned data-misfit Hessian,

which also plays a key role in the discussions that follow.

Lastly, we remark on the case when the forward operator is affine. This will be the case for inverse problems governed by linear PDEs with inhomogeneous source volume or boundary source terms. The model inverse problem considered in Section 5.2 is an example of such problems. In that case, the forward operator may be represented as the affine map $\mathcal{G}(m)=\mathcal{F}m+\boldsymbol{d}$, where $\mathcal{F}$ is a bounded linear transformation. Under the Gaussian assumption on the prior and noise, the posterior is a Gaussian with the same covariance operator as in (2.3) and with the mean given by $m_{\text{MAP}}^{\boldsymbol{y}}=\mathcal{C}_{\text{post}}\left(\sigma^{-2}\mathcal{F}^{*}(\boldsymbol{y}-d)+\mathcal{C}_{\text{pr}}^{-1}m_{\text{pr}}\right)$.

Discretization. We discretize the inverse problem using the continuous Galerkin finite element method. Consider a nodal finite element basis of compactly supported functions $\{\phi_{i}\}_{i=1}^{N}$. The discretized inversion parameter is represented as $m_{h}=\sum_{i=1}^{N}m_{i}\phi_{i}$. Following common practice, we identify $m_{h}$ with the vector of its finite element coefficients, $\boldsymbol{m}=[m_{1}\;m_{2}\;\cdots\;m_{N}]^{\top}$. The discretized inversion parameter space is thus $\mathbb{R}^{N}$ equipped with the mass-weighted inner product $\left\langle\boldsymbol{u},\boldsymbol{v}\right\rangle_{\mathbf{M}}\vcentcolon=\boldsymbol{u}^{\top}\mathbf{M}\boldsymbol{v}$. Here, $\mathbf{M}$ is the finite element mass matrix, $M_{ij}:=\int_{\Omega}\phi_{i},\phi_{j}\,d\boldsymbol{x}$, for $i,j\in\{1,\ldots,N\}$. Note that this mass-weighted inner product is the discretized $L^{2}(\Omega)$ inner product. Throughout the article, we use the notation $\mathbb{R}^{N}_{\mathbf{M}}$ for $\mathbb{R}^{N}$ equipped with the mass-weighted inner product $\left\langle\cdot,\cdot\right\rangle_{\mathbf{M}}$.

We use boldfaced symbols to represent the discretized versions of the operators appearing in the Bayesian inverse problem formulation. For details on obtaining such discretized operators, see Bui-ThanhGhattasMartinEtAl13 . The discretized solution, observation, and forward operators are denoted by $\mathbf{S}$, $\mathbf{B}$, and $\mathbf{F}$, respectively. Similarly, the discretized Hessian is presented as $\mathbf{H}$. We denote the discretized prior and posterior covariance operators by $\mathbf{\Gamma}_{\text{pr}}$ and $\mathbf{\Gamma}_{\text{post}}$, respectively. Note that $\mathbf{\Gamma}_{\text{pr}}$ and $\mathbf{\Gamma}_{\text{post}}$ are selfadjoint operators on $\mathbb{R}^{N}_{\mathbf{M}}$.

In the present work, an experimental design corresponds to an array of sensors selected from a set of candidate sensor locations, $\{x_{i}\}_{i=1}^{n}\subset\Omega$. In a classical OED problem, an experimental design is called optimal if it minimizes a notion of posterior uncertainty in the inversion parameter. This is different from a goal-oriented approach, where we seek designs that minimize the uncertainty in a goal quantity of interest.

To formulate an OED problem, it is helpful to parameterize sensor placements in some manner. A common approach is to assign weights to each sensor in the candidate sensor grid. That is, we assign a weight $w_{i}\geq 0$ to each $x_{i}$, $i\in\{1,\ldots,n\}$. This way, a sensor placement is identified with a vector $\boldsymbol{w}\in\mathbb{R}^{n}$. Each $w_{i}$ may be restricted to some subset of $\mathbb{R}$ depending on the optimization scheme. Here, we assume $w_{i}\in\{0,1\}$; a weight of zero means the corresponding sensor is inactive.

The vector $\boldsymbol{w}$ of the design weights is incorporated in the Bayesian inverse problem formulation through the data-likelihood Alexanderian21 . This yields a $\boldsymbol{w}$-dependent posterior measure. In particular, the posterior covariance operator is given by

There are several classical criteria in the OED literature. One example is the A-optimality criterion, which is defined as the trace of $\mathcal{C}_{\text{post}}(\boldsymbol{w})$. The corresponding discretized criterion is

The A-optimality criterion quantifies the average posterior variance of the inversion parameter field. To define a goal-oriented analogue of the A-optimality criterion, we need to consider the posterior variance of the goal-functional $\mathcal{Z}$ in (1.2). This is discussed in the next section.

In large-scale inverse problems, it is expensive to build the forward operator, the prior and posterior covariance operators, or the Hessian of the goal-functional, $\bar{\mathbf{H}}_{\text{z}}$ in (4.1). Therefore, matrix-free methods that only require applications of these operators on vectors are essential. A key challenge here is computation of the traces in (4.2). In this section, we present an approach for computing $\mathbf{\Psi}$ that relies on randomized trace estimation Avron2011 . As noted in AlexanderianPetraStadlerEtAl14 , the trace of a linear operator $\mathbf{T}$ on $\mathbb{R}^{N}_{\mathbf{M}}$ can be approximated via

This is known as a Monte Carlo trace estimator. The number $p$ of the required trace estimator vectors is problem-dependent. However, in practice, often a modest $p$ (in order of tens) is sufficient for the purpose of optimization.

We use Monte Carlo estimators to approximate the trace terms in (4.2). In particular, we use

where we have exploited the fact that $\mathbf{\Gamma}_{\text{pr}}$ and $\mathbf{\Gamma}_{\text{post}}$ are selfadjoint with respect to the mass-weighted inner product $\left\langle\cdot,\cdot\right\rangle_{\mathbf{M}}$. Thus, letting

we can estimate $\mathbf{\Psi}$ by

This enables a computationally tractable approach for approximating $\mathbf{\Psi}$. We outline the procedure for computing $\mathbf{\mathbf{\Psi}_{\text{rand,p}}}$ in Algorithm 1. The computational cost of this approach is discussed in Section 4.4.

The utility of methods based on Monte Carlo trace estimators in the context of OED for large-scale inverse problems has been demonstrated in previous studies such as HaberHoreshTenorio08 ; HaberMagnantLuceroEtAl12 ; AlexanderianPetraStadlerEtAl14 . A key advantage of the present approach is its simplicity. However, further accuracy and efficiency can be attained by exploiting the low-rank structures embedded in the inverse problem. This is discussed in the next section.

Here we present a structure-aware algorithm for estimating the $G_{q}$-optimality criterion that exploits low-rank components within the inverse problem. Namely, we leverage the often present low-rank structure in the (discretized) prior-preconditioned data misfit Hessian, $\tilde{\mathbf{H}}_{\text{mis}}\vcentcolon=\sigma^{-2}\mathbf{\Gamma}_{\text{pr}}^{1/2}\mathbf{F}^{*}\mathbf{F}\mathbf{\Gamma}_{\text{pr}}^{1/2}$.

Let us denote

Note that the posterior covariance operator can be represented as

As shown in Bui-ThanhGhattasMartinEtAl13 , we can obtain a computationally tractable approximation of $\mathbf{\Gamma}_{\text{post}}$ using a low-rank representation of $\tilde{\mathbf{H}}_{\text{mis}}$. Let $\{(\lambda_{i},\boldsymbol{v}_{i})\}_{i=1}^{k}$ be the dominant eigenpairs of $\tilde{\mathbf{H}}_{\text{mis}}$. We use

where $\mathbf{V}_{k}=[\boldsymbol{v}_{1}\;\boldsymbol{v}_{2}\;\cdots\boldsymbol{v}_{k}]$ and $\mathbf{\Lambda}_{k}=\text{diag}(\lambda_{1},\ldots,\lambda_{k})$. Now, define $\gamma_{i}\vcentcolon=\lambda_{i}/(\lambda_{i}+1)$ and $\mathbf{D}_{k}=\text{diag}(\gamma_{1},\gamma_{2},\cdots,\gamma_{k})$. We can approximate $\tilde{\mathbf{P}}$ using the Sherman-Woodbury-Morrison formula,

Substituting $\tilde{\mathbf{P}}$ by $\tilde{\mathbf{P}}_{\text{k}}$ in (4.6), yields the approximation

Subsequently, the $G_{q}$-optimality criterion (4.2) is approximated by

The following result provides a convenient expression for computing $\mathbf{\Psi}_{k}$.

Note that the second term in (4.10), $\mathrm{tr}(\tilde{\mathbf{H}}_{\text{z}}^{2})$, is a constant that does not depend on the experimental design (sensor placement). Therefore, when seeking to optimize $\mathbf{\Psi}_{k}$ as a function of $\boldsymbol{w}$, we can neglect that constant term and focus instead on minimizing the functional

The spectral approach for estimating the $G_{q}$-optimality criterion is outlined in Algorithm 2.

Note that the approximate posterior covariance operator $\mathbf{\Gamma}_{\text{post},\text{k}}$ can be used to estimate the classical A-optimality criterion as well. Namely, we can use $\mathrm{tr}(\mathbf{\Gamma}_{\text{post}})\approx\mathrm{tr}(\mathbf{\Gamma}_{\text{post},\text{k}})=\mathrm{tr}(\mathbf{\Gamma}_{\text{pr}})-\mathrm{tr}\left(\mathbf{\Gamma}_{\text{pr}}^{1/2}\mathbf{V}_{r}\mathbf{D}_{k}\mathbf{V}_{k}^{*}\mathbf{\Gamma}_{\text{pr}}^{1/2}\right)$. Since $\mathbf{\Gamma}_{\text{pr}}$ is independent of the experimental design, A-optimal designs can be obtained by minimizing

Furthermore, the present spectral approach can also be used for fast computation of the $G_{\ell}$-optimality criterion. In particular, it is straightforward to note that $G_{\ell}$-optimal designs can be computed by minimizing

This is accomplished by substituting $\mathbf{\Gamma}_{\text{post},\text{k}}$ into the discretized $G_{\ell}$-optimality criterion, given by (4.3), and performing some basic manipulations.

In this section, we present an algorithm for estimating the $G_{q}$-optimality criterion that relies on computing a low-rank SVD of prior-preconditioned forward operator. This approach relies on the specific form of the $\boldsymbol{w}$-dependent posterior covariance operator; see (2.7).

Before outlining our approach, we make the additional definitions

The following result enables a tractable representation for the $G_{q}$-optimality criterion.

Using Proposition 2, we can state the $G_{q}$-optimality criterion $\mathbf{\Psi}$ in (4.2) as

Note that the second term is independent of the design weights $\boldsymbol{w}$. Thus, we can ignore this term when minimizing $\mathbf{\Psi}$. In that case, we focus on

Computing the first term requires applications of $\mathbf{\Gamma}_{\text{post}}$ to vectors. We note,

This only requires a linear solve in the measurement space, when computing $\mathbf{D}_{\boldsymbol{w}}\tilde{\mathbf{F}}_{\boldsymbol{w}}\mathbf{\Gamma}_{\text{pr}}^{1/2}\boldsymbol{v}$. Once a low-rank SVD of $\tilde{\mathbf{F}}$ is available, this can be done without performing any PDE solves. The trace term in (4.15) can also be computed efficiently. First, we build

at the cost of $d$ applications of $\bar{\mathbf{H}}_{\text{z}}$ to vectors. The remaining part of computing $\hat{\mathbf{\Psi}}$ does not require any PDE solves. Let $\left\langle\cdot,\cdot\right\rangle_{2}$ denote the Euclidean inner product and $\{\boldsymbol{e}_{i}\}_{i=1}^{d}$ the standard basis in $\mathbb{R}^{d}$. We have

Computing this expression requires calculating $\mathbf{D}_{\boldsymbol{w}}\boldsymbol{e}_{i}$, for $i\in\{1,\ldots,d\}$, which amounts to building $\mathbf{D}_{\boldsymbol{w}}$. We are now in a position to present and an algorithm for computing the $G_{q}$-optimality criterion using a low-rank SVD of $\tilde{\mathbf{F}}$. This is summarized in Algorithm 3.

Here, we discuss the computational cost of the three algorithms presented above. We measure complexity in terms of applications of the operators $\mathbf{F}$ and $\mathbf{F}^{*}$, $\mathbf{\Gamma}_{\text{pr}}$, and $\bar{\mathbf{H}}_{\text{z}}$. Note that $\mathbf{F}$ and $\mathbf{F}^{*}$ correspond to forward and adjoint PDE solves. First we highlight the key computational considerations for each algorithm.

Computing $\mathbf{\Psi}_{\text{rand,p}}$:

The bottleneck in evaluating $\mathbf{\Psi}_{\text{rand,p}}$ is the need for $p$ applications of $\mathbf{\Gamma}_{\text{post}}$ and $2p$ applications of $\bar{\mathbf{H}}_{\text{z}}$. We assume that a Krylov iterative method is used to apply $\mathbf{\Gamma}_{\text{post}}$ to vectors, requiring $\mathcal{O}(r)$ iterations. In the present setting, $r$ is determined by the numerical rank of the prior-preconditioned data-misfit Hessian. Thus, each application of $\mathbf{\Gamma}_{\text{post}}$ requires $\mathcal{O}(r)$ forward and adjoint solves.

Computing $\mathbf{\Psi}_{\text{spec,k}}$:

In Algorithm 2, we need to compute the $k$ leading eigenpairs of $\tilde{\mathbf{H}}_{\text{mis}}$. In our implementation, we use the Lanczos method, costing $\mathcal{O}(k)$ applications of $\mathbf{F}$ and $\mathbf{F}^{*}$. Note also that Algorithm 2 requires $k+1$ applications of $\bar{\mathbf{H}}_{\text{z}}$ to vectors.

Computing $\mathbf{\Psi}_{\text{svd,r}}$:

This algorithm requires a low-rank SVD of $\tilde{\mathbf{F}}$ computed up-front. This can be done using a Krylov iterative method or randomized SVD HalkoMartinssonTropp11 . In this case, a rank $r$ approximation costs $\mathcal{O}(r)$ applications of $\mathbf{F}$ and $\mathbf{F}^{*}$. This algorithm also requires $d$ applications of $\bar{\mathbf{H}}_{\text{z}}$.

For readers’ convenience, we summarize the computational complexity of the methods in Table 1.

A few remarks are in order. Algorithm 2 and Algorithm 3 are more accurate than the randomized approach in Algorithm 1. When deciding between the spectral and low-rank SVD algorithms, several considerations must be accounted for. First, we note that the spectral approach is particularly cheap when the size of the desired design, i.e., the number of active sensors, is small. If the design vector $\boldsymbol{w}$ is such that $\|\boldsymbol{w}\|_{1}=k$, then $\text{rank}(\tilde{\mathbf{H}}_{\text{mis}})\leq k$. Thus, we only require the computation of the $k$ leading eigenvalues of $\tilde{\mathbf{H}}_{\text{mis}}$. The low-rank SVD approach is advantageous in the case when the forward model $\mathbf{F}$ is expensive and applications of $\bar{\mathbf{H}}_{\text{z}}$ are relatively cheap. This is due to the fact that no forward or adjoint solves are required in Algorithm 3, after precomputing the low-rank SVD of $\tilde{\mathbf{F}}$. However, the number of applications of $\bar{\mathbf{H}}_{\text{z}}$ is fixed at $d$, where $d$ is the number of candidate sensor locations. Lastly, we observe that all algorithms given in Section 4 may be modified to incorporate the low-rank approximation of $\bar{\mathbf{H}}_{\text{z}}$. Implementation of this is problem specific. Thus, methods presented are agnostic to the structure of $\bar{\mathbf{H}}_{\text{z}}$

Below, we first describe the model inverse problem under study and the goal-functional. Subsequently, we present our computational results.

In this section, we consider an example where the goal-functional depends nonlinearly on the inversion parameter.