Bayesian Optimization of Function Networks
Authors: Raul Astudillo
Peter I. Frazier
Venue: NeurIPS 2021 Poster (2021)
Keywords: Bayesian
optimization
network
composite
nested
grey
box

Abstract: We consider Bayesian optimization of the output of a network of functions, where each function takes as input the output of its parent nodes, and where the network takes significant time to evaluate. Such problems arise, for example, in reinforcement learning, engineering design, and manufacturing.  While the standard Bayesian optimization approach observes only the final output, our approach delivers greater query efficiency by leveraging information that the former ignores: intermediate output within the network. This is achieved by modeling the nodes of the network using Gaussian processes and choosing the points to evaluate using, as our acquisition function, the expected improvement computed with respect to the implied posterior on the objective. Although the non-Gaussian nature of this posterior prevents computing our acquisition function in closed form, we show that it can be efficiently maximized via sample average approximation. In addition, we prove that our method is asymptotically consistent, meaning that it finds a globally optimal solution as the number of evaluations grows to infinity, thus generalizing previously known convergence results for the expected improvement. Notably, this holds even though our method might not evaluate the domain densely, instead leveraging problem structure to leave regions unexplored. Finally, we show that our approach dramatically outperforms standard Bayesian optimization methods in several synthetic and real-world problems.

Bayesian Optimization of Function Networks
Raul Astudillo Cornell University 2a598@cornell.ede
Peter I. Frazier
Cornell University
#F98@cornell.ede
Abstract
We consider Bayesian optimization of the output of network functions, where each function takes input the output of its parent nodes, and where the net- work takes significant time to evaluate. Such problems arise, for example, in reinforcement learning, engineering design, and manufacturing. While the stan- dard Bayesian optimization approach observes only the final output, our approach delivers greater query efficiency by leveraging information that the former ignores: intermediate output within the network. This is achieved by modeling the nodes of the network using Gaussian processes and choosing the points evaluate using, as our acquisition function, the expected improvement computed with respect to the implied posterior on the objective. Although the non-Gaune non-Gaussian nature of this posterior prevents computing our acquisition function in closed form, we show that it can be efficiently maximized via sample average approximation. In addition, we prove that our method is asymptotically consistent, meaning that it finds a globally optimal solution as the number of evaluations grows infinity, thus generalizing previously known convergence results for the expected improvement. Notably, this holds even though our method might not evaluate domain densely, instead leveraging problem structure to leave regions unexplored. Finally, we show that our approach dramatically outperforms standard Bayesian optimization methods several synthetic and real-world problems.
Introduction
We consider Bayesian optimization (BO) of objective functions defined by a series of time-come consuming- to-evaluate functions, fi..... .......fk. fK. arranged in directed acyelic network, so that each function takes as input the output of its parent nodes. As we detail below, these problems arise in BO-based policy search reinforcement learning (Lizotte al. 2007), optimization of complex systems modeled simulation, and calibration of time-consuming physics-based models.
To illustrate, introduce running example of vaccine manufacturing (Sekhon and Saluja 2011 1), focusing on the portion of the manufacturing process that uses live cells to produce proteins needed in vaccine. begins with cell culture, which living cells grown and used "factories" to produce proteins. This process is controlled by vector, 31, containing the temperature, pH, and CO. content used when growing these cells. The output of this process is the quantity of the desired protein yı = fa (x1), i.c. the yield of this step, along with other byproducts. This output is passed into second process, purification, which removes byproducts and is controlled by a vector I2 comprising temperature, pressure, and flow rate. The yield of the desired protein from this second step 1/2 f2(x2, v1). This output enters a third step, formulation, in which we formulate the raw protein into form that can be distributed controlled by a third set of parameters. This determines the yield of the overall process ys f1 y2). We wish choose (æ1, 22, 233) maximize overall function network in Figure [
protein yield. This problem is summarized
The problem described above and other similar problems can be tackled with Bayesian optimization (BO), which has been shown to perform well compared to other derivative-free global optimization
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

methods for time-consumen consuming-tor torenng--o--vvlatte objective functions (Snock et al., 2012 Frazier 2018). A standard BO algorithm would fit Gaussian process (GP) (Rasmussen and Williams (006) model on the objective function (Vs, which depends on (I1, 22, F33) and use it, along with an acquisition function, to sequentially choose the points to evaluate. Under this standard approach, however, evaluations of the intermediate nodes, Ji.....Kk-1. would be ignored despite being available when computing the objective function. In the example above, this corresponds to looking only at the yield of the overall process, and not of each individual step.
y1
y/2
1/a
fi
f2
f3
I
22
T'3
Figure 1: Vaccine manufacturing as a function network. Protein = [(1) is created, then purified with yield y22 = f2(22, and formulated with yield ys fs( 1/2)- The goal is to find (21, 22, v3)) that maximizes Vs-
In this paper, we introduce a novel BO approach that leverages function network structure for substantially more efficient optimization. This approach models the individual nodes of the network using distinct GPs. This allows incorporating observations of each node's output recursively into non-Gaussian posterior on the network's overall output. Our approach then chooses the points to evaluate using the expected improvement (Jones et al. 1998) computed with respect to this implied posterior on the objective function. The non-Gaussian nature of this posterior prevents the expected improvement from having a closed form. However, we show that it can still be efficiently maximized via sample average approximation (Kleywegt et al. 2002).
Our approach can outperform standard BO by leveraging information from internal nodes unavailable to standard methods. We briefly explain one way this can happen, in the context of the example above. In vaccine manufacturing, each function fIe 1)- k = 2, 3. is bounded between 0 and U because new protein cannot created in the purification and formulation stages. Moreover, application experts have a prior on what values for fa (2 02-1)/ 1/yk-1 should be achievable if is set well. Thus, if we see that ya is unexpectedly poor, information from intermediate nodes can be extremely valuable: if yı and y/2 are as expected, then this suggests the problem is with z3; if y as expected but y/2 is poor, then the problem likely with 22; and if y/ is poor then the problem is likely with ø1- (If there is problem with Iı. there may also be a problem with Tw' k > k. but we can focus on fixing zg first.) Thus, by observing intermediate nodes, we can instantly reduce the effective dimensionality of the input space by factor of K = 3.
We show that our method asymptotically consistent, i.e., that it discovers the global optimum given sufficiently many samples. Remarkably, in contrast with most BO methods, it may do so without measuring densely over the feasible domain, instead leveraging function network structure to exclude regions as unnecessary explore. This indicates the power of function network structure to improve query efficiency.
We demonstrate through numerical experiments that access to additional information available in a problem formulated function network can dramatically accelerate optimization. We study four synthetic problems and four real-world problems: a manufacturing problem similar in spirit to the vaceine example above, an active learning problem with robotic arm, and two problems arising in epidemiology, one calibrating an epidemic model and the other designing a testing strategy to control the spread of COVID-19. Our method significantly outperforms competing methods that utilize less information, in some cases by ~5% and in other cases by several orders of magnitude.
2 Related Work
Our work occurs within BO, a framework for global optimization of expensive-to-evaluse black-boxx functions that originated with the work of Zhilinsks (1755) and Močkus (1975), and has recently become popular due to its remarkable performance in hyperparameter tuning of machine learning algorithms (Snoek et al. 2012) (wersky et al. 2013) Wu al. 2019). We refer the reader to (Shhhiiaii et al. 2016) and Frazier (2018) for modern introductions to BO.

---

Our approach can be catalogued as a grey-box BO method since it does not treat the objective function entirely as black box, and instead exploits known structure to improve sampling efficiency.
Other examples of grey-box BO approaches include
BO (Kandasamy et al. 2017 Wu
jet al. 2019), which leverages cheaper approximations of the objective function; BO of objective functions that are the integral of an expensive-to-ewe to-evlluute integrand (Williams et al.) 2000) Toscano Palmerin and Frazier 2018: Cakmak et al. 2020); BO of objective functions that are a sum of squared errors (Uhrenholt and Jensen 2019); and, more generally. BO of objective functions that are the composition of an sxpensive-to-evallee
inner function and a known inexpensive-to-evale function (Astudillo and Frazier 2019). of which our work can be seen significant generalization. We refer the reader to Astudillo and Frazier (021) for a tutorial on grey-box BO.
outer
Our work is also closely related to Marque-PPucher et al. (2019), which proposes a method for efficient sequential experimental design of nested computer codes, also using GPs. This work focuses on the case where there are only two node functions, and one takes as input the output of the other. In contrast with our work, the goal of the proposed method is to learn the output code as accurately as possible within limited evaluation budget, but optimization is not pursued.
Optimization of composite (a.k.a nested) functions has also been considered in the gradient-based optimization literature (Shapiro 2003 Drusvyatskiy and Paquette| 2019 Charisopoulos et al. 202T Balasubramanian et al. 2020). In contrast with ours, these works assume that evaluations are inexpensive, and often also some form of convexity, along with availability of gradients. However, this literature is similar in spirit to ours in that information from inner functions improves efficiency.
Function networks arise in many application areas. While these applications have not, to our knowledge, been previously formulated as specific instances of the general function network model we propose, their literatures are nonetheless relevant to our work.
In engincering and aerospace design, function networks arise in multidisciplinary opti- mization (Cramer et al. 1994] Amaral et al., 2014; Benaouali and Kachel 2019), where simulators focusing on different physical laws are coupled into a function network. For example, a simulation of airflow over a wing may output the forces on the wing to another simulation of mechanical stress on the wing structure.
In drug discovery and materials design, function networks arise from the sequence-structure function (Sadowski and Jones 2009) and cnapoostiioo-stuuctue--reeee
Sestrick-Simpers
et al. 2016) paradigms. Here, decision variables (composition, e.g., what fraction of what raw materials are used) determine the structure of the material (how the atoms combine together), which in turn determines how the material behaves (properties).
Function networks arise in the design of queuing networks (Fu and Henderson 2017 Arcelli 2020). This includes manufacturing systems (Ghasemi et al., 2018), where the partially- processed output of one workstation is input the next workstation, the design of service systems (Wang et al., (2020) such as hospitals and airport security checkpoints, and choosing traffic signal timings for a city's street network (Osorio and Bierlaire) 2013).
Finally, function networks arise in reinforcement learning (Sutton and Barto 2018) and Markov decision processes (Puterman 1990), where the transition kernel translorms the state variable the start of time range to another state variable at the end of this range. This outputted state variable is the input to the next time range. (5.2]s an example.
3 Problem Setting
We consider objective functions evaluated via a series of functions, fi, fk. arranged in directed network so that each function in this network takes as input the output of its parent nodes, and assume that evaluating this collection of functions is consuminn. The network structure is encoded as directed graph with nodes V = {1..... K} and directed edges E = { k) f& take as input the output of f,}. We assume that this graph is acyclic and has a single leaf node whose output is the objective to optimize.
Let J (k) denote the set of parent nodes of node We assume, without loss of generality, that the nodes are ordered such that j < k for all j E J (k)..
J. In addition to the output of its parent nodes, we
We allow ((k) 0. which corresponds to the root node(s).
3

---

assume that each function, fu. takes as input (potentially empty) subset, (lk) C {l..... D). of the components of the decision vector x E R."
Let ha (x)..... h (a:) denote the values of the K nodes in the function network when it evaluated xc These values defined recursively by
h₂(2)) fk (o (x)a)())) k k= l.=.......
ie.. by evaluating each function in the network as the values of its parent nodes become available. This recursion is well-defined by our assumption that the graph is acyelic. The objective function is then g hg, and the goal is to solve
max g(z), ZEX
where X C R.D is a simple compact set, such hyper-rectangle.
The standard BO approach to this problem models g using a GP prior distribution. This approach iteratively chooses the next point at which to evaluate as follows. Given n observations of the objective function, g(x1).... : ). it computes the posterior distribution on g. which itself GP. It then uses this posterior distribution to compute an acquisition function (Frazier 2018) that quantifies the value of the information that would result from a function evaluation at any given point. Finally, it chooses the next point to evaluate, Ix+1, as the point that maximizes this acquisition function. Notably, although the outputs of fi..... fk-1 would be observed as part of this approach. these evaluations would be ignored by standard BO when calculating the posterior distribution and corresponding acquisition function.
g(xx),
Bayesian Optimization with Full Network Observations
This section describes our approach. Like standard BO, it consists of a statistical model and an acquisition function. Unlike standard BO, however, our approach leverages the network structure of the problem by utilizing intermediate outputs within the network. As we describe below, this is achieved by modeling the node functions, f1..... .kk, as GPs, which in turn implies a non- Gaussian distribution on g ($4.1). Our acquisition function the expected improvement under this posterior distribution ( (.4.)) which no longer has a closed form and thus we maximize it via sample average approximation ($4.3). We end up this section by proving that our acquisition function is
asymptotically consistent despite not necessarily sampling X densely ($4.4)
4.1 Statistical Model
Instead of modeling g directly, we model f fK. as drawn from independent GP prior distribu- tions. Let & and Σg, denote the prior mean and covariance functions of fk, A 1,.... K, respec- tively. When g is evaluated at z, also get to observe the value of fx at (II(K). h, (I)) Thus, after querying at n points, Ze, f 1 I...., n, the observations of the values of fk. 1: 1,...,K, K, at (z), h()()) (z2)), e ... n, imply posterior distributions fi,..., fK. which (conditionally) independent GPs with mean and covariance functions fn.k and Em k I.....K K. These can be computed in closed form using the standard GP regression equations (see, e.g., Ras mussen and Williams 2006). For completeness, we include these equations in of the supplement.
again
The posterior distributions on fi...., described above imply a posterior distribution on This distribution is in general non-Gaussian. Thus, unlike in the standard setting, where g is modeled as GP, classical acquisition functions such as the expected improvement cannot be computed in closed form. However, as we describe next, drawing samples from this distribution is simple.
Thanks to the acyclic structure of the underlying network that defines g, a sample, ((x) 7u (a) from the posterior distribution on g at x can be obtained following the iterative process described next. In each iteration, h = 1...., K, we obtain a sample, hg (z), from the posterior distribution on (xx by drawing a sample from the normal distribution with mean jlnkk (zI(K), (x)())()))) and
For ease of presentation we assume that all the functions in the network are expensive- expensive-to-evaluse to-evaluate, and thus require to be modeled as GPs. However, if any of the functions, say ft. is not expensive-to taaeiee-----vvvatteee we can simply take μo fy and Σx = 0.

---

standard deviation
Jn,k (ZI(s). hoss)) En,i (oT()), i (2) (10)) (o))))))())))))))))))))))))))( 1/2
Supporting efficient calculation, the samples ha (a)) will have already been obtained in previous iterations of the for loop since j < k for all j (kk (Note that J(1) (). This procedure is summarized in Algorithm
Algorithm 1 Draw sample from the posterior g(z))
Require: x E X
1: for k = I....,Ay K do
2: (x(2) N Honk (erc)) on. (x10))))00)(()))
3: end for
4: return ((x) = lg((x)
We end this section by noting that, while the statistical model described above could be catalogued as a deep GP (Damianou and Lawrence 2013), in the sense that we have GP layers in an architecture described by directed acyclic graph, inference in our model is faster. Typically, deep GP inference requires marginalization over latent values of GP layers. In our setting, however, observation structure creates conditional independence across layers, avoiding the need to marginalize.
4.2 Expected Improvement for Function Networks
Our acquisition function is the expected improvement, computed with respect to the implied posterior distribution on g under the statistical model described in $4.1
EI-FN,,()) = E, [[((>)) g4
where 91 max;=l..... (Iin) is the best value observed so far, and E., is the expectation computed with respect to the GP time-n posterior distributions on fi....,Jx fk- To distinguish it from the classical expected improvement, refer for function networks (EI-FN) in the remainder of this work.
acquisition function the expected improvement
4.3 Maximization of EI-FN via Sample Average Approximation
Since the posterior distribution on g is non- non-Gaussian, in contrast with the classical expected improve- ment acquisition function, EI-FN,, does not admit a closed form expression. However, EI-FN,, (z) can be computed with arbitrary precision following simple Monte Carlo (MC) approach:
M
EI-FN. (z) R
Σ m=1
((x)((m)
(m)
In
M
where ()).,, (M) are samples drawn from the posterior distribution on g(z), which can be obtained following the approach described in $4.1
Following the above MC scheme to approximately compute EI-FN, one can aim to maximize this function using derivative-free global optimization algorithm for meexpesivve-to-evaluat functions, such CMA (Hansen 2016), by drawing samples, ())) ... (ar) independently for each at which EI-FN,, is evaluated. However, this approach slow since evaluations EI-FN are noisy and derivative information is not leveraged. Instead, propose to maximize EI-FN following
sample average approximation (SAA) approach (Kleywegt et al. 2002 Balandat et al. 2020).
Succinctly, the SAA approach works by building MC approximation of EI-FN that is deterministic given finite of random variables not depending on æ, and maximizing this approximation instead of EI-FN itself. The key observation for creating this approximation is that a sample, ÿ(r), can be
obtained (xx) (x), where h (x), ... Tog (x) are defined recursively by
Th (x; 2) pr.e (ZI(M), Prs (x; 2)) + on.k (er(s))T 2)) Zk,
5

---

where Z = (Z)..... Zk)
is standard normal random vector, and we write (x)))) as hy (æ; Z) to make its dependence on Z explicit. Analogously, we also write g(x) as g(:; Z). This can be seen as an extension of the so-called reparametrization trick for acquisition functions (Wilson et al. (0118).
If now fix M samples drawn from the K-dimensional standard normal distribution, z(1) ..... MMM, then
M Σ
EI-FN, (I; ((.M
1 M
{g xx Z(:)) g m=
is MC approximation of EI-FN that is deterministic given Z(1,, ..... Z(ar) as desired. Moreover, Proposition Proposition[]bellv [ellow shows that, under mild regularity conditions, any maximizer of the above SAA converges in probability exponentially fast to maximizer of EI-FN, M + oc, thus suggesting that in practice it is safe to use small values of M. This result is a generalization Theorem in Balandat et al. (2020). Its proof can be found in of the supplement.
Proposition 1. Suppose that the functions Pr,k and On,ko k = I,..,, K, are all Lipschitz continuous and let
E aargmax EI-FN, (x; (0.7))))
X, argmax EI-FN,,
.EX
xEX
Then, for every e > 0, there exist A, a > O such that P dist
( x.) > >() & Ae-aM, MEN.
Finally, we note that EI-FN, is differentiable with respect to x, provided that Hk,n and Ei,n h 1..... K. are all differentiable. Thus. EI-FN, can be maximized using gradient- based deterministic optimization algorithm such as L-BFGS-B (Byrd al. 1995). with multiple restarts.
4.4 Asymptotic Consistency of EI-EN without Dense Measurements
To shed light on the convergence properties of the expected improvement acquisition function under statistical model, we prove that, under suitable regularity conditions, EI-FN is asymptotically consistent, meaning that it finds the global optimum of the objective function as the number of evaluations grows to infinity. This builds on results shown for the classical expected improvement
(Vazquez and Bect 2010 Bect et al., 2019) and the expected improvement for composite functions (Astudillo and Frazier 2019). and we rely on similar assumptions. This result is stated Theorem below and its proof can be found in of
of the supplement.
Theorem 1. Suppose that En+1 E arggaaxexx (I- for all n. Then, under regularity condi- tions stated precisely in the supplement, In maxxex ((x) as - oc almost surely under the prior distribution on fi...., fK.
Significant novelty in our proof arises from the fact that EI-FN's measurements are not necessarily dense in X. In nearly all existing work, consistency of BO methods is shown by first showing that the measurements dense in the objective's domain (sce, e.g., Vazquez and Bect 2010). Surprisingly, the function network structure may allows us to exclude certain regions of X as suboptimal after only finitely many measurements, allowing our method to be consistent without measuring everywhere. Such ability gives insight into EI-FN's s strong empirical performance compared to methods ignoring function network structure. As stated in Proposition 2 below, we provide a function network where EI-FN provides a consistent estimate of the global optımum without sampling densely.
Proposition There exists a function network (detailed in of the supplement) in which EI-FN consistent but whose measurements are not necessarily dense in X.
5 Numerical Experiments
We compare the performance of our algorithm (EI-FN) against the classical expected improvement (EI), i.e., the expected improvement under a GP model over g. We also compare with the performance of two other standard algorithms: the algorithm that chooses the points to evaluate uniformly at random over X (Random); and the knowledge gradient (KG) (also under GP model over g), more sophisticated acquisition function that often delivers a better performance than the expected improvement (Wu and Frazier 2016). The problems in 5.2 and 5.3 fall within the framework

---

of Astudillo and Frazier (2019) and we include its proposed EI-CF algorithm as a benchmark. All algorithms were implemented in BoTorch (Balandat al., 2020).
We evaluate on 4 synthetic problems and real-world problems. These are described below or in $ of the supplement, with the supplement describing manufacturing problem building on I COVID-19 testing problem building on $5.3] and a robot control problem similar to the one described in $5.2 In all problems, a first stage of evaluations is performed using 2(d + 1) points chosen uniformly at random over X. A second stage (pictured in plots) is then performed using each of the algorithms. Error bars in Figure 2 show the mean of the best objective value observed so far, plus and minus 1.96 times the standard deviation divided by the square root of the number of replications. Since the difference in performance in some of our experiments is better appreciated in a logarithmic scale, we also include plots showing the logre-regret for such experiments in Figure 3 Each experiment was replicated 30 times. Experimental setup details and runtimes are available in the supplement. Code to reproduce our numerical experiments can be found at htttts::///gh
5.1 Synthetic Test Functions
We create synthetic test problems by arranging standard test functions from the global optimization literature (Jamil and Yang 2013) into function networks. These explore a variety of network structures in an sesy-to-reproduee form, and are named after the standard test function used to define the function network. We describe these briefly here and then in full detail in Of of the supplement.
Alpine2 and Rosenbrock both arrange K nodes in series, where each node except the first node takes the output of the previous node as input. Additionally, in Alpine2, each node takes a distinct dimension of the decision vector x as input. In Rosenbrock, 21 and 22 are inputs to the first node, I2 and 23 are inputs to the second, and so on. These network architectures arise in manufacturing problems like the example in $ as well as business operations with queues like boarding an aircraft or fulfilling drive-through orders. For Alpine2 we set K = 6, and for Rosenbrock we set K = 4.
Ackley has 3 nodes. The first two nodes each take the same 6-dimensional input. Their outputs are passed to the third node that produces the final output. This type of function network arises in algorithm design for two-sided markets (Li al. 202I). like Uber and AirBnB, where the first node simulates an intervention's 's effect on riders (or guests), the second simulates its effect on drivers (or hosts), and the third simulates the matching process where riders and drivers (or guests and hosts) interact to produce transactions.
Drop-Wave has two nodes. The first node takes a two-dimensional vector x as input. This node's s output is passed to the second node, which produces the objective value. This network architecture is representative of multidisciplinary engineering design (Benaouali and Kachel 2019), for example in aerospace. where a small number of distinct black-box simulators simulate processes governed by physical laws that affect each other through small number of channels, such as aircraft engine simulation (the first node) determining heat generated while flying, which is then inputted to a cemperature-depente simulation of mechanical stress on the aircraft's frame (the second node).
5.2 Fetch-and-Reace with Robotic Arm
This test problem is obtained by adapting the Fetch environment from OpenAl Gym (see Plappert jtt al. (2018)). The goal is to move the gripper of robotic arm to a target location with only three movements. We formulated this problem as function network with 3 nodes, each representing movement of the robotic arm. Each of these nodes takes as input the current location of the gripper along with a vector of forces to be applied to the robotic arm in that step, and produces output the location of the arm after this movement is complete. (Note that the output of each node is 3-dimensional and thus this can also be thought of as function network with single-output nodes). The objective to minimize is the distance between the gripper and the target in the final step. Figure 4 shows animation of this problem.
We formalize the above problem as follows. Let Zinit: zearget E R denote the object's initial and target locations, respectively. At each time step, I, we choose the vector of forces to be applied to the robotic arm I, E [-1, After this movement, the location of the object becomes 2411- The goal is to choose I; for I,...,T T to minimize || L|argee 2T 1|2- We set 2nnt (0,0,0).
7

---

drop wave
reserbreck. K-4
adkley
alpine2, K
0.9
bust dbjective value
bestt www.ccooeooooooooooooooooooooo waltow
-100
300
9.6
200
220
0.4
330
100
25
130
20
50
1.00
25
75
manuterturing lie
could-19 group lexing
1otch and reach
coldenic mode: calloration
4.4
0.0
4.2
3.50
best ob/ocht value
-9.75
-10.00
-0.6
Random
3.4
-0.8
KG
30
32
-10.50
-1.0
BIHN
number ovaluations
number Lacions
number of osaalaations
rumber evaluation
Figure 2: Top: Results on synthetic problems that adapt widely used synthetic test functions into function networks. Bottom: Results on realistic problems: manufacturing line, design of testing protocols for COVID-19, fetch-and-reach with a robotic arm, and calibration of an epidemic model. EI-FN substantially improves over benchmark methods, with larger improvements for problems with higher-dimensional decision vectors and more nodes.
rowembrock, K
covis-19 group testing
fetch-and--nach
model cullbnation
0.2
1.0
0.0
-0.2
0.0
Sogo (ww.et
-.44
Random
-1
-0.6
-1.0
-2
-0.8
BSCF
-3
-1.0
-2.0
BIFN
50 evalustions
number evaluations
evaluttoos
number evaliations
Results on four of our test problems. contrast with Figure above, which shows the best objective value found, here we plot the logre-regret.
Figure 3:
Starget = (-12, 0.2), and T 3. This can be interpreted as a function network by associating each time step with triplet of node functions fc f1, which take input and produce 1 = fxx))) as output.
A very similar experiment to the one described above can be found in $ of the supplement. It considers a variation of the active learning for robot pushing problem introduced by Waan and Jegelka (2017) whose goal is to teach robot to push an object to predetermined target location. As in the experiments here, EI-FN outperforms other methods significantly, including EI-CF.
5.3 Calibration of an Epidemic Model
Here we consider calibration of compartmental stochastic models to data, a widely-used tool in epidemiology, medical modeling, and ecology (Sandberg 1978). Function networks well-suited to exploit the structure of such models. We focus on calibration of a specific epidemic model for influenza, building toward COVID-19 mitigation benchmark in the next section. We first describe the epidemic model, then the calibration problem and, finally, formulation as a function network.
SIS Epidemiological Model: We calibrate to data widely used epidemiological model, the SIS model (see, e.g.. Garnett 2002), that models diseases like influenza capable of reinfecting individuals multiple times. In this model, individuals either do not have the discase and are "susceptible" (S) or
8

---

Figure 4: A sequence of screenshots showing three consecutive movements performed by the robotic arm described in $5.2
have the disease and are "infectious" (I). We consider SIS model of two interacting populations, where infections occur at population-ddependen rates.
The model is dynamic, with time indexed by - 0, ,. At the beginning of each time period I, the fraction of population E {0, }} that infectious J,,,- We assume each population is of equal size, N. During this time period, each person in group i comes into close physical contact with S.,.c people from group j. When this contact is between an infectious person and a susceptible one, infects the susceptible person. A fraction Lj, of the people from group j involved in such interactions are infectious and a fraction (1 - 45,t) from group are susceptible. A number of new infections in group result, N(1 T I As fraction of group i's population, this is (1 J., Summing across j. we have (1 - I,,e) E, (t) new group infections. At the same time, infections resolve at a rate of per period. This results in a decrease in the infectious population in group of 1 Putting this together, the number of infectious individuals in group i at the start of the next time period is 4,,+1 = 1,,(11 - 7) + (1 - 1.,) E; B.5,a
Calibration: The SIS model has parameters 14,0 7 and Bi.j where 0 S t < T and i, E (1,2). We calibrate the parameters 3 = (Ba,, i.j. while fixing 7 0.5 and I,,0 .01 (for both i). We up T 3, using a held-out value for 3. We let robs
simulate a trajectory of infections from denote the fraction of group i observed to be infected at time t in this trajectory. We then search for the vector 3 that, when passed to the SIS model, minimizes the mean squared error between this trajectory and the SIS model predictions. Letting (3) indicate this predicted value, the goal is to minimize the mean-squared error (MSE), MSE = Eiz (1obs lr)(3)2 do not include t since (,) the same for all B.
Formulation as a Function Network: We encode this as a function network using 2T + nodes, illustrated in Figure 5 For each time period and each group i, node takes input I (I j = 0, 1) and 3e = (5 j.j' E {0, 13) and produces output I,,+1- (For 0, the input Jo is not needed since this is the same for all F.) Then, one additional node takes the output of the other nodes (Li, : = 0, 1, 1, I........ T)) as its input and produces the sum squared errors Ei - ( output. We treat this final node as known (its GP prior has a kernel of 0).
EI-CF benchmark: The fact that the final node in this problem (denoted "MSE" in Figure 5)) has known structurre permits comparing against the EI-CF method for BO of composite functions Astudillo and Frazier 2019) benchmark. EI-CF substantially less general than our method (EI-FN): it restricted to settings with one time-consuming black-bor multi-output node that provides input to one -o-eluuue node with known structure. To apply EI-CF to this problem, the single black-box multi-output node takes as input and produces the vector (3),)(3) : i, t) as output. This output is then supplied to the "MSE" node. This approach ignores the fact that l does not depend on 3r, i' > t, and depends only indirectly on Be, < t through 1j,t-1, j = 1, 2.
5.4 Diseussion
Across the wide range of problems considered, EI-FN significantly improves query efficiency over standard BO methods that ignore the function network structure of evaluations. The benefits range

---

L1
L2
L3
E L1
1.2
1.3
13
31
82
MSE
R3
R2
R3
2,
12.2
12.3
I1
22
Figure 5: Function network for the (left) epidemic calibration problem in CÓVID-19 pooled testing optimization problem described the supplement.
$5.3 and (right) the
from a 5% improvement in the value of the best point found on the Drop- Wave and manufacturing problems to several orders of magnitude in the Rosenbrock and epidemic model calibration problems.
The largest benefit arise when the control imput is ligh-dimensional but the input to individual nodes is low-dimensional. On the epidemic model calibration problem, we see EI-CF (Astudillo and Frazier 2019) outperforming El and KG by several orders of magnitude, and EI-FN outperforming EI-CF by several additional orders of magnitude. As noted above, EI-CF can be seen as a special case of EI-FN using a less informative function network that hides observations from some nodes. This consistent with observations of function network structure allowing substantial improvement in query efficiency, and observing more of the internal function network structure providing more value.
6 Conclusion
We introduced novel BO approach for objective functions defined by a series of expensive-to- evaluate functions, arranged in network that each function takes as input the output of its parent nodes. These objective functions arise in a wide range of application domains. However, existing methods cannot leverage this structure. Our approach models the outputs of the functions this network instead of only the objective function, as is standard in BO. Our experiments show that, by doing so, this approach can dramatically outperform standard BO methods.
Though we see substantial benefits from our approach, there are limitations. First, it requires more computation than standard BO methods, explored in the supplement. (When the objective is iime-conumiing, the improved query efficiency more than makes up for the additional computation required.) Second, while we have demonstrated our method in problems with up nodes, and computational speed would support more, our method does not (yet) scale to hundreds of nodes. Third, while we show consistency, it would be instructive to complement our empirical results showing fast convergence with a theoretical understanding of convergence rates. Existing approaches to prove convergence rates for the classical expected improvement heavily rely on properties of its analytical expression (Bull) 2011 Ryzhovv] 2016)), and thus are not directly generalizable to our setting. This is, however, an exciting direction for future work.
As with any powerful new method for optimizing time-consuming-to-comm coocpputte black-box functions, ours can accelerate many applications. While this includes innovations that generally benefit society, such as improvements to public health and vaccine manufacturing, it also includes the design of weapons and other engineering systems that could harm individuals. Thus, is important that society enact guardrails that ensure proper use of our methodology.
Acknowledgments
The authors were partially supported by AFOSR FA9550-19-1-0223 and FA9550-20-1-0351.
10

---

References
Amaral, S., Allaire, D., and Willcox, K. (2014). A decomposition- based approach to uncertainty analysis of feed-forward multicomponent systems. International Journal for Numerical Methods in Engineering, 100(13);922-1055
Arcelli, D. (2020). Exploiting queuing networks to model and assess the performance of adlf-ada software systems: A survey. Procedia Computer Science, 170:498-505.
Astudillo, R. and Frazier, P. I. (2019). Bayesian optimization of composite functions. In Chaudhuri, K. and Salakhutdinov, R., editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 354-363. PMLR.
Astudillo, R. and Frazier, P. I. (2021). Thinking inside the box: A tutorial on grey-box Bayesian optimization. In Proceedings of the 2021 Winter Simulation Conference. Institute of Electrical and Electronics Engineers, Inc.
Balandat, M., Karrer, B., Jiang, D., Daulton, S., Letham, B., Wilson, A. G., and Bakshy, E. (2020). Botorch: A framework for efficient Monte-Carlo Bayesian optimization. In Larochelle, H.,, Ranzato, M., Hadsell, R., Balcan, M. F., and Lin, H., editors, Advances in Neural Information Processing Systems, volume 33, pages 21524-21538. Curran Associates, Inc.
Balasubramanian, K., Ghadimi, S., and Nguyen, A. (2020). Stochastic multi-level composition opti- mization algorithms with level-independent convergence rates. arXiv preprint arXiv:2008. 10526.
Bect, J.. Bachoc, F., and Ginsbourger, D. (2019). A supermartingale approach to Gaussian process based sequential design of experiments. Bernoulli, 25(4A)):883-2911.
Benaouali, A. and Kachel, S. (2019). Multidisciplinary design optimization of aircraft wing using commercial software integration. Aerospace Science and Technology. 92:766-776.
Bull, A. D. (2011). Convergence rates of efficient global optimization algorithms. Journal of Machine Learning Research, 12(88)2289--2904
Byrd, R. H., Lu, P., Nocedal, J., and Zhu, C. (1995). A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16(5):1190-125)
Cakmak, S., Astudillo Marban, R., Frazier, P., and Zhou, E. (2020). Bayesian optimization of
risk measures. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M. F,, and Lin, H., editors, Advances in Neural Information Processing Systems, volume 33, pages 20130-2014) Curran Associates, Inc.
Charisopoulos, V., Chen, Y., Davis, D., Díaz, M., Ding, L., and Drusvyatskiy, D. (2021). Low- rank matrix recovery with composite optimization: Good conditioning and rapid convergence. Foundations of Computational Mathematics, 1555):150-1
Cramer, E. J., Dennis, Jr, J. E., Frank, P. D., Lewis, R. M., and Shubin, G. R. (1994). Problem formulation for multidisciplinary optimization. SIAM Journal on Optimization, 4(4):754-776.
Damianou, and Lawrenee. N. D. (2013). Deep Gaussian processes. Carvalho, C. M. and Ravikumar, P., editors, Proceedings of the 16th International Conference on Artificial Intelligence and Statistics, volume of Proceedings of Machine Learning Research, pages 207-215. PMLR.
Drusvyatskiy, D. and Paquette, C. (2019). Efficiency of minimizing compositions of convex functions and smooth maps. Mathematical Programming, 178(1)50035558
Frazier, P. I. (2018). A tutorial on Bayesian optimization. arXiv preprint 107.0221.
Fu, M. C. and Henderson, S. G. (2017). History of seeking better solutions, AKA simulation optimization. In Chan, W. K. V., D' Ambrogi,, A., Zacharewicz, G., N. Mustafee, G. W., and Page, E., editors, Proceedings of the 2017 Winter Simulation Conference, pages 131-157. Institute of Electrical and Electronics Engineers, Inc.
Garnett, G. P. (2002). An introduction to mathematical models in sexually transmitted disease epidemiology. Sexually transmitted infections, 78(1):7-12.
11

---

Ghasemi, A., Heavey, C., and Laipple, G. (2018). A review of simulation simulation-optimzaatii optimization methods with upplications semiconductor operational problems. In Rabe, M., Juan, A. A.. Mustafee, N., Skoogh, A., Jain, S., and Johansson, B., editors, Proceedings of the 2018 Winter Simulation Conference, pages 3672-3683. Institute of Electrical and Flectronics Engineers, Inc.
Hansen, N. (2016). The CMA evolution strategy: A tutorial. arXiv preprini urXiv: 1604.00772.
-Simpic J. R., Gregoire, J. M., and Kusne, A. G. (2016). Perspective: Composition- ptructure-property mapping in high-throughuut experiments: Turning data into knowledge. APL Materials, 455):053211.
Jamil, M. and Yang, X.-S. (2013). literature survey of benchmark functions for global optimization problems. arXiv preprint arXiv:1308. 1308.4008.
Jones, D. R., Schonlau, M., and Welch, W. J. (1998). Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13(4):455-492
Kandasamy, K., Dasarathy, G., Schneider, J., and Póczos, B. (2017). Multi-fidelity Bayesian optimisation with continuous approximations. In Precup. D. and Teh, Y. W., editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 1799-1808. PMLR.
Kleywegt, A. J., Shapiro, A., and Homem-de Mello, T. (2002). The sample average approximation method for stochastic discrete optimization. SIAM Journal on Optimization, 12(2):479-502.
Li, H., Zhao, G., Johari, R., and Weintraub, G. Y. (2021). Interference, bias, and variance in two-sided marketplace experimentation: Guidance for platforms. arXiv preprint arXiv:2 /2222.
Lizotte, J.J., Wang. T., Bowling, M. H., and Schuurmans, D. (2007). Automatic gait optimization with Gaussian process regression. In Proceedings of the 20th International Joint Conference on Artificial Intelligence, pages 944-949.
Marque-Pucheu, S., Perrin, G., and Garnier, J. (2019). Efficient sequential experimental design for surrogate modeling of nested codes. ESAIM: PS, 23:245-270.
Močkus, J. (1975). On Bayesian methods for seeking the extremum. In Marchuk, G. I,, editor, Optimization Techniques IFIP Technical Conference, volume 27 Lecture Notes in Computer Science, pages 400-404. Springer- Verlag.
Osorio, C. and Bierlaire, M. (2013). A simulation- optimization framework for urban trans- portation problems. Operations Research, 61(6):1333-1344
Plappert, M., Andrychowicz, M., Ray, A., McGrew, B., Baker, B., Powell, G., Schneider, J., Tobin, J., Chociej, M., Welinder, P., et al. (2018). Multi-goal reinforcement learning: Challenging robotics environments and request for research. arXiv preprint arXiv: 1802.09464.
Puterman, M. L. (1990). Markov decision processes. In Stochastic Models, volume 2 of Handbooks in Operations Research and Management Science, pages 331 331-434. Elsevier.
Rasmussen, C. E. and Williams, C. K. (2006). Gaussian processes for machine learning. MIT Press.
Ryzhov, I. O. (2016). On the convergence rates of expected improvement methods. Operations Research, 6-55)11515-
Sadowski, M. and Jones, D. (2009). The sequence-structurg relationship and protein function prediction. Current Opinion in Structural Biology. 1933)3357-36..
Sandberg, I. (1978), On the mathematical foundations of compartmental analysis in biology, medicine, and ecology. IEEE transactions on Circuits and Systems, 2575):273-
Sekhon, B. and Saluja, V. (2. Biosimilars: An overview. Biosimilars, 1(1):1-11.
Shahriari, B., Swersky, K., Wang, Z.,, Adams, R. P., and De Freitas, N. (2016). Taking the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, (04(1):148-175.
12

---

Shapiro, A. (2003). On a class of nonsmooth composite functions. Mathematics of Operations Research, 2874
Snoek, J., Larochelle, H., and Adams, R. (2012). Practical Bayesian optimization of machine learn- ing algorithms. In Pereira, F., Burges, C. J. C., Bottou, L., and Weinberger, K. Q. editors, Advances in Neural Information Processing Systems, volume 25, pages 2951-2959. Curran Associates, Inc.
Sutton, R. S. and Barto, A. G. (2018). Reinforcement learning: An introduction. MIT press.
Swersky, K., Snoek, J., and Adams, R. P. (2013). Multi-task Bayesian optimization. In Burges, C. C., Bottou, L., Welling, M., Ghahramani, Z., and Weinberger, K. Q. editors, Advances in Neural Information Processing Systems, volume 26, pages 2004-2012. Curran Associates, Inc.
Toscano-Palmerin, S. and Frazier, P. I. (2018). Bayesian optimization with expensive integrands. arXiv preprini arXiv: 18800308666
Uhrenholt, A. K. and Jensen, B. S. (2019). Efficient Bayesian optimization for target vector estimation. In Chaudhuri. K. and Sugiyama, M., editors, Proceedings of the 22nd International Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings of Machine Learning Research, pages 2661-2670. PMLR.
Vazquez, E. and Bect, (2010). Convergence properties of the expected improvement algorithm with fixed mean and covariance functions. Journal of Statistical Plamning and Inference, 140( 3):3088- 3095.
Wang. X., Gong, X., Geng, N., Jiang, Z., and Zhou, L. (2020). Metamodel -based simulation optimisation for bed allocation. International Journal of Production Research, 58(20):6315 6335.
Wang, Z. and Jegelka, S. (2017). Max-value entropy search for efficient Bayesian optimization. In
Precup, D. and Teh, Y. W., editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 3627-3635. PMLR.
Williams, B. J., Santner, T. and Notz, W. I. (2000). Sequential design of computer experiments to minimize integrated response functions. Statistica Sinica, 10(4): 133-1152.
Wilson, J., Hutter, F., and Deisenroth, M. (2018). Maximizing acquisition functions for bayesian optimization. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 31, pages 9884-9895. Curran Associates, Inc.
Wu, J. and Frazier, (2016). The parallel knowledge gradient method for batch Bayesian optimization. In Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 29, pages 3134-3142, Curran Associates, Inc.
Wu, J., Dascano-Palmering S., Frazier, I., and Wilson, A. G. (2019). Practical multi-fidelity Bayesian optimization for hyperparameter tuning. In Proceedings of the 35th Conference Uncertainty in Artificial Intelligence.
Zhilinskas, A. (1975). Single-step Bayesian search method for an extremum of functions of a single variable. Cybernetics, 1 11(1):160-14
13


Reviews:
[{"review_id": "ITLWDHGXzY", "reviewer": "Reviewer_LhdJ", "summary": "This paper presents an extension of Bayesian optimization with Gaussian processes and expected improvement for target functions that are the results of composite functions from a DAG model. The paper introduces an extension of the expected improvement for the specific case which does not allow for a closed form solution. Instead the authors approximated the acquisition function with a sample approximation and optimized it using sample average approximation. The paper also includes some theoretical results in terms of consistency (although not convergence) and extensive results.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper extends on the work of Astudillo and Frazier, 2019 for general compositions of functions. In fact, even the structure of the paper is similar which simplifies the reading process if you are familiar with the previous paper. However, I think that this paper abuses the supplementary material. The paper should be self contained, but the text includes 14 times the word “supplement”. I understand that space constraints are an issue, especially after extensions from previous revisions. For example, in my opinion it would be better to fully describe and explain a couple of experiments in the main text, even including Figure 2 as is, instead of half describing all of them.\n\nThe extension and the method proposed seem original and, as pointed out in the related work, the problem can be motivated in many areas of science and engineering. I think the NeurIPS community might be quite interested in the RL setup.\n\nThe whole point of the statistical model and the SAA optimization is the result of the posterior distribution not being a Gaussian in the general case. However, the samples of any output h are drawn from Gaussian distributions, including hK = g. Furthermore, the authors recommend using SAA with deterministic samples, because the evaluations of the acquisition function are noisy. However, isn’t that the result of using a naive MC strategy? Having a DAG, wouldn’t be better to use belief propagation or exploit the structure with some MCMC or particle filter and have a more accurate representation of the acquisition function? I can see how BP can be much more expensive than the suggested approach, but this method is already quite expensive compared to regular BO. Furthermore, it seems that the use of SAA is not clearly motivated, being the default method for MC-based acquisitions in BoTorch.\n\nThe results are quite interesting. Clearly, the method reaches SOTA, even compared to simpler versions, such as EI-CF. However, the authors mix absolute value with log(regret). For instance, they claim that the improvement in the robot problem is quite large. However, any real robot application can be considered “solved” with log(regret) = -2 or -3. From a pure numerical/benchmark point of view, the result is still interesting, but then I would like to have a more realistic RL setup. For example, the suggested setup or learning a sequence of waypoints as a policy has been previously done in BO for robot navigation:\n\nRuben Martinez-Cantin, Nando de Freitas, Eric Brochu, Jose Castellanos and Arnaud Doucet (2009) A Bayesian Exploration-Exploitation Approach for Optimal Online Sensing and Planning with a Visually Guided Mobile Robot. Autonomous Robots - Special Issue on Robot Learning, Part B, 27(3):93-103.\n\nSmall comments:\n\n-Given that KG is sometimes superior to EI, would it be possible to have KG-FN or is there a limitation that could prevent it?\n\n-If the robot pushing and epidemic model are in the perfect setup for EI-CF, how would you explain the improvement of EI-FN?\n\n----\nNotes after author response: The authors have addressed most of my concerns. I'm glad that the new experiment worked so well. It clearly shows the advantages of the proposed method and I have updated my score as a result. \n\nAs a final suggestion, one of the advantages of the method is clearly computational cost. Thus, I would address that with some results or comments regarding that. The fact that \"vanilla\" KG is x4 more expensive than EI-FN is remarkable, as one would expect EI-FN to be the more expensive given the number of GPs to manage.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "8"}, "extra_text": {}}, {"review_id": "D7Pw8yzBtQ0", "reviewer": "Reviewer_HdZA", "summary": "This paper extends methods for Gaussian process \"Bayesian Optimization\" (BO)  of vectors of hyperparameters to structural networks of function, to include hyperparameters from functions of intermediate stages. Unlike a conventional approach, this method considers the output of intermediate stages.  The hypothesis addressed by the paper that including additional information from intermediate stages is a reasonable one, however it raises a challenge since despite the intermediate functions each being modelled as a GP, the final output g is not necessarily Gaussian.  To accommodate this, they apply sample average approximation to sample from the posterior on the final output.\n\nThe authors apply this method in the context of several application areas, having run experiments on 8 problems with results showing substantially faster convergence than competing methods.  \n\nThis paper offers a convincing empirical study together with theoretical convergence results. This is an interesting paper with useful learning to the general community.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper's approach is to provide an approximation to the posterior for the acceptance function, EI-FN, using a sampling method based on the 'so-called re-parameterization trick.\"  In Proposition 1 they show a convergence proof for the function.   The optimization efficiency they claim is demonstrated in the set of empirical problems on both simulated and actual data. \n\nAs for the evaluation criterion used, in the calculation of the approximated EI-FN function, the method generates M samples (even from the posterior distribution without actually evaluating the network) that can potentially reduce the variance of this evaluation. It appears that the posterior distribution is still estimated only using the n samples. When counting the number of iterations,  in Figure 2, is the x-axis  only counting n as the iteration number?  One wonders if the the total number of function evaluations would be a fairer basis for comparison , given the implication of this paper is that this “synthetic” M samples are the key to improve the convergence rate. This leads to a speculation if there exists other ways of generating the samples that can possibly lead to even better convergence rate (either based on posterior distribution or not).\n\nAs claimed, basic information-theoretic notions lead one to believe that the intermediate values used in the EI-FN function should lead to more efficient sampling estimates, hence faster optimization.  This property, although demonstrated in the empirical results, is not supported, and deserves further study since the analysis provided only hints in broad terms why it comes about. Convergence rates alone are only an indirect explanation for the observed efficiency.  Perhaps there are diagnostics derived from comparing conventional acquisition functions with the approximations generated?  \n\nTo speculate further, the example network structures resemble probabilistic graphical models, so one could imagine an inferential approach to compute distributions over the x's that would enable estimation of how much improvement is due to an x's additional information.   This is not to criticize the work to date, but it's surprising that thoughts about this did not get a mention. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors are cognizant that this technique can be used for ends that are both socially valuable or destructive, and state as much.  In it's current nascent form there is no reason to presume it's use would encourage abusive uses. ", "review_text": "", "extra_scores": {"ethical_concerns": "No ethical concerns arise in this current work. ", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3 hours"}, "extra_text": {}}, {"review_id": "55uzGPFC0w", "reviewer": "Reviewer_vzc8", "summary": "This paper focuses on Bayesian optimization of function networks, where the objective function has a known structure that is a composition of individual (node) functions that take the output of the previous function as input (along with potentially other parameters). Observations at each node function are obtained at each objective function evaluation. Such objective functions arise in a variety of applications including materials design and physical processes. Exploiting knowledge of the network structure can result in a lower effective input dimensionality of the model of each node function, depending on the network structure. By modeling each node function as distinct GP, the authors demonstrate improved sample efficiency using EI compared to optimizing the global objective function with a single GP or using compositional objectives.", "questions": "", "limitations": "", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths:\n-\tThe idea that accounting for this structure is beneficial is simple and intuitive\n- The empirical performance is quite compelling\n\nWeaknesses:\n-\tThere is not much technical novelty. Given the distinct GPs modeling the function network, the acquisition function and sampling procedure are not novel\n-\tThe theoretical guarantee is pretty weak (random search is asymptotically optimal).\n  - The discussion of not requiring dense coverage to prove the method is asymptotically consistent is interesting, but the utility of proposition 2 is not clear because although dense coverage is a consideration for proving consistency, it is not really a practical reality in sample-efficient optimization—typically BO would not have dense coverage.\n\nQuestions/comments:\n-\tThere is no discussion of observation noise, which is a practical concern in many of the real world use cases mentioned in the paper. The approach of using GPs to model nodes in function network can naturally handle noisy observations, so only the acquisition function would need to be adjusted to account for noisy observations since the best objective value would be unknown. I expect that the empirical performance would remain the same (e.g. using Noisy EI from Letham et al. 2019), but the computation would be much more expensive. It would be good to discuss and demonstrated performance under noisy observations.\n-\tHow does the number of MC samples affect performance, empirically? How does the network structure affect this?\n-\tIt would be interesting to see a head-to-head comparison with deep GPs. How different are the runtimes (including inference times) and empirical performances? \n- Since the core contribution is modeling each node in the function network with a distinct GP, it would be good to see more evaluation of the function network model's predictive performance compared to a alternative modeling choices (e.g. individual models with a compositional objective, vanilla global gp, deep gp)\n\nGrammar:\n-\tL238 “out method” -> “our method”\n-\tL335 “structurre” -> “structure”\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The discussion of the work's limitations is quite thorough, and it proposes interesting directions for future work. The authors have addressed potential negative societal impacts.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "-QcAYAZYDRy", "reviewer": "Reviewer_r59T", "summary": "This work investigates the problem of Bayesian optimization of the output of a network of functions. The network is organized as a directed acyclic graph, where each function takes as input the output of its parent nodes as well as (part of ) the input. Contrary to the standard setting that only observes the final output, this work assumes that intermediate nodes can also produce an output (which can be stochastic or deterministic). The authors model the nodes of the network using Gaussian processes, and employ an acquisition function that corresponds to the expected improvement with respect to the implied posterior on the final output under a proper statistical model. The authors optimize the acquisition function using sample average approximation. The authors show that their approach satisfies asymptotic convergence. Various experiments on synthetic and real problems demonstrate that the improved model and acquisition function result in visibly improved solution quality, at the expense of the runtime cost.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths\n- This work introduces a new setting for Bayesian Optimization (BO) that is highly relevant in practice, e.g., in engineering applications. This work seems to be the first to address this setting. For this reason, this work is fairly significant.\n- On the theory front, the work includes a convergence proof. Furthermore, it investigates the interesting property that the acquisition function for function networks does not need to measure the optimization domain densely. This suggests a possible mechanism via which the new acquisition function can lead to more efficient BO.\n- The paper is well written and easy to follow.\n- The experimental section is well written and demonstrates that the new model and acquisition function achieve superior result in terms of the objective value. The problems considered in the experimental study are diverse and quite realistic.\n\nConcerns\n- Using GPs to model the various nodes in the function network, and optimizing the acquisition function via sample average approximation are not particularly novel for that setting. Even though the problem statement in this work is novel, the authors have chosen a rather straightforward approach to tackle it. That being said, they do show consistency.\n- Convergence rates are not discussed, only consistency.\n- The method is quite time-consuming compared to simpler competitors, as illustrated in the runtimes. This problem can be exacerbated for large function networks with several nodes.\n\nAdditional questions to authors\n- How does the computational complexity depend on the number of nodes in the function network? Experimentally, function networks of up to how many nodes have the authors been able to handle?\n- The framework can deal with both deterministic and non-deterministic nodes. What about nodes whose output is hidden? Can the framework handle such cases for an arbitrary node (assuming of course the DAG structure for the function network)?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations and potential societal impact are discussion in the Conclusion.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper addresses Bayesian optimization which leverages the intermediate outputs as well as final outputs when they are available by evaluations in the function network. The problem setting is novel, extending the previous work [Astudillo, R. and Frazier, P., 2019]. A cascade of GPs is used as a surrogate model, leading to non-Gaussian process. A sample average approximation approach is used to optimize the EI. Strong empirical results are provided to demonstrate that the proposed method indeed achieves superior performance. The strength of this paper is in the novel problem setting, its usefulness in engineering applications mainly for practitioners. The downside is in the limited novelty in the proposed solution, since it is a direct employment of GPs and SAA. During the committee discussion period, I had a few communications with reviewers. The problem setting, which appears often in engineering applications, is valuable for further studying.  Thus, even if the proposed solution is not novel, this work studies an important and valuable problem, and could serve as a starting point for new works that could develop more sophisticated algorithms or deeper theory for the DAG-dependency setting.


Decision: Accept (Poster)