# Scalable method for Bayesian experimental design without integrating over
posterior distribution

Vinh Hoang  Corresponding author<EMAIL_ADDRESS>Luis Espath
Faculty of Science, University of Nottingham, United Kingdom.  Sebastian
Krumscheid  Steinbuch Center for Computing and Institute for Applied and
Numerical Mathematics, Karlsruhe Institute of Technology, 76131 Karlsruhe,
Germany.  Raúl Tempone  Chair of Mathematics for Uncertainty Quantification,
RWTH Aachen University Germany.  Computer, Electrical and Mathematical
Sciences and Engineering, KAUST, and Alexander von Humboldt professor in
Mathematics of Uncertainty Quantification, RWTH Aachen University.

###### Abstract

We address the computational efficiency in solving the A-optimal Bayesian
design of experiments problems for which the observational map is based on
partial differential equations and, consequently, is computationally expensive
to evaluate. A-optimality is a widely used and easy-to-interpret criterion for
Bayesian experimental design. This criterion seeks the optimal experimental
design by minimizing the expected conditional variance, which is also known as
the expected posterior variance. This study presents a novel likelihood-free
approach to the A-optimal experimental design that does not require sampling
or integrating the Bayesian posterior distribution. Our proposed approach is
developed based on two properties of the conditional expectation: the law of
total variance and the property of orthogonal projection. The expected
conditional variance is obtained via the variance of the conditional
expectation using the law of total variance, and we take advantage of the
orthogonal projection property to approximate the conditional expectation. We
derive an asymptotic error estimation for the proposed estimator of the
expected conditional variance and show that the intractability of the
posterior distribution does not affect the performance of our approach. We use
an artificial neural network (ANN) to approximate the nonlinear conditional
expectation in the implementation of our method. We then extend our approach
for dealing with the case that the domain of experimental design parameters is
continuous by integrating the training process of the ANN into minimizing the
expected conditional variance. Specifically, we propose a nonlocal
approximation of the conditional expectation and apply transfer learning to
reduce the number of evaluations of the observation model. Through numerical
experiments, we demonstrate that our method greatly reduces the number of
observation model evaluations compared with widely used importance sampling-
based approaches. This reduction is crucial, considering the high
computational cost of the observational models. Code is available at
https://github.com/vinh-tr-hoang/DOEviaPACE.

###### Contents

  1. 1 Introduction
  2. 2 Background
  3. 3 Estimation of the tECV using PACE
    1. 3.1 PACE-based MC estimator of the tECV
    2. 3.2 Estimation of errors
    3. 3.3 Reduced variance estimator
  4. 4 Numerical experiment: Linear-Gaussian setting
    1. 4.1 One-dimensional case
    2. 4.2 High-dimensional case
  5. 5 Stochastic optimization for the continuous design domain
    1. 5.1 Nonlocal approximation of the conditional expectation
    2. 5.2 Algorithm
  6. 6 Numerical experiment: Electrical impedance tomography
    1. 6.1 Formulation and finite element model of EIT experiment
    2. 6.2 Numerical results
      1. 6.2.1 Analysis of error in estimating the tECV
      2. 6.2.2 Minimization of the tECV
  7. 7 Conclusion
  8. A Computing tECV using the IS approach
  9. B Proofs
    1. B.1 Proof of Proposition 1
    2. B.2 Proof of Proposition 2
    3. B.3 Orthogonality property of the CE
  10. C Linear approximation of the CE
    1. C.1 Analytical formulation of the linear approximation
    2. C.2 Empirical linear approximation of the CE
  11. D Reduced variance estimators of MSE and tECV
  12. E Gradient of the weighted MSE w.r.t. design parameters
  13. F  Derivation of the finite element formulation

##  1 Introduction

The design of experiments (DOE) aims to systematically plan experiments to
collect data and make accurate conclusions about a particular process or
system. This methodology has wide-ranging applications, such as in engineering
[1, 2], pharmaceutical [3, 4], and biological fields [5, 6]. Fundamentally,
ill-posed problems and uncertainties naturally arise in the DOE. The Bayesian
experimental design provides a general probabilistic framework for dealing
with these challenges [7, 8]. There are several optimality criteria available
in the Bayesian experimental design framework, for example, the A-optimal
criterion minimizes the expected conditional variance (ECV), whereas the
information gain criterion maximizes the expected Kullback–Leibler divergence
between the prior and posterior distributions. This study focuses on the first
criterion, _i.e._ , the A-optimal DOE. Compared with alternative optimality
criteria, the A-optimality criterion is of great practical appeal due to its
straightforward interpretation. Indeed, a reduced posterior variance indicates
a reduction in uncertainty.

The Bayesian experimental design generally requires estimating the expected
conditional statistical quantities, such as the ECV or the expected
information gain (EIG). Because this task is computationally expensive, many
previous works have sought to develop efficient methods for its execution.
Alexanderian et al. [9] proposed a scalable method for solving the A-optimal
DOE problem, that was specifically tailored to scenarios in which the
observational map is linear. For a nonlinear observational map, the posterior
distribution is typically intractable. Although the Markov chain Monte Carlo
(MCMC) algorithm and the approximate Bayesian computation approach can be used
to sample an intractable posterior distribution, these approaches are too
computationally intensive for the Bayesian DOE because the problem requires
sampling many different posterior distributions for each design candidate
[10]. A popular approach to alleviate the computational cost of the Bayesian
DOE for nonlinear observational maps is to use a Laplace approximation of the
posterior distribution [11, 12, 13, 14]. In [15], Alexanderian et al. used the
Laplace approximation approach in the context of the A-optimal DOE. In [16],
Beck et al. used this Laplace approximation approach for estimating the EIG,
which was later incorporated with the multilevel method in [17] and the
stochastic gradient descent to find the optimal design in [18]. For a
nonlinear observation map, the approaches mentioned above estimate expected
conditional quantities using a three-step process: first, they sample the
posterior distribution or its approximation for each observational sample;
second, they estimate the posterior characteristic, _e.g._ , the posterior
variance or the Kullback–Leibler divergences with respect to (_w.r.t._) the
prior density; and finally, they repeat the previous step for many
observational samples to estimate the expected value of the posterior
characteristics. The first two steps are equivalent to solving multiple
inverse problems, each of which is computationally expensive, particularly for
high-dimensional cases, because the posterior distribution becomes
intractable.

This study presents a novel method for estimating the ECV using the
projection-based approximation of the conditional expectation (PACE), designed
to handle computationally expensive observational maps. The relationship
between the uncertainty of the quantities of interest, repre-sented by random
variable (RV) $Q$, and their corresponding experimental observations, denoted
as RV $Y_{d}$, is assumed to be $Y_{d}=h(Q,d)+\varXi$, where $h$ is the
observational map, $\varXi$ represents measurement error, and $d$ denotes the
design parameters. Using the law of total variance, the ECV of the RV $Q$
given RV $Y_{d}$, denoted as
$\operatorname{E}\left[\operatorname{Var}\left[Q\mid Y_{d}\right]\right]$, is
obtained by approximating the conditional expectation (CE) of the RV $Q$ given
RV $Y_{d}$, denoted as $\operatorname{E}\left[Q\mid Y_{d}\right]$.
Specifically, the ECV can be computed as the difference between the variance
of $Q$ and the variance of CE $\operatorname{E}\left[Q\mid Y_{d}\right]$. To
approximate $\operatorname{E}\left[Q\mid Y_{d}\right]$, we utilize the
orthogonality of CE and introduce the map $\phi(Y_{d})$, which minimizes the
mean square error $\operatorname{E}\left[\lVert
Q-\phi(Y_{d})\rVert_{2}^{2}\right]$ under the assumption of the finite
variance of both $Q$ and $Y_{d}$. The proposed approach prevents the need to
sample or approximate the posterior distribution and eliminates the evaluation
of the likelihood function. Notably, we show that the relative mean absolute
error (MAE) of the PACE-based estimator of the ECV is of the order
$\mathcal{O}(N^{-2})$, where $N$ represents the number of evaluations of the
observational map. Moreover, the computational efficiency of the proposed
approach remains unaffected by the intractability of the posterior. In a
related work, Hoang et al. [19] used the PACE approach to develop a machine
learning-based ensemble conditional mean filter (ML-EnCMF) for nonlinear data
assimilation problems and demonstrated that the ML-EnCMF outperforms linear
filters in terms of accuracy.

To implement our method, we use an artificial neural network (ANN) [20, 21] to
approximate the CE, owing to ANN’s versatility and demonstrated effectiveness
in handling high-dimensional regression problems. To deal with continuous
design domains, we present a novel algorithm that effectively minimizes the
ECV using the stochastic gradient descent method. The training process of the
ANN is integrated into the algorithm employed to optimize the design
parameters, thereby improving computational efficiency. This integration is
possible because the objective functions for optimizing the ECV and training
the ANN are identical. Moreover, we propose a nonlocal approximation of the CE
and apply transfer learning to reduce the number of evaluations of the
observation model.

The remainder of the paper is structured as follows. In Sec. 2, we summarize
the background of the Bayesian experimental design. Section 3 details the PACE
framework used for the A-optimal DOE and its error estimation. In Section 4,
we illustrate our method and verify its error estimation using the linear-
Gaussian setting of the DOE. Then, in Section 5, we examine the continuous
design domain scenario and develop a stochastic optimization algorithm to
solve the A-optimal DOE. In Section 6, we apply our approach for the optimal
design of electrical impedance tomography experiments used to identify the
properties of a composite material. Finally, in Section 7, we conclude the
paper with a summary and perspectives.

##  2 Background

In Bayesian experimental design, the uncertainty associated with the
quantities of interest is modeled as an RV, denoted here as $Q$ and takes
value in $\mathbb{R}^{n}$. We consider a deterministic observational map,
denoted as $h$, that is parameterized by a vector of design parameters $d$ in
a domain $\mathcal{D}\subset\mathbb{R}^{\delta}$. This map transforms each
vector $q\in\mathbb{R}^{n}$ into a corresponding noise-free observational
vector in $\mathbb{R}^{m}$. We consider the scenario in which map $h$ consists
of two components: a computationally expensive numerical model denoted as
$h_{m}$, which solves the partial differential equation that governs the
experiments, and a measurement operator $h_{o}$. Typically, $h=h_{o}\circ
h_{m}$. Assuming that the measurement error is additive, we model the
observational RV $Y_{d}$ for a given vector $d$ of the experimental design
parameters as

${Y_{d}}(\omega)=h(Q(\omega),d)+\varXi(\omega),$ (1)

where $\varXi$ is the observational error RV, and $\omega$ denotes an outcome
in the sample space $\varOmega$ of the underlying probability space
$(\varOmega,\mathfrak{A},\mathbb{P})$. Further, subscript $d$ indicates the
dependence of the observational RV $Y_{d}$ on the design parameter vector $d$.
The inverse problem involves updating the prior distribution of the quantities
of interest, given the specific measurement data, $y$.

Using the Bayesian identification framework is a standard approach to solve
the inverse problem owing to the mathematically well-posedness of the
framework and its ability to handle uncertainty. Let us assume that the
distributions of the RVs $Q$, ${Y_{d}}$, and $\varXi$ have finite variances
and are absolutely continuous, _i.e._ , their densities exist. For a fixed
vector $d$, given the prior PDF $\pi_{Q}$ of RV $Q$ and the observational data
$y$, the Bayesian posterior PDF $\pi_{Q\mid{Y_{d}}}(\cdot\mid y)$ is given as
follows:

$\pi_{Q\mid{Y_{d}}}(q\mid
y)=\dfrac{\pi_{Q}(q)\,\pi_{\varXi}(y-h(q,d))}{\pi_{{Y_{d}}}(y)},$ (2)

where $\pi_{\varXi}$ is the density function of RV $\varXi$, and
$\pi_{{Y_{d}}}(y)$ is the evidence given by

$\pi_{{Y_{d}}}(y)=\int_{\mathbb{R}^{n}}\pi_{Q}(q)\,\pi_{\varXi}(y-h(q,d)){}\mathrm{d}q.$
(3)

The Bayesian experimental design searches for the parameter vector $d$ that
minimizes the posterior uncertainty. Here, we consider the A-optimal DOE,
which essentially minimizes the expected posterior variance. The posterior
variance can be formulated via the Bayes’ rule, Eq. (2), as

$\operatorname{Var}\left[Q\mid y;d\right]=\int_{\mathbb{R}^{n}}q^{\odot
2}\;\pi_{Q|{Y_{d}}}(q|y)\mathrm{d}q-\left[\int_{\mathbb{R}^{n}}q\;\pi_{Q|{Y_{d}}}(q|y)\mathrm{d}q\right]^{\odot
2},$ (4)

where ⊙ denotes the Hadamard square operator, e.g., ${q}^{\odot
2}=[q_{1}^{2},\dots,q_{n}^{2}]^{\top}$. Modeling the measurement data as an
RV, the distribution of the posterior variance is represented by the
conditional variance $\operatorname{Var}\left[Q\mid{Y_{d}}\right]$, which is
an $\mathbb{R}^{n}$-valued RV defined as

$\operatorname{Var}\left[Q\mid{Y_{d}}\right](\omega)\equiv\operatorname{Var}\left[Q\mid{Y_{d}}(\omega);d\right]\bigr{]}.$
(5)

The A-optimal DOE approach seeks the experimental design parameters that
minimize the ECV of the RV $Q$ given RV $Y_{d}$, denoted as
$\operatorname{E}\left[\operatorname{Var}\left[Q\mid{Y_{d}}\right]\right]$. In
cases where $n>1$, implying a multidimensional scenario, the A-optimal DOE
approach uses the element-wise sum of the ECV as the objective function, which
is referred to as the _total ECV_ (tECV) in this study. We denote
$Q=[Q_{1},...,Q_{n}]^{\top}$ where $Q_{i}$ is the $i$-th component of $Q$. The
tECV $V$ for a given vector $d\in\mathbb{R}^{\delta}$ of the experimental
design parameters is given by

$V(d)=\sum_{i=1}^{n}\operatorname{E}\left[\operatorname{Var}\left[Q_{i}\mid{Y_{d}}\right]\right],$
(6)

where $\operatorname{Var}\left[Q_{i}\mid{Y_{d}}\right]$ is the conditional
variance of the $i$-th component of the RV $Q$. Alternatively, the tECV can be
expressed as the trace of the expected conditional covariance, represented as
$V(d)\equiv\operatorname{tr}\bigr{(}\textnormal{Cov}\left[Q\mid
Y_{d}\right]\bigr{)}$, with $\operatorname{Cov}$ and $\operatorname{tr}$
denoting the covariance and trace operators, respectively.

We obtain the A-optimal DOE by minimizing the tECV, _i.e._ , its parameter
vector $d_{\text{A}}$ satisfies

$d_{\text{A}}=\;\arg\;\min_{d}\;V(d).$ (7)

By minimizing the tECV, the A-optimal DOE approach attempts to reduce the
overall variability of the posterior.

Solving the optimization problem in Eq. (7) requires an efficient and accurate
estimation of the tECV. A double-loop algorithm is needed for using the Monte
Carlo (MC) method to estimate the tECV based on Eq. (6). The outer loop
samples the RV $Y_{d}$ , and the inner loop estimates the posterior variance
for each sample of the outer loop. Let $\\{y^{(i)}\\}_{i=1}^{N_{\text{o}}}$ be
$N_{\text{o}}$ samples from the outer loop. Estimating the posterior variance
$\operatorname{Var}\left[Q\mid y^{(i)};\,d\right]$ for each sample $y^{(i)}$
involves sampling the posterior distribution $\pi_{Q|{Y_{d}}}(\cdot\mid
y^{(i)})$. However, the sampling of an intractable posterior distribution
poses a significant challenge, _i.e._ , the posterior can be concentrated in
considerably smaller regions compared to the prior, especially in high-
dimensional problems where $n,m\gg 1$. Using the importance sampling (IS)
method or the MCMC method to compute the tECV becomes computationally
expensive in such cases. The formulation of the IS double-loop estimator can
be found in Appendix A. Enhancing the efficiency of the IS estimator can be
done using the Laplace approximation. However, this method assumes that the
underlying distribution closely resembles a Gaussian distribution and requires
that the maximum a posteriori estimator be determined.

Notably, estimating the tECV involves applying the IS method to compute the
variances of $N_{\text{o}}$ different posteriors, which results in a
computational cost of the multiplicative order $N_{\text{o}}\times
N_{\text{i}}$, where $N_{\text{i}}$ is the number of samples in the inner
loop. Similarly, other approaches that estimate tECV by sampling the
posterior, such as the MCMC method, suffer from the same multiplicative growth
of computational cost. Herein, we present a novel approach to approximate the
tECV via the CE taking into account the orthogonal projection of the CE. The
PACE-based approach, as discussed in Section 3, does not require approximating
or sampling the posterior distribution.

##  3 Estimation of the tECV using PACE

We demonstrate that the tECV ($V(d)$) can be evaluated by approximating the CE
rather than sampling the conditional variance. For a given experimental
parameter vector $d$, the CE of RV $Q$ given RV ${Y_{d}}$, which is denoted as
$\operatorname{E}\left[Q\mid{Y_{d}}\right]$, satisfies

$\int_{A}Q(\omega)\mathrm{d}\mathbb{P}(\omega)=\int_{A}\operatorname{E}\left[Q\mid
Y_{d}\right](\omega)\mathrm{d}\mathbb{P}(\omega),$ (8)

for every measurable set $A$ in the $\sigma$-algebra generated by the RV
${Y_{d}}$. According to Doob–Dynkin lemma, the CE is a composition of the map
$\phi_{d}\;:\;\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ and the RV ${Y_{d}}$
which is given as

$\operatorname{E}\left[Q\mid{Y_{d}}\right](\omega)=\phi_{d}({Y_{d}}(\omega)).$
(9)

The map $\phi_{d}$ can be defined using the posterior density as

$\phi_{d}(y)=\int_{\mathbb{R}^{n}}q\;\pi_{Q|{Y_{d}}}(q,y)\mathrm{d}q.$ (10)

Alternatively, map $\phi_{d}$ can be obtained from the orthogonal projection
property. For the finite-variance RVs $Q$ and $Y_{d}$, the CE
$\operatorname{E}\left[Q\mid{Y_{d}}\right]$ is the $L_{2}$ orthogonal
projection of RV $Q$ onto the $\sigma$-algebra generated by the RV ${Y_{d}}$,
which is formulated as

$\operatorname{E}\left[Q\mid{Y_{d}}\right]=\phi_{d}({Y_{d}}),\quad\text{where}\quad\phi_{d}=\arg\min_{f\in\mathcal{S}(\mathbb{R}^{m},\mathbb{R}^{n})}\operatorname{E}\left[\bigl{\lVert}Q-f\circ{Y_{d}}\bigr{\rVert}_{2}^{2}\right],$
(11)

where $\lVert\cdot\rVert_{2}$ is the $L_{2}$ norm. Here,
$\mathcal{S}(\mathbb{R}^{m},\mathbb{R}^{n})$ represents the set of all
functions $f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ for which the variance
of the RV $f({Y_{d}})$ is finite. Appendix B.3 provides a proof of the
orthogonal projection property. Theoretical properties of the CE can be found
in [22] and in [23, Chapter 4]. The following proposition provides a
formulation for the tECV calculated using the CE
$\operatorname{E}\left[Q\mid{Y_{d}}\right]$.

######  Proposition 1.

The tECV $V(d)$ given in Eq. (6) can be expressed using CE
$\operatorname{E}\left[Q\mid{Y_{d}}\right]$ as:

$V(d)=\operatorname{E}\left[\left\lVert
Q-\operatorname{E}\left[Q\mid{Y_{d}}\right]\right\rVert_{2}^{2}\right].$ (12)

A proof of Proposition 1 is given in Appendix B.1. Using Proposition 1, the
tECV can be estimated by approximating CE
$\operatorname{E}\left[Q\mid{Y_{d}}\right]$.

###  3.1 PACE-based MC estimator of the tECV

We leverage the orthogonal projection property described in Eq. (11) to
approximate the CE, eliminating the necessity to sample the posterior
distributions as in Eq. (10). For a nonlinear observational map $h$, the map
$\phi_{d}$ is typically nonlinear and does not allow a closed-form solution.
In such cases, we employ the nonlinear regression method to approximate the
CE, which involves three key components. First, we select a suitable
parameterized functional subspace of finite dimension
$\mathcal{S}^{\prime}(\mathbb{R}^{m},\mathbb{R}^{n})\subset\mathcal{S}(\mathbb{R}^{m},\mathbb{R}^{n})$.
This subspace represents a restricted set of functions that are used to
approximate the map $\phi_{d}$. Second, we solve the projection problem
described in Eq. (11) within this chosen subspace. Last, we employ an
appropriate MC estimator to estimate the MSE
$\operatorname{E}\left[\left\lVert Q-f\circ{Y_{d}}\right\rVert_{2}^{2}\right]$
for $f\in\mathcal{S}^{\prime}$.

By restricting $\mathcal{S}^{\prime}$ to a linear function space, we obtain a
linear approximation of the map $\phi_{d}$, as detailed in Appendix C.1.
However, this linear approximation of the CE is known to be biased and tends
to overestimate the tECV. The remainder of this section is dedicated to
developing a suitable method for approximating the nonlinear CE.

Our primary objective in this section is to approximate the CE given a fixed
design parameter vector $d$. We will expand our approach in Section 5 to
encompass the neighborhood surrounding a given vector $d$. Let
$\bm{\theta}\in\mathbb{R}^{\beta}$ denote the vector containing the
hyperparameters of functions $f$ in the selected subspace
$\mathcal{S}^{\prime}$. We use the orthogonal projection property (Eq. (11))
to approximate the map $\phi_{d}$ using a function
$f^{*}:=f(\cdot\;;\bm{\theta}^{*})$ in the subspace $\mathcal{S}^{\prime}$,
such that

$\bm{\theta}^{*}=\arg\;\min_{\bm{\theta}}\;\operatorname{E}\left[\bigl{\lVert}Q-f({Y_{d}};\bm{\theta})\bigr{\rVert}_{2}^{2}\right].$
(13)

To simplify the notation, we introduce the MSE function
$\mathcal{M}:\mathcal{S}\rightarrow\mathbb{R}_{+}$, which is defined as

$\mathcal{M}(f)\equiv\operatorname{E}\left[\bigl{\lVert}Q-f({Y_{d}})\bigr{\rVert}_{2}^{2}\right].$
(14)

With this notation in place, the relations between functions in $\mathcal{S}$
can be summarized as:

$\mathcal{M}(f)\geq\mathcal{M}(f^{*})\geq\mathcal{M}(\phi_{d})=V(d),\;\forall
f\in\mathcal{S^{\prime}}.$

A straightforward method for estimating the MSE $\mathcal{M}(f)$ is to use the
crude Monte Carlo (MC) estimator. Let
$D_{N}=\\{\left(q^{(i)},y^{(i)}\right)\\}_{i=1}^{N}$ be an $N$-sized dataset
of _i.i.d._ samples of the RV pairs $(Q,Y_{d})$, where
$\\{q^{(i)}\\}_{i=1}^{N}$ are the _i.i.d._ samples of RV $Q$, and
$\\{y^{(i)}\\}_{i=1}^{N}$ are the corresponding _i.i.d._ samples of RV
${Y_{d}}$, which are obtained as

$\\{y^{(i)}\\}_{i=1}^{N}=\left\\{y^{(i)}\mid
y^{(i)}=h(q^{(i)},d)+{\xi}^{(i)},\quad
q^{(i)}\in\\{q^{(i)}\\}_{i=1}^{N}\right\\}.$ (15)

Here, $\\{\xi^{(i)}\\}_{i=1}^{N}$ are the _i.i.d._ samples of $\varXi$.

Let $D_{M}=\\{\left(q^{(i)},y^{(i)}\right)\\}_{i=1}^{M}$ be an $M$-sized
dataset of _i.i.d._ samples of the RV pairs $(Q,Y_{d})$, which is
statistically independent from $D_{N}$. We denote the crude MC estimator of
the MSE $\mathcal{M}(f)$ using dataset $D\in\\{D_{N},D_{M}\\}$ as:

$\widehat{\mathcal{M}}(f\mid D)=\dfrac{1}{|D|}\sum_{(q,y)\in
D}\bigl{\lVert}q-f(y)\bigr{\rVert}_{2}^{2},$ (16)

where $|D|$ is the cardinality of dataset $D$. Owing to Proposition 1, we
estimate the tECV using the following _PACE-based MC estimator_
$\widehat{V}_{d}(D_{N},D_{M})$:

$\displaystyle\widehat{V}_{d}(D_{N},D_{M})$
$\displaystyle:=\widehat{\mathcal{M}}(f(\cdot;\bm{\theta}_{D_{N}})\mid
D_{M})\quad$ (17a) $\displaystyle\text{where}\quad\bm{\theta}_{D_{N}}$
$\displaystyle=\arg\;\min_{\bm{\theta}}\;\widehat{\mathcal{M}}(f(\cdot;\bm{\theta})\mid
D_{N}),$ (17b)

which requires a numerical solution of the optimization problem stated in Eq.
(17b). We observe that unlike the double-loop IS estimator, which incurs a
computational cost proportional to the product $(N_{\text{o}}\times
N_{\text{i}})$ as explained in Section 2, the PACE-based approach exhibits a
linear cost proportional to the sizes of datasets $D_{N}$ and $D_{M}$.

Although the PACE-based approach requires approximating the CE, it completely
alleviates the need to sample different posterior densities. This allows us to
bypass the complication associated with posterior intractability. Using the
regressor, we obtain approximate values of the CE samples with minimal
computing effort. In the following subsections, we analyze the statistical
error of the PACE-based MC estimator of the tECV. Moreover, we propose a
control-variate version for this PACE-based MC estimator.

###  3.2 Estimation of errors

In this subsection, we aim to analyze the MAE of the estimator
$\widehat{V}_{d}(D_{N},D_{M})$. Let $\epsilon_{S^{\prime}}$,
$\epsilon_{\text{opt}}$, and $\epsilon_{\text{MC}}$ denote the approximation
error due to the choice of the subspace $\mathcal{S}^{\prime}$, the
optimization error for the numerical solution to the minimization problem
stated in Eq. (17b) with the finite-size dataset $D_{M}$, and the estimating
error of the MC estimator due to the finite size of the dataset $D_{N}$,
respectively, _i.e._ ,

$\displaystyle\epsilon_{S^{\prime}}$
$\displaystyle:=\operatorname{E}\left[\bigl{|}\mathcal{M}(f^{*})-\mathcal{M}(\phi_{d})\bigr{|}\right],$
(18) $\displaystyle\epsilon_{\text{opt}}$
$\displaystyle:=\operatorname{E}\left[\bigl{\lvert}\mathcal{M}(f(\cdot;\bm{\theta}_{D_{N}}))-\mathcal{M}(f^{*})\bigr{\rvert}\right],$
(19) $\displaystyle\epsilon_{\text{MC}}$
$\displaystyle:=\operatorname{E}\left[\bigl{\lvert}\widehat{\mathcal{M}}(f(\cdot;\bm{\theta}_{D_{N}})\mid
D_{M})-\mathcal{M}(f(\cdot;\bm{\theta}_{D_{N}}))\bigr{\rvert}\right].$ (20)

The MAE between the estimator $\widehat{V}_{d}(D_{N},D_{M})$ and $V_{d}$ is
bounded using the triangular inequality as follows:

$\displaystyle\operatorname{E}\left[\bigl{|}\widehat{V}_{d}(D_{N},D_{M})-V_{d}\bigr{|}\right]=$
$\displaystyle\operatorname{E}\left[\bigl{|}\widehat{\mathcal{M}}(f(\cdot;\bm{\theta}_{D_{N}})\mid
D_{M})-\mathcal{M}(\phi_{d})\bigr{|}\right]$ (21) $\displaystyle\leq$
$\displaystyle\operatorname{E}\left[\bigl{|}\widehat{\mathcal{M}}(f(\cdot;\bm{\theta}_{D_{N}})\mid
D_{M})-\mathcal{M}(f(\cdot;\bm{\theta}_{D_{N}}))\bigr{|}\right]$
$\displaystyle+\operatorname{E}\left[\bigl{|}\mathcal{M}(f(\cdot;\bm{\theta}_{D_{N}}))-\mathcal{M}(f^{*})\bigr{|}\right]$
$\displaystyle+\operatorname{E}\left[\bigl{|}\mathcal{M}(f^{*})-\mathcal{M}(\phi_{d})\bigr{|}\right]$
$\displaystyle=$
$\displaystyle\epsilon_{\text{MC}}+\epsilon_{\text{opt}}+\epsilon_{S^{\prime}}.$

Aiming at deriving an asymptotic error estimator that captures the error
evolution _w.r.t._ the size of the datasets $D_{N}$ and $D_{M}$, we make use
of the following assumptions:

######  Assumption 1.

$\mathcal{S^{\prime}}$ is a convex set that contains constant functions.

######  Assumption 2.

$\operatorname{Var}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right]=\mathcal{O}\bigl{(}2\operatorname{E}\left[\left\lVert
Q-f^{*}({Y_{d}})\right\rVert_{2}^{2}\right]^{2}\bigr{)}\;<\infty.$ (22)

######  Assumption 3.

$\epsilon_{\text{opt}}=\mathcal{O}\left(\operatorname{E}\left[\left\lvert\widehat{\mathcal{M}}(f^{*}|D_{N})-\mathcal{M}(f^{*})\right\rvert\right]\right).$
(23)

The first assumption is a common constraint imposed on the choice of the
subspace $\mathcal{S^{\prime}}$. This assumption implies that
$\operatorname{E}\left[(Q-f^{*}({Y_{d}}))^{\top}g({Y_{d}}))\right]=0$ for
every $g\in\mathcal{S^{\prime}}$. By choosing $g$ as a constant function in
this expression, we deduce that the mean of the RV $Q-f^{*}({Y_{d}})$ is zero.
We justify the second assumption by noting that for every zero-mean Gaussian
RV $X$, we have
$\operatorname{Var}\left[X^{2}\right]=\operatorname{E}\left[X^{4}\right]-(\operatorname{E}\left[X^{2}\right])^{2}=2(\operatorname{E}\left[X^{2}\right])^{2}$.
The third assumption asserts that the error of the optimization problem under
the finite-sized dataset $D_{N}$ is on the same order as the statistical error
of the MC estimator $\widehat{\mathcal{M}}(f^{*}|D_{N})$.

######  Proposition 2 (Error estimation).

Assuming that the RVs $Q$ and ${Y_{d}}$ have finite variance, and that the
Assumptions 1, 2, and 3 hold, the MAE of the PACE-based MC estimator
$\widehat{V}_{d}(D_{N},D_{M})$ satisfies

$\operatorname{E}\left[\bigl{|}\widehat{V}_{d}(D_{N},D_{M})-V_{d}\bigr{|}\right]=\mathcal{O}\Biggl{(}\Bigl{(}\dfrac{2}{\sqrt{\pi
N}}+\dfrac{2}{\sqrt{\pi M}}\Bigr{)}\operatorname{E}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right]\Biggr{)}+\epsilon_{S^{\prime}}.$
(24)

as ${N,M\ \gg 1}$.

If we further assume the approximation error $\epsilon_{S^{\prime}}=0$, then

$\displaystyle\operatorname{E}\left[\bigl{\lvert}\widehat{V}_{d}(D_{N},D_{M})-V_{d}\bigr{\rvert}\right]$
$\displaystyle=\mathcal{O}\Biggl{(}\Bigl{(}\dfrac{2}{\sqrt{\pi
N}}+\dfrac{2}{\sqrt{\pi M}}\Bigr{)}\operatorname{E}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right]\Biggr{)}$ (25)
$\displaystyle=\mathcal{O}\Biggl{(}\Bigl{(}\dfrac{2}{\sqrt{\pi
N}}+\dfrac{2}{\sqrt{\pi M}}\Bigr{)}\operatorname{E}\left[\left\lVert
Q-\phi_{d}(Y)\right\rVert_{2}^{2}\right]\Biggr{)}$
$\displaystyle=\mathcal{O}\left(\dfrac{2}{\sqrt{\pi N}}+\dfrac{2}{\sqrt{\pi
M}}\right)V_{d}$

To evaluate the accuracy of an estimator $\widetilde{V}_{d}$, we use the
relative MAE defined as

$\text{relMAE}\;(\widetilde{V}_{d})\equiv\dfrac{\operatorname{E}\left[|\widetilde{V}_{d}-V_{d}|\right]}{V_{d}}.$
(26)

When $\epsilon_{S^{\prime}}=0$, the relative MAE of the estimator
$\widehat{V}_{d}(D_{N},D_{M})$ is given by

$\displaystyle\text{relMAE}\;(\widehat{V}_{d})$
$\displaystyle\equiv\dfrac{\operatorname{E}\left[\bigl{|}\widehat{V}_{d}(D_{N},D_{M})-V_{d}\bigr{|}\right]}{V_{d}}$
(27) $\displaystyle=\mathcal{O}\Bigl{(}\dfrac{2}{\sqrt{\pi
N}}+\dfrac{2}{\sqrt{\pi M}}\Bigr{)},$

for $N,M\gg 1$. When the upper bounds of the errors $\epsilon_{S^{\prime}}$,
$\epsilon_{\text{opt}}$, and $\epsilon_{\text{MC}}$ are available, it is
possible to derive more precise bound for the MAE of the PACE-based MC
estimator. However, those upper bounds are problem-dependent and beyond the
scope of this study.

###  3.3 Reduced variance estimator

We often model the measurement error RV $\varXi$ using simple parameterized
distributions such as the Gaussian distribution, Poisson distribution, or
binomial distribution. Consequently, sampling the RV $\varXi$ is
computationally inexpensive. Considering this observation, we augment datasets
$D_{N}$ and $D_{M}$ to reduce the statistical error of the PACE-based MC
estimators. The augmented datasets are obtained as follows:

$\displaystyle D_{N}^{\text{rv}}$
$\displaystyle=\left\\{\left(q^{(i)},y^{(i,j)}\right)\mid
y^{(i,j)}=h(q^{(i)},d)+{\xi}^{(i,j)}\right\\}_{i=1,\dots,N,\;j=1,\dots,a},$
(28) $\displaystyle D_{M}^{\text{rv}}$
$\displaystyle=\left\\{\left(q^{(i)},y^{(i,j)}\right)\mid
y^{(i,j)}=h(q^{(i)},d)+{\xi}^{(i,j)}\right\\}_{i=1,\dots,M,\;j=1,\dots,a},$
(29)

where $q^{(i)}$ and $\xi^{(i,j)}$ are the _i.i.d._ samples of RVs $Q$ and
$\varXi$, respectively, and $a\in\mathbb{N}_{+}$ represents the augmentation
multiplier. By substituting the augmented datasets $D_{N}^{\text{rv}}$ and
$D_{M}^{\text{rv}}$ for $D_{N}$ and $D_{M}$, respectively, in Eq. (16), we
obtain the reduced variance estimators of the MSEs
$\widehat{\mathcal{M}}(f\mid D_{N})$ and $\widehat{\mathcal{M}}(f\mid D_{M})$
as follows:

$\displaystyle\widehat{\mathcal{M}}^{\text{rv}}(f\mid D_{N})$
$\displaystyle=\dfrac{1}{|D_{N}|\times a}\sum_{(q,y)\in
D^{\text{rv}}_{N}}\left\lVert q-f(y)\right\rVert_{2}^{2}$ (30)
$\displaystyle\widehat{\mathcal{M}}^{\text{rv}}(f\mid D_{M})$
$\displaystyle=\dfrac{1}{|D_{M}|\times a}\sum_{(q,y)\in
D^{\text{rv}}_{M}}\left\lVert q-f(y)\right\rVert_{2}^{2}.$

Appendix D explains the reduction in the statistical error obtained with the
use of the estimator $\widehat{\mathcal{M}}^{\text{rv}}$ relative to the crude
one $\widehat{\mathcal{M}}$ .

##  4 Numerical experiment: Linear-Gaussian setting

In this section, we present the results of numerical experiments to
demonstrate the error estimation outlined in Section 3.2. We focus on a
linear-Gaussian scenario, where both distributions of the prior and the
measurement error are Gaussian and the observational map is linear. This
setting is advantageous for error analysis because it allows for an analytical
solution to the tECV. We will expand our numerical experiments to a broader
nonlinear setting in Section 6.

To assess the effectiveness of the PACE-based approach, we conduct numerical
experiments in two different scenarios. In Section 4.1, we apply our method to
a one-dimensional setting, _i.e._ , $q\in\mathbb{R}$, and in Section 4.2, we
evaluate our approach in a high-dimensional scenario. Furthermore, we compare
the performance of the PACE-based approach with that of the IS-based approach,
by examining the effect of two conditions: (1) increasing the dimension of the
inferred parameter vector $q$, and (2) reducing the measurement error
variance. These factors significantly exacerbate the intractable property of
the posterior distribution.

###  4.1 One-dimensional case

In this section, we consider the following setting:

$\displaystyle h(Q,d)=\dfrac{Q}{(d-0.5)^{2}+1},\quad d\in[0,1],$ (31a)
$\displaystyle
Q\sim\mathcal{N}(0,\sigma_{q}^{2}),\;\text{with}\quad\sigma_{q}=2,$ (31b)
$\displaystyle\varXi\sim\mathcal{N}(0,\sigma_{\xi}^{2}),$ (31c)

where $\mathcal{N}(0,\sigma^{2})$ is a Gaussian distribution with a mean of
zero and variance $\sigma^{2}$. We analyze two cases of variances in
measurement error: $\sigma^{2}_{\xi}=0.01^{2}$ and
$\sigma^{2}_{\xi}=0.001^{2}$. Given that the observational map is linear in
terms of $q$ and both the prior and measurement error distributions are
Gaussian, the CE map $\phi_{d}$ defined in Eq. (11) is linear and possesses a
closed-form expression. Thus, we obtain

$\phi_{d}(y)=\frac{a\sigma_{q}^{2}}{a^{2}\sigma_{q}^{2}+\sigma_{\xi}^{2}}\,y.\;\text{where}\;a=\dfrac{1}{(d-0.5)^{2}+1}.$
(32)

Appendix C.1 provides a detailed description of this linear approximation.

To implement our method in this setting, we employ linear regression to
empirically approxi-mate the CE by solving Eq. (11) within the subspace of
linear functions $S^{\prime}$, given the sample dataset
$D_{N}=\\{(q^{(i)},y^{(i)})\\}_{i=1}^{N}$. For detailed information, please
refer to Appendix C.2. We then utilize the obtained CE approximation to
estimate the tECV using Eq (30). To assess the accuracy of our estimation, we
compute the empirical relative MAE metric defined in Eq.(26). In particular,
we aim to compare the relative MAE with the estimation in Eq.(27), considering
the absence of bias error $\epsilon_{S^{\prime}}$ in the linear-Gaussian
setting. Additionally, we implement the IS-based approach described in
Appendix A to estimate the tECV and to evaluate its relative MAE.

Based on the result of Proposition 2, setting $N\equiv M$ is a straightforward
decision for the PACE-based approach. In the IS-based approach, the ECV is
estimated using a double-loop MC algorithm, which is detailed in Appendix A.
Because the posterior variance remains invariant _w.r.t_ $y$ in the linear-
Gaussian setting, we set the number of samples used by the outer loop to
$N_{\text{o}}=1$. In the PACE-based approach, the total number of samples is
$N+M$, while for the IS-based approach, it is $N_{\text{i}}+1$, where
$N_{\text{i}}$ represents the number of processed by the inner loop. Fig. 1
shows the relative MAE for $d=0.5$ for two different measurement variances,
_i.e._ , $\sigma_{\xi}^{2}=0.01^{2}$ and $\sigma_{\xi}^{2}=0.001^{2}$. Fig. 2
showcases the tECV for different values of $d$ and for
$\sigma_{\xi}^{2}=0.01^{2}$. To estimate the relative MAE, we perform 1000
statistically independent simulations.

(a)

(b)

Figure 1: Comparison of the PACE- and IS-based approaches for estimating the
tECV for $d=0.5$: (a) $\varXi\sim\mathcal{N}(0,0.01^{2})$ and (b)
$\varXi\sim\mathcal{N}(0,0.001^{2})$.

As can be observed in Fig. 1, the PACE-based method requires considerably
fewer samples than the IS-based approach for achieving a similar relative MAE
value. For example, for $\sigma_{\xi}^{2}=0.01^{2}$ without applying data
augmentation, the required number of samples is reduced by approximately 30
fold, and for $\sigma_{\xi}^{2}=0.001^{2}$, it is reduced by about 200 fold.
When data augmentation is used, our approach further reduces the number of
samples by nearly two additional orders of magnitude.

Overall, the evolution of the empirical relative MAE confirms the theoretical
error estimation obtained in Proposition 2, _i.e._ ,
$\mathcal{O}\bigl{(}\frac{2}{\sqrt{\pi N}}+\frac{2}{\sqrt{\pi M}}\bigr{)}$.
Particularly, we have
$\;\operatorname{Var}\left[\left(Q-f^{*}(Y_{d})\right)^{2}\right]=2\operatorname{E}\left[\left(Q-f^{*}({Y_{d}})\right)^{2}\right]^{2}$,
which validates Assumption 2. As the CE in this example is a linear function,
$\epsilon_{opt}$ is much smaller than the estimated order
$\Bigl{(}\operatorname{E}\left[\left\lvert\widehat{\mathcal{M}}(f^{*}|D_{N})-\mathcal{M}(f^{*})\right\lvert\right]\Bigr{)}$
that is assumed in the Assumption 3. Consequently, the obtained relative MAE
is consistently smaller than the theoretical error estimation of
$\bigl{(}\frac{2}{\sqrt{\pi N}}+\frac{2}{\sqrt{\pi M}}\bigr{)}$.

We also observe that the relative MAE in the PACE-based approach is
_invariant_ _w.r.t._ the measurement error variance. In contrast, the IS-based
method requires a considerably larger number of samples for a smaller
measurement error variance. Particularly, a reduction in the measurement error
variance decreases the evidence terms of Bayes’ formula. Consequently, the
statistical error inherent in the IS-based approach becomes amplified.

(a)

(b)

Figure 2: Comparison of the PACE- and IS-based approaches for estimating ECV
with $\varXi\sim\mathcal{N}(0,0.01^{2})$ using 10,000 samples: (a) PACE-based
approach without data augmentation and (b) IS-based approach. The error bar
plot depicts the mean and standard deviation of the empirical ECV values
estimated from 1000 statistically independent MC simulations.

Fig. 2 clearly illustrates that the PACE-based approach accurately estimates
the A-optimal DOE, _i.e._ , $d_{A}=0.5$, with 10,000 samples. However, the IS-
based approach results in significant errors in identifying the optimal design
due to its statistical errors in estimating the tECV.

###  4.2 High-dimensional case

In this section, we address the high-dimensional linear-Gaussian scenario. The
setting can be described as follows:

$\displaystyle h(q,d)$ $\displaystyle=\dfrac{1}{(d-0.5)^{2}+1}\;Q,\quad
d\in[0,1],$ (33) $\displaystyle Q$
$\displaystyle\sim\mathcal{N}(0_{n},\mathbf{I}_{n}),$ $\displaystyle\varXi$
$\displaystyle\sim\mathcal{N}(0_{n},0.1^{2}\mathbf{I}_{n}),$

where $0_{n}\in\mathbb{R}^{n}$ and $\mathbf{I}_{n}\in\mathbb{R}^{n\times n}$
represent a zero vector and an identity matrix, respectively. We perform an
analysis similar to the one presented in Section 4.1.

(a)

(b)

(c)

(d)

Figure 3: Comparison of the PACE- and IS-based approaches with high-
dimensional linear examples: (a) $n=2$, (b) $n=5$, (c) $n=10$, and (d) $n=20$.
For $n=20$, all 100 statistically independent MC simulations that uses the IS-
based approach return a _not-a-number_ error. Consequently, the numerical
result of the IS-based approach for $n=20$ is not available.

Fig. 3 illustrates the relative MAE for different values of
$n\in\\{2,5,10,20\\}$. The numerical results demonstrate that the IS-based
approach requires an extremely large number of samples for high-dimensional
cases. Specifically, for $n\geq 5$, we do not observe the convergence of the
IS-based approach. For $n=20$, the IS-based approach frequently returns a
_not-a-number_ error due to the decreasing evidence term in the Bayes’ theorem
(see Eq. (3)). In contrast, our PACE-based approach effectively mitigates the
curse of dimensionality.

##  5 Stochastic optimization for the continuous design domain

The PACE approach can be used directly to solve the A-ODE problem in a
discrete and finite design domain. By estimating the tECV for each
experimental design candidate using PACE, we identify the optimal solution as
the design with the minimum tECV. In this section, we shift our focus to the
continuous design domain. We assume that the gradient of the measurement map
_w.r.t._ the design parameters, $\nabla_{d}h$, exists and can be computed
numerically. To solve the A-ODE problem efficiently, we apply stochastic
gradient descent techniques on top of the PACE-based MC estimator of tECV.

When stochastic gradient optimization algorithms are used to solve the
A-optimal DOE prob-lem, the computation of the gradient of the tECV _w.r.t._
the design parameters, denoted as ${\nabla_{d}V_{d}}$, can pose a challenge.
To address this issue, we propose a nonlocal approximation of the CE. Our
approach involves approximating the CE within a neighborhood surrounding a
given design parameter vector to enable the evaluation of ${\nabla_{d}V_{d}}$.
We utilize ANNs as the approximation tool, although other regression
techniques can also be incorporated into our method.

Let
$f_{\text{N}}(\cdot;\bm{\theta}):\mathbb{R}^{m+\delta}\rightarrow\mathbb{R}^{n}$
be an ANN map, where $\bm{\theta}$ denotes the hyperparameters of the map. The
input of the proposed ANN is the concatenation of an $m$-dimensional vector
$y$ and a $\delta$-dimensional vector $d$. By combining Eqs. (7), (11), and
(12) the A-optimal design parameter vector, $d_{\text{A}}$, can be obtained by
solving the following optimization problem

$d_{\text{A}}=\arg_{d}\;\min_{d,\bm{\theta}}\;\operatorname{E}\left[\left\lVert
Q-f_{\text{N}}(Y_{d},d;\bm{\theta})\right\rVert_{2}^{2}\right]$ (34)

which minimizes the MSE _w.r.t._ the ANN’s hyperparameters and design
parameters. Notably, the optimization problem that approximates the CE and the
problem that seeks the A-optimal DOE minimize the same objective function.

During the course of the iterative optimization process to solve Eq. (34), the
ANN may need to be retrained when the algorithm provides a new candidate
vector for $d_{\text{A}}$. Here, we employ the transfer learning technique,
_i.e._ , the weights trained for the previous design candidate are used as the
starting point for the current training process. This simple transfer learning
technique efficiently reduces the training time and the amount of data
required. The reminder of this section provides a more detailed explanation of
the nonlocal approximation of the CE in Section 5.1 and discusses the
algorithms that are used for implementing our method in Section 5.2.

###  5.1 Nonlocal approximation of the conditional expectation

Given a design candidate $d_{\kappa}$, we use a PDF $w_{\kappa}$ defined over
the design domain $\mathcal{D}$ such that
$w_{\kappa}(d):=\frac{w(d,d_{\kappa})}{c}$. Here $w(d,d_{\kappa})$ is a weight
function such as a kernel, and $c$ is a normalization constant defined as
$c:=\int_{\mathcal{D}}{w(d,d_{\kappa})}\mathrm{d}d\,<\infty$. In our notation,
$\kappa\in\mathbb{N}_{+}$ denotes the iteration counter of the main procedure,
which will be discussed in Section 5.2.

To approximate the CE in the neighborhood of a design $d_{\kappa}$, we
determine the hyperparameters $\bm{\theta}_{\kappa}$ of the ANN $f_{\text{N}}$
by solving the following optimization problem:

$\bm{\theta}_{\kappa}=\arg\quad\min_{\bm{\theta}}\quad{L}_{\kappa}(\bm{\theta}),$
(35)

where $L_{\kappa}$ is the weighted MSE function defined as

$L_{\kappa}(\bm{\theta})=\int_{\mathcal{D}}\operatorname{E}\left[\left\lVert
Q-f_{\text{N}}(Y_{d},d;\bm{\theta})\right\rVert_{2}^{2}\right]w_{\kappa}(d)\mathrm{d}d.$
(36)

We use the variance reduction technique, similar to Eq. (30), to approximate
the function $L_{\kappa}$ as

$\widehat{L}^{\text{vr}}_{\kappa}(\bm{\theta})=\dfrac{1}{N\times
a}\sum_{i=1}^{N}\sum_{j=1}^{a}\left\rVert
q^{(i)}-f_{\text{N}}(y^{(i,j)},d^{(i)};\bm{\theta})\right\rVert_{2}^{2},$ (37)

where $\\{q^{(i)}\\}_{i=1}^{N}$ are the _i.i.d._ samples of RV $Q$,
$\\{d^{(i)}\\}_{i=1}^{N}$ are the _i.i.d_ samples generated according to the
PDF $w_{\kappa}$, and $y^{(i,j)}=h(q^{(i)},d^{(i)})+\xi^{(i,j)}$, where
$\xi^{(i,j)}$ is the _i.i.d._ sample of RV $\varXi$.

A simple choice for the kernel function is the density function of a
multivariate Gaussian distribution, where the marginal variances are used as
tuning parameters. Intuitively, increasing the marginal variances of PDF
$w_{\kappa}$ reduces the error of the CE approximation over the entire design
domain; however, this also requires a considerably large number of samples to
accurately estimate the weighted MSE $\widehat{L}^{\text{vr}}_{\kappa}$. In
our algorithm, which is discussed in more detail in Section 5.2, for the
$\kappa$-_th_ iteration, we only need to evaluate the tECV and its derivative
in the neighborhood of $d_{\kappa}$. Therefore, a required standard deviation
of PDF $\eta_{\kappa}$ is much smaller than the characteristic length of the
design domain, which improves the computational efficiency by reducing the
number of samples required to estimate the weighted MSE
$\widehat{L}^{\text{vr}}_{\kappa}$.

###  5.2 Algorithm

In this subsection, we present a stochastic gradient descent algorithm to
solve Eq. (34). The algorithm utilizes the nonlocal approximation of the CE to
estimate the tECV $V_{d}$ and its derivative $\nabla_{d}V_{d}$. A typical
$\kappa$-_th_ iteration of the proposed algorithm consists of two steps:

  * i)

Given a design parameter candidate $d_{\kappa-1}$, we apply _transfer
learning_ and solve Eq. (35) to update the hyper-parameters
$\bm{\theta}_{\kappa}$,

  * ii)

We update the design parameter candidate using the gradient of the MSE _w.r.t_
$d$, _i.e._ , $\nabla_{d}\,\operatorname{E}\left[\lVert
Q-f_{\text{N}}(Y_{d},d;\bm{\theta}_{\kappa})\rVert_{2}^{2}\right]$. This step
returns an updated candidate $d_{\kappa}$.

Algorithm 1 summarizes the main procedure, which calls the suboptimization
procedures that correspond to steps i) and ii) described in Algorithms 2 and
3, respectively. The suboptimization procedures in each step are executed
using the Adam algorithm [24].

Algorithm 1 Main

1:Observational model, its derivatives ($\nabla_{d}h$), and number of
iteration $K$

2:Initiate a value for $d_{0}$

3:for $\kappa$ from 1 to $K$ do

4: approximating the conditional expectation using Algorithm 2 $\rightarrow$
$\bm{\theta}_{\kappa}$

5: solving $\arg\;\min_{d}\operatorname{E}\left[\lVert
Q-f_{\text{N}}(Y_{d},d;\bm{\theta}_{\kappa})\rVert_{2}^{2}\right]$ using
Algorithm 3 $\rightarrow$ $d_{\kappa}$

6:end for

7:return $d_{\text{A}}\leftarrow d_{K}$

In Algorithm 2, we generate _i.i.d._ samples $\\{d^{(i)}\\}_{i=1}^{N}$,
$\\{q^{(i)}\\}_{i=1}^{N}$, and $\\{\xi^{(i)}\\}_{i=1}^{N}$ following PDFs
$w_{\kappa}$, $\pi_{Q}$, and $\pi_{\varXi}$, respectively. Next, we use those
samples to compute the samples of the observational RV,
$y^{(i,j)}=h(q^{(i)},d^{(i)})+\xi^{(i,j)}$ for $j=1,\dots,a$. Finally, we
employ the Adam algorithm to train the ANN. We use the dataset
$\left\\{\left(q^{(i)},d^{(i)},y^{(i,j)}\right)\right\\}_{i=1,\dots,N;\;j=1,\dots
a}$ and the loss function $\widehat{L}^{\text{vr}}_{\kappa}$ (defined in Eq.
(37)) for this training.

To reduce the numbers of epochs ($e_{1}$) and data samples ($N\times a$), we
employ the _transfer learning_ technique [25, 26], which reuses the
hyperparameters obtained from the previous iteration,
$\bm{\theta}_{\kappa-1}$, as the initial values for the current iteration
$\kappa$.

Algorithm 2 Nonlocal approximation of the CE operator

1:$d_{\kappa}$, weight function $w(\cdot;d_{\kappa})$, number of epoch
$e_{1}$, and learning rate $\alpha_{1}$

2:Generate $\\{d^{(i)}\\}_{i=1}^{N}\sim w_{\kappa}(d)$,
$\\{q^{(i)}\\}_{i=1}^{N}\sim\pi_{Q}$,
$\\{\xi^{(i,j)}\\}_{i=1,\dots,N;\;j=1,\dots a}\sim\pi_{\varXi}$

3:Evaluate $\\{h(q^{(i)},d^{(i)})\\}_{i=1}^{N}$

4:Compute samples $\\{y^{(i,j)}\\}_{i=1,\dots,N;\;j=1,\dots a}$ where
$y^{(i,j)}=h(q^{(i)},d^{(i)})+\xi^{(i,j)}$

5:Collect data into a dataset
$D^{\text{vr}}_{N}=\left\\{\left(q^{(i)},d^{(i)},y^{(i,j)}\right)\right\\}_{i=1,\dots,N;\;j=1,\dots
a}$

6:Train the ANN using the initial weights $\bm{\theta}_{\kappa-1}$ with
$e_{1}$ epoch, Adam algorithm, learning rate $\alpha_{1}$, loss function
$\widehat{L}^{\text{vr}}_{\kappa}$ defined in Eq. (37) and dataset
$D^{\text{vr}}_{N}$

7:return Hyper-parameters $\bm{\theta}_{\kappa}$ of ANN $f_{\text{N}}$

In Algorithm 3, we generate _i.i.d._ samples $\\{q^{(i)}\\}_{i=1}^{M}$
following PDF $\pi_{Q}$. The dataset $\\{q^{(i)}\\}_{i=1}^{M}$ is divided into
different batches of size $M_{b}$, denoted as $\\{B_{1},\dots,B_{\beta}\\}$.
For each batch $B$ and design candidate $d_{0}$, we compute the corresponding
samples of the observational RV as $y^{(i,j)}=h(q^{(i)},d_{0})+\xi^{(i,j)}$,
where $q^{(i)}\in B$, $j=1,\dots,a$, and $\xi^{(i,j)}$ are the _i.i.d._
samples of the RV $\varXi$. We then update the design parameter using

$d^{\prime}=d_{0}-\dfrac{\alpha_{2}}{M_{b}\times
a}\sum_{i=1}^{M_{b}}\sum_{j=1}^{a}\nabla_{d}\Bigl{[}\lVert
q-f_{\text{N}}(h(q,d)+\xi,d)\rVert_{2}^{2}\Bigr{]}(q^{(i)},y^{(i,j)},d_{0}),$
(38)

where $\alpha_{2}$ is a learning rate. Here, $\nabla_{d}\Bigl{[}\lVert
q-f_{\text{N}}(h(q,d)+\xi,d)\rVert_{2}^{2}\Bigr{]}(q^{(i)},y^{(i,j)},d_{0})$
represents the gradient of $\lVert q-f_{\text{N}}(h(q,d)+\xi,d)\rVert_{2}^{2}$
_w.r.t._ $d$ and is evaluated at $(q^{(i)},y^{(i,j)},d_{0})$. The value of
$d^{\prime}$ is then assigned to $d_{0}$ for the next batch.

Algorithm 3 Minimizing the MSE $\operatorname{E}\left[\lVert
Q-f_{\text{N}}(Y_{d},d;\bm{\theta}_{\kappa})\rVert_{2}^{2}\right]$ _w.r.t._
$d$

1:Observational model ($h$) and its derivatives ($\nabla_{d}h$), map
$f_{\text{N}}(\cdot;\bm{\theta}_{\kappa})$ from Algorithm 2, $d_{\kappa}$,
number of samples $M$, number of epoch $e_{2}$, batch-size $M_{b}$, and
learning rate $\alpha_{2}$

2:$d_{0}\leftarrow d_{\kappa}$

3:Generate samples $\\{q^{(i)}\\}_{i=1}^{M}$ following prior PDF $\pi_{Q}$

4:for $e=1$ to $e_{2}$ do

5: for $b=1$ to $M/M_{b}$ do

6: Generate samples $\\{\xi^{(i,j)}\\}_{i=1,\dots,M_{b};\;j=1,\dots,a}$ of the
measurement error $\varXi$

7: Evaluate $y^{(i,j)}=h(q^{(i)},d_{0})+\xi^{(i,j)}$

8: Updating the design candidate using Adam algorithm based on Eq. (38)

9: $d_{0}\leftarrow d^{\prime}$

10: end for

11:end for

12:return $d_{\kappa+1}\leftarrow d_{0}$

In Appendix E, we use the chain rule to derive the formulation of the gradient
$\nabla_{d}\Bigl{[}\lVert
q-f_{\text{N}}(h(q)+\xi,d)\rVert_{2}^{2}\Bigr{]}(q^{(i)},y^{(i,j)},d_{0})$ in
terms of the gradients of the observational map $h$ and the ANN $f_{\text{N}}$
_w.r.t_ $d$. Notably, the gradients of the ANN function $f_{\text{N}}$ are
obtained using the backpropagation [27].

##  6 Numerical experiment: Electrical impedance tomography

Electrical impedance tomography (EIT) aims to determine the conductivity of a
closed body by measuring the potential of electrodes placed on its surface in
response to applied electric currents. In this section, we revisit an
experimental design problem previously studied in [16, 17, 18]. The objective
of the experiment is to recover the fiber orientation in each orthotropic ply
of a composite laminate material. This is achieved by measuring the potential
when applying low-frequency electric currents through electrodes attached to
the body surface. The formulation of the EIT experiment is described in
Section 6.1, where we adopt the entire electrode model from [28].

Previous studies employed the maximization of EIG as the optimality criterion
and estimated the EIG using the IS technique combined with the Laplace
approximation of the posteriors. In Section 6.2, we analyze our PACE-based
approach for estimating the tECV and compare the obtained results with those
using the IS-based technique. Subsequently, we utilize the optimization
algorithms discussed in Section 5 to solve the DOE problem and present the
numerical results.

###  6.1 Formulation and finite element model of EIT experiment

We consider a rectangular body ${B}$ which has boundary $\partial B$ and
consists of two plies ${B}=\bigcup_{i=1}^{2}{B}_{i}$. The exterior boundary
surfaces are equipped with $N_{\text{el}}=10$ electrodes on the area
$E_{1},\dots,E_{N_{\text{el}}}$.

These electrodes are used to apply known currents, and they also allow us to
measure the electric potential. The configuration used in this study is
illustrated in Fig. 4.

Figure 4: Illustration of the EIT experiment for a two-ply composite sample.
The black rectangles indicate the electrode positions, and
$I_{1},\dots,I_{10}$ are their applied currents.

The fiber orientation over body ${B}$ is modeled as a random field
$\eta:{B}\times\varOmega\rightarrow[-\pi,\pi]$ as follows:

$\eta(x,\omega)=\begin{cases}&\eta_{1}(\omega)\quad\text{in}\;{B}_{1},\\\
&\eta_{2}(\omega)\quad\text{in}\;{B}_{2}.\end{cases}$ (39)

The RV $Q$ in this example is a random vector valued in $[-\pi,\pi]^{2}$ that
represents the prior distribution of the fiber orientations of the plies,
defined as

$Q(\omega):=[\eta_{1}(\omega),\eta_{2}(\omega)]^{\top}.$ (40)

We represent the quasistatic current flux and the quasi-static potential as
random fields as $\zeta:{B}\times\varOmega\rightarrow\mathbb{R}^{3}$ and
$u:{B}\times\varOmega\rightarrow\mathbb{R}$, respectively. The electrical
current field and the potential field satisfy the following PDE:

$\begin{aligned} \nabla\cdot\zeta&=0,\\\ \zeta&=\sigma(\eta)\;\nabla
u,\end{aligned}\;$ (41)

where $\sigma$ is the conductivity that depends on the fiber orientation of
the plies. The relation between the conductivity and the fiber orientations
$\eta$ is given as

$\sigma(\eta)=\mathcal{R}(\eta)^{\top}\cdot\bar{\sigma}\cdot\mathcal{R}(\eta),$
(42)

where $\mathcal{R}(\eta)$ is the rotation matrix,

$\mathcal{R}(\eta)=\begin{bmatrix}\cos(\eta)&0&\sin(\eta)\\\ 0&1&0\\\
-\sin(\eta)&0&\cos(\eta)\end{bmatrix}.$ (43)

Here, $\bar{\sigma}$ is a $3\times 3$ constant matrix, which we set to
$\bar{\sigma}=\operatorname{diag}(10^{-2},10^{-3},10^{-3})$. The boundary
conditions are given by

$\displaystyle\zeta\cdot n=0$
$\displaystyle,\quad\text{on}\,\partial{B}\setminus(\cup E_{l}),$ (44)
$\displaystyle\int_{E_{l}}\zeta\cdot n\mathrm{d}{}x=I_{l}$
$\displaystyle,\quad l=1,\dots,N_{\text{el}},$
$\displaystyle\frac{1}{E_{l}}\int_{E_{l}}u\mathrm{d}{}x+z_{l}\int_{E_{l}}\zeta\cdot{n}\mathrm{d}{}x=U_{l}$
$\displaystyle,\quad l=1,\dots,N_{\text{el}},$

where $n$ is the unit outward normal, and $I_{l}$ and $U_{l}$ are the applied
current and observable electric potential at electrode $E_{l}$, respectively.
The first boundary condition states the no-flux condition, and the second one
implies that the total injected current through each electrode is known. The
third boundary condition captures the surface impedance effect, meaning that
the shared interface between the electrode and material has an infinitesimally
thin layer with a surface impedance of $z_{l}$. In this study, we set to
$z_{l}=0.1$ for the surface impedance. We use the following two constraints
(Kirchhoff law of charge conservation and ground potential condition) to
guarantee the well-posedness:

$\sum_{l=1}^{N_{\text{el}}}I_{l}=0\quad\text{and}\quad\sum_{l=1}^{N_{\text{el}}}U_{l}=0.$
(45)

To solve the PDE stated in Eq. (41), we develop a 2D FE model using a
quadratic mesh of size $50\times 6$ elements. Details on the weak formulation
of the FE model are described in Appendix F.

###  6.2 Numerical results

##### Setting

Ten electrodes are placed on the surface of the rectangular domain
${B}=[0,20]\times[0,2]$, five on the top and five on the bottom of the
surface. The fiber orientations $\eta_{1}$ in ply ${B}_{1}$ and $\eta_{2}$ in
ply ${B}_{2}$ are the quantities of interest with the following assumed prior
distributions:

$\eta_{1}\sim\mathcal{U}\left(\frac{\pi}{4.5},\frac{\pi}{3.5}\right),\quad\eta_{2}\sim\mathcal{U}\left(-\frac{\pi}{3.5},-\frac{\pi}{4.5}\right).$
(46)

The observational model is given as

$\displaystyle Y_{d}$ $\displaystyle=h(Q,d)+\varXi,$ (47)
$\displaystyle\text{where}\quad h(Q,d)$
$\displaystyle:=[U_{1}(Q,d),\dots,U_{10}(Q,d)]^{\top},$

and $Q$ is defined in Eq. (40). We consider two cases of the observational
error distribution, which are $\varXi\sim\mathcal{N}(0,10^{2}\mathbf{I}_{10})$
and $\varXi\sim\mathcal{N}(0,3^{2}\mathbf{I}_{10})$, _i.e._ , the error
standard deviations are about 5% and 1.5% of the ground-truth value,
respectively. We treat the applied electrode currents as the design parameters
denoted as $d:=[I_{1},\dots,I_{9}]^{\top}$, where we apply the condition
$I_{10}=-\sum_{l=1}^{9}I_{l}$ in accordance with the Kirchhoff law of charge
conservation. We seek the A-optimal DOE parameterized by vector
$d\in[-1,1]^{9}$.

To accelerate our analysis, we construct a surrogate model using an ANN to
approximate the observational model $h$. The input of this surrogate model
consists of the fiber orientations $[\eta_{1},\eta_{2}]^{\top}$ and the
applied electrode currents $[I_{1},\dots,I_{9}]^{\top}$, and its output is the
vector of the potential at ten electrodes $[U_{1},\dots,U_{10}]^{\top}$. The
surrogate model is trained on a dataset obtained by running the FE model for a
large number of input samples such that
$[\eta_{1},\eta_{2}]^{\top}\sim\mathcal{U}\left([\frac{\pi}{4.5},\frac{\pi}{3.5}]\times[-\frac{\pi}{3.5},-\frac{\pi}{4.5}]\right)$
and $[I_{1},\dots,I_{9}]^{\top}\sim\mathcal{U}([-1,1]^{9})$ conditioned by
$\lvert\sum_{l=1}^{9}I_{l}\rvert\leq 1$.

####  6.2.1 Analysis of error in estimating the tECV

Here we study the errors in estimating the tECV for a specified design
parameter vector. For the PACE-based approach, we implement three models: i)
using the linear approximation for the CE (see Appendix C.1), ii) using the
ANN for approximating the CE without data augmentation, and iii) using the ANN
for approximating the CE and applying the data augmentation method described
in Section 3.3. The ANNs in cases ii and iii have two hidden layers of 100
neurons each. The input and output layers have 10 and 2 neurons, respectively,
which corresponds to the dimensions of vectors $y$ and $q$, respectively.
Dataset $D_{N}$ is split into two datasets for training and testing with a
size ratio of 1:1. We use the Adam algorithm [24] with a learning rate of
$0.0005$ and a batch size of $100$. The algorithm is set to run for a maximum
of $10,000$ epochs and features an early-stop mechanism that is activated when
the loss function shows no further improvement. To compute the tECV, we fix
$M=N$ based on the Proposal 2.

As a closed form of the tECV is not available in this example, we employ the
IS approach with $200,000$ _i.i.d_ samples, and we use the obtained result as
the reference, _i.e._ , $V$ in Eq. (26). Fig. 5 depicts the relative MAEs in
estimating the tECV for the optimal design found in [17]. To obtain these
estimated relative MAEs, we perform 50 statistically independent simulations
to empirically compute the expectation operator in Eq. (26).

(a)

(b)

Figure 5: Relative MAE for estimating the tECV for $d=d_{\text{optimal}}$ (a)
$\varXi\sim\mathcal{N}(0,10^{2}\mathbf{I}_{10})$, and (b)
$\varXi\sim\mathcal{N}(0,3^{2}\mathbf{I}_{10})$. For the PACE-based approach,
the number of samples is equal to $N+M$. For the IS-based approach, the number
of samples is equal to $(N_{\text{i}}+1)\times N_{\text{o}}$. For the sake of
simplicity, we fix $N=M$ and $N_{\text{i}}=N_{\text{o}}$.

Because the measurement map is nonlinear, the CE linear approximation contains
bias error, which results in a relative MAE of 20%. In contrast, we observe
negligible bias error when we use the ANN for approximating the CE. When
comparing models ii and iii, we observe that the data augmentation method
described in Section 3.3 plays a crucial role in reducing the estimation
error, particularly for a small number of samples, as shown in Fig.5. For the
case $N+M=200$, applying the data augmentation reduces the relative MAE by
more than three-fold for $\varXi\sim\mathcal{N}(0,10^{2}\mathbf{I}_{10})$.

With the IS-based approach, estimating the tECV requires a double-loop MC
simulation, where the numbers of samples for the outer and inner loops are
$N_{\text{o}}$ and $N_{\text{i}}$, respectively (see Appendix A). Here, we set
$N_{\text{o}}=N_{\text{i}}$ for simplicity. Optimizing the ratio
$N_{\text{o}}:N_{\text{i}}$ is problem dependent and beyond the scope of this
study. Compared with the IS-based approach, the PACE-based approach combined
with data augmentation significantly decreases the number of required samples.
For example in Fig. 5(a), to reach a relative MAE of 0.9%, the IS approach
requires more than 4000 samples, while the PACE-based approach using data
augmentation needs only 1000 samples. Notably, for the case
$\varXi\sim\mathcal{N}(0,3^{2}\mathbf{I}_{10})$ depicted in Fig. 5(b),
although the computational efficiency of the PACE approach does not exhibit a
significant change compared with the case in Fig. 5(a), the performance of the
IS approach deteriorates substantially. In this case, the statistical error of
the IS approach is even greater than that of the PACE approach using the
linear approximation. Moreover, the relative MAE of the PACE approach verifies
its theoretical estimation, $\mathcal{O}\bigl{(}\sqrt{\frac{2}{\pi
N}}+\sqrt{\frac{2}{\pi M}}\bigr{)}$ stated in Proposition 2.

####  6.2.2 Minimization of the tECV

In this section we apply the algorithms developed in Section 5 to determine
the A-optimal DOE. We approximate the CE nonlocally using an ANN having two
hidden layers. The input layers have 19 neurons that correspond to the sum of
the dimensions of vectors $d\in[-1,1]^{9}$ and $y\in\mathbb{R}^{10}$. The
number of neurons in the output layer and each hidden layer are 2 and 100,
respectively. We choose the density of the normal distribution
$\mathcal{N}(\text{0}_{9},0.2\times\mathbf{I}_{9})$ as the kernel function
used to evaluate the weighted loss function $L_{\kappa}$ in Eq. (36). From
this point onward, we keep $\varXi\sim\mathcal{N}(0,10^{2}\mathbf{I}_{10})$ as
was done [18].

The main procedure described in Algorithm 1 is implemented with 20 iterations.
In each iteration, we first use Algorithm 2 to approximate the CE over 1000
epochs. We use $N=500$ samples to estimate the weighted MSE $L_{\kappa}$ (see
Eq. (37)), and we augment the dataset by 30 times using the data augmentation
method described in Section 3.3. Other settings for training the ANN are kept
identical to those in Section 6.2.1. Given $\bm{\theta}$, we execute Algorithm
3 over 20 epochs using $N=25$ and a single batch. The dataset is also
augmented by 30-fold. In Algorithm 3, the learning rate is gradually reduced
from $0.1$ to $0.02$ to reduce the stochastic effect in the final epochs and
ensure the convergence of the algorithm. Algorithm 3 returns $d_{\kappa+1}$,
which is then utilized in the subsequent iteration.

For each iteration of the main procedure, Algorithms 2 and 3 must evaluate the
observational map 500 times. Additionally, the latter evaluates the gradient
$\nabla_{d}h$ 500 more times. In total, Algorithm 1 executes $20,000$
evaluations of the observational map and $10,000$ evaluations of the gradient
$\nabla_{d}h$. The total number of evaluations of the observational map and
its derivative is of the same order of magnitude as the estimations of the
tECV for a single design (see Fig. 5). Our method gains significant
computational efficiency owing to the nonlocal approximation of the CE and the
applied transfer learning technique. The performance of the algorithm in terms
of the design vector and the tECV is depicted in Fig. 6. Convergence can be
clearly observed after 15 iterations.

(a)

(b)

Figure 6: Convergence in terms of (a) the design parameter and (b) the tECV.
Figure 7: Potential fields solved using the FE model for the initial DOE
(above) and the A-optimal DOE (below).

##### A-optimal DOE

The A-optimal design found by the algorithm closely resembles the solution
$d_{\text{A}}=[1,\;1,\;1,\;-1,-1,\;1,\;1,-1,-1]^{\top}$. The potential fields
solved using the FE model are illustrated in Fig. 7 for both the initial DOE
and the A-optimal DOE. The magnitude of the potential field under the
A-optimal DOE, particularly at the electrodes, is considerably larger than
that of the initial DOE. Consequently, the effect of the measurement error on
the posterior variance is minimized. In this example, minimizing the tECV
leads to the DOE that is identical to the DOE maximizing the EIG, reported in
[18]. This observation suggests that the conditional variance and the
information gain are strongly correlated.

The typical posterior densities obtained from the initial DOE and the
A-optimal DOE are depicted in Fig. 8(a) and (b), respectively. The posterior
mean that is inferred using the A-optimal DOE predicts the ground-truth value
with a negligible error, whereas the posterior variance is significantly
reduced compared with that of the initial DOE.

(a)

(b)

Figure 8: Contour plots of the posterior densities and their mean values for
(a) initial DOE and (b) A-optimal DOE. The ground truth value is
$\eta_{1}=0.748$, $\eta_{2}=-0.848$.

##  7 Conclusion

This study presents an efficient computational method for estimating the tECV
within the A-optimal DOE framework using the PACE approach. The intractability
of the posterior distribution does not affect the computational efficiency of
our approach, as both the approximating and the sampling the posterior
distribution are excluded. Moreover, we derive an asymptotic error estimation
of our method and verify it through numerical experiments.

To address continuous design domains, we combine the PACE framework with
stochastic optimization algorithms to seek the A-optimal DOE. We demonstrate
our method by using the ANN to approximate the CE. The stochastic optimization
algorithms used for optimizing the ANN’s weights and for finding the optimal
design are integrated, as their loss functions are identical. We further
propose a nonlocal approximation of the CE to reduce the number of evaluations
of the observation map, which can be computationally demanding in practice.
Numerical experiments show that our approach requires significantly fewer
evaluations of the observational map compared to the crude IS method.

##  Appendix A Computing tECV using the IS approach

This appendix describes the computation of the ECV using the IS estimator of
CE. This approach requires a double-loop MC simulation. The inner one
estimates the posterior variances using Eq. (4), while the outer one estimates
the expected posterior variance using Eq. (6).

The outer loop for a given design setting $d$ is formulated as

$\widehat{V}^{\text{IS}}=\dfrac{1}{N_{\text{o}}}\sum_{i=1}^{N_{\text{o}}}\sum_{k=1}^{n}\widehat{\operatorname{Var}}(Q\mid
y^{(i)}),\quad\text{with}\quad
y^{(i)}=h(q^{(i)},d)+\xi^{(i)},\;i=1,\dots,N_{\text{o}},$ (48)

where $\\{q^{(i)}\\}_{1}^{N_{\text{o}}}$ and
$\\{\xi^{(i)}\\}_{1}^{N_{\text{o}}}$ are the $N_{\text{o}}$ _i.i.d._ samples
of the RV $Q$ and $\varXi$, respectively. For each sample pair
$(q^{(i)},y^{(i)})$, an inner loop is performed to estimate the posterior
variance as $\widehat{\operatorname{Var}}(Q\mid y^{(i)})$

$\widehat{\operatorname{Var}}\left[Q\mid
y^{(i)}\right]=\dfrac{\sum_{j=1}^{N_{\text{i}}}[q^{(i,j)}]^{\odot
2}\pi_{\varXi}\left(h(q^{(i,j)},d)-y^{(i)}\right)}{\sum_{j=1}^{N_{\text{i}}}\pi_{\varXi}\left(h(q^{(i,j)})-y^{(i)}\right)}-\Biggl{[}\dfrac{\sum_{j=1}^{N_{\text{i}}}q^{(i,j)}\pi_{\varXi}\left(h(q^{(i,j)},d)-y^{(i)}\right)}{\sum_{j=1}^{N_{\text{i}}}\pi_{\varXi}\left(h(q^{(i,j)},d)-y^{(i)}\right)}\Biggr{]}^{\odot
2},$ (49)

where $\\{q^{(i,j)}\\}_{i=1,j=1}^{N_{\text{o}},N_{\text{i}}}$ are
$N_{\text{i}}\times N_{\text{o}}$ _i.i.d._ samples of the RV $Q$. This method
requires $(N_{\text{i}}+1)\times N_{\text{o}}$ evaluations of the
observational map $h$.

##  Appendix B Proofs

###  B.1 Proof of Proposition 1

###### Proof.

By combining the laws of total mean and total variance of the conditional
expectation, _i.e._

$\displaystyle\operatorname{E}\left[\operatorname{E}\left[Q_{i}\mid{Y_{d}}\right]\right]$
$\displaystyle=\operatorname{E}\left[Q_{i}\right],$ (50)
$\displaystyle\operatorname{Var}\left[Q_{i}\right]$
$\displaystyle=\operatorname{E}\left[\operatorname{Var}\left[Q_{i}\mid{Y_{d}}\right]\right]+\operatorname{Var}\left[\operatorname{E}\left[Q_{i}\mid{Y_{d}}\right]\right]$

for $i=1,\dots,n$, we attain

$\displaystyle\operatorname{E}\left[\operatorname{Var}\left[Q_{i}|{Y_{d}}\right]\right]$
$\displaystyle=\operatorname{Var}\left[Q_{i}\right]-\operatorname{Var}\left[\operatorname{E}\left[Q_{i}\mid{Y_{d}}\right]\right]$
(51)
$\displaystyle=\operatorname{E}\left[Q_{i}^{2}-\operatorname{E}\left[Q_{i}\right]^{2}\right]-\operatorname{E}\left[\operatorname{E}\left[Q_{i}\mid{Y_{d}}\right]^{2}-\bigl{(}\operatorname{E}\left[\operatorname{E}\left[Q_{i}\mid{Y_{d}}\right]\right]\bigr{)}{}^{2}\right]$
(52)
$\displaystyle=\operatorname{E}\left[Q_{i}^{2}-\operatorname{E}\left[Q_{i}\mid{Y_{d}}\right]^{2}\right].$
(53)

Hence, the tECV can be formulated as

$\displaystyle V(d)$
$\displaystyle\equiv\sum_{i=1}^{n}\operatorname{E}\left[\operatorname{Var}\left[Q_{i}|{Y_{d}}\right]\right]$
(54)
$\displaystyle=\operatorname{E}\left[\sum_{i=1}^{n}Q_{i}^{2}-\operatorname{E}\left[Q_{i}\mid{Y_{d}}\right]^{2}\right]$
(55) $\displaystyle=\operatorname{E}\left[\left\lVert
Q-\operatorname{E}\left[Q\mid{Y_{d}}\right]\right\rVert_{2}^{2}\right].$ (56)

∎

###  B.2 Proof of Proposition 2

###### Proof.

We structure our proof of Proposition 2 in three steps.

_Step 1: Estimation of the optimization error $\epsilon_{\text{opt}}$._
Following the central limit theorem,
$\sqrt{N}\bigl{(}\widehat{\mathcal{M}}(f^{*}|D_{N})-\mathcal{M}(f^{*})\bigr{)}\rightarrow\mathcal{N}(0,\operatorname{Var}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right])$ as $N\gg 1$. For a RV
$X\sim\mathcal{N}(0,\sigma^{2})$, we have
$\operatorname{E}\left[|X|\right]=\sqrt{2/\pi}\;\sigma$. Consequently, the
statistical error of the MC estimator $\widehat{\mathcal{M}}(f^{*}|D_{N})$ can
be asymptotically estimated as

$\displaystyle\operatorname{E}\left[\left\lvert\widehat{\mathcal{M}}(f^{*}\mid
D_{N})-\mathcal{M}(f^{*})\right\rvert\right]$
$\displaystyle\approx\sqrt{\dfrac{2}{\pi
N}}\operatorname{Var}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right]^{1/2}$ (57)
$\displaystyle=\mathcal{O}\Bigl{(}\dfrac{2}{\sqrt{\pi
N}}\operatorname{E}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right]\Bigr{)},$

where we use Assumption 2 to obtain the last step. Combining Assumption 3 and
the result stated in Eq. (57) yields the following estimation of the
optimization error $\epsilon_{\text{opt}}$

$\epsilon_{\text{opt}}=\mathcal{O}\Bigl{(}\dfrac{2}{\sqrt{\pi
N}}\operatorname{E}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right]\Bigr{)}.$ (58)

Moreover, using the orthogonal property
$\operatorname{E}\left[(Q-f^{*}({Y_{d}}))^{\top}f({Y_{d}}))\right]=0$ for
every $f\in\mathcal{S^{\prime}}$, we obtain

$\displaystyle\epsilon_{\text{opt}}$
$\displaystyle=\operatorname{E}\left[\left\lVert\mathcal{M}(f(\cdot;\bm{\theta}_{D_{N}}))-\mathcal{M}(f^{*})\right|\right]$
(59)
$\displaystyle=\operatorname{E}\left[\Biggl{|}\operatorname{E}\left[\left\lVert
Q-f({Y_{d}};\bm{\theta}_{D_{N}})\right\rVert_{2}^{2}\;\Big{|}\;\bm{\theta}_{D_{N}}\right]-\operatorname{E}\left[\left\lVert
Q-f^{*}({Y_{d}})\right\rVert_{2}^{2}\right]\Biggr{|}\right]$
$\displaystyle=\operatorname{E}\left[\operatorname{E}\left[\left\lVert
f({Y_{d}};\bm{\theta}_{D_{N}})-f^{*}({Y_{d}})\right\rVert_{2}^{2}\;\Big{|}\;\bm{\theta}_{D_{N}}\right]\right]$
$\displaystyle=\operatorname{E}\left[\left\lVert
f({Y_{d}};\bm{\theta}_{D_{N}})-f^{*}({Y_{d}})\right\rVert_{2}^{2}\right].$

By combining Eqs. (58) and (59), we obtain

$\operatorname{E}\left[\left\lVert
f({Y_{d}};\bm{\theta}_{D_{N}})-f^{*}({Y_{d}})\right\rVert_{2}^{2}\right]=\mathcal{O}\Bigl{(}\dfrac{2}{\sqrt{\pi
N}}\operatorname{E}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right]\Bigr{)}.$ (60)

_Step 2: Estimation of the MC estimator error $\epsilon_{\text{MC}}$._
Following the central limit theorem, for fixed hyperparameters
$\bm{\theta}_{D_{N}}$ as $M\gg 1$, we obtain

$\sqrt{M}\left(\widehat{\mathcal{M}}(f(\cdot;\bm{\theta}_{D_{N}})\;\big{|}\;D_{M})-\mathcal{M}\left(f({Y_{d}};\bm{\theta}_{D_{N}}\right)\right)\rightarrow\mathcal{N}(0,\operatorname{Var}\left[\left\lVert
Q-f({Y_{d}};\bm{\theta}_{D_{N}})\right\rVert_{2}^{2}\;\Big{|}\;\bm{\theta}_{D_{N}}\right]^{1/2}).$
(61)

We estimate the error $\epsilon_{\text{MC}}$ as $M\to\infty$ using the central
limit theorem as

$\displaystyle\epsilon_{\text{MC}}$
$\displaystyle\equiv\operatorname{E}\left[\operatorname{E}\left[\bigl{|}\widehat{\mathcal{M}}(f(\cdot;\bm{\theta}_{D_{N}})|D_{M})-\mathcal{M}(f({Y_{d}};\bm{\theta}_{D_{N}}))\bigr{|}\;\Big{|}\;\bm{\theta}_{D_{N}}\right]\right]$
(62) $\displaystyle\approx\operatorname{E}\left[\sqrt{\dfrac{2}{\pi
M}}\operatorname{Var}\left[\left\lVert
Q-f({Y_{d}};\bm{\theta}_{D_{N}})\right\rVert_{2}^{2}\;\Big{|}\;\bm{\theta}_{D_{N}}\right]^{1/2}\right]$
$\displaystyle=\mathcal{O}\Bigl{(}\sqrt{\dfrac{2}{\pi
M}}\operatorname{Var}\left[\left\lVert
Q-f^{*}({Y_{d}})\right\rVert_{2}^{2}\right]^{1/2}\Bigr{)}$
$\displaystyle=\mathcal{O}\left(\dfrac{2}{\sqrt{\pi
M}}\operatorname{E}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right]\right),$

where the third step is obtained owing to $\operatorname{E}\left[\left\lVert
f({Y_{d}};\bm{\theta}_{D_{N}})-f^{*}({Y_{d}})\right\rVert_{2}^{2}\right]=\mathcal{O}(2/\sqrt{\pi
N})$ as stated in Eq. (60).

_Step 3._ Eventually, by combining Eqs. (21), (58), and (62), we can estimate
the error of the estimator $\widehat{V}_{d}(D_{N},D_{M})$ as

$\operatorname{E}\left[\left\lvert\widehat{\mathcal{M}}(f(\cdot;\bm{\theta}_{D_{N}})\;\big{|}\;D_{M})-\mathcal{M}(\phi_{d})\right\rvert\right]=\mathcal{O}\Biggl{(}\Bigl{(}\dfrac{2}{\sqrt{\pi
N}}+\dfrac{2}{\sqrt{\pi M}}\Bigr{)}\operatorname{E}\left[\left\lVert
Q-f^{*}(Y_{d})\right\rVert_{2}^{2}\right]\Biggr{)}+\epsilon_{S^{\prime}}.$
(63)

∎

###  B.3 Orthogonality property of the CE

######  Theorem 1 (Orthogonality property).

Let $Q$ and $Y$ be the finite-variance RVs valued in $\mathbb{R}^{n}$ and
$\mathbb{R}^{m}$, respectively . Let $L_{2}(\sigma_{Y})$ be the collection of
all random variables of type $g(Y)$, where
$g:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ is an arbitrary function
satisfying $\operatorname{E}\left[\lVert g(Y)\rVert_{2}^{2}\right]<\infty$.
Then,

  * A)

$\operatorname{E}\left[\;f(Y){}^{\top}\;\bigl{(}Q-\operatorname{E}\left[Q\mid
Y\right]\bigr{)}\;\right]=0,$ (64)

  * B)

the CE $\operatorname{E}\left[Q\mid Y\right]$ is the RV in $L_{2}(\sigma_{Y})$
that minimizes the MSE, i.e.,

$\displaystyle\operatorname{E}\left[Q\mid Y\right]$ $\displaystyle=\phi(Y),$
(65) $\displaystyle\text{where}\quad\phi$ $\displaystyle=\arg\min_{f(Y)\in
L_{2}(\sigma_{Y})}\operatorname{E}\left[\left\lVert
Q-g(Y)\right\rVert_{2}^{2}\right].$

###### Proof.

For any measurable set $B\subset\mathbb{R}^{m}$, let
$A=Y^{f\;-1}(B)\equiv\\{\omega\;:Y(\omega)\in B\\}$. As $A\in\sigma_{Y}$, we
achieve

$\displaystyle\operatorname{E}\left[\mathbf{1}_{B}(Y)\;\left(Q-\operatorname{E}\left[Q\mid
Y\right]\right)\right]$
$\displaystyle=\operatorname{E}\left[\mathbf{1}_{B}(Y)\;Q\right]-\operatorname{E}\left[\mathbf{1}_{B}(Y)\operatorname{E}\left[Q\mid
Y\right]\right]$ (66)
$\displaystyle=\int_{A}Q(\omega)\mathrm{d}\mathbb{P}(\omega)-\int_{A}\operatorname{E}\left[Q|Y\right](\omega)\mathrm{d}\mathbb{P}(\omega)$
$\displaystyle=0,$

where $\mathbf{1}_{B}(Y(\omega))=1$ if $Y(\omega)\in B$ and $0$ otherwise.

Let $X^{n}$ be a $\sigma_{Y}$-measurable simple RV that converges in terms of
distribution to $g\circ Y$. Using the result from Eq. (66), we prove part A)
ot Theorem 1.

Let $Z=\operatorname{E}\left[Q\mid Y\right]-\phi(Y)$, we obtain

$\displaystyle\operatorname{E}\left[\bigl{\lVert}Q-\phi(Y)\bigr{\rVert}_{2}^{2}\right]$
$\displaystyle=\operatorname{E}\left[\bigl{\lVert}Q-\operatorname{E}{\bigl{[}Q\mid
Y\bigr{]}}+Z\bigr{\lVert}^{2}\right]$ (67)
$\displaystyle=\operatorname{E}\left[\bigl{\lVert}Q-\operatorname{E}{\bigl{[}Q\mid
Y\bigr{]}}\bigr{\rVert}_{2}^{2}\right]+\operatorname{E}\left[\bigl{\lVert}Z\bigr{\rVert}_{2}^{2}\right],$

as the cross-product term vanishes. Eq. (67) directly shows that
$\operatorname{E}\left[\bigl{\lVert}Q-\phi(Y)\bigr{\rVert}_{2}^{2}\right]$ is
minimized when $Z$ is the $n$-dimensional vector of zeros, or
$\operatorname{E}\left[Q\mid Y\right]=\phi(Y)$ (part B). ∎

##  Appendix C Linear approximation of the CE

###  C.1 Analytical formulation of the linear approximation

The linear approximation of the CE $f_{\text{l}}(Y_{d})=\bm{A}Y_{d}+b$ is the
solution to the least mean-squares problem

$\bm{A},b=\quad\arg\min_{\bm{A}^{\prime}\in\mathbb{R}^{n\times
m},b^{\prime}\in\mathbb{R}^{n}}\operatorname{E}\left[\bigl{\lVert}Q-\bm{A}^{\prime}Y_{d}-b^{\prime}\bigr{\rVert}_{2}^{2}\right].$
(68)

Using the first order necessary conditions

$\displaystyle\operatorname{E}\left[Q-\bm{A}Y_{d}-b\right]$
$\displaystyle=0_{n},$ (69)
$\displaystyle\operatorname{E}\left[(Q-\bm{A}Y_{d}-b)Y_{d}^{\top}\right]$
$\displaystyle=\bm{0}_{n\times m},$

where $\bm{0}_{n\times m}$ denotes the $n\times m$ matrix of zeros, we obtain

$\displaystyle\bm{A}$
$\displaystyle=\operatorname{E}\left[\bigl{(}Q-\operatorname{E}[Q]\bigr{)}Y_{d}^{\top}\right]\operatorname{E}\left[\bigl{(}Y-\operatorname{E}{[Y_{d}]}\bigr{)}Y^{\top}\right]^{-1}$
(70)
$\displaystyle=\operatorname{E}\left[\bigl{(}Q-\operatorname{E}{[Q]}\bigr{)}{\bigl{(}Y_{d}-\operatorname{E}{[Y_{d}]}\bigr{)}^{\top}}\right]\operatorname{E}\left[\bigl{(}Y_{d}-\operatorname{E}{[Y_{d}]}\bigr{)}{\bigl{(}Y-\operatorname{E}{[Y_{d}]}\bigr{)}^{\top}}\right]^{-1}$
$\displaystyle=\textnormal{Cov}\left[Q,Y_{d}\right]\textnormal{Cov}\left[Y_{d}\right]^{-1},$
$\displaystyle b$ $\displaystyle=\operatorname{E}\left[Q-\bm{A}Y_{d}\right].$

###  C.2 Empirical linear approximation of the CE

Given dataset $D_{N}={(q^{(i)},y^{(i)})}_{i=1}^{N}$ of the $N$ _i.i.d._
samples of the RV pair $(Q,Y)$, the empirical linear approximation of the CE
is obtained as

$\widehat{\phi}(y)=\widehat{\operatorname{Cov}}[{Q,Y_{d}}]\;\bigl{[}\widehat{\operatorname{Cov}}[{Y_{d}}]\bigr{]}^{-1}y+\widehat{b}$
(71)

where

$\displaystyle\widehat{\operatorname{Cov}}[{Q,Y_{d}}]=\dfrac{1}{N}\sum_{i=1}^{N}[q^{(i)}-\overline{q}]\,[y^{(i)}-\overline{y}]^{\top},$
(72)
$\displaystyle\widehat{\operatorname{Cov}}[{Y_{d},Y_{d}}]=\dfrac{1}{N}\sum_{i=1}^{N}[y^{(i)}-\overline{y}]\,[y^{(i)}-\overline{y}]^{\top},$
$\displaystyle\widehat{b}=\overline{q}-\dfrac{1}{N}\sum_{i=1}^{N}\widehat{\operatorname{Cov}}[{Q,Y_{d}}]\;\bigl{[}\widehat{\operatorname{Cov}}[{Y_{d}}]\bigr{]}^{-1}y^{(i)}.$

Here, $\overline{q}$ and $\overline{y}$ are empirical means of the RVs $Q$ and
$Y_{d}$, respectively, computed as

$\overline{q}=\dfrac{1}{N}\sum_{i=1}^{N}q^{(i)},\quad\overline{y}=\dfrac{1}{N}\sum_{i=1}^{N}y^{(i)}.$
(73)

The reduced variance estimators corresponding to those in Eq. (72) are
straightforwardly obtained by replacing the role of the dataset $D_{N}$ with
its augmented version $D_{N}^{\textnormal{vr}}$ (see (Eq. 28)).

##  Appendix D Reduced variance estimators of MSE and tECV

We show that the estimator $\widehat{\mathcal{M}}^{\text{vr}}$ (see Eq. (30))
provides an estimation with reduced variance compared with the crude MC
estimator $\widehat{\mathcal{M}}$ (see Eq. (16)). For a given vector $d$, we
have

$\displaystyle\lim_{a\rightarrow\infty}\widehat{\mathcal{M}}^{\text{vr}}(f\mid
D_{N})$
$\displaystyle\;=_{\text{a.s.}}\;\dfrac{1}{N}\sum_{i=1,\dots,N}\operatorname{E}\Bigl{[}\bigl{\lVert}q^{(i)}-f\bigl{(}h(q^{(i)},d)+\varXi\bigr{)}\rVert_{2}^{2}\Bigr{]}$
(74)
$\displaystyle=\;\dfrac{1}{N}\sum_{i=1,\dots,N}\operatorname{E}\left[A(Q,\varXi)|Q=q^{(i)}\right],$

$\mathbb{P}$-almost surely, where

$A(q,\xi)\equiv\bigl{\lVert}q-f\bigl{(}h(q,d)+\xi\bigr{)}\bigr{\rVert}_{2}^{2}.$
(75)

Let $\widehat{\mathcal{M}}^{\text{vr}*}(f)$ denote the right-hand-side term in
Eq. (74). Because $\\{q^{(i)}\\}_{i=1}^{N}$ are _i.i.d._ samples, we
approximately quantify the statistical errors of the estimators
$\widehat{\mathcal{M}}$ and $\widehat{\mathcal{M}}^{\text{vr}*}(f)$,
respectively, as

$\displaystyle\operatorname{Var}\left[\widehat{\mathcal{M}}(f|D_{N})-\operatorname{E}\left[A(Q,\varXi)\right]\right]$
$\displaystyle\approx\dfrac{\operatorname{Var}\left[A(Q,\varXi)\right]}{N},$
(76a)
$\displaystyle\operatorname{Var}\left[\widehat{\mathcal{M}}^{\text{vr}*}(f)-\operatorname{E}\left[A(Q,\varXi)\right]\right]$
$\displaystyle\approx\dfrac{\operatorname{Var}\left[\operatorname{E}\left[A(Q,\varXi)\mid
Q\right]\right]}{N}.$ (76b)

Using the law of total variance,

$\operatorname{Var}\left[A(Q,\varXi)\right]=\operatorname{Var}\left[\operatorname{E}\left[A(Q,\varXi)\mid
Q\right]\right]+\operatorname{E}\left[\operatorname{Var}\left[A(Q,\varXi)\mid
Q\right]\right],$ (77)

we obtain

$\operatorname{Var}\left[\operatorname{E}\left[A(Q,\varXi)\mid
Q\right]\right]\leq\operatorname{Var}\left[A(Q,\varXi)\right].$ (78)

Therefore, for $a\gg 1$,
$\operatorname{Var}\left[\widehat{\mathcal{M}}^{\text{vr}}(f\mid
D_{N})-\operatorname{E}\left[A(Q,\varXi)\right]\right]\leq\operatorname{Var}\left[\widehat{\mathcal{M}}(f\mid
D_{N})-\operatorname{E}\left[A(Q,\varXi)\right]\right]$.

##  Appendix E Gradient of the weighted MSE w.r.t. design parameters

The gradient $\nabla_{d}\Bigl{[}\left\lVert
q-f_{\text{N}}(h(q)+\xi,d)\right\rVert_{2}^{2}\Bigr{]}$ is expanded using the
chain rule as

$\displaystyle\nabla_{d}\Bigl{[}\left\lVert
q-f_{\text{N}}(h(q,d)+\xi,d)\right\rVert_{2}^{2}\Bigr{]}=-2$
$\displaystyle\Bigl{[}\nabla_{d}\,f_{\text{N}}(y,d)+\nabla_{y}\,f_{\text{N}}(y,d)\,\nabla_{d}\,h(q,d)\Bigl{]}^{\top}$
(79) $\displaystyle\Bigl{[}q-f_{\text{N}}(h(q,d)+\xi,d)\Bigr{]},$

where $y=h(q,d)+\xi$. Notably, for ANN $f_{\text{N}}$, the Jacobian matrices
$\nabla_{d}f_{\text{N}}$ and $\nabla_{y}f_{\text{N}}$ can be directly obtained
numerically using the backpropagation [27]. We finally obtain the value of
$\nabla_{d}\Bigl{[}\left\lVert
q-f_{\text{N}}(h(q)+\xi,d)\right\rVert_{2}^{2}\Bigr{]}(q^{(i)},y^{(i,j)},d_{0})$
as

$\displaystyle\nabla_{d}\Bigl{[}\left\lVert
q-f_{\text{N}}(h(q,d)+\xi,d)\right\rVert_{2}^{2}\Bigr{]}$
$\displaystyle(q^{(i)},y^{(i,j)},d_{0})=$ (80) $\displaystyle-2$
$\displaystyle\Bigl{[}\nabla_{d}\,f_{\text{N}}(y^{(i,j)},d_{0})+\phantom{\Bigl{(}}\nabla_{y}\,f_{\text{N}}(y^{(i,j)},d_{0})\,\nabla_{d}\,h(q^{(i)},d_{0})\Bigr{]}^{\top}$
$\displaystyle\Bigl{[}q^{(i)}-f_{\text{N}}(y^{(i,j)},d_{0})\Bigr{]}.$

##  Appendix F  Derivation of the finite element formulation

We define $\mathcal{H}:=H^{1}({B}\times\mathbb{R}^{N_{\rm{el}}})$ as the space
of the solution for the potential field $(u(\omega),{U}(\omega))$ for
$\omega\in\Omega$ and the Bochner space
$L_{\mathbb{P}}^{2}(\Omega;\mathcal{H})$ as

$L_{\mathbb{P}}^{2}(\Omega;\mathcal{H}):=\left\\{(u,{U}):\Omega\to\mathcal{H}\quad\text{s.t.}\,\int_{\Omega}\lVert(u(\omega),{U}(\omega))\rVert_{\mathcal{H}}^{2}\mathrm{d}{}\mathbb{P}(\omega)<\infty\right\\},$
(81)

where $U(\omega)=[U_{1}(\omega),\dots,U_{N_{\rm{el}}}(\omega)]^{\top}$. For
the bilinear form $\mathcal{B}:\mathcal{H}\times\mathcal{H}\to\mathbb{R}$,
which is given as

$\mathcal{B}((u,{U}),(v,{V})):=\int_{\mathcal{D}}\sigma\,\nabla u\cdot\nabla
v\mathrm{d}{}{B}+\sum_{l=1}^{N_{\rm{el}}}\frac{1}{z_{l}}\int_{E_{l}}({U}_{l}-u)({V}_{l}-v)\mathrm{d}{}E_{l},$
(82)

we aim to determine $(u,{U})\in L_{\mathbb{P}}^{2}(\Omega;\mathcal{H})$ such
that the weak formulation

$\mathcal{B}\left((u(\omega),{U}(\omega)),(v,{V})\right)=\sum_{l=1}^{N_{\rm{el}}}I_{l}\,U_{l}(\omega)$
(83)

is fulfilled for all $(v,{V})\in L_{\mathbb{P}}^{2}(\Omega;\mathcal{H})$
almost surely.

## References

  * [1] Turab Lookman, Prasanna V Balachandran, Dezhen Xue, and Ruihao Yuan.  Active learning in materials science with emphasis on adaptive sampling using uncertainties for targeted design.  npj Computational Materials, 5(1):21, 2019.
  * [2] Laura Ilzarbe, María Jesús Álvarez, Elisabeth Viles, and Martín Tanco.  Practical applications of design of experiments in the field of engineering: a bibliographical review.  Quality and Reliability Engineering International, 24(4):417–428, 2008.
  * [3] Stavros N. Politis, Paolo Colombo, Gaia Colombo, and Dimitrios M. Rekkas.  Design of experiments (doe) in pharmaceutical development.  Drug development and industrial pharmacy, 43(6):889–901, 2017.
  * [4] Bhupinder Singh, Rajiv Kumar, and Naveen Ahuja.  Optimizing drug delivery systems using systematic “design of experiments”. Part I: fundamental aspects.  Critical Reviews™ in Therapeutic Drug Carrier Systems, 22(1), 2005\.
  * [5] Roger Mead.  Statistical methods in agriculture and experimental biology.  CRC Press, 2017.
  * [6] James H Bullard, Elizabeth Purdom, Kasper D Hansen, and Sandrine Dudoit.  Evaluation of statistical methods for normalization and differential expression in mrna-seq experiments.  BMC bioinformatics, 11(1):1–13, 2010.
  * [7] Kathryn Chaloner and Isabella Verdinelli.  Bayesian experimental design: A review.  Statistical science, pages 273–304, 1995.
  * [8] Xun Huan and Youssef M. Marzouk.  Simulation-based optimal Bayesian experimental design for nonlinear systems.  Journal of Computational Physics, 232(1):288–317, 2013.
  * [9] Alen Alexanderian, Noemi Petra, Georg Stadler, and Omar Ghattas.  A-optimal design of experiments for infinite-dimensional Bayesian linear inverse problems with regularized $l_{0}$-sparsification.  SIAM Journal on Scientific Computing, 36(5):A2122–A2148, 2014.
  * [10] Elizabeth G. Ryan, Christopher C. Drovandi, and Anthony N. Pettitt.  Simulation-based fully Bayesian experimental design for mixed effects models.  Computational Statistics & Data Analysis, 92:26–39, 2015.
  * [11] Håvard Rue, Sara Martino, and Nicolas Chopin.  Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(2):319–392, 2009.
  * [12] Elizabeth G. Ryan, Christopher C. Drovandi, James M. McGree, and Anthony N. Pettitt.  A review of modern computational algorithms for bayesian optimal design.  International Statistical Review, 84(1):128–154, 2016.
  * [13] Arved Bartuska, Luis Espath, and Raúl Tempone.  Small-noise approximation for Bayesian optimal experimental design with nuisance uncertainty.  Computer Methods in Applied Mechanics and Engineering, 399:115320, 2022.
  * [14] Arved Bartuska, André Gustavo Carlon, Luis Espath, Sebastian Krumscheid, and Raúl Tempone.  Double-loop quasi-monte carlo estimator for nested integration.  arXiv preprint arXiv:2302.14119, 2023.
  * [15] Alen Alexanderian, Noemi Petra, Georg Stadler, and Omar Ghattas.  A fast and scalable method for a-optimal design of experiments for infinite-dimensional Bayesian nonlinear inverse problems.  SIAM Journal on Scientific Computing, 38(1):A243–A272, 2016.
  * [16] Joakim Beck, Ben Mansour Dia, Luis F.R. Espath, Quan Long, and Raúl Tempone.  Fast Bayesian experimental design: Laplace-based importance sampling for the expected information gain.  Computer Methods in Applied Mechanics and Engineering, 334:523–553, 2018.
  * [17] Joakim Beck, Ben Mansour Dia, Luis Espath, and Raúl Tempone.  Multilevel double loop monte carlo and stochastic collocation methods with importance sampling for Bayesian optimal experimental design.  International Journal for Numerical Methods in Engineering, 121(15):3482–3503, 2020.
  * [18] André Gustavo Carlon, Ben Mansour Dia, Luis Espath, Rafael Holdorf Lopez, and Raúl Tempone.  Nesterov-aided stochastic gradient methods using Laplace approximation for Bayesian design optimization.  Computer Methods in Applied Mechanics and Engineering, 363:112909, 2020.
  * [19] Truong-Vinh Hoang, Sebastian Krumscheid, Hermann G. Matthies, and Raúl Tempone.  Machine learning-based conditional mean filter: A generalization of the ensemble kalman filter for nonlinear data assimilation.  Foundations of Data Science, 0:–, 2022.
  * [20] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  MIT press, 2016.
  * [21] Frank Rosenblatt.  The perceptron: a probabilistic model for information storage and organization in the brain.  Psychological review, 65(6):386, 1958.
  * [22] Adam Bobrowski.  Functional analysis for Probability and Stochastic Process.  Cambridge University press, 2005.
  * [23] Rick Durrett.  Probability: theory and examples.  Cambridge University press, 2019.
  * [24] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv:1412.6980, 2014.
  * [25] Lisa Torrey and Jude Shavlik.  Transfer learning.  In Handbook of research on machine learning applications and trends: algorithms, methods, and techniques, pages 242–264. IGI global, 2010\.
  * [26] Karl Weiss, Taghi M Khoshgoftaar, and DingDing Wang.  A survey of transfer learning.  Journal of Big data, 3(1):1–40, 2016.
  * [27] Andreas Griewank and Andrea Walther.  Evaluating Derivatives.  Society for Industrial and Applied Mathematics, second edition, 2008.
  * [28] Erkki Somersalo, Margaret Cheney, and David Isaacson.  Existence and uniqueness for electrode models for electric current computed tomography.  SIAM Journal on Applied Mathematics, 52(4):1023–1040, 1992.