# Fleming-Viot helps speed up variational quantum algorithms in the presence
of barren plateaus

Daniel Mastropietro  CNRS-IRIT, Université de Toulouse INP, 31071, Toulouse,
France  Georgios Korpas  HSBC Lab, Innovation & Ventures, London, HSBC, United
Kingdom  Department of Computer Science, Czech Technical University in Prague,
Czech Republic  Vyacheslav Kungurtsev  Department of Computer Science, Czech
Technical University in Prague, Czech Republic  Jakub Marecek  Department of
Computer Science, Czech Technical University in Prague, Czech Republic

###### Abstract

Inspired by the Fleming-Viot stochastic process, we propose a variant of
Variational Quantum Algorithms benefitting from a parallel implementation of
the classical step of learning parameters of a given variational form, with
the aim of avoiding regions of the parameter space known as barren plateaus.
In the Fleming-Viot tradition, parallel searches are called particles. In our
proposed approach, the search by a Fleming-Viot particle is stopped when it
encounters a region where the gradient is too small or noisy. The stopped
particle continues the search after being regenerated at another potentially
more interesting location of the parameter space, biasing the exploration away
from barren plateaus. We analyze the behavior of the Fleming-Viot particles
from a theoretical standpoint, backed up with numerical experiments on
synthetic problems as well as on selected instances of the Max-Cut problem on
graphs, which show that our method performs better than plain-vanilla variants
when there are large barren plateaus.

##  1 Introduction

In both academia and industry, there has been considerable interest in solving
optimization problems [1] using near-term quantum computers [2], albeit often
in the form of quantum annealers [3, 4, 5, 6, 7, 8, 9, e.g.], or in the form
of variational quantum algorithms (VQAs) on gate-based quantum computers [10,
11, 12, 13, 14, e.g.] or some combination of the two [15, e.g.]. VQAs [11, 13]
are sometimes suggested as potential candidates to demonstrate an early
quantum advantage in the near-term noisy quantum computers due to their
simplicity, robustness to noise, and because the resources required are shared
between a quantum computer and a classical optimizer. So far, VQAs have shown
some potential in applications in a number of domains, including chemistry
[16, 17, 18, 19, 20, 21], strongly correlated systems in condensed matter
physics [22, 23, 24, 25, 26], combinatorial optimization [10, 27, 28, 29, 30,
31, 32, 33, 34], solving linear and nonlinear equations [35, 36, 37],
supervised quantum machine learning [38, 39], generative models [40, 41, 42],
and quantum neural networks [43, 44, 45, 46], to name a few.

Figure 1: The VQA schema. In dotted square the classical-quantum feedback
loop. The circuit is initialized with parameters $\vartheta$ with the purpose
to obtain the empirical estimate of the expectation value of a set of
observables $\\{O_{k}\\}_{k\in\mathbb{N}}$. In the classical part, a classical
optimizer minimizes the cost function in order to obtain a new set of
parameters, $\vartheta^{\prime}$ to be fed back to the quantum circuit.

Variational quantum algorithms work by first running a parameterized quantum
circuit that encodes the problem of interest, measuring the expected value of
some observable relevant to this problem, classically learning the parameters
that minimize a cost function related to the expected value of the observable,
and repeating until the process converges to some satisfactory value. See Fig.
1.

###  1.1 Challenges encountered in VQAs

Despite the potential of variational quantum algorithms, there exist several
challenges in their deployment to practical problems. Concretely, the
optimization problems involved are usually highly non-convex, and in general
it is expected that they require super-exponential time (in the number of
variables) to find the global optimum. Furthermore, in their practical
implementation, the classical optimization subroutine is done with gradient-
based algorithms for which, unless some underlying tame structure exists [47,
48, 49, 50, 51], in general there is no guarantee that such algorithms
converge to a global minimum whatsoever. The algorithm may get trapped in
(suboptimal) local minima or saddle points. The iteration complexity as well
as the impact of the biased noise present in NISQ devices are studied in [52].
Further convergence results are presented in [53].

One major drawback of variational quantum algorithms is that they suffer from
the so-called barren plateau problem [54, 55, 56, 57, 58, 59] especially when
the problem is approached with hardware efficient ansätze (that is, initial
solution guesses $\vartheta$ of the problem). Barren plateaus are large
regions in the parameter space where the gradient of the cost function
vanishes or is noisy around zero, which, in turn, makes efficient optimization
very difficult. Concretely, most gradient descent algorithms work as

$\displaystyle\vartheta_{\ell}(t+1)-\vartheta_{\ell}(t)=-\eta\frac{\partial
C}{\partial\vartheta_{\ell}},$

where $t$ is the time step of gradient descent dynamics, $\eta$ is the
learning rate, and $C$ is the cost function. The observation of Ref. [54] is
that, if the chosen variational circuit is highly random, the gradient of the
cost function is generically suppressed by the dimension of the Hilbert space,
and the variation of the cost function during gradient descent is very small,
that is, $\delta C\equiv\frac{C(t+1)-C(t)}{C(t)}\ll 1$ for step $t$. Recently,
Ref. [60] introduced a complete Lie-algebraic characterization of barren
plateaus and proved the dependence of the vanishing gradients on the dimension
of the underlying Hilbert space, that is

$\displaystyle{\rm Var}(\partial_{\vartheta_{\ell}}C)=\mathcal{O}(2^{-n}),$

where $n$ is the number of qubits.

In general, the problem of barren plateaus has also been argued to be related
to the choice of cost function [61], circuit entanglement [62, 63], circuit
parameterization expressivity [64], the presence of Pauli noise [55] while it
has been suggested that gradient-free methods such as the finite difference
stochastic approximation algorithm (FDSA) and the simultaneous perturbation
stochastic approximation algorithm (SPSA) may also suffer from the problem of
barren plateaus [57]. It should be noted that, within quantum machine learning
applications, similar to deep learning, one wants to achieve a high level of
expressivity. Recent work in Ref. [65] has shown that highly expressive
parameterized circuits give rise both to standard as we all as noise-induced
barren plateaus. Therefore, we conclude that this problem can hinder the
possibility of observing near-term advantage with variational quantum
algorithms.

Several methods have been devised to avoid barren plateaus. For example, [66]
uses the technique of classical shadows [67]. Other approaches have been
proposed in several works, for example, [68, 69]. Worth mentioning is the
approach of Ref. [70] where the authors construct certain conic extensions of
the parameterized quantum circuits allowing for jumps out of barren plateaus.
Other methods include certain parameter initialization strategies [71],
gradients via parameter correlation [72], ansatz classical neural network pre-
training [73], and layer-wise training [74].

Before we dive into our proposal, we would like to motivate our ideas by
considering a problem that appears to be quite similar. Consider random walks
on graphs. Despite their natural appeal for graph exploration and numerous
applications in distributed systems, random walks exhibit a significant
limitation: slow mixing. For example, in a simple ring topology with $n$
nodes, a random walk requires an average of $\mathcal{O}(n^{2})$ steps to
cover the entire graph, whereas a deterministic traversal can achieve the same
in just $n$ steps. The cover time, which represents the expected duration for
a random walk to visit every node at least once, emerges as a critical metric,
emphasizing the challenges posed by random walks. Ref. [75] suggested a way to
address the latency of random walks by utilizing multiple simultaneous random
walks. The hypothesis is that by initiating $k$ random walks concurrently, the
graph coverage can be expedited, reducing the cover time. It has been
empirically shown that a spectrum of outcomes, ranging from logarithmic speed-
ups in cases like ring graphs to exponential speed-ups in structures like the
barbell graph, can be obtained, while for a wide range of interesting graphs,
a linear speed-up in cover time has been observed as long as $k$ remained
within a permissible range.

Interestingly, the above limitation is somewhat related to the training of a
variational quantum algorithm as seen from the point of view of a non-finite
unbounded graph. In particular, the parameter landscape is analogous to
intricate graph structures, where random walks struggle with the extra
difficulty of the vanishing gradients.

###  1.2 Our contribution

The usage of multiple parallel random walks to expedite coverage can inspire
similar potentially heuristic strategies in this context and to some extent
our proposal follows along this line; as we will explain in detail in Sec. 2
by employing multiple parallel gradient descent initializations, with a
specific killing and regeneration rule, the so-called Fleming-Viot inspired
algorithm, described in Alg. 1, seems to better navigate the parameter
landscape, efficiently circumventing barren plateaus by biasing the paths
towards other locations that avoid them.

The rest of the paper is organized as follows: in Sec. 2 we introduce the
Fleming-Viot process and describe how we can apply its principles to the
classical optimization step of VQAs, which gives rise to Alg. 1. In Sec. 3 we
provide and prove a theorem on the reduced expected hitting time of our method
against vanilla simulated annealing. In Sec. 4 we showcase an application of
our Fleming-Viot inspired algorithm for a synthetic problem with a large
percentage of barren plateaus in the parameter landscape as well as an
application on the QAOA algorithm for an instance of the Max-Cut problem with
8 nodes. We summarize and conclude in Sec. 5.

##  2 Fleming-Viot Process within VQAs

The Fleming-Viot stochastic process (FV) was initially introduced by Fleming
and Viot [76], to model biological evolution. In the 2000s, several works
studied the relationship between the Fleming-Viot process and Markov
processes, starting with Burdzy _et al._ [77] who proved it can be used to
approximate the eigenfunctions of the Laplacian (the generator of a Brownian
motion) in a bounded domain. More recently, [78] proved this approximation is
valid for diffusion processes that are more general than the Brownian motion,
in both bounded and unbounded domains. The FV process was originally defined
as a continuous-time process, but it can also be defined as a discrete-time
process, see [79], which is the context that will be utilized in our work.

In its original definition, the Fleming-Viot process is a continuous-time
stochastic process consisting of $N$ particles that initially evolve
independently of each other on a state space $\Lambda\cup\\{0\\}$. Each
particle follows the dynamics of an underlying Markov process as long as it
does not visit $0$, an absorption state, i.e., where the underlying process
stays forever or, in other words, is killed. The FV process is defined so that
when a particle is killed, one of the remaining particles, chosen uniformly at
random, is split into two particles, each at the same state in $\Lambda$ of
the particle before the split. After the split, each particle evolves
independently of each other, as before the split.

In the context of the present article, the Fleming-Viot process is used as a
driving idea of a parallelized optimization process. Each particle represents
a different instance of a stochastic optimizer of choice (e.g., SPSA or
Simulated Annealing) of a parameter vector of interest, $\vartheta$, whose
values evolve following the rules of the selected optimizer. Each particle
initially evolves independently of each other until it is considered killed,
i.e. when the corresponding $\vartheta$ reaches a subset of the parameter
space considered not to be interesting for further exploration by the
optimization process. At this point, two things can happen to the killed
particle, which correspond to two different regeneration strategies, called
exploitation and exploration:

  * •

Under the exploitation strategy, the particle is regenerated to the position
of another particle chosen uniformly at random among the particles not
currently killed –which are expected to be in a more interesting location than
the killed particle.

  * •

Under the exploration strategy, the particle is regenerated to a randomly
chosen point in the parameter space.

Figs. 3 and 4 illustrate this. For every killed particle, the choice of one of
these options is made with a fixed probability defined by the so-called
“exploration rate” parameter. Notice that in the first case, and thanks to the
stochastic nature of the optimizer, the two particles initially at the same
position evolve independently of each other after the regeneration.

Overall, in the VQA context, the FV-based parallel optimizer is used for the
classical optimization part with a “killing" rule on the average magnitude of
the gradient in a recent window of optimization steps.

Figure 2: Figure (a) denotes the standard quantum-classical loop of VQAs as
shown in more detail in Fig. 1. The output of the parameterized circuit
$U(\vartheta)$ is input to the cost function $C(\vartheta)$ which is optimized
with a classical optimizer, e.g., SPSA, by iterating $N$ times. The updated
output is fed into the circuit for a new run. Figure (b) denotes the
parallelization of the classical part of the VQA with the difference being
that the iterates of cost function, whose optimizers are called particles
$p=[p_{1},\ldots,p_{n}]$, form essentially a parallelized optimization
procedure where each $p_{i}$ iterates $M$ times such that $nM=N$.

The details of the FV-based parallel optimization algorithm are defined in
Alg. 1.

##  3 Theoretical Analysis

Our main result is a hitting-time analysis for Alg. 1. In the next section
(Sec. 4), the analysis is corroborated by numerical results. To present the
analysis, let us define several quantities that we will refer to:

  * •

$x$ will denote the optimization variable, i.e. $\vartheta$ in Figure 1. We
shall denote $X_{k}$ the random process associated with the iterates of the
two algorithms, Simulated Annealing with Fleming-Viot (SA+FV) and without
(SA).

  * •

$D$ is the (Lebesgue) measure of the domain for the optimization variable $x$.

  * •

$f(x)$ is the objective function with Lipschitz-smooth gradient.

We denote the equation corresponding to SA as:

$X_{k+1}=X_{k}-\eta_{k}g_{k}+\sqrt{\eta_{k}}\epsilon_{k}n_{k}$ (1)

where $g_{k}\approx\nabla f(x_{k})$ is a noisy estimate of the gradient,
$\eta_{k}$ is the learning rate, $\epsilon_{k}$ is the temperature defining
the noise level of the SA algorithm, and $n_{k}$ is a standard Gaussian random
variable. We are interested, in particular, with the impact of the barren
plateaus on the hitting time of SA (1), where:

######  Definition 1 (Barrren plateau)

A ball in the space of $x$ with center $\hat{x}$ and radius $R$ is an
$\alpha$-barrren plateau of function $f$ whenever $\|\nabla f(x)\|\leq\alpha$
for all $x\in\mathbb{B}(\hat{x},R)$.

Furthermore, we assume:

######  Assumption 1 (Stationary points and regions)

When we observe the stochastic process (1) arriving at
$\|X_{k+1}-X_{k}\|_{2}^{2}\leq\beta$, we assume that for some predefined small
positive constant $\beta$, this indicates that $X_{k}$ satisfies $\|\nabla
f(X_{k})\|\leq\alpha$. Furthermore we assume that this is a local minimizer
with probability $l/s$, a saddle point with probability $a/s$, and an
$\alpha$-barrren plateau with probability $b/s$. Simplistically, this could be
seen as having $s=l+b+a$ $\beta$-stationary points or regions, out of which
there are $l$ local minimizers, one of which is global, $b$ plateaus that are
$\alpha$-barrren, and $a$ saddle points. At the same time, consider that a
fraction $0<B<1$ of the space is barren, i.e. $BD$ is the total measure of a
barren plateau, and that each $\alpha$-barrren plateau is a ball of equal
size, i.e., of $BD/b$ measure, from which we can compute the radius to be
proportional to $(BD/b)^{1/d}$, with $d$ the dimension of $x$.

When $x$ is in an $\alpha$-barren plateau, the noise in the iteration
dominates the gradient term. Thus, we can treat the gradient as insignificant,
and so we simplistically model the processes.

######  Assumption 2 (Behaviour in an $\alpha$-barren plateau)

We assume that with FV, identification of barrenness is performed with full
accuracy. Thus, upon reaching an $\alpha$-barren plateau:

  1. 1.

SA+FV immediately identifies the plateau and performs a jump escaping the
plateau.

  2. 2.

SA behaves as a Gaussian random walk, being defined entirely by the noise
component of the iteration.

Then:

######  Theorem 1 (Hitting-time analysis)

Under Assumptions 1–2, the expected time of SA+FV to hit a global minimizer is
given by,

$C_{s}D\alpha^{-4}+(s-1)\left[C_{s}D\alpha^{-4}+\frac{aC_{a}}{s\alpha}+\frac{(l-1)C_{l}}{s\alpha}\right]$

and the expected hitting time of SA (1) without FV is,

$\begin{array}[]{l}C_{s}D\alpha^{-4}+(s-1)\left[C_{s}D\alpha^{-4}+\frac{aC_{a}}{s\alpha}+\frac{(l-1)C_{l}}{s\alpha}\right.\\\
\qquad\qquad\left.+\frac{C_{f}BD}{sb}\right],\end{array}$

where the constants $C_{s},C_{a},C_{l},C_{f}$ scale the escape time for the
stationary points of Assumptions 1, that is, stationary points, saddle points,
local minimizers, and barren plateaus.

Proof. Consider a process $X_{k}$ driven by SA with FV and $\hat{X}_{k}$
driven by SA without it. Consider the beginning of an event as the
initialization of SA, at $k=0$, at some random non-barren point.

With $\alpha$ the gradient threshold and a step size annealing rate of
$n^{1/2}$, the hitting time to a stationary point is $C_{s}D\alpha^{-4}$ [80,
Theorem 3.3]. We consider that generally there is an equal probability of the
iterate hitting any one of the $s$ stationary points. Let us consider each
type of point and its consequences on the two algorithms’ iterations.

With probability $1/s$, the process has reached a global minimum.

Otherwise, with probability $a/s$, both algorithms take $C_{a}\alpha^{-1}$
[80, Lemma 4.2] to escape a saddle point. Subsequently the next iterated event
begins, i.e., we have another SA path from an agnostic starting point
proceeding until it hits a stationary point.

Alternatively, with probability $(l-1)/s$, both algorithms take
$C_{l}\alpha^{-1}$ [81, Corollary 6] to escape a local minimum. Subsequently,
the next iterated event begins.

Finally, with probability $b/s$, the algorithm reaches a barren plateau. At
this point FV escapes, reinitializing the next event (to a randomly chosen
point). However, without FV, the time the algorithm remains at the barren
plateau has the expected value of the sojourn time of a Gaussian random walk
in a convex body. The expected sojourn time is given by $C_{f}R^{d}$ [82].
Subsequently, the event restarts.

Now we sum these quantities. The final run, that is, the event at which the
stationary point reached is a global minimizer, lasts for as long as the
algorithm iterates reach a stationary point (which in this case is the global
minimum).

The expected number of stationary points traversed before a global minimizer
is reached is a geometric random variable with parameter $1/s$. Thus, for an
expected $s-1$ time steps, we hit the other stationary points, and this
enables us to compute the expected hitting time as the $s-1$ sum of the
probabilities of hitting the different types of stationary points and their
expected contribution to the total iteration time.

This means that:

###### Corollary 1

Compared to the expected hitting time of SA, the expected time of SA+FV to hit
a global minimizer is sped up by a total expected number of iterations given
by

$\frac{C_{f}BD(s-1)}{bs}$

i.e., the speed up is an additive term that is linear with respect to both the
number of separated barren plateaus as well as the portion of the space that
is barren.

Note that this appears to be independent of $\alpha$. However, this is
misleading, as we expect the definition of the threshold to affect $B$ and
$b$, thus the more relaxed the condition, the larger $\alpha$, the larger the
region that is labeled “barren”.

Figure 3: Particle reinitialization. The particle is reinitialized to a random
point according to some probability distribution $\mathcal{D}$. Above,
particle $p_{3}$ is absorbed and reinitialized at any of the “red” spots in
the non-barren space according to $\mathcal{D}$ (e.g., uniform distribution).
Note that $\alpha$ is computed by a burning period. Figure 4: Particle
reactivation amounts to absorbing a particle that “arrived” at a region in
which the value of the gradient is below the threshold $\alpha$ and
regeneration, according to some probability distribution $\mathcal{D}$, at the
position of another active particle. Above, particle $p_{3}$ is being
regenerated at either $i_{n}^{p_{1}}$ or $i_{n}^{p_{2}}$. Our experiments
showed that this might not be the preferable method, though, since usually one
particle suffices to explore an interesting region.

##  4 Numerical Experiments

We analyzed the FV-based parallel optimizer using simulated annealing (SA)
[83, 84, 85] as each particle’s optimizer, because SA gives a good balance
between simplicity and power. The experiments were run against a benchmark
defined by using SA as a single optimizer on as many ansätze (initial
parameter values) as the number of particles used in the FV-based parallel
optimizer. The values of the ansätze in the SA run and in the FV-based run are
the same, so that both optimizers start at the same locations in the parameter
space.

We carried out two sets of experiments, each under a different optimization
scenario on a bounded domain; concretely:

  1. (_i_)

a two-dimensional optimization on synthetic functions generated with a
predefined percentage of barren areas;

  2. (_ii_)

a multidimensional optimization using VQA on the Max-Cut problem approximated
by a noiseless quantum circuit.

During the search for the optimum parameters, any parameter dimension whose
value falls outside the domain in the next iterated value is clipped to the
border of the domain in that dimension. The SA parameter in Eq. (1) defining
the temperature is set to $\epsilon_{k}=0.1$, $\forall k$; the learning rate
$\eta_{k}$ is initially selected from an analysis of the function’s landscape
in the neighbourhood of the anszatz and then decayed with a $\frac{1}{n}$
power law, in order to satisfy the convergence conditions of the stochastic
approximation algorithm on which SA is based.

Each set of experiments is defined by a parameterization specific to the
optimization problem to be solved. The results are collected and summarized
over 10-replication runs on each experiment.

Of particular interest is the analysis of the effect of the exploitation vs.
exploration strategies, which are controlled by the exploration rate
parameter. Thus, in each of the two optimization scenarios, two exploration
rates were considered: $50\%$ and $100\%$. Note that, in the latter case, all
killed particles are regenerated under the exploration strategy described in
Sec. 2.

###  4.1 Synthetic Experiments

Each synthetic experiment is defined by an expected percentage of barren
plateaus in the optimization landscape whose domain is set as the two-
dimensional rectangular area $[0,1]\times[0,1]$. The synthetic function is
constructed as the output of a smoothing spline interpolation on a set of
$(x,y,z)$ points, using the procedure described in App. A.

(a)

(b)

Figure 5: Violin plots on 10 replications of the relative error (log scale) of
the optimum value found for the synthetic experiments. Red: SA; green: FV. The
number of optimization steps run per particle is 50. Lower position of the
violins is better. More details are given in the text.

Parameter | Value Range | Value Used
---|---|---
Domain | Bounded in $\mathbb{R}^{2}$ | $[0,1]\times[0,1]$
Barren % | 0%-100% | 70%
| Gradient absorption
---
value
$\geq 0$ | 1
| Exploration
---
rate
0%-100% | 50% and 100%
| Number of
---
particles
$\geq 2$ | 10

Table 1: Main parameters used in the synthetic experiments

Tab. 1 lists the most important simulation parameters for the synthetic
experiments. Fig. 5 shows the results of the experiments carried out on 10
synthetic functions for two different exploration rates: $50\%$ and $100\%$.
Concretely, we have produced the violin plots on 10 replications of the
relative error of the optimum value found for 10 synthetic functions designed
with $\sim$ 70% of barren areas. To increase the ability to characterise the
aspects of the optimization functions that might be more favorable to a
successful application of the FV strategy, compared to the plain SA algorithm,
we have chosen the ansätze uniformly at random among the grid points of the
nominal barren areas (see the function synthesis process in App. A), as long
as their distance to the global optimum is larger than $0.6$, to try to avoid
any lucky choice of the ansatz that would help the algorithm find the global
optimum quickly. Each panel corresponds to a different exploration rate. In
each panel, violin plots are paired by the respective synthetic function, and
those pairs are sorted by decreasing observed advantage of FV over SA111The
advantage of FV over SA is indicated in parentheses below each group labels on
the horizontal axis. It is calculated as the difference of the paired medians
relative to the worst possible error, namely the error obtained by using the
function’s maximum value as the estimate of the global minimum. Such worst
possible error is marked as a horizontal gray line at the top of each pair..
The labels below each group (e.g. “3-SA, 3-FV”) indicate the function ID and
the method associated to the results, which enables a comparison of the impact
of the two exploration rates on the results for a particular function. In red,
we denote simulated annealing, while in green we denote our Fleming-Viot
inspired process. In these experiments we chose the number of optimization
steps run per particle to be 50. We observe that in general Fleming-Viot
performs favorably against vanilla SA, but its advantage strength over SA
depends on the function being optimised. We observe that the top-5 and
bottom-5 functions are the same in the $50\%$ and in the $100\%$ exploration
rates; but the $100\%$ exploration rate produces clearly larger advantage of
FV over SA than the $50\%$ exploration rate for the five top-ranked functions
(IDs $3,5,6,7,$ and $9$); on the contrary, the advantage is very similar for
the bottom-ranked functions, with IDs $1,2,4,8,$ and $10$. This suggests that
the exploration strategy has a larger impact on the algorithm’s performance
than the exploitation strategy.

(a)

(b)

Figure 6: The top-ranked (a) and bottom-ranked (b) synthetic functions in the
results of 5(b) corresponding to the case of the exploration rate 100%. In
each panel, the 3D landscape is shown on the left, and the corresponding
contour plot is shown on the right. Shaded squared areas are superimposed on
the right plot which correspond to the synthesized nominal (i.e. before
smoothing) barren plateaus, whose colour varies from blue (lower nominal
function value) to green (higher nominal function value); the red star
indicates the location of the global minimum, and the red points show the 100
ansätze ($=10$ particles $\times$ $10$ replications).

Fig. 6 shows the landscape of the top-ranked function and of the bottom-ranked
function in the results of the $100\%$ exploration rate case. Although the
function characteristics that are more favourable to the FV strategy are not
straightforward identifiable, we observe that the landscape of the function
where FV most outperforms SA (i.e. function ID 9 for which, in median, the
relative error of the estimated optimum is reduced by $33.9\%$ of the worst
possible error thanks to FV) corresponds to a landscape where the global
minimum is located at a quite steep valley with its minimum at $\sim-20$. On
the other end, FV performs similarly to SA ($0\%$ advantage) on a function (ID
2) that has a large positive peak of $\sim 20$ while the global minimum takes
a value close to 0, i.e. similar to the value of the function at the barren
plateaus.

###  4.2 QAOA Experiments

The Quantum Approximate Optimization Algorithm (QAOA) [10] (see also [86] for
a survey) is a VQA proposed to solve combinatorial optimization problems.
Specifically, in QAOA two variational parameter vectors are utilised,
$\beta,\gamma\in\mathbb{R}^{L}$, and $L$ layers of the products of the mixer-
unitary $U_{j}(\beta)=e^{-i\beta_{j}X}$, wth $\beta_{j}\in[0,\pi]$, and
problem unitary $U(\gamma)=e^{-i\gamma_{j}C}$, with $\gamma_{j}\in[0,2\pi]$
are applied per layer., where $i$ is the imaginary unit. Here
$X=\sum_{j}X_{j}$ and $X_{j}$ denotes the Pauli-$X$ gate on qubit $j$ and $C$
is the cost function of the encoded problem.

Parameter | Value Used
---|---
Number of qubits | 8
Level of connectedness | 100%
Weights on graph edges | Uniform distribution
Number of layers $L$ | 1
Gradient absorption value | 1.0
Exploration rate | 50% and 100%
Number of particles $p$ | 10

Table 2: Main parameters used in the Max-Cut experiments.

Then QAOA produces the state

$\left|\psi_{\gamma,\beta}\right\rangle=U_{L}(\beta_{L})U_{L}(\gamma_{L})\ldots
U_{1}(\beta_{1})U_{1}(\gamma_{1})|+\rangle.$

An empirical estimate $\widehat{C}$ of the expectation value of the cost
function $\langle\psi_{\gamma,\beta}|C|\psi_{\gamma,\beta}\rangle$ is then
measured on a specified number of shots run on the quantum circuit, and the
goal is to classically optimize the variational parameters and repeat so as to
converge to an optimal value. Due to the periodicy of the unitaries, the
search for the optimum $\beta$ and $\gamma$ values happens in the hyper-cube
with side $[0,\pi]$ on the $\beta$ dimensions and side $[0,2\pi]$ on the
$\gamma$ dimensions.

Fig. 7 shows the graph used in the problem and Fig. 8 presents the relative
error of the maximum value found by the two optimization algorithms under
comparison: SA and FV.

Figure 7: Fully-connected 8-node weighted graph of the analyzed Max-Cut
problem, where the edge width is proportional to the weight value which
follows a uniform distribution in $[0,1]$.

Concretely, just like in the previous example with the synthetic experiment,
we present violin plots on 10 replications of the relative error of the
optimum value found for the Max-Cut problem on a fully-connected 8-node graph
with edge weights chosen uniformly at random between 0 and 1. The results are
shown for 10 different test graphs and each panel corresponds to a different
exploration rate. In each panel, the violin plots are paired by the respective
optimization function and, as in the synthetic experiments presented in Sec.
4.1, those pairs are sorted by decreasing observed advantage of FV over SA, as
defined in footnote 1. The labels below each group (e.g. “3-SA, 3-FV”)
indicate the function ID associated to the optimization problem and the method
used, so that it is possible to compare the results on a particular function
obtained by the two exploration rates. As before, in red we denote the SA
results and in green the FV results. The number of shots for the evaluation of
the quantum circuit is 512 and the number of optimization steps run per
particle is 50.

As opposed to the synthetic experiments, in the Max-Cut experiments the
advantage of FV over SA takes similar values in the two analyzed exploration
rates across the analyzed functions, whose ranking is almost the same in both
cases. However, in the $100\%$ exploration rate, we observe a reduced spread
of the violin plots in the FV approach on the larger error end, suggesting
again that exploration is more important than exploitation.

(a)

(b)

Figure 8: Violin plots on 10 replications of the relative error of the optimum
value found for the Max-Cut problem on a fully-connected 8-node graph with
edge weights chosen uniformly at random between 0 and 1. Red: SA, green: FV.
Lower position of the violins is better.

##  5 Summary and conclusions

We have explored a novel approach to enhancing variational quantum algorithms,
while taking inspiration from the Fleming-Viot stochastic process,
traditionally applied in modeling biological evolution. Our primary objective
was to address the challenge of barren plateaus in variational quantum
algorithms – regions in the parameter space where gradients vanish or are
excessively noisy, impeding effective optimization.

We proposed a parallel implementation of the classical learning step in
variational quantum algorithms, utilizing the concept of multiple searches
(termed particles) that are terminated upon encountering barren plateaus and
then regenerated (using either an exploitation or an exploration strategy)
into potentially more fruitful areas of the parameter space. This strategy
biases the search for the global optimum away from barren plateaus, in order
to improve the overall efficiency of variational quantum algorithms.

Our methodology was tested through numerical experiments on synthetic problems
and specific instances of the Max-Cut problem on small graphs. The results
demonstrated that, in certain situations, our approach outperforms traditional
simulated annealing methods. They also suggest that the exploration strategy
for particle regeneration is more important than the exploitation strategy in
terms of algorithm performance. We believe that this approach not only
enhances the optimization process but also offers a fresh perspective on
tackling the inherent challenges in quantum computing optimizations.

Our research contributes to the ongoing efforts in quantum computing to
realize the potential of variational quantum algorithms for near-term
advantage, especially in the face of obstacles like barren plateaus. The
insights gained from our work open new research directions in near-term
quantum algorithms, particularly in developing more robust and efficient
classical optimization strategies for variational quantum algorithms.

#### Acknowledgements

We would like to thank Matthieu Jonchkeere, from CNRS-LAAS in Toulouse,
France, for the insights received during the development of the algorithm.
Jakub would like to acknowledge many discussions with Andrew J. Parkes up to
2012 concerning landscapes and algorithm analysis.

#### Disclaimer

This paper was prepared for information purposes, and is not a product of HSBC
Bank Plc. or its affiliates. Neither HSBC Bank Plc. nor any of its affiliates
make any explicit or implied representation or warranty and none of them
accept any liability in connection with this paper, including, but limited to,
the completeness, accuracy, reliability of information contained herein and
the potential legal, compliance, tax or accounting effects thereof. Copyright
HSBC Group 2023.

Algorithm 1 FV-based parallel optimization algorithm

1:$T$: the number of classical optimization iterations, $N$: the number of
Fleming-Viot particles, $B$: the burn-in time steps used to compute the
reference gradient norm involved in the particle’s absorption condition, $W$:
the size of the time window on which the average gradient is computed,
$\alpha$: the relative gradient threshold for absorption of a particle,
$\varepsilon$: the exploration rate.

2:Estimate of the parameter vector minimizing the system energy:

3:$\text{avg\\_abs\\_gradient}\leftarrow 0$

4:for $i=1$ to $B$ do

5: OptimizeAndRecordGradient

6: Update avg_abs_gradient

7:end for

8:$\text{ref\\_gradient}\leftarrow\alpha\times\text{avg\\_abs\\_gradient}$

9:for $i=1$ to $T$ do

10: OptimizeAndRecordGradient

11: CheckAbsorptionAndRegenerateParticles

12:end for

13:procedure OptimizeAndRecordGradient

14: for $j=1$ to $N$ do

15: Perform optimization step for particle $j$

16: Record observed gradient for particle $j$

17: end for

18:end procedure

19:procedure CheckAbsorptionAndRegenerateParticles

20: for $j=1$ to $N$ do

21: if avg gradient of particle $j$ over last $W$ steps < ref_gradient then

22: $\text{rnd}\sim\mathcal{U}[0,1]$

23: if $\text{rnd}<\epsilon$ then

24: ReinitializeParticle(j)

25: else

26: ReactivateParticle(j)

27: end if

28: end if

29: end for

30:end procedure

31:procedure ReinitializeParticle(j)

32: Regenerate particle $j$ to a randomly chosen point in the parameter space

33: Reset the history of recorded gradient values for particle $j$

34:end procedure

35:procedure ReactivateParticle(j)

36: Regenerate particle $j$ to a randomly chosen particle ($r$) among non-
killed particles

37: Assign the last $W$ gradient values of particle $r$ as the last $W$
gradient values of particle $j$

38: Reset the learning rate of particle $j$ to the initial learning rate of
particle $r$ with the goal of increasing the divergence between the two
particles, thus boosting exploration

39:end procedure

##  Appendix A Procedure used to generate the synthetic functions test bench

The following procedure was used to generate the synthetic functions used as
test bench of the proposed optimization algorithm. Given:

  * •

the 2D rectangular parameter domain on which the synthetic function will be
generated (e.g. $[0,1]\times[0,1]$),

  * •

the number of grid points to consider in each dimension on which the 3D
$(x,y,z)$ point values are proposed for function $f$ to be synthesized (e.g.
100 points),

  * •

the nominal percentage of barren plateaus,

  * •

the nominal width (number of grid points), in each parameter dimension, of
each barren plateau to be generated,

  * •

the width of the border margin, i.e. the margin to leave between the border of
the domain and the barren plateaus (e.g. 2 grid points),

  * •

the smoothing parameter for the smoothing spline interpolation that is used to
define the final value of function $f$ in each grid point (e.g. 1),

the synthetic process consists of roughly the following iterative steps:

  1. (1)

Set the value of the function at the grid points in the border margin to a
value a little larger than 1 (e.g. 2), so that the global minimum of the
function is less likely to be located at the border, which is an unwanted
situation222The situation is unwanted because it artificially helps the
optimization algorithm find the global optimum. This aid is due to the fact
that the parameter estimate at each optimization step is clipped to the
domain. Therefore, if after a particular step the parameter estimate happens
to fall outside the domain, its location might end up being close to the
global optimum just because of clipping..

  2. (2)

Randomly sample a point among the grid points in a subset of the domain that
guarantees there is no overlap between the barren plateaus and the border
margin; this point represents the center of a barren plateau.

  3. (3)

Construct a square centered at the sampled point, which represents the new
proposed barren plateau.

  4. (4)

Sample a value from a uniformly distributed random variable in [0, 1] that
represents the proposed constant value of the optimization function at each
grid point of the barren plateau.

  5. (5)

Check overlap between the newly constructed square and other barren plateaus
previously constructed. If there is overlap, regions are merged into a single
barren area. The function values assigned to each point in each of the regions
before merging is maintained.

  6. (6)

Compute the barren plateau percentage as the percent of grid points in the
barren areas constructed so far and the total number of grid points in the
domain. If this percentage is smaller than the nominal percentage of barren
plateau, go to step (2).

  7. (7)

Sample a value from a lognormal distribution with zero mean and a specified
scale parameter (e.g. 1). Set the negative of this value to the function value
at a randomly selected point in the grid, where the global minimum of the
function is expected to be located.

  8. (8)

Fit a 2D smoothing cubic spline (using the chosen smoothing parameter) to all
the sampled function values: (i) the values assigned to the grid points in all
generated barren plateaus, (ii) the value assigned to the grid points at the
four borders of the domain, (iii) the negative function value representing the
nominal global minimum. The smoothing spline fit is used to assign a function
value to each grid point in the domain. Note: due to the smoothing spline fit,
the location and value of the global minimum may not be the same as the
selected values in step (7). The size and location of the nominal barren areas
may also be affected by the smoothing interpolation. If the location of the
global minimum is at the border of the domain, in order to avoid the unwanted
situation described in footnote 2, the synthesized function is discarded and
the whole process described here is run again until the global minimum falls
in an interior point of the domain.

## References

  * [1] A. Abbas, A. Ambainis, B. Augustino, A. Bärtschi, H. Buhrman, C. Coffrin, G. Cortiana, V. Dunjko, D. J. Egger, B. G. Elmegreen, F. Fratini, B. Fuller, J. Gacon, C. Gonciulea, S. Gribling, S. Gupta, S. Hadfield, R. Heese, G. Kircher, T. Kleinert, T. Koch, G. Korpas, S. Lenk, J. Marecek, G. Mazzola, S. Mensa, N. Mohseni, G. Nannicini, C. Oâ€™Meara, E. P. Tapia, S. Pokutta, M. Proissl, P. Rebentrost, E. Sahin, B. C. B. Symons, S. Tornow, V. Valls, S. Woerner, M. L. Wolf-Bauwens, J. Yard, S. Yarkoni, D. Zechiel, S. Zhuk, and C. Zoufal, “Quantum optimization: Potential, challenges, and the path forward,” to appear, 2023.
  * [2] J. Preskill, “Quantum computing in the NISQ era and beyond,” Quantum, vol. 2, p. 79, aug 2018.
  * [3] T. Kadowaki, “Study of optimization problems by quantum annealing,” arXiv preprint quant-ph/0205020, 2002.
  * [4] A. Das and B. Chakrabarti, eds., Quantum Annealing and Other Optimization Methods.  Springer Berlin Heidelberg, 2005.
  * [5] S. Yarkoni, E. Raponi, T. Bäck, and S. Schmitt, “Quantum annealing for industry applications: Introduction and review,” Reports on Progress in Physics, 2022.
  * [6] V. Mehta, F. Jin, H. De Raedt, and K. Michielsen, “Quantum annealing with trigger hamiltonians: Application to 2-satisfiability and nonstoquastic problems,” Phys. Rev. A, vol. 104, no. 3, p. 032421, 2021.
  * [7] E. J. Crosson and D. A. Lidar, “Prospects for quantum enhancement with diabatic quantum annealing,” Nature Reviews Physics, vol. 3, p. 466â€“489, May 2021.
  * [8] B. K. Chakrabarti, H. Leschke, P. Ray, T. Shirai, and S. Tanaka, “Quantum annealing and computation: challenges and perspectives,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 381, Dec. 2022.
  * [9] S. Certo, G. Korpas, A. Vlasic, and P. Intallura, “Predicting Ising model performance on quantum annealers,” arXiv preprint arXiv:2311.07388, 2023\.
  * [10] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” arXiv preprint arXiv:1411.4028, 2014.
  * [11] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, “The theory of variational hybrid quantum-classical algorithms,” New Journal of Physics, vol. 18, no. 2, p. 023023, 2016.
  * [12] G. G. Guerreschi and A. Y. Matsuura, “Qaoa for max-cut requires hundreds of qubits for quantum speed-up,” Scientific Reports, vol. 9, May 2019.
  * [13] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, et al., “Variational quantum algorithms,” Nature Reviews Physics, vol. 3, no. 9, pp. 625–644, 2021\.
  * [14] D. J. Egger, J. Mareček, and S. Woerner, “Warm-starting quantum optimization,” Quantum, vol. 5, p. 479, 2021.
  * [15] S. H. Sack and M. Serbyn, “Quantum annealing initialization of the quantum approximate optimization algorithm,” Quantum, vol. 5, p. 491, 2021.
  * [16] P. J. Oâ€™Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, et al., “Scalable quantum simulation of molecular energies,” Phys. Rev. X, vol. 6, no. 3, p. 031007, 2016.
  * [17] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” Nature, vol. 549, no. 7671, pp. 242–246, 2017.
  * [18] J. Colless, V. Ramasesh, D. Dahlen, M. Blok, M. Kimchi-Schwartz, J. McClean, J. Carter, W. de Jong, and I. Siddiqi, “Computation of molecular spectra on a quantum processor with an error-resilient algorithm,” Phys. Rev. X, vol. 8, Feb. 2018.
  * [19] A. J. McCaskey, Z. P. Parks, J. Jakowski, S. V. Moore, T. D. Morris, T. S. Humble, and R. C. Pooser, “Quantum chemistry as a benchmark for near-term quantum computers,” npj Quantum Information, vol. 5, no. 1, p. 99, 2019\.
  * [20] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, “Quantum computational chemistry,” Reviews of Modern Physics, vol. 92, no. 1, p. 015003, 2020.
  * [21] Google AI Quantum and Collaborators, F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, et al., “Hartree-fock on a superconducting qubit quantum computer,” Science, vol. 369, no. 6507, pp. 1084–1089, 2020.
  * [22] C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, et al., “Self-verifying variational quantum simulation of lattice models,” Nature, vol. 569, no. 7756, pp. 355–360, 2019.
  * [23] C. Lyu, V. Montenegro, and A. Bayat, “Accelerated variational algorithms for digital quantum simulation of many-body ground states,” Quantum, vol. 4, p. 324, 2020.
  * [24] J. W. Z. Lau, K. Bharti, T. Haug, and L. C. Kwek, “Noisy intermediate scale quantum simulation of time dependent hamiltonians,” arXiv preprint arXiv:2101.07677, 2021.
  * [25] T. Haug and K. Bharti, “Generalized quantum assisted simulator,” Quantum Science and Technology, vol. 7, no. 4, p. 045019, 2022.
  * [26] M. Gong, H.-L. Huang, S. Wang, C. Guo, S. Li, Y. Wu, Q. Zhu, Y. Zhao, S. Guo, H. Qian, et al., “Quantum neuronal sensing of quantum many-body states on a 61-qubit programmable superconducting processor,” Science Bulletin, vol. 68, no. 9, pp. 906–912, 2023.
  * [27] E. Farhi and A. W. Harrow, “Quantum supremacy through the quantum approximate optimization algorithm,” arXiv preprint arXiv:1602.07674, 2016.
  * [28] Z. Wang, S. Hadfield, Z. Jiang, and E. G. Rieffel, “Quantum approximate optimization algorithm for maxcut: A fermionic view,” Phys. Rev. A, vol. 97, no. 2, p. 022304, 2018.
  * [29] G. E. Crooks, “Performance of the quantum approximate optimization algorithm on the maximum cut problem,” 2018.
  * [30] N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow, A. Cross, D. J. Egger, S. Filipp, A. Fuhrer, J. M. Gambetta, M. Ganzhorn, et al., “Quantum optimization using variational algorithms on near-term quantum devices,” Quantum Science and Technology, vol. 3, no. 3, p. 030503, 2018.
  * [31] M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, S. Boixo, et al., “Quantum approximate optimization of non-planar graph problems on a planar superconducting processor,” Nature Physics, vol. 17, no. 3, pp. 332–336, 2021.
  * [32] S. Brandhofer, D. Braun, V. Dehn, G. Hellstern, M. Hüls, Y. Ji, I. Polian, A. S. Bhatia, and T. Wellens, “Benchmarking the performance of portfolio optimization with qaoa,” Quantum Information Processing, vol. 22, no. 1, p. 25, 2022.
  * [33] S. Boulebnane and A. Montanaro, “Solving boolean satisfiability problems with the quantum approximate optimization algorithm,” 2022.  Presented at QIP 2023.
  * [34] L. Zhu, H. L. Tang, G. S. Barron, F. A. Calderon-Vargas, N. J. Mayhall, E. Barnes, and S. E. Economou, “Adaptive quantum approximate optimization algorithm for solving combinatorial problems on a quantum computer,” Phys. Rev. Res., vol. 4, p. 033029, Jul 2022.
  * [35] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, and P. J. Coles, “Variational quantum linear solver,” arXiv preprint arXiv:1909.05820, 2019\.
  * [36] X. Xu, J. Sun, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, “Variational algorithms for linear algebra,” Science Bulletin, vol. 66, no. 21, pp. 2181–2188, 2021.
  * [37] M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch, “Variational quantum algorithms for nonlinear problems,” Phys. Rev. A, vol. 101, p. 010301, Jan 2020.
  * [38] M. Schuld and N. Killoran, “Quantum machine learning in feature hilbert spaces,” Phys. Rev. Lett., vol. 122, no. 4, p. 040504, 2019.
  * [39] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, “Supervised learning with quantum-enhanced feature spaces,” Nature, vol. 567, no. 7747, pp. 209–212, 2019.
  * [40] P.-L. Dallaire-Demers and N. Killoran, “Quantum generative adversarial networks,” Phys. Rev. A, vol. 98, no. 1, p. 012324, 2018.
  * [41] M. Benedetti, D. Garcia-Pintos, O. Perdomo, V. Leyton-Ortega, Y. Nam, and A. Perdomo-Ortiz, “A generative modeling approach for benchmarking and training shallow quantum circuits,” npj Quantum Information, vol. 5, no. 1, p. 45, 2019.
  * [42] E. Y. Zhu, S. Johri, D. Bacon, M. Esencan, J. Kim, M. Muir, N. Murgai, J. Nguyen, N. Pisenti, A. Schouela, K. Sosnova, and K. Wright, “Generative quantum learning of joint probability distribution functions,” Phys. Rev. Res., vol. 4, p. 043092, Nov 2022.
  * [43] Y. Du, M.-H. Hsieh, T. Liu, and D. Tao, “Expressive power of parametrized quantum circuits,” Phys. Rev. Res., vol. 2, no. 3, p. 033125, 2020.
  * [44] J. Bausch, “Recurrent quantum neural networks,” in Advances in Neural Information Processing Systems (H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, eds.), vol. 33, pp. 1368–1379, Curran Associates, Inc., 2020.
  * [45] E. A. Cherrat, I. Kerenidis, N. Mathur, J. Landman, M. Strahm, and Y. Y. Li, “Quantum vision transformers,” arXiv preprint arXiv:2209.08167, 2022.
  * [46] X. Hou, G. Zhou, Q. Li, S. Jin, and X. Wang, “A duplication-free quantum neural network for universal approximation,” Science China Physics, Mechanics & Astronomy, vol. 66, no. 7, p. 270362, 2023.
  * [47] J. Bolte and E. Pauwels, “Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning,” Mathematical Programming, vol. 188, pp. 19–51, 2021.
  * [48] J. Bolte, T. Le, E. Pauwels, and T. Silveti-Falls, “Nonsmooth implicit differentiation for machine-learning and optimization,” in Advances in Neural Information Processing Systems (M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, eds.), vol. 34, pp. 13537–13549, Curran Associates, Inc., 2021.
  * [49] C. Aravanis, J. Aspman, G. Korpas, and J. Marecek, “Polynomial matrix inequalities within tame geometry,” arXiv preprint arXiv:2206.03941, 2022\.
  * [50] J. Aspmann, V. Kungurtsev, and R. R. Seraji, “Riemannian stochastic approximation for minimizing tame nonsmooth objective functions,” arXiv preprint arXiv:2302.00709, 2023.
  * [51] J. Aspman, G. Korpas, and J. Marecek, “Taming binarized neural networks and mixed-integer programs,” arXiv preprint arXiv:2310.04469, 2023.
  * [52] V. Kungurtsev, G. Korpas, J. Marecek, and E. Y. Zhu, “Iteration complexity of variational quantum algorithms,” arXiv preprint arXiv:2209.10615, 2022\.
  * [53] L. Binkowski, G. Koßmann, T. Ziegler, and R. Schwonnek, “Elementary proof of qaoa convergence,” arXiv preprint arXiv:2302.04968, 2023.
  * [54] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, “Barren plateaus in quantum neural network training landscapes,” Nature communications, vol. 9, no. 1, p. 4812, 2018.
  * [55] S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, “Noise-induced barren plateaus in variational quantum algorithms,” Nature communications, vol. 12, no. 1, pp. 1–11, 2021.
  * [56] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, “Cost function dependent barren plateaus in shallow parametrized quantum circuits,” Nature Communications, vol. 12, Mar. 2021.
  * [57] A. Arrasmith, M. Cerezo, P. Czarnik, L. Cincio, and P. J. Coles, “Effect of barren plateaus on gradient-free optimization,” Quantum, vol. 5, p. 558, 2021.
  * [58] Z. Holmes, A. Arrasmith, B. Yan, P. J. Coles, A. Albrecht, and A. T. Sornborger, “Barren plateaus preclude learning scramblers,” Phys. Rev. Letters, vol. 126, no. 19, p. 190501, 2021.
  * [59] M. Larocca, P. Czarnik, K. Sharma, G. Muraleedharan, P. J. Coles, and M. Cerezo, “Diagnosing barren plateaus with tools from quantum optimal control,” Quantum, vol. 6, p. 824, 2022.
  * [60] E. Fontana, D. Herman, S. Chakrabarti, N. Kumar, R. Yalovetzky, J. Heredge, S. H. Sureshbabu, and M. Pistoia, “The adjoint is all you need: Characterizing barren plateaus in quantum ansätze,” arXiv preprint arXiv:2309.07902, 2023.
  * [61] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, “Cost function dependent barren plateaus in shallow parametrized quantum circuits,” Nature communications, vol. 12, no. 1, p. 1791, 2021.
  * [62] T. L. Patti, K. Najafi, X. Gao, and S. F. Yelin, “Entanglement devised barren plateau mitigation,” Phys. Rev. Res., vol. 3, no. 3, p. 033090, 2021.
  * [63] C. O. Marrero, M. Kieferová, and N. Wiebe, “Entanglement-induced barren plateaus,” PRX Quantum, vol. 2, no. 4, p. 040316, 2021.
  * [64] Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles, “Connecting ansatz expressibility to gradient magnitudes and barren plateaus,” PRX Quantum, vol. 3, no. 1, p. 010313, 2022.
  * [65] M. Schumann, F. K. Wilhelm, and A. Ciani, “Emergence of noise-induced barren plateaus in arbitrary layered noise models,” 2023.
  * [66] S. H. Sack, R. A. Medina, A. A. Michailidis, R. Kueng, and M. Serbyn, “Avoiding barren plateaus using classical shadows,” PRX Quantum, vol. 3, p. 020365, Jun 2022.
  * [67] H.-Y. Huang, R. Kueng, and J. Preskill, “Predicting many properties of a quantum system from very few measurements,” Nature Physics, vol. 16, no. 10, pp. 1050–1057, 2020.
  * [68] K. Zhang, L. Liu, M.-H. Hsieh, and D. Tao, “Escaping from the barren plateau via gaussian initializations in deep variational quantum circuits,” Advances in Neural Information Processing Systems, vol. 35, pp. 18612–18627, 2022.
  * [69] L. Friedrich and J. Maziero, “Avoiding barren plateaus with classical deep neural networks,” Phys. Rev. A, vol. 106, no. 4, p. 042433, 2022.
  * [70] L. Binkowski, G. KoÃŸmann, T. J. Osborne, R. Schwonnek, and T. Ziegler, “From barren plateaus through fertile valleys: Conic extensions of parameterised quantum circuits,” 2023.
  * [71] E. Grant, L. Wossnig, M. Ostaszewski, and M. Benedetti, “An initialization strategy for addressing barren plateaus in parametrized quantum circuits,” Quantum, vol. 3, p. 214, 2019.
  * [72] T. Volkoff and P. J. Coles, “Large gradients via correlation in random parameterized quantum circuits,” Quantum Science and Technology, vol. 6, no. 2, p. 025008, 2021.
  * [73] G. Verdon, M. Broughton, J. R. McClean, K. J. Sung, R. Babbush, Z. Jiang, H. Neven, and M. Mohseni, “Learning to learn with quantum neural networks via classical neural networks,” arXiv preprint arXiv:1907.05415, 2019.
  * [74] A. Skolik, J. R. McClean, M. Mohseni, P. van der Smagt, and M. Leib, “Layerwise learning for quantum neural networks,” Quantum Machine Intelligence, vol. 3, pp. 1–11, 2021.
  * [75] N. Alon, C. Avin, M. Koucky, G. Kozma, Z. Lotker, and M. R. Tuttle, “Many random walks are faster than one,” in Proceedings of the twentieth annual symposium on parallelism in algorithms and architectures, pp. 119–128, 2008.
  * [76] W. Fleming and M. Viot, “Some measure-valued Markov processes in population genetics theory,” Indiana Univ. Math. J., vol. 28, no. 5, pp. 817–843, 1979.
  * [77] K. Burdzy, R. Holyst, and P. March, “A Fleming-Viot particle representation of the Dirichlet Laplacian,” Comm. Math. Phys., vol. 214, no. 3, pp. 679–703, 2000.
  * [78] I. Grigorescu and M. Kang, “Hydrodynamic limit for a fleming–viot type system,” Stochastic processes and their applications, vol. 110, no. 1, pp. 111–143, 2004.
  * [79] A. Budhiraja, N. Fraiman, and A. Waterbury, “Approximating quasi-stationary distributions with interacting reinforced random walks,” ESAIM: Probability and Statistics, vol. 26, pp. 69–125, 2022.
  * [80] X. Chen, S. S. Du, and X. T. Tong, “On stationary-point hitting time and ergodicity of stochastic gradient langevin dynamics,” Journal of Machine Learning Research, 2020.
  * [81] Y. Zhang, P. Liang, and M. Charikar, “A hitting time analysis of stochastic gradient langevin dynamics,” in Conference on Learning Theory, pp. 1980–2022, PMLR, 2017.
  * [82] Z. Ciesielski and S. Taylor, “First passage times and sojourn times for Brownian motion in space and the exact hausdorff measure of the sample path,” Transactions of the American Mathematical Society, vol. 103, no. 3, pp. 434–450, 1962.
  * [83] M. Pincus, “Letter to the editor—a monte carlo method for the approximate solution of certain types of constrained optimization problems,” Operations Research, vol. 18, pp. 1225–1228, Dec. 1970.
  * [84] A. Khachaturyan, S. Semenovsovskaya, and B. Vainshtein, “The thermodynamic approach to the structure analysis of crystals,” Acta Crystallographica Section A, vol. 37, pp. 742–754, Sept. 1981.
  * [85] P. J. M. van Laarhoven and E. H. L. Aarts, Simulated Annealing: Theory and Applications.  Springer Netherlands, 1987.
  * [86] 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, vol. 10, no. 2, p. 021067, 2020.