Bayesian Adaptive Calibration and Optimal Design
Authors: Rafael Oliveira
Dino Sejdinovic
David Howard
Edwin V. Bonilla
Venue: NeurIPS 2024 poster (2024)
Keywords: Gaussian processes
Bayesian inference
variational inference
experimental design
active learning
calibration
Primary area: active_learning
TL;DR: We propose a method which jointly estimates posteriors and informative designs to calibrate computer simulations of physical processes.

Abstract: The process of calibrating computer models of natural phenomena is essential for applications in the physical sciences, where plenty of domain knowledge can be embedded into simulations and then calibrated against real observations. Current machine learning approaches, however, mostly rely on rerunning simulations over a fixed set of designs available in the observed data, potentially neglecting informative correlations across the design space and requiring a large amount of simulations. Instead, we consider the calibration process from the perspective of Bayesian adaptive experimental design and propose a data-efficient algorithm to run maximally informative simulations within a batch-sequential process. At each round, the algorithm jointly estimates the parameters posterior distribution and optimal designs by maximising a variational lower bound of the expected information gain. The simulator is modelled as a sample from a Gaussian process, which allows us to correlate simulations and real data with the unknown calibration parameters. We show the benefits of our method when compared to related approaches across synthetic and real-data problems.

Bayesian Adaptive Calibration and Optimal Design
Rafael Oliveira" CSIRO's Data61
Dino Sejdinovic University of Adelaide Adelaide, Australia
David Howard CSIRO's Data61 Brisbane, Australia
Edwin V. Bonilla CSIRO's Data61 Sydney, Australia
Sydney, Australia
Abstract
The process of calibrating computer models of natural phenomena is essential for applications in the physical sciences, where plenty of domain knowledge can be embedded into simulations and then calibrated against real observations. Current machine learning approaches, however, mostly rely on rerunning simulations over fixed set of designs available in the observed data, potentially neglecting informative correlations across the design space and requiring large amount simulations. Instead, we consider the calibration process from the perspective of Bayesian adaptive experimental design and propose a data-efficient algorithm to run maximally informative simulations within batch-sequential process. At each round, the algorithm jointly estimates the parameters of the posterior distribution and optimal designs by maximising variational lower bound of the expected information gain. The simulator is modelled as a sample from a Gaussian process, which allows to correlate simulations and observed data with the unknown calibration parameters. We show the benefits of our method when compared to related approaches across synthetic and real-data problems.
Introduction
In many scientific and engineering disciplines, computer simulation models form an essential part of the process of predicting and reasoning about complex phenomena, especially when real data scarce. These simulation models depend on the inputs set by the user, commonly referred to as designs, and on number of parameters representing unknown physical quantities, known as calibration parameters. The problem of setting these parameters so as to closely match observations of the real phenomenon is known the calibration of computer models [1].
The seminal work by Kennedy and O' Haaaan [l] introduces the Bayesian framework for calibration of simulation models, using Gaussian processes [2] to account for the differences between the model and reality, as well as for uncertainty the calibration parameters. While the simulator is an essential tool when obtaining real data is expensive or unfeasible, each run of simulator may itself involve significant computational resources, especially in applications such climate science or complex engineering systems. In this situation, it is imperative to run simulations at carefully chosen settings of designs as well of calibration inputs, using current knowledge to optimise resource use [3-5].
In this contribution, we bridge Bayesian calibration with adaptive experimental design [6] and use information-theone
criteria [7] to guide the selection of simulation settings so that they are most informative about the true value of the calibration parameters. We refer to our approach as BACON (Bayesian Adaptive Calibration and Optimal desigN). BACON allows computational resources to be focused on simulations that provide the most value in terms of reducing epistemic uncertainty. Importantly, in contrast prior work, it optimises designs jointly with calibration inputs in order to capture informative correlations across both spaces. Experimental results on synthetic experiments and robotie gripper design problem demonstrate the benefits of BACON compared to competitive
Conresponding author; rafael. www.aattaadeelliiiiraaaaaaaaaeed cairo.au
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

baselines in terms of computational savings and the quality of the estimated posterior under similar computational constraints.
2 Problem formulation
Let X' - Y represent mapping of experimental designs X to the outcomes of a physical process f(x) E y c R. We are given a set of observed outcomes Y = (3/....., yR] and their associated designs XR {x C X. Observations are coruppted by noise y/ f(x)) + 14+ where V4 ~ (V, o.) is zero-mean Gaussian noise, for E (1,..., R). addition, we have access to the output of computer model h X' @ IK given design input and simulation parameters. Given an optimal setting for the calibration parameters 0 E O, the simulator 1((x, 0²), can be used approximate the outcomes of the real physical process f(x). However, 0" is unknown, and evaluations of the simulator h are costly, though cheaper than executing real experiments evaluating f. Our task is to optimally estimate O' given the real data y outputs of the simulator h and prior distribution p(0"), representing initial assumptions about 0
More concretely, let g := [4(8,, 0,1 represent simulated outcomes for a set of designs Rs := { ) C X' and simulation parameters θ {0. O. Given the cost of rumning simulations, we will associate the simulator with latent function (usually referred to as emulator) drawn from Gaussian process (GP) prior and assume simulation outputs and real data follow joint probability distribution p(y R 8 )")
In this setting, the Bayesian experimental design objective is to propose sequence of simulations which will maximise the expected information gain (EIG) about 0":
(Is θs) = H(p(0" (RR)
E p H
=Ep(s Ng. Rs. [KL "s
= =1(0":9 $5 | YR.
6s),
where HI(-) represents the entropy of probability distribution, (-I(- denotes the Kullback-Leiblee divergence, and 1(0"; S.s y') is the mutual information between 0" and the simulator output ý.s given the real observations y and the simulator inputs to be optimized. We note here that, in our setting, the real observations y a always fixed. Therefore, intuitively, the EIG above captures the reduction in uncertainty that will be obtained when selecting (A's, θs) averaged over all the possible outcomes y.
3 Related work
Our work consists of deriving a Bayesian adaptive experimental design (BAED) approach to the problem of calibration. Therefore, in the following, will briefly discuss current literature these two main research areas.
3.1 Adaptive experimental design
The problem of experimental design has long history [8], spanning from classical fixed design patterns to modern adaptive approaches [9]. Optimal experimental design consists of selecting experiments which will maximise some form of criterion involving measure of utility of the experiment and its associated costs [10]. Under the Bayesian formulation, uncertainty in the outcomes of the process is considered, and the optimality of design is measured in terms of its expected utility [11]. Information theory then allows us to quantify information gain as a utility function, which is commonly applied in modern approaches to Bayesian experimental design [12].
The estimation of posterior distributions becomes computational bottleneck for informatiionthhortticc Bayesian frameworks. Recent work has focused on addressing the difficulties in estimating the expected information gain by means of, e.g., variational inference [13], density-ratio estimation [14], importance sampling [15], and the learning of efficient policies to propose designs [16, 17]. These methods, however, usually assume that the simulator known and inexpensive to evaluate. In contrast, the simulations themselves are modelled as expensive experiments for us, and apply

---

Gaussian process models as emulators to capture uncertainty over the black box simulator. In addition, traditional BAED approaches assume that the prior is trivial to sample from and evaluate densities of, while in our case the starting prior is p(0")y g), which is likely non-trivial. We refer the reader to the recent review on modern Bayesian methods for experimental design by Rainforth et al. [18] for further details on BAED.
3.2 Active learning for calibration
Experimental design approaches generally aim towards the selection of designs for physical ex- periments, whereas we are concerned with the problem of running optimal simulated experiments for model calibration in the presence of real data. When simulations are resource-intensive, a few methods have been derived based on the Bayesian calibration framework proposed by Kennedy and O' Hagan [1]. Busby and Feraille [19] present an algorithm learn GP emulators for simulator which can then be combined with Bayesian inference algorithms, such as Markov chain Monte Carlo [20], to provide a posterior distribution over parameters, In their approach, the optimised variables are solely the calibration parameters, and the selection criterion is based on minimising the integrated mean-square error of the GP predictions. Many other approaches can be applied to this setting by modelling the simulator or its associated likelihood function as a GP, including Bayesian optimisation [3, 21, 22] and methods for adaptive Bayesian quadrature [23, 24). Besides GPs, other algorithms focusing on the selection of calibration parameters have been derived using ensembles of neural networks [25] and deep reinforcement learning [26]. These frameworks, however, do not allow for the selection of simulation design points, usually keeping them co-located with the real data.
Allowing for design point decisions to be included, Leatherman al. [4] presented approaches for combined simulation and physical experimental design following geometric and prediction-error based criteria, though using an offlin,, -on-euuential framework. More recently, Marmin and Filippone [5] derived a deep Gaussian process [27] framework for Bayesian calibration problems and discussed an application with experimental design among other examples. Their experimental design approach to calibration was based on choosing simulations that maximally reduce the variational posterior variance over the calibration parameters, as measured by the derivatives of the evidence lower bound with respect to (w.rtt) variance parameters. In contrast, we aim to directly maximise the information gain w.r.t. the unknown calibration parameters.
4 Gaussian processes for Bayesian calibration
To estimate information gain, we need a probabilistic model which can correlate simulations with real data and the unknown parameters 0* Ideally, the model needs allow for a computationally tractable conditioning on the parameters 0* and account for the discrepancy between real and simulated data. Hence, we follow the Bayesian calibration approach in Kennedy and O' O'Hagan [1] and model:
fxx = p 0*) + =xx), x 0 p(0)),
where E : X' E represents the error (or discrepancy) between simulations and real outcomes, and E IR accounts for possible differences in scale. We place Gaussian process priors on the simulator
h ~ (O, it) and on the error function E ke).
4.1 Bi-fidelity exact Gaussian process model
Since both h and are GPs, simulations and real outcomes can be jointly modelled single Gaussian process. In fact, both the simulator h and the true function f can be seen as different levels of fidelity of the same underlying process, with h representing coarser version of f. Namely, let s E S = {0, 1} denote a fidelity parameter. The combined model is then given by:
[1(xx0) 0), f(x,0, 0, s) (ph(x, 0) + Exx).
8 0
s s=1.
such that
0* 1) and 1(x 0) = f(x, 0,0), 0), for any x, & E X and Ā E O. As a result, for arbitrary points in the joint space 2 E 2 := X' x θ x S. the following covariance function parameterises the combined GP model f GP(O, k):
k(2,, x') == kp(8, 8/ ( 0). (x,, 0')) + s8 kx
3

---

where kp(s, := (1 ((p- (p 1)), z := (x, and z' := (x'.6'ss'). s'). Therefore, any set of real and simulated evaluations are joint normally distributed under combined GP model.
4.2 Joint probabilistic model and predictions
Let Z == ZR(0") == (xx, 0" .1)] represent the set of partially observed inputs for real data YR. and let Zg := [(&,, 0, o)] denote the current set of simulation imputs for the observations s- Under the GP prior, the joint probability model p(S s. Y 0²) can be decomposed as:
s.V. R " 1000)'('''>))'000000))))))'))1 1 pl9 po(sslDpis Dry alf, 0r ("(((("))((" df,
p(Ss.YR: 0-) p(Ss.
(5)
where f. (Z))))) E I R+5 and Z(0") == (Z=²(0), 2 corresponds the full set of inputs. The GP prior then allows us to model real and simulated outcomes jointly as Gaussian random vector f:
f ~ ((0)))
(6)
where K(0") k (Z)())) Z(0)) = (kz,z') e denotes the prior covariance matrix. As- suming Gaussian noise model for the observations y [(x, 0*) + +x(x + v, with (0 ož), the marginal distribution over the observations y = [y RR. 3]] is available closed form as:
pl s. y = (y; 0, K(0") + Ey),
where By denotes the covariance matrix of the observation noise, i.e., Eyla = of for any z, with 8: = 1, and [Eybs = elsewhere.
Under the GP assumptions, we can make predictions about y 1(8, 0) at any pair of 8, Ô E X x O. Conditioning on 0 and dataset D, == {xRR y a: &, S. J.J. let Z, (0*) == {Zx(0*) Z,} denote the set of inputs up to time & conditional on 0" and y, the corresponding outputs. then have that:
p 0 D.) = (( (2
(8)
for 2 = (x. 0), where:
(e(2; 0¹) = k2,2;:0 0-) (K,(0)) + Ey,)
(9)
k. (,, ; 0*) 2(2,2)) kz 0-) "(K((0)) -z,';
(10)
of(2;07 == 1x 2; ))
(11)
with (2 0") (0", and K, (0')) := k (0") (0")). We next describe how apply this model to derive Bayesian adaptive calibration algorithm.
Bayesian adaptive calibration and optimal design
In this section, we describe an approach design experiments for calibration of computer models that incorporates information gathered during the experiments iteratively. We refer to these types of designs adaptive. Thus, we consider the sequential design of experiments setting, where each iteration E N, we optimise:
EIG, 0) = T(0"; : x 3,0,D.-1)
p DD (p 0 10, 0, D.-1))]
(12)
log p00") |1. x, θ D.
EE 1 (8.0.D...)
p(0 (D...)
given the dataset D,-1 == (Xa. y &o ..11, of observations. Given that the expected information gain is submodular [28]. sequential approach allows us to get close enough (usually factor of least - 1/e [29]) to the optimal EIG over the whole experiment, while also allowing algorithmic decisions to adapt to current estimates for ((0((D,)
In practice add small nugger term to the diagonal of the noise covariance matrix for numerical stability.

---

In general, computing the full EIG objective (1), or its sequential version (12), is intractable, as that requires estimating the true posterior and its density conditioned on sampled data. Note that both p(0"ló 1ú, &, D1 and p(g. 0- st. 0 D:-1) depend on the posterior p(0 ( 1D
p(0 (ū. & 0, D. D) p(ù, 0 18, Ó, D2-1) p(0 & 0, D:-1)
(13)
p(ù, 0 8 0, pp-j)- 8, 0, Dr (p00))p(01PP-1
(14)
where the conditional predictive density p(ğ)0". x 0, D.-1) Gaussian and available in closed form (Eq. 8). Clearly, in general, the true posterior is intractable, since p(0" (D,) = pD 0 1p(0") and
p(2))
p(P)) la ((), )p(0)) -00)) d0* involves integration over the entire parameter space θ which can be high dimensional and involve highly non-linear operations, such as computing inverse covariances. In addition, the marginal predictive density p(u]&, 0, Dz-1) Se P(. 0- 02 D) d0 is also
usually intractable for the same reasons.
5.1 Variational EIG lower bound
Following Foster et al. [13], we replace the EIG by a variational objective which does not require the true posterior density over 0" This formulation allows jointly estimate approximation the posterior and select optimal design points & and simulation parameters Ó. Applying the variational lower bound by Barber and Agakov [30] to Eq. 12 yields the following alternative to the EIG:
g(0") 9. 8, x) log p(0-1D--1)
(&G. Ō == E 100.0018.00000
S (X,
(15)
where q(0" g. S, Q) is any conditional probability density model. The gap is given by the expected Kullback-Leiberr (KL) divergence between the true and the variational posterior [13, Sec. A.1]:²
EIG θ) EIG,(&, 0,g) = =E 2.0.D.......... (Dx. 12,- 9) (00 (6)) 2 0. (16)
Maximising the variational EIG lower bound w.r.t. the variational distribution g then provides us with approximation to p(0" jū, &, Ô, D,-1). Therefore, we can simultancously obtain maximally informative designs and optimal variational posteriors by jointly optimising the EIG lower bound w.r.t. the simulator inputs and the variational distribution as:
& θt, Qt E argmax
FIG, 0,9) argmax "... [log y (0)))
(17)
REX.0E0.0E
*x REx.0E0,0E0
given a suitable variational family Q of conditional distributions. Note that, in this formulation, only need samples from the posterior (001D(,-) to estimate the expectation above, which can be approximated via Monte Carlo, without requiring densities other than that of the variational model q.
5.2 Algorithm
Algorithm 1 summarises the method we propose, which we name Bayesian Adaptive Calibration and Optimal desigN (B (BACON). The algorithm starts with an initial dataset Do containing the real data (and possibly previously available simulation data) and an estimate of the posterior given the initial data p(0" (Do). Posterior estimates in BACON can be represented by samples obtained via Markov chain Monte Carlo (MCMC) or variational inference over the GP model and the currently available data D, Note that we only need samples from the previous posterior to estimate the expectation in Eq. 17, with no need to directly evaluate its probability densities. Each iteration starts by optimising the variational EIG lower bound using the objective in Eq. 17 to jointly select an optimal design & simulation parameters Ô, and variational posterior dt- Given the new design &, we run the simulation with the chosen parameters θ,, observing the GP model are then updated with the new data, potentially including a re-estimation of the GP hyper-parameters via, for example, maximum likelihood estimation. The process then repeats given the updated GP and posterior for up to a given number of iterations T. At the end, a final posterior pr (0") = p(0" y R- Vr) and conditional density model qr obtained.
outcome 9t- The calibration posterior Pr pr(0") and
'W will times write (00 (j) to denote (0" g &; 0) to avoid notation clutter, the dependence the inputs (&, 0) remains implicit through the conditioning on g.
5

---

Algorithm BACON
input Do == R. nn::
[Real data}
input po(0") = p00 (Do)
(MCMC or prior distribution
for 1
, ] do
& 01 argmax 2.0 E P P1
[log (0" [)]
{Optimise EIG,)
Ét h(;;; 0.) D. = 1 1
(Run simulation
(Update GP model)
(( p(0 0 (D--1)
(Update posterior via MCMC or VI)
end for
output pT")
[Final posterior}
5.3 Variational posteriors
Any conditional probability density model (0))))) estimating probability densities over the parameter space O given an observation 1 could suit our method. In the following, we describe two possible parameterisations for this model. The first facilitates marginalising latent inputs in GP regression [31, 32], while the second better captures multi-modality in the posterior.
Conditional Gaussian models. Assuming construct a variational conditional density model as:
approximate p(0" (D:) as a Gaussian, we can
(0" jj. 0) == A'"0": mo (Ú. x2 0), (0,
(18)
where my and Eg are given by parametric models. such as neural networks, with parameters Φ. To ensure ΣΦ (-) positive-definite. it can be parameterised by its Cholesky decomposition Ex(- - ( (.),, where Lø (-) is a bower-triangular matrix with positive diagonal entries.
Conditional normalising flows Normalising flows [33] apply the change- of ofee--o--vriira variable formula to derive composable, invertible transformations Sw of fixed base distribution po:
(E)) g(k) o... gw (Eo), Eo- Po
(19)
The log-probability density of point & gw (Eo) under this model can be calculated as:
K
log (E: w) = log po (En)
log Σ 2=1
J (5; (E)-1)
where == g= (E), E, = g (E)-1). and ( is the Jacobian matrix of the jth transform gw for j E (1,...,K), Several invertible flow architectures have been proposed in the literature, including radial and planar flows [33], autoregressive models [34-36] and models based on splines [37].
To derive a conditional density model qa (0" [M). conditional normalising flows map the original flow parameters via neural network model r 9 w [38, 39]. The resulting variational conditional density model then given by:
log gol (0
|0. &, 0) = log pr (0"; (0 & 0))
(20)
5.4 Batch parallel evaluations
Often simulations can be run in parallel by spawning multiple processes in a single machine or over a high-performance computing cluster. In this case, proposing batches with multiple simulation inputs can be more effective than running single simulations a sequence. Optimising the EIG w.r.t. a batch of inputs B == {&,, ) instead of single points, we obtain batch version of Algorithm 1. seeking batch that maximises the mutual information between the parameters
In this case, 0" and the resulting simulation outcomes, i.e.:
9 (o. B
(B), I(0": ( D,-1) EMCL,E,
log p(0") (D:1)
(21)
We optimise this objective with variational models that accept multiple conditioning observations q("" -.,,..,, (D). In practice, this simply amounts to replacing the single conditioning entries the models in Sec. 5.3 by the concatenated batch cormmutation-invariann deep set encoding [16, 40].

---

Thus
MAP
Predictior RMGT
Reverse
25
30-
043
15
04
042
38
02
65
100
150
200
300 120
200
300
250
MAP error
(b) RMSE
(c) DEL lp.)
(d) Posterior
Experimental results on synthetic data where the target posterior p' is unimodal. The first 3 plots show estimates for performance metrics as function of the number of simulations run (not including the initial data). Estimates were computed based on the posterior estimates for each method available during their run, with random using p(0") D-optimal and BACON using MCMC posteriors, and IMSPE using Dirac delta (reverse KL undefined, not shown) the MAP estimate as posterior estimates. Results are averaged over 10 trials, and shaded areas indicate #1 standard deviation. The rightmost plot shows the target posterior, with the true O" indicated by a star.
Figure I:
HAP
Redittior RMSE
feverse
05-
08%
-0.5
SACO's
0:0
-1.2-
Us
125
(a) MAP error
(b) RMSE
(c) DKL (Pc (p.)
(d) Posterior
Figure 2: Experimental results on synthetic data where the target posterior p" bimodal. See Fig. 1
for details, with the exception that the rightmost plot now shows the bimodal target posterior.
Experiments
In this section, we present experimental results on synthetic and real-data problems evaluating the proposed variational Bayesian adaptive calibration framework against baselines. Further experimental details can be found in Appendix A and in our code repository."
Performance metries. We evaluated each method against set of performance metrics, which we now describe. The caciimum-a--psttriiie (MAP) eтror measures the distance between the mode of the variational distribution and the true parameters 0" To measure the quality of the learnt model in predicting real outcomes, we also evaluated the root mean square error (RMSE) between the expected GP predictions under the learnt variational distribution and real outcomes:
RMSE == N ΣM (Eyco) [u(x;, [xxx 0" 0)] ) where y = =[x;)) + v,," are observations of the real process over a set of test points (x7 C X' placed on a uniform grid over the design space.
Information gain. Lastly, we also evaluated two sample-based estimates of the KL divergence [41]. Namely, D (Prllpo) corresponds to the KL divergence between the final MCMC posterior (given all simulations and real data) and the initial one (given only the real data and initial set of randomised simulations) both estimated over the learnt GP model. The column D (px (prllp") indicates the KL divergence between the final MCMC posterior Pr and the posterior p" with full knowledge of the simulator, which can be cheaply evaluated in this synthetic scenario. The average of (KI(PPT Po) is indicator for the expected information gain (1) of algorithm, given that it the expected relative entropy across the possible trajectories of observations. In contrast, DRL. (Rr(ppllp"" indicates how far the estimates are from the best possible posterior obtainable with a model that given the available real data and (a potentially infinite amount of) simulations.
'Code available at: https: : // thub. com/csiro- http//cioo-fumm/bacee funml/bacon
7

---

DKL
(pT p PO) T
(PT
pp]
DKL (PT PP) T
DKL
pp'
BACON
1.00 = 0.06
0.76 = 0.13
BACON
0.40 = 0.03
0.45 = 0.06
IMSPE
0.89 + 0.11
1.05 + 0.19
IMSPE
0.19 + 0.04
0.70 0.07
D-optim.
0.42 + 0.11
1.09 0.15
D-optim.
0.07 + 0.02
0.94 + 0.03
Random
0.62 + 0.07
1.18 0.13
Random
0.28 + 0.07
0.54 0.07
VBMC
0.53 + 0.02
VBMC
0.49 + 0.13
(a) Unimodal posterior
(b) Bimodal posterior
Table 1: Results for 2+2D synthetic problem after T 50 iterations (batch of B 4). Here DKL (prl Po) corresponds to the KL divergence between the final posterior (estimated after each algorithm's 's with all the data it collected) and the starting one (higher better), while (pr pppl"")) the KL between the final posterior and the posterior with full knowledge of the simulator p" (lower is better). All posteriors deviations were estimated from 10 independent runs.
sampled via MCMC using 4000 samples. Averages and standard
6.1 Baselines
Our algorithmic baselines were chosen to illustrate the main approaches currently available in the literature. All baselines are implemented as sequential methods. in the sense that their GP models are updated with the latest batch of observations before proceeding to the next iteration.
Random search. This baseline samples simulation designs & ((X) from uniform distribution over the design space X' and calibration parameters from the prior è p(0").
IMSPE with MAP estimates. The integrated mean squared prediction error (IMSPE) [42] criterion chooses designs & and calibration 0, parameters by minimising the GP prediction error:
IMSPE,
[[(f(z) (z;))(;0))) f(2), Di)
(2; Dr, f())) dz. (22)
The posterior MAP estimate 02 E argmaxp p(0 D) is used as a point estimate for the true 0" The integral is approximated sum uniform grid of designs and samples from the calibration prior, making IMSPE equivalent to active learning Cohn [43] and also a form of A-optimality [28].
D-optimal designs. We provide experimental results with additional baseline following D- optimality criterion, a classic experimental design objective. Optimal candidate designs according this criterion points of maximum uncertainty according to the model [28]. If model the simulator as the unknown variable of interest, this corresponds to selecting designs where we have maximum entropy of the Gaussian predictive distribution p( & 0 D,-1). This approach, therefore, simply attempts to collect an informative set of simulations according to the GP prior over the simulator h only, without considering the information in the real data. Running D-optimality on 0*, instead, would lead back the EIG criterion we use.
Variational Bayesian Monte Carlo (VBMC). Acerbi [44] presents an adaptive Bayesian quadra- ture method learn posterior distributions over models with black-box likelihood functions. The method estimates the posterior p(0 ly h) by modelling the log-joint log n: 0" 0'(h) as sample from Gaussian process. VBMC then learns variational posterior approximation by maximising lower-confidence bound over the ELBO given by the GP estimates. Calibration parameter queries 0, are obtained by optimising quadrature-base
acquisition functions. Regarding design points, simulations are always run on the set of real design points XxR in the observed data, which is fixed.
6.2 Synthetic experiments
For this experiment, we sampled function f ~ CP, k.) to use our simulator and compared different algorithms. Following a sparse GP approach [45], function sampled from GP can
'The original paper proposed analytic solutions to Eq. 22 tailored for specific kernels. However, we decided keep codebase generic work with different kernels, and therefore opted for numerical approximation.

---

DKL
PP)) t
DKL (PT
1
BACON
0.37 = 0.09
0.07 = 0.06
IMSPE
0.22 + 0.11
0.45 + 0.21
D-optimal
0.21 + 0.08
0.23 + 0.10
Random
0.32 0.09
0.20 + 0.14
VBMC
5.48 + 1.66
Table 2: Results on the location finding problem after 7= 30 iterations with B = 4, R = 20 "real"
data points and an initial set of 20 simulations. Estimates were averaged over 10 independent runs.
be approximated as f(z) x. k(z,, Zm) )K / uM, where uM Nuū.: Ex) is a sample from an M M-dimensional Gaussian, Zaz {z, M C X x θ x f0, 10, for a given M. As the number of points M + 30. if the pseudo-inputs ZM form a dense set, the approximate f should converge in distribution to sample from the Gaussian process GP(O, k). In our case, to sample Zas. we sample designs from a uniform distribution over the design space, calibration parameters from the prior, and fidelities from Bernoulli distribution with parameter set to 0.5. We also set üи := O and Σns- == ks (Z a. We repeatedly run loop of T iterations for each algorithm, with different random seeds.
We run each algorithm for T =- 50 iterations using batch of B =- designs per iteration. Each of the methods using GP approximations for the simulator are initialised with 20 observations and R = 5 real data points. To configure VBMC, we allow it to run an equivalent maximum amount of objective function evaluations. The design space is set as the 2- 2-dimensional unit box X' == [0, 1] and the "true" parameters are sampled from a standard normal prior p(0*) := NO00; 0, I) also over a 2D space, totalling 4-dimensional problem space.
Results are presented in Fig. 1 and 2. Fig. shows case where the GP Ga-sampled simulator led to unimodal target posterior. In this case, we see that BACON is able to achieve fast convergence in terms of MAP estimates and KL divergence towards the target posterior, while IMSPE dominates in terms of simulator approximation error as measured by the RMSE. As the posterior is unimodal and quite concentrated around the true parameter, it is natural that a method relying MAP estimates, such as RMSE, would perform well. In contrast, when the posterior is multimodal, as shown in the bimodal case in Fig. 2, MAP estimates are not necessarily reliable any more, as they might get stuck on a non-informattvve mode, leading to biased estimates for IMSPE and a significant drop in performance. Lastly, note that D-optimal and random designs can also lead to RMSE approaching the lowest (as determined by the noise level with S2 0.5) in some circumstances. However, these approaches do not directly provide posterior approximations and may fail in more complex scenarios.
In terms of final posterior estimates, Table 1 shows that VBMC estimates reach the closest to the full-knowledge target posterior p" in the unimodal case, while BACON is able surpass the other GP-emulation approaches in terms of information gain. For the bimodal case, however, we see that BACON gains an advantage over VBMC. Recall that VBMC relies on variational mixture of Gaussian distributions, while BACON applies conditional normalising flows for its posterior approximations, which lead to increased flexibility. In addition, despite the slightly worse performance than VBMC, BACON also provides a GP model that can be used as an emulator for the simulator (and to approximate the real process), while VBMC's focus is on approximating the likeliool..oo.
6.3 Finding the location of hidden sources
We consider the problem of finding the location of 2 hidden sources in a 2D environment following the setting in Foster et al. [16]. We are provided with R 20 initial measurements and an initial set of S = 20 randomised simulations without knowledge of the true parameters which the data was generated with. Sources are sampled from a standard normal, the design space limited to the unit box, and noise is sampled with o, 0.5. Our results are presented in Table 2, showing similar tendency in higher information gain for our method, and a very low KL w.r.t. p². Note that a higher information gain indicates more informative posterior, whose entropy will be much lower relative to the starting distribution, compared to the other methods. In addition, the ideal p², which a GP-based posterior should converge to in the limit of infinite data, is not known by the methods, only Po- Therefore, besides obtaining maximally informative data, we have shown that BACON is also efficient in approximating posteriors over black-box simulators, while also learning GP emulator.

---

SOFA
SarA
(a) Platform
(b) Real grasp
(c) Simulation
Figure 3: Soft-robotics grasping experiment. We calibrate soft materials simulator against real data from physical grasping from an automated experimentation platform
DKL (PT p
BACON
1.32 + 0.05
IMSPE
1.56 + 0.08
D-optimal
1.50 + 0.05
Random
1.48 + 0.07
Table 3: Soft-robotics simulator calibration final results after T 10 with B = 16 points per batch.
The target posterior p" was inferred using a large set of 1024 random simulations uniformly covering the design and parameter space. Performance averaged over 4 independent runs.
6.4 Soft-robotic grasping simulator calibration
For this experiment, we are provided with dataset containing R = 10 real measurements of the peak grasping force of soft robotic gripper designs on range of testing objects (see Fig. 3). The gripper designs follow fin-ray pattern parameterised by geometric parameters [46], and are interested in estimating unknown physics parameters, the Young's modulus of clasticity and the coefficient of static friction with the objects. To simulate the gripper designs, use the SOFA framework [47] to reproduce the grasping scenario and provide an estimate of the peak grasping force. In particular, for this paper, focus on the grasping of a spherical object, which provides simpler geometry and lower discrepancy with respect to real data measurements compared to more complex objects. This experiment provides with benchmark where simulations are expensive run, taking from minutes to a few hours to run (depending on mesh resolution) on a high-performance computing platform. Therefore, it is important choose minimum amount of informative simulations.
Our results are shown in Table 3. Each algorithm was initialised with a set of 123 random simulations and run for T = 10 iterations. The results show that BACON achieves the closest approximation to the target posterior. IMSPE highly concentrated its parameter choices around its posterior mode estimate, while other baselines were too spread, both leading inferior posterior approximations (see Fig. 4 in the appendix) and showing the advantage of BACON's joint optimisation and inference.
7 Conclusion, limitations and future work
We have developed BACON, a Bayesian approach that carries out parameter calibration of computer models and optimal design of experiments jointly. does by optimizing information-theorett criterion so that input designs and calibration parameters are selected to be maximally informa- tive about the optimal parameters. Our method provides full posterior over optimal calibration parameters as well as an accurate Gaussian process based estimation of the computer model (i.c., an emulator). One of the main limitations of the presented framework, however, scalability to large datasets, due to the cubic computational complexity of exact inference with GPs. A potential extension with scalable sparse variational GP models [48] using a conditional distribution model for the inducing points is discussed in Sec. B.2. We emphasize that our proposed method is still applicable to many real practical settings, where the problem constraints do not demand a very large number of simulation samples. Lastly, we also note that the method can be adapted to work with rector-valued observations by the use of multi-output GP models [49]. Further discussions on limitations and future work can be found in appendix (see Appendix B and C).
10

---

Acknowledgements
This project was supported by resources and expertise provided by CSIRO IMT Scientific Computing. We are also grateful for the support of CSIRO's Data61 soft-robotics team, especially Josh Pinskier, Xing Wang, Lois Liow, Sarah Baldwin, James Brett and Vinoth Viswanathan, in the experimental data collection and simulations setup for the soft-robotics calibration problem.
References
[1] Mare C. Kennedy and Anthony O' O'Hagan. Bayesian calibration of computer models. Journal of the Royal Statistical Sociery: Series B (Statistical Methodology), 63(3):425 464, 2001.
[2] Carl E. Rasmussen and Christopher K. Williams. Gaussian Processes for Machine Learning. The MIT Press, Cambridge, MA, 2006.
[3] Michael U. Gutmann and Jukka Corander. Bayesian optimization for likelihood-free inference of simulator-based statistical models. Journal of Machine Learning Research, 17, 2016.
[4] Erin R. Leatherman, Angela M. Dean, and Thomas J. Santner. Designing combined physical and computer experiments to maximize prediction accuracy. Computational Statisties and Data Analysis, 1 113:346-366 3622663662222 2017.
[5] Sébastien Marmin and Maurizio Filippone. Deep Gaussian processes for calibration of computer models (with discussion). Bayesian Analysis, 17(4):1301-1350, 2022.
[6] Tom Rainforth, Adam Foster, Desi R. Ivanova, and Freddie Bickford Smith. Modern Bayesian Experimental Design. Statistical Science, 39(1): 114, 2024.
[7] Thomas M. Cover and Joy A. Thomas. Elements of Information Theory. John Wiley Sons, 2005.
[8] Ronald A. Fisher. The design of experiments. Oliver & Boyd, Oxford, England, 1935.
[9] Stewart Greenhill, Santu Rana, Sunil Gupta, Pratibha Vellanki, and Svetha Venkatesh. Bayesian optimization for adaptive experimental design: review. IEEE Access, 8: 13937-113948, 2020.
[10] J. Kiefer. Optimum experimental designs. Journal of the Royal Statistical Society. Series B (Methodological), 21(2):272-319 1959.
[11] Kathryn Chaloner and Isabella Verdinelli. Bayesian experimental design: A review. Statistical Science, 10(3):273-304, 1995.
[12] Elizabeth G. Ryan, Christopher C. Drovandi, James M. Megree, and Anthony N. Pettitt. A review of modern computational algorithms for Bayesian optimal design. International Statistical Review, 84(1):128-154, 2016.
[13] Adam Foster, Martin Jankowiak, Eli Bingham, Paul Horsfall, Yee Whye Teh, Tom Rainforth. and Noah Goodman. Variational Bayesian optimal experimental design. In 33rd Conference on Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada, 2019.
[14] Steven Kleinegesse and Michael Gutmann. Efficient Bayesian experimental design for implicit models. In Kamalika Chaudhuri and Masashi Sugiyama, editors, Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics (AISTATS), Proceedings of Machine Learning Research, Naha, Okinawa, Japan, 2019. PMLR.
[15] Joakim Beck, Ben Mansour Dia, Luis FR Espath, Quan Long, and Raul Tempone. Fast Bayesian experimental design: Laplace- importance sampling for the expected information gain. Computer Methods in Applied Mechanics and Engineering. 334:523-553, 2018.
[16] Adam Foster, Desi R. Ivanova, Ilyas Malik, and Tom Rainforth. Deep Adaptive Design: Amortizing sequential Bayesian experimental design. In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning (ICMI. 2021), I), volume 139 of Proceedings of Machine Learning Research, pages 3384-3395. PMLR, 2021.
11

---

[17] Tom Blau, Edwin V. Bonilla, Iadine Chades, and Amir Dezfouli. Optimizing sequential experimental design with deep reinforcement learning. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, Baltimore, Maryland, USA, 2022. PMLR.
[18] Tom Rainforth, Adam Foster, Desi R Ivanova, and Freddie Bickford Smith. Modern bayesian experimental design. arXiv preprint arXiv:2302. /4545, 2023.
[19] D. Busby and M. Feraille. Adaptive design of experiments for calibration of complex simulators An application to uncertainty quantification of mature oil field. Journal of Physics: Conference Series, 135, 2008.
[20] Christophe Andrieu, Nando De Freitas, Arnaud Doucet, and Michael I. Jordan. An introduction to MCMC for machine learning. Machine Learning, 50(1-2):5-44 -2:,544334450 2003.
[21] Soumalya Sarkar, Sudeepta Mondal, Michael Joly, Matthew E. Lynch, Shaunak D. Bopardikar, Ranadip Acharya, and Paris Perdikaris. Multifidelity and multiscale Bayesian framework for high-dimensional engineering design and calibration. Journal of Mechanical Design, 141(12), 2019.
[22] Rafael Oliveira, Lionel Ott, and Fabio Ramos. No-regret approximate inference via Bayesian optimisation. In 37th Conference on Uncertainty in Artificial Intelligence (UAI). PMLR, 2021.
[23] Luigi Acerbi. Variational Bayesian Monte Carlo with noisy likelihoods. In Larochelle, M Ranzato, R Hadsell, M F Balcan, and H Lin, editors, 34th Conference on Neural Information Processing Systems (NeurIPS 2020), volume 33, pages 821 1-2222, 2020.
[24] Marko Järvenpää, Michael U. Gutmann, Aki Vehtari, and Pekka Marttinen. Parallel Gaussian process surrogate Bayesian inference with noisy likelihood evaluations. Bayesian Analysis, 2020.
[25] Mucahit Cevik, Mehmet Ali Ergun, Natasha K Stout, Amy Trentham- Dietz, Mark Craven, and Oguzhan Alagoz. Using Active Learning for Speeding up Calibration in Simulation Models. Medical decision making: an international journal of the Society for Medical Decision Making, 36(5):581-593, 2016.
[26] Yuan Tian, Manuel Arias Chao, Chetan Kulkarni, Kai Goebel, and Olga Fink. Real-time model calibration with deep reinforcement learning. Mechanical Systems and Signal Processing, 165 (July 2021):108288 2022.
[27] Andreas C. Damianou and Neil D. Lawrence. Deep Gaussian processes. Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS), 31, 2013.
[28] Andreas Krause, Ajit Singh, and Carlos Guestrin. Near-optimal sensor placements in Gaussian processes: Theory, efficient algorithms and empirical studies. Journal of Machine Learning Research, 9:235-284 2008
[29] Daniel Golovin and Andreas Krause. Adaptive submodularity: Theory and applications in active learning and stochastic optimization. Journal of Artificial Intelligence Research, 42: 427-486, 2011.
[30] David Barber and Felix Agakov. The im algorithm: A variational approach to information maximization. In Proceedings of the 16th International Conference on Neural Information Processing Systems, NIPS 03, page 201-208 Cambridge, MA, USA, 2003. MIT Press.
[31] Patrick Dallaire, Camille Besse, and Brahim Chaib-Draa. An approximate inference with Gaussian process to latent functions from uncertain data. Neurocomputing, 74:1945-1955. 2011.
[32] Andreas C. Damianou, Michalis K Titsias, and Neil D Lawrence. Variational inference for latent variables and uncertain inputs in Gaussian processes. Journal of Machine Learning Research, 17(1):1-62. 12, 2016.
12

---

[33] Danilo Jimenez Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International Conference on Machine Learning (ICML 2015), volume 2, pages 1530-1538, Lille, France, 2015.
[34] Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, llya Sutskever, and Max Welling. Improved variational inference with inverse autoregressive flow. In D. Lee, M. Sugiyama, U. Luxburg, Guyon, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 29, 2016.
[35] George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation. In Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30, 2017.
[36] Chin Wei Huang, David Krueger, Alexandre Lacoste, and Aaron Courville. Neural autoregres- sive flows. In 35th International Conference on Machine Learning (ICML 2018), volume 5, pages 3309-3324, Stockholm, Sweden, 2018.
[37] Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows, In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32, 2019.
[38] Christina Winkler, Daniel Worrall, Emiel Hoogeboom, and Max Welling. Learning likelihoods with conditional normalizing flows. arXiv preprint arXiv:1912.00042 2019.
[39] AbdeIrahman Abdelhamed, Marcus A Brubaker, and Michael S Brown. Noise flow: Noise modeling with conditional normalizing flows. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 3163-3173, 2019.
[40] Manzil Zaheer, Satwik Kottur, Siamak Ravanbhakhsh, Barnabás Póczos, Ruslan Salakhutdinov, and Alexander Smola. Deep sets. In Advances in Neural Information Processing Systems, volume 30, pages 3392-3402 2017.
[41] Zoltán Szabó. Information theoretical estimators toolbox. Journal of Machine Learning Research, 15:283-287, 2014.
[42] Scott Koermer, Justin Loda, Aaron Noble, and Robert B. Gramacy. Active Learning for Simulator Calibration. arXiv, 2023. URL http: //arxiv...vv//// org/ aa/aabaa/23310 10228.
[43] Annie Sauer, Robert B. Gramacy, and David Higdon. Active Learning for Deep Gaussian Process Surrogates. Technometrics, 65(1):1-39, 2022. ISSN 15372723, doi: 10. 10.1080/0440100 2021.2008505. URL httpp::///oot org/ /tt. http://.....r.1....11010//////// 1080/00401706. 2000000055.....
[44] Luigi Acerbi. Variational Bayesian Monte Carlo. In 32nd Conference on Neural Information Processing Systems (NeurIPS 2018), Montréal, Canada, 2018
[45] Edward Snelson and Zoubin Ghahramani. Sparse gaussian processes using pseudo-inpus In Y. Weiss, B. Schölkopf, and J. Platt, editors, Advances in Neural Information Processing Systems, volume 18, 2005.
[46] Xing Wang. Bing Wang, Joshua Pinskier, Yue Xie, James Brett, Richard Scalzo, and David Howard. Fin-bayes: A multi-objective bayesian optimization framework for soft robotic fingers. Soft Robotics, 2024.
[47] François Faure, Christian Duriez, Hervé Delingette, Jérémie Allard, Benjamin Gilles, Stéphanie Marchesseau, Hugo Talbot, Hadrien Courtecuisse, Guillaume Bousquet, Igor Peterlik, and Stéphane Cotin. SOFA: A Multi-Model Framework for Interactive Physical Simulation. In Yohan Payan, editor, Soft Tissue Biomechanical Modeling for Computer Assisted Surgery, volume 11 of Studies in Mechanobiology, Tissue Engineering and Biomaterials, pages 283-321. Springer, June 2012. URL https: inria a. hal cal.ssieen sci wecceenneeallll81155551
[48] Michalis K. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In International Conference on Artificial Intelligence and Statistics (AISTATS), Clearwater Beach, Florida, USA, 2009.
13

---

[49] Mauricio A Alvarez, Lorenzo Rosasco, Neil D Lawrence, et al. Kernels for vector- valued functions: review. Foundations and Trends@ in Machine Learning, 4(3):195-266, 2012.
[50] Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D. Goodman. Pyro: Deep Universal Probabilistic Programming. Journal of Machine Learning Research, 2018.
[51] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy. BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization. In Advances in Neural Information Processing Systems 33, 2020. URL http://aree orgg /9100 06403.
[52] Matthew D. Hoffman and Andrew Gelman. The no-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15:1593 15:1593-1623, 2014.
[53] Michalis Titsias and Neil Lawrence. Bayesian Gaussian process latent variable model. In International Conference on Artificial Intelligence and Statistics (AISTATS), pages 844-851, 2010.
[54] Vidhi Lalchand, Aditya Ravuri, and Neil D. Lawrence. Generalised GPLVM with stochastic variational inference. In Gustau Camps- Valls, Francisco J. R. Ruiz, and Isabel Valera, editors, Proceedings of The 25th International Conference on Artificial Intelligence and Statisties, volume 151 of Proceedings of Machine Learning Research, pages 7841-7864. PMLR, 2022.
[55] James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, UAI' page 282-290, Arlington, Virginia, USA, 2013. AUAI Press.
14

---

Algorithm 2 BACON (split training)
Do = {Kn, RR);
for 1 E (1, (1,..., ; do
Mi-1, < UpdateGP(D,-i)
10 MCMC (D-1D= p(
p(0" (D-1)) Pd-1 == SA -i=1% s0;
q: - TrainFlow (ft-1. D1-1)
{ 1- t OptimiseDesignss(A qt. 1, D(-1) 0r.: == h(R de,i) (parallel) for E (1,.... B}
(Run batch of simulations] {Update GP dataset)
D. = D (Re,i, Or. 01,
end for
A Additional details on the experiments
For all experiments, we use conditional normalising flows as the variational model for BACON. Our implementation for BACON and most of the baselines, except for VBMC, is based on Pyro probabilistic programming models [50]. Gaussian process modelling code is based on BoTorch" [51]. The flow architecture is chosen for each synthetic-data problem by running hyper-parameter tuning with a simplified version of the problem. Most Gaussian process models are parameterised with Matérn kernels [2, Ch. 4] and constant or zero mean functions. Pyro's MCMC with its default no-U-turn (NUTS) sampler [52] was applied to obtain samples from ((0'|Dr-1 at each iteration t. KL divergences are computed from samples using nearest neighbours estimator implemented in the information theoretical estimators (ITE) package" [41].
A.1 Synthetic GP problem
The GP prior was set with given by a squared exponential kernel and kg given by a Matérn kernel with smoothness parameter set to 2.5 [2]. The conditional normalising flow was configured with 2 layers of neural spline flows [37]. Batches of arbitrary size are used for conditioning via permutation invariant set encoder, similar to Blau al. [17]. with 2-layer, 32-units-wide fully- connected hyperbolic tangent neural network passing through a summation at the end. Gradient- based optimisation is run using Adam with learning rate 10-3 for the flow parameters and 0.05 for the simulation design points, both using cosine annealing with warm restarts as a learning rate scheduler. 256 samples subsampled from the MCMC posterior estimate expectations for both this and the Sraation-inding problem.
Algorithm with split training. For the synthetic GP problem, provide more detailed pseudo- code of our algorithmic implementation using an option for training the conditional normalising flow and optimising the designs separately. Specifically, had flexible optimisation loop, where we had the option to separate the training of the conditional normalising flow model from the optimisation of the design points, as shown Algorithm 2. This approach can make the algorithm more stable, though at the cost of a longer runtime. This option was only applied to the Ga-based synthetic experiments, while for the other experiments we ran the full joint optimisation over both the simulation inputs (x, 0) and the variational parameters of the conditional model g.
applied MCMC to estimate our posteriors and
A.2 Location finding problem
For this experiment we used more up-to-date Zuko" implementations of the conditional normalising flow models. which were again set as neural spline flows [37] combined with set encoder to condition on arbitrary batch sizes. Further architectural details can be found in our code repository. 256 samples
For VBMC, we used its author's Python implementation at: https: httpb//ttttttttt///g www./cerbilabb/pyeeeaeacooeooc
Boooorch: httpe://botorch.. org
s IITE package: http::://bttiieee...r////a/////e
Zuko: attps://zule how.dttheeoo...///taabbb.e.
15

---

Algorithm 3 TrainFlow
input Pt. Dr
for {1, ,..., N do
(Ó. |=| (1
c)pe + cp
[ U(())
[0: [=1 it
(ds. (ij=1 (N 1 (0:) Ex( P (0))))))
S.B
Φ - Φ Σニ V log Us
0: Dr U [8, 0ル
end for output 94
Algorithm 4 OptimiseDesigns
input Qt. 0.-1 Di-1
θ = - ~ (1 c) + ep
X [8 ((()
for for, do
(0) ~Pt-
19 ~N(H,,11 (X, (; E,- (x
(X,0) t -(X,6) + Ei Xx log (0) D. D.
end for output {2, m
were subsampled from the MCMC posterior at each iteration to estimate expectations for EIG lower bound computations. The simulations kernel & was Matérn 2.5 kernel. For this experiment we did not model the error term, leaving it with a zero kernel, since data is generated directly from the simulator with no further error component, only Gaussian noise with standard deviation of 0.5. Final KL estimates were computed using the maximum-a-poster model learnt with the random search approach to minimise biases in the estimate of DEL (pz|lpo) due to differing GP hyper-parament parameters across baselines.
dyper-parameter of the GP
A.3 Soft-robotics simulation problem
The prior for the calibration parameters p(0") in this experiment consisted a 2-dimensional standard normal transformed through sigmoid and affine transform composition to provide a smooth uniform distribution over a pre-specified range for the calibration parameters. Such smooth approximation allows gradients to be computed near the edges of the parameter space while not allowing optimisation to take the calibration parameter candidates outside the uniform prior boundaries, since these would be placed at infinity under the normalised space. The conditional normalising flow model used Zuko's 's implementation of neural spline flows with 10 transform layers. The set encoder consisted of a 2-layer fully connected 32-unit-wide neural network encoding each input into an 8-dimensional output which was then summed and passed through as the context input condition the flow. Adam again was used for optimisation with learning rate of 0.001 for the flow and 0.05 for the simulation inputs. Monte Carlo expectation estimates used 256 samples from the current MCMC posterior each joint optimisation step.
A.4 Hyper-parameter tuning
Besides the GP hyperparameters (e.g., lengthscales, noise variance, etc.), which had be tuned for the non-GP-basee problems, there are optimisation settings (i.e., step sizes, scheduling rates, etc.), conditional density model hyper-parameters (i.e., normalising flow architecture), and other algorithmic settings, e.g., the designs batch size B. The latter is dependent on the available computing resources (e.g., number of CPU cores or compute nodes for simulations in a high-peeformance computing system). We tuned optimisation settings and architectural parameters for the conditional
16

---

Reference posterior
0.25
-1
0.10
0.05
+3
0.00
(a) Reference posterior p*
BACON
IMSPE
0.25
025
0.15
0.15
0.05
0.05
(b) BACON final posterior
IMSPE baseline final posterior
D-optimal
Random
2
0.20
0.20
1
0
-1
0.10
-1
0.10
-2
0.00
-3
0.00
(d) D-optimal baseline final posterior
(c) Random baseline final posterior
Figure 4: Final posterior approximations p(0")D7) and simulation parameter (red crosses) choices by each method for the soft-robotics simulator calibration problem after one of the runs. The target/reference posterior (a) was inferred using large number (1024) of simulations following a Latin hypercube pattern over the combined design X and calibration parameters space and uniform prior p(0) over the same range as the smooth uniform prior the algorithms used. The posteriors are plotted as a 2D histogram over the normalised range (after an affine and sigmoid transform). which the algorithms used for optimisation. The KL divergences in Table are computed with respect to this reference posterior. Also note that the simulation parameters Ô in the plot correspond to different algorithmic choices for design inputs &, which are 9-dimensional variables that are not plotted here.
17

---

normalising flows via Bayesian optimisation with short runs (e.g., 10-20 iterations) on the synthetic problem. However, depending on the number of parameters, simpler approach. like grid search, might be enough. GP hyper-parameters were optimised online via maximum a posteriori estima- tion after each iteration's batch update. Further implementation details can be found in our code repository." 10
B Extensions of the proposed approach
In the following, present two extensions to deal with limitations of the current approach. Namely, we can amortise inference over the calibration posterior by reutilising the learnt conditional distribu- tion models priors, instead of having run, for example, MCMC. Secondly, present derivations for a scalable sparse GP version of our method.
B.1 Amortisation
We use conditional variational distribution model for ((0"))). The main advantage of training conditional model is that, once new data is observed, we readily obtain approximation to the new posterior p(0" (D;) p(0" lr, St, 0t, Dr-1) (0cc) There is, therefore, potential to reuse the variational posterior as the prior for the next iteration, and all the optimisation concentrated within a single loop.
Approximate objective. We are still left with terms dependent on the posterior from the previous iteration p00 (D:-1) in Eq. 15. Firstly, however, note that the denominator inside the expectation is constant w.r.t. the optimisation variables, not affecting the maximiser. Secondly, we may replace the joint predictive distribution p(Ù. 0" & D. by approximation using the previous optimal variational posterior Qt-1 as:
p(ù. 0" &R, 0 D. - (1) 0" & 0) p(0 0" x D,,,D.1100 0
(23)
where (0*) == gr- (0"10r1)))) R p DD-I. The following objective then approximately shares the same set of maximisers as the variational lower bound EIG,(8,)) 0, q):
& 0, q: E
argmax E 10.0000)))886) [log ((0") [0]
(24)
In practice, reusing the variational conditional posterior may tend to degenerate the approximation over time. However, that can be corrected by rerumning MCMC variational inference scheme over the data to obtain a fresh new posterior at every few iterations.
B.2 Conditional sparse models for large datasets
Computing the variational EIG requires evaluating expectations with respect to the posterior predictive distribution p( 0" x, Ô Dy). Note, however, that, as 0" appears inside a matrix inversion in the GP predictive (Eq. 8), each sample of p(ū x, 0, D:) requires a (NN3) computation cost, where N₂ R+ is the number of data points iteration N. This cost may quickly become prohibitive for reasonably large datasets, which are easily obtainable in batch settings (Sec. 5.4), rendering EIG computations infeasible. To scale our method handle large amounts of data, then need GP models that can reduce this computational complexity, while still allowing us to obtain reasonable EIG estimates.
B.2.1 Variational sparse GP approximation
We consider an augmentation the original GP model which allows to sparsify its covariance matrix, reducing the computational complexity of GP predictions. Following the variational sparse GP approach [48], let == fZ..) e R-M denote a vector of M inducing variables representing unknown function values given set of pseudo-inputs Z., The joint distribution between observations y.
Codd available at: https: //gi thub. wwmccccco////eeeeoooooooeeccc.ce fow.rrooofummmbbonnn
18

---

function values f: := f(Z(6")). inducing variables u and the unknown parameters 0 can be written as:
U. u1 = ply, (lf",
(25)
where ply = N,,;;,, Ey),
p(f|u, 0*) = ju (o') ","",, and puu) = (u; 0, K.n), using notation shortcuts Kun := iZz,, Zn), K ju (0)) k(Z(0"), Z..)
(26)
and Kj ,j₃00) :- k(Z) (0))). We may now formulate an evidence lower bound (ELBO) based on the joint variational density q(f, u, 0²) as:
log ply, F.f,uu. u, 0*)
logply) E Fg((...0))
+ DKL u, (")ll((f u. 0yl))
g(f. u. 0¹)
(27)
log fly,f,u, 0')
2 E (f...r))
g(f. u. 0*)
Since (qL(,((, u, 8¹) )))llp(( u, 0) 2 O. and if and only if q(f,, u, 0*) = p(f. u, 0 ly), max- imising the ELBO above w.n.. provides us with an approximation to the joint posterior. Choosing q(f, 0*) == p(f u, 0 )g( 0*) simplifies the ELBO to [53]:
ply viyl))p(u)p(00) (uu
logp(y) 2 E vft.u..-) log
(28)
Sparse variational GP approaches can reduce the computational complexity of Bayesian inference on GPs to ((NM) or even ((M3) [48, 54), where N is the number of data points.
B.2.2 Structure of the joint variational posterior
If would take mean-field approach setting quu, 0¹) 0() further simplify, leading to a few computational advantages, as explored by Bayesian GP-LVM methods [53, 32, 54]. However, in our experimental design context, this approach leads to a few issues. Firstly, using the mean-field posterior as a replacement for our joint posterior breaks the dependence between and 0*, leading their mutual information (a.k.a. EIG) to be zero regardless of the design inputs & and Ó. Secondly, although u and 0 are independent according to their priors (Eq. 25), they become dependent when conditioned on the data. In fact, the true posterior over given the data and the true parameters 0" is exactly Gaussian:
((") the ELBO above would
p(u[D,, 0-) = Nu;;; H4 (Z,; 0²), (Z,,, ZZ.
(29)
where (e(: 0*) and by (,,;;;; :; 0*) are given by Eq. 9 and Eq. 10, respectively. Note, however, that the posterior over 0" should not be Gaussian for general non-linear kernel k. Therefore, it makes more sense for us to model q(u, 0*) == q ((]("") ). Moreover, learning a Gaussian conditional model over and flexible variational distribution over 0* should be enough to allow us recover the true posterior, since puu, 0" D.) = p(u DD 0²) )p(6") (D:).
Optimal variational inducing- cointt distribution. Given 0* E O, we have a standard sparse GP model. The optimal variational inducing- point distribution available in closed form following standard results [48] as:
g" (ul00) = NNu: (4, E.,))),
(30)
where the distribution parameters are:
(0) == Knn (Knu + (0))
1 (
(31)
(0) := K.n. (Km + (1)( Knn.
(32)
and the conditional V matrices are given by:
4 (0) K (0)2
(33)
22( (0)=
(34)
for 0 E O. The computational cost of sampling predictions with this model then reduces from (NN))
to ONN
19

---

Parametric variational inducing distribution. To further reduce the computational cost of pre- dictions, we may accept a sub-optimal conditional variational inducing-pons distribution given by parametric model:
gc (<(100) == Nu; mc(0"), ((())))(((
(35)
following the architecture in Sec. 5.3. This formulation allows approximate the evidence lower bound Eq. 28 w.r.t. quu (0)) via mini hini-batching [see 55]. To do so, we approximate fi := f(z,) via conditionally independent samples given for E in Eq. 28 decomposes as a sum which is amenable mini-batching:
NN.. As result, the data-dependent term
N
(36)
where (c(f. ul0") p(fil 0" (a<(u(("), The variational parameters Ç need be optimised within a second optimisation loop after the data update in Algorithm I w.r.t.:
(<C) :- (0-)
Σ
E
("e (37)
Although the GP update is no longer available in closed form, we gain computational efficiency for large volumes of data. Applying mini-batches of size I. < N to Eq. 37 results in a computational cost (LM)) (or ((M)), if M > L),, which is smaller than the cost O(NM2) of the optimal variational distribution q" (u(0")
C Further discussion on limitations
High-dimensional settings. The dimensionality of our search space consists of the combined dimensionality of the designs X' and calibration parameters space θ which can be large in practical applications. In general, in higher dimensions, one is to expect that the algorithm will require larger number of iterations to find suitable posterior approximations due to the possible increase in complexity of the posterior. The analysis of such complexity, however, is problemn-deenndet and outside the scope of this work. In addition, note that we do not mean that the per-iteration runtime is directly affected, since what dominates the cost of inference sampling from the GP, whose runtime complexity is dominated by the cube of the number of data points due to matrix inversion operation, while being only linear in dimensionality.
Gaussian assumptions. We make Gaussian assumptions when modelling the simulator and the approximation errors, which can be seen as restrictive for some applications. However, the errors are Sub-Gaussian (i.e., its tail probabilities decay faster than that of a Gaussian), as is the case for bounded errors, we conjecture that GP model can still be a suitable surrogate, as it would not underestimate the етrоr uncertainty. If the error function is sampled from some form of heavy-tailed stochastic process (e.g., a Student-T process), the GP would, however, tend under estimate uncertainty and lead to possibly optimistic EIG estimates that make the algorithm under-explore the search space. Changing from GP model to another type of stochastic process model that can capture heavier tails would be possible, though require significant changes to the algorithm's 's predictive equations. We, however, believe that most real-world cases would present errors which are least bounded (and therefore sub- uub-Gaussian) with respect to the simulations.
20

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
Experimental results confirm the main claims in the introduction and abstract.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Limitations are discussed in the conclusion section and Appendix C.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question:
For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
21

---

Justification:
Theoretical results have not been derived for this paper.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: Code is included, though soft-roboties experiment data is protected.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architccture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reprodueibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
22

---

Answer: [Yes]
Justification:
Code (included with submission) will be made public for most of the results, except soft-robotics data, which is subject to internal restrictions.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (uttps://nips.c cc/
cooiisssssiiicyyyccccoooc oolicy) for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: wttpc//guiieee//oooeeeeeeeee
/ /nipss.. cc/
ississiinns asionPolicyy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Main details are provided, and the code has been made publicly available.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Standard deviations are reported with every performance plot and table.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
23

---

It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [No]
Justification: Full experiment details will be provided for camera-ready version.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https http///eerrpp.eeeee cc wwtppcccc/ccaeeeeeeeee wtw.tiiiiieeeeeeeeeeeeeeeeeeeeee
Answer: [Yes]
Justification:
No data subject to the NeurIPS Code of Ethics has been used in this work.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: This work is of a theoretical nature introduce new methods for general class applications, potentially science and engineering.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.

---

The conference expects that many papers will be foundational research and not tied particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [No]
Justification:
There are currently no plans to release any dataset other than synthetic data
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Lieenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: Mostly open-source code has been used to base this project
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
25

---

If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [No]
Justification: Except for code to reproduce experiments,
assets introduced.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mooel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This paper does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
26


Reviews:
[{"review_id": "rjcrwqc3ws", "reviewer": "Reviewer_bRKR", "summary": "This paper considers the problem of calibration of computer models as an active learning problem. Given the objective being maximizing the expected information gain (EIG) about calibration parameters, and based on the assumption of linear dependency between simulator outcome and true observation, this work proposes a Gaussian process model that jointly models true observations, calibration parameters and simulator outcomes, and maximization of EIG is performed based on this model. Due to the intractability of the EIG objective, the authors further propose to use variational objective in replace of the original EIG in finding the next design parameters and calibration parameters to sample.", "questions": "- The proposed approach certainly makes sense when both spaces of calibration parameters and design parameters are low-dimensional. It may be helpful to include a discussion about choosing appropriate method when the number of calibration parameters is bigger or the number of design parameters is bigger, or both.", "limitations": "Please also include limitation of the dimensionality issue as pointed out above.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- The problem setting is interesting and is indeed important in engineering and physical sciences as computer simulators are often used in those areas.\n- The approach of using a single GP to model calibration parameters and design parameters jointly is novel.", "weaknesses": "- Optimizing calibration parameters and design parameters jointly seems to make the problem harder because dimension of the search space is the sum of both spaces, therefore, the applicability of this method to real-world problems may be limited.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their comments and feedback. Indeed, the dimensionality of the problem consists of the sum of the dimensionalities of the design space and the calibration parameters space. The purpose of this paper, however, was to propose a general method for Bayesian calibration via adaptive experimental design and some of the aforementioned limitations were out of our current scope. In high-dimensional settings, there are two main issues to address.\n\n* GP models with conventional translation-invariant radial kernels start to lose predictive performance in high-dimensional spaces due to the disappearing volume of the $n$-ball at higher dimensions.\n* Finding optimal design points in higher dimensions may become increasingly difficult due to the vast spaces.\n\nThe first issue may be addressed by incorporating more specialised, possibly non-stationary, kernels, which better capture correlations in higher-dimensional spaces. For example, Li et al. (2024, below) have recently shown that infinite-width Bayesian neural networks (BNNs) maintain satisfying modelling capabilities in higher-dimensional active learning settings. Due to their equivalence with GPs (Lee et al., 2018, below), infinite-width BNNs can be implemented as a GP model with a kernel given by the infinite-width limit covariance function of a BNN of a given architecture. A few results are known in the literature and software packages are available for such. The main point, however, is that this approach would not lead to any modifications to our proposed algorithm, since it is still GP-based. We will make these clarifications along with the corresponding set of reference in the revised manuscript.\n\nThe second issue may result in slower, though still computationally feasible, convergence of optimisation approaches based on stochastic gradient descent (see Arous et al, 2021, below). To mitigate possible issues with finding good local optima, multiple random restarts can be applied.\n\nIn general, in higher dimensions, one is to expect that the algorithm will require a larger number of iterations $T$ to find suitable posterior approximations due to the possible increase in complexity of the posterior $p(\\boldsymbol{\\theta}^*|\\mathcal{D}_t)$. The analysis of such complexity, however, is problem dependent and outside the scope of this work. In addition, note that we do not mean that the per-iteration runtime is directly affected, since what dominates the cost of inference is sampling from the GP, whose runtime complexity $\\mathcal{O}(t^3 + dt^2)$ is dominated by the cube of the number of data points due to a matrix inversion operation, while being only linear in dimensionality $d$ (the term $dt^2$ is the cost of computing the $t^2$ entries in the kernel matrix $\\mathbf{K}_t$). We will expand on these limitations in the final version of the paper.\n\n**References**\n* Arous, G. Ben, Gheissari, R., and Jagannath, A. (2021). Online stochastic gradient descent on non-convex losses from high-dimensional inference. *Journal of Machine Learning Research*, 22.\n * Lee, J., Bahri, Y., Novak, R., Schoenholz, S. S., Pennington, J., and Sohl-dickstein, J. (2018). *Deep neural networks as Gaussian processes*. In International Conference on Learning Representations (ICLR).\n * Li, Y. L., Rudner, T. G. J., and Wilson, A. G. (2024). *A Study of Bayesian Neural Network Surrogates for Bayesian Optimization*. 2024 International Conference on Learning Representations (ICLR 2024).", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "i8vTtL3nKr", "reviewer": "Reviewer_Z11Z", "summary": "This paper addresses the challenge of calibrating expensive-to-evaluate computer models using Bayesian adaptive experimental design. The novelty of the proposed method (BACON) lies in using the expected information gain (EIG), which is a principled information theoretic criterion for active learning, to perform calibration of models. Another point of departure from the existing literature is the fact that BACON performs active learning in the joint space of design and parameters.", "questions": "* Can you say something about the tightness of the bound in Section 5.1? What does it depend on? (I suppose these would have been discussed in ref [13] but it would be nice if mentioned here as well)\n* Can you explain a bit more why VBMC performs better in the synthetic experiments but fares poorly in the other experiments? Does it have anything to do with the dimensionality of the problem?\n* Reporting the KL divergence between the prior and the posterior after T iterations tells us how much the posterior has changed, but that does not mean we are converging to the true posterior, right (we may be confidently biased)? Is the KL divergence between $p_T$ and $p^*$ a better indicator for accuracy?", "limitations": "The authors have discussed the limitations adequately.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This is a technically solid and well-written paper which I really enjoyed reading. The idea of using the EIG criterion over the joint space of $\\theta$ and $x$ may sound simple and straightforward once you see it written down that it is almost surprising that no has done this before.  Many good papers seem obvious once you read them, and I feel like this paper belongs in that category. Despite being a notation heavy paper, the authors did a good job of making the writing clear and concise.", "weaknesses": "I do not have any major concerns. Some questions/comments that might help improve the paper are as follows:\n* It would have been nice to see an ablation study where the design and parameters are not jointly optimized over, in order to ascertain the benefit of doing so.\n* The error function $\\varepsilon$ is also modellled as a GP. I wonder how does BACON's performance get affected if this assumption does not hold (i.e. when this error model is misspecified).\n* Perhaps a discussion of the hyperparameters/settings of the proposed method would be nice to have, in terms of how to set them and how sensitive the performance of BACON is to their values.\n* MMD is not mentioned in the text despite being plotted in Figure 1(d).\n* Table 3 is not referenced in Section 6.4.\n* Please include up/down arrows in the tables next to the columns so that it is easy to read them.\n* Appendix B seems incomplete.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their positive and insightful feedback. We will revise the text addressing the issues raised. In addition to our global response in the main author rebuttal, our response to the reviewer's specific questions follows below.\n\n### Weaknesses\nRegarding the issues raised as weaknesses, we elaborate on answers below.\n> It would have been nice to see an ablation study where the design and parameters are not jointly optimized over, in order to ascertain the benefit of doing so.\n\nTo optimise only the calibration parameters while keeping the designs fixed usually involves setting them to the designs in the real data. However, the choice of designs to use for each chosen calibration parameter in the batch, especially when the batch size and number of real data points do not match, leads to some ambiguities and not a very obvious pathway to implementation. Our current  proxy for the non-joint-optimisation case is VBMC, whose performance is presented in Table 1 to 3.\n\n> The error function $\\epsilon$ is also modelled as a GP. I wonder how does BACON's performance get affected if this assumption does not hold (i.e. when this error model is misspecified).\n\nIf the errors are still sub-Gaussian (i.e., its tail probabilities decay faster than that of a Gaussian), as is the case for bounded errors, we conjecture that a GP model can still be a suitable surrogate, as it would not underestimate the error uncertainty. If the error function is sampled from some form of heavy-tailed stochastic process (e.g., a Student-T process), the GP would, however, tend to under estimate uncertainty and lead to possibly optimistic EIG estimates that make the algorithm under-explore the search space. Changing from a GP model to another type of stochastic process model that can capture heavier tails would be possible, though require significant changes to the algorithm's predictive equations. We, however, believe that most real-world cases would present errors which are at least bounded (and therefore sub-Gaussian) with respect to the simulations. We will add a discussion to the revision.\n\n>Perhaps a discussion of the hyperparameters/settings of the proposed method would be nice to have, in terms of how to set them and how sensitive the performance of BACON is to their values.\n\nBesides the GP hyperparameters (e.g., lengthscales, noise variance, etc.), there are optimisation settings (i.e., step sizes, scheduling rates, etc.), conditional density model hyper-parameters (i.e., normalising flow architecture), and algorithm's settings, e.g., the designs batch size $B$. The latter is dependent on the available computing resources (e.g., number of CPU cores or compute nodes for simulations in a high-performance computing system). We tuned optimisation settings and architectural parameters for the conditional normalising flows via Bayesian optimisation with short runs (e.g., 10-30 iterations) on the synthetic problem. However, depending on the number of parameters, a simpler approach, like grid search, might be enough. The GP hyper-parameters were optimised via maximum a posteriori estimation after each iteration's batch update. We will elaborate on the details of our hyper-parameter optimisation setup in the revision.\n\nThanks for noticing the other minor issues, which we will make sure to address in the revised text. MMD stands for maximum mean discrepancy, an integral probability metric that quantifies the distance between two probability distributions based on their embeddings in a reproducing kernel Hilbert space. We used the implementation available in the ITE package [39]. We will add a proper definition and reference in the revision.\n\n### Answers to questions\n**Tightness of the variational EIG bound.**\nThe bound in Sec. 5.1, quantifying the difference between the variational EIG and the true EIG, is exact. A full derivation can be found in Foster et al. [13, Appendix A.1], but we will add it to the revision, as suggested.\n\n**VBMC's performance** VBMC relies on a mixture of Gaussians as the variational approximation to the posterior, which tends to be quite smooth. Therefore, the algorithm may struggle to approximate non-smooth posteriors, as is the case for the location finding problem, for example. Another potential point of struggle for VBMC is in problems where the posterior's mode is not concentrated, but spread along a path, e.g., a circle, as shown in Oliveira et al. [22]. The dimensionality of the calibration parameters space may also affect VBMC's posterior approximation, since Gaussian components can be seen as a form of radial basis functions, whose representation power deteriorates with the increase in dimensionality. Another key difference to note is that VBMC does not produce an emulator approximating the simulator, but an approximation to the log-likelihood function, in contrast to the other GP-emulation-based methods we assess.\n\n**KL divergence performance criteria.** The expected information gain is equivalent to the expected KL divergence between the posterior and the prior (Eq. 1) and quantifies the expected reduction in uncertainty (i.e., entropy) on the posterior after incorporating observations produced by a given set of designs. The average KL in the tables, therefore, can be seen as an approximation to the final EIG. In our setup, we are approximating a latent function representing the simulator, and as the number of observations grows to infinity, that approximation should asymptotically concentrate at the true simulator. The posterior obtained with the true simulator is simply $p^*$. However, there might be a long gap between an optimal EIG and its asymptotic limit, which might not be possible to quantify. Computing the KL divergence between the final posterior approximation $p_T$ and $p^*$ is then a more direct measure of how close the algorithm has got to the asymptotic posterior.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "XoP8VQyts1", "reviewer": "Reviewer_bPdt", "summary": "The paper proposes a more data-efficient algorithm inspired by Bayesian adaptive experimental design. This algorithm runs maximally informative simulations in a batch-sequential process, estimating posterior distribution parameters and optimal designs by maximizing a variational lower bound of the expected information gain. The algorithm is validated on both synthetic and real-data problems.", "questions": "1. In Figure 1:  (a) It's challenging to ascertain if stability has been achieved, which affects the credibility of the conclusion that BACON achieves rapid convergence in terms of MAP estimates. (b). It's difficult to determine if EIG statistically outperforms others in terms of RMSE.\n2. The paper lacks clarity in specifying its contributions and novelty. Could you please elaborate on the main distinction between BACON and other Bayesian optimal design methods that incorporate variational inference (e.g. Variational Bayesian Optimal Experimental Design)?\n3. There are numerous instances where variational inference (VI) can fail and produce poor approximations of the target posterior distribution. It would be beneficial to investigate the performance of BACON in scenarios where VI struggles to capture the characteristics of the target distribution, such as when dealing with highly correlated coordinates or multimodal targets.", "limitations": "Variational inference (VI) is typically used to approximate posterior distributions in large-scale datasets and high-dimensional parameter spaces. However, the paper indicates that BACON, which uses VI in Bayesian optimal design, still struggles to scale to large datasets. This raises the question of why VI is used in BACON instead of traditional MCMC sampling methods like HMC. The rationale for choosing VI over these traditional methods remains unclear.", "rating": 3, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The paper is well-written, offering comprehensive background information and a thorough review of the literature. The method is rigorously compared with other related approaches across multiple metrics on both synthetic and real-data datasets.\n2. The paper focuses on a well-motivated and challenging problem.", "weaknesses": "The paper can be improved in the following ways:\n1. My biggest concern for this paper is its novelty: replacing EIG by a variational lower bound has been explored and well studied by many literature. (e.g. [18], [29])\n2. The method is currently compared only against \"Random\" and \"IMSPE\". However, there exist numerous other variants of Bayesian optimal design and frequentist approaches with diverse optimality criteria. It would be valuable to compare the proposed method with these state-of-the-art alternatives.\n3. The variational inference is usually applied to sampling from posterior with large-scale datasets or high-dimensional parameter spaces. But the paper only present results on small datasets with low-dimensional parameter spaces, which makes it less convincing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for the detailed comments and insightful feedback. We provide a global response clarifying the contrast with the state-of-the-art and details on an additional baseline, but further elaborate on specific details relevant to the reviewer's comments below.\n\n**Novelty.** Concerning novelty, we highlight that it is not trivial to apply existing Bayesian adaptive experimental design (BAED) methods to the calibration setting with an expensive black-box simulator. Existing adaptive experimental design approaches (e.g., Foster et al. [13, 16]) rely on \\emph{unrestricted} access to the simulator $h$ to generate simulated observations $\\hat{y}$, while in our case we have a limited budget due to the cost of simulations (e.g., long runtime). An alternative would be to sample observations from the GP model, instead, and perhaps learn an adaptive design policy based on GP-based observation histories [16, 17]. The GP then would serve as a proxy for the simulator and the actual simulations would serve as the ``real'' experiments within a typical BAED framework. However, sampling from a GP breaks important conditional independence assumptions which these methods rely on, since simulation outcomes $\\hat{y}_1, \\dots, \\hat{y}_t$ sampled from a GP would no longer be conditionally independent given the true calibration parameters $\\boldsymbol{\\theta}^*$, as $p(\\hat{y}_t | \\boldsymbol{\\theta}^*, \\mathbf{\\hat{x}}_t, \\boldsymbol{\\hat{\\theta}}_t)$ $\\neq$ $p(\\hat{y}_t | \\boldsymbol{\\theta}^*,  \\mathbf{\\hat{x}}_t, \\boldsymbol{\\hat{\\theta}}_t, \\mathcal{D}\\_{t-1})$\nunder a GP prior for a latent simulator $h$. In addition, the most appropriate prior for these methods would not be the unconditional prior $p(\\boldsymbol{\\theta}^*)$, but $p(\\boldsymbol{\\theta}^*|\\mathbf{y}_R)$, instead, which is non-trivial.\n\n**Contrast with the state of the art.** Besides the reasons outlined above, one could think of Bayesian optimal experimental design methods that consider black-box simulators [14, 19]. However, these approaches often only consider the calibration space [19], keeping the designs fixed, or are focused on the problem of finding a single optimal design for a real experiment [14]. We chose VBMC as a representative for approaches which focus on the choice of calibration parameters, while keeping the design points fixed. VBMC's performance is presented in Table 1 to 3 in the manuscript in addition to the other two baselines.\n\n**Additional baseline.** For the rebuttal, we are including experimental results with an additional baseline representing a classic experimental design criterion. We implemented a strategy which selects designs of maximal predictive entropy for a GP model with only the simulation data. Such strategy can be related to the D-optimality criterion when considering the stochastic process representing the unknown simulator as the random variable of interest. We include a table with the new results in the PDF as part of our author rebuttal.In additional, we also show above that the IMSPE is A-optimal and equivalent to the Active Learning Cohn (ALC) criterion used in other active learning methods (Sauer et al., 2022, in Author Rebuttal above). Our new results show that BACON is still superior in terms of maximising the expected information gain (EIG), as measured by the KL divergence between the final (MCMC) posterior estimate $p_T$ and the true posterior $p^*$ (with perfect knowledge of the simulator).\n\n**Reason for variational inference.** We chose a variational approximation as it allows us to formulate a lower bound on the expected information gain (EIG) and perform *joint optimisation and inference* with a single well-defined objective. This is perhaps more difficult to achieve through sampling.  VI for us is then simply a framework to reformulate both inference and design problems as optimisation. Despite the similarities, however, note that our formulation works by minimising the *forward* KL divergence $KL(p||q)$, instead of the usual reverse $KL(q||p)$ which VI works with. The forward KL tends to match the moments of the target distribution better, while reverse KL mainly focuses on the modes. As such, it would behave differently in the situations which the reviewer described. On multimodal targets, for example, an optimal forward-KL approximation would tend to spread its mass across the detected modes. Therefore, even if the approximation is limited due to the representation power of the chosen parametric family (i.e., conditional normalising flows in our case), the uncertainty on the parameters estimate can still be captured, which is the main driver for an information-theoretic criterion, like the EIG. In addition, a reverse KL objective could lead to certain pathologies that cause discontinuities in the EIG estimation, due to the variational posterior possibly locking to very different modes (see Appendix G in Foster et al. [13]), which we avoid with the forward KL.\n\n**Performance in Figure 1.** A large part of the variance of the results in Figure 1 and Table 1 is actually due to the variability of the GP samples used as simulators for each individual run. The average behaviour of the metrics (as presented by the solid mean curves), therefore, is more indicative of the actual performance of the algorithms. In that case, we see that the MAP error drops quicker for BACON than for the other algorithms, as well as the RMSE values.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "4MF4n2QNvL", "reviewer": "Reviewer_CjgP", "summary": "This paper addresses the problem of calibrating simulation models. Simulation models depend on inputs set by the user, referred to as designs, and parameters representing unknown physical quantities, called calibration parameters. The task is to find calibration parameters such that simulations match real observations. To that end, the authors propose an active learning scheme in which maximally informative design and calibration parameters are iteratively chosen to construct the training set. They also propose a Gaussian process structure adapted to this setting.", "questions": "I have no questions.", "limitations": "I do not see any unaddressed limitations.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "### Originality\n* This piece of work is new to me.\n\n### Quality\n* The method is sound, and I did not identify any flaws.\n* Related works are discussed.\n* The experiments support the claims.\n\n### Clarity\n* The paper is well articulated.\n* The problem is clearly introduced, and the intuition behind the proposed solution is provided early on.\n* Figures are clean. \n\n### Significance\n* Although my knowledge of the field is too limited to have a strong opinion, the problem addressed seems important.", "weaknesses": "### Originality\n* I have no concerns regarding originality.\n\n### Quality\n* I have no concerns regarding the quality.\n\n### Clarity\n* Calibration parameters are also referred to as simulation parameters, which is confusing.\n* Some minor comments for the camera-ready version. There is a typo in equation 6 (missing parenthesis). In algorithm 1, in the \"update posterior\" line, I believe this should be $\\mathcal{D}_t$. Table 3 is not referred to in section 6.4.\n\n### Significance\n* I have no concerns regarding significance.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their positive feedback and for letting us know about the clarity issues. Table 3 consists of the results discussed in Sec. 6.4, which unfortunately missed a direct reference to the table. We will make sure to address these and the other issues in the revision.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We would like to thank all the reviewers for their constructive feedback and the time and effort applied in reviewing our manuscript. We provide individual responses to each review, but we also address some of the main common points here. In addition, we have *new results* with additional baselines available in the PDF attached to this rebuttal.

# Applicability of existing BAED methods
In contrast to traditional Bayesian adaptive experimental design (BAED) approaches, the Bayesian calibration problem presents a few key differences and challenges.
* **The experiment *is* the simulation.** Our problem is to select informative simulation designs $(\mathbf{\hat{x}}, \boldsymbol{\hat{\theta}})$ to run in order to reduce uncertainty about calibration parameters $\boldsymbol{\theta}^*$ (e.g., physical properties) that influence real observed data. The possibly expensive simulations, therefore, are the experiments we run, and the role of the real data is to form our prior beliefs on the unknown parameters of interest, which are then updated using simulations.
* **A model for the simulator (a.k.a. an emulator) is needed.** To run state-of-the-art BAED approaches, we would need an emulator for the simulator that is able to correlate the simulation outcomes with the unknown parameter $\boldsymbol{\theta}^*$ for the real data. In our case, we use a GP model, which is a natural choice for Bayesian emulation, given its well known capabilities and guarantees. However, a GP also models correlations across the outcomes, which breaks the conditional independence assumption across simulation outcomes (e.g., $\hat{y}_i, \hat{y}_j$), since: $$p(\hat{y}_i, \hat{y}_j|\boldsymbol{\theta}^*, \mathbf{y}_R) = p(\hat{y}_i| \hat{y}_j, \boldsymbol{\theta}^*, \mathbf{y}_R) p(\hat{y}_j|\boldsymbol{\theta}^*, \mathbf{y}_R)\neq p(\hat{y}_i|\boldsymbol{\theta}^*, \mathbf{y}_R)p(\hat{y}_j|\boldsymbol{\theta}^*, \mathbf{y}_R),$$ an important assumption which existing methods rely on. Ignoring it would lead to an over-estimation of the predictive variances of $\hat{y}$ by using the GP predictive marginal $p(\hat{y}_i|\boldsymbol{\theta}^*, \mathbf{y}_R)$, instead of the conditional $p(\hat{y}_i| \hat{y}_j, \boldsymbol{\theta}^*, \mathbf{y}_R)$, besides also shifting the predictive mean, leading to biases in the EIG and misguiding the algorithm.
* **The "prior" is non-trivial.** The prior $p(\boldsymbol{\theta}^*|\mathbf{y}_R)$ for the calibration problem is actually a posterior over the parameters $\boldsymbol{\theta}^*$ given the real data $\mathbf{y}_R$. In contrast, standard BAED approaches assume the prior is relatively simple to sample from and evaluate densities of, possibly with a closed-form expression. So we also cannot rely on easily evaluating prior densities or on sampling from them, except for approximations provided by, e.g., MCMC or VI.

Considering the points above, the application of existing BAED methods to our setting is not straightforward and led us to the development of a dedicated algorithmic framework to tackle these issues.

# Additional baseline
We provide experimental results with an additional baseline following a *D-optimality* criterion, a classic experimental design objective, in the rebuttal PDF. Optimal candidate designs according to this criterion are points of maximum uncertainty according to the model [28]. If we model the simulator $h$ as the unknown variable of interest, this corresponds to selecting designs where we have maximum entropy of the Gaussian predictive distribution $p(\hat{y}|\mathbf{\hat{x}}, \boldsymbol{\hat{\theta}}, \mathbf{\hat{y}}_t)$. This approach, therefore, simply attempts to collect an informative set of simulations according to the GP prior only, without considering the real data. Running D-optimality on $\boldsymbol{\theta}^*$, in contrast, would lead back to the EIG criterion we use.

We would also like to highlight that the IMSPE criterion, which we use as a baseline, is equivalent to the Active Learning Cohn (ALC) criterion (Sauer et al., 2022, below) and a form of *A-optimality* [28], another classic experimental design objective. The ALC maximises the reduction in predictive variance of the model across a set of reference designs. Since the current predictive variance is constant w.r.t. the new data, maximising the ALC is equivalent to minimising the IMSPE. In addition, an integral over the predictive variance is also evaluating the trace of the posterior covariance operator of the simulations GP. Namely, letting $\Sigma := \mathbb{E}[h \otimes h ]$ denote the covariance operator of $h\sim \mathcal{GP}(0, k)$, assuming $h$ is realised in an $L_2(\Xi)$ space of square-integrable functions of $\xi := (\mathbf{\hat{x}}, \boldsymbol{\hat{\theta}}) \in \Xi := \mathcal{X} \times \Theta$, we have that:
$$\mathrm{Tr}(\Sigma) = \mathrm{Tr}(\mathbb{E}[h \otimes h ]) = \mathbb{E}[\mathrm{Tr}(h \otimes h) ] = \mathbb{E}[\lVert h \rVert_2] = \mathbb{E}\left[\int_{\Xi} h(\xi)^2 d\xi \right] = \int_{\Xi} \mathbb{E}\left[ h(\xi)^2 \right] d\xi = \int_{\Xi} \sigma^2(\xi) d\xi,$$
where we applied the linearity of the trace and of the expectation. A-optimality seeks designs $\xi$ which maximally reduce the trace of the posterior covariance of the unknown variable, i.e., maximising the difference $\mathrm{Tr}(\Sigma) - \mathbb{E}[\mathrm{Tr}(\Sigma)|\xi, \hat{y}]$. Therefore, maximising the ALC (or minimising the IMSPE) is equivalent to selecting A-optimal points. We revised our IMSPE implementation based on the methodology presented in Sauer et al. (2022) and included the results in the PDF.

Our new results show that BACON can still reach superior performance when compared to classic experimental design criteria in terms of expected information gain (EIG), while keeping competitive performance in terms of posterior approximation.

**References**

* Sauer, A., Gramacy, R. B., \& Higdon, D. (2022). Active Learning for Deep Gaussian Process Surrogates. *Technometrics, 65*(1).

Meta-review: This paper introduces a method for calibrating expensive computer models while jointly optimizing design parameters using Bayesian adaptive experimental design. The approach employs expected information gain for calibration and performs active learning in the joint space of design and parameters.
The reviewers generally agree that the paper addresses an important problem in engineering and physical sciences. Two reviewers recommend acceptance. One reviewer gives a borderline accept, while another recommends rejection. In the final discussion, however, reviewers highlighted the paper as technically solid and well-written as well as novel. 
The authors are encouraged to incorporate their rebuttal explanations and additional results into the final version of the paper.
Decision: Accept (poster)