# The promising path of evolutionary optimization to avoid barren plateaus

Jakab Nádori,1 Gregory Morse,2 Zita Majnay-Takács,2 Zoltán Zimborás,2,3,4
Péter Rakyta,1,3∗

1Department of Physics of Complex Systems, Eötvös Loránd University, Budapest,
Hungary
2Department of Programming Languages and Compilers, Eötvös Loránd University,
Budapest, Hungary
3Wigner Research Center for Physics, P.O. Box 49, 1525 Budapest, Hungary
4Algorithmiq Ltd, Kanavakatu 3C 00160 Helsinki, Finland
∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
> Variational quantum algorithms are viewed as promising candidates for
> demonstrating quantum advantage on near-term devices. These approaches
> typically involve the training of parameterized quantum circuits through a
> classical optimization loop. However, they often encounter challenges
> attributed to the exponentially diminishing gradient components, known as
> the barren plateau (BP) problem. This work introduces a novel optimization
> method designed to alleviate the adverse effects of BPs during circuit
> training. Our approach to select the optimization search direction relies on
> the distant features of the cost-function landscape. This enables the
> optimization path to navigate around barren plateaus without the need for
> external control mechanisms. We have successfully applied our optimization
> strategy to quantum circuits comprising $16$ qubits and $15000$ entangling
> gates, demonstrating robust resistance against BPs. Additionally, we have
> extended our optimization strategy by incorporating an evolutionary
> selection framework, enhancing its ability to avoid local minima in the
> landscape. The modified algorithm has been successfully utilized in quantum
> gate synthesis applications, showcasing a significantly improved efficiency
> in generating highly compressed quantum circuits compared to traditional
> gradient-based optimization approaches.

## Introduction

The Variational Quantum Eigensolver (VQE) originally formulated by Peruzzo et
al. (?) and McCLean et al. (?), has garnered considerable attention within the
research community over the past two decades (?). VQE uses the variational
principle to calculate the ground state energy of a Hamiltonian, a fundamental
problem in quantum chemistry and condensed matter physics. Conventional
computing methods face limitations in accuracy and efficiency, due to the
exponentially growing complexity of the quantum state used to describe many-
body systems. VQE, on the other hand, is expected to be efficient to model
these complex wave-functions in polynomial time, making it one of the most
promising near-term applications for quantum computing. In addition,
variational quantum algorithms have demonstrated some degree of resilience to
noise in the quantum hardware, making them a good candidate to deliver
reasonable calculation results even on near-term quantum computers.
Consequently, VQE stands as a prominent example of NISQ algorithms. In its
most general form, it seeks to compute an upper bound for the ground-state
energy of a Hamiltonian, typically provide the first step to determine the
energetic properties of molecules and materials. The potential applications of
VQE are therefore extensive, encompassing areas such as drug discovery (?, ?),
material science (?), chemical engineering (?, ?, ?), quantum optimization (?,
?), and quantum machine learning (?, ?).

While the training process in variational quantum algorithms is proven to be
NP-Hard (?), the literature highlights four primary areas of research aimed at
bringing VQE closer to practical real-world applications. These areas
encompass: (i) The development of optimal measurement techniques to minimize
the required number of circuit repetitions for VQE execution. (ii) The
establishment of extensive parallelization methods across multiple quantum
computers. (iii) The creation of strategies to address potential challenges
associated with vanishing gradients during the optimization process for larger
systems. (iv) The exploration of tailored error mitigation approaches
specifically designed for the VQE algorithm. The solutions to these open
research questions will show the pathways for VQEs to attain quantum
advantage, particularly as quantum computing hardware scales up and noise
levels hopefully diminish.

In this study, we address a common pitfall encountered in variational quantum
algorithms during parameter optimization, known as the barren plateau (BP)
problem, as described in prior research (?, ?, ?). The BP arises due to the
phenomenon where the gradients of the cost function vanish exponentially as
the number of qubits in the quantum model increases. For cost functions that
exhibit BPs, exponentially many measurements may be required to determine the
minimization direction in gradient-based optimization. This challenge in
scaling applies not only to gradient-based optimization but also to
derivative-free optimization methods, as discussed in Arrasmith et al. (?),
and remains unresolved even with the application of higher-order optimization
techniques (?). Research efforts have linked the BP problem to the overall
expressiveness of the variational quantum circuit ansatz (?). Also, the degree
of entanglement of the wave-function (?, ?, ?, ?), and the non-locality of the
cost function (?, ?, ?) have also high influence on encountering BPs during
the optimization process. These findings imply that BPs are a fundamental
property of variational quantum algorithms at scale, needing extra measures to
mitigate their impact.

Numerous methods have been proposed to address this issue, but it remains
uncertain whether these approaches can resolve the BP problem. Early research
initiatives focused on strategies to circumvent BPs during the initialization
phase of algorithms. For instance, Grant et al. (?) suggested initializing the
circuit with blocks containing identity gates. Skolik et al. (?) proposed
optimizing the circuit layer by layer, while Dborin et al. (?) recommended
employing a matrix product state ansatz for the optimization to reduce
entanglement. An alternative strategy was formulated in the study of (?) to
take advantage of correlated or restricted single-qubit rotation parameters
(or direction). An alternative approach utilizing the method of natural
gradient (?, ?) was also investigated in the literature and has been found to
be a promising strategy compared to conventional gradient-based optimizers.
Though, the individual optimization steps turned to be quite expensive and
it’s tolerance against BPs was not addressed.

Taking a different approach, the relationship between the occurrence of BPs
and the structure of the cost function has been explored in references (?) and
(?). Additionally, the introduction of concepts like entanglement-induced BPs
(?) and noise-induced BPs (?) have added new routes to the discussion. The
relation between BPs and entanglement has lead to various proposals that
suggest controlling entanglement to mitigate BPs (?, ?, ?, ?). Furthermore,
the idea of weak BPs, as defined for a subset of qubits (?), has been employed
to detect early indicators of BPs and prevent them from occurring. In
addition, (?) has demonstrated that entanglement-induced BPs and BPs arising
from local cost functions are essentially the same. Therefore, strategies to
avoid entanglement-induced BPs are equivalent to strategies for avoiding BPs
associated with local cost functions. It’s important to highlight that all of
the suggested methods for mitigating BPs primarily focus on avoiding entering
onto a BP in the first place. There is currently no known methodology to
follow when the optimization process is already in the midst of a barren
plateau. One might even argue that barren plateaus can be utilized to
characterize the trainability of variational quantum algorithms: an algorithm
shows better trainability when experiencing less impact from barren plateaus.

In this study, we turn over this reasoning and investigate the possibility to
engineer a BP resistant optimization technique, showing lower degree of
dependence on BP indicators like the vanishing gradient components or high
entanglement entropy. We present a novel evolutionary strategy that has proven
to be remarkably effective in optimization, even when initiated from a point
in the parameter space associated with high entanglement entropy. This
optimization approach seamlessly integrates into the standard VQE framework,
as it only requires the evaluation of the target cost function to be
optimized. The idea to use evolutionary machine learning practices has been
proven in several scenarios, trying to increase the training capabilities of
variational quantum algorithms in the context of quantum architecture search
(?, ?, ?, ?, ?, ?, ?) Though, none of these studies delivered a viable
approach to mitigate the consequences of BPs on the optimization process. Our
novel evolutionary approach implemented in the SQUANDER circuit optimization
package (?) attempts to fill this gap. We demonstrate the efficiency of our
approach through various applications, including optimizing the ground state
of the model like the XXX Heisenberg or the SYC model (?, ?, ?), and quantum
gate compilation of deep circuits. Among the key advantages of our algorithm,
we can also mention the computations cost of a single iteration do not exceed
the complexity of a conventional gradient based approach optimizing over the
same number of parameters. By evaluating distant features of the cost-function
landscape, the optimization strategy can avoid flat areas and not even miss
narrow valleys leading to acceptable solutions. While our studies are focusing
on VQEs, they readily extend to other variational hybrid algorithms as well,
such as quantum machine learning, quantum optimization or variational time
evolution. The developed optimization algorithm can be also ported to
realistic experiments performed on quantum processors, since it relies on the
evaluation of the cost function with simple phase shifts applied on the
parameters. We performed our numerical experiments using quantum computer
simulations implemented in the SQUANDER package (?).

The structure of the paper is organized as follows. In the first section
following the introduction we overview the concept of the VQE algorithm.
Subsequently, we present the technical background of our novel optimization
method exploiting long-range properties of the optimization landscape to
determine the search direction and performing evolutionary selection strategy
to further increase the success rate to find the global minimum. In the third
section we outline the concept of the entanglement entropy – BP link, using
the second-Rényi entropy for BP inference. In the fourth section we present
our numerical results on training quantum to approximate the ground states of
the Heisenberg and SYC models. Finally, we conclude our work in the last
section.

## Brief description of the VQE algorithm

As implied by its name, VQE aims to approximate the ground state of a quantum
system that is challenging to simulate using classical hardware. Thus, the
first step in the VQE is to define the system for which we want to find the
ground state. By providing a Hamiltonian $\hat{H}$ and a trial wave-function
$|\Psi\rangle$, the ground state energy associated with this Hamiltonian,
$E_{0}$, bounded by

$E_{0}\leq\frac{\langle\Psi|\hat{H}|\Psi\rangle}{\langle\Psi|\Psi\rangle}$ (1)

The primary objective of VQE is to determine a parameterization for
$|\Psi\rangle$ that minimizes the expectation value of $\hat{H}$. This
expectation value serves as an upper bound for the ground state energy (?)
and, in an ideal scenario, should be practically identical to it, meeting the
desired level of precision. Therefore, VQE operates by variational
optimization of the parameters within the quantum circuit ansatz.

Figure 1: Schematic process of the VQE algorithm. The expectation value sof
the Hamiltonian is evaluated using the quantum processor parameterized with a
parameter vector $\boldsymbol{\theta}$. The VQE energy is optimized via
classical computer by updating the parameter vector $\boldsymbol{\theta}$
provided for the quantum circuit in each iteration.

To transform the minimization task into a quantum computer-executable problem,
the first step is to define an ansatz quantum circuit used to prepare the
trial wave-function on a quantum processor. Therefore, we represent
$|\Psi\rangle$ as the result of applying a generic parameterized unitary
operation $U(\boldsymbol{\theta})$ to an initial state of $N$ qubits, with the
parameter set denoted as $\boldsymbol{\theta}$, taking values in the range
$\left(0,2\pi\right]$. For simplicity, the qubit register is typically
initialized as $|0\rangle^{\otimes N}$, denoted as $|\mathbf{0}\rangle$. The
VQE procedure is illustrated in Fig. 1. Given that $|\Psi\rangle$ must be a
normalized wave-function due to the unitary evolution, we can now formulate
the VQE optimization problem as follows, omitting the normalization factor
from Eq. (1):

$E_{VQE}=\min\limits_{\theta}\langle\mathbf{0}|U^{\dagger}(\theta)\hat{H}|U(\theta)\mathbf{0}\rangle$
(2)

Equation (2) is referred to as the cost function in the context of the VQE
optimization problem. This terminology is adopted from the fields of machine
learning and mathematical optimization literature.

To implement the VQE algorithm on real quantum processors, a typical approach
is to expand the Hamiltonian in terms of quantities that can be directly
measured on a quantum device. This representation is often achieved by
expressing the Hamiltonian as a weighted sum of Pauli spin operators.
Observables suited for immediate measurement on a quantum device can be formed
by tensor products of spin operators. These observables can be defined as
Pauli strings denoted as $\hat{P}_{\alpha}\in{I,X,Y,Z}^{\otimes N}$, where $N$
represents the number of qubits utilized in the quantum register. The
Hamiltonian can then be reformulated as follows:

$\hat{H}=\sum\limits_{\alpha}^{\mathcal{P}}{w_{\alpha}\hat{P}_{\alpha}}$ (3)

with $w_{\alpha}$ being a set of weights, and $\mathcal{P}$ the number of
Pauli strings in the Hamiltonian. Equation (2) then becomes

$E_{VQE}=\min\limits_{\theta}\sum\limits_{\alpha}^{\mathcal{P}}{w_{\alpha}\langle\mathbf{0}|U^{\dagger}(\theta)\hat{P}_{\alpha}|U(\theta)\mathbf{0}\rangle}$
(4)

where the terms
$E_{P_{a}}=\langle\mathbf{0}|U^{\dagger}(\theta)\hat{P}_{\alpha}|U(\theta)\mathbf{0}\rangle$
corresponds to the expectation value of a Pauli string $\hat{P}_{\alpha}$ and
can be measured by a quantum device, while the summation of these terms are
calculated using a classical machine. Generally, to evaluate $E_{VQE}$ with
the aid of a quantum chip it is required to sample the expectation values of
the Pauli strings several times. Typically, in order to attain a precision of
$\epsilon$ for the expectation value of an operator, it is necessary to
conduct approximately $\mathcal{O}(1/\epsilon^{2})$ iterations of the circuit
execution, each terminated by a measurement at the end [38]. Therefore, to
increase the efficiency of the VQE execution measurement strategies are
studied to minimize the number of the required repetitions [82,120,121].

Once the parameters $\mathbf{\theta}$ have been successfully optimized, the
trial state $|\Psi\rangle$ serves as a model for the system’s ground state
wave-function under investigation. Throughout the optimization process, the
ansatz parameters must undergo iterative updates until convergence is
achieved. Therefore, the selection of a classical optimization algorithm plays
a crucial role. Firstly, it directly influences the number of measurements
needed for each optimization step and the total number of iterations required
to reach convergence. Secondly, specific optimization strategies have the
potential to partially mitigate particular optimization challenges, such as
the BPs problem [131–135]. Though, a full-fledged solution to handle BPs
during the optimization process still forms an open challenge for researchers,
as the proposed strategies rather aim to avoid BPs and not overcome them.

## Evolutionary optimization technique to train quantum circuits

Evolutionary algorithms (?, ?) apply the principles of evolution found in
nature to find the most optimal configuration of properties to accomplish
specific tasks. They are often applied in scenarios where solving problem in
polynomial time is not feasible, like in the case of combinatorial
optimization problems, leading to NP-hard tasks at scale. In this sense,
evolutionary algorithms optimize among different possibilities in the space of
solution candidates, which are encoded in specific properties of the
individuals, referred to as _agents_ in the followings. Evolution is therefore
based on recombination over a group of agents, accompanied with random
mutations. Natural selection between the agents is performed according to the
optimization criteria, encoded in the cost function. This process is repeated
until the individuals can not improve any more in solving the optimization
problem.

The agents explore the parameter space to find a favorable route that
eventually decreases the cost function the most. A possible and numerically
efficient scenario might be to follow a trajectory traced by the optimization
of the cost function with respect to a small subset of the free parameters. In
other words, the agents perform in each iteration an exact line search along
the direction determined in a randomly selected parameter-subspace while
keeping the rest of the parameters constant. Studies of Refs. (?) and (?)
showed that by meeting certain conditions in the parameterization of the
quantum gates, the single parameter dependence of the energy landscape in a
VQA algorithm can be described by a simple sinusoidal function with a period
of $2\pi$. This finding has inspired the creation of two high performance
optimization algorithm demonstrated in this work. First of all we should
mention that merely three calls to the cost function are sufficient to find
the exact minimum with respect to a single parameter. Secondly, global-scaled
properties on the optimization landscape can be collected by sampling the cost
function along single-parameter directions, leading to the formulation of a
novel strategy on deciding the search direction during the optimization
iterations. These mathematical findings enable a monotone step-by-step
reduction of the cost function by iteratively selecting a new subset of
parameters during the optimization process. The result of the procedure is
determined by both the starting point in the parameter space and the sequence
of the chosen parameters. Unfortunately, neither the optimal starting point
nor the most optimal parameter subset sequence is known in advance. While a
predefined sequence of parameters, such as an iterative parameter sweep
following a periodic pattern, might prematurely lead to convergence in most
cases, a well-chosen sequence of optimized parameters could significantly
enhance the overall optimization efficiency. In this work we show via
numerical experiments, that probability based approaches promise a high
success rate in finding a sub-optimal sequence of the optimized parameters.

### Evolutionary selections

The algorithm described in this study is based on a highly-efficient
parameter-wise optimization routine and on an evolutionary selection of the
sequence of the parameters over which the cost function is minimized. In our
algorithm, each agent randomly selects the next parameters to be optimized by
sampling from a uniform distribution. After optimizing the cost function with
respect to the chosen parameters, agents pick up a next subset to optimize,
again via uniform distribution sampling.

This repetitive process charts a path in the parameter space, over which the
cost function is being optimized. Upon completing a specified number of
iterations (set to a value between $200$ and $20000$ in our reference
implementation within the SQUANDER package), the evolution of the agents are
compared to each other and the most successful agent (returning the most
favorable cost function) is selected. We refer to this step in the algorithm
as _agent synchronization protocol_.

At this stage, the individual agents are randomly assigned to either continue
their optimization path with a probability of $p_{exploration}$ or start a new
route from the actual position of the most successful agent. The latter option
has a probability of $1-p_{exploration}$. In the reference implementation of
the algorithm $p_{exploration}=0.2$ was chosen as default.

Figure 2: Illustration of the evolutionary algorithm in a corner case of a
single-parameter optimization per iteration per step. In each episode the
agents perform parameter-wise optimization along a path defined by a randomly
chosen sequence of parameters. At the end of an epoch the most successful
agent is selected and its state is spawned over the remaining agents for the
next episode. The remaining agents are either terminated (indicated with red
crosses), or randomly selected to keep their state and resume their evolution
to further explore the parameter space (see the red arrows indicating the
chosen agents). Some randomly chosen to implement a randomization on their
stat e(see the blue arrows indicating the chosen agents). The figure
illustrate a specific example with $11$ agents, resulting in a sequence of
best solutions $f_{2}>f_{11}>f_{8}>\dots>f_{1}$, with the final result $f_{1}$
after $L$ episodes.

We split the optimization algorithm into episodes: each episode terminates
with the described agent synchronization protocol and encompasses the
evolution of the agents between two synchronization points. To further enhance
the success rate through avoiding local minima of the energy landscape, the
parameter vector associated with a given agent undergoes a perturbation with a
probability $p_{randomization}$ right after the agent-synchronization
protocol. This action causes an adjustment in the elements of the parameter
vector by a randomly selected value from the range $[-r,r]$. In our reference
implementation of the algorithm $p_{randomization}=0.2$, and $r=0.3$ was
chosen as default. The described algorithm is schematically visualized in Fig.
2.

### Single parameter optimization

The algorithm’s performance is highly influenced by the efficiency of the line
search that is commonly enhanced by performing inexact line search (?).
Instead of finding the exact minimum in the search directions, an inexact line
search terminates after finding an acceptably good point providing sufficient
decrease in the cost function. In our case, however, the sinusoidal structure
of the landscape enables to design an efficient method to find the exact
minimum with respect to the chosen parameter just in few evaluation of the
cost function. Firstly, the circuit can be decomposed into a sequence of
single-parameter unitaries, each having a form
$U_{i}=\textrm{exp}(\theta_{i}\hat{H}_{i})$, where $\hat{H}_{i}$ is a
hermitian matrix generating a rotation. Scaling this matrix allows each
unitary $U_{i}$ to be adjusted to have a period of $2\pi$. Secondly, if
$\hat{H}_{i}$s are full-rank matrices, the unitaries $U_{i}$s can’t facilitate
controlled rotations. Consequently, parameter-free entangling gates are added
to the circuit which do not affect the sinusoidal nature of the single-
parameter dependence of the cost function which reads as

$E_{VQE}=\kappa\cdot\sin(2p_{i}+\xi)+C\;.$ (5)

The derivation of this equation is straightforward as the state generated by
the circuit can be given in the form

$|\Psi\rangle=\cos(p_{i})|a\rangle+\sin(p_{i})|b\rangle$ (6)

with respect to the parameter $p_{i}$, where $|a\rangle$ and $|b\rangle$ are
normalized quantum states. By evaluating the expectation value
$\langle\Psi|\hat{H}|\Psi\rangle$ we end up with Eq. (5) by introducing
suitable constants $\kappa$, $\xi$ and $C$. In this work we generalize this
equation to cases when parametric controlled rotation gates are also
incorporated in the circuit, particularly significant in certain approaches of
quantum gate compilations (?, ?). In this scenario the general form of the
cost function combines two frequencies:

$E_{VQE}=\kappa\cdot\sin(2p_{i}+\xi)+\gamma\cdot\sin(p_{i}+\phi)+C\;.$ (7)

The extra sinusoidal term in this expression can be explained by the
additional constant term $|c\rangle$ in the quantum state originating from
controlled unitary gate operations. Namely, the state generated by the circuit
as a function of the parameter $p_{i}$ can be written in the form

$|\Psi\rangle=\cos(p_{i})|a\rangle+\sin(p_{i})|b\rangle+|c\rangle$ (8)

By evaluating the expectation value of the Hamiltonian, Eq. (7) can be
recovered. Unlike the original setup with $3$ cost function calls, in this
case $5$ samples are needed to fully recover the cost function and to solve
the minimization problem with respect to a single parameter. For more details
see the Sec. S2 of the supplementary material. Computationally, the increased
number of the required function evaluations might question the application of
this generalized approach, however, the enhanced expressibility of the circuit
with parametric controlled rotations might improve the trainability of the
circuit, which has been already proven in quantum gate compilation experiments
(?, ?). In VQE experiments controlled rotation gates can be considered as an
abstraction, encompassing two parameter-free controlled two-qubit gates and
two single-qubit rotations with correlated angles of rotation. As was shown in
the study of Refs. (?) and (?), correlated parameters in the circuits might
increase their trainability.

### Optimization over gradient-free multi-parameter search direction

To further generalize the optimization approach and potentially improve its
efficiency, we extend the concept of the parameter-wise optimization to a
finite subset of the parameters. While this extension adds complexity to the
line search phase of the optimization, it has the potential to result in much
faster convergence of the solution.

As mentioned in the preceding section, it requires $3$ function calls to
determine the exact minimum of the cost function along a single parameter
while keeping the others fixed. However, if the function value at the initial
point $\boldsymbol{\theta}$ is known from the previous optimization iteration,
the cost of determining the minimum with respect to a single parameter reduces
to $2$ function calls. Let’s assume that at the start of the iteration step, a
subset $\Lambda\subset{1,\dots,L}$ of the parameters was randomly selected
(where $L$ is the total number of parameters). For each parameter $\theta_{i}$
with $i\in\Lambda$, the parameter-wise minimum $\theta_{i}^{*}$ was determined
using Eq. (5). Then we define the search direction $\mathbf{d}$ in the
parameter space with components

$d_{i}=\left\\{\begin{matrix}\theta_{i}^{*}-\theta_{i}&\textrm{if
}i\in\Lambda\\\ 0&\textrm{otherwise.}\end{matrix}\right.$ (9)

The method for determining the gradient-free search direction is illustrated
schematically in Fig. 3. In each iteration, the search direction $\mathbf{d}$
is established with respect to the randomly selected subset of parameters.
Subsequently, a line search is conducted over the points of
$\boldsymbol{\theta}+t\cdot\mathbf{d}$ with $0\leq t\leq 1$ to identify the
point with the lowest cost function, which serves as the starting point for
the next iteration.

Following the best practices applied in conventional gradient-based
optimization methods, we incorporate an (gradient-free) inexact line search
into the algorithm by taking a predefined number of equidistant samples from
the cost function and selecting the point for the next iteration with the
lowest value.

Figure 3: A typical $E_{VQE}$ energy landscape as a function of two selected
parameters ”$i$” and ”$j$”, while the other parameters are kept fixed. The
blue colored arrows $d_{1}$ and $d_{2}$ label the two components of the search
vector given by Eq. (9). The search vector determines the path of the line
search between the initial (red colored) and minimal point (black colored). In
the next iteration step a new line search is performed originating from the
black point and optimizing over a new subset of parameters.

We assert that the described procedure will converge to either a global or
local minimum, and this is ensured by two key properties of the algorithm. (i)
Firstly, the cost function consistently decreases with each iteration. (ii)
Secondly, any local minimum in the landscape is an attracting fixed-point of
the optimization process. Therefore, like other optimization strategies, the
described algorithm might get trapped in a local minimum as well. An avenue
for significant enhancement could involve executing the optimization with
multiple agents and implementing evolutionary selection among them, as
detailed in a previous section. However, our numerical experiments have
indicated that even with a single agent, the described algorithm proves to be
quite efficient in training quantum circuits up to sizes accessible on
classical simulators, as getting trapped in local minima showed to be highly
unlikely.

## Discerning BPs via monitoring the entropy of subsystems

As mentioned in the introduction, a high vale of entanglement (or von Neumann)
entropy (?) is often accompanied with the occurrence of BPs in the
optimization landscape. The rapid growing of the entanglement entropy server
as an indicator of BP formation, making it a valuable monitoring tool during
optimization. Inspired by these findings, we also use the concept of entropy
monitoring during the optimization process to characterize the path in the
parameter landscape passed during the optimization process.

Moreover, the study of ref. (?) introduces the concept of weak BPs (WBPs) as
precursors of BPs. WBPs emerge when the entanglement of a local subsystem
exceeds a certain threshold identified by the entanglement of a fully
scrambled state. In contrast to BPs, WBPs can be efficiently diagnosed using
the few-body density matrices via evaluating the second Rényi entropy, given
by

$S_{2}=-\ln\textrm{Tr}\rho_{A}^{2}\;.$ (10)

Here $\rho_{A}=\textrm{tr}_{B}\rho$ is the reduced density matrix where $A$
denotes the subset of qubits that are measured and $B$ is the rest of the
system. While in a simulation the evaluation of the entropy is
straightforward, it’s value can also be inferred experimentally. The algorithm
introduced in (?), called classical shadow estimation, efficiently estimates
the expectation value not only of the cost function and its gradients but also
the second Rényi entropy of small-sized subsystems. For practical reasons, we
will always normalize the second Rényi entropy by the Page entropy $S_{Page}$,
which is the upper bound of the von Neumann entropy of a Haar random state
(?):

$S_{Page}(k,N)\approx k\ln 2-\frac{1}{2^{N-2k+1}}\;.$ (11)

Here $k$ labels the number of qubits in the chosen subspace in the $N$-qubit
register. Study of Ref. (?) showed that in the limit of (nearly) maximal
possible entanglement of a small sub-region A, the expected value of the
second Rényi entropy comes close to the Page entropy.

## Results

In this section, we present the outcomes of our numerical simulations
conducted on variational problems, specifically the estimation of the ground
state of an XXX Heisenberg Hamiltonian and the SYC model. In both of these
models, the ground state satisfies volume law, indicating that the
entanglement entropy of any bipartition of the state scales with the volume,
i.e., $(\rho_{A})\sim|A|$ (refer to Ref. (?) for a comprehensive review of
these concepts). Additionally, we share numerical results obtained from
quantum gate decomposition problem scenarios. Our findings demonstrate that
the developed evolutionary algorithm excels in deep quantum circuit synthesis
tasks, outperforming traditional optimization methods such as the well-known
Broyden–Fletcher–Goldfarb–Shanno (BFGS) (?) or the ADAM (?) optimization
methods. These traditional methods prove ineffective in producing meaningful
results for the gate decomposition.

We conducted the training of the targeted quantum circuits utilizing the high-
performance quantum computer simulator embedded in the SQUANDER package. To
assess the capabilities of the new optimization technique, we deliberately
excluded any noise in the simulation of the circuits and assumed perfect
measurements on the qubits. This approach allowed us to concentrate
specifically on the performance of the evolutionary optimizer across various
applications and leave the impact of noise and other imperfections of a VQE
experiment for future studies.

In the training of quantum circuits, we considered the circuit of the form of
”hardware-efficient” ansatz described by a unitary

$U\left(\boldsymbol{\theta}\right)=\prod\limits_{i=1}^{P}\left(W_{i}\prod\limits_{j=1}^{N}U_{i}^{j}(\theta_{i}^{j})\right)\;.$
(12)

The $N$-qubit circuit is composed from periodic layers of parameter-free
entangling layers $W_{i}$ and parametric single-qubit rotation layers composed
from independent single-qubit rotations $U_{i}^{j}$ acting on qubit $j$ in the
layer $i$. The entangling layers are composed from controlled-not (CNOT) gates
connecting nearest neighbor qubits, as shown in Fig.4 for the case of $5$
qubits.

Figure 4: The structure of the ”hardware-efficient” circuit ansatz used in our
numerical experiments. The single qubit rotations are expressed by three
consecutive rotations around axis $Z$, $Y$ and $Z$ again. The circuit ansatz
is built from identical layers, denoted with the parameter $p$.

Our initial objective in these numerical experiments is to demonstrate the
efficient utilization of the gradient-free multi-parameter search direction in
circumventing BPs during circuit training. To achieve this, we focus on a
specific scenario in which only a single agent operates on a randomly chosen
subset of the parameters. In this form, the strategy closely resembles
traditional gradient descend optimization algorithms. The key difference is
the incorporated routine to determine the search direction. While conventional
gradient-based optimizers choose the search direction based on local
properties of the cost-function landscape, such as the gradient or curvature,
our proposed alternative approach suggests that determining the search
direction via long-range properties of the optimization landscape yields
superior results.

Our second objective was to enhance the optimization strategy by introducing
mechanisms to prevent or significantly reduce the likelihood of getting
trapped in local minima within the optimization landscape. To achieve this, we
employed a multi-agent approach and implemented evolutionary selection among
them.

We showcase the effectiveness of the evolutionary approach on quantum circuit
synthesis problems, particularly in scenarios where our previous findings
indicated the presence of local minima leading to the termination of most
gradient-based optimizer processes. The gate synthesis problem emerged as a
promising candidate to investigate the impact of evolutionary selection on the
optimization’s success rate while maintaining a low level of complexity in the
cost-function evaluation. This aspect of our numerical experiments is crucial,
as optimization with evolutionary selection involves a significant amount of
redundant calculations that need to be discarded.

### Heisenberg XXX Model

The Heisenberg model, composed in terms of the $\hat{X}$, $\hat{Y}$ and
$\hat{Z}$ Pauli spin operators, encompasses the most general form of
Hamiltonian that can be simulated on a qubit-based quantum processor. When
scaled up, the Heisenberg model with randomly chosen coupled qubit pairs,
becomes hard to simulate by classical means. The reason for this is attributed
to the non-local interactions between the sites resulting in a volume-law
entanglement scaling for a typical bipartite cut. Specifically, here we focus
our study on $k$-local Hamiltonians $\hat{H}$, defined as sum of terms each
containing at most $k$ Pauli matrices. In our numerical simulations we take
$k$ to be finite and fixed, while the number of qubits $N\gg k$. A typical
example describing such a model for $k=2$ and including a bias term related to
the effect of an external magnetic field is given by equation

$\hat{H}_{XXX}=\sum\limits_{i,j\in
V_{\mathcal{G}}}J\left(\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{j}+\hat{\sigma}^{y}_{i}\hat{\sigma}^{y}_{j}+\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}\right)+h_{z}\sum\limits_{N}\hat{\sigma}^{z}_{i}\;.$
(13)

Despite of the local $k$-terms in the formula, according to the reasoning of
Refs. (?, ?) the VQE formulation of such models will still exhibit BPs in a
limit of $\mathcal{O}(poly(N))$ deep quantum circuits. Therefore, with a deep
enough quantum circuit the VQE algorithm will also suffer from a BP, making
this model to be a suitable test case to show the efficiency of the developed
evolutionary optimization technique in training the circuit.

Figure 5: VQE simulations implementing the strategy based on the gradient-free
multi-parameter search direction with a single agent. The training of the
quantum circuits was conducted to find the ground state for the $N=10$-qubit
XXX Heisenberg model on a $3$-regular random graph depicted in the subset of
sub-figure (a) and for a $10$-qubit SYC model. Panels (a) and (b) show the
evolution of the VQE energy and the normalized second-Rényi entropy evaluated
on the sites $0$ and $1$ (labeled by blue in the insets) for the Heisenberg
model. Sub-figures (c) and (d) show the same kind of results for the SYC
model. The energy and the normalized entropy of the exact ground state for the
Heisenberg system are $E_{0}=-21.28$ eV and $S_{2}/S_{Page}=0.667$. The exact
ground state of the SYC model is $E_{0}=-3.82$ eV with normalized entropy
$S_{2}/S_{Page}=0.980$. For the initialization we used random parameters.

### Sachdev–Ye–Kitaev model

From theoretical point of view, the Sachdev–Ye–Kitaev (SYK) model- (?, ?, ?),
with it’s non-local Hamiltonian structure in the Puali basis, is a frequently
visited representative of variational problems that are proven to be hard to
solve in VQE context. This is because the SYC model features a ground state
that is nearly maximally entangled (?, ?), so crossing a BP becomes inevitable
during the VQE process (?). From solid state physics point of view, the SYK
model can be used to describe $2N$ spinless Majorana fermions $\chi_{i}$
satisfying the following anticommutation relations
$\\{\chi_{i},\chi_{j}\\}=\delta_{i,j}$. The Hamiltonian of the model reads as

$\hat{H}_{SYC}=\sum\limits_{1\leq i<j<k<l\leq
2N}J_{i,j,k,l}\;\chi_{i}\chi_{j}\chi_{k}\chi_{l}\;,$ (14)

where the couplings $J_{i,j,k,l}$ are taken randomly from a Gaussian
distribution with zero mean and a variance of

$\textrm{var}\left(J_{i,j,k,l}\right)=\frac{3!}{(N-3)(N-2)(N-1)}J^{2}\;.$ (15)

This fermionic Hamiltonian can be adopted to a quantum hardware by mapping the
couplings between the Majorana fermions onto spin-chain variables by the
nonlocal Jordan-Wigner transformation:

$\chi_{2i}=\sigma^{x}_{1}\sigma^{x}_{2}\dots\sigma^{x}_{i-1}\sigma^{y}_{i}\qquad\chi_{2i-1}=\sigma^{x}_{1}\sigma^{x}_{2}\dots\sigma^{x}_{i-1}\sigma^{z}_{i}\;.$
(16)

such that $\\{\chi_{i},\chi_{j}\\}=\delta_{i,j}$ will hold on. This encoding
enables one to model $2N$ Majorana fermions with only $N$ logical qubits.
Generally speaking, we can consider the SYC model as a corner case of the
Heisenberg model, when no restriction is made on the number of coupled qubits
in a single term in the Hamiltonian. For our studies, we set $J=1$ and
consider a system of $N=10$ qubits.

### Numerical results

In Fig. 5, we present the results of our numerical simulations conducted on
$10$-qubit Heisenberg XXX and SYC models. These systems have been recognized
as challenging for VQE in numerous studies, and our tests with conventional
gradient-based optimizers affirmed these observations. To address the
optimization problem (2), we explored four strategies available within the
SQUANDER package. The Powell-variant of the BFGS optimizer (?) and the
Gradient Descent methods incorporate an inexact line search (?) in each
iteration, proving to be quite efficient for smaller problems (e.g., in the
synthesis of $3-5$ qubit circuits). The handcrafted implementations of these
optimizers also consider the periodicity of the optimization landscape by
appropriately bounding the line search during optimization. Additionally, we
included the ADAM optimization algorithm in our tests, although we did not
implement the adaptive learning rate strategy proposed in (?). Our findings
indicate that, despite exhaustive optimization and adaptation to quantum
circuit training circumstances, these optimization engines struggle to solve
these optimization problems at scale.

In contrary, our developed optimization algorithm demonstrated the ability to
mitigate the effects of BPs in situations where other optimization methods
fail. Notably, our approach avoids BPs without the need for external
interventions or monitoring of specific quantities. As anticipated, in our
numerical simulations, gradient-based optimization approaches became trapped
near the energy level of the initial point.

The left side of Fig.5 illustrates our VQE experiments on a $10$-qubit
Heisenberg model identical to the system studied in (?). In each iteration
step, we performed optimization by randomly selecting $64$ parameters of the
circuit and conducted a line search along the direction determined by Eq.(9).
Throughout the iterations, we utilized a single agent, and no evolutionary
steps were applied during the simulations. Still, the optimization problem
could be solved, while significantly reducing the computational complexity.

The outcome of this optimization process displayed a steady decrease in the
VQE energy during the iterations [see Fig.5.(a)], while the evolution of the
second-Rényi entropy demonstrated a rather non-trivial pattern [see
Fig.5.(b)]. As the VQE instances in our simulations were initiated to random
parameters, the entropy at the beginning of the optimization picked up large
values implying large entanglement in the resulting state vector. However, as
the optimization progressed, the entropy rapidly decreased, even falling below
the value calculated from the exact ground state. Given that a sharp decrease
in the evolution of the cost function generally correlates with low values of
the second-Rényi entropy, this observation is expected. As the search
direction and the bounds of the line search are determined via globally-scaled
properties of the optimization landscape, the optimization process can
navigate away from flat areas and avoid overstepping narrow valleys, guiding
toward acceptable solutions.

Upon further analysis of the numerical results, it is valuable to explore the
resulting properties of the quantum state $|\Psi(\boldsymbol{\theta})\rangle$
obtained via circuits of varying depths. In addition to examining the energy
and entropy associated with the resulting state, a direct comparison between
the state $|\Psi(\boldsymbol{\theta})\rangle$ and the exact ground state can
be characterized by their overlap integral

$M=\left|\langle\Psi_{0}|\Psi(\boldsymbol{\theta})\rangle\right|^{2}$ (17)

Values of $M$ close to $1$ signify that the VQE wave function faithfully
reproduces the exact ground state. On the left side of Fig.5, we compare a VQE
optimization performed with circuits of depths $p=50$ and $100$ layers (refer
to Fig.4 for the circuit structure). As depicted, with $100$ layers, both the
VQE energy and the second-Rényi entropy closely approach the theoretical
results of the exact ground state. The overlap integral $M$ calculated for the
produced VQE state is also very close to $1$, indicating that the VQE
experiment successfully reproduced the exact ground state with high accuracy.
(The resulting VQE circuit with parameters can be accessed within the SQUANDER
package (?)).

For a circuit with $50$ layers, on the other hand, we could not approximate
the ground state with such high quality. Neither the energy nor the entropy
[plotted in blue in Fig. 5.(a)-(c)] fully converged to the values
corresponding to the exact ground state. The overlap integral at the end of
the optimization (terminated at $\sim 3\times 10^{6}$ cost function calls) is
$M=0.97$. The results obtained for even shallower circuits rapidly become
inaccurate. The explanation behind this observation is related to the
expressiveness of the circuits. Considering $10$ qubits, a general wave-
function can be described by $2^{10}=1024$ complex amplitudes, resulting in
$2\times 2^{10}-2=2046$ free parameters (two parameters are constrained since
the state is normalized). The $50$-layered VQE circuit has $2700$ free
parameters, while the $100$-layered variant can be tuned with $5400$ free
parameters. As one can see, the $50$-layered circuits contain enough
parameters to encode any quantum state, but the convergence slows down at the
end of the optimization. The $100$-layered circuit converges much faster to
the solution, as it is likely over-parameterized.

Our results obtained for the SYC model further validate these considerations.
As depicted in the right panels of Fig. 5, the VQE optimization for the SYC
model exhibits similar behavior to the Heisenberg model on the left. With
$p=100$ layers in the circuit ansatz, we achieved a decent approximation of
the ground state energy. In contrast, at $p=50$ layers, the approximation
becomes less accurate after $6.5\times 10^{-6}$ evaluations of the cost
function. The second-Rényi entropy remains high throughout the optimization
process, though it does not hinder the optimizer from obtaining high-quality
results.

Regarding the exact solution of the SYC model, we need to mention that the
diagonalization of the SYC Hamiltonian revealed a $2$-fold degeneracy of the
states. In addition, the $2$-fold ground state and the nearest excited states
(also $2$-fold degenerate) were close to each other in energy. While the
energy of the ground state was $E_{0,SYC}=-3.819$ eV, the first excited energy
was $E_{1,SYC}=-3.799$ eV. The VQE energy obtained with the training of the
$p=100$ layered circuit turned out to be $E_{VQE}=3.814$ eV, a value that
falls between the ground state energy and the first excited energy. Thus, the
VQE state is expected to have a finite component orthogonal to the subspace of
the $2$-fold ground state. Fortunately, this orthogonal component has only an
overlap $M=0.098$ with the subspace of the first excited states, while the
overlap of the VQE wave function with the subspace of the $2$-fold ground
state is $M=0.902$. Hence, the components of the VQE wave function
predominantly come from the ground state and the first excited state.

At this point we would like to also mention that the optimization Ref. (?) was
not able to solve neither the same Heisenberg model nor the SYC model of
identical size with an accuracy comparable to our results. (We also notice
that (?) used slightly different circuit ansatz compared to ours.)

Figure 6: VQE simulations implementing the strategy based on the gradient-free
multi-parameter search direction with a single agent. The training of the
quantum circuits was conducted to find the ground state for the $N=16$-qubit
XXX Heisenberg model on a $3$-regular random graph depicted in the subset of
subfigure (c). Panels (a) and (b) show the evolution of the VQE energy and the
normalized second-Rényi entropy on sites $0$ and $1$ (labeled by blue in the
insets) for the Heisenberg model evaluated with shallower circuits consisting
of $p=50$ and $100$ layers, while . Subfigures (c) and (d) show the same kind
of results obtained with deep circuits consisting of $p=500$ and $1000$
layers. The black numbers and dots represent the evolution of the of the
overlap integral between the VQE wave function and the exact ground state
defined by Eq.(17). The energy and the normalized entropy of the exact ground
state for the Heisenberg system are $E_{0}=-32.84$ eV and
$S_{2}/S_{Page}=0.693$. The optimization was initialized with parameters equal
to $0$.

After examining $10$-qubit systems, we raise the stakes and turn our attention
to even more challenging problems. While it is not a significant challenge for
quantum computer simulators to execute circuits beyond $20$ qubits, the VQE
experiments become quite expensive at this scale due to the large number of
cost function evaluations. Therefore, in this work, we limit our
considerations to problems up to $16$ qubits. The goal of our VQE experiments
is to reproduce the ground state of a randomly generated $16$-qubit XXX
Heisenberg model (with $k=3$), as visualized in the inset of Fig. 6.(c).
Furthermore, to increase the efficiency of the optimization process, we
initiate the training of the circuits with parameter values equal to $0$.
Since the circuits are applied to the initial state $|\mathbf{0}\rangle$ and
initially, all the single-qubit rotations will be identities with this choice,
the two-qubit entangling gates will lose their effect as well. This is because
all the control qubits are initially in state $0$, so the controlled two-qubit
gates incorporated into the circuit will act as an identity gate as well.
Consequently, the resulting state will have zero initial entanglement entropy,
promising a steep initial descent in the VQE energy. Despite of this
simplification, our method remains autonomous in avoiding BPs on its own,
there is still no need for any external intervention during the optimization.

A general $16$-qubit wave function can be characterized by $2^{16}-2=131070$
parameters. Consequently, circuits with either $p=50$ or $p=100$ layers are
expected to provide much less accurate VQE wave functions than in the
previously addressed $10$-qubit case, where the number of trainable parameters
in the circuit was sufficiently large. In the $16$-qubit case, the
$100$-layered circuit ansatz contains only $9000$ free parameters, which is
much less than required to approximate the target state. The left panels of
Fig.6 show the numerical results of the conducted VQE experiments with
circuits consisting of $p=50$ and $p=100$ layers. As expected, it was not
possible to precisely approximate the ground state energy in these instances
(VQE energy of $-31.1$ eV with the $50$-layered circuit and $-31.6$ eV with
the $100$-layered circuit). The $100$-layered circuit gave significantly
better results than those achieved with the $50$-layered circuit. However,
neither of these circuits could get close to the exact ground state energy of
the model. After more than $10^{8}$ evaluations of the cost function, the VQE
wave function still showed a significant difference compared to the exact
ground state, and the overlap calculated with the wave function generated by
the $100$-layered circuit did not exceed the value of $M=0.75$. The black
numbers in Fig.6.(a) indicate the evolution of the overlap with the number of
cost function evaluations. As mentioned, these findings were not surprising,
as the number of trainable parameters in these circuits was only $4500$ and
$9000$ for the $50$\- and $100$-layered circuits, respectively. The
expressiveness of these circuits lags behind the numbers needed for the
accurate approximation of a general $16$-qubit wave function. Therefore, we
increased the number of layers in our circuit ansatz by an order.

The right panels of Fig.6 show our numerical results obtained with these deep
circuits composed of $p=500$ and $p=1000$ layers. As the plotted results
indicate, the developed optimization strategy had no difficulties proceeding
with the training of these deep circuits. We still executed only a single
agent in these numerical experiments, and optimized $128$ parameters in each
iteration. We observed that the efficiency of the optimization varies from run
to run. Fig.6 shows the most successful executions, while in several
experiments, we observed signs of convergence to a local minimum around the
energy level of $\sim-29$ eV. We argue that this issue could be resolved by
increasing the number of agents in the process due to the evolutionary
selection of the most successful agents. However, due to the overwhelming
numerical requirements of the performed VQE simulations, we did not try out
this approach in these simulations. We revisit this aspect of the algorithm in
the next section by addressing quantum gate synthesis problems, where it
becomes essential to avoid local minima to retrieve meaningful results. In
these experiments, evolutionary selection will play a crucial role in solving
the addressed problems.

In the conducted $16$-qubit VQE optimization our numerical simulation yielded
a VQE energy of $E_{VQE}=-31.99$ eV with $p=500$ layers and $E_{VQE}=-32.00$
eV with $p=1000$ layers after $8\times 10^{7}$ cost function evaluation, at
which point the converge towards the exact ground state energy has
dramatically slowed down. The largest value of the overlap integral (17) with
the exact ground state turned to be $M=0.86$.

Our numerical results highlight two main conclusions: (i) Firstly, the
probabilistic algorithm developed in this work, making use of the global
properties of the optimization landscape, is a promising candidate capable of
solving variational quantum problems at scale, even beyond the capabilities of
classical quantum computer simulators. (ii) Secondly, when using a general
circuit ansatz to perform VQE optimization (like the circuit structure used in
our experiments), highly accurate results require an exponentially scaling
number of layers in the circuit ansatz. We acknowledge that the circuit ansatz
used in our calculations may not be the most efficient one to model the ground
states of the studied models. Ansatz structures accounting for the symmetry of
the addressed physical model might significantly reduce the circuit depth for
practical applications (?, ?, ?). Our aim in this work was to examine the
numerical capabilities of the new solver strategy. The non-efficient
representation of the circuit ansatz did not impose limitations on us in
achieving this goal. On the contrary, we could execute quite deep circuits and
demonstrate that our strategy can be successfully applied even under extreme
circumstances.

### Quantum compilation of deep circuits

A commonly employed approach in the compilation and optimization of large
quantum circuits involves partitioning the circuit into smaller blocks that
can be individually optimized or re-synthesized. This allows for the treatment
of wide quantum circuits, as demonstrated in the work of Ref. (?), where
circuits ranging from $4$ to $400$ qubits were addressed. Previous approaches
following this partitioning strategy (?) were shown to scale to similar
circuit sizes but were practically limited to using only $3$-qubit sized
partitions due to inefficient numerical optimization. In the experiments of
Ref. (?), this limitation was overcome, extending the partition size to $6$
qubits. The choice of the number of qubits per partition is crucial, as more
qubits are likely to capture additional domain-level physical interactions,
leading to improved optimization opportunities.

On the contrary, for partitions containing more qubits, the optimization
problem becomes exceptionally challenging for traditional gradient-based
solvers due to the dense presence of local minima in the landscape (?, ?).
Here, we demonstrate that the evolutionary optimization algorithm described in
this study enables the application of the adaptive gate synthesis approach of
Ref. (?) to problems that were found to be intractable with conventional
optimization engines, such as the BFGS method (?). In order to measure the
distance between the initial unitary $U$ of the quantum program and the
synthesized unitary $V$ (obtained from the synthesized circuit) we used the
frobenius norm given by

$f(U,V)=\frac{1}{2}\left\|V-U\right\|_{F}^{2}=d-\textrm{Re}\left[{\textrm{Tr}}(U^{\dagger}V)\right],$
(18)

with $d$ being the dimension of the Hilbert-space. For further details of the
decomposing algorithm please refer to either Ref. (?) or to Ref. (?).

Circuit name | qubits | Initial | QISKIT | Ref (?, ?) | SQUANDER (?)
---|---|---|---|---|---
|  | $CNOT$ | $CNOT$ | $CNOT$ | $\overline{T}\,[s]$ | $CNOT$ | $\overline{T}\,[s]$
4gt5_76 | $5$ | $46$ | $529$ | $24$ | $1711$ | $22$ | $79$
4gt10-v1_81 | $5$ | $66$ | $372$ | $39$ | $65737$ | $37$ | $2245$
one-two-three-v1_99 | $5$ | $59$ | $302$ | $45$ | $80390$ | $35$ | $4693$
one-two-three-v0_98 | $5$ | $65$ | $213$ | $61$ | $175994$ | $42$ | $20313$
one-two-three-v2_100 | $5$ | $32$ | $502$ | $37$ | $5141$ | $27$ | $2687$
4mod7-v1_96 | $5$ | $72$ | $150$ | $33$ | $10255$ | $26$ | $2951$
aj_e11_165 | $5$ | $69$ | $337$ | $36$ | $15585$ | $23$ | $2823$
alu-v2_32 | $5$ | $72$ | $469$ | $41$ | $33820$ | $31$ | $1639$
alu-v4_36 | $5$ | $51$ | $193$ | $40$ | $11090$ | $28$ | $2106$
4gt5_77 | $5$ | $58$ | $338$ | $19$ | $2855$ | $17$ | $483$
4gt12-v0_86 | $6$ | $116$ | $1837$ | - | - | $40$ | $8174$
4gt12-v0_87 | $6$ | $112$ | $1837$ | $47$ | - | $36$ | $8598$
4gt12-v0_88 | $6$ | $86$ | $1837$ | $44$ | - | $36$ | $6577$
Table 1: $CNOT$ gate count comparison of deep $4-6$-qubit circuits obtained by
decomposing and optimizing unitaries taken from the online database (?). The
individual columns label the name of the circuit, the number of qubits in the
circuit, the initial $CNOT$ gate count of the circuit, the $CNOT$ gate count
of the quantum circuit synthesized by QISKIT transpile function. The last four
columns show the $CNOT$ count and average execution times obtained from Refs.
(?) and (?), and the same kind of results obtained with the new evolutionary
algorithm with the SQUANDER package. The final error $f$ of the approximation
calculated via Eq. (18) was less than $10^{-8}$ in each decomposition
corresponding. We performed the benchmark calculations on a computing server
equipped with $32$-Core AMD EPYC 7542 Processor (providing $64$ threads with
multi-threading) and with $128$GB of memory.

Our numerical experiences revealed that quantum gate synthesis provides a
perfect opportunity to explore problems where the presence of local minima in
the optimization landscape causes traditional optimizers to fail, even at
smaller scales. In our experiments, we observed that employing the developed
evolutionary algorithm with multiple agents is indeed necessary to synthesize
deep circuits while achieving a significant reduction in gate count. We
executed $64$ concurrent agents and synchronizing their evolution after
$200-2000$ optimization steps. To reduce the computational footprint, we
limited the randomly chosen set of parameters to contain only a single
parameter. Nevertheless, the developed evolutionary algorithm demonstrated
exceptional efficiency in solving the addressed problems. For optimal
performance, we combined the evolutionary algorithm with the BFGS strategy,
switching to the latter when the cost function was reduced down to $10^{-3}$
by the evolutionary algorithm. The results of our numerical experiments are
presented in Table. 1. Comparing these results with previous findings, we can
draw two conclusions. Firstly, the average execution time was reduced by an
order of magnitude in all addressed use cases where comparisons could be made.
Secondly, the evolutionary algorithm exhibited remarkable stability and
efficiency in solving the optimization problem with a precision of $10^{-8}$.

### Complexity considerations

It is worth to study the computational complexity of the developed
evolutionary algorithm described in this study. There are two key components
in relation with our optimization technique characterizing it’s cost. These
are the utilization of the agents undergoing evolutionary selections and the
size of the parameter subset over which the cost function gets optimized. Both
of these dimensions, though, fall back to the count of sampling from the
single parameter cross sections of the cost function. Since the terminating
point of the line search can be reused in the next iteration, only $2$ cost
function evaluations are needed to reconstruct the parameter dependence of the
cost function, if Eq. (5) is used to characterize the optimization landscape.
In the more general case of Eq. (7) $4$ cost function evaluations are required
to determine the minimum with respect to a single parameter. Therefore, $2LM$
(or $4LM$) cost function evaluations are required to determine the search
direction, and additional constant number of evaluations to perform the
(inexact) line search. These numbers compare close to the case of regular
gradient descend method. In particular, on a quantum computer the general rule
of parameter shift (?) can be used to evaluate a gradient component of the
cost function with $2$ evaluations of the cost function. By carefully planning
the number or agents and the number of optimized parameter components per
iteration, the overall complexity of our optimization algorithm can be set not
to exceed the complexity of conventional gradient based methods.

## Discussion

In this manuscript, we have presented a novel optimization method for training
quantum circuits. A prevalent challenge in variational quantum algorithms is
the hindering impact of barren plateaus (BPs) on circuit training. The
vanishing gradient problem associated with BPs significantly limits the
applicability of near-term quantum algorithms. Clearly speaking, the
application of a given variational algorithm becomes questionable when the
training process fails at large scaled problem instances. Therefore,
developing optimization strategies resistant to BPs is highly desirable. We
demonstrated that the long-range properties of the optimization landscape can
be leveraged to determine a search direction and optimization bounds along
this direction. This method introduces a novel line-search approach for
training quantum circuits, autonomously maintaining low values of the second-
Rényi entropy during optimization. Given the inherent connection between
entropy and BPs, this property of the developed algorithm implies a natural
tendency to avoid BPs during training without requiring external control
mechanisms.

We applied the optimization strategy to VQE simulations on randomly generated
Heisenberg XXX and Sachdev–Ye–Kitaev (SYC) models, all of which were
previously identified to exhibit barren plateaus (BPs) during VQE process. We
increased the depth of the training circuits to $p=1000$ entangling layers
with $16$ qubits. The developed optimization strategy, utilizing the gradient-
free search direction method, consistently decreased the VQE energy over
thousands of iterations without getting trapped by BPs. We demonstrated these
outstanding capabilities of the solver by training quantum circuits containing
up to $15000$ entangling gates.

In this study we also addressed another common issue of optimization
approaches: the local minima scattered over the optimization landscape causes
early termination of the optimization and prevents from solving the addressed
problem. This issue is particularly critical for quantum gate synthesis, where
only results close to the global minimum are acceptable. Our numerical
experiments demonstrated a significant increase in efficiency to avoid local
minima and converge to the global minimum when employing evolutionary
selection of agents. In most of the studied circuits, we achieved a
substantial reduction in gate count and decreased the computational time by an
order of magnitude. For these experiments, we utilized single-parameter
optimization to keep computational requirements low.

We believe that the optimal training performance for large-scale variational
quantum calculations can be achieved through the combination of these
practices. However, due to the exponential scaling of quantum circuit
simulation with the number of qubits, the limitations of classical computers
have prevented us from conducting exhaustive research by applying the
gradient-free multi-parameter search direction strategy and evolutionary
selection in combination.

Combining the described probabilistic optimization strategy with machine
learning practices has the potential to further improve its performance.
Replacing the probability-based selection with machine learning inference
could significantly reduce the number of redundant evaluations of the cost
function during circuit training. We leave the exploration of this possibility
for future research projects.

## Acknowledgments

Founding: This research was supported by the Ministry of Culture and
Innovation and the National Research, Development and Innovation Office within
the Quantum Information National Laboratory of Hungary (Grant No.
2022-2.1.1-NL-2022-00004), by the ÚNKP-22-5 New National Excellence Program of
the Ministry for Culture and Innovation from the source of the National
Research, Development and Innovation Fund, and by the Hungarian Scientific
Research Fund (OTKA) Grants No. K134437 and FK135220. RP. acknowledge support
from the Hungarian Academy of Sciences through the Bolyai János Stipendium
(BO/00571/22/11) as well. We acknowledge the computational resources provided
by the Wigner Scientific Computational Laboratory (WSCLAB) (the former Wigner
GPU Laboratory).

Data and materials availability: The compiled quantum circuits reported in
this work and the VQE trained circuits are publicly available within the
SQUANDER package in QASM and other binary formats at GitHub repository (?)

## References

  * 1. A. Peruzzo, et al., Nature Communications 5, 4213 (2014).
  * 2. J. R. McClean, J. Romero, R. Babbush, A. Aspuru-Guzik, New Journal of Physics 18, 023023 (2016).
  * 3. J. Tilly, et al., Physics Reports 986, 1 (2022). The Variational Quantum Eigensolver: a review of methods and best practices.
  * 4. Y. Cao, J. Romero, A. Aspuru-Guzik, IBM Journal of Research and Development 62, 6:1 (2018).
  * 5. N. S. Blunt, et al., Journal of Chemical Theory and Computation 18, 7001 (2022).
  * 6. V. Lordi, J. M. Nichol, MRS Bulletin 46, 589 (2021).
  * 7. A. Kandala, et al., Nature 549, 242 (2017).
  * 8. G. A. Quantum, et al., Science 369, 1084 (2020).
  * 9. Y. Cao, et al., Chemical Reviews 119, 10856 (2019).
  * 10. M. P. Harrigan, et al., Nature Physics 17, 332 (2021).
  * 11. N. Lacroix, et al., PRX Quantum 1, 020304 (2020).
  * 12. V. Havlíček, et al., Nature 567, 209 (2019).
  * 13. S. Johri, et al., npj Quantum Information 7, 122 (2021).
  * 14. L. Bittel, M. Kliesch, Phys. Rev. Lett. 127, 120502 (2021).
  * 15. D. Wecker, M. B. Hastings, M. Troyer, Phys. Rev. A 92, 042303 (2015).
  * 16. J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, H. Neven, Nature Communications 9, 4812 (2018).
  * 17. M. Cerezo, A. Sone, T. Volkoff, L. Cincio, P. J. Coles, Nature Communications 12, 1791 (2021).
  * 18. A. Arrasmith, M. Cerezo, P. Czarnik, L. Cincio, P. J. Coles, Quantum 5, 558 (2021).
  * 19. M. Cerezo, P. J. Coles, Quantum Science and Technology 6, 035006 (2021).
  * 20. A. A. Holmes, H. J. Changlani, C. J. Umrigar, Journal of Chemical Theory and Computation 12, 1561 (2016).
  * 21. C. Ortiz Marrero, M. Kieferová, N. Wiebe, PRX Quantum 2, 040316 (2021).
  * 22. T. L. Patti, K. Najafi, X. Gao, S. F. Yelin, Phys. Rev. Res. 3, 033090 (2021).
  * 23. Z. Holmes, et al., Phys. Rev. Lett. 126, 190501 (2021).
  * 24. K. Wang, Z. Song, X. Zhao, Z. Wang, X. Wang, npj Quantum Information 8, 52 (2022).
  * 25. A. V. Uvarov, J. D. Biamonte, Journal of Physics A: Mathematical and Theoretical 54, 245301 (2021).
  * 26. K. Sharma, M. Cerezo, L. Cincio, P. J. Coles, Phys. Rev. Lett. 128, 180505 (2022).
  * 27. E. Grant, L. Wossnig, M. Ostaszewski, M. Benedetti, Quantum 3, 214 (2019).
  * 28. A. Skolik, J. R. McClean, M. Mohseni, P. van der Smagt, M. Leib, Quantum Machine Intelligence 3, 5 (2021).
  * 29. J. Dborin, F. Barratt, V. Wimalaweera, L. Wright, A. G. Green, Quantum Science and Technology 7, 035014 (2022).
  * 30. Z. Holmes, K. Sharma, M. Cerezo, P. J. Coles, PRX Quantum 3, 010313 (2022).
  * 31. T. Haug, M. S. Kim, Optimal training of variational quantum algorithms without barren plateaus (2021).
  * 32. D. Wierichs, C. Gogolin, M. Kastoryano, Phys. Rev. Res. 2, 043246 (2020).
  * 33. S. Wang, et al., Nature Communications 12, 6961 (2021).
  * 34. J. Kim, Y. Oz, Journal of Statistical Mechanics: Theory and Experiment 2022, 073101 (2022).
  * 35. J. Kim, Y. Oz, Phys. Rev. A 106, 052424 (2022).
  * 36. R. Wiersema, C. Zhou, J. F. Carrasquilla, Y. B. Kim, SciPost Phys. 14, 147 (2023).
  * 37. S. H. Sack, R. A. Medina, A. A. Michailidis, R. Kueng, M. Serbyn, PRX Quantum 3, 020365 (2022).
  * 38. Y. Du, T. Huang, S. You, M.-H. Hsieh, D. Tao, npj Quantum Information 8, 62 (2022).
  * 39. S. Y.-C. Chen, Proceedings of the 2023 International Workshop on Quantum Classical Cooperative, QCCC ’23 (Association for Computing Machinery, New York, NY, USA, 2023), p. 17–20.
  * 40. U. Las Heras, U. Alvarez-Rodriguez, E. Solano, M. Sanz, Phys. Rev. Lett. 116, 230504 (2016).
  * 41. Z. Lu, P.-X. Shen, D.-L. Deng, Phys. Rev. Appl. 16, 044039 (2021).
  * 42. H. L. Tang, et al., PRX Quantum 2, 020310 (2021).
  * 43. L. Ding, L. Spector, Proceedings of the Genetic and Evolutionary Computation Conference Companion, GECCO ’22 (Association for Computing Machinery, New York, NY, USA, 2022), p. 2190–2195.
  * 44. Sequential quantum gate decomposer, https://github.com/rakytap/sequential-quantum-gate-decomposer (2021).
  * 45. A. Kitaev, A simple model of quantum holography (2015).
  * 46. S. Sachdev, J. Ye, Phys. Rev. Lett. 70, 3339 (1993).
  * 47. J. Maldacena, D. Stanford, Phys. Rev. D 94, 106002 (2016).
  * 48. G. Arfken, H. Weber, Mathematical Methods for Physicists pp. 957–961 (1985).
  * 49. D. Fogel, IEEE Transactions on Neural Networks 5, 3 (1994).
  * 50. T. H. W. Bäck, et al., Evolutionary Computation 31, 81 (2023).
  * 51. K. M. Nakanishi, K. Fujii, S. Todo, Phys. Rev. Res. 2, 043158 (2020).
  * 52. M. Ostaszewski, E. Grant, M. Benedetti, Quantum 5, 391 (2021).
  * 53. Line Search Methods (Springer New York, New York, NY, 2006), pp. 30–65.
  * 54. P. Rakyta, Z. Zimborás, Efficient quantum gate decomposition via adaptive circuit compression (2022).
  * 55. N. A. Nemkov, E. O. Kiktenko, I. A. Luchnikov, A. K. Fedorov, Quantum 7, 993 (2023).
  * 56. T. Volkoff, P. J. Coles, Quantum Science and Technology 6, 025008 (2021).
  * 57. I. Bengtsson, K. Zyczkowski, Geometry of Quantum States: An Introduction to Quantum Entanglement (Cambridge University Press, 2006).
  * 58. H.-Y. Huang, R. Kueng, J. Preskill, Nature Physics 16, 1050 (2020).
  * 59. S. K. Foong, S. Kanno, Phys. Rev. Lett. 72, 1148 (1994).
  * 60. C. Kelley, Iterative Methods for Optimization, Frontiers in Applied Mathematics (Society for Industrial and Applied Mathematics, 1999).
  * 61. D. P. Kingma, J. Ba, Adam: A method for stochastic optimization (2014).
  * 62. Y. Huang, Y. Gu, Phys. Rev. D 100, 041901 (2019).
  * 63. M. J. D. Powell, Mathematical Programming 38, 29 (1987).
  * 64. B. O’Gorman, W. J. Huggins, E. G. Rieffel, K. B. Whaley, Generalized swap networks for near-term quantum computing (2019).
  * 65. G.-L. R. Anselmetti, D. Wierichs, C. Gogolin, R. M. Parrish, New Journal of Physics 23, 113010 (2021).
  * 66. J. S. Kottmann, Quantum 7, 1073 (2023).
  * 67. A. Kukliansky, E. Younis, L. Cincio, C. Iancu, Qfactor: A domain-specific optimizer for quantum circuit instantiation (2023).
  * 68. X.-C. Wu, M. G. Davis, F. T. Chong, C. Iancu, 2021 International Conference on Rebooting Computing (ICRC) (2021), pp. 35–46.
  * 69. P. Rakyta, et al., Journal of Computational Physics 500, 112756 (2024).
  * 70. A. P. A. Zulehner, R. Wille, https://github.com/iic-jku/ibm_qx_mapping (2019).
  * 71. D. Wierichs, J. Izaac, C. Wang, C. Y.-Y. Lin, Quantum 6, 677 (2022).