Bayesian Adaptive Calibration and Optimal Design

Authors: rafael  oliveira
dino  sejdinovic
david  howard
edwin v bonilla
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 05:28:04

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.

Main paper: 
 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 is scarce. These simulation models depend on the inputs set by the user, commonly referred to as designs, and on a 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 as the calibration of computer models.The seminal work by Kennedy and O'Hagan [1] 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 in the calibration parameters. While the simulator is an essential tool when obtaining real data is expensive or unfeasible, each run of a simulator may itself involve significant computational resources, especially in applications such as climate science or complex engineering systems. In this situation, it is imperative to run simulations at carefully chosen settings of designs as well as of calibration inputs, using current knowledge to optimise resource use [3] [4] [5] .In this contribution, we bridge Bayesian calibration with adaptive experimental design [6] and use information-theoretic 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 to prior work, it optimises designs jointly with calibration inputs in order to capture informative correlations across both spaces. Experimental results on synthetic experiments and a robotic gripper design problem demonstrate the benefits of BACON compared to competitive baselines in terms of computational savings and the quality of the estimated posterior under similar computational constraints. 
problem formulation: 
 Let f : X → Y represent a mapping of experimental designs x ∈ X to the outcomes of a physical process f (x) ∈ Y ⊂ R. We are given a set of observed outcomes y R = [y 1 , . . . , y R ] T and their associated designs X R := {x i } R i=1 ⊂ X . Observations are corrupted by noise as y i = f (x i ) + ν i , where ν i ∼ N (0, σ 2 ν ) is zero-mean Gaussian noise, for i ∈ {1, . . . , R}. In addition, we have access to the output of a computer model h : X × Θ → R given a design input and simulation parameters. Given an optimal setting for the calibration parameters θ * ∈ Θ, the simulator h(x, θ * ), can be used to approximate the outcomes of the real physical process f (x). However, θ * is unknown, and evaluations of the simulator h are costly, though cheaper than executing real experiments evaluating f . Our task is to optimally estimate θ * given the real data y R , outputs of the simulator h and a prior distribution p(θ * ), representing initial assumptions about θ * . More concretely, let ŷS := [h(x i , θi )] S i=1 represent simulated outcomes for a set of designs X S := {x i } S i=1 ⊂ X and simulation parameters Θ S := { θi } S i=1 ⊂ Θ. Given the cost of running simulations, we will associate the simulator h with a latent function (usually referred to as emulator) drawn from a Gaussian process (GP) prior and assume simulation outputs and real data follow a joint probability distribution p(y R , ŷS , θ * ).In this setting, the Bayesian experimental design objective is to propose a sequence of simulations which will maximise the expected information gain (EIG) about θ * :EQUATIONwhere H(•) represents the entropy of a probability distribution, D KL (•||•) denotes the Kullback-Leibler divergence, and I(θ * ; ŷS | y R ) is the mutual information between θ * and the simulator output ŷS given the real observations y R and the simulator inputs to be optimized. We note here that, in our setting, the real observations y R are always fixed. Therefore, intuitively, the EIG above captures the reduction in uncertainty that will be obtained when selecting ( X S , Θ S ) averaged over all the possible outcomes ŷS . 
related work: 
 Our work consists of deriving a Bayesian adaptive experimental design (BAED) approach to the problem of calibration. Therefore, in the following, we will briefly discuss current literature on these two main research areas.The problem of experimental design has a 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 a 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 a 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 a computational bottleneck for information-theoretic 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 is known and inexpensive to evaluate. In contrast, the simulations themselves are modelled as expensive experiments for us, and we apply Gaussian process models as emulators to capture uncertainty over the black-box simulator. In addition, traditional BAED approaches assume the prior is trivial to sample from and evaluate densities of, while in our case the starting prior is p(θ * |y R ), 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.Experimental design approaches generally aim towards the selection of designs for physical experiments, 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 to learn GP emulators for a 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 based on selecting 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 design points, keeping them fixed.Allowing for design point decisions to be included, Leatherman et al. [4] presented approaches for combined simulation and physical experimental design based on geometric and prediction-error-based criteria, but using an offline, non-sequential framework. More recently, Marmin and Filippone [5] derived a deep Gaussian process [27] framework for Bayesian calibration problems and discussed an application to 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.r.t.) variance parameters. In contrast, we aim to directly maximise the information gain w.r.t. the unknown calibration parameters. 
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 θ * . Ideally, the model needs to allow for a computationally tractable conditioning on the parameters θ * and account for the differences between real and simulated data. Hence, we follow the Bayesian calibration approach in Kennedy and O'Hagan [1] and model:EQUATIONwhere ε : X → R represents the error (or discrepancy) between simulations and real outcomes, and ρ ∈ R accounts for possible differences in scale. We place Gaussian process priors on the simulator h ∼ GP(0, k) and on the error function ε ∼ GP(0, k ε ).Since both h and ε are GPs, simulations and real outcomes can be jointly modelled as a 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 a coarser version of f . Let s ∈ S := {0, 1} denote a fidelity parameter. The combined model is then given by:EQUATIONsuch that f (x) = f (x, θ * , 1) and h(x, θ) = f (x, θ, 0), for any x, x ∈ X and θ ∈ Θ. As a result, for arbitrary points in the joint space z, z ′ ∈ Z := X × Θ × S, the following covariance function parameterises the combined GP model f ∼ GP(0, k):EQUATIONwhere k ρ (s, s ′ ) := (1 + s(ρ -1))(1 + s ′ (ρ -1)), z := (x, θ, s), and z ′ := (x ′ , θ ′ , s ′ ). Therefore, any set of real and simulated evaluations are joint normally distributed under a combined GP model.Let Z R := Z R (θ * ) := [(x i , θ * , 1)] R i=1 represent the set of partially observed inputs for real data y R , and let Z S := [(x i , θ, 0)] S i=1 denote the current set of simulation inputs for the observations ŷS . Under the GP prior, the joint probability model p(ŷ S , y R , θ * ) can be decomposed as:EQUATIONwhere f := f (Z(θ * )) ∈ R R+S , and Z(θ * ) := {Z R (θ * ), Z S } corresponds to the full set of inputs.The GP prior then allows us to model real and simulated outcomes jointly as a Gaussian random vector f :EQUATIONK(θ * ) := k(Z(θ * ), Z(θ * )) = [k(z, z ′ )] z,z ′ ∈Z(θ * )denotes the prior covariance matrix. Assuming a Gaussian noise model for the observations y = f (x, θ * ) + ε(x) + ν, with ν ∼ N (0, σ2 ν ), the marginal distribution over the observations y := [y T R , ŷT S ] T is available in closed form as:EQUATIONwhere Σ y denotes the covariance matrix of the observation noise, i.e., [Σ y ] ii = σ 2 ν for any z i with s i = 1, and [Σ y ] ij = 0 elsewhere. 2 Under the GP assumptions, we can make predictions about ŷ = h(x, θ) at any pair of x, θ ∈ X × Θ. Conditioning on θ * and a datasetD t := {X R , y R , X t , Θ t , ŷt }, let Z t (θ * ) := {Z R (θ * ), Z t }denote the set of inputs up to time t conditional on θ * , and y t the corresponding outputs. We then have that:EQUATIONfor ẑ := (x, θ), where:EQUATIONEQUATIONwith k t (ẑ; θ * ) := k(Z t (θ * ), ẑ) and K t (θ * ) := k(Z t (θ * ), Z t (θ * )). We next describe how to apply this model to derive a Bayesian adaptive calibration algorithm.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: Code is included, though soft-robotics 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 a 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 it 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 to 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 a 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 submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or 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 reproducibility.In the case of closed-source models, it 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 to have some path to reproducing or verifying the results. 
bayesian adaptive calibration and optimal design: 
 In this section, we describe an approach to design experiments for calibration of computer models that incorporates information gathered during the experiments iteratively. We refer to these types of designs as adaptive. Thus, we consider the sequential design of experiments setting, where at each iteration t ∈ N, we optimise:EQUATIONgiven the dataset D t-1 := {X R , y R , X t-1 , Θ t-1 , ŷt-1 } of observations. Given that the expected information gain is submodular [28] , a sequential approach allows us to get close enough (usually a factor of at least 1 -1/e [29] ) to the optimal EIG over the whole experiment, while also allowing algorithmic decisions to adapt to current estimates for p(θ * |D t ).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(θ * |ŷ, x, θ, D t-1 ) and p(ŷ, θ * |x, θ, D t-1 ) depend on the posterior p(θ * |D t-1 ), as:EQUATIONEQUATIONwhere the conditional predictive density p(ŷ|θ * , x, θ, D t-1 ) is Gaussian and available in closed form (Eq. 8). Clearly, in general, the true posterior is intractable, since p(θ * |D t ) = p(Dt|θ * )p(θ * ) p(Dt) and p(D t ) = Θ p(D t |θ * )p(θ * ) dθ * involves integration over the entire parameter space Θ, which can be high dimensional and passed through highly non-linear operations such as inverse covariances. In addition, the marginal predictive p(ŷ|x, θ, D t-1 ) = Θ p(ŷ, θ * |x, θ, D t-1 ) dθ * is usually also intractable for the same reasons.Following Foster et al. [13] , we replace the EIG by a variational objective which does not directly involve the true posterior over θ * . This formulation allows us to jointly estimate an approximation to the posterior and select optimal design points x and simulation parameters θ. Applying the variational lower bound by Barber and Agakov [30] to Eq. 12 yields the following alternative to the EIG:EQUATIONwhere q(θ * |ŷ, x, θ) is any conditional probability density model. The gap is given by the expected Kullback-Leibler (KL) divergence between the true and the variational posterior [13]:3 EQUATIONMaximising the variational EIG lower bound w.r.t. the variational distribution q then provides us with an approximation to p(θ * |ŷ, x, θ, D t-1 ). Therefore, we can simultaneously 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:EQUATIONgiven a suitable variational family Q of conditional distributions.Algorithm 1 summarises the method we propose for Bayesian adaptive calibration and optimal design (BACON). The algorithm is given an initial dataset D 0 containing the real data (and possibly previously available simulation data) and an estimate of the posterior given the real datap 0 p(θ * |X R , y R ).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 real dataD 0 := {X R , y R }.Note that we only need samples from the previous posterior to estimate the expectation in Eq. 17, with no need to directly evaluate their probability densities. Each iteration starts by optimising the variational EIG lower bound using the objective in Eq. 17 to jointly select an optimal design xt , simulation parameters θt and variational posterior q t . Given the new design xt , we run the simulation with the chosen parameters θt , observing a new outcome ŷt . The calibration posterior p t (θ * ) and 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 p T (θ * ) = p(θ * |y R , ŷT ) and a conditional density model q T are obtained.Algorithm 1 BACON input D 0 := {X R , y R }; {Real data} input p 0 (θ * ) := p(θ * |D 0 ) {MCMC or VI prior distribution} for t ∈ {1, . . . , T } do xt , θt , q t ∈ argmax x, θ,q E pt-1(ŷ,θ * |x, θ) [log q(θ * |ŷ)] ŷt := h(x t , θt ) {Run simulation} D t := D t-1 ∪ {x t , θt , ŷt } {Update GP model} p t (θ * ) = p(θ * |D t-1 ){Update posterior via MCMC or VI} end for output p T (θ * ) {Final posterior}Any conditional probability density model q(θ * |ŷ) estimating probability densities over the parameter space Θ given an observation ŷ 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 we can approximate p(θ * |D t ) as a Gaussian, we may set:EQUATIONwhere m ϕ and Σ ϕ are given by parametric models, such as neural networks, with parameters ϕ. To ensure Σ ϕ (•) is positive-definite, it can be parameterised by its Cholesky decomposition T , where L ϕ (•) is a lower-triangular matrix with positive diagonal entries.Σ ϕ (•) = L ϕ (•)L ϕ (•)Conditional normalising flows Normalising flows [33] apply the change-of-variable formula to derive composable, invertible transformations g w of a fixed base distribution p 0 :EQUATIONThe log-probability density of a point ξ = g w (ξ 0 ) under this model can be calculated as:log p K (ξ; w) = log p 0 (ξ 0 ) - K j=1 log J (j) w (ξ j-1 ),where ξ 0 := g -1 w (ξ), ξ j := g (j)w (ξ j-1 ), and J (j) w is the Jacobian matrix of the jth transform g(j)w , for j ∈ {1, . . . , K}. Several invertible flow architectures have been proposed in the literature, including radial and planar flows [33] , autoregressive models [34] [35] [36] and models based on splines [37] .To derive a conditional density model q ϕ (θ * |ŷ), conditional normalising flows map the original flow parameters w via a neural network model r ϕ : ŷ → w [38, 39] . In this case, the variational model is given by:EQUATIONOften simulations can be run in parallel by spawning multiple processes in a single machine or over a compute cluster. In this case, proposing batches of simulation inputs can be more effective than running single simulations in a sequence. Optimising the EIG w.r.t. a batch of inputs B := {x i , θi } B i=1 , instead of single points, we obtain a batch version of Algorithm 1. In this case, we are seeking a batch that maximises the mutual information between the parameters θ * and the resulting observations, i.e.:EQUATIONWe approximate this objective by using variational models that accept multiple conditioning observations q(θ * |ŷ 1 , . . . , ŷB ). In the case of scalar observations, this simply amounts to replacing the scalar inputs to the conditional models in Sec. 5.3 by vector-valued inputs. Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? 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. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
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 experiment details can be found in Appendix A and in our code repository. 4Performance metrics. We evaluated each method against a set of performance metrics, which we now describe. The maximum-a-posteriori (MAP) approximation error measures the distance between the mode of the variational distribution and the true parameters θ * . 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 := 1 N N i=1 (E q(θ) [µ(x * i , θ * ; θ)] -y * i ) 2, where y * i = f (x * i ) + ν * i are observations of the true function over a set of designs {x * i } N i=1 ⊂ X placed on a uniform grid the design space.Information gain. Lastly, we also evaluated two sample-based estimates of the KL divergence [40] . Namely, D KL (p T ||p 0 ) 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 an initial set of randomised simulations) both estimated over the learnt GP model. ) corresponds to the KL divergence between the final posterior (estimated after each algorithm's run with all the data it collected) and the starting one (higher is better), while D KL (p T ||p * ) is the KL between the final posterior and the posterior with full knowledge of the simulator p * (lower is better). All posteriors were sampled via MCMC using 4000 samples. Averages and standard deviations were estimated from 10 independent runs.Our algorithmic baselines were chosen to illustrate the main approaches currently available in the literature. Both are employed as adaptive baselines, in the sense that their GP models are updated with the latest observations before proceeding to the next iteration.Random search. This baseline samples simulation designs xt ∼ U(X ) from a uniform distribution over the design space X and calibration parameters from the prior θt ∼ p(θ * ).IMSPE with MAP estimates. The integrated mean squared prediction error (IMSPE) [41] criterion chooses designs xt and calibration θt parameters by minimising the GP prediction error:EQUATIONThe posterior's MAP estimate θ * t ∈ argmax θ p(θ|D t-1 ) is used as a point estimate for the true θ * . The integral is approximated as a sum over a uniform grid of designs and samples from the calibration prior, 5 making IMSPE equivalent to active learning Cohn [42] and a form of A-optimality [28] . D-optimal designs. We provide experimental results with an additional baseline following a Doptimality criterion, a classic experimental design objective. Optimal candidate designs according to this criterion are points of maximum uncertainty according to the model [28] . If we 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(ŷ|x, θ, D t-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 θ * , in contrast, would lead back to the EIG criterion we use.Variational Bayesian Monte Carlo (VBMC). Acerbi [43] presents an adaptive Bayesian quadrature method to learn posterior distributions over models with black-box likelihood functions. The method estimates the posterior p(θ * |y R , h) by modelling the log-joint log p(y R , θ * |h) as a sample from a Gaussian process. VBMC then learns a variational posterior approximation by maximising a lower-confidence bound over the ELBO given by the GP estimates. Calibration parameter queries θt are obtained by optimising quadrature-based acquisition functions. Regarding design points, simulations are always run on the set of real design points X R in the observed data, which is fixed.For this experiment, we sampled a function f ∼ GP(0, k) to use as our simulator and compared different algorithms. Following a sparse GP approach, a function sampled from a GP can be 2 : Results on the location finding problem after T = 30 iterations with B = 4, R = 20 "real" data points and an initial set of 20 simulations. Estimates were averaged over 10 independent runs.D KL (p T ||p 0 ) ↑ D KL (p T ||p * ) ↓ BACON 0.approximated as f (z) ≈ k(z, Z M )K -1 M u M , where u M ∼ N (û M , Σ M ) is a sample from an M - dimensional Gaussian, Z M := {z i } M i=1 ⊂ X × Θ × {0, 1}, for a given M . As the number of points M → ∞, if the pseudo-inputs Z M form a dense set, we have that f converges in distribution to a sample from the Gaussian process GP(0, k). In our case, to sample Z M , we sample designs from a uniform distribution over the design space, calibration parameters from the prior, and fidelities from a Bernoulli distribution with parameter set to 0.5. We also set ûM := 0 and Σ M := K M = k(Z M , Z M ). We repeatedly run a loop of T iterations for each algorithm, with each repetition running on independent f sampled from the same GP prior.We run each algorithm for T := 50 iterations using a batch of B := 4 designs per iteration. Each of the methods using GP approximations for the simulator are initialised with 20 observations and R = 5 real data. To configure VBMC, we allow it to run an equivalent maximum amount of objective function evaluations. The design space is set as the 2-dimensional unit box X := [0, 1] 2 and the "true" parameters for each run are sampled from a standard normal prior p(θ * ) := N (θ * ; 0, I) also over a 2D space, totalling a 4-dimensional problem space.Fig. 1 shows a case where the GP-sampled simulator led to a 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 on 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-informative 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 σ ν = 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 to 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 a mixture of Gaussians, while BACON applies conditional normalising flows for its posterior approximations, which lead to increased generalisability. 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 approximates the real process), while VBMC's focus is on approximating the log-likelihood.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. Our results are presented in Table 2 , which show a similar tendency in higher information gain for our method, and a very low KL w.r.t. p * . Note that a high information gain indicates a 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 p 0 . Therefore, besides obtaining maximally informative data, we have shown that BACON is also efficient in approximating posteriors over black-box simulators, while also learning a GP emulator model. 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 was averaged over 4 independent runs.For this experiment, we are provided with a dataset containing R = 10 real measurements of the peak grasping force of soft-robotic gripper designs on a range of testing objects (see Fig. 3 ). The gripper designs follow a fin-ray pattern parameterised by 9 geometric parameters [44] , and we are interested in estimating 2 unknown physics parameters, the Young's modulus of elasticity and the coefficient of static friction with the objects. To simulate the gripper designs, we use the SOFA framework [45] to reproduce the grasping scenario and provide an estimate of the peak grasping force. In particular, for this paper, we focus on the grasping of a spherical object, which provides a simpler geometry and lower discrepancy with respect to real data measurements compared to more complex objects. This experiment provides us with a benchmark where simulations are expensive to run, taking from minutes to a few hours to run (depending on mesh resolution) on a high-performance computing platform. Therefore, it is important to choose a 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 to inferior posterior approximations (see Fig. 4 in the appendix) and showing the advantage of BACON's joint optimisation and inference.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: Main details are provided and code is included.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a 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, or as supplemental material. 
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. It does so by optimizing an information-theoretic criterion so that input designs and calibration parameters are selected to be maximally informative about the optimal parameters. Our method provides a full posterior over optimal calibration parameters as well as an accurate Gaussian process based estimation of the computer model (i.e., an emulator). One of the main limitations of the presented framework is scalability to large datasets, due to the cubic computational complexity of exact inference with GPs. However, our method can be extended to work with scalable sparse variational GP models [46] by using a conditional distribution model for the inducing points (see 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 multi-output observations by the use of multi-output GP models [47] . Further discussions on limitations and future work can be found in our appendix (see Appendix B and C).Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: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, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run 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). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It 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 a 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 of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 • 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. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS 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. 
 
 Appendix: D 0 := {X R , y R }; for t ∈ {1, . . . , T } do µ t-1 , k t-1 ← UpdateGP(D t-1 ) {θ * i } S A i=1 MCMC ∼ p(θ * |D t-1 ) ∝ p(θ * )N (y R ; µ t-1 (Z R (θ * ); θ * ), Σ t-1 (Z R (θ * ); θ * ) + σ 2 ν I) p(θ * |D t-1 ) ≈ pt-1 := 1 S A S A i=1 δ θ * i q t ← TrainFlow(p t-1 , D t-1 ) {x t,i , θt,i } B i=1 ← OptimiseDesigns(q t , pt-1 , D t-1 ) ŷt,i := h(x t,i , θt,i ) (parallel) for i ∈ {1, . . . , B} {Run batch of simulations} D t := D t-1 ∪ {x t,i , θt,i , ŷt,i } B i=1 {Update GP dataset} end forFor 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, 6 is based on Pyro probabilistic programming models [48] . Gaussian process modelling code is based on BoTorch7 [49] . 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 [50] was applied to obtain samples from p(θ * |D t-1 ) at each iteration t. KL divergences are computed from samples using a nearest neighbours estimator implemented in the information theoretical estimators (ITE) package8 [40] .The GP prior was set with k given by a squared exponential kernel and k ε 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 a permutation invariant set encoder, similar to Blau et al. [17] , with a 2-layer, 32-units-wide fullyconnected hyperbolic tangent neural network passing through a summation at the end. Gradient-based optimisation is run using Adam with a 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 were subsampled from the MCMC posterior to estimate expectations for both this and the location-finding problem.Algorithm with split training. For the synthetic GP problem, we provide a more detailed pseudocode of our algorithmic implementation using an option for training the conditional normalising flow and optimising the designs separately. Specifically, we applied MCMC to estimate our posteriors and had a 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 in 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 GP-based synthetic experiments, while for the other experiments we ran the full joint optimisation over both the simulation inputs (x, θ) and the variational parameters of the conditional model q.For this experiment we used more up-to-date Zuko9 implementations of the conditional normalising flow models, which were again set as neural spline flows [37] combined with a set encoder to condition on arbitrary batch sizes. Further architectural details can be found in our code repository. 256 samplesAlgorithm 3 TrainFlow input pt , D t for n ∈ {1, . . . , N } do { θi } B i=1 ∼ (1 -ϵ)p t + ϵp {x i } B i=1 ∼ U(X ) {θ * i } S i=1 ∼ pt {ŷ i,j } S,B i,j=1 ∼ N (µ t ({x i , θi } B i=1 ; {θ * i } S i=1 ), Σ t ({x i , θi } B i=1 ; {θ * i } S i=1 )) ϕ ← ϕ + η S S i=1 ∇ ϕ log q ϕ θ * i D t ∪ {x j , θj , ŷi,j } B j=1 end for output q ϕ Algorithm 4 OptimiseDesigns input q t , pt-1 , D t-1 Θ = { θi } B i=1 ∼ (1 -ϵ)p t-1 + ϵp X = {x i } B i=1 ∼ U(X ) for n ∈ {1, . . . , N } do {θ * i } S i=1 ∼ pt-1 {ŷ i } S i=1 ∼ N (µ t-1 ( X, Θ; {θ * i } S i=1 ), Σ t-1 ( X, Θ; {θ * i } S i=1 )) ( X, Θ) ← ( X, Θ) + η S S i=1 ∇ X, Θ log q t θ * i D t ∪ { X, Θ, ŷi } end for output {x i , θi } B i=1were subsampled from the MCMC posterior at each iteration to estimate expectations for EIG lower bound computations. The simulations kernel k was a 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 a standard deviation of 0.5. Final KL estimates were computed using the maximum-a-posteriori hyper-parameters of the GP model learnt with the random search approach to minimise biases in the estimate of D KL (p T ||p 0 ) due to differing GP hyper-parameters across baselines.The prior for the calibration parameters p(θ * ) in this experiment consisted of a 2-dimensional standard normal transformed through a sigmoid and an 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 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 to condition the flow. Adam again was used for optimisation with a 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 at each joint optimisation step.Besides the GP hyperparameters (e.g., lengthscales, noise variance, etc.), which had to be tuned for the non-GP-based 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-performance computing system). We tuned optimisation settings and architectural parameters for the conditional 3 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 x, which are 9-dimensional variables that are not plotted here.normalising flows via Bayesian optimisation with short runs (e.g., 10-20 iterations) on the synthetic problem. However, depending on the number of parameters, a simpler approach, like grid search, might be enough. GP hyper-parameters were optimised online via maximum a posteriori estimation after each iteration's batch update. Further implementation details can be found in our code repository. 10In the following, we 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 distribution models as priors, instead of having to run, for example, MCMC. Secondly, we present derivations for a scalable sparse GP version of our method.We use a conditional variational distribution model for q(θ * |ŷ). The main advantage of training a conditional model is that, once new data ŷt is observed, we readily obtain an approximation to the new posterior as p(θ* |D t ) = p(θ * |ŷ t , xt , θt , D t-1 ) ≈ q t (θ * |ŷ t ).There is, therefore, potential to reuse the variational posterior as the prior for the next iteration, and all the optimisation is concentrated within a single loop.Approximate objective. We are still left with terms dependent on the posterior from the previous iteration p(θ * |D t-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(ŷ, θ * |x, θ, D t-1 ) by an approximation using the previous optimal variational posterior q t-1 as:EQUATIONwhere q t-1 (θ * ) := q t-1 (θ * |ŷ t-1 ) ≈ p(θ * |D t-1 ). The following objective then approximately shares the same set of maximisers as the variational lower bound EIG t (x, θ, q):EQUATIONIn practice, reusing the variational conditional posterior may tend to degenerate the approximation over time. However, that can be corrected by rerunning MCMC or a variational inference scheme over the data to obtain a fresh new posterior at every few iterations.Computing the variational EIG requires evaluating expectations with respect to the posterior predictive distribution p(ŷ|θ * , x, θ, D t ). Note, however, that, as θ * appears inside a matrix inversion in the GP predictive (Eq. 8), each sample of p(ŷ|θ * , x, θ, D t ) requires a O(N 3 t ) computation cost, where N t := R + t is the number of data points at iteration t ∈ 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 to handle large amounts of data, we then need GP models that can reduce this computational complexity, while still allowing us to obtain reasonable EIG estimates.We consider an augmentation to the original GP model which allows us to sparsify its covariance matrix, reducing the computational complexity of GP predictions. Following the variational sparse GP approach [46] , let u := f (Z u ) ∈ R M denote a vector of M inducing variables representing unknown function values at a given set of pseudo-inputs Z u . The joint distribution between observations y, function values f := f (Z(θ * )), inducing variables u and the unknown parameters θ * can be written as:EQUATIONwhere p(y| f ) = N (y; f , Σ y ),EQUATIONand p(u) = N (u; 0, K uu ), using notation shortcutsEQUATIONK f f (θ * ) := k(Z(θ * ), Z(θ * )).We may now formulate an evidence lower bound (ELBO) based on the joint variational density q( f , u, θ * ) as:EQUATIONSince D KL (q( f , u, θ * )||p( f , u, θ * |y)) ≥ 0, and 0 if and only if q( f , u, θ * ) = p( f , u, θ * |y), maximising the ELBO above w.r.t. q provides us with an approximation to the joint posterior. Choosing q( f , u, θ * ) := p( f |u, θ * )q(u, θ * ) simplifies the ELBO to [51] :EQUATION)Sparse variational GP approaches can reduce the computational complexity of Bayesian inference on GPs to O(N M 2 ) or even O(M 3 ) [46, 52] , where N is the number of data points.If we would take a mean-field approach setting q(u, θ * ) := q(u)q(θ * ), the ELBO above would further simplify, leading to a few computational advantages, as explored by Bayesian GP-LVM methods [51, 32, 52] . 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 θ * , leading their mutual information (a.k.a. EIG) to be zero regardless of the design inputs x and θ. Secondly, although u and θ * are independent according to their priors (Eq. 25), they become dependent when conditioned on the data. In fact, the true posterior over u given the data and the true parameters θ * is exactly Gaussian:EQUATIONwhere µ t (•; θ * ) and k t (•, •; θ * ) are given by Eq. 9 and Eq. 10, respectively. Note, however, that the posterior over θ * should not be Gaussian for a general non-linear kernel k. Therefore, it makes more sense for us to model q(u, θ * ) := q(u|θ * )q(θ * ). Moreover, learning a Gaussian conditional model over u and a flexible variational distribution over θ * should be enough to allow us to recover the true posterior, since p(u, θ * |D t ) = p(u|D t , θ * )p(θ * |D t ).Optimal variational inducing-point distribution. Given θ * ∈ Θ, we have a standard sparse GP model. The optimal variational inducing-point distribution is available in closed form following standard results [46] as:EQUATIONwhere the distribution parameters are:EQUATIONand the conditional Ψ matrices are given by:EQUATIONfor θ ∈ Θ. The computational cost of sampling predictions with this model then reduces from O(N 3 ) to O(N M 2 ).Parametric variational inducing distribution. To further reduce the computational cost of predictions, we may accept a sub-optimal conditional variational inducing-point distribution given by a parametric model:EQUATIONfollowing the architecture in Sec. 5.3. This formulation allows us to approximate the evidence lower bound in Eq. 28 w.r.t. q(u|θ * ) via mini-batching [see 53]. To do so, we approximate fi := f (z i ) via conditionally independent samples given u, for i ∈ {1, . . . , N }. As a result, the data-dependent term in Eq. 28 decomposes as a sum which is amenable to mini-batching:EQUATIONwhere q ζ ( fi , u|θ * ) = p( fi |u, θ * )q ζ (u|θ * ). The variational parameters ζ need to be optimised within a second optimisation loop after the data update in Algorithm 1 w.r.t.:EQUATION)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 L ≪ N to Eq. 37 results in a computational cost O(LM 2 ) (or O(M 3 ), if M > L), which is smaller than the cost O(N M 2 ) of the optimal variational distribution q * (u|θ * ).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 a 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 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 is dominated by the cube of the number of data points due to a 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, if 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 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.• Theorems and Lemmas that the proof relies upon should be properly referenced.In practice, we add a small nugget term to the diagonal of the noise covariance matrix for numerical stability.We will at times write q(θ * |ŷ) to denote q(θ * |ŷ, x, θ) to avoid notation clutter, as it is implicit the dependence on the inputs (x, θ) through ŷ.Code available at: https://github.com/csiro-funml/baconThe original paper proposed analytic solutions to Eq. 22 tailored for specific kernels. However, we decided to keep our codebase generic to work with different kernels, and therefore opted for a numerical approximation.For VBMC, we used the author's Python implementation at: https://github.com/acerbilab/pyvbmcBoTorch: https://botorch.orgITE package: https://bitbucket.org/szzoli/ite-in-pythonZuko: https://zuko.readthedocs.io/stable/Code available at: https://github.com/csiro-funml/bacon 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.
questions: 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: limitations: Please also include limitation of the dimensionality issue as pointed out above.
paper_summary: paper_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.
strength_weakness: strength_weakness: I do not have any major concerns. Some questions/comments that might help improve the paper are as follows:
* 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.
* 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).
* 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.
* MMD is not mentioned in the text despite being plotted in Figure 1(d).
* Table 3 is not referenced in Section 6.4.
* Please include up/down arrows in the tables next to the columns so that it is easy to read them.
* Appendix B seems incomplete.
questions: 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)
* 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?
* 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: limitations: The authors have discussed the limitations adequately.
paper_summary: paper_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.
strength_weakness: strength_weakness: The paper can be improved in the following ways:
1. 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])
2. 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.
3. 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.
questions: 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.
2. 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)?
3. 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: 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: ### Originality
* I have no concerns regarding originality.

### Quality
* I have no concerns regarding the quality.

### Clarity
* Calibration parameters are also referred to as simulation parameters, which is confusing.
* 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.

### Significance
* I have no concerns regarding significance.
questions: questions: I have no questions.
limitations: limitations: I do not see any unaddressed limitations.

Comments: comment: Thanks for the feedback. We understand the reviewer's concern, as we perhaps have not selected the most appropriate performance metrics to show, and the discussion is lacking a few key insights. We would like to point out that most of the posteriors in the synthetic calibration problem of Sec. 6.2 are multimodal, since the simulators are simply random functions drawn from a GP prior and only 5 "real" data points were provided, alongside 20 initial simulation points. In these multimodal problems, MAP estimates, which show most of the confidence interval overlaps in Fig. 1, should not be considered a primary metric of performance, since the mode of the posterior will often not match the true calibration parameter, given the low amount of data.  We, however, have been able to show that our proposed method (BACON) achieves its intended goal, which is to maximise the EIG, as measured by the expected KL divergence from final to initial posterior $\mathbb{D}_{\mathrm{KL}}(p_T||p_0)$ (see Eq. 1 for the equivalence), when compared to the baselines in the paper across all experimental benchmarks. We will clarify these points in the revision and include a few examples of some of the posterior distributions we find in the synthetic and real data problems to better illustrate such challenges.
Thanks for the response. I have read other comments as well and decide to keep my scoring unchanged.
Thank you for your response. However, I still find Figure 1 unclear. The confidence intervals of the different methods overlap, suggesting that the proposed method may not be significantly better than the alternatives. As a result, my concerns remain, and I would like to retain my current score.
Thanks for the update.

I have read the rebuttal and the other reviews and keep my score unchanged.
I thank the authors for their clarifying responses to my questions. The general response outlining the difference between experimental design and Bayesian calibration makes the utility of their approach clearer. I think this paper is a useful contribution when calibrating parameters of a computationally expensive mechanistic model that depends on a design variable (which is slightly different to the kind of models calibrated using simulation-based inference). Hence, I am happy to recommend accept.