Non-intrusive optimal experimental design for large-scale nonlinear Bayesian inverse problems using a Bayesian approximation error approach

In this article we are interested in parameters which take values in possibly infinite-dimensional Hilbert spaces. For a Hilbert space $\mathscr{H}$, the corresponding inner product is denoted by $\langle\cdot,\cdot\rangle_{\mathscr{H}}$ and the associated norm by $\left\|\cdot\right\|_{\mathscr{H}}$. For $\mathscr{H}_{1}$ and $\mathscr{H}_{2}$ Hilbert spaces, we let $\mathcal{L}(\mathscr{H}_{1},\mathscr{H}_{2})$ denote the space of bounded linear operators from $\mathscr{H}_{1}$ to $\mathscr{H}_{2}$.

In this article we are particularly interested in Gaussian measures on Hilbert spaces. To this end, throughout the article we use $\mathcal{N}(v_{0},\mathcal{C}_{vv})$ to denote a Gaussian measure with mean $v_{0}$ and covariance operator $\mathcal{C}_{vv}$ . In the infinite dimensional case, the covariance operator is required to satisfy certain regularity assumptions to ensure the Bayesian inverse problem is well-defined [53]. As such, we assume $\mathcal{C}_{vv}$ is a strictly positive self-adjoint trace-class operator. We recall the Gaussian measure $\mathcal{N}(v_{0},\mathcal{C}_{vv})$ on $\mathscr{H}$ induces the Cameron-Martin space $\mathscr{E}\coloneqq\text{range}(\mathcal{C}_{vv}^{1/2})$, which is endowed with the inner product $\langle u,w\rangle_{\mathscr{E}}\coloneqq\langle\mathcal{C}_{vv}^{-1/2}u,\mathcal{C}_{vv}^{-1/2}w\rangle_{\mathscr{H}}$ for all $u,w\in\mathscr{E}$.

For $\mathscr{H}_{1}$ and $\mathscr{H}_{2}$ real Hilbert spaces, and for a linear operator $\mathcal{A}\in\mathcal{L}(\mathscr{H}_{1},\mathscr{H}_{2})$, we denote by $\mathcal{A}^{\ast}\in\mathcal{L}(\mathscr{H}_{2},\mathscr{H}_{1})$ the adjoint. Moreover, for $u\in\mathscr{H}_{1}$ and $w\in\mathscr{H}_{2}$, we define the (outer product) operator $u\otimes w$ by $(u\otimes w)v=\langle w,v\rangle_{\mathcal{H}_{2}}u$, for any $v\in\mathscr{H}_{2}$. The (cross-)covariance operator (between $u$ and $w$) is then defined

The discussion presented here is carried out in the possibly infinite-dimensional Hilbert space setting [53]. We consider the problem of inferring an unknown parameter $v\in\mathscr{H}$ given observations $\mathbf{d}\in\mathbb{R}^{d}$ that are related to $v$ as

where $\mathbf{e}\in\mathbb{R}^{d}$ denotes measurement noise. In our target applications, the (potentially) nonlinear parameter-to-observable map $\mathcal{G}:\mathscr{H}\rightarrow\mathbb{R}^{d}$ can be decomposed into the composition of two operators, $\mathcal{G}(\cdot)\coloneqq(\mathcal{V}\circ\mathcal{S})(\cdot)$, where the parameter-to-state map $\mathcal{S}:\mathscr{H}\rightarrow\mathcal{X}$ (for some suitably-chosen function space $\mathcal{X}$) requires solving a partial differential equation (PDE), and $\mathcal{V}:\mathcal{X}\rightarrow\mathbb{R}^{d}$ denotes a linear spatio-temporal observation operator.

In the Bayesian approach for parameter estimation, $v$ is treated as a random variable and is endowed with a prior probability law ($\mu_{v}$) that encompasses any knowledge we have about $v$ prior to data collection. Here we use a Gaussian prior $\mathcal{N}(v_{0},\mathcal{C}_{vv})$ with prior mean $v_{0}\in\mathscr{E}$ and a trace-class covariance operator $\mathcal{C}_{vv}$. In practice it is often assumed that the data is corrupted by Gaussian noise which is independent of the parameter, i.e., $\mathbf{e}\sim\mathcal{N}(\mathbf{e}_{0},\mathbf{\Gamma_{e}})$, where $\mathbf{\Gamma_{e}}\in\mathbb{R}^{d\times d}$ is a symmetric-positive-definite (SPD) covariance matrix and $\mathbf{e}\perp v$.

Within the Bayesian framework, the solution to the inverse problem is the (data-informed) posterior law of $v|\mathbf{d}$, which we denote by $\mu_{v|d}$. The posterior law is established as a Radon-Nikodym derivative via Bayes’ law and takes the following form under our additive Gaussian noise model,

Fully characterizing the posterior in the case of PDE-constrained inverse problems is generally infeasible. As a computationally practical alternative it is common to compute a Gaussian approximation to the posterior measure. A commonly utilized approximation is the so-called Laplace approximation [55], which approximates $\mu_{v|d}$ as a Gaussian, $\mathscr{L}_{\mu_{v|d}}\coloneqq\mathcal{N}(v_{\rm MAP},\mathcal{C}_{\rm post})$, centered around the maximum a posteriori (MAP) estimate [27],

The approximate posterior covariance operator of the Laplace approximation is given by $\mathcal{C}_{\rm post}=\mathcal{H}^{-1}_{v_{\rm MAP}}$, where $\mathcal{H}_{v_{\rm MAP}}$ denotes the Hessian of $\mathcal{J}$ evaluated at the MAP estimate.

Construction of the Laplace approximation presents many challenges. While the MAP estimator is guaranteed to exist under relatively mild assumptions on the PTO map and data misfit, uniqueness of the estimator can not be guaranteed in general for nonlinear Bayesian inverse problems (see [19] for details). Additional theoretical difficulties arise if the map $\mathcal{G}$ is not differentiable since the Laplace approximation requires access to derivatives of the PTO map. Constructing the Laplace approximation is also computationally costly. Finding the MAP estimator via minimization of (4) can involve many applications of the PTO map and its adjoint. Additionally, the inverse Hessian defining the posterior covariance operators typically not available in explicit form and needs to be estimated with additional forward and adjoint solves using, e.g., optimal low-rank approximations as described in [52]. This computational cost is exacerbated when one seeks optimal designs for nonlinear Bayesian inverse problems. As discussed in Section 3, OED in the nonlinear setting typically requires solving many Bayesian inverse problems. Thus, typical approaches (e.g., [4]) involve building many Gaussian approximations.

An alternative approach to obtaining a Gaussian approximation to the posterior distribution that circumvents some of the aforementioned challenges is to employ a linear surrogate to the PTO map. We leave the precise form of the linear (affine) surrogate arbitrary for now, and denote the linear map with $\mathcal{F}\in\mathcal{L}(\mathscr{H},\mathbb{R}^{d})$. Naturally, replacing the nonlinear map with some linear approximation leads to errors that need to be accounted for to avoid erroneous overconfidence and bias in the resulting Gaussian posterior. These model approximation errors can be incorporated using the Bayesian approximation error (BAE) approach [34, 36].

Before outlining the BAE approach to account for the resulting approximation errors, we recall the explicit results of the linear Gaussian case.

Neglecting the model errors and uncertainties induced by the use of a surrogate to the PTO map in general leads to overconfidence in biased estimates [36, 34]. To avoid this, we employ the BAE approach, which has been used for incorporating various types of model uncertainties into the Bayesian inverse problem as well as the OED problem (see, e.g., [6, 48, 26, 2, 2]). Here we outline the BAE approach as it pertains to our method.

As alluded to in the preceding section, the typical approach to OED within the Bayesian framework requires repeatedly solving the optimization problem (4) and constructing the (local) Laplace approximation^***An alternative (approximate) approach is proposed in [56] (see Section 3.4) where the Laplace approximation is carried out at samples from the prior, thus negating the need to compute any MAP estimates. However, this approximate approach is not well-suited for highly ill-posed inverse problems.. While using a surrogate (or reduced order model approximation) to the PTO map can reduce the computational costs associated with computing the MAP estimate, the cost for accurate Monte Carlo (MC) approximation of the optimality criterion can still be significant since it requires approximating the MAP estimate for many data samples.

In what follows, we let $\mathcal{F}:\mathscr{H}\to\mathbb{R}^{d}$ denote a surrogate model. The starting point for the BAE approach is to rewrite the accurate relationship between the data and parameters (2) as

where $\mathbb{R}^{d}\ni\boldsymbol{\boldsymbol{\varepsilon}}(v)\coloneqq\mathcal{G}(v)-\mathcal{F}(v)$ is the approximation error, and $\mathbb{R}^{d}\ni\boldsymbol{\eta}(v)\coloneqq\boldsymbol{\boldsymbol{\varepsilon}}(v)+\mathbf{e}$ is the total error. At this point, a conditional Gaussian approximation for the approximation error is made, i.e.,

with (where we denote $\boldsymbol{\varepsilon}_{0}\coloneqq\mathbb{E}\boldsymbol{\varepsilon}$)

where $\mathcal{C}_{v\varepsilon}$ and $\mathcal{C}_{\varepsilon v}$ defined as in (1). As a direct consequence of (7) we have

and we approximate Bayes’ law (3) by

In some applications of the BAE the further approximation $\mathcal{C}_{v\varepsilon}=0$ is employed, which is often referred to as the enhanced error model [6, 34].

A particularly straight-forward choice of surrogate model is a (affine) linear model, i.e., take $\mathcal{F}\in\mathcal{L}(\mathscr{H},\mathbb{R}^{d})$. The standard choice for such a surrogate model is the (generalized) derivative of $\mathcal{G}$ evaluated at some nominal point $v_{\ast}\in\mathscr{H}$, i.e., $\mathcal{F}v=\mathcal{G}(v_{\ast})+\mathsf{D}_{m}\mathcal{G}(v_{\ast})(v-v_{\ast})$. However, as long as the approximation error is taken into account using the BAE approach, the resulting (approximate) posterior is independent of the particular choice of linear surrogate, as stated in the following theorem from [49].

As an immediate consequence, using (5), we have the following corollary.

In various inverse problems there are additional uncertain model parameters which are not estimated [34]. This is usually because $a$) these parameters are not of interest, or $b$) estimation of these parameters is too costly (or impossible). We will refer to these additional (to the primary parameters of interest) parameters as auxiliary parameters, though in the literature the terms nuisance parameters, secondary parameters, and latent parameters are also common. Neglecting the uncertainty in the auxiliary parameters while inferring the primary parameters typically results in misleading estimates and significantly underestimated uncertainty. The same can also be said for the OED process; neglecting the uncertainty in the auxiliary parameters while carrying out OED generally yields sub-optimal designs [2].

The BAE approach can be applied as a means to account for uncertainty in auxiliary parameters during both the inference and OED stages. Specifically, let $v=(m,\xi)$, with $m\in\mathscr{M}$ the primary parameter $\xi\in\mathscr{X}$ the auxiliary parameter. Then assuming the accurate PTO model is $\mathcal{G}:\mathscr{M}\times\mathscr{X}\to\mathbb{R}^{d}$, we can rewrite the relationship between the data and the primary parameters as (cf. (6))

In keeping with the theme of the paper, we will be particularly interested in cases where the surrogate model is linear, i.e., $\mathcal{F}\in\mathcal{L}(\mathscr{M},\mathbb{R}^{d})$. For such cases the following corollary, outlining the independence of the approximate posterior to the choice of linear(-ized) model can be useful.

Here we focus on inverse problems where data is collected at a set of sensors. In this setting, the OED goal is finding an optimal set of locations for sensor deployment in some prescribed measurement domain $\mathscr{D}$. To this end, we fix a candidate set of locations $\{\mathbf{s}_{i}\}_{i=1}^{s}$ (with each $\mathbf{s}_{i}\in\mathscr{D}$) and define the design problem as that of finding an optimal subset of locations from this set^†††This discrete approach to sensor placement is common in OED literature [1], however a continuous formulation is also possible [47].. To distinguish between different designs or sensor arrangements, a binary weight $w_{i}\in\{0,1\}$ is assigned to each sensor location $\mathbf{s}_{i}$. If $w_{i}=1$, then data is collected using the sensor at location $\mathbf{s}_{i}$, whereas a weight of $0$ implies no measurement is conducted at $\mathbf{s}_{i}$. For general Bayesian inverse problems, the optimal arrangement of $k$ sensors is then defined through the optimal weight vector $\mathbf{w}^{*}\in\{0,1\}^{s}$, which minimizes an expected utility function, i.e.,

The utility function $U$ assesses the effectiveness of using any sensor combination (defined via the weight vector $\mathbf{w}\in\{0,1\}^{s}$) in solving the Bayesian inverse problem with measurement data $\mathbf{d}\in\mathbb{R}^{d}$. Since the posterior measure depends on the data, minimizing the expectation of the utility function ensures that the chosen sensor placement works well on average for all possible realizations of the data. The utility function is typically problem-specific and some common choices for infinite-dimensional Bayesian OED are described in [1]. For illustrative purposes, we focus on the A-optimality criterion, though our approach could be extended to other utility functions. The $k$-sensor A-optimal design minimizes the expected value of the average posterior pointwise variance. Letting $\mathcal{C}_{v|d,w}(\mathbf{w},\mathbf{d})$ denote the posterior covariance operator at a fixed design $\mathbf{w}$ and measurement data $\mathbf{d}$, the A-optimal design can be obtained by solving the minimization problem (17) with $U(\mathbf{w},\mathbf{d})\coloneqq\operatorname*{trace}(\mathcal{C}_{v|d,w}(\mathbf{w},\mathbf{d}))$.

For Bayesian inverse problems governed by nonlinear PTO maps or involving non-Gaussian prior measures, there is no closed-form expression for the data-dependent posterior covariance operator. Thus efficient techniques for approximating the expected utility are required. In the infinite-dimensional setting, the approximation is often carried out by combining a Gaussian approximation to the posterior measure with a sample average approximation (SAA) to the expectation, e.g., using a Laplace or Gauss-Newton approximation (see [4, 43]). As mentioned in Section 2.2, employing these techniques to approximate the SAA to the expected utility requires computing the MAP estimator for many data samples.

In this section, we outline our approach to optimal sensor placement for Bayesian inverse problems where the accurate (but computationally costly) PTO map has been replaced with some linear surrogate. As in Section 2.3, with $\mathcal{F}\in\mathcal{L}(\mathscr{H},\mathbb{R}^{d})$ we denote the surrogate map that approximates the accurate nonlinear operator $\mathcal{G}:\mathscr{H}\rightarrow\mathbb{R}^{d}$.

To incorporate model error into the design-dependent Bayesian inverse problem, we follow the BAE approach outlined in Section 2.3. The data measured at any subset of the candidate locations can be modeled by

where $\varepsilon(v,\mathbf{w})=\mathbf{W}\left(\mathcal{G}(v)-\mathcal{F}v\right)$, and the random variable $\boldsymbol{\eta}|v(\mathbf{w})\sim\mathcal{N}(\mathbf{W}\boldsymbol{\eta}_{0|v},\mathbf{W}\mathbf{\Gamma}_{\eta|v}\mathbf{W}^{T})$ is used to approximate the total error at the chosen sensor locations. The total error mean $\boldsymbol{\eta}_{0|v}$ and covariance $\mathbf{\Gamma}_{\eta|v}$ are defined in (2.3) and (7) and can be estimated using Monte Carlo sampling as described in the end of Section 2.3. Since the design matrix enters linearly into the model, the design-dependent approximate posterior is also independent of the particular choice of linearization. That is, we can extend the results of 1 and 1 to design-dependent Bayesian inverse problems. This is done in the following corollary.

An immediate consequence of the invariance of the A-optimality criterion to the specific linearization is that the uncertainty-aware A-optimal designs are also independent of the specific choice of surrogate map $\mathcal{F}$. Thus the trivial map, the zero operator $\mathsf{O}\colon v\mapsto\boldsymbol{0}$, leads to the same optimal sensor placements as a more complicated tailored surrogate when model error is incorporated using the BAE approach. This is the key insight that facilitates our greedy, derivative-free approach to optimal experimental design, which we describe in Section 3.3.

In this section, we outline a greedy procedure for computing A-optimal sensor placements for both design problems (19) and (20). Although we emphasize that our approach could be used with the zero map surrogate (and this is what we will use for the numerical examples in Section 4 and Section 5), for generality we present the procedure using an arbitrary linear surrogate map $\mathcal{F}$.

For the remainder of this section, let $\textbf{v}\in\mathbb{R}^{n}$, $\mathbf{G}\colon\mathbb{R}^{n}\rightarrow\mathbb{R}^{d}$, and $\mathbf{F}\in\mathcal{L}(\mathbb{R}^{n},\mathbb{R}^{d})$ denote the discretized parameter $v$, accurate PTO map $\mathcal{G}$, and linear surrogate PTO map $\mathcal{F}$, respectively. The discretized uncertainty-aware model is thus

where $\boldsymbol{\eta}|\textbf{v}\sim\mathcal{N}(\boldsymbol{\eta}_{0|v},\mathbf{\Gamma}_{\eta|v})$, with $\boldsymbol{\eta}_{0|v}=\mathbf{e}_{0}+\boldsymbol{\varepsilon}_{0}-\mathbf{C}_{\varepsilon v}\mathbf{C}_{vv}^{-1}\left(\textbf{v}-\textbf{v}_{0}\right)$ and $\mathbf{\Gamma}_{\eta|v}=\mathbf{\Gamma_{e}}+\mathbf{\Gamma}_{\varepsilon\varepsilon}-\mathbf{C}_{\varepsilon v}\mathbf{C}_{vv}^{-1}\mathbf{C}_{v\varepsilon}$. Here, the discretized statistics $\textbf{v}_{0},\,\boldsymbol{\varepsilon}_{0},\,\mathbf{C}_{vv},\,\mathbf{C}_{v\varepsilon}=\mathbf{C}^{\ast}_{\varepsilon v}$, and $\mathbf{\Gamma}_{\varepsilon\varepsilon}$ are computed via Monte Carlo (as in (14)) using the discretized operators $\mathbf{F}$ and $\mathbf{G}$. Letting $\widetilde{\mathbf{F}}=\mathbf{F}+\mathbf{C}_{\varepsilon v}\mathbf{C}_{vv}^{-1}$, the corresponding discretized posterior covariance operator is then

where $\widetilde{\mathbf{F}}^{\ast}$ denotes the adjoint of $\widetilde{\mathbf{F}}$, $\mathbf{Q}(\mathbf{W})\coloneqq\left[\mathbf{W}\left(\mathbf{\Gamma}_{\eta|v}+\widetilde{\mathbf{F}}\mathbf{C}_{vv}\widetilde{\mathbf{F}}^{\ast}\right)\mathbf{W}^{T}\right]^{-1}$ and the last equation follows from the Woodbury matrix identity and $\widetilde{\mathbf{F}}\mathbf{C}_{vv}=\mathbf{F}\mathbf{C}_{vv}+\mathbf{C}_{\varepsilon v}$.

Using the last equality and the cyclic property of the trace operator, we have

Thus, letting $\mathbf{K}(\mathbf{w})\coloneqq\mathbf{Q}(\mathbf{W})\mathbf{W}\left(\mathbf{F}\mathbf{C}_{vv}+\mathbf{C}_{\varepsilon v}\right)\left(\mathbf{C}_{vv}{\mathbf{F}}^{\ast}+\mathbf{C}_{v\varepsilon}\right)\mathbf{W}^{T}$ the A-optimal design satisfies

This reformulation of the A-optimality criterion is computationally advantageous if $n\gg d$, which is often the case for inverse problems, particularly when dealing with discretizations of infinite-dimensional parameters $v$.