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.

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.

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 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})$ under 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.

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.

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).

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.

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.

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.

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.

Weaknesses

Regarding the issues raised as weaknesses, we elaborate on answers below.

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.

To 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.

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).

If 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.

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.

Besides 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.

Thanks 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.

Answers to questions

Tightness of the variational EIG bound. The 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.

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.

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.

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.

• 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.
• Finding optimal design points in higher dimensions may become increasingly difficult due to the vast spaces.

The 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.

The 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.

In 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.

References

• 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.
• 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).
• 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).