Large-scale Bayesian optimal experimental design with derivative-informed projected neural network

Let $\mathcal{D}\subset\mathbb{R}^{n_{x}}$ denote a physical domain in dimension $n_{x}=1,2,3$. We consider the problem of inferring an uncertain parameter field $m$ defined in the physical domain $\mathcal{D}$ from noisy data $\mathbf{y}$ and a complex model represented by PDEs. Let $\mathbf{y}\in{\mathbb{R}}^{n_{y}}$ denote the noisy data vector of dimension $n_{y}\in{\mathbb{N}}$, given by

which is contaminated by the additive Gaussian noise $\boldsymbol{\varepsilon}\sim\mathcal{N}(\boldsymbol{0},\mathbf{\Gamma}_{\text{n}})$ with zero mean and covariance $\mathbf{\Gamma}_{\text{n}}\subset{\mathbb{R}}^{n_{y}\times n_{y}}$. Specifically, $\mathbf{y}$ is obtained from observation of the solution of the PDE model at $n_{y}$ sensor locations. $\mathcal{F}$ is the parameter-to-observable map which depends on the solution of the PDE and an observation operator that extracts the solution values at the $n_{y}$ locations.

We consider the above inverse problem in a Bayesian framework. First, we assume that $m$ lies in an infinite-dimensional real separable Hilbert space $\mathcal{M}$, e.g., $\mathcal{M}=\mathcal{L}^{2}(\mathcal{D})$ of square integrable functions defined in $\mathcal{D}$. Moreover, we assume that $m$ follows a Gaussian prior measure $\mu_{\text{pr}}=\mathcal{N}({m_{\text{pr}}},\mathcal{C}_{\text{pr}})$ with mean ${m_{\text{pr}}}\in\mathcal{M}$ and covariance operator $\mathcal{C}_{\text{pr}}$, a strictly positive, self-adjoint, and trace-class operator. As one example, we consider $\mathcal{C}_{\text{pr}}=\mathcal{A}^{-2}$, where $\mathcal{A}=-\gamma\Delta+\delta I$ is a Laplacian-like operator with prescribed homogeneous Neumann boundary condition, with Laplacian $\Delta$, identity $I$, and positive constants $\gamma,\delta>0$; see Stuart10 ; Bui-ThanhGhattasMartinEtAl13 ; PetraMartinStadlerEtAl14 for more details. Given the Gaussian observation noise, the likelihood of the data $\mathbf{y}$ for the parameter $m\in\mathcal{M}$ satisfies

where

is known as a potential function. By Bayes’ rule, the posterior measure, denoted as $\mu_{\text{post}}(m|\mathbf{y})$, is given by the Radon-Nikodym derivative as

where $\pi(\mathbf{y})$ is the so-called normalization constant or model evidence, given by

This expression is often computationally intractable because of the infinite-dimensional integral, which involves a (possibly large-scale) PDE solve for each realization $m$.

To measure the information gained from the data $\mathbf{y}$ in the inference of the parameter $m$, we consider a Kullback–Leibler (KL) divergence between the posterior and the prior, defined as

which is random since the data $\mathbf{y}$ is random. We consider a widely used optimality criterion, expected information gain (EIG), which is the KL divergence averaged over all realizations of the data, defined as

where the last equality follows Bayes’ rule (4) and the Fubini theorem under the assumption of proper integrability.

We consider the OED problem for optimally acquiring data to maximize the expected information gained in the parameter inference. The experimental of design seeks to choose $r$ sensor locations out of $d$ candidates $\{\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{d}\}$ represented by a design matrix ${\mathbf{W}}\in{\mathbb{R}}^{r\times d}\in\mathcal{W}$, namely, if the $i$-th sensor is placed at $\boldsymbol{x}_{j}$, then ${\mathbf{W}}_{ij}=1$, otherwise ${\mathbf{W}}_{ij}=0$:

Let ${\mathcal{F}}_{d}:\mathcal{M}\mapsto{\mathbb{R}}^{d}$ denote the parameter-to-observable map and $\boldsymbol{\varepsilon}_{d}\in{\mathbb{R}}^{d}$ denote the additive noise, both using all $d$ candidate sensors; then we have

Then the likelihood (2) for a specific design ${\mathbf{W}}$ is given by

and the normalization constant also depends on ${\mathbf{W}}$ as

From Section 2.2, we can see that the EIG $\Psi$ depends on the design matrix ${\mathbf{W}}$ through the likelihood function $\pi_{\text{like}}(\mathbf{y}|m,{\mathbf{W}})$. To this end, we formulate the OED problem to find an optimal design matrix ${\mathbf{W}}^{*}$ such that

with the ${\mathbf{W}}$-dependent EIG given by

To facilitate the presentation of our computational methods, we make a finite-dimensional approximation of the parameter field by using a finite element discretization. Let $\mathcal{M}_{n}\subset\mathcal{M}$ denote a subspace of $\mathcal{M}$ spanned by $n$ piecewise continuous Lagrange polynomial basis functions $\{\psi_{j}\}^{n}_{j=1}$ over a mesh with elements of size $h$. Then the discrete parameter $m_{h}\in\mathcal{M}_{n}$ is given by

The Bayesian inverse problem is then stated for the finite-dimensional coefficient vector $\mathbf{m}=(m_{1},\dots,m_{n})^{T}$ of $m_{h}$, with $n$ possibly very large. The prior distribution of the discrete parameter $\mathbf{m}$ is Gaussian $\mathcal{N}(\mathbf{m}_{\text{pr}},\mathbf{\Gamma}_{\text{pr}})$ with $\mathbf{m}_{\text{pr}}$ representing the coefficient vector of the discretized prior mean of ${m_{\text{pr}}}$, and $\mathbf{\Gamma}_{\text{pr}}$ representing the covariance matrix corresponding to $\mathcal{C}_{\text{pr}}=\mathcal{A}^{-2}$, given by

where $\mathbf{A}$ is the finite element matrix of the Laplacian-like operator $\mathcal{A}$, and $\mathbf{M}$ is the mass matrix. Moreover, let $\mathbf{F}_{d}:{\mathbb{R}}^{n}\to{\mathbb{R}}^{d}$ denote the discretized parameter-to-observable map corresponding to $\mathcal{F}_{d}$, we have $\mathbf{F}={\mathbf{W}}\mathbf{F}_{d}$ as in (9). Then the likelihood function corresponding to (10) for the discrete parameter $\mathbf{m}$ is given by

To solve the OED problem (12), we need to evaluate the EIG repeatedly for each given design ${\mathbf{W}}$. The double integrals in the EIG expression can be computed by a double-loop Monte Carlo (DLMC) estimator $\Psi^{dl}$ defined as

where $\mathbf{m}_{i}$, $i=1,\dots,n_{\text{out}}$, are i.i.d. samples from prior $\mathcal{N}(\mathbf{m}_{\text{pr}},\mathbf{\Gamma}_{\text{pr}})$ in the outer loop and $\mathbf{y}_{i}=\mathbf{F}(\mathbf{m}_{i})+\boldsymbol{\varepsilon}_{i}$ are the realizations of the data with i.i.d. noise $\boldsymbol{\varepsilon}_{i}\sim\mathcal{N}(\boldsymbol{0},\mathbf{\Gamma}_{\text{n}})$. $\hat{\pi}(\mathbf{y}_{i},{\mathbf{W}})$ is a Monte Carlo estimator of the normalization constant $\pi(\mathbf{y}_{i},{\mathbf{W}})$ with $n_{\text{in}}$ samples in the inner loop, given by

where ${\mathbf{m}}_{i,j}$, $j=1,\dots,n_{\text{in}}$, are i.i.d. samples from the prior $\mathcal{N}(\mathbf{m}_{\text{pr}},\mathbf{\Gamma}_{\text{pr}})$.

For complex posterior distributions, e.g., high-dimensional, locally supported, multi-modal, non-Gaussian, etc., evaluation of the normalization constant is often intractable, i.e., a prohibitively large number of samples $n_{\text{in}}$ is needed. As one particular example, when the posterior of $\mathbf{m}$ for data $\mathbf{y}_{i}$ generated at sample $\mathbf{m}_{i}$ concentrates in a very small region far away from the mean of the prior, the likelihood $\pi_{\text{like}}(\mathbf{y}_{i}|\mathbf{m}_{i,j},{\mathbf{W}})$ is extremely small for most samples $\mathbf{m}_{i,j}$, which which leads to a requirement of a large number of samples to evaluate $\hat{\pi}(\mathbf{y}_{i},{\mathbf{W}})$ with relatively small estimation error. This is usually prohibitive, since one evaluation of the parameter-to-observable map, and thus one solution of the large-scale PDE model, is required for each of $n_{\text{out}}\times n_{\text{in}}$ samples. This $n_{\text{out}}\times n_{\text{in}}$ PDE solves are required for each design matrix ${\mathbf{W}}$ at each optimization iteration.

Recent research has motivated the deployment of neural networks as surrogates for parametric PDE mappings BhattacharyaHosseiniKovachki2020 ; FrescaManzoni2022 ; KovachkiLiLiuEtAl2021 ; LiKovachkiAzizzadenesheliEtAl2020b ; LuJinKarniadakis2019 ; OLearyRoseberryChenVillaEtAl2022 ; NelsenStuart21 ; NguyenBui2021model ; OLearyRoseberryVillaChenEtAl2022 . These surrogates can be used to accelerate the computation of the EIG within OED problems. Specifically, to reduce the prohibitive computational cost, we build a surrogate for the parameter-to-observable map $\mathbf{F}_{d}:{\mathbb{R}}^{n}\mapsto{\mathbb{R}}^{d}$ at all candidate sensor locations by a derivative-informed projected neural network (DIPNet) OLearyRoseberry2020 ; OLearyRoseberryDuChaudhuriEtAl2021 ; OLearyRoseberryVillaChenEtAl2022 . Often, PDE-constrained high-dimensional parametric maps, such as the parameter-to-observable map $\mathbf{F}_{d}$, admit low-dimensional structure due to the correlation of the high-dimensional parameters, the regularizing property of the underlying PDE solution, and/or redundant information in the data from all candidate sensors. When this is the case, the DIPNet can exploit this low-dimensional structure and parametrize a parsimonious map between the most informed subspaces of the input parameter and the output observables. The dimensions of the input and output subspaces are referred to as the “information dimension” of the map, which is often significantly smaller than the parameter and data dimensions. The architectural strategy that we employ exploits compressibility of the map, by first reducing the input and output dimensionality via projection to informed reduced bases of the inputs and outputs. A neural network is then used to construct a low-dimensional nonlinear mapping between the two reduced bases. Error bounds for the effects of basis truncation, and parametrization by neural network are studied in BhattacharyaHosseiniKovachki2020 ; OLearyRoseberryDuChaudhuriEtAl2021 ; OLearyRoseberryVillaChenEtAl2022 .

For the input parameter dimension reduction, we use a vector generalization of an active subspace (AS) ZahmConstantinePrieurEtAl2020 , which is spanned by the generalized eigenvectors (input reduced basis) corresponding to the $r_{M}$ largest eigenvalues of the eigenvalue problem

where the eigenvectors $\mathbf{v}_{i}$ are ordered by the decreasing generalized eigenvalues $\lambda_{i}^{\text{AS}}$, $i=1,\dots,r_{M}$. For the output data dimension reduction, we use a proper orthogonal decomposition (POD)ManzoniNegriQuarteroni2016 ; QuarteroniManzoniNegri2015 , which uses the eigenvectors (output reduced basis) corresponding to the first $r_{F}$ eigenvalues of the expected observable outer product matrix,

where the eigenvectors $\mathbf{\phi}_{i}$ are ordered by the decreasing eigenvalues $\lambda_{i}^{\text{POD}}$, $i=1,\dots,r_{F}$. When the eigenvalues of AS and POD both decay quickly, the mapping $\mathbf{m}\mapsto\mathbf{F}_{d}(\mathbf{m})$ can be well approximated when $\mathbf{m}$ and $\mathbf{F}_{d}$ are projected to the corresponding subspaces with small $r_{M}$ and $r_{F}$; in this case approximation error bounds for reduced basis representation of the mapping are given by the trailing eigenvalues of the systems (19),(20). This allows one to detect appropriate “breadth” for the neural network via the direct computation of the associated eigenvalue problems, removing the need for ad-hoc neural network hyperparameter search for appropriate breadth. The neural network surrogate $\tilde{\mathbf{F}}_{d}$ of the map $\mathbf{F}_{d}$ then has the form

where $\mathbf{\Phi}_{r_{F}}\in{\mathbb{R}}^{d\times r_{F}}$ represents the POD reduced basis for the output, $\mathbf{V}_{r_{M}}\in{\mathbb{R}}^{n\times r_{M}}$ represents the AS reduced basis for the input, $\mathbf{f}_{r}\in{\mathbb{R}}^{r_{F}}$ is the neural network mapping between the two bases parametrized by weights $\mathbf{w}$ and bias $\mathbf{b}$. Since the reduced basis dimensions $r_{F},r_{M}$ are chosen based on spectral decay of the AS and POD operators, we can choose them to be the same; for convenience we denote the reduced basis dimension instead by $r$. The remaining difficulty is how to properly parametrize and train the neural network mapping. While the use of the reduced basis representation for the approximating map allows one to detect appropriate breadth for the neural network by avoiding complex neural network hyperparameter searches, and associated nonconvex neural network trainings, how to choose appropriate depth for the network is still an open question. While neural network approximation theory suggests deeper networks have richer representative capacities, in practice, for many architectures, adding depth eventually diminishes performance in what is known as the “peaking phenomenon”Hughes1968 . In general finding appropriate depth for e.g., fully-connected feedforward neural networks requires re-training from scratch different networks with differing depths. In order to avoid this issue we employ an adaptively constructed residual network (ResNet) neural network parametrization of the mapping between the two reduced bases. This adaptive construction procedure is motivated by recent approximation theory that conceives of ResNets as discretizations of sequentially minimizing control flows LiLinShen22 , where such maps are proven to be universal approximators of $L^{p}$ functions on compact sets. A schematic for our neural network architecture can be seen in Fig. 1.

This strategy adaptively constructs and trains a sequence of low-rank ResNet layers, where for convenience we take $r=r_{M}=r_{F}$ or otherwise employ a restriction or prolongation layer to enforce dimensional compatibility. The ResNet hidden neurons at layer $i+1$ have the form

with $z_{i+1},z_{i},b_{i}\in\mathbb{R}^{r}$, $w_{i,2},w_{i,1}^{T}\in\mathbb{R}^{r\times k}$, where the parameter $k<r$ is referred to as the layer rank, and it is chosen to be smaller than $r$ in order to impose a compressed representation of the ResNet latent space update (22). This choice is guided by the “well function” property in LiLinShen22 . The ResNet weights $\mathbf{w}=[(w_{i,2},w_{i,1},b_{i})]_{i=0}^{\text{depth}}$ consist of all of the coefficient arrays in each layer. Given appropriate reduced bases with dimension $r$, the ResNet mapping between the reduced bases is trained adaptively, one layer at a time, until over-fitting is detected in training validation metrics. When this is the case, a final global end-to-end training is employed using a stochastic Newton optimizer OLearyRoseberryAlgerGhattas2020 . This architectural strategy is able to achieve high generalizability for few and (input-output) dimension-independent data, for more information on this strategy, please see OLearyRoseberryDuChaudhuriEtAl2021 .

By removing the dependence of the input-to-output map on their high-dimensional and uninformed subspaces (complement to the low-dimensional and informed subspaces), we can construct a neural network of small input and output size that requires few training data. Since these architectures are able to achieve high generalization accuracy with limited training data for parametric PDE maps, they are especially well suited to scalable EIG approximation, since they can be efficiently queried many times at no cost in PDE solves, and require few high fidelity PDE solutions for their construction.

We propose to train a DIPNet surrogate $\tilde{\mathbf{F}}_{d}$ for the parameter-to-observable map ${\mathbf{F}}_{d}$, so that $\hat{\pi}(\mathbf{y}_{i},{\mathbf{W}})$ can be approximated as

where we omitted a constant $\frac{1}{(2\pi)^{n_{\text{out}}/2}\det(\mathbf{\Gamma}_{\text{n}})^{1/2}}$ since it appears in both the numerator and denominator of the argument of the log in the expression for the EIG. To this end, we can formulate the approximate EIG with the DIPNet surrogate as

where $\boldsymbol{\varepsilon}_{d,i}$ are i.i.d. observation noise. Thanks to the DIPNet surrogate, the EIG can be evaluated at negligible cost (relating to PDE solver cost) for each given ${\mathbf{W}}$, and does not require any PDE solves.

By the definition of the EIGs (17) and (24), we have

For sufficiently large $n_{\text{in}}$, we have that the normalization constants $\hat{\pi}(\mathbf{y}_{i},{\mathbf{W}})$ and $\tilde{\pi}(\mathbf{y}_{i},{\mathbf{W}})$ are bounded away from zero, i.e.,

for some constant $c_{i}>0$. Then we have

where we used the bound (29), which implies the bound (30) with constant

With the DIPNet surrogate, the evaluation of the DLMC EIG $\Psi^{nn}$ defined in (26) does not involve any PDE solves. Thus to solve the optimization problem

we can directly use a greedy algorithm that requires evaluations of the EIG, not its derivative w.r.t. ${\mathbf{W}}$. Let $S_{d}$ denote the set of all $d$ candidate sensors; we wish to choose $r$ sensors from $S_{d}$ that maximize the approximate EIG (approximated with the DIPNet surrogate). At the first step with $t=1$, we select the sensor $v^{*}\in S_{d}$ corresponding to the maximum approximate EIG and set $S^{*}=\{v^{*}\}$. Then at step $t=2,\dots,r$, with $t-1$ sensors selected in $S^{*}$, we choose the $t$-th sensor $v^{*}\in S_{d}\setminus S^{*}$ corresponding to the maximum approximate EIG evaluated with $t$ sensors $S^{*}\cup\{v^{*}\}$; see Algorithm 1 for the greedy optimization process, which can achieve (quasi-optimal) experimental designs with an approximation guarantee under suitable assumptions on the incremental information gain of an additional sensor; see JagalurMohanMarzouk21 and references therein. Note that at each step the approximate EIG can be evaluated in parallel for each sensor choice $S^{*}\cup\{v\}$ with $v\in S_{d}\setminus S^{*}$.

For the first numerical experiment, we consider an inverse wave scattering problem modeled by the Helmholtz equation with uncertain medium in the two-dimensional physical domain $\Omega=(0,3)^{2}$

where $u$ is the total wave field, $k\approx 4.55$ is the wavenumber, and $e^{2m}$ models the uncertainty of the medium, with the parameter $m$ a log-prefactor of the squared wavenumber. The right hand side $f$ is a point source located at $\mathbf{x}=(0.775,2.5)$. The perfectly matched layer (PML) boundary condition approximates a semi-infinite domain. The candidate sensor locations $\boldsymbol{x}_{i}$ are linearly spaced in the line segment between the edge points $(0.1,2.9)$ and $(2.9,2.9)$, with coordinates $\{(0.1+2i/35,2.9)\}_{i=0}^{49}$, as shown in Fig. 2. The prior distribution for the uncertain parameter $m$ is Gaussian $\mu_{\text{pr}}=\mathcal{N}({m_{\text{pr}}},\mathcal{C}_{\text{pr}})$ with zero mean ${m_{\text{pr}}}=0$ and covariance $\mathcal{C}_{\text{pr}}=(5.0I-\Delta)^{-2}$. The mesh used for this problem is uniform of size $128\times 128$. We use quadratic elements for the discretization of the wave field $u$ and linear elements for the parameter $m$, leading to a discrete parameter $\mathbf{m}$ of dimension $16,641$. The dimension of the wave field $u$ is $66049$, the wave is more than sufficiently resolved in regards to the Nyquist sampling criteria for wave problems. A sample of the parameter field $m$ and the PDE solution $u$ is shown in Fig. 2 with all $50$ candidate sensor locations marked in circles.

The network has $10$ low-rank residual layers, each with a layer rank of $10$. For this numerical example we demonstrate the effects of using different breadths in the neural network representation, in each case the ResNet learns a mapping from the first $r_{M}$ basis vectors for active subspace to the first $r_{M}$ basis vectors of POD. In the case that $r_{M}>50$ we use a linear restriction layer to reduce the ResNet latent representation to the $50$ dimensional output. For the majority of the numerical results, we employ a $r_{M}=50$ dimensional network. The neural network is trained adaptively using $4915$ training samples, and $1228$ validation samples. Using $512$ independent testing samples, the DIPNet surrogate was $81.56\%$ $\ell^{2}$ accurate measured as a percentage by

For more details on this neural network architecture and training, see OLearyRoseberryDuChaudhuriEtAl2021 . The computational cost of the $50$ dimensional active subspace projector using $128$ samples is equivalent to the cost of $842$ additional training data; since the problem is linear the additional linear adjoint computations are comparable to the costs of the training data generation. As we will see in the next example, when the PDE is nonlinear the additional active subspace computations are marginal. Thus for fair comparison with the same computational cost of PDE solves, we use $4915+842=5757$ samples for simple MC.

To test the efficiency and accuracy of our DIPNet surrogate, we first compute the log normalization constant $\log\pi(\mathbf{y})$ with our DIPNet surrogate for given observation data $\mathbf{y}$ generated by $\mathbf{y}={\mathbf{W}}\mathbf{F}_{d}(\mathbf{m})+\boldsymbol{\varepsilon}$. $\mathbf{m}$ is a random sample from the prior. We use in total $60000$ random samples for Monte Carlo (MC) to compute the normalization constant as the ground truth. Fig. 3 shows the $\log\pi(\mathbf{y})$ comparison for three random designs ${\mathbf{W}}$ that choose $15$ sensors out of $50$ candidates. We can see that DIPNet surrogate converges to a value close to the ground truth MC reference, while for the (simple) MC with $5757$ samples, the approximated value indicated by the green star is much less accurate than the DIPNet surrogate using $60000$ samples with similar cost in PDE solves. Note that the DIPNet training and evaluation cost for this small size neural network is negligible compared to PDE solves.

The left and middle figures in Fig. 4 illustrate the sample distributions of the relative approximation errors for the log normalization constant $\log\pi(\mathbf{y})$ and the EIG $\Psi$ ($n_{\text{out}}=200$) with (DIPNet MC with $60000$ samples) and without (simple MC with $5757$ samples) the DIPNet surrogate, compared to the true MC with $60000$ samples. These results show that using DIPNet surrogate we can approximate $\log\pi(\mathbf{y})$ and $\Psi$ much more accurately with less bias compared to the simple MC. The sample distributions of EIG at $200$ random designs compared to the design chosen by the greedy optimization using DIPNet surrogate for different number of sensors are shown in the right figure, from which we can see that our method can always chose better designs with larger EIG values than all the random designs.

Fig. 5 gives approximate values of the EIG with increasing number of outer loop samples $n_{\text{out}}$ using: the DIPNet surrogate with inner loop $n_{\text{in}}=60000$, simple DLMC with inner loop $n_{\text{in}}=5757$, and the truth computed with inner loop $n_{\text{in}}=60000$. We can see that the approximate values by the DIPNet surrogates are almost the same as the truth, while simple DLMC results are very inaccurate.

To show the effectiveness of truncated rank (breadth) for DIPNet surrogate, we evaluate the log normalization constant $\log\pi(\mathbf{y})$ and EIG $\Psi$ with breadth $=10,25,50,100$ and compared with true MC and simple MC in Fig. 6. We can see that with increasing breadth, the relative error is decreasing, but gets worse when breadth reaches $100$. With breadth $=100$, the difficulties of neural network training start to dominate and diminish the accuracy. We can also see it in the right part of the figure. The relative error of EIG approximation reduces (close to) linearly with respect to the generalization error of the DIPNet approximation of the observables with network breadth $=10,25,50$, which confirms the error analysis in Theorem 1. However, when the breadth increases to $100$, the neural network becomes less accurate (without using more training data), leading to less accurate EIG approximation.

For the second numerical experiment we consider an OED problem for an advection-diffusion-reaction equation with a cubic nonlinear reaction term. The uncertain parameter $m$ appears as a log-coefficient of the cubic nonlinear reaction term. The PDE is defined in a domain $\Omega=(0,1)^{2}$ as

Here $k=0.01$ is the diffusion coefficient. The volumetric forcing function f is a smoothed Gaussian bump located at $\boldsymbol{x}=(0.7,0.7)$,

The velocity field $\mathbf{v}$ is a solution of a steady-state Navier Stokes equation with shearing boundary conditions driving the flow (see the Appendix in OLearyRoseberryVillaChenEtAl2022 for more information on the flow). The candidate sensor locations are located in a linearly spaced mesh-grid of points in $(0.1,0.9)\times(0.1,0.9)$, with coordinates $\{(0.1i,0.1j),i,j=1,2,\dots,9\}$. The prior distribution for the uncertain parameter $m$ is a mean zero Gaussian with covariance $\mathcal{C}_{\text{pr}}=(I-0.1\Delta)^{-2}$. The mesh used for this problem is uniform of size $128\times 128$. We use linear elements for both $u$ and $m$, leading to a discrete parameter of dimension $16,641$. Fig. 7 gives a prior sample of the parameter field $m$ and the solution $u$ with all $100$ candidate sensor locations in white circles.

The neural network surrogate is trained adaptively using $409$ training samples and $102$ validation samples. Using $512$ independent testing samples the DIPNet network was $97.13\%\ell^{2}$ accurate (see equation 38). The network has $20$ low-rank residual layers, each with a layer rank of $10$, the breadth of the network is $r_{M}=r_{F}=25$. The computational cost of the $25$ dimensional active subspace projector using $256$ samples is equivalent to the cost of $34$ additional training data. As was noted before when the PDE is nonlinear the linear adjoint-based derivative computations become much less of a computational burden. Thus we use $409+34=443$ samples for simple MC for fair comparison.

We first examine the log normalization constant $\log\pi(\mathbf{y})$ computed with our $409$ PDE-solve-based DIPNet surrogate compared against the truth computed with $60000$ MC samples and the simple MC computed with $443$ samples. Fig. 8 shows the $\log\pi(\mathbf{y})$ comparison for three random designs that select $15$ sensors out of $100$ candidates. We can see that DIPNet MC converges to a value close to the true MC curves, while the simple MC’s green star computed with the same number of PDE solves ($443$) as DIPNet MC, has much worse accuracy.

The left figure of Fig. 9 shows the relative errors for $\log\pi(\mathbf{y})$ computed with the DIPNet surrogate and simple MC using $443$ samples based on the true MC with $60000$ samples, for $200$ random designs. We see again that DIPNet gives better accuracy with less bias than simple MC.

Fig. 10 shows the EIG approximations of three random designs with increasing number of outer loop samples $n_{\text{out}}$ using the DIPNet MC with $60000$ (inner loop) samples, simple MC with $443$ samples, and true MC with $60000$ samples. We can see that the values of DIPNet MC are quite close to the true MC, while simple MC is far off. Relative errors of the EIG $\Psi$ ($n_{\text{out}}=100$) computed with the DIPNet surrogate, the simple MC for $200$ random designs is given in the middle figure of Fig. 9. With the DIPNet DLMC $\Psi^{nn}$, we can use the greedy algorithm to find the optimal designs. The DIPNet greedy designs are presented as the pink crosses in the right figure of Fig. 9. We can see that the designs chosen by the greedy algorithm have much larger EIG values than all $200$ random designs.