# Bayesian Optimization for QAOA

Simone Tibaldi<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia dell’Università di Bologna, I-40127 Bologna, Italy  INFN, Sezione di Bologna, I-40127 Bologna, Italy  Davide Vodola  Dipartimento di Fisica e Astronomia dell’Università di Bologna, I-40127 Bologna, Italy  INFN, Sezione di Bologna, I-40127 Bologna, Italy  Edoardo Tignone  Leithà S.r.l. | Unipol Group, Bologna, Italy  Elisa Ercolessi  Dipartimento di Fisica e Astronomia dell’Università di Bologna, I-40127 Bologna, Italy  INFN, Sezione di Bologna, I-40127 Bologna, Italy

###### Abstract

The Quantum Approximate Optimization Algorithm (QAOA) adopts a hybrid quantum-
classical approach to find approximate solutions to variational optimization
problems. In fact, it relies on a classical subroutine to optimize the
parameters of a quantum circuit. In this work we present a Bayesian
optimization procedure to fulfil this optimization task, and we investigate
its performance in comparison with other global optimizers. We show that our
approach allows for a significant reduction in the number of calls to the
quantum circuit, which is typically the most expensive part of the QAOA. We
demonstrate that our method works well also in the regime of slow circuit
repetition rates, and that few measurements of the quantum ansatz would
already suffice to achieve a good estimate of the energy. In addition, we
study the performance of our method in the presence of noise at gate level,
and we find that for low circuit depths it is robust against noise. Our
results suggest that the method proposed here is a promising framework to
leverage the hybrid nature of QAOA on the noisy intermediate-scale quantum
devices.

##  1 Introduction

Hybrid quantum-classical variational algorithms [1, 2, 3] play a central role
in the current research on noisy intermediate-scale quantum (NISQ) devices
[4]. In a hybrid variatonal setting, a classical computer is entrusted with
the non-trivial task of optimizing the parameters of a quantum state. These
algorithms implement a heuristic protocol to approximately solve variational
problems including combinatorial optimization tasks, which are ubiquitous and
have a great practical importance [5], and are indeed one of the main drivers
of the industry interest towards quantum computing applications.
Unfortunately, the problems belonging to this class are hard to solve with
classical methods [6]. In this paper we focus on the Max-Cut and the Max
Indepedent Set problems defined on specific graph instances.

Among the hybrid variational algorithms, the Quantum Approximate Optimization
Algorithm (QAOA) [7] is extensively studied [8] as a promising algorithm to
investigate quantum speedups on NISQ devices and has been implemented on
several experimental platforms, such as Rydberg atom arrays [9],
superconducting processors [10], trapped-ions simulators [11], as well as
simulated on classical devices [12].

Similarly to other hybrid variational algorithms, QAOA consists in a sequence
of parametrized quantum gates applied to a wavefunction, on which an
expectation value of some operator, tipycally the Hamiltonian, is
reconstructed from measurements. The task of the classical subroutine is to
optimize the gate parameters in order to minimize such expectation value.
Every variational quantum algorithm therefore requires the estimation of the
wavefunction, which is a sub-optimal problem that requires resources that grow
exponentially with the system size [13]. Furthermore, the interplay between
the classical and quantum parts of the algorithm entails to run the quantum
circuit a large number of times, thus being expensive in term of resources.

To overcome these issues, an efficient classical optimization routine is
crucial. Different techniques have been proposed for optimizing variational
quantum circuits, e.g. Nelder-Mead [14], Monte-Carlo [15], gradient descent
[16], iterative schemes [8, 17, 18] and also Bayesian methods [19, 20, 21, 22,
23]. Here, we focus on a Bayesian optimization framework, which is suitable
for gradient-free global optimization of black-box functions [24, 25]. We
explore its behaviour in comparison with other global optimizers and we show
that the convergence rate to a minimum is faster. We also demonstrate that the
Bayesian approach is efficient in terms of number of circuit runs and is
robust against noise sources.

The paper is structured as follows. In Section 2 we introduce the QAOA
algorithm. In Section 3 we give a detailed presentation of the Bayesian
algorithm. In Section 4 we present the result of applying this method to QAOA,
compare it to other global optimization methods, evaluate its performance with
a limited number of circuit runs and against simulated quantum noise. Finally,
in Section 5 we draw our conclusions.

##  2 QAOA for combinatorial problems

The QAOA is a variational quantum algorithm that performs hybrid quantum-
classical optimization [7]. Given a cost function $C(z)$ with
$z=(z_{1},\dots,z_{i},\dots,z_{N})$ with $z_{i}\in\\{0,1\\}$, QAOA aims at
finding the bitstring $z^{\star}$ that minimizes the cost. In order to do so,
the cost function is translated into a quantum operator $H_{C}$ that is
diagonal in the computational basis $\ket{z}=\ket{z_{1}\dots z_{i}\dots
z_{N}}$ for $N$ qubits i.e. $H_{C}\ket{z}=C(z)\ket{z}$.

The QAOA circuit consists in preparing an initial state of $N$ qubits, usually
$|+\rangle=\sum_{z}|z\rangle/\sqrt{2^{N}}$, and then applying two unitary
operators alternatively: one generated by $H_{C}$, the other generated by
$H_{M}=\sum_{i}X_{i}$, where $X_{i}$ is the flip (NOT) operator on the $i$-th
qubit for a number of layers $p$ which is called the depth of the circuit. In
this way the QAOA circuit prepares the state

$|\bm{\theta}\rangle=\prod_{l=1}^{p}e^{-i\beta_{l}H_{M}}e^{-i\gamma_{l}H_{C}}|+\rangle,$
(1)

where $\bm{\theta}=(\bm{\gamma},\bm{\beta})$ are $2p$ parameters. By measuring
the state $\ket{\bm{\theta}}$ in the computational basis an estimate of the
energy $E(\bm{\theta})=\langle\bm{\theta}|H_{C}|\bm{\theta}\rangle$ is
obtained. This energy is fed to a classical routine which looks for the set of
angles $\bm{\theta}^{\star}=(\bm{\gamma}^{\star},\bm{\beta}^{\star})$ that
minimizes $E(\bm{\theta})$. Several strategies have been proposed for finding
the optimal parameters $\bm{\theta}^{\star}$. In this work we rely on Bayesian
optimization.

##  3 Bayesian optimization

Bayesian optimization is a global optimization strategy which allows one to
find within relatively few evaluations the minimum of a noisy, black-box
objective function $f(\bm{\theta})$ that is in general expensive to evaluate
[26]. The algorithm for the optimization method can be summarized as follows:
(i) it treats the objective function $f$ as a random function by choosing a
prior (also called surrogate model) for $f$. Several choices for the surrogate
model are possible [24], in this work we adopt the so-called Gaussian process
[27]. (ii) The prior is then updated through the likelihood function by
gathering observations of $f$ and therefore forming the posterior
distribution. (iii) The posterior distribution is finally used to construct an
auxiliary function, called acquisition function, that is in general cheap to
evaluate. The point where the acquisition function is maximized gives the next
point where $f$ will be evaluated [25].

Since Bayesian optimization requires no previous knowledge on $f$, it appears
to be a well-suited technique for optimizing the parameters of a variational
circuit running on NISQ devices.

In the following sections we describe the Gaussian process, the optimization
routine and the acquisition function in detail.

###  3.1 Gaussian process

Since the function $f(\bm{\theta})$ ($\bm{\theta}\in A\subset\mathbb{R}^{d}$)
to be optimized is unknown, we may think of it as belonging to a random
process, i.e an infinite collection of random variables defined for every
point $\bm{\theta}\in A$. A random process is called Gaussian if the joint
distribution of any finite collection of those random variables is a
multivariate normal distribution defined by a mean function $\mu(\bm{\theta})$
and covariance (or kernel) function $k(\bm{\theta},\bm{\theta}^{\prime})$
[27]. The mean function is related to the expected value of the function $f$,
while the kernel estimates the deviations of the mean function from the value
of $f$:

$\displaystyle\mu(\bm{\theta})$ $\displaystyle=$
$\displaystyle\mathbb{E}[f(\bm{\theta})],$ (2) $\displaystyle
k(\bm{\theta},\bm{\theta}^{\prime})$ $\displaystyle=$
$\displaystyle\mathbb{E}[(f(\bm{\theta})-\mu(\bm{\theta}))(f(\bm{\theta}^{\prime})-\mu(\bm{\theta}^{\prime}))],$
(3)

where $\mathbb{E}$ denotes the expectation value w.r.t. the infinite
collection of functions belonging to the random process. Conceptually, the
mean encloses the knowledge of the function $f$ to reconstruct while $k$
represents the uncertainty we have on such reconstruction.

Since we assume $f$ to be smooth, we choose for $k$ the Matérn kernel, a
stationary kernel [27] that depends on the distance between the points
$\bm{\theta}$ and $\bm{\theta}^{\prime}$, defined as

$k(\bm{\theta},\bm{\theta}^{\prime})=\sigma^{2}\bigg{(}1+\frac{\sqrt{3}\lVert\bm{\theta}-\bm{\theta}^{\prime}\rVert_{2}}{\ell}\bigg{)}e^{-\frac{\sqrt{3}\lVert\bm{\theta}-\bm{\theta}^{\prime}\rVert_{2}}{\ell}},$
(4)

where $\lVert\cdot\rVert_{2}$ is the 2-norm and $\sigma^{2}$ and $\ell$ are
two hyperparameters characterizing the Gaussian process. The hyperparameter
$\sigma^{2}$ defines the variance of the random variables whereas $\ell$ is a
characteristic length-scale which regulates the decay of the correlation
between points: in the limit of $\ell\rightarrow\infty$ all points are equally
correlated, for $\ell\rightarrow 0$ all points are uncorrelated.

###  3.2 Bayesian optimization

The main steps of the algorithm for Bayesian optimization can be summarized in
the pseudocode in Algorithm 1 (see also Appendix A for details).

Set the prior on $f$ as a Gaussian process;

Evaluate $f$ at $N_{W}$ different points $\bm{\theta}_{i}$;

Define the initial training set
$\mathcal{D}=\\{(\bm{\theta}_{i},f(\theta_{i}))\\}_{i=1}^{N_{W}}$;

Compute the hyperparameters $\sigma,\ell$ based on $\mathcal{D}$;

Set the guess for the minimum of $f$ to
$f_{m}=\min[\\{f(\theta_{i})\\}_{i=1}^{N_{W}}$];

while  _$n\leq N_{\mathrm{BAYES}}$_ do

Update the posterior distribution on $f$ using the training set $\mathcal{D}$;

Compute the acquisition function with the updated posterior;

Find $\tilde{\bm{\theta}}$ that maximes the acquisition function;

Evaluate $f(\tilde{\bm{\theta}})$;

if  _$f(\tilde{\bm{\theta}}) <f_{m}$_ then

Set the guess for the minimum of $f$ to $f_{m}=f(\tilde{\bm{\theta}})$;

end if

Append $(\tilde{\bm{\theta}},f(\tilde{\bm{\theta}}))$ to $\mathcal{D}$;

Compute the new hyperparameters $\sigma,\ell$;

Increment $n$;

end while

return $f_{m}$

Algorithm 1 Pseudo-code for Bayesian optimization

The optimization starts by a warmup phase where a number $N_{W}$ of
evaluations of the objective function $f$ is performed. These evaluations take
place at random chosen values of the parameters $\bm{\theta}_{i}$ and are
collected in the training set
$\mathcal{D}=\\{(\bm{\theta}_{i},y_{i}=f(\theta_{i}))\\}_{i=1}^{N_{W}}$ of the
optimization. Given the set $\mathcal{D}$, we define the design matrix
$\bm{\Theta}=(\bm{\theta}_{1},\dots,\bm{\theta}_{N_{W}})$ with the parameters,
the vector $\bm{y}\in\mathbb{R}^{N_{W}}$ with the observations via
$\bm{y}=(y_{1},\dots,y_{N_{W}})$ and the covariance matrix
$\bm{K}\in\mathbb{R}^{N_{W}\times N_{W}}$ formed by evaluating the covariance
function in Eq. (3) for each pair of points
$\bm{\theta}_{i},\bm{\theta}_{j}\in\bm{\Theta}$ via

$\bm{K}_{i,j}=k(\bm{\theta}_{i},\bm{\theta}_{j}),$ (5)

where $\bm{K}_{i,j}$ denotes the $(i,j)$ element of the matrix $\bm{K}$. The
hyperparameters entering the kernel function (Eq.(2)) are optimized as
explained in Sec. 3.4.

The training set will be used at each step of the optimization for
incorporating the acquired knowledge in the Gaussian process, by updating the
prior distribution and thus for predicting the values of $f$ at any point
$\bm{\theta}\in A$. This prediction is done by conditioning the prior Gaussian
process on the observations in $\mathcal{D}$ [27]. This operation returns a
posterior distribution still described by a Gaussian process multinomial
distribution characterized by a posterior mean $\mu^{\prime}$ and covariance
$k^{\prime}$ given by

$\displaystyle\mu^{\prime}$
$\displaystyle=\bm{\kappa}^{T}\cdot\bm{K}^{-1}\cdot\bm{y}$ (6) $\displaystyle
k^{\prime}$
$\displaystyle=k(\bm{\theta},\bm{\theta})-\bm{\kappa}^{T}\cdot\bm{K}^{-1}\cdot\bm{\kappa},$
(7)

where $\bm{\theta}$ is a generic point in $A$, $\bm{\kappa}$ is a column
vector formed by evaluating the covariance function $k$ between the generic
point $\bm{\theta}$ and all the points in $\bm{\Theta}$, i.e. its $j$-th
element is $\kappa_{j}=k(\bm{\theta},\bm{\theta}_{j})$. Eq. (6) shows that the
new mean is a linear combination of the observations $\bm{y}$.

###  3.3 Acquisition Function

The next step in the Bayesian optimization is computing the acquisition
function whose maximum gives the next point at which to evaluate the objective
function. A common choice of acquisition function is the expected improvement
that suggests to evaluate $f$ in a point that, on average, improves on $f_{m}$
the most [25]. This choice corresponds to defining the acquisition function
$\textrm{EI}(\bm{\theta})=\mathbb{E}[u(\bm{\theta})]$ as the average value of
the utility function $u(\bm{\theta})=\max[0,f_{m}-f(\bm{\theta})]$ such that
the lower $f(\bm{\theta})$ is with respect to the current minimum, the larger
the utility $u(\bm{\theta})$ will be.

By considering that $f(\bm{\theta})$ is a Gaussian process, we can obtain an
analytical expression for $\mathrm{EI}(\bm{\theta})$ as

$\mathrm{EI}(\bm{\theta})=\Phi(z)(f_{m}-\mu^{\prime})+\phi(z)k^{\prime},$ (8)

where $\mu^{\prime}$ and $k^{\prime}$ are obtained for the point $\bm{\theta}$
by using Eqs. (6) and (7); $\Phi(\cdot)$ and $\phi(\cdot)$ are respectively
the cumulative distribution function and the probability density function of
the standard normal distribution and the quantity $z$ is defined as
$z=(f_{m}-\mu^{\prime})/k^{\prime}$. The two terms in Eq. (8) resume the
trade-off between exploitation and exploration: the first term, being
proportional to the difference between the current minimum and the mean value
of the posterior, brings the optimization towards points with lower
$\mu^{\prime}$ whereas the second one promotes points with larger
$k^{\prime}$, i.e. with higher uncertainty.

In this work, we compute the point $\tilde{\bm{\theta}}$ that maximizes the
acquisition function by using the differential evolution algorithm [28], a
population-based metaheuristic search algorithm (see Appendix C for details).

Figure 1: Energy distributions of graphs. (a) MIS: Distribution of the
energies of the possible bitstrings for the graph of 6 nodes (shown in the
inset, nodes in red correspond to one solution). The red bar to the left
highlights the two solution bitstrings of the MIS problem on such graph. (b)
Max-Cut: Distribution of the energies of the possible bitstrings for the graph
of 10 nodes (shown in the inset, nodes in red correspond to one solution). The
red bar highlights the two solution bitstrings of the Max-Cut problem on such
graph.

###  3.4 Hyperparameters

We are now only left with the task of picking the best hyperparameters
$\sigma,\ell$ for the Matérn kernel. This is typically done by considering the
marginal likelihood [27]

$p(\bm{y}|\bm{\Theta})=\int
p(\bm{y}|\bm{f},\bm{\Theta})p(\bm{f}|\bm{\Theta})d\bm{f},$ (9)

where the prior $p(\bm{f}|\bm{\Theta})$ and the likelihood
$p(\bm{y}|\bm{f},\bm{\Theta})$ are Gaussian and the marginalization is done
over the function values $\bm{f}$. Given the Gaussian nature of the likelihood
and the prior, a closed form of the log marginal likelihood can be obtained
(for more details see [27]):

$\displaystyle\begin{split}\log
p(\bm{y}|\bm{\Theta})=&-\frac{1}{2}\bm{y}^{T}\cdot\bm{K}^{-1}\cdot\bm{y}\\\
&-\frac{1}{2}\log\det\bm{K}-\frac{N}{2}\log 2\pi,\end{split}$ (10)

where $N$ is the number of observations in the design matrix $\bm{\Theta}$. In
Eq. (10), the first term specifies how well the process fits the data, the
second term instead acts as a regularization factor on the elements of the
kernel matrix. When fitting the Gaussian process to a new set of points, the
best hyperparameters $(\tilde{\sigma}^{2},\tilde{\ell})$ can be found by
maximizing the log marginal likelihood in Eq. (10). For the optimization of
$\log p(\bm{y}|\bm{\Theta})$, we use the quasi-Newton method L-BFGS [29] with
multiple restarting points which proved to be efficient on the flat landscape
of the likelihood (see Appendix A for details).

##  4 Results

In this section we apply the Bayesian optimization to the QAOA parameters. We
consider two well-known combinatorial problems defined on graphs, the Max-Cut
and the MIS.

Max-Cut – Given a graph $G=(V,E)$ where $V$ is the set of nodes and $E$ the
set of edges, the Max-Cut problem consists in finding a partition of the
graph’s vertices $V$, $P=\\{V_{0},V_{1}\\}$, such that the number of edges
between $V_{0}$ and its complement $V_{1}$ is as large as possible. It is
known to be a NP-hard problem [30]. We can define the assignment of the nodes
to the sets $V_{0}$ and $V_{1}$ by labelling with label “0” the nodes $v\in
V_{0}$ and with label “1” the nodes $v\in V_{1}$. In these terms, the Max-Cut
consists in finding the largest number of edges connecting the bits labelled
with “0” to the bits labelled with “1”. On a quantum computer, the labels 0
and 1 are replaced by the computational basis states $\ket{0}$ and $\ket{1}$
and the cost Hamiltonian can be written as:

$H_{C}^{\text{MC}}=-\sum_{(i,j)\in E}(1-Z_{i}Z_{j})/2.$ (11)

The cost function has minimum eigenvalue on the states $\ket{z_{i}z_{j}}$ for
which $z_{i}\neq z_{j}$ which represent separate partitions.

Figure 2: QAOA at $\bm{p=1}$. (a) Landscape of the energy $E(\bm{\theta}_{1})$
obtained on the 6 nodes graph solving the MIS problem. The red cross indicates
the angles corresponding to the final state $\ket{\bm{\theta}_{F}}$ with
largest fidelity. (b) Landscape of the fidelity $F(\bm{\theta}_{1})$. (c)
Squared amplitudes of the two states $\ket{\bm{\theta}_{E}}$,
$\ket{\bm{\theta}_{F}}$. The solution bitstrings are highlighted in red.
Values of the qubits are given in the order shown in the inset of Fig. 1(a).

MIS – The Max Independent Set problem consists in finding the largest number
of graph nodes which are not adjacent and the corresponding cost quantum
Hamiltonian is:

$H_{C}^{\text{MIS}}=\sum_{i}Z_{i}+\omega\sum_{(i,j)\in E}Z_{i}Z_{j},$ (12)

with $\omega$ a parameter that balances the effect of the first term (which
maximizes the number of qubits in $\ket{1}$) and the second one (which
prevents neighbour bits to be equal). We set $\omega=2$ for the rest of the
work.

During the optimization process we monitor the approximation ratio
$R=E(\bm{\theta})/E_{GS}$ where $E_{GS}$ is the energy of the solution
bitstring. We also look at the fidelity defined as
$F=|\braket{\bm{\theta}}{z^{\star}}|^{2}$ where $\ket{z^{\star}}$ is the state
which encodes the solution. The following results are obtained on two
3-regular graphs of 6 and 10 nodes which are plotted in the insets of Figs.
1(a)-(b).

QAOA at low vs. large depth – We start by looking at the QAOA at depth $p=1$
for the 6 nodes graph. It corresponds to a shallow circuit that depends only
on two parameters $\bm{\theta_{1}}=(\gamma_{1},\beta_{1})$. We consider the
MIS problem, and we plot the landscapes of both the energy
$E(\bm{\theta}_{1})$ and the fidelity $F(\bm{\theta}_{1})$ (Fig. 2(a) and (b),
respectively) for values of $\gamma_{1},\beta_{1}\in[0,\pi]$ due to the
simmetry of the problem. We see that the landscape of the energy, which is the
function to minimize, is rather flat with two global maxima and minima,
corresponding to the best solutions possible at $p=1$.

Interestingly, we find that the QAOA state
$\ket{\bm{\theta_{E}}}=\sum_{z}\alpha_{z,E}\ket{z}$, corresponding to the
parameters which minimize the energy, is not the state
$\ket{\bm{\theta_{F}}}=\sum_{z}\alpha_{z,F}\ket{z}$ with largest fidelity. To
see how they differ we plot the squared amplitudes
$\lvert\alpha_{z,E}\rvert^{2}$ and $\lvert\alpha_{z,F}\rvert^{2}$ of both
states in Fig. 2(c) as histograms. The fidelity of $\ket{\bm{\theta_{F}}}$
w.r.t. the solution $\ket{z^{\star}}$ is, as expected, much larger than that
one of $\ket{\bm{\theta_{E}}}$, yet the latter has a lower energy because it
has many non zero amplitudes along excited states with low energy. This
unravels the problem of optimizing the QAOA parameters by only looking at the
energy $E(\bm{\theta})$. There is a large concentration of excited states with
energy comparable to the energy of the ground state, as shown in the
histograms of panels (a) and (b) of Fig. 1. It is difficult to increase the
amplitude corresponding to the solution when many other states can contribute
with low values of the energy.

The divergence between lowest energy and highest fidelity is guaranteed to
disappear theoretically for $p\rightarrow\infty$. For this reason, we apply
Bayesian optimization to the problem and we show in Fig. 3 that the
approximation ratio and fidelity both tend to 1 for $p\sim 12$. Yet, we
already see a good performance at $p=4$ where $R\sim 0.7$ and we have $F\sim
0.5$ meaning about a $50\%$ chance of measuring the solution on the state
obtained with QAOA.

Figure 3: Results increasing depth. Average approximation ratio (plotted as
$1-R$) and fidelity $F$ for increasing values of circuit depth from 1 to 12
over 50 runs. Shaded areas correspond to one standard deviation. Results were
obtained on the 6 nodes graph of Fig. 1(a).  Figure 4: Comparison among
optimizers. (a) The plot shows the average number of calls $N_{C}$ to the
quantum circuit of each optimizer in order to obtain the same approximation
ratio as the Bayesian optimization. (b) - (c) Average approximation ratio
(plotted as $1-R$) and fidelity during the optimization with the different
methods at $p=7$ over 30 runs. Shaded areas correspond to one standard
deviation. Results are obtained on the 10 nodes graph of Fig. 1(b).

Comparing resources –  Increasing the depth of a variational circuit increases
the number of parameters that must be optimized. In turn, this is expected to
increase the number of calls to the quantum circuit needed to reach a good
approximate solution, which is a problem in the current NISQ era, since
running a quantum circuit can be costly due to both state preparation routines
and recalibrations of the device.

Bayesian optimization mitigates such problem, and allows to achieve a good
approximate solution within relatively fewer calls to the circuit compared to
other global optimization methods. To show this, we ran differential
evolution, basin-hopping and dual annealing (see Appendix C and D for details)
on the 10 nodes graph of Fig. 1(b) for the Max-Cut problem at different
depths. Fig. 4(a) shows the average number of calls to the circuit that the
other optimizers need in order to reach the same approximation ratio of
Bayesian optimization. To gain more insight, we plot in panels (b) and (c) of
Fig. 4 the average approximation ratio and fidelity during the run of the
algorithm for each method at $p=7$. We see that Bayesian optimization stops at
a lower $R$ than basin-hopping and dual annealing, but it reaches a value of
$R=95\%$ with $\sim 500$ runs of the circuit compared to the other two methods
which take, in order, 1400 and 10800.

Slow Measurements – The energy $E(\bm{\theta})$ is obtained by measuring the
QAOA state after running the circuit: we refer to these two operations
combined together as a “shot”. By measuring on the $Z$ basis at each shot we
get a bitstring, and we calculate its classical energy associated to the
combinatorial problem. The precision in the reconstruction of $E(\bm{\theta})$
depends on the number of shots $N_{S}$. Since we consider this as a
multinomial sampling problem we expect the variance of the reconstructed
energy to depend on $N_{S}^{-1}$. In many scenarios of NISQ devices it is
necessary to balance $N_{S}$ with the desired standard deviation. For this
reason, we compare the average approximation ratio obtained with the exact
energy (simulated) with the energy reconstructed with a limited number of
shots.

We show in Fig. 5 such comparison with a number of shots $N_{S}$ equal to
1024, 128, 64, 16, 4. Looking at the approximation ratio $R$ (Fig. 5(a)) we
see that taking $N_{S}=128$ shots reduces $R$ by $5\%$ w.r.t. $N_{S}=1024$ and
going to $N_{S}=64$ reduces it by a further $5\%$. This behaviour then stops
and even reverses its trend. In fact, we even see an average increase going
from $N_{S}=16$ to $4$. This is understandable since the reconstruction of the
energy with as little as 4 shots is not indicative of the real energy of the
state. Specifically, from a final QAOA state we might sample the solution
bitstring 2, 3 or even 4 times out of 4 and the expectation of the energy on
these three samplings would be very different. This behaviour is indeed
confirmed also by the fidelity in Fig. 5(a) which follows the same trend as
the approximation ratio.

To have a better understanding of how the algorithm adapts to the sampling
noise, we look at the kernel noise parameter $\sigma^{2}_{N}$ which is learned
by the Gaussian process during the fitting at each step of the optimization
(see Appendix B for details on the noise hyperparameter). The plot in Fig.
5(b) shows that, after an initial phase, the kernel noise sets at a specific
value at around 400 steps. In addition to that, the lower the number of shots
the larger the noise parameter learned. In fact, by fitting the average kernel
noise found at the end of the training (Fig. 5(c)) we obtain that
$\sigma^{2}_{N}$ follows a power law with $N_{S}^{-1.1}$. This trend is
comparable to the expected trend for the variance $N_{S}^{-1}$ of the
reconstruction of the energy. This shows that the Gaussian process adapts to
sampling noise.

Figure 5: Slow measurements. (a) Average approximation ratio (plotted as
$1-R$) and fidelity ($F$) as a function of the number of shots $N_{S}$. The
$1/N_{S}=0$ points indicates the exact evaluation of the energy
$E(\bm{\theta})$. Shaded areas correspond to one standard deviation. (b)
Kernel noise $\sigma^{2}_{N}$ learned by fitting the data at each step of the
optimization for different number of shots $N_{S}$. (c) Average kernel noise
learned by the Gaussian process as a function of $N_{S}$ (blue circles). The
plot also shows a linear fit ($\sim 1/N^{1.1}$, orange line) of the logarithm
of the data.

Quantum noise – Another relevant issue in the state-of-the-art NISQ devices
are the sources of quantum noise which can interfere with the quantum circuit.
Every device has different sources of noise depending on the underlying
technology. In order to simulate it without specifying the device technology
we add random noise on every QAOA parameter. In this way Eq. (1) for the Max-
Cut problem is modified as:

$\displaystyle\begin{split}|\bm{\theta}\rangle=\prod_{l=1}^{p}e^{-i\sum_{i}\beta_{l}^{i}X_{i}}\
e^{i\sum_{\langle
i,j\rangle}\gamma_{l}^{(i,j)}Z_{i}Z_{j}}|+\rangle,\end{split}$ (13)

where $\beta_{l}^{i}$ and $\gamma_{l}^{(i,j)}$ act differently on every
qubit/edge of the graph at every layer because they are affected by Gaussian
random noise with mean zero and standard deviation $\sigma_{QN}$. In Fig. 6 we
plot $R$ as a function of depth at different values of $\sigma_{QN}$ on the 10
nodes graph for the Max-Cut problem.

By increasing $\sigma_{QN}$ and $p$, we expect to obtain a worse approximation
ratio $R$ because the error accumulates along the circuit as the number of
parameters grows. Indeed, as we can see in the figure, from $p\geq 5$ the
obtained $R$ decreases w.r.t. the noiseless case, decreasing even by $20\%$
for $p=9$ with $\sigma_{QN}=0.1$. For $p\leq 3$, however, $R$ increases with
$p$, which indicates that, for shallow circuits, Bayesian optimization is
robust against noise.

Figure 6: Approximation ratio $R$ for different values of the quantum noise
$\sigma_{QN}$. The noise is simulated by adding random Gaussian noise with
mean zero and standard deviation $\sigma_{QN}$ to the variational parameters
$(\bm{\gamma},\bm{\beta})$. The plot shows the effects of the noise
$\sigma_{QN}$ on the final obtained approximation ratio $R$ as a function of
the QAOA depth $p$ for different $\sigma_{QN}$. Shaded areas correspond to one
standard deviation. The results are obtained on the graph with 10 nodes in
Fig. 1(b).

##  5 Conclusions

In this work we have presented the Bayesian optimization algorithm as a
subroutine to optimize the variational parameters of the QAOA. We have applied
it to find the solutions of two combinatorial optimization problems, the Max-
Cut and the MIS on two graph instances.

After introducing the QAOA and the details of the Bayesian optimization
algorithm, we have focused on its capability to adapt to the data and to
predict new possible optimal points by exploiting both the accumulated
knowledge from the previous observations and the uncertainty with respect to
the optimization landscape.

We have analyzed some details of the QAOA at low circuit depth with the
purpose of presenting some of the issues related to the optimization of a
variational quantum algorithm. These include the flatness of the energy
landscape and the limited information that we can retain from the energy of
the QAOA state compared to its overlap with the ground state. After that, we
have compared the Bayesian optimization with other global optimization
methods, and we have showed that they require more calls, in the order of tens
or hundreds, to the quantum circuit with respect to the Bayesian optimizer.
This is a first sign that this method responds more efficiently to the
requirements from the quantum part of the QAOA.

Secondly, we have analyzed the effects of a finite number of measurements for
the reconstruction of the energy landscape. We have showed that the results
are slightly altered by a $5\%$ decrease in the approximation ratio by using
1024 measurements compared to the optimization with the exact energy. A lower
number of measurements will results in a decreasing approximation ratio. We
have also showed that the Gaussian process learns to add a noise
hyperparameter which is proportional to the variance expected from the
reconstruction of the energy. This can be seen as a further example of
adaptation of the Bayesian algorithm to the data.

Finally, we have simulated a noisy algorithm and we have showed that for
shallow circuits, with depth $p\in[1,3]$, the approximation ratio is improved
even for large value of the noise. For deeper circuit, $p\geq 5$, the
intensity of noise sensibly affects the final approximation ratio, as
expected.

These findings show that Bayesian optimization is a robust method which can
account for both quantum and sampling noise. For this reason it represents a
valid tool for solving optimization problems via hybrid algorithms to be run
on a NISQ device.

We thank C. Sanavio and F. Dell’Anna for useful discussions. This research is
funded by the International Foundation Big Data and Artificial Intelligence
for Human Development (IFAB) through the project “Quantum Computing for
Applications”. E. Ercolessi, S. Tibaldi and D. Vodola are partially supported
by Istituto Nazionale di Fisica Nucleare (INFN) through the project “QUANTUM”
and the QuantERA 2020 Project “QuantHEP”.

## References

  * McClean et al. [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18, 023023 (2016).
  * Peruzzo et al. [2014] A. Peruzzo, J. R. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5, 4213 (2014).
  * Moll et al. [2018] N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow, A. Cross, D. J. Egger, S. Filipp, A. Fuhrer, J. M. Gambetta, M. Ganzhorn, A. Kandala, A. Mezzacapo, P. Müller, W. Riess, G. Salis, J. Smolin, I. Tavernelli, and K. Temme, Quantum optimization using variational algorithms on near-term quantum devices, Quantum Science and Technology 3, 030503 (2018).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Du and Pardalos [1999] D. Du and P. M. Pardalos, eds., Handbook of Combinatorial Optimization (Springer, 1999).
  * Garey and Johnson [1982] M. R. Garey and D. S. Johnson, Computers and Intractability: A Guide to the Theory of NP-Completeness (W. H. Freeman and Co, 1982).
  * Farhi et al. [2014] E. Farhi, J. Goldstone, and S. Gutmann, A Quantum Approximate Optimization Algorithm, (2014), arXiv:1411.4028.
  * Zhou et al. [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Quantum Approximate Optimization Algorithm: Performance, Mechanism, and Implementation on Near-Term Devices, Phys. Rev. X 10, 021067 (2020).
  * Ebadi et al. [2022] S. Ebadi et al., Quantum optimization of maximum independent set using Rydberg atom arrays, Science 376, 1209 (2022).
  * Harrigan et al. [2021] M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, S. Boixo, and et al., Quantum approximate optimization of non-planar graph problems on a planar superconducting processor, Nature Physics 17, 332 (2021).
  * Pagano et al. [2020] G. Pagano, A. Bapat, P. Becker, K. S. Collins, A. De, P. W. Hess, H. B. Kaplan, A. Kyprianidis, W. L. Tan, C. Baldwin, L. T. Brady, A. Deshpande, F. Liu, S. Jordan, A. V. Gorshkov, and C. Monroe, Quantum approximate optimization of the long-range Ising model with a trapped-ion quantum simulator, Proceedings of the National Academy of Sciences 117, 25396 (2020).
  * Medvidović and Carleo [2021] M. Medvidović and G. Carleo, Classical variational simulation of the Quantum Approximate Optimization Algorithm, npj Quantum Information 7, 101 (2021).
  * Cramer et al. [2010] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, Efficient quantum state tomography, Nature Communications 1, 149 (2010).
  * Guerreschi and Smelyanskiy [2017] G. G. Guerreschi and M. Smelyanskiy, Practical optimization for hybrid quantum-classical algorithms, (2017), arXiv:1701.01450 [quant-ph].
  * Wecker et al. [2016] D. Wecker, M. B. Hastings, and M. Troyer, Training a quantum optimizer, Phys. Rev. A 94, 022309 (2016).
  * Wang et al. [2018] Z. Wang, S. Hadfield, Z. Jiang, and E. G. Rieffel, Quantum approximate optimization algorithm for MaxCut: A fermionic view, Phys. Rev. A 97, 022304 (2018).
  * Mele et al. [2022] A. A. Mele, G. B. Mbeng, G. E. Santoro, M. Collura, and P. Torta, Avoiding barren plateaus via transferability of smooth solutions in Hamiltonian Variational Ansatz, arXiv: Quantum Physics (2022), arXiv:2206.01982 [quant-ph].
  * Mbeng et al. [2019] G. B. Mbeng, R. Fazio, and G. E. Santoro, Quantum Annealing: a journey through Digitalization, Control, and hybrid Quantum Variational schemes, arXiv: Quantum Physics (2019), arXiv:1906.08948 [quant-ph].
  * Otterbach et al. [2017] J. S. Otterbach et al., Unsupervised Machine Learning on a Hybrid Quantum Computer, (2017), arXiv:1712.05771 [quant-ph].
  * Zhu et al. [2019] D. Zhu, N. M. Linke, M. Benedetti, K. A. Landsman, N. H. Nguyen, C. H. Alderete, A. Perdomo-Ortiz, N. Korda, A. Garfoot, C. Brecque, L. Egan, O. Perdomo, and C. Monroe, Training of quantum circuits on a hybrid quantum computer, Science Advances 5, eaaw9918 (2019).
  * Self et al. [2021] C. N. Self, K. E. Khosla, A. W. R. Smith, F. Sauvage, P. D. Haynes, J. Knolle, F. Mintert, and M. S. Kim, Variational quantum algorithm with information sharing, npj Quantum Information 7, 116 (2021).
  * Wang et al. [2021] G. Wang, D. E. Koh, P. D. Johnson, and Y. Cao, Minimizing Estimation Runtime on Noisy Quantum Computers, PRX Quantum 2, 010346 (2021).
  * Tamiya and Yamasaki [2022] S. Tamiya and H. Yamasaki, Stochastic gradient line Bayesian optimization for efficient noise-robust optimization of parameterized quantum circuits, npj Quantum Information 8, 90 (2022).
  * Shahriari et al. [2016] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. de Freitas, Taking the Human Out of the Loop: A Review of Bayesian Optimization, Proceedings of the IEEE 104, 148 (2016).
  * Frazier [2018] P. I. Frazier, A Tutorial on Bayesian Optimization, (2018), arXiv:1807.02811 [stat.ML].
  * Snoek et al. [2012] J. Snoek, H. Larochelle, and R. P. Adams, Practical Bayesian Optimization of Machine Learning Algorithms, in Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 2, NIPS’12 (Curran Associates Inc., Red Hook, NY, USA, 2012) p. 2951–2959.
  * Rasmussen and Williams [2005] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) (The MIT Press, 2005).
  * Price et al. [2006] K. Price, R. M. Storn, and J. A. Lampinen, Differential evolution: a practical approach to global optimization (Springer Science & Business Media, 2006).
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, On the limited memory BFGS method for large scale optimization, Mathematical Programming 45, 503 (1989).
  * Edwards [1973] C. S. Edwards, Some Extremal Properties of Bipartite Subgraphs, Canadian Journal of Mathematics 25, 475 (1973).
  * Wales and Doye [1997] D. J. Wales and J. P. K. Doye, Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms, The Journal of Physical Chemistry A 101, 5111 (1997).
  * Xiang et al. [1997] Y. Xiang, D. Sun, W. Fan, and X. Gong, Generalized simulated annealing algorithm and its application to the Thomson model, Physics Letters A 233, 216 (1997).

##  Appendix A Details on the Bayesian optimization

In this appendix we give a detailed review of the algorithm we used in this
work for performing the Bayesian optimization for the QAOA.

Figure 7: Posterior and acquisition function during optimization. Here we plot
the posterior and the acquisition at three different steps of the
optimization: (from left to right) at warmup, after 10 steps and after 20
steps. In rows (a) and (b) we plot respectively the mean
$\mu^{\prime}(\bm{\theta})$ and the variance $k^{\prime}(\bm{\theta})$ of the
posterior while in row (c) the acquisition function $\mathrm{EI}(\bm{\theta})$
(see Eqs. (6) - (8) in the text). These data, obtained running the Bayesian
optimization on QAOA with $p=1$ on the graph of Fig. 1(a), are shown as a
function of the two variational parameters
$\bm{\theta}=(\gamma_{1},\beta_{1})$. The red dots indicate the warmup points
($N_{W}=5$, in this case) while the crosses are the points that have been
selected by subsequents Bayesian optimization steps as new candidate solutions
(with a color scale from first (orange) to last (pink)). At every Bayesian
optimization step, $\mu^{\prime}(\bm{\theta})$ encodes the knowledge of the
landscape of the function to optimize, while $k^{\prime}(\bm{\theta})$
contains the uncertainty. The acquisition function $\mathrm{EI}(\bm{\theta})$
combines the information from $\mu^{\prime}$ and $k^{\prime}$ and the position
of its maximum proposes the next possible optimal point. From these
considerations, we see that after 20 steps the mean
$\mu^{\prime}(\bm{\theta})$ (row (a)) recreates the landscape (our knowledge
of $E(\bm{\theta}_{1})$) with more precision than at warmup (compare with the
real energy landscape in Fig. 2(a)). In addition to that,
$\mathrm{EI}(\bm{\theta})$ becomes more and more flat in all the landscape
except for a small area in the top right corner on the right panel of (c).
This means that the Gaussian process has acquired enough knowledge from the
data to converge to the energy minimum.

###  A.1 Bayesian optimization algorithm

This algorithm is made out of three phases: warmup, kernel optimization,
acquisition function maximization and is summarized in Algorithm 1.

Warmup – In the warmup phase we start with a set of $N_{W}=10$ points
$X=\\{\bm{\theta}_{j}\\}_{j=1}^{N_{W}}$ with
$\bm{\theta}_{j}\in\mathbb{R}^{d}$ where $d=2p$ and $p$ is the depth of the
QAOA circuit. Each point $\bm{\theta}_{j}$ is a set of angles
$(\bm{\gamma},\bm{\beta})$. These points are sampled from the latin hypercube
of bounds $[0,\pi]^{d}$. For each set of angles we estimate the energy of QAOA
$y_{j}$ and we create the design matrix
$\bm{\Theta}=(\bm{\theta}_{1},\dots,\bm{\theta}_{N_{W}})$ and the observation
vector $\bm{y}$ with the energies of the points $\bm{\theta}_{j}$. We also
store the point with lowest energy as $(\bm{\theta}_{m},f_{m})$.

Kernel optimization – In this part we look for the hyperparameters
$(\tilde{\sigma},\tilde{\ell})$ which maximize the likelihood function
$p(\bm{y}|\bm{\Theta})$ (Eq. (10)). This optimization is performed by
repeating the L-BFGS [29] minimization on $-p(\bm{y}|\bm{\Theta})$ for 10
times and selecting the best parameters found. The parameters found at every
step of the optimization are plotted in panels (c) and (d) of Fig. 8.

Acquisition function maximization – Once the hyperparameters are set the
algorithm exploits its knowledge and uncertainty of the data to propose a new
point $\bm{\theta}^{\prime}$ with the new parameters where we evaluate the
QAOA circuit. This is done by maximizing the expected improvement in Eq. (8).

The new point $\bm{\theta}^{\prime}$ maximizing the expected improvement is
then added to the dataset $\mathcal{D}$ and the algorithm is repeated. The
procedure stops after $N_{\mathrm{BAYES}}$ iterations. Fig. 7 provides an
illustrative example (with $p=1$, $N_{W}=5$, and $N_{\mathrm{BAYES}}=20$) of
how the Bayesian method operates and moves through the landscapes.

##  Appendix B Bayesian optimization with noise

In this paper we consider two scenarios in which the energy $E(\bm{\theta})$
is affected by a source of noise: the finite number of samplings and the
quantum noise at the gate level. In the context of Bayesian optimization we
can account for noise modifying the kernel function by adding a term
$\sigma^{2}_{N}\mathbb{I}$ like so:

$k(\bm{\theta},\bm{\theta}^{\prime})=\sigma^{2}\bigg{(}1+\frac{\sqrt{3}\lVert\bm{\theta}-\bm{\theta}^{\prime}\rVert_{2}}{\ell}\bigg{)}e^{-\frac{\sqrt{3}\lVert\bm{\theta}-\bm{\theta}^{\prime}\rVert_{2}}{\ell}}+\sigma_{N}^{2}\mathbb{I}$
(14)

This is usually called a white kernel and it is a parameter added to the
diagonal to account for random fluctuations around the true value of
$f(\bm{\theta})$. In this way the new predicted mean and covariance for a set
of data $\bm{\Theta}$ can be easily calculated to be:

$\displaystyle\mu^{\prime}$
$\displaystyle=\bm{\kappa}^{T}\cdot(\bm{K}+\sigma^{2}_{N}\mathbb{I})^{-1}\cdot\bm{y}$
(15) $\displaystyle k^{\prime}$
$\displaystyle=k(\bm{\theta},\bm{\theta})-\bm{\kappa}^{T}\cdot(\bm{K}+\sigma^{2}_{N}\mathbb{I})^{-1}\cdot\bm{\kappa}$
(16)

The constant $\sigma_{N}^{2}$ belongs to the list of hyperparameters (along
with $\sigma^{2}$ and $\ell$) that are optimized with the log marginal
likelihood which now takes the form

$\displaystyle\log p(\bm{y}|\bm{\Theta})=-$
$\displaystyle\frac{1}{2}\bm{y}^{T}\cdot(\bm{K}+\sigma^{2}_{N}\mathbb{I})^{-1}\cdot\bm{y}$
$\displaystyle-$
$\displaystyle\frac{1}{2}\log\det(\bm{K}+\sigma^{2}_{N}\mathbb{I})$
$\displaystyle-$ $\displaystyle\frac{N}{2}\log 2\pi.$

We show how the hyperparameter $\sigma^{2}_{N}$ is learned during training in
the text in Fig. 5(b).

Figure 8: Parameters of Bayesian optimization at p = 2, 7. Plots of the
parameters changing during Bayesian optimization for two runs with
$N_{BAYES}=600$ steps. (a) Approximation ratio $R$. (b) Fidelity $F$. (c)
Kernel constant $\sigma^{2}$. (d) Kernel correlation length $\ell$. (e)
Standard deviation of the expected improvement and (f) average distance of the
points of differential evolution at the last generation $N_{T}$.

##  Appendix C Differential evolution

Finding the point $\tilde{\bm{\theta}}$ of the parameter space that maximizes
the expected improvement $\mathrm{EI}(\bm{\theta})$ (Eq. 8) is not an easy
task since $\mathrm{EI}(\bm{\theta})$ can show a fairly flat landscape [26],
in particular after many optimization steps (for example, see Fig. 7(c)).

To compute the maximum of $\mathrm{EI}(\bm{\theta})$ in this work we use the
differential evolution algorithm [28]. This is an evolutionary method in which
populations of points $\\{\bm{\theta}\\}$, called generations, are iteratively
obtained from the previous ones until convergence. The algorithm starts by
initializing a generation (a) and then the population is updated following
three main steps: (b) mutation, (c) cross-over and (d) selection.

(a) We choose as starting population $N_{P}=15\cdot 2p$ points
$\\{\bm{\theta}_{i,1}\\}$ where the index $i\in\\{1,\dots,N_{P}\\}$ uniquely
identifies the point within the belonging population, while the index $1$
indicates that the point belongs to the first generation $g=1$. These points
are randomly generated on the latin $2p$-cube of bounds $[0,\pi]^{2p}$ and to
each point there is an associated expected improvement
$\mathrm{EI}(\bm{\theta}_{i,1})$.

(b) For each $\bm{\theta}_{i,g}$ (called parent point) in the population, the
differential evolution picks three random points, different from
$\bm{\theta}_{i,g}$, labelled by $r_{0},r_{1},r_{2}$ within the corresponding
population, and creates a new point as:

$\bm{v}_{i,g}=\bm{\theta}_{r_{0},g}+M(\bm{\theta}_{r_{1},g}-\bm{\theta}_{r_{2},g}),$
(17)

where $M\in(0.5,1)$ is a hyperparameter that is selected randomly at every
generation. Through Eq. 17 the differential evolution mutates and recombines
the population to create another set of parent points $\bm{v}_{i,g}$.

(c) A new point $\bm{u}_{i,g}$ (offspring point) is created from
$\bm{\theta}_{i,g}=(\theta_{i,g}^{1},\dots,\theta_{i,g}^{2p})$ and
$\bm{v}_{i,g}=(v_{i,g}^{1},\dots,v_{i,g}^{2p})$ choosing randomly between
their coordinates $\theta_{i,g}^{j}$ and $v_{i,g}^{j}$ for every
$j=1,\dots,2p$.

(d) Finally, if $\mathrm{EI}(\bm{u}_{i,g})\geq\mathrm{EI}(\bm{\theta}_{i,g})$
the algorithm replaces $\bm{\theta}_{i,g}$ with $\bm{u}_{i,g}$ in the next
generation, otherwise $\bm{\theta}_{i,g}$ is kept.

Steps (b) through (d) are repeated for $N_{T}$ generations $g$. The algorithm
stops when two convergence criteria are fulfilled: the standard deviation
$\sigma_{EN}$ of the population’s expected improvement (see Fig. 8(e)) and the
average distance $D$ among the population points (see Fig. 8(f)) are below a
certain threshold (we set $\sigma_{EN}=D=10^{-3}$). When this happens, the
point in the population with the maximum expected improvement is selected as
$\tilde{\bm{\theta}}$. We notice that the criterion on the distance guarantees
that in a flat landscape like the one of $\mathrm{EI}(\bm{\theta})$ the points
do not get stuck on a plateau and concentrate closer to a unique candidate
maximum. Although the algorithm requires many evaluations of
$\mathrm{EI}(\bm{u}_{i,g})$ (as shown also in the main text) it is a valid
algorithm for finding the maximum within the flat landscape of the acquisition
function.

##  Appendix D Other optimizers

Basin-hopping – Basin-hopping is a global stochastic optimization algorithm
[31]. It combines two steps: (i) a local optimization which proposes a
candidate solution and (ii) a perturbation of such candidate in order to make
it hop to other basins which might contain a global optimal point. The new
point is accepted or rejected according to a probability which depends on a
“temperature” parameter. The “temperature” parameter decreases with the
iteration number so that, at the beginning, new proposals are easily accepted
while, at larger iterations, the algorithms becomes more and more selective.
The algorithm runs for a fixed number of iterations.

Dual annealing – This global optimization algorithm is the generalized form of
the simulated annealing and it is paired with a local optimization which is
performed at the end of the annealing to refine the solution [32]. It is a
variation of a hill climbing algorithm in which a solution is randomly
perturbed and the new proposed point is accepted with a probability that
depends on the difference in energy between the two points. This probability
also depends on a “temperature” parameter that, like in the basin-hopping
case, decreases with the number of iterations in order to converge to a
candidate solution.