# Collaborative and Distributed Bayesian Optimization via Consensus:
Showcasing the Power of Collaboration for Optimal Design

Xubo Yue  Northeastern University, Boston  Yang Liu  Virginia Tech, Blacksburg
Albert S. Berahas  University of Michigan, Ann Arbor  Blake N. Johnson
Virginia Tech, Blacksburg  Raed Al Kontar  University of Michigan, Ann Arbor
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Optimal design is a critical yet challenging task within many applications.
This challenge arises from the need for extensive trial and error, often done
through simulations or running field experiments. Fortunately, sequential
optimal design, also referred to as Bayesian optimization when using
surrogates with a Bayesian flavor, has played a key role in accelerating the
design process through efficient sequential sampling strategies. However, a
key opportunity exists nowadays. The increased connectivity of edge devices
sets forth a new collaborative paradigm for Bayesian optimization. A paradigm
whereby different clients collaboratively borrow strength from each other by
effectively distributing their experimentation efforts to improve and fast-
track their optimal design process. To this end, we bring the notion of
consensus to Bayesian optimization, where clients agree (i.e., reach a
consensus) on their next-to-sample designs. Our approach provides a generic
and flexible framework that can incorporate different collaboration
mechanisms. In lieu of this, we propose transitional collaborative mechanisms
where clients initially rely more on each other to maneuver through the early
stages with scant data, then, at the late stages, focus on their own
objectives to get client-specific solutions. Theoretically, we show the sub-
linear growth in regret for our proposed framework. Empirically, through
simulated datasets and a real-world collaborative sensor design experiment, we
show that our framework can effectively accelerate and improve the optimal
design process and benefit all participants.

##  1 Introduction

The success of many real-world applications critically depends on trial &
error. Often the goal is to manipulate a set of variables, called designs, to
achieve a desired outcome or system response. For instance, material
scientists perform time-consuming and expensive experiments to determine
optimal compositions (Zhang et al., 2023a, ) that produce a material with
desired properties. Similarly, additive manufacturers must calibrate many
design parameters, including laser power, beam diameter, and hatching pattern
through trial & error (Buchanan and Gardner,, 2019) so that their product
matches its intended shape. Those studies, whether through experimentation or
simulations like finite element models, consume resources (e.g., time and
budget) that can significantly limit progress.

Fortunately, sequential design has played a key role in accelerating the
design process (Le Gratiet and Cannamela,, 2015; Kusne et al.,, 2020; Huang et
al.,, 2021). In the engineering and statistical literature, sequential design
can be categorized into two types (Ezzat et al.,, 2018): (i) sequential design
for exploration, that aims to explore a response surface, (ii) sequential
design for optimization, that aims to find designs that optimize a target
response. Our work focuses on sequential design for optimization, known as
sequential optimal design. This field has also been coined more recently as
Bayesian optimization (BO) when using surrogates with a Bayesian flavor (e.g.,
Gaussian processes - $\mathcal{GP}$s). Hereon, we will use the term Bayesian
optimization since, without loss of generality, we exploit $\mathcal{GP}$s in
the process.

Rather than running exhaustive brute-force experiments over a dense grid, BO
employs a sequential strategy to conduct experiments and observe new samples.
It first fits a surrogate that estimates the design-response relationship from
existing data. Afterward, a utility, also known as acquisition, function
(Gramacy and Lee,, 2009; Overstall and Woods,, 2017; Lee et al.,, 2018) is
defined to hint at the benefits of sampling/experimenting at new design
points. Based on this utility metric, the next-to-sample design point is
selected. The procedure is sequentially iterated over several rounds till
budget constraints or exit conditions are met. Needless to say, sequential
optimal design/BO has been extensively studied (Antony,, 2014; Shahriari et
al.,, 2015; Kleijnen,, 2018) and has found success in a wide variety of
disciplines across physics, chemistry, mechanical engineering (Ekström et
al.,, 2019; Kusne et al.,, 2020; Zhang et al.,, 2020; Rahman et al.,, 2021;
Wang and Dowling,, 2022), amongst many others.

In lieu of the aforementioned successes, this work aims to bring BO to a
collaborative paradigm. The main question we ask is: How can multiple clients
collaborate to improve and fast-track their design processes? With today’s
advances in computation and communication power at edge devices (Kontar et
al.,, 2021), it has become more plausible for potentially dispersed clients to
share information, distribute trial & error efforts, and fast-track the design
process so that all participants gain benefit. Here, collaborating clients can
be scientists, robots, multiple finite element simulations, etc.

However, to enable collaboration, key challenges exist. First, and most
importantly, is how to distribute the sequential optimization process. While
BO has been extensively studied in the past decades, few literature exists on
collaborative BO. The second challenge is heterogeneity. Despite trying to
optimize similar processes, clients may operate under different external
factors and conditions. As such, retaining client-specific optimal solutions
is of importance. The third challenge is privacy. To encourage participation,
a collaborative process should refrain from sharing client-specific outcomes.

In an effort to address the opportunity and stated challenges above, we
summarize our contributions below:

  * •

Bayesian optimization via Consensus: We bring the notion of “consensus” to BO
where clients perform experiments locally and agree (i.e., reach a consensus)
on their individualized next-to-sample designs. Our approach provides a
generic and flexible framework that can incorporate different collaboration
mechanisms. It hinges upon a consensus matrix that evolves with iterations and
flexibly determines the next-to-sample designs across all clients. Some
interesting by-products of our framework are:

    1. 1.

Amenability to full decentralization that respects communication restrictions
and allows clients to perform the collaborative process without the
orchestration of a central entity.

    2. 2.

Ability to only share designs and not client-specific outcomes to preserve
some privacy.

    3. 3.

Amenability to any utility/acquisition function of choice.

  * •

Regret minimization: We support our proposed collaborative framework with a
theoretical justification. We show that, under some mild conditions, with a
high probability, the cumulative regret for each client has a sublinear growth
rate. This implies that our collaborative algorithm can bring clients within
regions of optimal designs.

  * •

The power of collaboration for optimal sensor design: We perform an experiment
for accelerated optimization of sensor design through collaborative finite
element analysis (FEA) workflows. The experiment showcases the ability of
collaboration to significantly outperform its non-collaborative counterpart.

The remainder of this paper is organized as follows. In Sec. 2, we present a
motivating example that drives this work. We then briefly discuss related
literature in Sec. 3. In Sec. 4, we formulate the problem mathematically and
provide our proposed collaborative BO framework. Sec. 5 and Sec. 6 showcase
the benefits of our approach through several simulation studies and an
experiment on accelerated sensor design. We conclude our paper with a brief
discussion about open problems in this relatively new field in Sec. 7.

##  2 Motivation

Our work is motivated by optimal sensor design. Biosensors are critical
bioanalytical technologies that enable the selective detection of target
analytes and have broad applications ranging from medical diagnostics, health
monitoring, food and water safety, and environmental monitoring (Cesewski and
Johnson,, 2020). A biosensor is defined as a device that is based on an
integrated biorecognition element and transducer (Thévenot et al.,, 2001).
There are two main categories of biosensors: (1) device-based biosensors, and
(2) methods-based biosensors (e.g., nanobiosensors) (Johnson and Mutharasan,,
2014). Device-based biosensors are often integrated with microfluidics and
exhibit form factors that can be physically integrated with and removed from
processes (e.g., thin-film, dip-stick) and preserve the characteristics of the
sample to be analyzed. Alternatively, methods-based biosensors based on
solutions and suspensions of functionalized particles provide detection by
mixing with a sample, and thus, are relatively destructive with respect to the
sample. Thus, device-based biosensors are commonly used in a “flow-and-
measure” or “dip-and-measure” format. In contrast, methods-based biosensors,
such as nanobiosensors, are typically used in a “mix-and-measure” format.

While biosensors have been created for sensitive and selective detection of
many target analytes ranging from small molecules, proteins, nucleic acids,
cells, biomarkers, and pathogens, there remains a need to further optimize
biosensor design (e.g., form factor, functionalization) and performance (e.g.,
sensitivity, measurement confidence, and speed) to meet the constraints and
requirements of industrial and commercial applications (Carpenter et al.,,
2018). Current thrust areas in biosensor design and performance optimization
can be categorized as driven by experimentation, simulation, data analytics,
or combinations thereof (Selmi et al.,, 2017; Polatoğlu et al.,, 2020; Koveal
et al.,, 2022; Zhang et al., 2023b, ; Zhang et al.,, 2024). For example, high-
throughput experimentation, sensor arrays, FEA, and data-driven biosensing
have been leveraged to improve the understanding, design, and performance of
biosensors, particularly device-based biosensors whose design, utility, and
performance are often linked to the characteristics of experimental
measurement formats, such as microfluidic channel design, flow field
parameters (e.g., flow rate), and other parameters of the measurement format
(e.g., sample injection time). However, it remains a challenge to optimize
biosensor performance, particularly in maximizing the amount of target analyte
detected, given the high-dimensionality of the design space associated with
biosensor design, functionalization (e.g., concentration of immobilized
biorecognition elements), and measurement format parameters. In particular,
there remains a need to establish closed-loop self-driving workflows for
engineering high-performance biosensors, such as by optimizing the biosensor
design, functionalization protocols, and measurement formats that synergize
with experimentation, simulation, and machine learning. Given the current
limitation, our central hypothesis is:

Hypothesis: Collaboration can accelerate the pace of optimal sensor design and
yield optimal design with minimal resource expenditures.

Fig. 1 provides a microcosm of our collaborative solution, which we have
tested in Sec. 6. Our experiment features multiple collaborating agents (FEA
simulations) that perform biosensor design. Then they will coordinate to
decide on their next simulations. Our goal is to rapidly discover the
biosensor design and measurement format parameters that find the the maximum
amount of captured target analyte.

Figure 1: Illustration of closed-loop biosensor and measurement format design
optimization driven by consensus BO-driven simulation.

##  3 Related Work

Optimal design has a longstanding history. Traditionally, the major trend for
optimal design was fixed-sample designs, where the best design is chosen from
a pre-determined set. Perhaps the earliest departures from fixed-sample
designs date back to a few years following World War II, where Friedman and
Savage, (1947); Box and Wilson, (1951); Robbins and Monro, (1951); Bellman,
(1956); Chernoff, (1959) proposed schemes to sequentially decide on next-to-
sample designs to eventually reach the optimal design. Amongst them were the
experiments done by E.P Box and co-workers Box and Wilson, (1951) at the
Imperial Chemical Industries, aimed at improving the yield of chemical
processes. Afterward, two seminal papers from Kushner, (1964); Sacks and
Ylvisaker, (1970) proposed deciding on the next-to-sample design points
through optimizing utility functions driven by surrogates that can quantify
uncertainty. This led to an era of development in BO that remains very active
till this day.

Recent advances in BO include but are not limited to: (i) Alternative
surrogate models: instead of using $\mathcal{GP}$ surrogates, Marmin and
Filippone, (2022); Ming et al., (2022); Sauer et al., (2022) suggest using
Deep $\mathcal{GP}$s as an alternative. This credits to many desirable
advantages of Deep $\mathcal{GP}$s, such as flexibility to non-stationarity
and robustness in handling abrupt regime changes in the training data. Besides
deep $\mathcal{GP}$s, Bayesian neural networks (Snoek et al.,, 2015;
Springenberg et al.,, 2016) have also been widely adopted as surrogates. (ii)
Multi-objective BO: where the goal is to simultaneously optimize multiple,
often competing, responses (Konomi et al.,, 2014; Svendsen et al.,, 2020).
Here surrogates that simultaneously learn the multiple outputs, such as the
Multi-output $\mathcal{GP}$ (Kontar et al.,, 2018), are often exploited and
developed. (iii) Multi-fidelity BO: when data is collected across multiple
fidelities, it becomes crucial to choose the fidelity to use when running an
experiment. Along this line, Le Gratiet and Cannamela, (2015); He et al.,
(2017); Stroh et al., (2022) have developed various sequential strategies to
tackle this challenge. (iv) New utility functions: developing new utility
functions remains one of the most active areas in optimal design. Along this
line, recent literature has investigated look-ahead utilities that chose a
design based on its utility over a rolling future horizon (Jiang et al.,,
2020; Lee et al.,, 2020; Yue and Kontar,, 2020).

Though BO has caught major attention over the past century, to our knowledge,
little literature exists on collaborative BO. Perhaps the two closest fields
are batch sequential design (or batch BO) and federated BO. In batch BO,
multiple designs are chosen from a common surrogate and utility. These designs
are then distributed across multiple compute nodes to be evaluated in parallel
(Azimi et al.,, 2010, 2012; Wu and Frazier,, 2016; Duan et al.,, 2017; Hunt,,
2020). Unfortunately, batch BO is designed for parallel computing, where it
focuses on optimizing a single objective using a single surrogate learned from
a centralized dataset. Therefore, it cannot handle cases where data comes from
diverse and potentially heterogeneous sources, nor does it preserve privacy.
All data needs to be agglomerated in one place to learn the surrogate and
optimize the utility. In a similar fashion, federated BO (Dai et al.,, 2020,
2021) tries to distribute experiments decided from a single objective.
However, they do so while preserving privacy. For instance, Dai et al., (2020)
share function realizations of the posterior belief using random Fourier
features to learn the common surrogate. While preserving privacy, federated BO
tries to circumvent its inherent homogeneity assumptions by reducing
collaboration after some time, yet such a method is limited to a Thomson-
sampling utility.

##  4 The Collaborative BO Framework

###  4.1 Setting the Stage

We start by describing our problem setting and introducing notation. Assume
there are $K\geq 2$ clients, and each client has a budget of $T$ experiments
across $T$ iterations. Denote by $t\in\\{0,1\ldots,T-1\\}$ the iteration
index. Clients can communicate with each other either via a central
orchestrator or direct communication (See Sec. 4.2).

Each client $k\in[K]\coloneqq\\{1,\ldots K\\}$ has an initial dataset
$\mathcal{D}^{(0)}_{k}=\\{\bm{X}^{(0)}_{k},\bm{y}^{(0)}_{k}\\}$ with
$N^{(0)}_{k}$ observations, where
$\bm{X}^{(0)}_{k}=\\{\bm{x}_{k,1},\ldots,\bm{x}_{k,N^{(0)}_{k}}\\}$ is a
$D\times N^{(0)}_{k}$ design matrix that contains the initials designs
$\bm{x}_{k,\cdot}\in\mathbb{R}^{D}$ and
$\bm{y}^{(0)}_{k}=(y_{k,1},\ldots,y_{k,N^{(0)}_{k}})^{\intercal}$ is an
$N^{(0)}_{k}\times 1$ vector that contains the corresponding observed
responses. The goal of each client is to find a set of client-specific optimal
designs

$\displaystyle\bm{x}^{*}_{k}=\operatorname*{argmax}_{\bm{x}\in\mathcal{X}\subseteq\mathbb{R}^{D}}f_{k}(\bm{x})\,,$

where $\mathcal{X}$ is a subset in $\mathbb{R}^{D}$, and
$f_{k}:\mathbb{R}^{D}\rightarrow\mathbb{R}$ is the true unknown design
function each client aims to optimize. Clearly, using first or second-order
optimization algorithms is not feasible since $f_{k}$ is a black-box. To
observe $f_{k}(\cdot)$ at a new design point $\bm{x}_{k}^{\text{new}}$ one
needs to run an experiment and observe $y_{k}(\bm{x}_{k}^{\text{new}})$, that
is possibly a noisy representation of $f_{k}(\bm{x}_{k}^{\text{new}})$;
$y_{k}=f_{k}+\epsilon_{k}$, where $\epsilon_{k}$ is an additive noise.
Therefore, the goal is to carefully decide on the next-to-sample design
$\bm{x}_{k}^{\text{new}}$ so that an optimal design is reached with the fewest
experiments possible.

To do so, at any time $t$, BO resorts to a utility function $U$,
$U(y_{k}(\bm{x});\mathcal{D}^{(t)}_{k}):\mathbb{R}^{D}\to\mathbb{R}$, that
quantifies the benefits gained if one were to conduct an experiment at a new
design, $\bm{x}_{k}^{(t)\text{new}}$ (Gramacy and Lee,, 2009; Konomi et al.,,
2014; Amo-Salas et al.,, 2016; Overstall and Woods,, 2017; Lee et al.,, 2018).

In a non-collaborative environment, each client $k$ chooses the next-to-sample
design by maximizing their own utility. However, since the utility is
dependent on the response $y_{k}$, one cannot calculate the utility except at
previously observed designs. Here BO resorts to a surrogate $\hat{y}_{k}$ that
estimates $\hat{y}_{k}(\bm{x})$ for any $\bm{x}\in\mathbb{R}^{D}$ using the
dataset $\mathcal{D}_{k}^{(t)}$. Such surrogates are often probabilistic
(e.g., $\mathcal{GP}$ (Sacks et al.,, 1989; Currin et al.,, 1991)) and are
capable of providing a predictive distribution
$\mathbb{P}_{\hat{y}_{k}(\bm{x})|\mathcal{D}^{(t)}_{k}}$ over
$\hat{y}_{k}(\bm{x})$. Equipped by the surrogate, client $k$, now chooses the
next-to-sample design by maximizing the expected utility

$\displaystyle\bm{x}_{k}^{(t)\text{new}}=\operatorname*{argmax}_{\bm{x}}\mathbb{E}_{\mathbb{P}_{\hat{y}_{k}|\mathcal{D}^{(t)}_{k}}}\left[U(\hat{y}_{k}(\bm{x});\mathcal{D}^{(t)}_{k})\right].$
(1)

Hereon, for notation brevity, we write the expected utility in (1) as
$\mathbb{E}_{{\hat{y}_{k}|\mathcal{D}^{(t)}_{k}}}\left[U(\hat{y}_{k}(\bm{x}))\right]$.
It is also worth noting that expectation in (1) with respect to $\hat{y}_{k}$
is sometimes replaced with $\hat{f}_{k}$ where
$\hat{y}_{k}=\hat{f}_{k}+\epsilon_{k}$, depending on the type of utility
function. See Sec. 4.4 for some examples.

###  4.2 Collaborative BO via Consensus

Now in a collaborative framework, given an arbitrary $t$, we aim to allow
clients to collaboratively decide on
$\\{\bm{x}_{k}^{(t)\text{new}}\\}_{k=1}^{K}$. A natural idea is to maximize
utility across all clients. This translates to the following problem:

$\displaystyle\max_{\\{\bm{x}_{k}\\}_{k=1}^{K}}\mathbb{E}_{k}\left[\mathbb{E}_{\hat{y}_{k}|\mathcal{D}^{(t)}_{k}}\left[U(\hat{y}_{k}(\bm{x}_{k}))\right]\right]=\max_{\\{\bm{x}_{k}\\}_{k=1}^{K}}\sum_{k=1}^{K}p_{k}\left[\mathbb{E}_{\hat{y}_{k}|\mathcal{D}^{(t)}_{k}}\left[U(\hat{y}_{k}(\bm{x}_{k}))\right]\right],$
(2)

where $p_{k}$ is some weight coefficient for client $k$ with
$\sum_{k=1}^{K}p_{k}=1$. Without loss of generality, hereon we assume
$p_{k}=\frac{1}{K}$. As shown, the key difference of (1) from (2) is taking
the expectation $\mathbb{E}_{k}$ over all participants. However, in itself,
(2) does not allow entities to borrow strength from each other as it can be
fully decoupled across the $K$ entities.

In order to enable collaboration and distribute experiments across the $K$
clients, we bring the notion of “consensus” (Shi et al.,, 2017; Berahas et
al.,, 2018) to BO. In the context of BO, consensus allows clients to agree
(reach a consensus) on their individualized next-to-sample designs through a
consensus matrix. Specifically, we modify (2) to

$\displaystyle\bm{x}^{(t)}_{k}=[\bm{x}^{(t)}_{\mathcal{C}}]_{k}=\operatorname*{argmax}_{\bm{x}_{k}}\left[\mathbb{E}_{\hat{y}_{k}|\mathcal{D}^{(t)}_{k}}\left[U(\hat{y}_{k}(\bm{x}_{k}))\right]\right]$
$\displaystyle\quad\text{and}\quad\bm{x}_{k}^{(t)\text{new}}=\left[(\bm{W}^{(t)}\otimes\bm{I}_{D})\bm{x}^{(t)}_{\mathcal{C}}\right]_{k}\,,$
(3)

where $[\cdot]_{k}$ represents the $k^{th}$ block of a vector,
$\bm{x}^{(t)}_{\mathcal{C}}=[\bm{x}^{(t)\top}_{1},\bm{x}^{(t)\top}_{2},\cdots,\bm{x}^{(t)\top}_{K}]^{\top}$
is concatenation of the designs across all clients, $\bm{W}^{(t)}$ is a
consensus matrix of size $K\times K$,

$\displaystyle\bm{W}^{(t)}=\begin{bmatrix}w^{(t)}_{11}&w^{(t)}_{12}&\cdots&w^{(t)}_{1K}\\\
w^{(t)}_{21}&w^{(t)}_{22}&\cdots&w^{(t)}_{2K}\\\
\vdots&\vdots&\ddots&\vdots\\\
w^{(t)}_{K1}&w^{(t)}_{K2}&\cdots&w^{(t)}_{KK}\end{bmatrix},$

$\bm{I}_{D}$ is a $D\times D$ identity matrix, and $\otimes$ denotes the
Kronecker product operation, i.e., $\bm{W}^{(t)}\otimes\bm{I}_{D}$ results in
a matrix of size $DK\times DK$. The matrix $\bm{W}^{(t)}$ is a symmetric,
doubly stochastic matrix (i.e., $\sum_{k}w^{(t)}_{kj}=\sum_{j}w^{(t)}_{kj}=1$
for $j,k\in[K]$) with non-negative elements.

The new objective (3) has several interesting features. For the sake of
compactness, unless necessary, we drop the superscript $t$ in the subsequent
discussion. First, the formulation presented in (3) is indeed reminiscent of
recent optimization approaches coined as consensus optimization. To see this,
notice that a doubly stochastic matrix $\bm{W}$ has the property that
$(\bm{W}\otimes\bm{I}_{D})\bm{x}_{\mathcal{C}}=\bm{x}_{\mathcal{C}}$ if and
only if $\bm{x}_{k}=\bm{x}_{j}$ for all $k,j\in[K]$ (Nedic and Ozdaglar,,
2009). As such, if we enforce
$(\bm{W}\otimes\bm{I}_{D})\bm{x}_{\mathcal{C}}=\bm{x}_{\mathcal{C}}$ as a
constraint, we can solve

$\displaystyle\max_{\bm{x}_{k}}\frac{1}{K}\sum_{k=1}^{K}\left[\mathbb{E}_{\hat{y}_{k}|\mathcal{D}_{k}}\left[U(\hat{y}_{k}(\bm{x}_{k}))\right]\right]\quad\text{subject
to}\quad\bm{x}_{k}=\bm{x}_{j},\,\,\,\forall k,j\in[K].$ (4)

The equality constraint (often referred to as the consensus constraint) is
imposed to enforce that local copies at every client are equal. That said, our
formulation and setting have several distinguishing features that
differentiate it from consensus optimization. The differences pertain to the
goals of the two approaches and are related to the consensus matrix $\bm{W}$.
In consensus optimization, the goal is for all clients to eventually agree on
a common decision variable. However, in our collaborative BO paradigm, we do
not want all clients to make the same decisions. Rather, we want clients to
borrow strength from each other while at the same time allowing for
personalized (per client) solutions. Even when clients are homogeneous,
enforcing the constraint will significantly reduce exploration as everyone
runs an experiment at the same location. Therefore, we do not explicitly
enforce the consensus constraint. Instead, we allow $\bm{W}$ to play an
aggregation role that decides on the next-to-sample designs
$\bm{x}_{k}^{(t)\text{new}}$ given the current utility maximizers
$\bm{x}_{\mathcal{C}}^{(t)}$. More importantly, $\bm{W}^{(t)}$ is time-varying
and, in the limit, converges to the identity matrix (see below for more
details and examples of the consensus matrices used). This allows clients to
borrow strength in the initial stages of the optimization, yet, eventually,
make personalized decisions.

Second, the consensus matrix $\bm{W}$ is doubly stochastic. By the
Birkhoff–von Neumann theorem (Marshall et al.,, 1979), any doubly stochastic
matrix can be expressed as a convex combination of multiple permutation
matrices. Mathematically, there exists $L$ non-negative scalars
$\\{\eta_{l}\\}_{l=1}^{L}$ such that $\sum_{l=1}^{L}\eta_{l}=1$ and
$\bm{W}=\sum_{l=1}^{L}\eta_{l}\bm{P}_{l}$. For example,

$\displaystyle\bm{W}=\begin{bmatrix}0.2&0.3&0.5\\\ 0.6&0.2&0.2\\\
0.2&0.5&0.3\end{bmatrix}=0.2\begin{bmatrix}0&1&0\\\ 0&0&1\\\
1&0&0\end{bmatrix}+0.2\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&1\end{bmatrix}+0.1\begin{bmatrix}0&1&0\\\ 1&0&0\\\
0&0&1\end{bmatrix}+0.5\begin{bmatrix}0&0&1\\\ 1&0&0\\\ 0&1&0\end{bmatrix},$

where each permutation matrix $\bm{P}_{l}$ can be viewed as an allocator that
assigns one design solution to one client. In this example, we can see that

$\displaystyle(\bm{P}_{1}\otimes\bm{I}_{D})\bm{x}_{\mathcal{C}}=\left(\begin{bmatrix}0&1&0\\\
0&0&1\\\
1&0&0\end{bmatrix}\otimes\bm{I}_{D}\right)\begin{bmatrix}\bm{x}_{1}\\\
\bm{x}_{2}\\\ \bm{x}_{3}\end{bmatrix}=\begin{bmatrix}\bm{x}_{2}\\\
\bm{x}_{3}\\\ \bm{x}_{1}\end{bmatrix}.$

This implies that the first allocator assigns the solution from client 1 to
client 2, and so forth. This is an interesting phenomenon that is related to
batch BO, where a batch of candidate design points is selected, and then
multiple experiments run in parallel, each corresponding to one candidate
solution. Now, instead of sticking to one allocation strategy, the consensus
approach can be viewed as a weighted average of all possible assignments,
where the designed consensus matrix dictates the weights, hence the impact of
clients on each other.

Third, over the course of the optimization, the consensus step
$(\bm{W}\otimes\bm{I}_{D})\bm{x}_{\mathcal{C}}$ naturally yields $K$ solutions
to all clients (i.e., one for each) so that clients can explore and exploit
the solution space independently and in a distributed manner. For example,
suppose $D=1$ and consider a consensus matrix
$\bm{W}=\begin{bmatrix}0.7&0.3\\\ 0.3&0.7\end{bmatrix}$ and two local
solutions $\bm{x}_{1}=5,\bm{x}_{2}=7$. The consensus step will yield
$\bm{W}\bm{x}_{\mathcal{C}}=(5.6,6.4)^{\intercal}$ such that client 1 will
take solution $5.6$, and client 2 will take solution $6.4$.

Fourth, from the previous example, it can be seen that the consensus matrix
$\bm{W}$ controls how much one client will affect the design choices for other
clients. As a result, the matrix $\bm{W}$ adds a layer of flexibility in
optimizing (3) and allows for heterogeneous clients. More specifically, in
Sec. 4.3, we will show two approaches to design $\bm{W}^{(t)}$ at each
iteration $t$ such that $\bm{W}^{(t)}\to\bm{I}$. The intuition is as follows.
In the early stages, client $k$ may not have enough observations to obtain a
high-quality surrogate model and therefore needs to borrow information from
other clients. In the late stages, as client $k$ accumulates sufficient data
and can construct a high-quality local surrogate model, it will focus more on
its local design problem to find client-specific optimal design points.

Finally, the consensus constraint provides a naturally distributed approach to
solve (3) sequentially. To see this, recall that each client can run $T$
experiments across $T$ iterations and the initial dataset is
$\mathcal{D}_{k}^{(0)}$. Then a natural algorithm would take the following
form: At iteration $t\in\\{0,\ldots,T-1\\}$, given
$\bm{x}_{\mathcal{C}}^{(t)}$, all clients calculate
$\bm{x}_{k}^{(t)\text{new}}=[(\bm{W}^{(t)}\otimes\bm{I}_{D})\bm{x}_{\mathcal{C}}^{(t)}]_{k}$.
Now, each client $k$ conducts an experiment at $\bm{x}_{k}^{(t)\text{new}}$,
and augments its dataset with the new observation. Each client then fits a
surrogate (e.g., $\mathcal{GP}$) using the current dataset
$\mathcal{D}_{k}^{(t+1)}$ and accordingly constructs the utility function
$U(\hat{y}_{k}(\bm{x}_{k}))$. Then the clients maximize their local utility
$\max_{\bm{x}_{k}}\mathbb{E}_{\hat{y}_{k}|\mathcal{D}^{(t+1)}_{k}}\left[U(\hat{y}_{k}(\bm{x}_{k}))\right]$
to obtain a new candidate solution $\bm{x}^{(t+1)}_{k}$. Finally, all clients
send their solutions to each other or a central orchestrator, and the process
repeats.

Algorithm 1 summarizes our Collaborative Bayesian Optimization via Consensus
(CBOC) framework. Fig. 2 presents a flowchart illustrating the collaborative
BO. Notably, our collaborative framework enjoys some nice theoretical
properties (See Sec. 4.5).

Data: $T$ iterations, $K$ number of clients, initial data
$\\{\mathcal{D}^{(0)}_{k}\\}_{k=1}^{K}$, initial consensus matrix
$\bm{W}^{(0)}$, initial designs to share $\\{\bm{x}_{k}^{(0)}\\}_{k=1}^{K}$

for  _$t$ $=0,1,\cdots T-1$_ do

for  _$k=1,\cdots,K$_ do

Consensus: client $k$ computes the consensus point
$\bm{x}_{k}^{(t)\text{new}}=[(\bm{W}^{(t)}\otimes\bm{I}_{D})\bm{x}_{\mathcal{C}}^{(t)}]_{k}$
;

Experiment: client $k$ conducts experiment using $\bm{x}_{k}^{(t)\text{new}}$
and observes $y_{k}(\bm{x}_{k}^{(t)\text{new}})$;

Data Augmentation: client $k$ augments dataset by
$\mathcal{D}_{k}^{(t+1)}=\mathcal{D}_{k}^{(t)}\cup\\{\bm{x}_{k}^{(t)\text{new}},y_{k}(\bm{x}_{k}^{(t)\text{new}})\\}$
;

Surrogate Modeling: client $k$ updates its surrogate model ;

Optimization: client $k$ finds their local utility maximizer
$\bm{x}^{(t+1)}_{k}=\operatorname*{argmax}_{\bm{x}_{k}}\mathbb{E}_{\hat{y}_{k}|\mathcal{D}^{(t+1)}_{k}}\left[U(\hat{y}_{k}(\bm{x}_{k}))\right]$
and shares this candidate sample ;

end for

end for

Algorithm 1 CBOC: Collaborative BO via Consensus

Figure 2: Flowchart of CBOC.

It is worthwhile noting that the collaborative process can be done in a
centralized or decentralized manner where clients share $\bm{x}^{(t)}_{k}$
either with a central orchestrator or with each other directly.

###  4.3 Designing the Consensus Matrix

From Algorithm 1, it can be seen that one of the key components in (3) is the
consensus matrix $\bm{W}$. In essence, as we have discussed earlier, the
consensus matrix $\bm{W}$ controls how much one client will affect the design
choices for other clients. Therefore, one needs to carefully design the
consensus matrix.

In this section, we propose two approaches for designing the consensus matrix.
The first approach assumes one does not have any prior information on
different clients and uniformly adjusts all entries in $\bm{W}^{(t)}$ at each
iteration $t$. The second approach carefully modifies the weights for each
client based on the “leader”, where the “leader” is defined as the client that
has observed the best improvement (e.g., the most significant improvement in
the utility). Below, we will detail both approaches.

##### Uniform Transitional Matrix

We propose to assign equal weights to all clients at iteration $0$ and then
gradually decay off-diagonal elements

$\displaystyle\bm{W}^{(0)}=\begin{bmatrix}\frac{1}{K}&\frac{1}{K}&\ldots&\frac{1}{K}\\\
\vdots&\vdots&\vdots&\vdots\\\
\frac{1}{K}&\frac{1}{K}&\ldots&\frac{1}{K}\end{bmatrix}.$

For simplicity, we assume all clients are connected. In practice, if clients
$i$ and $j$ are not connected, we can set $w_{ij}=w_{ji}=0$ and reweigh the
other components of the $\bm{W}$ matrix to ensure it is doubly stochastic.

Figure 3: Illustration of the collaborative process. Colors represent client-
specific black-box functions.

As the optimization progresses, we gradually increase diagonal elements
$w_{kk}$ to 1 and decrease the off-diagonal elements $w_{kj}$ for $k\neq j$ to
0. For example, if we adjust weights linearly, then

$\displaystyle\bm{W}^{(t+1)}=\bm{W}^{(t)}+\begin{bmatrix}\frac{K-1}{TK}&-\frac{1}{TK}&\ldots&-\frac{1}{TK}\\\
\vdots&\vdots&\vdots&\vdots\\\
-\frac{1}{TK}&-\frac{1}{TK}&\ldots&\frac{K-1}{TK}\end{bmatrix}.$ (5)

Ultimately, $\bm{W}^{(t)}$ will converge to the identity matrix $\bm{I}$. The
intuition behind this design is illustrated in Fig. 3. In this figure, each
colored line represents a black-box design function for a client. In the early
stages, client $k$ does not have enough observations to obtain a high-quality
utility and therefore needs to borrow more information from other clients. As
client $k$ has more data in the late stages, it will focus more on its own
design problem to find client-specific optimal designs.

Note that $\bm{W}^{(t)}\rightarrow\bm{I}$, is mandatory if clients have some
heterogeneity. To see this, assume all clients have recovered their optimal
design $\bm{x}_{k}^{*}$, then if $\bm{W}\neq\bm{I}$, the consensus step will
always move the experiment away from $\bm{x}_{k}^{*}$.

##### Leader-driven Matrix

Another approach is to adjust the weights of $\bm{W}$ dynamically based on the
“leader”, i.e., the client that observed the best improvement. As a result,
all other clients will follow the leader and explore the region that
potentially contains the best solutions. Initially, we create two matrices
$\bm{W}^{(0)}_{1}$ and $\bm{W}^{(0)}_{2}$, and assign equal weights to all
entries of both matrices. Matrix $\bm{W}_{1}$ will be used as a baseline, and
is updated using (5). Matrix $\bm{W}_{2}$ will be modified based on
$\bm{W}_{1}$ and by incorporating the leader information. The purpose of using
$\bm{W}_{1}$ is to ensure that $\bm{W}_{2}$ still converges to $\bm{I}$ so
that every single client will eventually focus on their own design objective.

At iteration $t$, each client $k$ shares a pair
$(\bm{x}^{(t)}_{k},S^{(t)}_{k})$, where $S_{k}$ is a reward that quantifies
the gained benefit at client $k$. As a simple example, we could set
$S^{(t)}_{k}=\max\mathbb{E}_{\hat{y}_{k}|\mathcal{D}^{(t)}_{k}}\left[U(\hat{y}_{k}(\bm{x}_{k}))\right]$.
We then sort $S^{(t)}_{k}$ and find
$k^{*(t)}=\operatorname*{argmax}_{k}\\{S^{(t)}_{1},S^{(t)}_{2},\ldots,S^{(t)}_{K}\\}$,
and treat client $k^{*(t)}$ as a leader at the current iteration. In essence,
client $k^{*(t)}$ can hint to others that the neighborhood of
$\bm{x}_{k^{*(t)}}^{(t)}$ is an area of potential improvement. To do so,
during the consensus step, for the non-lead clients $k\neq k^{*(t)}$, we
decrease their weights by $\frac{1}{TK}$ and increase their off-diagonal
weights with the leader $k^{*(t)}$ by $\frac{K-1}{TK}$. To this end, we adjust
blocks of $\bm{W}^{(t)}_{2}$ as follows:

$\displaystyle[\bm{W}^{(t)}_{2}]_{\neq
k^{*(t)},k^{*(t)}}=[\bm{W}^{(t)}_{1}]_{\neq
k^{*(t)},k^{*(t)}}+\frac{K-1}{TK},$
$\displaystyle[\bm{W}^{(t)}_{2}]_{k^{*(t)},\neq
k^{*(t)}}=[\bm{W}^{(t)}_{1}]_{k^{*(t)},\neq k^{*(t)}}+\frac{K-1}{TK},$
$\displaystyle[\bm{W}^{(t)}_{2}]_{\neq k^{*(t)},\neq
k^{*(t)}}=[\bm{W}^{(t)}_{1}]_{\neq k^{*(t)},\neq k^{*(t)}}-\frac{1}{TK},$
$\displaystyle[\bm{W}^{(t)}_{2}]_{k^{*(t)},k^{*(t)}}=[\bm{W}^{(t)}_{1}]_{k^{*(t)},k^{*(t)}}-\frac{(K-1)^{2}}{TK},$

where $[\bm{W}_{2}^{(t)}]_{\neq i,j}$ represents all elements that are in the
$j^{th}$ column but not in the $i^{th}$ row of $\bm{W}_{2}^{(t)}$.
Interestingly, the weight $[\bm{W}^{(t)}_{2}]_{k^{*(t)},k^{*(t)}}$ shrinks to
maintain double stochasticity. Such an assignment allows the leader $k^{*(t)}$
to explore new regions.

Here we also suggest two heuristics: (i) To avoid the situation where the same
client is selected in succession (this will cause the same client to keep
exploring rather than exploiting), we propose to select the second largest
index from $\\{S^{(t+1)}_{1},S^{(t+1)}_{2},\ldots,S^{(t+1)}_{K}\\}$ if
$k^{*(t)}=k^{*(t+1)}$. (ii) In the extreme case where
$[\bm{W}^{(t)}_{2}]_{k^{*(t)},k^{*(t)}}$ is negative, set
$[\bm{W}^{(t)}_{2}]_{k^{*(t)},k^{*(t)}}$ to zero and reweigh other components
accordingly.

To see an example, consider a scenario where $D=1$, $K=3$ and $T=10$. We
initialize both $\bm{W}^{(0)}_{1}$ and $\bm{W}^{(0)}_{2}$ as

$\displaystyle\bm{W}^{(0)}_{1}=\bm{W}^{(0)}_{2}=\begin{bmatrix}\frac{1}{3}&\frac{1}{3}&\frac{1}{3}\\\
\frac{1}{3}&\frac{1}{3}&\frac{1}{3}\\\
\frac{1}{3}&\frac{1}{3}&\frac{1}{3}\end{bmatrix}.$

Suppose each client received
$(\bm{x}^{(0)}_{1},1),(\bm{x}^{(0)}_{2},5),(\bm{x}^{(0)}_{3},4)$. Then we
adjust $\bm{W}^{(0)}_{2}$ as

$\displaystyle\bm{W}^{(0)}_{2}=\begin{bmatrix}\frac{1}{3}-\frac{1}{30}&\frac{1}{3}+\frac{2}{30}&\frac{1}{3}-\frac{1}{30}\\\
\frac{1}{3}+\frac{2}{30}&\frac{1}{3}-\frac{4}{30}&\frac{1}{3}+\frac{2}{30}\\\
\frac{1}{3}-\frac{1}{30}&\frac{1}{3}+\frac{2}{30}&\frac{1}{3}-\frac{1}{30}\end{bmatrix}.$

In this example, client 2 observed that sampling at $\bm{x}^{(0)}_{2}$ yields
the highest benefits, and therefore we put more weight on client 2. As a
result, clients 1 and 3 will explore towards $\bm{x}^{(0)}_{2}$. On the other
hand, we decrease $[\bm{W}_{2}^{(0)}]_{2,2}$ to make sure that client 2 does
not over-explore the region that contains $\bm{x}^{(0)}_{2}$.

###  4.4 Contextualization under a specific surrogate and utility

Now, given the generic framework presented in 4.2 and 4.3, we will
contextualize Algorithm 1 and provide a concrete example. The iteration index
$t$ is dropped for simplicity unless necessary.

####  4.4.1 Gaussian Process Surrogate

At each iteration, we place a $\mathcal{GP}$ prior on the surrogate
$\hat{f}_{k}$. We model the additive noise term $\epsilon_{k}$ as independent
and identically distributed ($i.i.d.$) noise that follows a normal
distribution with zero mean and $v_{k}^{2}$ variance. Now, given a new input
location $\bm{x}^{\text{test}}$, the posterior predictive distribution of
$\hat{f}_{k}(\bm{x}^{\text{test}})\sim\mathbb{P}_{\hat{f}_{k}(\bm{x}^{\text{test}})|\mathcal{D}^{(t)}_{k}}\coloneqq\mathcal{N}(\mu_{k}(\bm{x}^{\text{test}};\mathcal{D}_{k}),\sigma^{2}_{k}(\bm{x}^{\text{test}};\mathcal{D}_{k}))$,
where

$\begin{split}&\mu_{k}(\bm{x}^{\text{test}};\mathcal{D}_{k})=\bm{K}(\bm{x}^{\text{test}},\bm{X}_{k})\left(\bm{K}(\bm{X}_{k},\bm{X}_{k})+v_{k}^{2}\bm{I}\right)^{-1}\bm{y}_{k},\\\
&\sigma^{2}_{k}(\bm{x}^{\text{test}};\mathcal{D}_{k})=\bm{K}(\bm{x}^{\text{test}},\bm{x}^{\text{test}})-\bm{K}(\bm{x}^{\text{test}},\bm{X}_{k})\left(\bm{K}(\bm{X}_{k},\bm{X}_{k})+v_{k}^{2}\bm{I}\right)^{-1}\bm{K}(\bm{X}_{k},\bm{x}^{\text{test}}),\end{split}$
(6)

and $\bm{K}(\cdot,\cdot):\mathbb{R}^{D}\times\mathbb{R}^{D}\to\mathbb{R}$ is a
covariance matrix whose entries are determined by some kernel function
$\mathcal{K}(\cdot,\cdot)$. Similarly,
$\mathbb{P}_{\hat{y}_{k}|\mathcal{D}^{(t)}_{k}}$ is derived by simply adding
$v_{k}^{2}$ to $\sigma^{2}_{k}(\bm{x}^{\text{test}})$.

####  4.4.2 Utility Function

Given the $\mathcal{GP}$ surrogate, one can build a utility function that
measures the benefits of conducting an experiment using a set of new design
points. One common example is the expected improvement (EI) utility expressed
as (Jones et al.,, 1998)

$\displaystyle\mathbb{E}_{\hat{f}_{k}|\mathcal{D}^{(t)}_{k}}\left[U(\hat{f}_{k}(\bm{x}))\right]$
$\displaystyle=\mathbb{E}_{\hat{f}_{k}|\mathcal{D}^{(t)}_{k}}\left[(\hat{f}_{k}(\bm{x})-y_{k}^{*(t)})^{+}\right]=\texttt{EI}^{(t)}_{k}(\bm{x})$
$\displaystyle=\sigma^{(t)}_{k}(\bm{x};\mathcal{D}^{(t)}_{k})\phi(z^{(t)}_{k}(\bm{x}))+(\mu^{(t)}_{k}(\bm{x};\mathcal{D}^{(t)}_{k})-y^{*(t)}_{k})\Phi(z^{(t)}_{k}(\bm{x})),$

where $a^{+}=\max(a,0)$, $y_{k}^{*(t)}=\max\bm{y}^{(t)}_{k}$ is the current
best response, $\phi(\cdot)$ (or $\Phi(\cdot)$) is a probability density
function (PDF) (or cumulative distribution function (CDF)) of a standard
normal random variable, and
$z^{(t)}_{k}(\bm{x})=\frac{\mu^{(t)}_{k}(\bm{x};\mathcal{D}^{(t)}_{k})-y^{*(t)}_{k}}{\sigma^{(t)}_{k}(\bm{x};\mathcal{D}^{(t)}_{k})}$.
Here note that the expectation is taken with respect to $\hat{f}_{k}$ rather
than $\hat{y}_{k}$. Another example is the knowledge gradient (KG) (Wu and
Frazier,, 2016) defined as
$\mathbb{E}_{\hat{y}_{k}|\mathcal{D}^{(t)}_{k}}\left[U(\hat{y}_{k}(\bm{x}))\right]=\texttt{KG}^{(t)}_{k}(\bm{x})=\mathbb{E}_{\hat{y}_{k}|\mathcal{D}^{(t)}_{k}}\left[\mu^{*{(t)}}_{k}(\bm{x};\mathcal{D}^{(t)}_{k}\cup\\{\bm{x},\hat{y}_{k}(\bm{x})\\})\right]-\max_{\bm{x}^{\prime}}(\mu^{(t)}_{k}(\bm{x}^{\prime};\mathcal{D}^{(t)}_{k}))$,
where
$\mu^{*(t)}_{k}(\bm{x};\mathcal{D}^{(t)}_{k}\cup\\{\bm{x},\hat{y}_{k}(\bm{x})\\})$
is the maximum of the updated posterior mean of the $\mathcal{GP}$ surrogate
by taking one more sample at location $(\bm{x},\hat{y}_{k}(\bm{x}))$. KG can
be interpreted as finding the new sampling location $\bm{x}$ that potentially
increases the maximum updated posterior mean. Hereon, in the remainder of this
paper, we drop $\mathcal{D}_{k}$ in $\mu_{k},\sigma_{k}$ for the sake of
compactness.

###  4.5 Theoretical Analysis

Despite its immense success, BO theory is still in its infancy due to many
fundamental challenges. First and foremost, the black-box nature of $f_{k}$
renders theory hard to derive due to the lack of known structure. Second, we
still have a limited understanding of the properties, such as Lipschitz
continuity, concavity, or smoothness, of many commonly employed utility
functions. For example, even the EI utility in general is not Lipschitz
continuous or concave. Third, despite recent advances in laying the
theoretical foundations for understanding the generalization error bounds of
$\mathcal{GP}$s (Lederer et al.,, 2019; Wang et al.,, 2020), bridging the gap
between these bounds and errors incurred in the utility function remains an
open problem and a rather challenging one.

To circumvent these open problems while providing a theoretical proof of
concept, we derive some theoretical insights that are confined to the EI
utility in conjunction with a smooth $\mathcal{GP}$ kernel and a homogeneity
assumption. While we believe that our results extend beyond these settings, we
leave this analysis as an enticing challenge for future research.

We focus on regret defined as
$r^{(t)}_{k}=f_{k}(\bm{x}_{k}^{*})-f_{k}(\bm{x}^{(t)\text{new}}_{k})$ for
client $k$ at iteration $t$. Intuitively, regret measures the gap between the
design function evaluated at the optimal solution $\bm{x}_{k}^{*}$ and the one
evaluated at the consensus solution $\bm{x}^{(t)\text{new}}_{k}$, at iteration
$t$. By definition, $r^{(t)}_{k}=0$ if CBOC recovers the global optimal
solution. Our theoretical guarantee studies cumulative regret, more
specifically, $R_{k,T}=\sum_{t=1}^{T}r^{(t)}_{k}$.

By definition,
$\texttt{EI}_{k}^{(t)}(\bm{x})=\mathbb{E}_{\hat{f}_{k}|\mathcal{D}^{(t)}_{k}}\left[(\hat{f}_{k}(\bm{x})-y_{k}^{*(t)})^{+}\right]\geq
0,\forall\bm{x}\in\mathbb{R}^{D}$. Therefore, we define a small positive
constant $\kappa$ such that when $\texttt{EI}_{k}^{(t)}(\bm{x})<\kappa$,
client $k$ stops its algorithm at iteration $t$. This stopping criterion is
only used for theoretical development. In practice, we will run our algorithm
for $T$ iterations or until all budgets are exhausted.

Below, we present our main Theorem and the sketch of our proof. Detailed
information and all supporting Lemmas and required assumptions can be found in
Appendices 1-2.

######  Theorem 1.

(Homogeneous Clients) Suppose $f_{1}=f_{2}=\cdots=f_{K}$, and suppose a
squared exponential kernel function
$\mathcal{K}_{k}(\bm{x},\bm{x}^{\prime})=u_{k}^{2}\exp\left(\frac{\left\lVert\bm{x}-\bm{x}^{\prime}\right\rVert^{2}}{2\ell_{k}^{2}}\right)$
is used for the $\mathcal{GP}$ surrogate, where $u_{k}$ is the variance scale
parameter and $\ell_{k}$ is the length parameter, and each client uses the EI
utility. Without loss of generality, assume
$\mathcal{K}_{k}(\bm{x},\bm{x}^{\prime})\leq 1$ and the initial sample size
for each client is $2$. Under some assumptions (Appendix 1), given any doubly
stochastic $\bm{W}^{(t)}$ with non-negative elements, for
$\epsilon>0,\delta_{1}\in(0,T)$, with probability at least
$(1-\frac{\delta_{1}}{T})^{T}$, the cumulative regret after $T>1$ iterations
is

$\displaystyle R_{k,T}$
$\displaystyle=\sum_{t=1}^{T}r_{k}^{(t)}\leq\sqrt{\frac{6T\left[(\log
T)^{3}+1+C\right](\log T)^{D+1}}{\log(1+v_{k}^{-2})}}$
$\displaystyle\qquad+\sqrt{\frac{2T(\log
T)^{D+4}}{\log(1+v_{k}^{-2})}}+\sum_{t=1}^{T}\mathcal{O}\left(\frac{1}{(\log(2+t))^{0.5+\epsilon}}\right)$
$\displaystyle\sim\mathcal{O}(\sqrt{T\times(\log T)^{D+4}}),$

where $C=\log[\frac{1}{2\pi\kappa^{2}}]$.

Theorem 1 shows that the cumulative regret of Algorithm 1 has a sublinear
growth rate in terms of the number of iterations. This implies that as the
algorithm proceeds with more iterations, the consensus solution for client $k$
will be close to the optimal solution $\bm{x}_{1}^{*}=\cdots=\bm{x}_{k}^{*}$.

Theorem 1 relies on the assumption of homogeneity. When heterogeneity is
present, where at least one pair of design functions has some differences,
providing convergence results is challenging. To our knowledge, the theory
here is vacant. Even if one were to provide a regularity assumption on the
heterogeneity, say $f_{k}$’s are uniformly bounded and have optimal designs
restricted within a small ball, it remains an open problem to understand how
this regularity will propagate to the utility where we are able to understand
the structural similarities and differences across
$\\{\texttt{EI}_{k}\\}_{k=1}^{K}$. Indeed, this stems from the aforementioned
second and third challenges in traditional BO. That said, Theorem 1 may serve
as a proof of concept, showcasing that our approach will drive clients to a
region that contains optimal designs. Afterward, each client will focus on
their own objective to obtain a client-specific optimal design (recall Fig.
3). As shown in Sec. 4.3 and Sections 5-6, our empirical results highlight the
advantageous properties of our approach under heterogeneous settings.

##### Proof Sketch

We hereby provide a proof sketch for Theorem 1. We first decompose the regret
$r_{k}^{(t)}$ as

$\displaystyle
r_{k}^{(t)}=\bigg{\\{}\underbrace{f_{k}(\bm{x}_{k}^{*})-y^{*(t)}_{k}}_{\text{A}}+\underbrace{y^{*(t)}_{k}-f_{k}(\bm{x}^{(t)\text{new}}_{k})}_{\text{B}}\bigg{\\}}.$

We let $\text{A}=f_{k}(\bm{x}_{k}^{*})-y^{*(t)}_{k}$ and
$\text{B}=y^{*(t)}_{k}-f_{k}(\bm{x}^{(t)\text{new}}_{k})$. Next, we aim to
bound terms A and B.

  * •

Term A defines the difference between the optimal value of the design function
and the current best-observed output value. Here, we first use a concentration
inequality from Srinivas et al., (2009) to bound the difference between the
mean function $\mu_{k}^{(t)}(\cdot)$ of the $\mathcal{GP}$ surrogate and the
truth $f_{k}(\cdot)$. We then show that, with a high probability, term A is
upper bounded by the EI value evaluated at the optimal design $\bm{x}_{k}^{*}$
plus a scaled predictive standard deviation term. Mathematically, with
probability $1-\delta$, $\delta\in(0,1)$,

$f_{k}(\bm{x}_{k}^{*})-y^{*(t)}_{k}\leq\texttt{EI}_{k}^{(t)}(\bm{x}_{k}^{*})+\sqrt{\beta_{k}^{(t)}}\sigma_{k}^{(t)}(\bm{x}_{k}^{*}),$

where $\\{\beta_{k}^{(t)}\\}_{t}$ is a non-decreasing sequence such that
$\beta_{k}^{(t)}\sim\mathcal{O}((\log\frac{t}{\delta})^{3})$.

  * •

Term B defines the difference between the current best-observed output value
and the underlying true function $f_{k}(\cdot)$ evaluated at the consensus
solution. To proceed, we further expand B as

B
$\displaystyle=y^{*(t)}_{k}-f_{k}(\bm{x}^{(t)\text{new}}_{k})=y^{*(t)}_{k}-\mu_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})+\mu_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})-f_{k}(\bm{x}^{(t)\text{new}}_{k}).$

We then bound B using a similar strategy in bounding term A and obtain

$\displaystyle\text{B}\leq\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\left(\tau(-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))-\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))+\sqrt{\beta_{k}^{(t)}}\right),$

where $\tau(z_{k}^{(t)}(\bm{x}))\coloneqq
z_{k}^{(t)}(\bm{x})\Phi(z_{k}^{(t)}(\bm{x}))+\phi(z_{k}^{(t)}(\bm{x}))$.

Our next goal is to show that the summation of $r_{k}^{(t)}$ over $T$
iterations (i.e., the cumulative regret) is bounded. Here, note that there are
two key components that appear in A and B: $\sigma^{(t)}_{k}(\cdot)$ and
$\tau(\cdot)$. First, we analyze the behavior of the cumulative predictive
variance $\sum_{t=1}^{T}\sigma^{2(t)}_{k}(\cdot)$. We show that, with a
squared exponential kernel,
$\sum_{t=1}^{T}\sigma^{2(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\leq\frac{2}{\log(1+v_{k}^{-2})}\mathcal{O}((\log
T)^{D+1})$. Second, we show that
$\tau(-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))\leq 1+\sqrt{C}$, where
$C=\log[\frac{1}{2\pi\kappa^{2}}]$. After some algebraic manipulation, we
obtain the upper bound stated in Theorem 1.

∎

##  5 Simulation Studies

In this section, we validate CBOC on several simulation datasets. We consider
a range of simulation functions from the Virtual Library of Simulation
Experiments (Surjanovic and Bingham,, 2013) as the underlying black-box design
functions for the clients. Our goal is to showcase the benefit of
collaboration in finding client-specific optimal designs.

We set the number of iterations for each testing function to $T=20D$. At
iteration $0$, $\bm{W}^{(0)}$ is initialized as a uniform matrix where each
entry has a weight equal to $\frac{1}{K}$. We use the uniform transitional
approach (CBOC-U) or the leader-driven approach (CBOC-L) to adjust
$\bm{W}^{(t)}$ at every iteration $t$. The initial dataset for each client
$\\{\mathcal{D}^{(0)}_{k}\\}_{k=1}^{K}$ contains $5D$ randomly chosen design
points. The performance of each client is evaluated using a Gap metric $G_{k}$
(Jiang et al.,, 2020) defined as

$\displaystyle
G_{k}=\frac{|y^{*(0)}_{k}-y^{*(T)}_{k}|}{|y^{*(0)}_{k}-y^{*}_{k}|},$

where $y^{*(0)}_{k}$ (or $y^{*(T)}_{k}$) is the best observed response at
iteration $0$ (or $T$), and $y^{*}_{k}=f_{k}(\bm{x}^{*}_{k})$ is the true
optimal response. A larger $G_{k}$ implies better performance. If the optimal
solution is recovered (i.e., $y^{*(T)}_{k}=y^{*}_{k}$), then $G_{k}=1$. We
compare CBOC using the EI utility with the following benchmark models: (1)
Individual: each client $k$ conducts BO without collaboration; (2) FedBO: the
state-of-the-art federated BO algorithm that builds upon Thompson sampling
(Dai et al.,, 2021). For both benchmarks, we collect $5D$ initial samples and
conduct experiments for $20D$ iterations. Our code is available at this GitHub
link.

###  5.1 Testing Function 1: Levy-$D$

Levy-$D$ is an $D$-dimensional function in the form of

$\displaystyle
f(\bm{x})=\sin^{2}(\pi\omega_{1})+\sum_{d=1}^{D-1}(\omega_{d}-1)^{2}[1+10\sin^{2}(\pi\omega_{d}+1)]+(\omega_{D}-1)^{2}[1+\sin^{2}(2\pi\omega_{D})],$

where $\omega_{d}=1+\frac{\tensor[^{d}]{x}{}-1}{4},\forall d\in[D]$ and
$\bm{x}=(\tensor[^{1}]{x}{},\ldots,\tensor[^{d}]{x}{}\ldots,\tensor[^{D}]{x}{})^{\intercal}\in[-10,10]^{D}$.
We first consider the homogeneous scenario where
$f_{1}(\bm{x})=\cdots=f_{K}(\bm{x})=f(\bm{x})$. We focus on maximizing the
design function $-f_{k}(\bm{x})$ (i.e., minimizing $f_{k}(\bm{x})$), $\forall
k\in[K]$. In Table 1, we report the average Gap across $K=5$ clients over 30
independent runs, defined as follows:

$\displaystyle\bar{G}=\frac{1}{30}\sum_{i=1}^{30}\frac{1}{K}\sum_{k=1}^{K}G^{(i)}_{k},$

where $G^{(i)}_{k}$ is the Gap metric for client $k$ at the $i$-th run.

Functions | CBOC-L | Individual | FedBO
---|---|---|---
Levy-$2$ | $0.993(\pm 0.002)$ | $0.931(\pm 0.004)$ | $0.990(\pm 0.002)$
Levy-$4$ | $0.987(\pm 0.003)$ | $0.926(\pm 0.007)$ | $0.951(\pm 0.006)$
Levy-$8$ | $0.969(\pm 0.010)$ | $0.925(\pm 0.013)$ | $0.938(\pm 0.009)$
Table 1: The average Gap across $K=5$ clients over 30 independent runs under a
homogeneous setting. We report standard deviations over 30 runs in brackets.

Second, we consider the heterogeneous scenario where each client has a
different underlying truth. To do so, for each client $k$, we transform the
Levy function $f(\bm{x})$ to

$\displaystyle f_{k}(\bm{x})=a_{1}f(\bm{x}+\text{vec}(a_{3}))+a_{2}\,,$

where $a_{1}\sim\texttt{Uniform}(0.5,1)$, $a_{2},a_{3}\sim\mathcal{N}(0,1)$,
and $\text{vec}(a_{3})$ is a $D$-dimensional vector whose elements are all
equal to $a_{3}$. This transformation will shift and re-scale the original
function, creating heterogeneous functional forms. The homogeneous scenario
can be viewed as a special case when $a_{1}=1,a_{2}=a_{3}=0$. We set $K=10$.
Other settings remain unchanged. Results are shown in Table 2.

Functions | CBOC-L | Individual | FedBO
---|---|---|---
Levy-$2$ | $0.990(\pm 0.001)$ | $0.942(\pm 0.005)$ | $0.958(\pm 0.003)$
Levy-$4$ | $0.984(\pm 0.002)$ | $0.933(\pm 0.012)$ | $0.940(\pm 0.009)$
Levy-$8$ | $0.949(\pm 0.008)$ | $0.917(\pm 0.008)$ | 0.$903(\pm 0.011)$
Table 2: The average Gap across $K=10$ clients over 30 independent runs under
a heterogeneous setting.

From Tables 1-2, we can derive two key insights. First, collaborative methods
yield superior performance than non-collaborative competitors through higher
average Gap metrics. This evidences the importance of collaboration in
improving the optimal design process. Second, CBOC outperforms all benchmarks.
This credits to CBOC’s ability to address heterogeneity through a flexible
consensus framework that allows clients to collaboratively explore and exploit
the design space and eventually obtain client-specific solutions.

###  5.2 An illustrative example

To visualize the performance of our method, we provide an illustration using
the Levy-$2$ function, $K=2$, and the same heterogeneity structure as in Sec.
5.1. Specifically, the two design functions are set to
$f_{1}(\bm{x})=f(\bm{x}+\text{vec}(1))+1,f_{2}(\bm{x})=2f(\bm{x}+\text{vec}(2))+2$.
Here, each client starts with five two-dimensional designs, and then both
CBOC-L and Individual are run for $T=40$. Fig. 4 demonstrates the landscape of
the original Levy-$2$ function and shows contour plots of $f_{1},f_{2}$. The
selected design points for each client are marked as red for CBOC-L and green
for Individual, and the iteration indices are labeled next to those points. We
did not label all points for better visualization.

From Fig. 4, we can see that through collaboration CBOC-L samples more
frequently near optimal regions, as evident by the larger number of red points
close to the optimal. This allows CBOC-L to reach the region that contains the
optimal design (white-most region) for both clients in 25 iterations. On the
other hand, Individual required around 38 iterations for client 1, while for
client 2, even at $T=40$, the optimal region was still not explored yet. This
again highlights the benefit of collaboration. More importantly, Fig. 4
highlights the ability of our method to operate under heterogeneity as both
clients were able to reach their distinct optimal design neighborhood, and
they do so much faster than operating in isolation.

Figure 4: Contour plot of $f_{1},f_{2}$ and selected design points. The x-axis
shows the first design input (i.e., $\tensor[^{1}]{x}{}$), and the y-axis
shows the second design input (i.e., $\tensor[^{2}]{x}{}$).

###  5.3 Testing Function 2: Shekel-$10$

Shekel-$10$ is a four-dimensional ($D=4$) with $10$ local minima. It has the
functional form:

$\displaystyle
f(\bm{x})=-\sum_{i=1}^{10}\left(\sum_{d=1}^{4}(\tensor[^{d}]{x}{}-F_{di})^{2}+\xi_{i}\right)^{-1}.$

We defer the specification of $F_{di}$ and $\xi_{i}$ to Appendix 3.

Similar to the heterogeneous scenario in Sec. 5.1, for each client $k$, we
transform the Shekel-$10$ function $f(\bm{x})$ to
$f_{k}(\bm{x})=a_{1}f(\bm{x}+\text{vec}(a_{3}))+a_{2}$, where
$a_{1}\sim\texttt{Uniform}(0.5,1)$ and
$a_{2}\sim\mathcal{N}(0,2),a_{3}\sim\mathcal{N}(0,1)$. We test all benchmarks
using $K=5,10,15$ and $20$ clients.

Functions | CBOC-L | CBOC-U | Individual | FedBO
---|---|---|---|---
$K=5$ | $0.475(\pm 0.053)$ | $0.462(\pm 0.052)$ | $0.350(\pm 0.055)$ | $0.370(\pm 0.047)$
$K=10$ | $0.516(\pm 0.049)$ | $0.501(\pm 0.029)$ | $0.364(\pm 0.035)$ | $0.422(\pm 0.040)$
$K=15$ | $0.577(\pm 0.033)$ | $0.553(\pm 0.036)$ | $0.356(\pm 0.022)$ | $0.496(\pm 0.053)$
$K=20$ | $0.592(\pm 0.036)$ | $0.572(\pm 0.015)$ | $0.335(\pm 0.028)$ | $0.535(\pm 0.051)$
Table 3: The average Gap across $K$ clients over 30 independent runs under a
heterogeneous setting.

Table 3 shows that CBOC outperforms both the non-collaborative method and
state-of-the-art FedBO. Interestingly, we observe that the performance of CBOC
improves as more clients participate in the collaboration process. For
example, the average Gap for CBOC when $K=5$ is 0.475. As we increase $K$ to
$20$, the average Gap becomes 0.592. Fig. 5 plots the evolution of the average
Gap with respect to iterations. This result further demonstrates the benefits
of collaboration.

Figure 5: The evolution of the average Gap with respect to iterations and
collaborators.

###  5.4 Other Testing Functions

Finally, we test our methods on three other functions: Branin, Ackley-$D$, and
Harmann-$6$. We use $K=10$. Experimental details and function specifications
are deferred to the appendix.

Functions | CBOC-L | CBOC-U | Individual | FedBO
---|---|---|---|---
Branin | $0.992(\pm 0.000)$ | $0.990(\pm 0.001)$ | $0.975(\pm 0.001)$ | $0.986(\pm 0.001)$
Ackley-$5$ | $\bm{0.656(\pm 0.057)}$ | $0.641(\pm 0.053)$ | $0.501(\pm 0.039)$ | $0.632(\pm 0.041)$
Hartmann-$6$ | $0.968(\pm 0.005)$ | $0.959(\pm 0.002)$ | $0.941(\pm 0.001)$ | $0.955(\pm 0.003)$
Table 4: The average Gap across $K=10$ clients over 30 independent runs under
a heterogeneous setting.

Similar to our previous analysis, Table 4 shows the superiority of our
approach.

##  6 Case Study

A case study on consensus BO-driven closed-loop biosensor design optimization
was performed (see Fig. 6). The closed-loop workflow was based on simulation
(FEA) guided by a decentralized computing process for collaboration (i.e.,
consensus BO). Given the demand for device-based biosensors for industrial
process analytical technology (e.g., for bioprocess monitoring and control)
and health monitoring (e.g., via wearable sensors), the case study was focused
on closed-loop optimization of a device-based biosensor. Specifically, our
case study focused on optimizing a device-based biosensor with a milli-scale
transducer used in a continuous flow format, which is consistent with several
types of device-based biosensors, including electrochemical, mechanical, and
electromechanical biosensors.

Figure 6:  Case study schematic.

Biosensor performance, specifically, transient binding of the target analyte,
was calculated using commercially available FEA software (COMSOL Multiphysics,
COMSOL). The transient fractional surface coverage of the bound target analyte
on the biosensor surface was calculated by numerical solution of a coupled
convection-diffusion-reaction model in 2D using a time-dependent study. The
convection-diffusion-reaction model was constructed by coupling laminar flow
and transport of diluted species physics with a surface reaction. The
computational domain consisted of a 6 $\times$ 1 mm2 (length × width)
rectangular fluidic channel that encompassed the domain $x\in[-3,3]$ mm and
$y\in[0,1]$ mm. The biosensor surface (i.e., surface 1) on which the binding
reaction between immobilized biorecognition elements and target analyte
occurred encompassed the domain $x\in[0,1]$ mm at $y=0$. The computational
domain contained five additional surfaces on which boundary and initial
conditions were applied. The coupled partial differential equations associated
with the convection-diffusion-reaction model were solved subject to the
following boundary and initial conditions:

For laminar flow physics:

  * •

Surface 2 (inlet): normal inflow velocity = $u_{\text{in}}$.

  * •

Surface 3 (outlet): static pressure = 0 Pa.

  * •

Surfaces 1 and 4-6 (top and bottom walls, including the sensing surface): no
slip.

Here, no slip represents a traditional no slip boundary condition. The fluid
properties were defined by water and obtained from the FEA software’s material
library. The initial velocity and pressure fields were zero.

For transport of diluted species physics:

  * •

Surface 1: $J_{0,c}=-r_{\text{ads}}+r_{\text{des}}$.

  * •

Surface 2: $c=c_{0}$.

  * •

Surface 3: $J_{0,c}=\bm{n}\cdot D_{i}\nabla c_{i}$.

  * •

Suraces 4-6: no flux.

In the above equations, $J_{0,c}$ is the mass flux,
$r_{\text{ads}}=k_{\text{on}}cc_{s}$ is the adsorption (binding) rate of the
target analyte, $c$ is the concentration of the target analyte in solution,
$c_{s}$ is the concentration of immobilized biorecognition element,
$r_{\text{des}}=k_{\text{off}}c_{s}^{*}$ is the desorption rate, $c_{s}^{*}$
is the concentration of occupied sites, $c_{0}=c_{\text{in}}\text{GP}(t)$,
$c_{\text{in}}=1$ nM is the concentration of the target analyte in the
injected sample, $\text{GP}(t)$ is a time-dependent Gaussian pulse function
with integral normalization and integral value of unity, and no flux
represents a traditional no flux boundary condition. The following values of
transport properties and rate constants were used in the simulation:
$D_{i}=1\times 10^{-11}$ m2/s is the diffusivity of the target analyte in
solution, the binding (forward) rate constant ($k_{\text{on}})=1\times 10^{6}$
M-1s-1, and the unbinding (reverse) rate constant ($k_{\text{off}}$) $=1\times
10^{-3}$ s-1. The initial concentration field was zero.

For surface reaction chemistry:

  * •

Surface 1: reaction rate $=r_{\text{ads}}-r_{\text{des}}$.

  * •

Surfaces 2-6: no reaction.

The following surface properties were used: $\rho_{s}=$ density of surface
sites and site occupancy number $=1$. The surface diffusion of the bound
target analyte was assumed to be zero. The initial concentration of bound
target was zero.

The convection-diffusion-reaction model was discretized and solved subject to
the aforementioned boundary and initial conditions using a physics-controlled
adaptive mesh that contained 373004 (plus 351428 degrees of freedom) in the
final mesh for a given combination of inputs
$[x_{1},x_{2},x_{3}]=[u_{\text{in}},\rho_{s},\text{GP}_{\text{std}}]$, where
$\text{GP}_{\text{std}}$ is the standard deviation of the Gaussian pulse
associated with the injected sample. The specific values of
$u_{\text{in}},\rho_{s},\text{GP}_{\text{std}}$ examined were systematically
selected by traditional BO and the consensus BO model, which are referred to
as “non-collaborative” and “collaborative” learning, respectively.

The simulation’s output of interest (i.e., the quantity to be optimized) was
the maximum fractional surface coverage of the bound target analyte
($\theta_{\text{max}}$) during the transient binding response, which is a
fundamental characteristic of biosensor function that is associated with the
performance characteristics of sensitivity, detection limit, dynamic range,
and speed.

The objective of this case study was to discover (i.e., learn) the biosensor
design and measurement format parameters that maximize $\theta_{\text{max}}$
(i.e., the maximum amount of captured target analyte). An optimal solution
(i.e., design) was sought within the design space of $x_{1}\in[1\times
10^{-4},1\times 10^{-2}]$ m/s, $x_{2}\in[1.8\times 10^{-8},3.3\times 10^{-8}]$
mol/m2, and $x_{3}\in[0.5,3]$ min. The range and limits of the design space
were selected based on practical values used in previous biosensing studies
(Selmi et al.,, 2017; Squires et al.,, 2008; Baronas et al.,, 2021; Johnson
and Mutharasan,, 2012). In summary, the objective of this case study was to
learn the optimal combination of
$u_{\text{in}},\rho_{s},\text{GP}_{\text{std}}$ that maximized the amount of
captured target analyte and to compare the learning performance achieved by
traditional BO performed by a group of independent clients that do not share
information with that achieved by collaborative learning performed by a group
of clients that share information and make decisions for next-to-test values
using consensus BO.

The case study involved two groups: a non-collaborative group and a
collaborative group, each consisting of three clients. The non-collaborative
group used conventional centralized BO algorithms, where each client solves
its own problem and there is no communication among clients. In contrast, the
collaborative group utilized the CBOC algorithm (Algorithm 1) to
collaboratively find optimal solutions. For CBOC, we adopted the leader-driven
matrix mentioned in Sec. 4.3. Fig. 7 shows the trend of biosensor performance,
specifically maximum amount of captured target species, throughout the
iterative closed-loop workflow. As shown in Fig. 7, the self-driving workflow
driven by consensus BO was capable of learning an optimized biosensor design
and measurement format parameter selection after ten rounds of experimentation
as evidenced by the trend of $\theta_{\text{max}}$. Additionally, the average
biosensor performance was higher and exhibited lower variance in the
collaborative group relative to the non-collaborative group, which highlights
the value of leveraging collaboration in closed-loop high-throughput
experimentation.

Figure 7: Trend in the biosensor’s maximum fractional surface coverage (i.e.,
the maximum amount of captured target analyte) vs. iteration number (i.e.,
sequential rounds of simulation) via BO without and with collaboration after
ten iterations.

This case study serves as a proof-of-concept showcasing that collaboration in
optimal design can reap benefits for the participating entities.

##  7 Discussion and Conclusion

In this paper, we propose a collaborative Bayesian optimization framework
built upon a consensus mechanism. Our experiments on simulated data and a
real-world sensor design case study show that collaboration through our
proposed framework can accelerate and improve the optimal design process.

Collaborative BO is still in its infancy stage, and few prior studies exist
along this line. Indeed, there are many avenues of potential improvement.
These avenues include: (i) extending the consensus framework to multi-
objective and multi-fidelity settings, (ii) developing a resource-aware
approach that allows some clients to perform more or fewer experiments
depending on their resources, (iii) understanding the theoretical conditions
needed for the collaborative design to outperform its non-collaborative
counterpart. If successful, this may lead to targeted algorithms that exploit
these conditions to improve performance. We hope this research may inspire
future work along these avenues and beyond in the future.

## Appendix

##  Appendix A Technical Results

Denote by $\bm{x}_{k}^{*}=\operatorname*{argmax}_{\bm{x}}f_{k}(\bm{x})$ the
global optimal solution of the design function for client $k$. Recall that, in
Algorithm 1,
$\bm{x}^{(t)\text{new}}_{k}=[(\bm{W}^{(t)}\otimes\bm{I}_{D})\bm{x}_{\mathcal{C}}^{(t)}]_{k}$.
Let $r^{(t)}_{k}=f_{k}(\bm{x}_{k}^{*})-f_{k}(\bm{x}^{(t)\text{new}}_{k})$ be
the regret for client $k$ at iteration $t$. Our goal is to derive an upper
bound on the cumulative regret $R_{k,T}=\sum_{t=1}^{T}r^{(t)}_{k}$ and show
that $R_{k,T}$ has a sublinear growth rate. We first define a stopping
criterion for the expected improvement utility function
$\texttt{EI}_{k}^{(t)}$. By definition,
$\texttt{EI}_{k}^{(t)}(\bm{x})\coloneqq\mathbb{E}_{\hat{f}_{k}|\mathcal{D}^{(t)}_{k}}\left[U(\hat{f}_{k}(\bm{x}))\right]=\mathbb{E}_{\hat{f}_{k}|\mathcal{D}^{(t)}_{k}}\left[(\hat{f}_{k}(\bm{x})-y_{k}^{*(t)})^{+}\right]\geq
0,\forall\bm{x}\in\mathbb{R}^{D}$. Therefore, we define a small positive
constant $\kappa$ such that when $\texttt{EI}_{k}^{(t)}(\bm{x})<\kappa$, we
stop our algorithm at iteration $t$. Next, we will present the assumptions we
made for our main Theorem.

######  Assumption 1.

(Homogeneity) Design functions are homogeneous (i.e.,
$f_{1}(\bm{x})=f_{2}(\bm{x})=\cdots=f_{K}(\bm{x})$ for all
$\bm{x}\in\mathbb{R}^{D}$).

######  Assumption 2.

(Initial Sample Size) Without loss of generality, we assume the initial sample
size is 2 such that $N^{(t)}_{k}=2+t$, for all $k\in[K]$, at iteration $t$.

######  Assumption 3.

(Bounded Utility Solutions) Assume
$\left\lVert\bm{x}^{(t)}_{k}-\bm{x}^{(t)}_{j}\right\rVert_{2}\leq\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$
for all $k,j\in[K]$, where $\epsilon$ is a positive constant.

######  Assumption 4.

(Bounded Variance) Assume
$|\sigma_{k}(\bm{x}^{(t)}_{k})-\sigma_{k}(\bm{x}^{(t)\textup{new}}_{k})|\leq\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$,
for all $k\in[K]$, where $\epsilon$ is a positive constant.

Assumption 3 states that the distance between design points returned by any
two clients is bounded by $\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$
at every iteration $t$. Recall that
$\bm{x}^{(t)}_{k}=\operatorname*{argmax}_{\bm{x}}\texttt{EI}^{(t)}_{k}(\bm{x})$.
Furthermore, Assumption 4 assumes that the difference between the
$\mathcal{GP}$ posterior variance evaluated at $\bm{x}^{(t)}_{k}$ and that
evaluated at the consensus solution is bounded. These two assumptions are
intuitively understandable: at the early stages, since clients have very
little data, their surrogate functions will differ even under a homogeneous
setting. Therefore, the error term
$\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$ is large (i.e., due to
small $t$). As the collaborative process proceeds, the error term decreases,
and gradually the solutions returned by all clients will be close since the
underlying true design functions are the same. Eventually if
$\bm{x}_{1}^{(t)}=\bm{x}_{2}^{(t)}=\cdots=\bm{x}_{K}^{(t)}$, then
$\bm{W}^{(t)}\bm{x}^{(t)}_{\mathcal{C}}=\bm{x}^{(t)}_{\mathcal{C}}$, since
$\bm{x}^{(t)}_{\mathcal{C}}$ concatenates all
$\\{\bm{x}_{k}^{(t)}\\}_{k=1}^{K}$. The variance difference
$|\sigma_{k}(\bm{x}^{(t)}_{k})-\sigma_{k}(\bm{x}^{(t)\text{new}}_{k})|$ will
be zero accordingly.

One natural question is: why do we believe the error term is decreasing at a
rate of $\mathcal{O}\left(\frac{1}{(\log N^{(t)}_{k})^{0.5+\epsilon}}\right)$?
Fortunately, this argument can be justified using the following Corollary.

######  Corollary 1.

(Lederer et al.,, 2019) Suppose $f_{k}(\cdot)$ is sampled from a zero mean
$\mathcal{GP}$ defined through a Lipschitz continuous kernel function
$\mathcal{K}_{k}(\cdot,\cdot)$. Given $N_{k}$ collected observations, if there
exists an $\epsilon$ such that
$\sigma_{k}(\bm{x})\in\mathcal{O}(\frac{1}{(\log N_{k})^{0.5+\epsilon}})$,
$\forall\bm{x}\in\mathbb{R}^{D}$, then it holds for every $\delta_{0}\in(0,1)$
that

$\displaystyle\mathbb{P}\left(\sup_{\bm{x}}\left\lVert\mu_{k}(\bm{x})-f_{k}(\bm{x})\right\rVert\in\mathcal{O}\left(\frac{1}{(\log
N_{k})^{\epsilon}}-\log(\delta_{0})\frac{1}{(\log
N_{k})^{0.5+\epsilon}}\right)\right)\geq 1-\delta_{0}.$

This Lemma implies that the $\mathcal{GP}$ surrogate has an approximation
error that is decreasing at a rate of $\mathcal{O}\left(\frac{1}{(\log
N^{(t)}_{k})^{0.5+\epsilon}}\right)$. This error bound decreases as more
samples are collected (larger $N^{(t)}_{k}$).

To this end, we state our main Theorem and its proof below. All supporting
lemmas are in Sec. B.

######  Theorem.

Suppose a squared exponential kernel function
$\mathcal{K}_{k}(\bm{x}_{1},\bm{x}_{2})=u_{k}^{2}\exp\left(\frac{\left\lVert\bm{x}_{1}-\bm{x}_{2}\right\rVert^{2}}{2\ell_{k}^{2}}\right)$
is used for the $\mathcal{GP}$ surrogate and each client uses the EI utility.
Without loss of generality, assume $\mathcal{K}_{k}(\bm{x}_{1},\bm{x}_{2})\leq
1$. Suppose assumptions 1-4 hold, and the algorithm will stop when
$\texttt{EI}_{k}^{(t)}<\kappa$, then given any doubly stochastic
$\bm{W}^{(t)}$ with non-negative elements, with probability at least
$(1-\frac{\delta_{1}}{T})^{T}$, where $\delta_{1}\in(0,T)$, the cumulative
regret after $T>1$ iterations is

$\displaystyle R_{k,T}$
$\displaystyle=\sum_{t=1}^{T}r_{k}^{(t)}\leq\sqrt{\frac{6T\left[(\log
T)^{3}+1+C\right](\log T)^{D+1}}{\log(1+v_{k}^{-2})}}+\sqrt{\frac{2T(\log
T)^{D+4}}{\log(1+v_{k}^{-2})}}+\sum_{t=1}^{T}\mathcal{O}\left(\frac{1}{(\log(2+t))^{0.5+\epsilon}}\right)$
$\displaystyle\sim\mathcal{O}(\sqrt{T\times(\log T)^{D+4}}),$

where $C=\log[\frac{1}{2\pi\kappa^{2}}]$.

###### Proof.

Recall that we define the regret as
$r_{k}^{(t)}=f_{k}(\bm{x}_{k}^{*})-f_{k}(\bm{x}^{(t)\text{new}}_{k})$. We
decompose the regret $r_{k}^{(t)}$ as

$\displaystyle
r_{k}^{(t)}=\bigg{\\{}\underbrace{f_{k}(\bm{x}_{k}^{*})-y^{*(t)}_{k}}_{\text{A}}+\underbrace{y^{*(t)}_{k}-f_{k}(\bm{x}^{(t)\text{new}}_{k})}_{\text{B}}\bigg{\\}},$

where $y^{*(t)}_{k}=\max\bm{y}^{(t)}_{k}$ is the largest observed output value
from iteration $t$.

By Lemma 2, we have, with probability $1-\delta$ ($\delta\in(0,1)$),
${\texttt{EI}}^{(t)}_{k}(\bm{x})\geq
I_{k}^{(t)}(\bm{x})-\sqrt{\beta_{k}^{(t)}}\sigma_{k}^{(t)}(\bm{x})$, where
$I_{k}^{(t)}(\bm{x})=(f_{k}(\bm{x})-y^{*(t)}_{k})^{+}$ and $\beta^{(t)}_{k}$
is defined in Lemma 1. Therefore, we obtain

$\displaystyle{\texttt{EI}}^{(t)}_{k}(\bm{x}_{k}^{*})\geq
I_{k}^{(t)}(\bm{x}_{k}^{*})-\sqrt{\beta_{k}^{(t)}}\sigma_{k}^{(t)}(\bm{x}_{k}^{*}),$
$\displaystyle\Rightarrow
I_{k}^{(t)}(\bm{x}_{k}^{*})\leq{\texttt{EI}}^{(t)}_{k}(\bm{x}_{k}^{*})+\sqrt{\beta_{k}^{(t)}}\sigma_{k}^{(t)}(\bm{x}_{k}^{*}).$

If $f_{k}(\bm{x})>y^{*(t)}_{k}$, we have
$I_{k}^{(t)}(\bm{x}_{k}^{*})=f_{k}(\bm{x}_{k}^{*})-y^{*(t)}_{k}$. Otherwise,
we have $I_{k}^{(t)}(\bm{x}_{k}^{*})=0\geq
f_{k}(\bm{x}_{k}^{*})-y^{*(t)}_{k}$. Under both situations, we bound term A as

A
$\displaystyle\leq\bigg{\\{}\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{*})+\sqrt{\beta^{(t)}_{k}}\sigma^{(t)}_{k}(\bm{x}_{k}^{*})\bigg{\\}},$

with probability $1-\delta$. Next, we want to study the relationship between
$\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{*})$ and
$\texttt{EI}_{k}^{(t)}(\bm{x}^{(t)}_{k})$ to further bound term A.

##### Relationship between $\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{*})$ and
$\texttt{EI}_{k}^{(t)}(\bm{x}^{(t)}_{k})$

By definition and Lemma 2, we have

$\displaystyle\texttt{EI}_{k}^{(t)}(\bm{x}^{(t)}_{k})$
$\displaystyle=\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\phi(z_{k}^{(t)}(\bm{x}^{(t)}_{k}))+(\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})-y^{*(t)}_{k})\Phi(z_{k}^{(t)}(\bm{x}^{(t)}_{k}))$
$\displaystyle=\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\tau(z_{k}^{(t)}(\bm{x}^{(t)}_{k})),$

where $z_{k}^{(t)}(\cdot)$ and $\tau(\cdot)$ have been defined in Lemma 2.
Therefore,

$\displaystyle|\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{(t)})-\texttt{EI}^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})|$
$\displaystyle=|\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\tau(z_{k}^{(t)}(\bm{x}^{(t)}_{k}))-\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))|.$

By Assumption 4, we have
$|\sigma_{k}(\bm{x}^{(t)}_{k})-\sigma_{k}(\bm{x}^{(t)\textup{new}}_{k})|\leq\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$.
This implies

$\displaystyle\sigma_{k}(\bm{x}^{(t)}_{k})-\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})\leq\sigma_{k}(\bm{x}^{(t)\textup{new}}_{k})\leq\sigma_{k}(\bm{x}^{(t)}_{k})+\mathcal{O}(\frac{1}{(\log(2+t)^{0.5+\epsilon}}).$

Since
$\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{(t)})\geq\texttt{EI}^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})$
($\bm{x}_{k}^{(t)}$ is the maximizer of $\texttt{EI}_{k}$), we take the lower
bound of $\sigma_{k}(\bm{x}^{(t)\textup{new}}_{k})$ and get

$\displaystyle|\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{(t)})-\texttt{EI}^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})|$
$\displaystyle\leq|\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\tau(z_{k}^{(t)}(\bm{x}^{(t)}_{k}))-(\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})-\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}}))\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))|$
$\displaystyle=|\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})(\tau(z_{k}^{(t)}(\bm{x}^{(t)}_{k}))-\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})))+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))|.$

We know $\tau(\cdot)$ is Lipschitz continuous on a bounded domain. Assume the
Lipschitz constant is $L_{\tau}$, we obtain

$\displaystyle|\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{(t)})-\texttt{EI}^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})|$
$\displaystyle\leq|\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})(\tau(z_{k}^{(t)}(\bm{x}^{(t)}_{k}))-\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})))+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))|$
$\displaystyle\leq
L_{\tau}|\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})||z_{k}^{(t)}(\bm{x}^{(t)}_{k})-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})|+|\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))|$
(7)

Recall that
$z_{k}^{(t)}(\bm{x})=\frac{\mu^{(t)}_{k}(\bm{x})-y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x})}$
for any $\bm{x}\in\mathbb{R}^{D}$, we then obtain

$\displaystyle|z_{k}^{(t)}(\bm{x}^{(t)}_{k})-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})|$
$\displaystyle=\left|\frac{\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})-y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})}-\frac{\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})-y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})}\right|$
$\displaystyle=\left|\frac{\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})-\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})}-\frac{\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})-\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})}\right|$
$\displaystyle=\left|\frac{\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})-\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})y^{*(t)}_{k}-\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})+\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})}\right|.$

Now, to bound
$|z_{k}^{(t)}(\bm{x}^{(t)}_{k})-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})|$, we
need to study the properties of $\mu_{k}$ and $\sigma_{k}$.

Lederer et al., (2019) show that $\mathcal{GP}$ with a squared exponential
kernel has a Lipschitz continuous posterior mean function. Assume the
Lipschitz constant is $L_{\mu}$. Therefore, we have
$|\mu^{(t)}_{k}(\bm{x}_{k}^{(t)\text{new}})-\mu^{(t)}_{k}(\bm{x}_{k}^{(t)})|\leq
L_{\mu}\left\lVert\bm{x}_{k}^{(t)\text{new}}-\bm{x}_{k}^{(t)}\right\rVert_{2}$.
By Assumption 3, we know
$\left\lVert\bm{x}_{k}^{(t)}-\bm{x}_{j}^{(t)}\right\rVert_{2}\leq\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$.
Therefore, we can show that

$\displaystyle\left\lVert\bm{x}_{k}^{(t)\text{new}}-\bm{x}_{k}^{(t)}\right\rVert_{2}$
$\displaystyle=\left\lVert\sum_{j=1}^{K}w^{(t)}_{kj}\bm{x}_{j}^{(t)}-\bm{x}_{k}^{(t)}\right\rVert_{2}$
$\displaystyle=\left\lVert\sum_{j=1}^{K}w^{(t)}_{kj}\bm{x}_{j}^{(t)}-\sum_{j=1}^{K}w^{(t)}_{kj}\bm{x}_{k}^{(t)}\right\rVert_{2}$
Since $\bm{W}^{(t)}$ is doubly stochastic
$\displaystyle\leq\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}}).$

Therefore, we know
$|\mu^{(t)}_{k}(\bm{x}_{k}^{(t)\text{new}})-\mu^{(t)}_{k}(\bm{x}_{k}^{(t)})|$
is of order $\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}}).$ Without loss
of generality, we assume $\mu_{k}^{(t)}\geq 0$ and the negative case follows a
similar proof. By Assumption 4, we know
$|\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})-\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})|\leq\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$.
Therefore $\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})$ is upper bounded by
$\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})+\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})$.
When
$\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})\geq\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})$,
we take the upper bound of $\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})$ and
derive that

$\displaystyle\left|\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})-\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\right|$
$\displaystyle=\left|(\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}}))\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})-\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\right|$
$\displaystyle=\left|\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})(\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})-\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k}))+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})\right|.$

By Assumption 4, we know $\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})$ is finite.
Otherwise, the bound on
$|\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})-\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})|$
is infinite. Also, $\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\neq 0$. Therefore, we
know
$\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})(\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})-\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k}))$
is of order $\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$. By Corollary
1, we know $\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})$ is finite. As a result,
$\left|\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})-\sigma_{k}^{(t)}(\bm{x}^{(t)}_{k})\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\right|$
is of order $\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$.

When
$\mu^{(t)}_{k}(\bm{x}^{(t)}_{k})\leq\mu^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})$,
we take the lower bound of $\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})$ and
the same conclusion holds.

As a result, we show that
$|z_{k}^{(t)}(\bm{x}^{(t)}_{k})-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})|$ is
of order $\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$. Therefore, Eq.
(A) gives
$\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{(t)})\leq\texttt{EI}^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$.
Additionally, we know
$\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{*})\leq\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{(t)})$
since $\bm{x}_{k}^{(t)}$ is the maximizer of $\texttt{EI}_{k}$. This gives
$\texttt{EI}^{(t)}_{k}(\bm{x}_{k}^{*})\leq\texttt{EI}^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})$.
Therefore, with probability $1-\delta$,

A
$\displaystyle\leq\texttt{EI}^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}})+\sqrt{\beta^{(t)}_{k}}\sigma^{(t)}_{k}(\bm{x}_{k}^{*})$
$\displaystyle=\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))+\sqrt{\beta^{(t)}_{k}}\sigma^{(t)}_{k}(\bm{x}_{k}^{*})+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}}),$

where the last equality holds by Lemma 2.

##### Bounding Term B

We first decompose B as

B
$\displaystyle=y^{*(t)}_{k}-f_{k}(\bm{x}^{(t)\text{new}}_{k})=y^{*(t)}_{k}-\mu_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})+\mu_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})-f_{k}(\bm{x}^{(t)\text{new}}_{k}).$

Using Lemma 1, we know
$\mu_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})-f_{k}(\bm{x}^{(t)\text{new}}_{k})\leq\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\sqrt{\beta_{k}^{(t)}}$.
Therefore, we can bound term B as

B
$\displaystyle=y^{*(t)}_{k}-f_{k}(\bm{x}^{(t)\text{new}}_{k})=y^{*(t)}_{k}-\mu_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})+\mu_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})-f_{k}(\bm{x}^{(t)\text{new}}_{k})$
$\displaystyle\leq\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})(-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))+\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\sqrt{\beta_{k}^{(t)}}$
$\displaystyle\leq\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\left(\tau(-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))-\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))+\sqrt{\beta_{k}^{(t)}}\right),$

where the last inequality uses the fact that $z=\tau(z)-\tau(-z)$.

##### Compiling A and B

Given all the previous information, we reach

$\displaystyle r_{k}^{(t)}$ $\displaystyle=\text{A}+\text{B}$
$\displaystyle\leq\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\left(\tau(-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))+\sqrt{\beta_{k}^{(t)}}\right)+\sqrt{\beta^{(t)}_{k}}\sigma^{(t)}_{k}(\bm{x}_{k}^{*})+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}}).$

Now, by Lemma 3, we have, with probability $1-\delta$,

$\displaystyle
r_{k}^{(t)}\leq\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\left[\sqrt{\beta_{k}^{(t)}}+1+\sqrt{C}\right]+\sqrt{\beta^{(t)}_{k}}\sigma^{(t)}_{k}(\bm{x}_{k}^{*})+\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}}),$

where $C$ is defined in Lemma 3. We next aim to further bound
$\sum_{t=1}^{T}r_{k}^{(t)}$ using Lemma 4.

Let $\delta=\frac{\delta_{1}}{T}$, where $\delta_{1}\in(0,T)$. By Lemma 4, and
using the fact that $\left[\sqrt{\beta_{k}^{(t)}}+1+\sqrt{C}\right]^{2}\geq
3(\beta_{k}^{(t)}+1+C)$, we have

$\displaystyle\sum_{t=1}^{T}\left(\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\left[\sqrt{\beta_{k}^{(t)}}+1+\sqrt{C}\right]\right)^{2}$
$\displaystyle\leq
3\sum_{t=1}^{T}(\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k}))^{2}\left[\beta_{k}^{(t)}+1+C\right]$
$\displaystyle\leq
3\sum_{t=1}^{T}(\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k}))^{2}\left[\beta_{k}^{(T)}+1+C\right]$
$\beta_{k}^{(t)}$ is non-decreasing
$\displaystyle\leq\frac{6\left[\beta_{k}^{(T)}+1+C\right]\gamma_{k,T}}{\log(1+v_{k}^{-2})}$
$\displaystyle\text{by Lemma \ref{lemma:5}},$

with probability $(1-\frac{\delta_{1}}{T})^{T}$, where $\gamma_{k,T}$ is
defined in Lemma 4. Therefore, by the Cauchy-Schwartz inequality, we have

$\displaystyle\sum_{t=1}^{T}\left(\sigma^{(t)}_{k}(\bm{x}^{(t)new}_{k})\left[\sqrt{\beta_{k}^{(t)}}+1+\sqrt{C}\right]\right)$
$\displaystyle\leq\sqrt{T}\sqrt{\sum_{t=1}^{T}\left(\sigma^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\left[\sqrt{\beta_{k}^{(t)}}+1+\sqrt{C}\right]\right)^{2}}$
$\displaystyle\leq\sqrt{T}\sqrt{\frac{6\left[\beta_{k}^{(T)}+1+C\right]\gamma_{k,T}}{\log(1+v_{k}^{-2})}}.$
by Lemma 4

Similarly, we can show that, with probability $(1-\frac{\delta_{1}}{T})^{T}$,

$\displaystyle\sum_{t=1}^{T}(\sqrt{\beta^{(t)}_{k}}\sigma^{(t)}_{k}(\bm{x}_{k}^{(t)\text{new}}))\leq\sqrt{T}\sqrt{\beta^{(T)}_{k}\sum_{t=1}^{T}\sigma^{2(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})}\leq\sqrt{\frac{2T\beta_{k}^{(T)}\gamma_{k,T}}{\log(1+v_{k}^{-2})}}.$

Finally, we obtain, with probability $(1-\frac{\delta_{1}}{T})^{T}$,

$\displaystyle\sum_{t=1}^{T}r_{k}^{(t)}\leq\sqrt{T}\sqrt{\frac{6\left[\beta_{k}^{(T)}+1+C\right]\gamma_{k,T}}{\log(1+v_{k}^{-2})}}+\sqrt{\frac{2T\beta_{k}^{(T)}\gamma_{k,T}}{\log(1+v_{k}^{-2})}}+\sum_{t=1}^{T}\mathcal{O}(\frac{1}{(\log(2+t))^{0.5+\epsilon}}).$

From Lemma 1, we know $\beta_{k}^{(T)}\sim\mathcal{O}((\log(T/\delta))^{3})$,
where $\delta=\frac{\delta_{1}}{T}$, $\forall\delta_{1}\in(0,T)$. Furthermore,
Contal et al., (2014) have shown that $\gamma_{k,T}\sim\mathcal{O}((\log
T)^{D+1})$ if a squared exponential kernel is used. Therefore, we can show
that $\sum_{t=1}^{T}r_{k}^{(t)}\sim\mathcal{O}(\sqrt{T\times(\log T)^{D+4}})$.
As a result, $\lim_{T\to\infty}\frac{R_{k,T}}{T}=0$.

∎

##  Appendix B Supporting Lemmas and Some Useful Equations

Before proving Lemma 2, we first need a Lemma from Srinivas et al., (2009).

######  Lemma 1.

(Srinivas et al.,, 2009) Let $\delta\in(0,1)$. Then, $\forall t$ and
$\forall\bm{x}\in\mathbb{R}^{D}$, there exists a non-decreasing positive
sequence $\\{\beta_{k}^{(t)}\\}_{t}$ such that
$\beta_{k}^{(t)}\sim\mathcal{O}((\log\frac{t}{\delta})^{3})$ and

$\displaystyle\mathbb{P}\left(|\mu^{(t)}_{k}(\bm{x})-f_{k}(\bm{x})|\leq\sqrt{\beta_{k}^{(t)}}\sigma_{k}^{(t)}(\bm{x})\right)\geq
1-\delta.$

######  Lemma 2.

Let $I_{k}^{(t)}(\bm{x})=(f_{k}(\bm{x})-y^{*(t)}_{k})^{+}$. We have
${\texttt{EI}}^{(t)}_{k}(\bm{x})=\sigma_{k}^{(t)}(\bm{x})\tau(z_{k}^{(t)}(\bm{x}))$,
where $\tau(z_{k}^{(t)}(\bm{x}))\coloneqq
z_{k}^{(t)}(\bm{x})\Phi(z_{k}^{(t)}(\bm{x}))+\phi(z_{k}^{(t)}(\bm{x}))$ and
$z_{k}^{(t)}(\bm{x})=\frac{\mu^{(t)}_{k}(\bm{x})-y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x})}$.
Additionally, with probability $1-\delta$ ($\delta\in(0,1)$),
${\texttt{EI}}^{(t)}_{k}(\bm{x})\geq
I_{k}^{(t)}(\bm{x})-\sqrt{\beta_{k}^{(t)}}\sigma_{k}^{(t)}(\bm{x})$.

###### Proof.

We have, with probability at least $1-\delta$,

$\displaystyle\texttt{EI}^{(t)}_{k}(\bm{x})$
$\displaystyle=\sigma_{k}^{(t)}(\bm{x})\phi(z_{k}^{(t)}(\bm{x}))+(\mu^{(t)}_{k}(\bm{x})-y^{*(t)}_{k})\Phi(z_{k}^{(t)}(\bm{x}))$
$\displaystyle=\sigma_{k}^{(t)}(\bm{x})\left[z_{k}^{(t)}(\bm{x})\Phi(z_{k}^{(t)}(\bm{x}))+\phi(z_{k}^{(t)}(\bm{x}))\right]$
$\displaystyle=\sigma_{k}^{(t)}(\bm{x})\tau(z_{k}^{(t)}(\bm{x})).$

Therefore, we have shown that
${\texttt{EI}}^{(t)}_{k}(\bm{x})=\sigma_{k}^{(t)}(\bm{x})\tau(z_{k}^{(t)}(\bm{x}))$.

Now, let
$q_{k}^{(t)}(\bm{x})=\frac{f_{k}(\bm{x})-y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x})}$,
we know

$\displaystyle z_{k}^{(t)}(\bm{x})-q_{k}^{(t)}(\bm{x})$
$\displaystyle=\frac{\mu^{(t)}_{k}(\bm{x})-y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x})}-\frac{f_{k}(\bm{x})-y^{*(t)}_{k}}{\sigma_{k}^{(t)}(\bm{x})}$
$\displaystyle=\frac{\mu^{(t)}_{k}(\bm{x})-f_{k}(\bm{x})}{\sigma_{k}^{(t)}(\bm{x})}$
$\displaystyle\geq-\frac{\sqrt{\beta_{k}^{(t)}}\sigma_{k}^{(t)}(\bm{x})}{\sigma_{k}^{(t)}(\bm{x})}.$
By Lemma 1

Therefore, $z_{k}^{(t)}(\bm{x})\geq
q_{k}^{(t)}(\bm{x})-\sqrt{\beta_{k}^{(t)}}$. We know $\tau(\cdot)$ is a non-
decreasing function. Therefore, we have

$\displaystyle\texttt{EI}^{(t)}_{k}(\bm{x})$
$\displaystyle=\sigma_{k}^{(t)}(\bm{x})\tau(z_{k}^{(t)}(\bm{x}))$
$\displaystyle\geq\sigma_{k}^{(t)}(\bm{x})\tau\left(q_{k}^{(t)}(\bm{x})-\sqrt{\beta_{k}^{(t)}}\right)$
$\displaystyle\geq\sigma_{k}^{(t)}(\bm{x})\left(q_{k}^{(t)}(\bm{x})-\sqrt{\beta_{k}^{(t)}}\right).$

The last inequality holds due to $\tau(z)\geq z,\forall z$. If
$I_{k}^{(t)}(\bm{x})=0,$ then we know $q_{k}^{(t)}(\bm{x})\leq 0$ and
$\sigma_{k}^{(t)}(\bm{x})\left(q_{k}^{(t)}(\bm{x})-\sqrt{\beta_{k}^{(t)}}\right)<0$.
Therefore, we have $\texttt{EI}^{(t)}_{k}(\bm{x})\geq
I_{k}^{(t)}(\bm{x})-\sqrt{\beta_{k}^{(t)}}\sigma_{k}^{(t)}(\bm{x})$. If
$I_{k}^{(t)}(\bm{x})>0$, then we also have

$\displaystyle\texttt{EI}^{(t)}_{k}(\bm{x})$
$\displaystyle\geq\sigma_{k}^{(t)}(\bm{x})\left(q_{k}^{(t)}(\bm{x})-\sqrt{\beta_{k}^{(t)}}\right)$
$\displaystyle=I_{k}^{(t)}(\bm{x})-\sqrt{\beta_{k}^{(t)}}\sigma_{k}^{(t)}(\bm{x}).$

∎

######  Lemma 3.

Given a constant $\kappa>0$, if
$y^{*(t)}_{k}-\mu_{k}^{(t)}(\bm{x}^{(t)\textup{new}}_{k})\geq 0$, then we have
$y^{*(t)}_{k}-\mu_{k}^{(t)}(\bm{x}^{(t)\textup{new}}_{k})\leq\sigma^{(t)}_{k}(\bm{x}^{(t)\textup{new}}_{k})\sqrt{C}$,
where $C=\log[\frac{1}{2\pi\kappa^{2}}]$. Furthermore, We have
$\tau(-z_{k}^{(t)}(\bm{x}^{(t)new}_{k}))\leq 1+\sqrt{C}$.

###### Proof.

Suppose the algorithm is not stopped. By Lemma 2, we know

$\displaystyle\kappa\leq\texttt{EI}^{(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})=\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\tau(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})).$

When $y^{*(t)}_{k}-\mu_{k}^{(t)}(\bm{x}^{(t)\textup{new}}_{k})\geq 0$ (i.e.,
$z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\leq 0$), we have

$\displaystyle\tau(z_{k}^{(t)}(\bm{x}_{k}^{(t)\text{new}}))\leq\phi(z_{k}^{(t)}(\bm{x}_{k}^{(t)\text{new}})).$

This implies

$\displaystyle\kappa$
$\displaystyle\leq\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\frac{1}{\sqrt{2\pi}}\exp\left[-\frac{1}{2}(z_{k}^{(t)}(\bm{x}_{k}^{(t)\text{new}}))^{2}\right]$
$\displaystyle\Rightarrow(z_{k}^{(t)}(\bm{x}_{k}^{(t)\text{new}}))^{2}\leq
2\log\frac{\sigma_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})}{\sqrt{2\pi}\kappa}\leq\log[\frac{1}{2\pi\kappa^{2}}]=C.$
Since $\mathcal{K}_{k}(\bm{x}_{1},\bm{x}_{2})\leq 1$

Therefore, $(z_{k}^{(t)}(\bm{x}^{(t)\text{new}}))^{2}\leq C$. By the
definition of $z_{k}^{(t)}(\bm{x}^{(t)\text{new}})$, we have
$y^{*(t)}_{k}-\mu_{k}^{(t)}(\bm{x}^{(t)\textup{new}}_{k})\leq\sigma^{(t)}_{k}(\bm{x}^{(t)\textup{new}}_{k})\sqrt{C}$.

Now, to prove $\tau(-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))\leq 1+\sqrt{C}$,
we consider two situations. When
$y^{*(t)}_{k}-\mu_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\geq 0$, we have

$\displaystyle\tau(-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))\leq 1+\sqrt{C}$

using the fact that $\tau(z)\leq 1+z,\forall z\geq 0$. When
$y^{*(t)}_{k}-\mu_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k})\leq 0$, we know
$\forall z\leq 0,\tau(z)\leq\phi(z)\leq 1$. Therefore,

$\displaystyle\tau(-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))\leq 1.$

As a result, we conclude that
$\tau(-z_{k}^{(t)}(\bm{x}^{(t)\text{new}}_{k}))\leq 1+\sqrt{C}$. ∎

######  Lemma 4.

Let $\text{J}^{(t)}_{k}$ be the mutual information at iteration $t$ and denote
by $\gamma_{k,T}=\max\text{J}^{(t)}_{k}$ the maximum information gain after
$T$ iterations. We have
$\sum_{t=1}^{T}\sigma^{2(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\leq\frac{2}{\log(1+\sigma_{k}^{-2})}\gamma_{k,T}$.

###### Proof.

By definition,
$\text{J}^{(t)}_{k}=\frac{1}{2}\log(1+v_{k}^{-2}\sigma^{2(t)}_{k}(\bm{x}^{(t)\text{new}}_{k}))$
(Contal et al.,, 2014). Given a constant
$s^{2(t)}=v_{k}^{-2}\sigma^{2(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})\leq
v_{k}^{-2}$, we have
$1\leq\frac{s^{(t)}}{\log(1+s^{(t)})}\leq\frac{1}{v_{k}^{2}\log(1+v_{k}^{-2})}$.

Therefore,

$\displaystyle\sum_{t=1}^{T}\sigma^{2(t)}_{k}(\bm{x}^{(t)\text{new}}_{k})$
$\displaystyle=\sum_{t=1}^{T}v_{k}^{2}s^{(t)}$
$\displaystyle\leq\sum_{t=1}^{T}v_{k}^{2}\left[\frac{\log(1+s^{(t)})}{v_{k}^{2}\log(1+v_{k}^{-2})}\right]$
$\displaystyle=\frac{2}{\log(1+v_{k}^{-2})}\frac{1}{2}\sum_{t=1}^{T}\log(1+v_{k}^{-2}\sigma^{2(t)}_{k}(\bm{x}^{(t)\text{new}}_{k}))$
$\displaystyle=\frac{2}{\log(1+v_{k}^{-2})}\text{J}^{(t)}_{k}$
$\displaystyle\leq\frac{2}{\log(1+v_{k}^{-2})}\gamma_{k,T}.$

∎

##  Appendix C Details on Testing Functions

In this section, we provide more details on the testing function.

##### Shekel-$10$

Shekel-$10$ is a four-dimensional ($D=4$) with $10$ local minima. It has a
functional form

$\displaystyle
f(\bm{x})=-\sum_{i=1}^{10}\left(\sum_{d=1}^{4}(\tensor[^{d}]{x}{}-F_{di})^{2}+\xi_{i}\right)^{-1},$

where $F_{di}$ is the $(d,i)$-th component of a matrix $\bm{F}$ and $\xi_{i}$
is the $i$-th component of a vector $\bm{\xi}$, and

$\displaystyle\bm{\xi}=\frac{1}{10}(1,2,2,4,4,6,3,7,5,5)^{\intercal},\bm{F}=\begin{bmatrix}4&1&8&6&3&2&5&8&6&7\\\
4&1&8&6&7&9&3&1&2&3.6\\\ 4&1&8&6&3&2&5&8&6&7\\\
4&1&8&6&7&9&3&1&2&3.6\end{bmatrix},$

$\bm{x}=(\tensor[^{1}]{x}{},\ldots,\tensor[^{d}]{x}{}\ldots,\tensor[^{10}]{x}{})^{\intercal}\in[0,10]^{D}$.

##### Harmann-$6$

Hartmann-$6$ is a six-dimensional ($D=6$) function in a form of

$\displaystyle
f(\bm{x})=\sum_{i=1}^{4}\alpha_{i}\exp\left(-\sum_{d=1}^{6}A_{id}(\tensor[^{d}]{x}{}-P_{id})^{2}\right),$

where

$\displaystyle A$ $\displaystyle=\begin{bmatrix}10&3&17&3.5&1.7&8\\\
0.05&10&17&0.1&8&14\\\ 3&3.5&1.7&10&17&8\\\
17&8&0.05&10&0.1&14\end{bmatrix},P=\begin{bmatrix}0.1312&0.1696&0.5569&0.0124&0.8283&0.5886\\\
0.2329&0.4135&0.8307&0.3736&0.1004&0.9991\\\
0.2348&0.1451&0.3522&0.2883&0.3047&0.6650\\\
0.4047&0.8828&0.8732&0.5743&0.1091&0.0381\end{bmatrix},$
$\displaystyle\alpha_{1}=1.0,\alpha_{2}=1.2,\alpha_{3}=3.0,\alpha_{4}=3.2,\bm{x}\in[0,1]^{6}.$

For each client, we generate $a_{1}\sim\texttt{Uniform}(0.5,2)$ and
$a_{2},a_{3}\sim\mathcal{N}(0,1)$.

##### Branin

Branin is a two-dimensional ($D=2$) function that has the form

$\displaystyle
f(\bm{x})=(\tensor[^{2}]{x}{}-\frac{5.1}{4\pi^{2}}\tensor[^{1}]{x}{}^{2}+\frac{5}{\pi}\tensor[^{1}]{x}{}-6)^{2}+10(1-\frac{1}{8\pi})\cos(\tensor[^{1}]{x}{})+10,$

where $\bm{x}=(\tensor[^{1}]{x}{},\tensor[^{2}]{x}{})^{\intercal}$ and
$\tensor[^{1}]{x}{}\in[-5,10],\tensor[^{2}]{x}{}\in[0,15]$. For each client,
we set $a_{1}\sim\texttt{Uniform}(0.5,1)$ and
$a_{2},a_{3}\sim\mathcal{N}(0,1)$.

##### Ackley-$D$

Ackley-$D$ is a $D$-dimensional function. When $D=2$, this function is
characterized by a nearly flat outer region and a large valley at the center.
This property renders the Ackley function very challenging to be optimized.
Given an arbitrary dimension $D$, Ackley-$D$ has the following function form:

$\displaystyle
f(\bm{x})=-20\exp\left(-0.2\sqrt{\frac{1}{D}\sum_{d=1}^{D}\tensor[^{d}]{x}{}^{2}}\right)-\exp\left(\frac{1}{D}\sum_{i=1}^{D}\cos(2\pi\tensor[^{d}]{x}{})\right)+20+\exp(1),$

where $\tensor[^{d}]{x}{}\in[-32.768,32.768],$ for $d=1,\ldots,D$. For each
client, we set $a_{1}\sim\texttt{Uniform}(1,2)$ and
$a_{2},a_{3}\sim\mathcal{N}(0.5,1)$.

## References

  * Amo-Salas et al., (2016) Amo-Salas, M., Delgado-Márquez, E., and López-Fidalgo, J. (2016).  Optimal experimental designs in the flow rate of particles.  Technometrics, 58(2):269–276.
  * Antony, (2014) Antony, J. (2014).  Design of experiments for engineers and scientists.  Elsevier.
  * Azimi et al., (2010) Azimi, J., Fern, A., and Fern, X. (2010).  Batch bayesian optimization via simulation matching.  Advances in Neural Information Processing Systems, 23.
  * Azimi et al., (2012) Azimi, J., Jalali, A., and Fern, X. (2012).  Hybrid batch bayesian optimization.  arXiv preprint arXiv:1202.5597.
  * Baronas et al., (2021) Baronas, R., Ivanauskas, F., Kulys, J., Baronas, R., Ivanauskas, F., and Kulys, J. (2021).  Biosensors acting in injection mode.  Mathematical Modeling of Biosensors, pages 183–205.
  * Bellman, (1956) Bellman, R. (1956).  A problem in the sequential design of experiments.  Sankhyā: The Indian Journal of Statistics (1933-1960), 16(3/4):221–229.
  * Berahas et al., (2018) Berahas, A. S., Bollapragada, R., Keskar, N. S., and Wei, E. (2018).  Balancing communication and computation in distributed optimization.  IEEE Transactions on Automatic Control, 64(8):3141–3155.
  * Box and Wilson, (1951) Box, G. and Wilson, K. (1951).  On the experimental designs for exploring response surfaces.  Ann Math Stat, 13:1–45.
  * Buchanan and Gardner, (2019) Buchanan, C. and Gardner, L. (2019).  Metal 3d printing in construction: A review of methods, research, applications, opportunities and challenges.  Engineering Structures, 180:332–348.
  * Carpenter et al., (2018) Carpenter, A. C., Paulsen, I. T., and Williams, T. C. (2018).  Blueprints for biosensors: design, limitations, and applications.  Genes, 9(8):375.
  * Cesewski and Johnson, (2020) Cesewski, E. and Johnson, B. N. (2020).  Electrochemical biosensors for pathogen detection.  Biosensors and Bioelectronics, 159:112214.
  * Chernoff, (1959) Chernoff, H. (1959).  Sequential design of experiments.  The Annals of Mathematical Statistics, 30(3):755–770.
  * Contal et al., (2014) Contal, E., Perchet, V., and Vayatis, N. (2014).  Gaussian process optimization with mutual information.  In International Conference on Machine Learning, pages 253–261. PMLR.
  * Currin et al., (1991) Currin, C., Mitchell, T., Morris, M., and Ylvisaker, D. (1991).  Bayesian prediction of deterministic functions, with applications to the design and analysis of computer experiments.  Journal of the American Statistical Association, 86(416):953–963.
  * Dai et al., (2020) Dai, Z., Low, B. K. H., and Jaillet, P. (2020).  Federated bayesian optimization via thompson sampling.  Advances in Neural Information Processing Systems, 33:9687–9699.
  * Dai et al., (2021) Dai, Z., Low, B. K. H., and Jaillet, P. (2021).  Differentially private federated bayesian optimization with distributed exploration.  Advances in Neural Information Processing Systems, 34:9125–9139.
  * Duan et al., (2017) Duan, W., Ankenman, B. E., Sanchez, S. M., and Sanchez, P. J. (2017).  Sliced full factorial-based latin hypercube designs as a framework for a batch sequential design algorithm.  Technometrics, 59(1):11–22.
  * Ekström et al., (2019) Ekström, A., Forssén, C., Dimitrakakis, C., Dubhashi, D., Johansson, H. T., Muhammad, A. S., Salomonsson, H., and Schliep, A. (2019).  Bayesian optimization in ab initio nuclear physics.  Journal of Physics G: Nuclear and Particle Physics, 46(9):095101.
  * Ezzat et al., (2018) Ezzat, A. A., Pourhabib, A., and Ding, Y. (2018).  Sequential design for functional calibration of computer models.  Technometrics, 60(3):286–296.
  * Friedman and Savage, (1947) Friedman, M. and Savage, L. J. (1947).  Planning experiments seeking maxima.  Techniques of statistical analysis, pages 365–372.
  * Gramacy and Lee, (2009) Gramacy, R. B. and Lee, H. K. (2009).  Adaptive design and analysis of supercomputer experiments.  Technometrics, 51(2):130–145.
  * He et al., (2017) He, X., Tuo, R., and Wu, C. J. (2017).  Optimization of multi-fidelity computer experiments via the eqie criterion.  Technometrics, 59(1):58–68.
  * Huang et al., (2021) Huang, C., Ren, Y., McGuinness, E. K., Losego, M. D., Lively, R. P., and Joseph, V. R. (2021).  Bayesian optimization of functional output in inverse problems.  Optimization and Engineering, 22:2553–2574.
  * Hunt, (2020) Hunt, N. (2020).  Batch Bayesian optimization.  PhD thesis, Massachusetts Institute of Technology.
  * Jiang et al., (2020) Jiang, S., Chai, H., Gonzalez, J., and Garnett, R. (2020).  Binoculars for efficient, nonmyopic sequential experimental design.  In International Conference on Machine Learning, pages 4794–4803. PMLR.
  * Johnson and Mutharasan, (2012) Johnson, B. N. and Mutharasan, R. (2012).  Sample preparation-free, real-time detection of microrna in human serum using piezoelectric cantilever biosensors at attomole level.  Analytical chemistry, 84(23):10426–10436.
  * Johnson and Mutharasan, (2014) Johnson, B. N. and Mutharasan, R. (2014).  Biosensor-based microrna detection: techniques, design, performance, and challenges.  Analyst, 139(7):1576–1588.
  * Jones et al., (1998) 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.
  * Kleijnen, (2018) Kleijnen, J. P. (2018).  Design and analysis of simulation experiments.  Springer.
  * Konomi et al., (2014) Konomi, B., Karagiannis, G., Sarkar, A., Sun, X., and Lin, G. (2014).  Bayesian treed multivariate gaussian process with adaptive design: Application to a carbon capture unit.  Technometrics, 56(2):145–158.
  * Kontar et al., (2021) Kontar, R., Shi, N., Yue, X., Chung, S., Byon, E., Chowdhury, M., Jin, J., Kontar, W., Masoud, N., Nouiehed, M., et al. (2021).  The internet of federated things (ioft).  IEEE Access, 9:156071–156113.
  * Kontar et al., (2018) Kontar, R., Zhou, S., Sankavaram, C., Du, X., and Zhang, Y. (2018).  Nonparametric modeling and prognosis of condition monitoring signals using multivariate gaussian convolution processes.  Technometrics, 60(4):484–496.
  * Koveal et al., (2022) Koveal, D., Rosen, P. C., Meyer, D. J., Díaz-García, C. M., Wang, Y., Cai, L.-H., Chou, P. J., Weitz, D. A., and Yellen, G. (2022).  A high-throughput multiparameter screen for accelerated development and optimization of soluble genetically encoded fluorescent biosensors.  Nature Communications, 13(1):2919.
  * Kushner, (1964) Kushner, H. J. (1964).  A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise.
  * Kusne et al., (2020) Kusne, A. G., Yu, H., Wu, C., Zhang, H., Hattrick-Simpers, J., DeCost, B., Sarker, S., Oses, C., Toher, C., Curtarolo, S., et al. (2020).  On-the-fly closed-loop materials discovery via bayesian active learning.  Nature communications, 11(1):5966.
  * Le Gratiet and Cannamela, (2015) Le Gratiet, L. and Cannamela, C. (2015).  Cokriging-based sequential design strategies using fast cross-validation techniques for multi-fidelity computer codes.  Technometrics, 57(3):418–427.
  * Lederer et al., (2019) Lederer, A., Umlauft, J., and Hirche, S. (2019).  Uniform error bounds for gaussian process regression with application to safe control.  Advances in Neural Information Processing Systems, 32.
  * Lee et al., (2020) Lee, E., Eriksson, D., Bindel, D., Cheng, B., and Mccourt, M. (2020).  Efficient rollout strategies for bayesian optimization.  In Conference on Uncertainty in Artificial Intelligence, pages 260–269. PMLR.
  * Lee et al., (2018) Lee, I.-C., Hong, Y., Tseng, S.-T., and Dasgupta, T. (2018).  Sequential bayesian design for accelerated life tests.  Technometrics, 60(4):472–483.
  * Marmin and Filippone, (2022) Marmin, S. and Filippone, M. (2022).  Deep gaussian processes for calibration of computer models (with discussion).  Bayesian Analysis, 17(4):1301–1350.
  * Marshall et al., (1979) Marshall, A. W., Olkin, I., and Arnold, B. C. (1979).  Inequalities: theory of majorization and its applications.
  * Ming et al., (2022) Ming, D., Williamson, D., and Guillas, S. (2022).  Deep gaussian process emulation using stochastic imputation.  Technometrics, pages 1–12.
  * Nedic and Ozdaglar, (2009) Nedic, A. and Ozdaglar, A. (2009).  Distributed subgradient methods for multi-agent optimization.  IEEE Transactions on Automatic Control, 54(1):48–61.
  * Overstall and Woods, (2017) Overstall, A. M. and Woods, D. C. (2017).  Bayesian design of experiments using approximate coordinate exchange.  Technometrics, 59(4):458–470.
  * Polatoğlu et al., (2020) Polatoğlu, İ., Aydın, L., Nevruz, B. Ç., and Özer, S. (2020).  A novel approach for the optimal design of a biosensor.  Analytical Letters, 53(9):1428–1445.
  * Rahman et al., (2021) Rahman, M. H., Xie, C., and Sha, Z. (2021).  Predicting sequential design decisions using the function-behavior-structure design process model and recurrent neural networks.  Journal of Mechanical Design, 143(8).
  * Robbins and Monro, (1951) Robbins, H. and Monro, S. (1951).  A stochastic approximation method.  The annals of mathematical statistics, pages 400–407.
  * Sacks et al., (1989) Sacks, J., Welch, W. J., Mitchell, T. J., and Wynn, H. P. (1989).  Design and analysis of computer experiments.  Statistical science, 4(4):409–423.
  * Sacks and Ylvisaker, (1970) Sacks, J. and Ylvisaker, D. (1970).  Statistical designs and integral approximation.  Proc. 12th Bienn. Semin. Can. Math. Congr.(R. Pyke, ed.), pages 115–136.
  * Sauer et al., (2022) Sauer, A., Gramacy, R. B., and Higdon, D. (2022).  Active learning for deep gaussian process surrogates.  Technometrics, pages 1–15.
  * Selmi et al., (2017) Selmi, M., Gazzah, M. H., and Belmabrouk, H. (2017).  Optimization of microfluidic biosensor efficiency by means of fluid flow engineering.  Scientific Reports, 7(1):5721.
  * Shahriari et al., (2015) Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., and De Freitas, N. (2015).  Taking the human out of the loop: A review of bayesian optimization.  Proceedings of the IEEE, 104(1):148–175.
  * Shi et al., (2017) Shi, X., Cao, J., and Huang, W. (2017).  Distributed parametric consensus optimization with an application to model predictive consensus problem.  IEEE Transactions on Cybernetics, 48(7):2024–2035.
  * Snoek et al., (2015) Snoek, J., Rippel, O., Swersky, K., Kiros, R., Satish, N., Sundaram, N., Patwary, M., Prabhat, M., and Adams, R. (2015).  Scalable bayesian optimization using deep neural networks.  In International conference on machine learning, pages 2171–2180. PMLR.
  * Springenberg et al., (2016) Springenberg, J. T., Klein, A., Falkner, S., and Hutter, F. (2016).  Bayesian optimization with robust bayesian neural networks.  Advances in neural information processing systems, 29.
  * Squires et al., (2008) Squires, T. M., Messinger, R. J., and Manalis, S. R. (2008).  Making it stick: convection, reaction and diffusion in surface-based biosensors.  Nature biotechnology, 26(4):417–426.
  * Srinivas et al., (2009) Srinivas, N., Krause, A., Kakade, S. M., and Seeger, M. (2009).  Gaussian process optimization in the bandit setting: No regret and experimental design.  arXiv preprint arXiv:0912.3995.
  * Stroh et al., (2022) Stroh, R., Bect, J., Demeyer, S., Fischer, N., Marquis, D., and Vazquez, E. (2022).  Sequential design of multi-fidelity computer experiments: maximizing the rate of stepwise uncertainty reduction.  Technometrics, 64(2):199–209.
  * Surjanovic and Bingham, (2013) Surjanovic, S. and Bingham, D. (2013).  Virtual library of simulation experiments: test functions and datasets.  Simon Fraser University, Burnaby, BC, Canada, accessed May, 13:2015.
  * Svendsen et al., (2020) Svendsen, D. H., Martino, L., and Camps-Valls, G. (2020).  Active emulation of computer codes with gaussian processes–application to remote sensing.  Pattern Recognition, 100:107103.
  * Thévenot et al., (2001) Thévenot, D. R., Toth, K., Durst, R. A., and Wilson, G. S. (2001).  Electrochemical biosensors: recommended definitions and classification.  Analytical Letters, 34(5):635–659.
  * Wang and Dowling, (2022) Wang, K. and Dowling, A. W. (2022).  Bayesian optimization for chemical products and functional materials.  Current Opinion in Chemical Engineering, 36:100728.
  * Wang et al., (2020) Wang, W., Tuo, R., and Jeff Wu, C. (2020).  On prediction properties of kriging: Uniform error bounds and robustness.  Journal of the American Statistical Association, 115(530):920–930.
  * Wu and Frazier, (2016) Wu, J. and Frazier, P. (2016).  The parallel knowledge gradient method for batch bayesian optimization.  Advances in neural information processing systems, 29.
  * Yue and Kontar, (2020) Yue, X. and Kontar, R. A. (2020).  Why non-myopic bayesian optimization is promising and how far should we look-ahead? a study via rollout.  In International Conference on Artificial Intelligence and Statistics, pages 2808–2818. PMLR.
  * (66) Zhang, J., Liu, Y., Singh, M., Tong, Y., Kucukdeger, E., Yoon, H. Y., Haring, A. P., Roman, M., Kong, Z. J., Johnson, B. N., et al. (2023a).  Rapid, autonomous high-throughput characterization of hydrogel rheological properties via automated sensing and physics-guided machine learning.  Applied Materials Today, 30:101720.
  * (67) Zhang, J., Srivatsa, P., Ahmadzai, F. H., Liu, Y., Song, X., Karpatne, A., Kong, Z., and Johnson, B. N. (2023b).  Reduction of biosensor false responses and time delay using dynamic response and theory-guided machine learning.  ACS sensors, 8(11):4079–4090.
  * Zhang et al., (2024) Zhang, J., Srivatsa, P., Ahmadzai, F. H., Liu, Y., Song, X., Karpatne, A., Kong, Z. J., and Johnson, B. N. (2024).  Improving biosensor accuracy and speed using dynamic signal change and theory-guided deep learning.  Biosensors and Bioelectronics, 246:115829.
  * Zhang et al., (2020) Zhang, Y., Apley, D. W., and Chen, W. (2020).  Bayesian optimization for materials design with mixed quantitative and qualitative variables.  Scientific reports, 10(1):1–13.