Simulating, Fast and Slow: Learning Policies for Black-Box Optimization

## Abstract

Abstract In recent years, solving optimization problems involving black-box simulators has become a point of focus for the machine learning community due to their ubiquity in science and engineering. The simulators describe a forward process f sim : ( 𝝍 , 𝒙 ) → 𝒚 : subscript 𝑓 sim → 𝝍 𝒙 𝒚 f_{\mathrm{sim}}:({\bm{\psi}},{\bm{x}})\rightarrow{\bm{y}} italic_f start_POSTSUBSCRIPT roman_sim end_POSTSUBSCRIPT : ( bold_italic_ψ , bold_italic_x ) → bold_italic_y from simulation parameters 𝝍 𝝍 {\bm{\psi}} bold_italic_ψ and input data 𝒙 𝒙 {\bm{x}} bold_italic_x to observations 𝒚 𝒚 {\bm{y}} bold_italic_y , and the goal of the optimization problem is to find parameters 𝝍 𝝍 {\bm{\psi}} bold_italic_ψ that minimize a desired loss function. Sophisticated optimization algorithms typically require gradient information regarding the forward process, f sim subscript 𝑓 sim f_{\mathrm{sim}} italic_f start_POSTSUBSCRIPT roman_sim end_POSTSUBSCRIPT , with respect to the parameters 𝝍 𝝍 {\bm{\psi}} bold_italic_ψ . However, obtaining gradients from black-box simulators can often be prohibitively expensive or, in some cases, impossible.
Furthermore, in many applications, practitioners aim to solve a set of related problems. Thus, starting the optimization ab initio , i.e. from scratch, each time might be inefficient if the forward model is expensive to evaluate.
To address those challenges, this paper introduces a novel method for solving classes of similar black-box optimization problems by learning an active learning policy that guides a differentiable surrogate’s training and uses the surrogate’s gradients to optimize the simulation parameters with gradient descent. After training the policy, downstream optimization of problems involving black-box simulators requires up to ∼ similar-to \sim ∼ 90% fewer expensive simulator calls compared to baselines such as local surrogate-based approaches, numerical optimization, and Bayesian methods.

## 1 Introduction

Across many science and engineering fields, such as medical imaging [, , ], wireless applications [, , , , ], particle physics [, , , ], and molecular dynamics [] and design [], practitioners often face the challenge of inferring unknown properties of a system from observed data [].
In most cases, these settings involve a forward physics-based black-box simulator(^1^11In our study, we consider stochastic and non-stochastic simulators. Our method applies to both types of simulators without requiring any modifications.) $f_{\mathrm{sim}}:(\bm{\psi},\bm{x})\rightarrow\bm{y}$ which maps simulation parameters $\bm{\psi}$ and input data $\bm{x}$ to observations $\bm{y}$ [].
For instance, in particle physics, optimizing the detector geometry to reduce the number of detected events from certain types of particles is crucial to designing an experiment. To that aim, $f_{\mathrm{sim}}$ simulates the detection of particles $\bm{y}$ given their properties, $\bm{x}$, and detector settings, $\bm{\psi}$.
Similarly, in wireless communication, placing a transmitter such that it provides best coverage to a set of users is a recurring problem, for large-scale network planning as well as wifi router placement.
To tackle this problem, $f_{\mathrm{sim}}$ models the received signal $\bm{y}$ at specific locations, given environment conditions $\bm{x}$ (e.g., receiver positions, material properties of the surroundings) and the transmitter position $\bm{\psi}$. All these examples require solving an optimization problem involving a black-box simulator to find the experiment parameters $\bm{\psi}$ that provide the desired observations.

Although these simulators faithfully model the known physical behavior of a system, they are often computationally expensive to run. Therefore, solving
black-box optimization problems with a minimal number of simulator calls is desirable.
When objective functions are (approximately) differentiable, we can use their gradients to guide the optimization process. For appropriate loss landscapes, notably those that are convex, this can achieve strong optimization performance [, , ]. However, many applications involve non-differentiable or, from the perspective of the practitioner, complete black-box simulators [, , , ]. In such cases, gradient-free optimization can be employed, for instance, by using evolutionary strategies [, ], or Bayesian optimization [, ]. To utilize gradient-based optimization, one relies on numerical differentiation [, ] or stochastic gradient estimation methods [, , , , ] to obtain approximate gradients.

Figure: Figure 1: Schematic view of our approach. (a) We study black-box optimization problem (over parameters $\bm{\psi}$), with an emphasis on using gradient information from a fast differentiable surrogate $f_{\phi}$ (b) To optimize $\bm{\psi}$ sample-efficiently, we employ a policy $\pi_{\theta}$ to actively determine whether retraining the surrogate is necessary before using the gradient information.
Refer to caption: extracted/5649542/figures/fig_approach.png

In this paper, we focus on stochastic gradient estimation techniques for black-box optimization. Inspired by Shirobokov et al. [], our approach involves leveraging gradients from a surrogate model trained to (locally) mimic the black-box simulator(^2^22The simulator can be either stochastic or deterministic.).
Gradient-based methods typically perform multiple simulator calls to estimate the gradients, thus making these approaches computationally demanding. To mitigate such a demand, we aim to minimize the number of required simulator calls by proposing to learn a *policy* to guide the optimization. The policy determines whether the current surrogate model (fast) can be used or instead a simulator call is necessary to update the surrogate (slow), see Figure 1. Furthermore, by drawing inspiration from the literature on active learning [, , , , , , , ] we also let our policy learn how to sample new data for training the local surrogate model. This offers additional control, which the policy may learn to exploit.

Our contribution can be summarized as follows: i) We introduce a Reinforcement Learning (RL) framework to learn a policy to reduce the number of computationally expensive calls to a black-box simulator required to solve an optimization problem; ii) We propose to learn a policy that determines when a simulator call is necessary to update the surrogate and when the current surrogate model can be used instead; iii) We implement a policy that also learns how to sample new data for training the surrogate model during the optimization process; iv) We assess the benefits of our RL-based approach on low- and high-dimensional global optimization benchmark functions and two real-world black-box simulators and show that our policy reduces the number of simulator calls up to $\sim$90% , compared to the baselines.

## 2 Related work

#### Simulation-based Inference

Our work lies at the intersection of black-box simulator-based optimization and active learning. Black-box optimization problems are ubiquitous in science and engineering, encompassing scenarios where unknown parameters must be deduced from observational data. These parameters can entail anatomies in MRI [], molecular structures [], particle properties [, ], and cosmological model parameters [], among others. The forward process is often a complex physical process that can be modelled by a simulator but does not provide a likelihood for easy inference. Simulation-based inference techniques aim to infer posterior distributions over these simulation parameters in such likelihood-free settings [, , ]. Other solutions may involve supervised learning on observation-parameter pairs or imitation learning [, ].
Our simulator-based optimization setting is a variation on these problems. Here, the objective is to find the optimal parameters of the simulator, where optimality is typically formulated in terms of desired observations. This methodology can be applied in various fields, such as MRI [, , ], particle physics [, , , ] and molecular design []. When the simulators are differentiable, direct gradient-based optimization can perform well [, , ]. However, a different approach is necessary in cases where the simulators are non-differentiable. Well-known gradient-free methods that may be employed in such settings include evolutionary strategies [, ] and Bayesian optimization [, , , ]. Nevertheless, these methods often require additional assumptions to make the optimization scalable in high dimensional parameter spaces [, ].

#### Approximate-Gradient Optimization

With the rise of deep learning, there has been a surge of interest in approximate-gradient optimization methods. While some authors consider numerical differentiation [, ], many others have focused on methods for efficiently obtaining approximate stochastic gradients [, , , , , ]. Another strategy involves training differentiable surrogate models to mimic the simulator and assuming that the gradients of the surrogate model are similar enough to those of the simulator []. Surrogate models have been trained for many applications, including wireless propagation modeling [, ], space weather prediction [], material discovery [], and fluid dynamics simulation []. This trend provides an opportunity for surrogate-based optimization of simulators, as surrogate models are readily available. Additionally, it has been observed by Shirobokov et al. [] that using (local) surrogate gradients is more efficient than many alternatives. Our work generalizes this setup by introducing a policy that guides the optimization by suggesting when and, optionally, how the surrogate should be updated during the optimization process.

#### Active Learning

When the policy decides how the surrogate should be updated, it does so using information provided by the surrogate. This is an example of active learning [], where the current instance of a task model (the surrogate) affects the data it sees in future training iterations. In particular, our policies are instances of learning active learning, where a separate model (our policy) is trained to suggest the data that the task model should be trained on [, , , , , , , ].

## 3 Background

We aim to optimize the simulation parameters of a black-box simulator using stochastic gradient descent.
The black-box simulator, $f_{\mathrm{sim}}$, describes a stochastic process(^3^33A non-stochastic simulator can be considered as a special case where $f_{\mathrm{sim}}$ places a delta distribution over observations.), $p(\bm{y}|\bm{\psi},\bm{x})$, from which we obtain the observations as $\bm{y}\ =\ f_{\mathrm{sim}}({\bm{\psi}},{\bm{x}})\sim p(\bm{y}|\bm{\psi},\bm{x})$, where $\bm{x}\sim q(\bm{x})$ is a stochastic input and ${\bm{\psi}}$ is the vector of simulation parameters.
Since, these simulators are typically not differentiable, we train a surrogate neural network to locally (in $\bm{\psi}$) approximate the simulator []. Gradients of these local surrogates, obtained through automatic differentiation, might then be used to perform the optimization over $\bm{\psi}$.
The goal is now to minimize an expected loss $\mathcal{L}$ over the space of the simulation parameters $\bm{\psi}$. As the functional form of the simulator is generally unknown, this expectation cannot be evaluated exactly and is instead estimated using $N$ Monte Carlo samples:

$$ $\displaystyle\bm{\psi}^{*}=\operatorname*{arg\,min}_{\bm{\psi}}\mathbb{E}\left [\mathcal{L}(\bm{y})\right]$ $\displaystyle=\operatorname*{arg\,min}_{\bm{\psi}}\int\mathcal{L}(\bm{y})\,p( \bm{y}|\bm{\psi},\bm{x})\,q(\bm{x})\,d\bm{x}\,d\bm{y},$ (1) $\displaystyle\approx\operatorname*{arg\,min}_{\bm{\psi}}\,\frac{1}{N}\,\sum_{i =1}^{N}\mathcal{L}(f_{\mathrm{sim}}(\bm{\psi},\bm{x}_{i})).$ (2) $$

After training a neural network surrogate $f_{\phi}:(\bm{\psi},\bm{x},\bm{z})\rightarrow\bm{y}$ on data generated with $f_{\mathrm{sim}}$, the optimization might be performed following gradients of the surrogate. Here, $\bm{z}$ is a randomly sampled latent variable that accounts for the stochasticity of the simulator. Gradients are then estimated as:

$$ $\nabla_{\bm{\psi}}\mathbb{E}\left[\mathcal{L}(\bm{y})\right]\approx\frac{1}{N} \,\sum_{i=1}^{N}\,\nabla_{\bm{\psi}}\,\mathcal{L}\left(f_{\phi}(\bm{\psi},\bm{ x}_{i},\bm{z}_{i})\right).$ (3) $$

Since running the forward process $f_{\mathrm{sim}}$, is often an expensive procedure, our goal is to minimize the number of simulator calls required to solve the optimization problem at hand.

## 4 Policy-based Black-Box Optimization

Following Shirobokov et al. [], we perform an iterative optimization based on the gradients obtained in  Eq. 3. At each point during the optimization, new values $\bm{\psi}_{j}$ are sampled within a box of fixed size $\epsilon$, centered around the current $\bm{\psi}$: $U_{\epsilon}^{\bm{\psi}}=\{\bm{\psi}_{j};|\bm{\psi}_{j}-\bm{\psi}|\leq\epsilon\}$. Then, input samples are obtained from $q({\bm{x}})$, and the simulator is called to obtain the corresponding $\bm{y}$ values. The resulting samples are stored in a history buffer $H$, from which the surrogate is trained from scratch.
Specifically, the surrogate is trained on samples extracted from $H$ that satisfies the condition that $\bm{\psi}_{j}$ lies within $U_{\epsilon}^{\bm{\psi}}$. The overall process required to generate new samples from the black-box simulator is what we refer to as a “*simulator call*”.

#### Policy-based Approach

We propose further reducing the number of simulator calls required for an optimization run with an RL-based approach. Our method involves utilizing a learned policy $\pi_{\theta}$, with learnable parameters $\theta$ to: i) decide whether a simulator call should be performed to retrain the local surrogate; and ii) define how to sample from the black-box simulator.

#### Sampling Strategy

To investigate the question concerning how to perform a simulator call, we train policies to additionally output the $\epsilon$ for constructing the sampling neighbour $U_{\epsilon}^{\bm{\psi}}$, which serves as our data acquisition function.
As $\epsilon$ parameterizes this acquisition function, such policies are an example of active learning []. In particular, these policies are instances of learning active learning [, , , , , , , ], as they learn a distribution over $\epsilon$. See Appendices B.1 and C.1 for a more detailed description concerning the policy implementation and training.

#### State Definition

We formalize the sequential optimization process as an episodic Markov Decision Process (MDP). The state $s_{t}$ (at timestep $t$) is given by the tuple $(\bm{\psi}_{t},t,l_{t},\sigma_{t})$, where $\bm{\psi}_{t}$ is the current parameter value, $l_{t}$ is the number of simulator calls already performed in the episode, and $\sigma_{t}$ is some measure of uncertainty produced by the surrogate.

#### Action Definition

Actions $a_{t}$ consist of binary valued variables $b\in\{0,1\}$, sampled from a Bernoulli distribution, where $1$ represents the decision to perform a simulator call. Additionally, we also train policies to determine, as part of the action, the trust region size for sampling new values for $\bm{\psi}$. The dynamics of the MDP is represented by means of the Adam optimizer [] which updates the current state by performing a single optimization step in the direction of the gradients of ${\bm{\psi}}$.

#### Reward Design

Episodes come to an end under three conditions: $\mathcal{A)}$ when the optimization reaches a parameter for which $\mathbb{E}\left[\mathcal{L}(\bm{y})\right]$ is below the target value $\tau$ (we call this termination - see Appendix D.4 for details concerning the choice of $\tau$); $\mathcal{B)}$ when the maximum number of timesteps $T$ is reached; or $\mathcal{C)}$ when the policy hits the available budget for simulator calls $L$. To incentivize reducing the number of simulator calls, rewards $r(s_{t},a_{t},s_{t+1})$ are 0 0 if $b=0$ and $-1$ if $b=1$. Additionally, a reward penalty is added when $\mathcal{B)}$ or $\mathcal{C)}$ occur to promote termination. The penalty is $-(L-l_{t})-1$ when $\mathcal{B)}$ occurs and $-1$ when $\mathcal{C)}$ occurs. This ensures the sum-of-rewards for non-terminating episodes is $-L-1$. We have observed that using reward penalties based on $l_{t}$ rather than $t$ improves training stability. We refer the reader to Appendix D.3 for further details concerning the reward design.

#### Local Surrogate

The decision to perform a simulator call should rely on the quality of the local surrogate. A well-fitted surrogate to the simulator at the current $\bm{\psi}$ will presumably provide useful gradients, so gathering additional data and retraining is unnecessary. Vice-versa, a badly fitted surrogate will likely not provide useful gradients and may be worth retraining, even if a simulator call is expensive. We use the uncertainty feature $\sigma$ to provide this information.

#### Local Surrogate Ensemble

To construct $\sigma$, we replace the local surrogate with an ensemble of local surrogates, all trained on and applied to the same input data. The use of an ensemble empowers our approach with the ability to estimate uncertainties without the overhead of training a posterior network as in the case of Bayesian models [, ]. Each surrogate is implemented as a two-layer Multi-Layer Perceptron (MLP) with Rectified Linear Unit (ReLU) activation function. Therefore, the additional resource requirement for training an ensemble instead of a single surrogate is negligible.
As the input, we use the tuple $({\bm{\psi}},{\bm{x}},{\bm{z}})$, where ${\bm{z}}$ is sampled from diagonal Normal distribution.

#### Uncertainty Feature

We compute the prediction mean per surrogate on $D$ samples as $\bm{\bar{y}}=\frac{1}{D}\,\sum_{i=1}^{D}\left[f_{\phi}(\bm{\psi},\bm{x}_{i},
\bm{z}_{i})\right]$, and construct $\sigma$ as the standard deviation over these mean predictions.
Specifically, ${\bm{z}}$ accounts for the stochasticity of $f_{\mathrm{sim}}$. Such an idea allows us to dramatically simplify the surrogate architecture compared to Shirobokov et al. []. Training GANs [] is notoriously more challenging than training a shallow MLP due to instabilities and mode collapse. Nonetheless, our “simpler” surrogate has enough capacity to locally approximate highly complex stochastic, and non-stochastic, simulators. See Appendix B.2 for further details concerning models implementation.

## 5 Experimental Results

Figure: Figure 2: Loss landscape and learned optimization trajectory for the Probabilistic Three Hump problem. The yellow region denotes $\bm{\psi}$ values that lead to termination. The $\epsilon=0.5$ neighbour around $\bm{\psi}_{0}$ (black cross) is visualized as the red box. Light green and blue arrows represent gradients from the surrogate or after a simulator call, respectively.
Refer to caption: extracted/5649542/figures/probabilistic_three_hump_opt_traj_1.png

To assess the performance of our method, we test it on two different types of experiments.

First, we consider stochastic versions of benchmark functions available in the optimization literature []. We consider the Probabilistic Three Hump, the Rosenbrock, and the Nonlinear Submanifold Hump problems.
These benchmark functions are relevant for two reasons:
(i) they allow us to compare our models against baselines on similar settings as in []; and
(ii) they allow to easily gain insights into models performance. Especially, since the Probabilistic Three Hump problem is a two dimensional problem, i.e. ${\bm{\psi}}\in\mathbb{R}^{2}$, we are able to conveniently visualize the objective landscape as well as the optimization trajectories. Furthermore, the Rosenbrock and Nonlinear Submanifold Hump problems allow us to test our approach on high-dimensional, more complex, problems before moving to real-world black-box simulators.

The second type of experiments concerns real-world black-box simulators. We consider applications from two different scientific fields, namely the Indoor Antenna Placement problem for wireless communications and the Muon Background Reduction problem for high energy physics.

#### Baselines

We compare our method against three baselines. We consider Bayesian optimization using Gaussian processes with cylindrical kernels [], numerical differentiation with gradient descent, and local surrogate-based methods (L-GSO) []. Furthermore, to guarantee a fair comparison against our models, we formulated an ensemble version(^4^44L-GSO-E averages the gradients over the ensemble. The model does not leverage any uncertainty since it always calls the black-box simulator.) (L-GSO-E) of the local surrogate for L-GSO.

#### Our Models

Policy methods are split into those that only output when to perform a simulator call, $\pi$-E, and those that also output how to sample ${\bm{\psi}}$ values by providing the neighbour size, $\epsilon$, that parameterizes the acquisition function, $\pi_{AL}$-E. L-GSO, its ensemble version (L-GSO-E), and $\pi$-E use a fixed value for $\epsilon$ that depends on the problem at hand (see Appendix D.1).
Finally, $\pi_{AL}^{G}$-E is a version of $\pi_{AL}$-E where the surrogate ensemble is always warm-started from the previous training step, such that the surrogate is continuously improved along the observed trajectories through $\bm{\psi}$-space (see Appendix B.1 for more details).

Figure: Figure 3: Benchmark function results. Top row: Probabilistic Three Hump problem. Middle row: Rosenbrock problem. Bottom row: Nonlinear Submanifold Hump Problem. AMO (the lower the better) concerning (a) a fixed and (c) a parameterized ${\bm{x}}$ distribution. ANC (the lower the better) regarding (b) a fixed and (d) a parameterized ${\bm{x}}$ distribution. Uncertainties are quantified over evaluation episodes and different random seeds.
Refer to caption: x1.png

#### Metrics

We report experimental results by using two different metrics: the Average Minimum of the Objective function (AMO) for a specific budget of simulator calls and the Average Number of simulator Calls (ANC) required to terminate an episode. The first quantity answers the question: *What is the lowest value for the objective function achievable for a given budget of simulator calls?*; that might be used as an indicator of the efficacy of each simulator call. The second quantity answers the question: *What is the simulator call budget required, on average, to solve a black-box optimization problem?*; that might indicate how good the policy is at leveraging the surrogate and understanding its reliability. Therefore, those two metrics allow us to benchmark our approach against others by looking at relevant quantities (see Appendix D.5 for more details).

### 5.1 Benchmark Functions

We consider a fixed and a parameterized input distribution for each benchmark function. Specifically, the latter setup corresponds to solving an entire family of related optimization problems each characterized by a different input distribution, $q_{i}({\bm{x}})$.
During training and evaluation of the policy, each episode is characterized by a different input distribution. In what follows we report the definition for each benchmark function only. A more detailed description can be found in Appendix D.1.

#### Probabilistic Three Hump Problem

As mentioned in the introduction to the section, the Probabilistic Three Hump problem concerns the optimization of a 2-dimensional vector ${\bm{\psi}}$. Specifically, the goal is to find ${\bm{\psi}}^{*}$ such that:
${\bm{\psi}}^{*}=\mathrm{arg\ min}_{\bm{\psi}}\mathbb{E}[\mathcal{L}(y)]=
\mathrm{arg\ min}_{\bm{\psi}}\ \mathbb{E}[\sigma(y-10)-\sigma(y)]$, where $\sigma$ is the sigmoid function and ${\bm{y}}$, the observations vector, is given by: $y\sim\mathcal{N}(y;\ \mu_{i},1),\,i\in\{1,2\}$.

Being ${\bm{\psi}}$ a 2-dimensional vector, its optimization trajectory is amenable to visualization. Figure 2 illustrates that a fully trained policy can exploit the local-surrogate as much as possible and only perform a simulator call when the model is far from the initial training location (red square in Figure 2).
Intuitively, such behaviour is foreseen. The surrogate model is expected to provide meaningful gradients in proximity to the ${\bm{\psi}}$ region where it was previously trained.
However, as we move away from that region, we expect the quality of the gradients to decline until a simulator call is triggered and the local-surrogate re-trained. However, moving away from the last training region is not the sole condition that might trigger a simulator call. For instance, towards the end of the trajectory, the policy decides to call the simulator twice to gather more data to train the surrogate and then calls the simulator again before ending the episode, indicating that a rapidly changing loss landscape may also trigger a simulator call.

#### Rosenbrock Problem

In the Rosenbrock problem, we aim to optimize ${\bm{\psi}}\in\mathbb{R}^{10}$ such that: ${\bm{\psi}}^{*}=\mathrm{arg\ min}_{\bm{\psi}}\ \mathbb{E}[\mathcal{L}(y)]=
\mathrm{arg\ min}_{\bm{\psi}}\ \mathbb{E}[y]$; where $y$ is given by: $y\sim\mathcal{N}(y;\ \gamma+x,\ 1)$, where $\gamma=\sum^{n-1}_{i=1}\left[(\psi_{i}-\psi_{i+1})^{2}+(1-\psi_{i})^{2}\right]$.

#### Nonlinear Submanifold Hump Problem

This problem share a similar formulation to the Probabilistic Three Hump problem. However, the optimization is realized by considering the embedding $\hat{{\bm{\psi}}}={\bm{B}}\ \mathrm{tanh}({\bm{A}}{\bm{\psi}})$, where ${\bm{A}}\in\mathbb{R}^{16\times 40}$ and ${\bm{B}}\in\mathbb{R}^{2\times 16}$, of the vector ${\bm{\psi}}\in\mathbb{R}^{40}$. Subsequently, $\hat{{\bm{\psi}}}$ is used in place of ${\bm{\psi}}$ in the Probabilistic Three Hump problem definition.

### 5.2 Real-world Simulators

We now focus on real-world optimization problems involving computationally expensive, non-differentiable black-box simulators. First, we look at the field of wireless communications considering two settings with a (non-stochastic) wireless ray tracer []. Then, we move to the world of subatomic particles and solve a detector optimization problem for which we use the high energy physics toolkits Geant4 [] (stochastic simulator) and FairRoot [].

#### Wireless Communication: Indoor Transmitting Antenna Placement

We study the problem of optimally placing a transmitting antenna in indoor environments to maximize the signal strength at multiple receiver locations.
Determining the signal strength in such a scenario typically requires a wireless ray tracer ([] in our case), which takes as input the transmit location candidate ${\bm{\psi}}\in{\mathbb{R}}^{3}$, alongside other parameters (e.g., receive locations, 3D mesh of scene).
To predict the signal strength for a particular link (i.e., a transmit-receive antenna pair), the ray tracer exhaustively identifies multiple propagation paths between the two antennas and calculates various attributes of each path (e.g., complex gains, time-of-flight).
The signal strength is computed from the coherent sum of the complex-valued gains of each path impinging on the receive antenna and is represented in log-scale (specifically, dBm).
Optimally placing the transmit antenna is typically slow, as it amounts to naively and slowly sweeping over transmit location choices ${\bm{\psi}}$ and observing the simulated signal strengths.
Instead, we employ our approach to “backpropagate” through the the surrogate and perform gradient descent steps on the location ${\bm{\psi}}$.
Specifically, we consider two indoor scenes for this experiment and investigate how to use our approach to find an optimal transmit location that maximizes signal strength in the 3d scene (column (a) in Figure 4).
The end goal in both cases is to find an optimal transmit antenna location ${\bm{\psi}}$ that maximizes the median signal strength calculated over a distribution of receive locations ${\bm{x}}\sim q({\bm{x}})$ (see Appendix D.2 for more details concerning simulations).

Figure: Figure 4: Wireless ray-tracing results. (a) Rendering of the indoor environment, (b) AMO (the lower the better), and (c) ANC (the lower the better). Uncertainties are quantified over evaluation episodes and different random seeds. Top row: Conference room environment. Bottom row: Office room environment.
Refer to caption: x2.png

#### Physics: Muon Background Reduction

We consider the optimization of the active muon shield for the SHiP experiment [].
Typically, optimizing a detector is a crucial step in designing an experiment for particle physics. For instance, the geometrical shape, the intensity and orientation of magnetic fields, and the materials used to build the detector play a crucial role in defining the detector’s “sensitivity” to specific types of particle interactions, i.e. events.
Observed events are usually divided into signal, i.e., interactions physicists are interested in studying, and “background”, i.e., events that are not of any interest and that might reduce the detector’s sensitivity. Concerning the SHiP experiment, muons represent a significant source of background; therefore, it is necessary to shield the detector against those particles.
The shield comprises six magnets, left image in Figure 7, each described by seven parameters. Hence, ${\bm{\psi}}\in\mathbb{R}^{42}$. To run the simulations, we use the Geant4 [] and FairRoot [] toolkits. The input distribution ${\bm{x}}$ describes the properties of incoming muons(^5^55Concerning the muons distribution, we use the same dataset as in Shirobokov et al. []. The dataset is available for research purposes.). Specifically, as in [], we consider the momentum ($P$), the azimuthal ($\phi$) and polar ($\theta$) angles with respect to the incoming $z$-axis, the charge $Q$, and $(x,y,z)$ coordinates. The goal is to minimize the expected value of the following objective function:

$$ $\displaystyle\mathcal{L}({\bm{y}};{\bm{\alpha}})=\sum^{N}_{i=1}\mathbb{I}_{Q_{ i}=1}\sqrt{(\alpha_{1}-({\bm{y}}_{i}+\alpha_{2}))/\alpha_{1}}+\mathbb{I}_{Q_{i }=-1}\sqrt{(\alpha_{1}+({\bm{y}}_{i}-\alpha_{2}))/\alpha_{1}}$ (4) $$

where $\mathbb{I}$ is the indicator function, $\alpha_{1}$ and $\alpha_{2}$ are known parameters defining the sensitive region of the detector, and $Q$ and ${\bm{y}}$ represent the electric charge and the coordinates of the observed muons, respectively. Minimizing $\mathcal{L}({\bm{y}};{\bm{\alpha}})$ corresponds to minimize the number of muons hitting the sensitive region of the detector.

### 5.3 Results & Discussion

Concerning the problems involving benchmark functions, policy-based methods achieve the highest performance in both scenarios with fixed and parameterized ${\bm{x}}$-distributions, as shown in Figure 3, with the $\pi$-E model scoring best especially concerning the average number of simulator calls required to terminate an episode (bar plots in the figures). Notably, we observe a significant reduction in the number of simulator calls, up to $\sim 90\%$. While the $\pi_{AL}$-E and $\pi_{AL}^{G}$-E models outperform all the baselines as well, there is no clear advantage compared to $\pi$-E. Those results can be interpreted based on the argument that, for the benchmark functions, we set the trust region size to the optimal value reported in Shirobokov et al. [], and therefore simplifying the problem for the $\pi$-E model (the L-GSO baselines use by default the “optimal” trust region size). However, it should also be noted that using a warm-started surrogate, $\pi_{AL}^{G}$-E, helps mitigate such a difficulty by improving the results compared to $\pi_{AL}$-E.
Similarly, the average minimum objective function value shows that our policies are typically better, with the only exception of the Nonlinear Submanifold Hump problem for which our models agree, within the error, with the baselines.
As noted in [], the BOCK baseline struggles in solving the Rosenbrock problem (Figure 3, middle row), likely due to the high curvature of the objective function under analysis. On the other hand, numerical differentiation appears to be less affected by this issue, thus reporting acceptable results for all the problems involving benchmark functions.
However, given that local surrogate methods (L-GSO) have been shown to outperform the other baselines or be on par in the worst case, we only consider that in the experiments involving real-world black-box simulators.

Figure: Figure 5: Physics experiments results. (a) Schematic view (^7^77Image from “*Optimising the active muon shield for the SHiP experiment at CERN*" In Journal of Physics: Conference Series, vol. 934, no. 1, p. 012050. IOP Publishing, 2017, by Baranov, A., et al. Licensed under CC BY 3.0) [] of the active muon shield baseline configuration. The detail concerning the “Target/Magnet hadron absorber” is not relevant to the current discussion. However, we reported that only for completeness. See [] for more details. (b) AMO (the lower the better), (c) and ANC (the lower the better). Uncertainties are quantified over evaluation episodes and different random seeds.
Refer to caption: x3.png

Regarding the real-world black-box simulators, we draw similar conclusions to the benchmark functions. Policy-based methods achieve the best performance in terms of AMO and ANC. However, in contrast to the problems involving the benchmark functions, the performance of the three different policies typically agrees within the uncertainties (Figure 4) with, in some cases, $\pi_{AL}$-E and $\pi_{AL}^{G}$-E performing better, on average, than $\pi$-E (Figure 7). It may be possible, therefore, that the true advantage of learning to adapt the trust region size, $\pi_{AL}$-E, combined with a warm-started surrogate, $\pi_{AL}^{G}$-E, is only revealed in highly complex optimization problems, such as the optimization of a detector for high energy physics experiments (Figure 7). We leave it to future works to conduct further investigations to assess the actual advantages that will result from learning the sampling strategy. We refer to Appendix A for a discussion concerning limitations and future works.

## 6 Conclusions

We proposed a novel method to minimize the number of simulator calls required to solve optimization problems involving black-box simulators using (local) surrogates. The core idea of our approach is to reinforce an active learning policy that controls when the black-box simulator is used and how to sample data to train the local surrogates. We trained three policy models and compared them against baselines, including local surrogate methods [], numerical differentiation, and Bayesian approaches []. We tested our approach on two different experimental setups: benchmark functions and real-world black-box simulators.

Our policy-based approach showed the best performance across all scenarios. In particular, compared to the baselines, we observed a significant reduction in the number of simulator calls, up to $\sim 90\%$. Our results suggest that local surrogate-based optimization of problems involving black-box forward processes benefits from the guidance of both simple policies and learned sampling strategies.

## Appendix A Broader Impact, Limitations and Future Workds

#### Broader Impact

This paper proposes a novel policy-based approach to guide local surrogate-based problem optimization with black-box simulators.
We believe the potential societal consequences of our work are chiefly positive, as it has the potential to promote the use of policy-based approaches in various scientific domains, particularly concerning optimization procedures involving black-box, non-differentiable, forward processes.
However, it is crucial to exercise caution and thoroughly comprehend the behaviour of the models to obtain tangible benefits.

#### Limitations and Future Works

Gradient-based optimization may get stuck in local optima of the loss surface $\mathbb{E}_{p(\bm{y}|\bm{\psi},\bm{x})}\left[\mathcal{L}(\bm{y})\right]$. Investigating whether introducing a policy into the optimization can help avoid such local minima, is an interesting direction of future research. The Probabilistic Three Hump problem has no local minima but does contain a few flat regions, where gradient-based optimization is more challenging. Exploratory experiments have provided weak evidence that the policy may learn to avoid such regions.

Hyperparameter tuning has mostly involved reducing training variance through tuning the number of episodes used for a PPO iteration, as well as setting learning rates and the KL-threshold. Little effort has been spent optimizing the policy or surrogate architectures; we expect doing so to further improve performance. Similarly, while PPO with a value function critic is a widely used algorithm, more recent algorithms may offer additional advantages, such as improved planning and off-policy learning for more data-efficient training [, ].

## Appendix B Implementation details

### B.1 Policy

Figure: (a) Actor.
Refer to caption: x4.png

The policy $\pi_{\theta}$ is composed of two separate neural networks: an Actor and a Critic. Both networks are ReLU MLPs with a single hidden layer of 256 neurons, schematically depicted in Figure 6. The input to both networks is the tuple: $(\bm{\psi}_{t},t,l_{t},\sigma_{t})$, where $\bm{\psi}_{t}$ is the current parameter value (at timestep $t$), $l_{t}$ is the number of simulator calls already performed this episode, and $\sigma_{t}$ is the standard deviation over the average surrogate predictions in the ensemble.

The Actor outputs either one or three values. The first value is passed through a sigmoid activation and treated as a Bernoulli random variable, from which we sample $b$, representing the decision to perform a simulator call or not. If the policy outputs three values, the second and third values are treated as the mean and standard deviation of a lognormal distribution from which we sample $\epsilon$, the trust region size, for the current timestep. The standard deviation value is passed through a softplus activation function to ensure it is positive.

The Critic outputs a value-function estimate $V_{\theta}(s)$, where $\theta$ are policy parameters. We use this estimate to compute advantage estimates in PPO, as explained in detail in section C.1. Since rewards have unity order of magnitude, we expect return values to be anywhere in $[-T,0]$. To prevent scaling issues, we multiply the Critic output values by $T$ before using them for advantage estimation.

#### The π A ⁢ L G superscript subscript 𝜋 𝐴 𝐿 𝐺 \pi_{AL}^{G} italic_π start_POSTSUBSCRIPT italic_A italic_L end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_G end_POSTSUPERSCRIPT method

When training a policy for downstream optimization of many related black-box optimization problems, it may be helpful to train a global surrogate simultaneously for such a problem setting. Such a global surrogate might provide better gradients for problem optimization, especially if it has been jointly optimized with the policy. We have implemented the $\pi_{AL}^{G}$ method to test this. Here, the policy outputs both the decision to perform a simulator call and the trust-region size, $\epsilon$, just as in $\pi_{AL}$. However, the surrogate ensemble is “warm-started” from the previous training step every time a retraining decision is made. This results in a continuously optimized surrogate ensemble for the training trajectories. To prevent the surrogate from forgetting old experiences too quickly, we employ a replay buffer that undersamples data from earlier iterations geometrically. Specifically, when training the surrogate with trust-region $U_{\epsilon}^{\bm{\psi}}$, we include all data inside $U_{\epsilon}^{\bm{\psi}}$ for the current episode, half of the data inside $U_{\epsilon}^{\bm{\psi}}$ from the previous episode, a quarter of the data seen two episodes ago, and so on.

### B.2 Surrogate

Figure: Figure 7: Schematic for the surrogate architecture. The surrogate is an MLP trained to mimic the simulator. It takes $(\bm{\psi},\bm{x},\bm{z})$ as input and outputs $\bm{y}$.
Refer to caption: x6.png

The surrogate consists of a ReLU MLP with two hidden layers of 256 neurons that takes as input $(\bm{\psi},\bm{x},\bm{z})$ and outputs $\bm{y}$. $\bm{z}$ is sampled from a 100-dimensional diagonal unit Normal distribution. The surrogate architecture is schematically depicted in Figure 7.

Surrogates are trained on data generated from $f_{\textrm{sim}}$. Following the approach outlined in Shirobokov et al. [], we sample $M$ values $\bm{\psi}_{j}$ inside the box $U_{\epsilon}^{\bm{\psi}}$ around the current parameter value using an adapted Latin Hypercube sampling algorithm. For each of those $\bm{\psi}_{j}$, we then sample $N=3\cdot 10^{3}$ $x$-values. We use $M=5$ for the Probabilistic Three Hump problem, $M=16$ for the Rosenbrock problem, and $M=40$ for the Nonlinear Submanifold Hump problem. As in Shirobokov et al. [], this means a single “simulator call” consists of $1.5\cdot 10^{4}$ function evaluations for Probabilistic Three Hump, $4.8\cdot 10^{4}$ for Rosenbrock, and $6.0\cdot 10^{4}$ for the Nonlinear Submanifold Hump.

To train the surrogates, we use the Adam optimizer for two epochs with a learning rate of $10^{-3}$ and a batch size of $512$. Each ensemble comprises three surrogates, each trained on identical data but a different random seed.

## Appendix C Training Details

### C.1 Training

We train our policy in an episodic manner by accumulating sequential optimization episodes and updating the policy using Proximal Policy Optimization (PPO) [] with Generalised Advantage Estimation (GAE) advantages [] (discount factor $\gamma=1.0$, GAE $\lambda=0.95$). Episodes terminate once any of the following conditions is met: A) the target value for the loss, $\tau$, has been reached, B) the number of timesteps $T=1000$ has been reached, or C) the number of simulation calls $L=50$ has been reached. For every training iteration, before doing PPO updates, we accumulate: 10 episodes for the Nonlinear Sub. Hump and 16 episodes for the Rosenbrock and Prob. Three Hump problems, 10 episodes concerning the wireless simulations and 5 for the high energy physics experiments. The different choices in the number of episodes to accumulate are mainly dictated by the time required to complete one episode.

We use the PPO-clip objective (with clip value $0.2$) on full trajectories with no entropy regularization to perform Actor updates.
We perform multiple Actor updates with the same experience until either the empirical KL-divergence between the old and new policy reaches a threshold ($3\cdot 10^{-3}$ for simulator-call decision actions, $10^{-2}$ for trust-region size $\epsilon$ actions), or 20 updates have been performed. In practice, we rarely perform the full 20 updates. Updates use the Adam optimizer with learning rate $3\cdot 10^{-4}$.

Similarly, we perform multiple Critic updates using the Mean-Squared Error (MSE) between the estimated values $V_{\theta}(s_{t})$ and the observed return (sum of rewards, as $\gamma=1.0$) $R_{t}$ at every timestep. We keep updating until either $\text{MSE}\leq 30.0$ or ten updates have been done. This approach helps the critic learn quickly initially and after seeing surprising episodes but prevents it from over-updating on similar experiences (as MSE will be low for those iterations). Updates use Adam with learning rate $10^{-4}$. See Algorithm 1 for the training pseudo-code and Algorithm 2 for the evaluation procedure pseudo-code.

To assess the performance of our models, we run 32 evaluation episodes for the benchmark functions and 20 and 5 evaluation episodes for the wireless and physics experiments, respectively. Moreover, we consider three random seeds for L-GSO and policy models, while we used ten random seeds for the BOCK and Num. Diff. baselines.

Figure: Algorithm 1 Training the (active learning) policy.

Figure: Algorithm 2 Inference with the (active learning) policy.

Figure 8 shows that the policy is actually able to learn when to call the simulator. Initially, during the first stages of the training, the policy generates completely random actions, resulting in an average probability of calling the simulator close to 0.5 (bottom plot in Figure 8). However, as the training progresses, such a probability gradually decreases, leading to a reduction in the number of simulator calls (top plot in Figure 8).

Figure: Figure 8: Top: Average number of simulator calls (purple curve) and average sum of rewards (blue curve) after each PPO iteration. Values are averaged across evaluation episodes. The upper and lower bound for the shadowed areas represent the max. and min. for each of the two mentioned metrics, respectively. Bottom: Average value of the probability of calling the black-box simulator (orange curve) and average critic loss (blue curve). Values are averages across evaluation episodes for each PPO iteration.
Refer to caption: extracted/5649542/figures/probabilistic_three_hump_policy_iterations.png

### C.2 Objective Landscape and Optimization Trajectory

Experiments with low-dimensional functions, such as the Probabilistic Three Hump problem (${\bm{\psi}}\in\mathbb{R}^{2}$), allow us to easily visualize optimization trajectories to gain insights into the models behaviour.

As mentioned in the main corpus of the paper, practitioners in many scientific fields may need to solve a set of related balck-box optimization problems that can become costly if each optimization process has to begin ab initio. Therefore, we investigated the robustness of the policy trained on a given setup, i.e. input $x$-distributions, and then tested on different ones. To mimic such a scenario, we consider a parameterized input ${\bm{x}}$-distribution. In real-world experiments, such a variation could correspond to different properties of the input data used to run the simulations. We already report the results concerning such tests in Section 5. In Figure 9, we show the optimization landscape for different ${\bm{x}}$-distributions for the Prob. Three Hump problem.
It is worth noticing that, although the minima generally correspond to similar neighbours of the $\psi$ values, the landscape dramatically changes from one distribution to another.

Figure: Figure 9: Loss landscape for the Probabilistic Three Hump Problem for different bounds on the $x$-distribution. Shown is $\frac{1}{n}\,\sum_{i=0}^{n-1}\mathcal{L}(f_{s}(\bm{\psi},\bm{x}_{i}))$ for a grid of $\bm{\psi}$ values ($n=100$). The yellow region denotes $\bm{\psi}$ values that lead to termination.
Refer to caption: extracted/5649542/figures/probabilistic_three_hump_psi_surface.png

### C.3 Experiments Compute Resources

Performing a single optimization for the benchmark functions and the wireless experiments does not require a significant amount of computational resources and can be conducted using any commercially available NVIDIA GPU. A single optimization can be easily fitted on a single GPU.
On the other hand, physics experiments require extensive computing resources for running simulations. While it is still feasible to run the entire optimization on a single machine, it might take a consistent amount of time when simulating thousands of particles. The primary bottleneck for such experiments stems from the Geant4 [] simulator, which is highly CPU-demanding.
Since the simulations of individual particles are independent of each other, they can be run in parallel without communication between processes.
In our experiments, we split up each simulation into chunks of 2000 particles which resulted in run times of 5-15 minutes per simulation on single CPU core, depending on the exact hardware.

## Appendix D Experimental Details

### D.1 Benchmark Functions

Our tests with benchmark functions employ a probabilistic version of three benchmark functions from the optimization literature: Probabilistic Three Hump, Rosenbrock, and Nonlinear Submanifold Hump. The first one is a two-dimensional problem that lends itself well to visualization. Instead, the $N$-dimensional Rosenbrock (with $N=10$) and the Nonlinear Submanifold Hump problems are used to test our method on higher-dimensional settings.

#### Probabilistic Three Hump Problem

The goal is to find the 2-dimensional $\bm{\psi}$ that optimizes:(^8^88Here the upper bound of $x_{1}$ and lower bound of $x_{2}$ are switched compared to the notation in Equation (3) of []. These bounds match the official implementation of L-GSO as of August 2023.)

$$ $\displaystyle\bm{\psi}^{*}=\operatorname*{arg\,min}_{\bm{\psi}}\,\mathbb{E} \left[\mathcal{L}(y)\right]=\operatorname*{arg\,min}_{\bm{\psi}}\,\mathbb{E} \left[\sigma(y-10)-\sigma(y)\right],\text{s.t.}$ $\displaystyle y\sim\mathcal{N}(y|\mu_{i},1),\,i\in\{1,2\},\;\mu_{i}\sim \mathcal{N}(x_{i}h(\bm{\psi}),1),\;x_{1}\sim U[-2,2],\;x_{2}\sim U[0,5],$ (5) $\displaystyle P(i=1)=\frac{\psi_{1}}{||\bm{\psi}||_{2}}=1-P(i=2),\;h(\bm{\psi} )=2\psi_{1}^{2}-1.05\psi_{1}^{4}+\psi_{1}^{6}/6+\psi_{1}\psi_{2}+\psi_{2}^{2}.$ $$

We consider an episode terminated when $\mathbb{E}\left[\mathcal{L}(y)\right]=\frac{1}{N}\,\sum_{i=1}^{N}\mathcal{L}(f
_{\textrm{sim}}(\bm{\psi},\bm{x}_{i}))\leq\tau=-0.8$, which we evaluate after every optimization step using $N=10^{4}$ samples. Following [], we use $\epsilon=0.5$ as the trust-region size. The optimization is initialized at $\bm{\psi}_{0}=[2.0,0.0]$; this is a symmetry point in the Three Hump function such that optimization with stochastic gradients can fall into either of the two wells around the two minima of the function. Such a procedure requires our methods to learn good paths to both optima, making the task more interesting. In principle, optimization could be initialized at any $\bm{\psi}_{0}$.

#### Rosenbrock Problem

The goal for this problem is to find the 10-dimensional $\bm{\psi}$ that optimizes:

$$ $\displaystyle\bm{\psi}^{*}=\operatorname*{arg\,min}_{\bm{\psi}}\,\mathbb{E} \left[\mathcal{L}(y)\right]=\operatorname*{arg\,min}_{\bm{\psi}}\,\mathbb{E} \left[y\right]$ (6) $\displaystyle y\sim\mathcal{N}\Bigg{(}y;\ \sum^{n-1}_{i=1}\left[(\psi_{i+1}- \psi_{i}^{2})^{2}+(\psi_{i}-1)^{2}\right]+x,\ 1\Bigg{)},\ x\sim\mathcal{N}(x; \mu,1);\ \mu\sim\mathrm{U}\left[-10,10\right]$ (7) $$

We consider an episode terminated when $\mathbb{E}\left[\mathcal{L}(y)\right]=\frac{1}{N}\,\sum_{i=1}^{N}\mathcal{L}(f
_{\textrm{sim}}(\bm{\psi},\bm{x}_{i}))\leq\tau=3.0$, which we evaluate after every optimization step using $N=10^{4}$ samples. Following [], we use $\epsilon=0.2$ as the trust-region size and $\bm{\psi}_{0}=[\overrightarrow{2.0}]\in\mathbb{R}^{10}$ to initialize the optimization.

#### Nonlinear Submanifold Hump Problem

In this problem, we seek to find the optimal parameters vector ${\bm{\psi}}$ in $\mathbb{R}^{40}$ by utilizing a non-linear submanifold embedding represented by $\hat{{\bm{\psi}}}={\bm{B}}\mathrm{tanh}({\bm{A}}{\bm{\psi}})$, where ${\bm{A}}\in\mathbb{R}^{16\times 40}$ and ${\bm{B}}\in\mathbb{R}^{2\times 16}$.
To achieve this, we use $\hat{{\bm{\psi}}}$ in place of ${\bm{\psi}}$ in the Probabilistic Three Hump problem definition. Also, for the current setup, we follow similar settings as in []: the orthogonal matrices ${\bm{A}}$ and ${\bm{B}}$ are generated via a QR-decomposition of a random matrix sampled from the normal distribution; we use $\epsilon=0.5$ as the trust-region size and initialize the optimization at $\bm{\psi}_{0}=[2.0,\overrightarrow{0.0}]\in\mathbb{R}^{40}$.

#### Parameterized x-dist.

In order to evaluate the generalization capabilities of our method, we further parameterize each target function by placing distributions on the bounds of the Uniform distributions from which $x_{1}$ and $x_{2}$ are sampled. We randomly sample new bounds in every episode to ensure that the policy is exposed to multiple related but distinct simulators during training and evaluation. Concerning the Hump problems, we sample the lower and upper bounds of $x_{1}$ from $\mathcal{N}(-2,0.5)$ and $\mathcal{N}(2,0.5)$, respectively. For $x_{2}$, we instead use $\mathcal{N}(0,1)$ and $\mathcal{N}(5,1)$. For the Rosenbrock problem, we sample the lower and upper bounds of $x$ from $\mathcal{N}(0,2)$ and $\mathcal{N}(10,2)$, respectively.
Occasionally, an episode may not terminate as the specified termination value $\tau$ is below the minimum loss value for some samplings.

### D.2 Real-world Simulators

#### Wireless Communication: Indoor Transmitting Antenna Placement

The goal in this scenario is to find an optimal transmit antenna location ${\bm{\psi}}$ that maximises the signal strength over multiple receiving antenna locations ${\bm{x}}\sim q({\bm{x}})$.
Now, we detail aspects on the experimental setup for the experiments.
We run wireless simulations using Matlab’s Antenna Toolbox [], by evaluating the received signal strength (sigstrength function).
The simulations are run in two 3d scenes (‘conferenceroom’ and ‘office’), both of which are available by default and we additionally let Matlab automatically determine the surface materials.
We use the ‘raytracing‘ propagation model with a maximum of two reflections and by disabling diffraction.
The end-objective is to find a transmit antenna location ${\bm{\psi}}$ maximize the received signal strength over locations ${\bm{x}}\sim q({\bm{x}})$.
We constrain the locations in a 3d volume spanning the entire XY area of the two scenes: 3$\times$3m in conferenceroom and 8$\times$5m in office.
The transmit elevations $\psi$ are constrained between 2.2-2.5m and 3.0-3.2m per scene, and the receive locations between 1.3-1.5m (identical for both scenes).
The end-objective is to identify a transmit location ${\bm{\psi}}$ such that the median receive signal strength is maximized over a uniform distribution of receive antenna locations $q({\bm{x}})$.

### D.3 Reward Design

The reward function is chosen to incentivize the policy to reduce the number of simulator calls. This is achieved by giving a reward of -1 every time the policy opts to call the simulator, contrasting with a reward of 0 when it does not. However, with this reward function the policy could achieve maximum return (of zero) by never calling the simulator even if this leads to non-terminating episodes. An extra term is required to make any non-terminating episodes worse than any terminating one. Since the minimum return is $-L$, corresponding to the maximum number of simulator calls for an episode, that is achieved by setting a reward penalty of $-(L-l_{t})-1$ whenever the episode ends for reasons other than reaching the target value $\tau$: if the simulator call budget has been exhausted, then $l_{t}=L$ and the penalty is $-1$; if the timestep budget has been exhausted, then we have accumulated $-l_{t}$ return already. In both cases, adding this penalty leads to a total return of $-L-1<-L$.

### D.4 Termination Value

Termination values $\tau$ for the Probabilistic Three Hump and Rosenbrock problems are chosen to trade-off episode length and optimisation precision. Selecting values very close to the exact minimum value of the objective function $\mathcal{L}$ leads to extremely long episodes, due to the stochastic nature of the optimisation process. Moreover, parameterizing the distribution of the ${\bm{x}}$ variables changes the (expected) objective value minimum, such that choosing a too low value for $\tau$ leads to episodes that cannot terminate even in theory. Computing the minimum of $\mathcal{L}$ on the fly for the various parameterizations of ${\bm{x}}$ is not trivial, and so we opted for choosing a $\tau$ that generally suffices for good performance across parameterizations of a given problem. These values are chosen by manually inspecting L-GSO runs.

### D.5 Metrics

As we mentioned in Section 5, we use two metrics to compare our models against the baselines: the Average Minimum of the Objective function (AMO) for a specific budget of simulator calls and the Average Number of simulator Calls (ANC) required to terminate an episode. We now delve deeper into both of them. The meaning of the latter is quite straightforward. We consider the average number of simulator calls to solve the problem. We compute the average across evaluation episodes and random seeds. In contrast, the AMO is slightly less intuitive to interpret. One might question whether the value of the ANC should align with the maximum value on the $x$-axes for the AMO. In other words, assuming that for a given model, the ANC is equal to, e.g. 10, *should one expect that at a value $x=10$, the AMO will be equal to the termination value?*
Generally speaking, the answer is *no*. To explain why that is the case, we can report the following example. Let us assume that, for a given model, we have the following three episodes, each characterized by a specific length and value of the objective function at each simulator call:

- •
Episode 1: [20, 12, 7, 5, 3, 1]
- •
Episode 2: [18, 6, 1]
- •
Episode 3: [15, 5, 1]

We assumed the target value, $\tau$, to be 1. For simplicity, we used integers for the objective value. As we can see from the example, we have $\mathrm{ANC}=4$. Now, if we examine the AMO for $x=4$, we find that it is equal to 5 since only the first episode contributes to it, which is greater than $\tau$. Therefore, one cannot directly map the $x$-axis from the AMO to the $y$-axis of the ANC. Such a one-to-one mapping would exist only when all episodes always require the same number of simulator calls, which is not the case. We hope that our explanation has clarified the interpretation of the results we reported in the main corpus of the paper.

## Peer Reviews

confidence: 4
contribution: 1
initial_rating: 5
presentation: 3
questions: The decision to train on only local samples makes sense since one only requires local gradients, but there may be common features of the simulator's input-output relations across parameter space. Would it make sense to somehow incorporate older samples as well? Are there obvious cases where one or the other strategy would be preferred? How should one set/explore the size of the memory buffer when applying this approach to a new task?

The study seems to use the same policy network throughout, so we're assuming that the simulator's input/output mapping changes quickly, but the conditions under which we should trust our surrogates changes more slowly? Perhaps this contrast could be discussed in the context of the assumptions we should make about the simulator's local and global properties in order for the proposed method to work well.

Would it make sense to provide the observations $y$ to the policy network?

What are the failure modes? If I try this method out on a new task and the policy or surrogate training fails to adequately learn, how can this be identified? What if I don't know what loss value can be achieved/expected in my system a priori?

Is there any reason to use only 2-layer networks, especially for the 42-dimensional parameter estimation task?
soundness: 2
strengths: The approach is straightforward and intuitively appealing.

The presentation and structure and generally clear and informative (with a few exceptions noted below).

The ability to efficiently estimate simulation parameters is clearly important to many fields, and using local information is bound to be a winning strategy on some problems.
summary: The authors introduce an strategy combining RL and neural network simulation surrogates to optimize simulation parameters. The basic idea is to train local surrogates to estimate simulator parameter gradients, while training a policy network that tells us whether to trust the surrogates or to retrain them on a new batch of simulations. Using this approach, they aim to accurately estimate simulation parameters while minimizing the number of calls to the simulator. They test this approach on a combination of benchmark and real-world tasks.
weaknesses: The fact that a single "simulator call" means calling $f_\text{sim}$ thousands of times should be prominently stated early on, and not hidden in the appendix. The figures give a false impression that the algorithms used here can work with dozens of $f_\text{sim}$ evaluations, while in fact they appear to require thousands, similar to SBI methods that do not rely on local gradient estimates. The total number of simulation calls required to train the policy network should also be reported. If we simply trained a network to map from observations to parameters using a least squares loss and this number of samples from the parameter prior, then evaluated the loss, would it actually be worse than the reported results for inference (while using zero simulator calls at inference)?

More generally because the simulation costs are comparable, this method should be directly compared to SBI methods relying on neural density estimations, for which there are mature implementations (Tejero-Cantero 2020) and benchmarks (Lueckmann 2021). Can the proposed method solve these problems, and how sample efficient is it compare to SNPE, SNLE etc.? Because the software is readily available, this should be straightforward to verify. Total simulation counts, including those required to train the policy network, should be included. It is also worth discussing the fact that an amortized, trained SBI method can carry out inference on the full posterior without any additional simulator calls. In principle there must be some important scenarios where using local information is more efficient, but this does not seem to have been demonstrated here.

The strength of using local gradients in $\psi$ space over other approaches such as neural density estimation or Bayesian optimization is that this should scale to higher-dimensional parameter spaces. But the highest we see here is 42, well within the reach of existing SBI methods. Furthermore, we never get to see the joint distribution of true and estimated parameters, nor do we get to compare simulation results from true vs. estimated parameters. Therefore, it's hard to say whether the tasks (especially the real world problems) have actually been solved, or whether the current method is simply failing to solve it in a more sample-efficient way than baselines.

It's not so clear why the latent variable $z$ is included, when only the mean output of each surrogate is used to estimate sigma. Do the distributions from individual surrogates show variation, and does this capture uncertainty? Or does including $z$ somehow regularize the surrogates? Do they actually use $z$ or discard it after training?

The section on metrics isn't very clear. Does "objective function" refer to the reward function for RL or the loss for local training? It's not clear in either case how closely either of these correspond to accurate estimation of $\psi$, which is the real aim here. Why not show accuracy of $\psi$ vs. total simulation count? While it's generally true that lowering the objective function means we've gotten closer to the distribution of observables under the ground truth parameters, exactly how close is difficult to determine from the reported metrics. Also, it's not clear to me whether AMO/ANC are calculated over the whole training process, or for a new episode when the policy has already been learned.

Overall the precise sequence of steps is not clearly stated in the main text, and inclusion of a simplified algorithm in the main text would be welcome. Giving a concrete example of what the surrogate training loss might be early on in the manuscript could help readability.

Some relevant citations are missing. Papamakarios and Murray should be cited for SBI. Nonnenmacher (2021) should be cited for emulator-based gradient estimation. Lueckmann (2019) should be cited for simulation-based inference with active learning and ensemble predictions.

References

Lueckmann, Jan-Matthis, et al. "Likelihood-free inference with emulator networks." Symposium on Advances in Approximate Bayesian Inference. PMLR, 2019.

Lueckmann JM, Boelts J, Greenberg D, Goncalves P, Macke J. Benchmarking simulation-based inference. InInternational conference on artificial intelligence and statistics 2021 Mar 18 (pp. 343-351). PMLR.

Nonnenmacher M, Greenberg DS. Deep emulators for differentiation, forecasting, and parametrization in Earth science simulators. Journal of Advances in Modeling Earth Systems. 2021 Jun;13(7):e2021MS002554.

Tejero-Cantero A, Boelts J, Deistler M, Lueckmann JM, Durkan C, Gonçalves PJ, Greenberg DS, Macke JH. SBI--A toolkit for simulation-based inference. arXiv preprint arXiv:2007.09114. 2020 Jul 17.

Papamakarios G, Murray I. Fast ε-free inference of simulation models with bayesian conditional density estimation. Advances in neural information processing systems. 2016;29.

---

confidence: 3
contribution: 1
initial_rating: 3
presentation: 2
questions: * What is the motivation for using an ensemble to build the $\sigma$? It does not seem to be explained whether it is effective compared to not using it.
soundness: 1
strengths: * The experiment relating to the antenna placement is interesting and could be the focus of an applications paper if expanded on.
* The motivation to work on black-box optimization and gradient-free approaches makes sense.
* The introduction is written well.
* It is good that the comparison is made with an ensemble of L-GSO models.
summary: The paper focuses on black-box optimization problems using a learnt RL policy. The policy is then used to guide the training of differentiable surrogate models, where the gradients of these surrogate models are used to perform black-box optimization. The paper compares to a local-surrogate based approach and a Bayesian optimization approach.
weaknesses: * The main weakness of the paper seems to be highlighted in the appendix, where algorithm 1 shows that the policy is required to be trained on simulated data prior to the actual black-box optimization in algorithm 2. However, this is not made clear in the main paper. This means that the results in the paper do not include the large number of simulations required to train a policy and seems to defeat the purpose of reducing the cost of simulation. Specifically, to train the policy, the simulator is queried in a triple nested for-loop, which implies the requirement of many queries before performing the black-box optimization. 
* The section “Policy-Based Black-Box Optimization” does not highlight the technical approach of the paper. It continuously refers to the appendix instead of introducing the technical content. It should be the case that this is the most important section of the paper and therefore information regarding training and policy implementation (when it is key to the algorithm) needs to be made explicit in this section.  
* The uncertainty feature is introduced and then not mentioned again. Is it necessary? Especially since the appendix mentions using an ensemble of 100 models. Therefore there is the potential for large overhead. 
* It is not clear from the paper what the definition of a “parameterized input distribution” is, and why it might even be needed. Therefore, it is hard to interpret the difference between columns (a) and (c) in Figure 3.
* The results in Figure 3 do not match Shirobokov et al. (2020) which is a concern. For example at 100 function calls in Figure 2a of the prior paper, no algorithm has reached -0.8, but all algorithms have reached -0.8 in this paper except numerical differentiation. The difference can also be seen for the 10D Rosenbrock and the nonlinear submanifold hump problem. This needs to be addressed in the paper as to why there is this discrepancy.   There is also a similar discrepancy for the Physics experiment where L-GSO gets to an objective value of about 70 after 20 calls in Shirobokov et al. (2020). 

Minor:
* Is “learning active learning” an accepted phrase in the literature? Maybe meta active learning as in Ravi and Larochelle (2018)?

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: - how are the surrogate neural networks trained?
- what are the values employed for K and G?
- Do the plot reporting the number of simulator calls incorporates the cost of training? This is key for a fair comparison in my opinion.

Minor:

- Typo l. 403: "receive locations" -> "receiving locations"
soundness: 2
strengths: - I am not an expert in the area, but the problem of amortizing an active learning strategy over multiple optimization is an interesting problem (although I am not sure how relevant this problem is in practice).
- The method is tested on multiple real-world experiments, and is able to handle stochastic and non-stochastic simulators, suggesting a potential high impact for practitioners.
- Aside from the issue related to amortization (discussed above), the presentation is very clear.
summary: This submission proposes a method to perform optimization of a loss function that takes the form of an expectation under a sampleable simulator, but whose (stochastic) gradient gradient is not analytically available.

To do so, the method follows an approach already explored in the field approximately minimizing the loss function using Stochastic Gradient Descent, where the gradient is obtained using a differentiable surrogate of the simulator. To produce a gradient estimate, at each gradient step, a policy network  determines whether to train a new simulator surrogate.
This policy network is trained in a separate algorithm prior to the optimization using a reward model which incentivises both reaching satisfying objective function values and minimizing the number of surrogate retraining operation, as the latter involves sampling from the simulator which may be costly in practice.

Variants of the method are constructed where the policy is also trained to estimate the width of the parameter neighborhood used to retrain the surrogate.

The method is showcased on a series of benchmarks and on real world experiments. In the case of the benchmark model, an additional "parameterized input distributions" setting is considered, where the goal is to be able to perform optimization of multiple loss functions, each of them with the same simulator function, but differing in which input distribution is used. The training of the policy network is amortized over all parameterized distributions.
weaknesses: If I understood correctly, I see two ways in which the method could be beneficial:

1. either the total number of simulator calls over the 2 phases ((1) training the policy network and (2) optimizing the function) is smaller than the number of simulator calls in  other methods.
2. or, as training the policy network can be amortized over multiple input distributions q_i, if one needs to solve multiple such optimization problems using different input distributions, the cost of training ends up being negligible compared to the gains at inference (compared to other methods) that the training policy network provides.


However, both of these benefits are neither clearly described nor demonstrated in the main body:
- regarding point 1: when reporting the metrics, it is not mentioned whether training is accounted for
- regarding point 2: the use of multiple input distributions is mentioned, but once (l 326), and without stating out clearly why it is an important use case for the method at hand. If point 2 ends up being the main advantage of the method, I would expect the introduction to provide more motivation as of why parameterized input distributions matter in practice. Moreover, I could not find a mention of parameterized distributions in real-world experiments.


In my opinion, the paper needs to demonstrate clearly that the additional training step is justified. Due to the reasons described above, this demonstration remains limited, and thus I will recommend rejection until these issues are addressed. I may have missed additional benefits, in which case I will change my recommendation.

Some important training and evaluation details are missing from the main body, see questions.