# Towards large-scale quantum optimization solvers with few qubits

Marco Sciorilli<EMAIL_ADDRESS>Quantum Research Center, Technology
Innovation Institute, Abu Dhabi, UAE  Lucas Borges<EMAIL_ADDRESS>Quantum Research Center, Technology Innovation Institute, Abu Dhabi, UAE
Federal University of Rio de Janeiro, Caixa Postal 652, Rio de Janeiro, RJ
21941-972, Brazil  Taylor L. Patti  NVIDIA, Santa Clara, California 95051, USA
Diego García-Martín  Information Sciences, Los Alamos National Laboratory, Los
Alamos, NM 87545, USA  Quantum Research Center, Technology Innovation
Institute, Abu Dhabi, UAE  Giancarlo Camilo  Quantum Research Center,
Technology Innovation Institute, Abu Dhabi, UAE  Anima Anandkumar  Department
of Computing + Mathematical Sciences (CMS), California Institute of Technology
(Caltech), Pasadena, CA, 91125 USA  NVIDIA, Santa Clara, California 95051, USA
Leandro Aolita  Quantum Research Center, Technology Innovation Institute, Abu
Dhabi, UAE

###### Abstract

We introduce a variational quantum solver for combinatorial optimizations over
$m=\order{n^{k}}$ binary variables using only $n$ qubits, with tunable $k>1$.
The number of parameters and circuit depth display mild linear and sublinear
scalings in $m$, respectively. Moreover, we analytically prove that the
specific qubit-efficient encoding brings in a super-polynomial mitigation of
barren plateaus as a built-in feature. This leads to unprecedented quantum-
solver performances. For $m=7000$, numerical simulations produce solutions
competitive in quality with state-of-the-art classical solvers. In turn, for
$m=2000$, an experiment with $n=17$ trapped-ion qubits featured MaxCut
approximation ratios estimated to be beyond the hardness threshold $0.941$. To
our knowledge, this is the highest quality attained experimentally on such
sizes. Our findings offer a novel heuristics for quantum-inspired solvers as
well as a promising route towards solving commercially-relevant problems on
near term quantum devices.

††preprint: APS/123-QED

Combinatorial optimizations are ubiquitous in industry and technology [1]. The
potential of quantum computers for these problems has been extensively studied
[2, 3, 4, 5, 6, 7, 8, 9, 10]. However, it is unclear whether they will deliver
advantages in practice before fault-tolerant devices appear. With only
quadratic asymptotic runtime speed-ups expected in general [11, 12, 13] and
low clock-speeds [14, 15], a major challenge is the number of qubits required
for quantum solvers to become competive with classical ones. Current
implementations are restricted to noisy intermediate-scale quantum devices
[16], with variational quantum algorithms [17, 18] as a promising alternative.
These are heuristic models – based on parameterized quantum circuits – that,
although conceptually powerful, face inherent practical challenges [19, 20,
21, 22, 23, 24, 25]. Among them, hardware noise is particularly serious, since
its detrimental effect rapidly grows with the number of qubits. This can
flatten out the optimization landscape – causing exponentially-small gradients
(barren plateaus) [24] or underparametrization [25] – or render the algorithm
classically simulable [20]. Hence, near-term quantum optimization solvers are
unavoidably restricted to problem sizes that fit within a limited number of
qubits.

In view of this, interesting qubit-efficient schemes have been explored [26,
27, 28, 29, 30, 31, 32]. In Refs. [26, 27], two or three variables are encoded
into the (three-dimensional) Bloch vector of each qubit, allowing for a linear
space compression. In contrast, the schemes of [28, 29, 30, 31, 32] encode the
$m$ variables into a quantum register of size
$\mathcal{O}\big{(}\log(m)\big{)}$. However, such exponential compressions
both render the scheme classically simulable efficiently and seriously limit
the expressivity of the models [28, 31]. Moreover, in Refs. [28, 29, 30, 31,
32], binary problems are relaxed to quadratic programs. These simplifications
strongly affect the quality of the solutions. In addition, the measurements
required by those methods can be statistically demanding. For instance, in a
deployment with $m=3964$ variables [30], most of the measurement outcomes
needed did not occur and were replaced by classical random bits, leading to a
low quality solution compared to state-of-the-art solvers. To the best of our
knowledge, no experimental quantum solver has so far produced non-trivial
solutions to problems with $m$ beyond a few hundreds [8, 9, 10, 33].
Furthermore, the interplay between qubit-number compression, loss function
non-linearity, trainability, and solver performance in general is mostly
unknown.

Figure 1:  Quantum optimization solvers with polynomial space compression.
Encoding: An exemplary MaxCut (or weighted MaxCut) problem of $m=9$ vertices
(graph on the left) is encoded into 2-body Pauli-matrix correlations across
$n=3$ qubits (Q1, Q2, Q3). The colour code indicates which Pauli string
encodes which vertex. For instance, the binary variable $x_{1}$ of vertex 1 is
encoded in the expectation value of $Z_{1}\otimes Z_{2}\otimes\openone_{3}$,
supported on qubits 1 and 2, while $x_{9}$ is encoded in
$Y_{1}\otimes\openone_{2}\otimes Y_{3}$, over qubits 1 and 3 (see Eq. (1)).
This corresponds to a quadratic space compression of $m$ variables into
$n=\mathcal{O}(m^{1/2})$ qubits. More generally, $k$-body correlations can be
used to attain polynomial compressions of order $k$. The Pauli set chosen is
composed of three subsets of mutually-commuting Pauli strings. This allows one
to experimentally estimate all $m$ correlations using only 3 measurement
settings throughout. Quantum-classical optimization: we train a quantum
circuit parametrized by gate parameters $\bm{\theta}$ using a loss function
$\mathcal{L}$ of the Pauli expectation values that mimics the MaxCut (or
weighted MaxCut) objective function (see Eq. (2)). The variational Ansatz is a
brickwork circuit with the number of 2-qubit gates and variational parameters
scaling very mildly with $m$ (see Fig. 2), and circuit depth sublinear in $m$.
This makes both experimentally- and training-friendly (see Fig. 3). Solution:
once the circuit is trained, we read-out its output $\bm{x}$ from the
correlations across single-qubit measurement outcomes on its output state.
Finally, we perform an efficient classical bit-swap search around $\bm{x}$ to
find potential better solutions nearby. The result of that search,
$\bm{x}^{*}$, is the final output of our solver.

Here we explore this territory. We introduce a hybrid quantum-classical solver
for binary optimization problems of size $m$ polynomially larger than the
number of qubits $n$ used. This is an interesting regime in that the scheme is
highly qubit-efficient while at the same time preserving the classical
intractability in $m$, leaving room for potential quantum advantages. We
encode the $m$ variables into Pauli correlations across $k$ qubits, for $k$ an
integer of our choice. A parameterized quantum circuit is trained so that its
output correlations minimize a non-linear loss function suitable for gradient
descent. The solution bit string is then obtained via a simple classical post-
processing of the measurement outcomes, which includes an efficient local bit-
swap search to further enhance the solution’s quality. Moreover, a beneficial,
intrinsic by-product of our qubit-number compression is a super-polynomial
suppression of the decay of gradients, from barren plateaus of heights
$2^{-\Theta(m)}$ with single-qubit encodings to $2^{-\Theta(m^{1/k})}$ with
Pauli-correlation encodings. In turn, the circuit depth scales sublinearly in
$m$, as $\mathcal{O}(m^{1/2})$ for quadratic ($k=2$) compressions and
$\mathcal{O}(m^{2/3})$ for cubic ($k=3$) ones. All these features make our
scheme more experimentally- and training-friendly than previous quantum
schemes, leading to significantly higher quality of solutions.
For example, for $m=2000$ and $m=7000$ MaxCut instances, our numerical
solutions are competitive with those of semi-definite program relaxations,
including the powerful Burer-Monteiro algorithm. This is relevant as a basis
for quantum-inspired classical solvers. In addition, we deploy our solver on
IonQ hardware, observing an impressive performance even without quantum error
mitigation. For example, for a MaxCut instance with $m=2000$ vertices encoded
into $n=17$ trapped-ion qubits, we obtain estimated approximation ratios above
the hardness threshold $r\approx 0.941$. This is the highest quality reported
by an experimental quantum solver on sizes beyond a few tens for MaxCut [8,
33] and a few hundreds for combinatorial optimizations in general [10, 9]. Our
results open up a promising framework to develop competitive solvers for
large-scale problems with small quantum devices.

##  I Results

###  I.1 Quantum solvers with polynomial space compression

We solve combinatorial optimizations over $m=\order{n^{k}}$ binary variables
using only $n$ qubits, for $k$ a suitable integer of our choice. Such a
compression is achieved by encoding the variables into $m$ Pauli-matrix
correlations across multiple qubits. More precisely, with the short-hand
notation $[m]:=\\{1,2,\ldots m\\}$, let $\bm{x}:=\\{x_{i}\\}_{i\in[m]}$ denote
the string of optimization variables and choose a specific subset
$\Pi:=\\{\Pi_{i}\\}_{i\in[m]}$ of $m\leq 4^{n}-1$ traceless Pauli strings
$\Pi_{i}$, i.e. of $n$-fold tensor products of identity ($\openone$) or Pauli
($X$, $Y$, and $Z$) matrices, excluding the $n$-qubit identity matrix
$\openone^{\otimes n}$. We define a Pauli-correlation encoding (PCE) relative
to $\Pi$ as

$\displaystyle x_{i}:=\text{sgn}\big{(}\langle\Pi_{i}\rangle\big{)}\text{ for
all }i\in[m],$ (1)

where sgn is the sign function and
$\langle\Pi_{i}\rangle:=\bra{\Psi}\Pi_{i}\ket{\Psi}$ is the expectation value
of $\Pi_{i}$ over a quantum state $\ket{\Psi}$. We focus on strings with $k$
single-qubit traceless Pauli matrices. In particular, we consider encodings
$\Pi^{(k)}:=\big{\\{}\Pi^{(k)}_{1},\ldots,\Pi^{(k)}_{m}\big{\\}}$ where each
$\Pi^{(k)}_{i}$ is a permutation of either $X^{\otimes
k}\otimes\openone^{\otimes n-k}$, $Y^{\otimes k}\otimes\openone^{\otimes
n-k}$, or $Z^{\otimes k}\otimes\openone^{\otimes n-k}$ (see left panel of Fig.
1 for an example with $k=2$). That is, $\Pi^{(k)}$ is the union of 3 sets of
mutually-commuting strings. This is experimentally convenient, since only
three measurement settings are required throughout. Using all possible
permutations for the encoding yields $m=3\binom{n}{k}$. In this work, we deal
mostly with $k=2$ and $k=3$, corresponding to $m=\frac{3}{2}n(n-1)$ and
$m=\frac{1}{2}n(n-1)(n-2)$, respectively. The single-qubit encodings of [26,
27], in turn, correspond to PCEs with $k=1$.

The specific problem we solve is weighted MaxCut, a paradigmatic NP-hard
optimization problem over a weighted graph $G$, defined by a (symmetric)
adjacency matrix $W\in\mathbb{R}^{m\times m}$. Each entry $W_{ij}$ contains
the weight of an edge $(i,j)$ in $G$. The set $E$ of edges of $G$ consists of
all $(i,j)$ such that $W_{ij}\neq 0$. We denote by $|E|$ the cardinality of
$E$. The special case where all weights are either zero or one defines the
(still NP-hard) MaxCut problem, where each instance is fully specified by $E$
(see Sec. III.1). The goal of these problems is to maximize the total weight
of edges cut over all possible bipartitions of $G$. This is done by maximizing
the quadratic objective function $\mathcal{V}(\bm{x}):=\sum_{(i,j)\in
E}W_{ij}(1-\,x_{i}\,x_{j})$ (the cut value).

We parameterize the state in Eq. (1) as the output of a quantum circuit with
parameters $\bm{\theta}$, $\ket{\Psi}=\ket{\Psi(\bm{\theta})}$, and optimize
over $\bm{\theta}$ using a variational approach [17, 18]. As circuit Ansatz,
we use the brickwork architecture shown in Fig. 1 (see Sec. III.3 for details
on the variational Ansatz). The goal of the parameter optimization is to
minimize the non-linear loss function

$\displaystyle\mathcal{L}$ $\displaystyle=\\!\\!\sum_{(i,j)\in
E}\\!\\!\\!W_{ij}\tanh\big(\alpha\,\langle\Pi_{i}\rangle\big{missing})\tanh\big(\alpha\,\langle\Pi_{j}\rangle\big{missing})+\mathcal{L}^{(\text{reg})}.$
(2)

The first term corresponds to a relaxation of the binary problem where the
sign functions in Eq. (1) are replaced by smooth hyperbolic tangents, better-
suited for gradient-descent methods [27]. The second term,
$\mathcal{L}^{(\text{reg})}$ (see Sec. III.2), forces all correlators to go
towards zero, which is observed to improve the solver’s performance (see App.
A in the Supplementary Information). One way to interpret the latter is in
view of frustration effects, by virtue of which the correlator magnitudes
cannot all be arbitrarily close to one. For instance, in the case of two-body
correlators ($k=2$), we prove in App. B that, if one imposes
$|\langle\Pi_{i}\rangle|=c$ for all $i\in[m]$, then $c=\mathcal{O}(m^{-1/2})$
is required in order to be able to encode arbitrary bit strings. However, this
also restricts the $\tanh$ to its linear regime ($\tanh(z)\approx z$ for small
$z$), which is inconvenient for the training. Hence, to restore the non-linear
behaviour, we introduce the rescaling factor $\alpha=\mathcal{O}({c}^{-1})$
(see App. C).

Once the training is complete, the circuit output state is measured and a bit-
string $\bm{x}$ is obtained via Eq. (1). Then, as a classical post-processing
step, we perform one round of single-bit swap search (of complexity
$\mathcal{O}(\absolutevalue{E})$) around $\bm{x}$ in order to find potential
better solutions nearby (see Sec. III.3). The result of the search,
$\bm{x}^{*}$, with cut value $\mathcal{V}(\bm{x}^{*})$, is the final output of
our solver.

Our work differs from [28, 29, 30, 31, 32] in fundamental ways. First, as
mentioned, those studies focus mainly on exponential compressions in qubit
number. These are also possible with PCEs, since there are $4^{n}-1$ traceless
operators available. However, apart from automatically rendering the schemes
classically simulable efficiently [28, 31], exponential compressions strongly
limit the expressivity of the model, since $L$-depth circuits contain
$\mathcal{O}(L\times\log(m))$ parameters. This affects the quality of the
solutions. In contrast, our method operates manifestly in the regime of
classically intractable quantum circuits (see App. D). Second, as for
experimental feasibility, the previous schemes require the measurement of
probabilities that are (at best) of order $m^{-1}$. In contrast, our solver
requires expectation values of order $m^{-\lfloor k/2\rfloor/k}$ (see App. C).
For $k=2$ and $k=3$, e.g., this is respectively quadratically and cubically
better than $m^{-1}$. Third, while in [28, 29, 30, 31] the problems are
relaxed to quadratic programs [28], Eq. (2) defines a highly non-linear
optimization. These combined features lead to solutions notably superior to
those of previous schemes.

###  I.2 Circuit complexities and approximation ratios

Figure 2: Gate complexity and performance. Left: Number of two-qubits gates
needed for achieving an average estimated approximation ratio
$\overline{r}\geq 16/17\approx 0.941$ (over 250 non-trivial random MaxCut
instances and 5 random initializations per instance) without the local bit-
swap search (quantum-circuit’s output $\bm{x}$ alone) versus $m$, both for
quadratic and cubic compressions. A linear scaling is observed in both cases.
Right: Maximum $r$ (now including the local bit-swap search step) over random
initializations for three specific MaxCut instances of different sizes as
functions of the compression degree $k$ ($10$ random initializations were used
for $m=800$ and $m=2000$, and $5$ for $m=7000$). For a fair comparison, the
total number of parameters is kept the same for all $k$. The horizontal lines
denote the reported results of the leading gradient-based SDP solver [34]
(dotted lines) and the powerful Burer-Monteiro algorithm [35, 36] (dashed
lines). Remarkably, our solver outperforms the former in all cases and even
the latter for the $m=2000$ instance at $k=6$ and 7.

Here, we investigate the quantum resources (circuit depth, two-qubit gate
count, and number of variational parameters) required by our scheme. Due to
the strong reduction in qubit number, an increase in required circuit depth is
expected to maintain the same expressivity. We benchmark on graph instances
whose exact solution
$\mathcal{V}_{\text{max}}:=\max_{\bm{x}}\mathcal{V}(\bm{x})$ is unknown in
general. Therefore, we denote by
$r_{\text{exact}}:=\mathcal{V}(\bm{x}^{*})/\mathcal{V}_{\text{max}}$ the exact
approximation ratio and by
$r:=\mathcal{V}(\bm{x}^{*})/\mathcal{V}_{\text{best}}$ the estimated
approximation ratio based on the best known solution
$\mathcal{V}_{\text{best}}$ available (see Sec. III.3).
In Fig. 2 (left panel), we plot the gate complexity required to reach
$\overline{r}=16/17\approx 0.941$ without doing the final local search step
(to capture the resource scaling exclusively due to the quantum subroutine) on
non-trivial random MaxCut instances of increasing sizes, for the encodings
$\Pi^{(2)}$ and $\Pi^{(3)}$. For $r_{\text{exact}}$, this value gives the
threshold for worst-case computational hardness. By non-trivial instances we
mean instances post-selected to discard easy ones (see Sec. III.3). The
results suggest that the number of gates scales approximately linearly with
$m$. The same holds also for the number of variational parameters, which is
proportional to the number of gates. In turn, the number of circuit layers
scales as $\mathcal{O}(m/n)$. For quadratic and cubic compressions, e.g., this
corresponds to $\mathcal{O}(m^{1/2})$ and $\mathcal{O}(m^{2/3})$,
respectively. These surprisingly mild scalings translate directly into
experimental feasibility and model-training ease. In fact, we observe (see
App. E) that the number of epochs needed for training also scale linearly with
$m$. Moreover, in App. F, we prove upper bounds on the numbers of measurements
required to estimate $\mathcal{L}(\bm{\theta})$. For $k=2$ and $k=3$, e.g.,
these bounds coincide and give
$\tilde{\mathcal{O}}\left(m\,(6|E|+m)^{2}\right)$.

In Fig. 2 (right), in turn, we plot solution qualities versus $k$, for three
MaxCut instances from the benchmark set `Gset` [37] (see Sec. III.3). The
total number of variational parameters is fixed by $m$ (or as close to $m$ as
allowed by the circuit ansatz) for a fair comparison, with the circuit depths
adjusted accordingly for each $k$. In all cases, $r$ increases with $k$ up to
a maximum, after which the performance degrades. This is consistent with a
limit in compression capability before compromising the model’s expressivity,
as expected. Remarkably, the results indicate that our solutions are
competitive with those of state-of-the-art classical solvers, such as the
leading gradient-based SDP solver [34], based on the interior points method,
and even the Burer-Monteiro algorithm [35, 36], based on non-linear
programming. Importantly, while our solver performs a single optimization
followed by a single-bit swap search, the Burer-Monteiro algorithm includes
multiple re-optimizations and two-bit swap searches (see App. G). This
highlights the potential for further improvements of our scheme. All in all,
the impressive performance seen in Fig. 2 is not only relevant for quantum
solvers, but also suggests our scheme as an interesting heuristic for quantum-
inspired classical solvers.

Figure 3:  Loss-function variance decay. Main: Average sample variance
$\overline{\text{Var}(\mathcal{L})}$ of $\mathcal{L}$, normalized by
$\alpha^{4}\sum_{(i.j)\in E}w_{ij}^{2}$, after plateauing, for the encodings
$\Pi^{(1)}$, $\Pi^{(2)}$, and $\Pi^{(3)}$, as a function of $m$ (in log-linear
scale). The agreement with the analytical expression in Eq. (3) is excellent
(the dashed, blue curve corresponds to the first term of the equation, which
decreases as $2^{-2n}$). Since $n=\mathcal{O}(m^{1/k})$, this translates into
a super-polynomial suppression in $m$ of the decay speed of
$\overline{\text{Var}(\mathcal{L})}$ for $k>1$. Inset: average variances of
the entries of the gradient of $\mathcal{L}$ as functions of the number of
layers, for quadratic compression ($\Pi^{(2)}$). Each curve corresponds to a
different $n$. Note the decay of the plateau (rightmost) values with $n$. The
black crosses indicate the depths needed to reach average approximation ratios
$>0.941$ (computational-hardness threshold) with the quantum-circuit’s output
$\bm{x}$ alone, i.e. excluding the final bit-swap search. Remarkably, in all
cases such ratios are attained before the variances have converged to their
asymptotic, steady values.

###  I.3 Intrinsic mitigation of barren plateaus

Another appealing feature of our solver emerging from the qubit-number
reduction is an intrinsic mitigation of the infamous barren plateau (BP)
problem [23, 38, 39, 40, 41], which constitutes one of the main challenges for
training variational quantum algorithms. BPs are characterized by a vanishing
expectation value of $\nabla\mathcal{L}$ over random parameter initializations
and an exponential decay (in $n$) of its variance. This jeopardizes the
applicability of variational quantum algorithms in general [19]. For instance,
the gradient variances of a two-body Pauli correlator on the output of
universal 1D brickwork circuits are known to plateau at levels exponentially
small in $n$ for circuit depths of about $10\times n$ [23]. Alternatively, BPs
can equivalently be defined in terms of an exponentially vanishing variance of
$\mathcal{L}$ itself (instead of its gradient) [42]. This is often more
convenient for analytical manipulations.

In App. I we prove that, if the random parameter initializations make the
circuits sufficiently random (namely, induce a Haar measure over the special
unitary group), the variance of $\mathcal{L}$ is given by

${\rm Var}(\mathcal{L})=\frac{\alpha^{4}}{d^{2}}\sum_{(i,j)\in
E}w^{2}_{ij}+\mathcal{O}\left(\frac{\alpha^{6}}{d^{3}}\right)\,,$ (3)

where $d=2^{n}$ is the Hilbert-space dimension. Interestingly, the leading
term in Eq. (3) appears also if one only assumes the circuits to form a
$4$-design, but it is then not clear how to bound the higher-order terms
without the full Haar-randomness assumption. However, we suspect that the
latter is indeed not necessary. In practice, the unitary-design assumption is
always approximately satisfied for sufficiently deep universal quantum
circuits. In line with that, for our loss function, we empirically observe
convergence to Eq. (3) at circuit depths of about $8.5\times n$. This is
illustrated in Fig. 3 for linear, quadratic, and cubic compressions, where we
plot the average sample variance $\overline{\text{Var}(\mathcal{L})}$ of
$\mathcal{L}$ over 100 non-trivial random MaxCut instances and 100 random
parameter initializations per instance, as a function of $m$. In contrast, the
depth needed to reach $r>0.941$ on average with the circuit’s output alone is
about $1.05\times n$ (see figure inset).

One observes an excellent agreement between
$\overline{\text{Var}(\mathcal{L})}$ and the first term of Eq. (3) for large
$m$. As $m$ decreases, small discrepancies appear, specially for $k=2$ and
$k=3$. This can be explained by noting that $\alpha\sim 1.5$ for $k=1$ whereas
$\alpha\sim 1.5\times n$ for $k=2$ and $k=3$ (see App. C), so that the second
term in (3) scales as $\frac{1}{2^{3n}}$ for the former but as
$\frac{n^{6}}{2^{3n}}$ for the latter. Hence, as $m$ (and so $n$) decreases,
that term requires smaller $m$ to become non-negligible for the former than
for the latter. Remarkably, the scaling
$\overline{\text{Var}(\mathcal{L})}\in\Theta(\alpha^{4}\,2^{-2n})$ in $n$
translates into a super-polynomial suppression of the decay speed in $m$ when
compared to single-qubit (linear) encodings. This means, for instance, that
quadratic encodings feature
$\overline{\text{Var}(\mathcal{L})}\in\Theta(\alpha^{4}\,2^{-2\sqrt{m}})$,
instead of
$\overline{\text{Var}(\mathcal{L})}\in\Theta(\alpha^{4}\,2^{-2\,m})$ displayed
by linear encodings. Importantly, the scaling obtained still represents a
super-polynomial decay in $m$. Yet, the enhancement obtained makes a
tremendous difference in practice, as shown in the figure by the orders of
magnitude separating the three curves.

Figure 4: Trapped-ion experimental implementation. Main: estimated
approximation ratios given by our scheme deployed on IonQ’s device Aria-1 as
functions of the number of measurements (per each of the three measurement
settings). Three problem instances (see main text for details) are shown: one
weighted MaxCut instance of $m=512$, solved with quadratic compression using
$n=19$ qubits (purple pentagons), one MaxCut instance of $m=800$, solved with
quadratic compression on $n=23$ qubits (red circles) and cubic compression on
$n=13$ qubits (yellow triangles), and another MaxCut instance of $m=2000$,
solved with cubic compression on $n=17$ qubits (blue squares). The two
horizontal lines indicate the Goemans-Williamson threshold (black dotted), at
$r\approx 0.878$, and the worst-case computational-hardness threshold (black
dashed), at $r=16/17$. Inset: experimental and numerically-simulated loss
functions $\mathcal{L}$ (at fixed, optimized parameters) versus number of
measurement shots per measurement setting (same shot range as in main figure),
for the $m=800$ instance with quadratic and cubic compressions.

###  I.4 Experimental deployment on quantum hardware

We experimentally demonstrate our quantum solver on IonQ’s trapped-ion device
Aria-1, for two MaxCut instances of $m=800$ and $2000$ vertices and a weighted
MaxCut instance of $m=512$ vertices. Details on the hardware and model
training are provided in Sec. III.4, while the choice of instances is detailed
in Sec. III.3 (see Table 1). We optimize the circuit parameters offline via
classical simulations and experimentally deploy the pre-trained circuit. Fig.
4 depicts the obtained approximation ratios for each instance as a function of
the number of measurement shots, employing both the quadratic and cubic Pauli-
correlation encodings, $\Pi^{(2)}$ and $\Pi^{(3)}$, respectively. For each
instance, we collected enough statistics for the approximation ratio to
converge (see figure inset). Gate infidelities strongly limit the circuit
size. We found a good trade-off at $95$ two-qubit gates altogether. This is
below the number required for these instance sizes according to Fig. 2 (left),
especially for $k=2$. However, remarkably, our solver still returns solutions
of higher quality than the Göemans-Williamson bound in all cases and even than
the worst-case hardness threshold in three out of the four experiments
performed. This is the first-ever quantum experiment to produce such high-
quality solutions for these sizes. As a reference, the largest MaxCut instance
experimentally solved with QAOA [2] has size $m=414$ and average and maximal
approximation ratios 0.57 and 0.69, respectively (see Table IV in Ref. [33]).

##  II Conclusions and Discussion

We introduced a scheme for solving binary optimizations of size $m$
polynomially larger than the number of qubits used. Pauli correlations across
few qubits encode each binary variable. The circuit depth is sublinear in $m$
while the numbers of parameters and training epochs approximately linear in
$m$. Moreover, the qubit-number compression brings in the remarkable by-
product of significantly suppressing the decay in $m$ of the variances of the
loss function (and its gradient, for that matter), as we have analytically
proven. These features, together with an educated choice of non-linear loss
function, allow us to solve large, computationally non-trivial instances with
unprecedentedly-high quality. Numerically, our solutions for $m=2000$ and
$m=7000$ MaxCut instances are competitive with those of semi-definite program
relaxations, including the powerful Burer-Monteiro algorithm. Experimentally,
in turn, for a deployment on $17$ trapped-ion qubits, we estimated
approximation ratios beyond the worst-case computational hardness threshold
$0.941$ for a non-trival MaxCut instance with $m=2000$ vertices. To our
knowledge, this the highest solution quality ever reported experimentally on
such instance sizes. We stress that these results are based on raw
experimental data, without any quantum error mitigation procedure to the
observables measured. Yet, our method is indeed well-suited for standard error
mitigation techniques [43, 44, 45], the use of which can enhance the solver’s
performance even further. In turn, although we have focused on quadratic
unconstrained binary optimization (QUBO) problems, the technique can be
straightforwardly extended to generic polynomial unconstrained binary
optimizations (PUBOs) [46] without any increase in qubit numbers. This is in
contrast to conventional PUBO-to-QUBO reformulations, which incur in expensive
overheads in extra qubits [47]. Finally, as with most variational quantum
algorithms, an important open question is the time complexity of
experimentally training the model. Particularly challenging is the number of
measurement runs required for estimating gradients with enough statistical
precision. In the frequent case where $|E|$ is linear in $m$, our analysis
gives a worst-case upper bound $\tilde{\mathcal{O}}(m^{3})$ to the sample
complexity of estimating the loss function. However, the actual scaling in
practice may be significantly better; and further improvement is even possible
through optimization of the hyperparameters $\alpha$ and $\beta$ on a case-by-
case basis. Another important benchmark for our method is fully-classical
schemes where the quantum circuit in our solver is replaced by a generative
neural network (see App. D). However, the fact that the quantum circuit in
question cannot be classically simulated efficiently in $m$ offers promising
prospects for future explorations. All in all, our framework offers an
interesting arena to develop competitive solvers for large-scale problems,
both with small quantum devices in the near term or with quantum-inspired
classical techniques.

##  III Methods

###  III.1 MaxCut problems

The weighted MaxCut problem is an ubiquitous combinatorial optimization
problem. It is a graph partitioning problem defined on weighted undirected
graphs $G=(V,E)$ whose goal is to divide the $m$ vertices in $V$ into two
disjoint subsets in a way that maximizes the sum of edge weights $W_{ij}$
shared by the two subsets – the so-called cut value. If the graph $G$ is
unweighted, that is, if $W_{ij}=1$ or $W_{ij}=0$ for every edge $(i,j)\in E$,
the problem is referred to simply as MaxCut. By assigning a binary label
$x_{i}$ to each vertex $i\in V$, the problem can be mathematically formulated
as the binary optimization

$\underset{\bm{x}\in\\{-1,1\\}^{m}}{\text{maximize}}\hskip
5.69046pt\sum_{i,j\in[m]}W_{ij}(1-\,x_{i}\,x_{j})\,.$ (4)

Since $\sum_{i,j\in[m]}W_{ij}$ is constant over $\bm{x}$, Eq. (4) can be
rephrased as a minimization of the objective function
$\bm{x}^{\text{T}}W\bm{x}$. This specific format is known as a quadratic
unconstrained binary optimization (QUBO). For generic graphs, solving MaxCut
exactly is NP-hard [48]. Moreover, even approximating the maximum cut to a
ratio $r_{\text{exact}}>\frac{16}{17}\approx 0.941$ is NP-hard [49, 50]. In
turn, the best-known polynomial-time approximation scheme is the Goemans-
Williamson (GW) algorithm [51], with a worst-case ratio
$r_{\text{exact}}\approx 0.878$. Under the Unique Games Conjecture, this is
the optimal achievable by an efficient classical algorithm with worst-case
performance guarantees. If, however, one does not require performance
guarantees, there exist powerful heuristics that in practice produce cut
values often higher than those of the GW algorithm. Two examples are discussed
in Best solutions known in Sec. III.3.

###  III.2 Regularization term

The regularization term in Eq. (2) penalizes large correlator values, thereby
forcing the optimizer to remain in the correlator domain where all possible
bit string solutions are expressible. Its explicit form is

$\displaystyle\mathcal{L}^{(\text{reg})}=\beta\,\nu\left[\frac{1}{m}\sum_{i\in
V}\tanh\big(\alpha\,\langle\Pi_{i}\rangle\big{missing})^{2}\right]^{2}.$ (5)

The factor $1/m$ normalizes the term in square brackets to $\order{1}$. The
parameter $\nu$ is an estimate of the maximum cut value: it sets the overall
scale of $\mathcal{L}^{(\text{reg})}$ so that it becomes comparable to the
first term in Eq. (2). For weighted MaxCut, we use the Poljak-Turzík lower
bound $\nu=w(G)/2+w(T_{\text{min}})/4$ [52], where $w(G)$ and
$w(T_{\text{min}})$ are the weights of the graph and of its minimum spanning
tree, respectively. For MaxCut, this reduces to the Edwards-Erdös bound [53]
$\nu=|E|/2+(m-1)/4$. Finally, $\beta$ is a free hyperparameter of the model,
which we optimize over random graphs to get $\beta=1/2$. Such optimizations
systematically show increased approximation ratios due to the presence of
$\mathcal{L}^{(\text{reg})}$ in Eq. (2) (see App. A).

###  III.3 Numerical details

Choice of instances. The numerical simulations of Figs. 2 (left) and Fig. 3
were performed on random MaxCut instances generated with the well-known `rudy`
graph-generator [54] post-selected so as to filter out easy instances. The
post-selection consisted in discarding graphs with less than 3 edges per node
on average or those for which a random cut gives an approximation ratio
$r>0.82$. The latter is sufficiently far from the Goemans-Williamson ratio
$0.878$ while still allowing efficient generation. For the numerics in Fig. 2
(right) and the experimental deployment in Fig. 4 we used 6 graphs from
standard benchmarking sets: the former used the G14, G23, and G60 MaxCut
instances from the `Gset` repository [37], while the latter used G1 and G35
from `Gset` and the weighted MaxCut instance pm3-8-50 from the `DIMACS`
library [55] (recently employed also in [27]). Their features are summarized
in Table 1.

Best solutions known. For the generated instances, the best solution is taken
as the one with the highest cut value between the (often coinciding) solutions
produced by two classical heuristics, namely the Burer-Monteiro [35] and the
Breakout Local Search [56] algorithms. For the instances from benchmarking
sets, we considered instead the best known documented solution. The
corresponding cut value, $\mathcal{V}_{\text{best}}$, is used to define the
approximation ratio achieved by the quantum solution $\bm{x}^{*}$, namely
$r=\mathcal{V}(\bm{x}^{*})/\mathcal{V}_{\text{best}}$.

Graph | $m$ | $|E|$ | $W_{ij}$ | Type | Use
---|---|---|---|---|---
pm3-8-50 | 512 | 1536 | $\pm 1$ |  $3D$ torus grid | Experiment
G1 | 800 | 19176 | $1$ | random | Experiment
G14 | 800 | 4694 | $1$ | planar | Numerics
G23 | 2000 | 19990 | $1$ | random | Numerics
G35 | 2000 | 11778 | $1$ | planar | Experiment
G60 | 7000 | 17148 | $1$ | random | Numerics
Table 1: Benchmark instances used in this work. Apart from the the number of
vertices, edges, and edge weights, we also include the type of graph as well
as its use.

Variational Ansatz. As circuit Ansatz, we used the brickwork architecture
shown in Fig. 1, with layers of single-qubit rotations, parameterized by a
single angle, followed by a layer of Mølmer-Sørensen (MS) two-qubit gates,
each with three variational parameters. Each single-qubit gate layer contains
rotations around a single direction (X, or Y, or Z), one at a time,
sequentially. Furthermore, we observed that many of the other commonly used
parameterized gate displays the same numerical scalings up to a constant.

Quantum-circuit simulations. The classical simulations of quantum circuits
have been done using two libraries: `Qibo` [57, 58] for exact state-vector
simulations of systems up to 23 qubits, and `Tensorly-Quantum` [59] for
tensor-network simulations of larger qubit systems.

Optimization of circuit parameters. Two optimizers were used for the model
training. SLSQP from the `scipy` library was used for systems small enough to
calculate the gradient using finite differences. In all other cases we used
Adam from the `torch`/`tensorflow` libraries, leveraging automatic
differentiation to speed up computational time. As a stopping criterion for
Adam, we halted the training after $50$ steps whose cumulative improvement to
the loss function was less then $0.01$. For both optimizers, the default
optimization parameters were used.

Classical bit-swap search as post-processing step. As mentioned, a single
round of local bit-swap search is performed on the bit string $\bm{x}$ output
by the trained quantum circuit. This consists of sequentially swapping each
bit of $\bm{x}$ and computing the cut value of the resulting bit string. If
the cut value improves, we retain the change. Else, the local bit flip is
reverted. There are altogether $\Theta(m)$ local bit flips. A bit flip on
vertex $i$ affects $\Theta(d(i))$ edges, with $d(i)$ the degree of the vertex.
Hence, an update of only $\Theta(d(i))$ terms in $\mathcal{V}(\bm{x})$ is
required per bit flip. The total complexity of the entire round is thus
$\Theta(\absolutevalue{E})$.

###  III.4 Experimental details

|  |  |  |  |  |  $r$ |
---|---|---|---|---|---|---|---
Graph | $k$ | $n$ | 1-q | 2-q | Epochs | Sim. | Exp.
pm3-8-50 | 2 | 19 | 199 | 90 | 13485 | 0.967 | 0.921
G1 | 2 | 24 | 192 | 36 | 4027 | 0.954 | 0.957
G1 | 3 | 13 | 170 | 36 | 2022 | 0.940 | 0.965
G35 | 3 | 17 | 193 | 88 | 4100 | 0.935 | 0.943
Table 2: Details about the experimentally deployed instances. For each
instance, we display $k$, $n$, the 1-qubit and 2-qubit gate counts, and number
of optimization epochs used during classical training. The last two columns
report the approximation ratios given by the classical simulation of the
noiseless circuit (sim.) and the best one observed in the experiment (exp.).
We note that the latter ratios lie $3$ standard deviations away from the
average solution obtained via a single-bit search over a randomly picked bit
string (see Sec. H in the SI).

Hardware details.  The experiments were deployed in the IonQ Aria-1 quantum
hardware, a 25-qubit device with all-to-all connectivity based on trapped
ytterbium ions. The VQA architecture of Fig. 1 was adapted accordingly to
hardware native gates. We used alternating layers of partially entangling
Mølmer-Sørensen (MS) gates and, depending on the experiment, rotation layers
composed of one or two native single-qubit rotations GPI and GPI2 (see Table
2). Since the $z$ rotation is done virtually on the IonQ Aria chip,
parameterized RZ rotation were also added at the end of every rotation layer
without any extra gate cost.

Resource analysis. We run a total of four experimental deployments. The three
selected instances were trained using exact classical simulation with Adam
optimizer, as detailed in Sec. III.3. In an attempt to get the best possible
solution within the limited depth constraints of the hardware, the stopping
criteria was relaxed to allow $150$ non-improving steps. This resulted in a
total number of training epochs considerably larger then the average case
scenario (see App. G). Table 2 reports the precise quantum (number of qubits
and gate count) and classical (number of epochs) resources, as well as the
observed results.

##  IV Acknowledgments

The authors would like to thank Marco Cerezo and Tobias Haug for helpful
conversations. D.G.M. is supported by Laboratory Directed Research and
Development (LDRD) program of LANL under project numbers 20230527ECR and
20230049DR. At CalTech, A.A. is supported in part by the Bren endowed chair,
and Microsoft, Google, Adobe faculty fellowships.

## References

  * Korte  _et al._ [2011] B. H. Korte, J. Vygen, B. Korte, and J. Vygen,  _Combinatorial optimization_ , Vol. 1 (Springer, 2011).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm, arXiv (2014), arXiv:1411.4028 .
  * Guerreschi and Matsuura [2019] G. G. Guerreschi and A. Y. Matsuura, QAOA for max-cut requires hundreds of qubits for quantum speed-up, Scientific Reports 9, 6903 (2019).
  * Akshay  _et al._ [2020] V. Akshay, H. Philathong, M. Morales, and J. Biamonte, Reachability deficits in quantum approximate optimization, Physical Review Letters 124, 10.1103/physrevlett.124.090504 (2020).
  * Akshay  _et al._ [2021] V. Akshay, D. Rabinovich, E. Campos, and J. Biamonte, Parameter concentrations in quantum approximate optimization, Physical Review A 104, 10.1103/physreva.104.l010401 (2021).
  * Wurtz and Love [2021] J. Wurtz and P. Love, Maxcut quantum approximate optimization algorithm performance guarantees for $p>$ 1, Physical Review A 103, 042612 (2021).
  * Pagano  _et al._ [2020] G. Pagano, A. Bapat, P. Becker, K. S. Collins, A. De, P. W. Hess, H. B. Kaplan, A. Kyprianidis, W. L. Tan, C. Baldwin, L. T. Brady, A. Deshpande, F. Liu, S. Jordan, A. V. Gorshkov, and C. Monroe, Quantum approximate optimization of the long-range ising model with a trapped-ion quantum simulator, Proceedings of the National Academy of Sciences 117, 25396 (2020).
  * Harrigan  _et al._ [2021] M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, _et al._ , Quantum approximate optimization of non-planar graph problems on a planar superconducting processor, Nature Physics 3, 1745 (2021).
  * Ebadi  _et al._ [2022] S. Ebadi, A. Keesling, M. Cain, T. T. Wang, H. Levine, D. Bluvstein, G. Semeghini, A. Omran, J. Liu, R. Samajdar, X.-Z. Luo, B. Nash, X. Gao, B. Barak, E. Farhi, S. Sachdev, N. Gemelke, L. Zhou, S. Choi, H. Pichler, S. Wang, M. Greiner, V. Vuletic, and M. D. Lukin, Quantum optimization of maximum independent set using rydberg atom arrays (2022).
  * Yarkoni  _et al._ [2022] S. Yarkoni, E. Raponi, T. Bäck, and S. Schmitt, Quantum annealing for industry applications: introduction and review, Reports on Progress in Physics 85, 104001 (2022).
  * Dürr and Hoyer [1996] C. Dürr and P. Hoyer, A quantum algorithm for finding the minimum, CoRR quant-ph/9607014 (1996).
  * Ambainis [2004] A. Ambainis, Quantum search algorithms, ACM SIGACT News 35, 22 (2004).
  * Montanaro [2020] A. Montanaro, Quantum speedup of branch-and-bound algorithms, Physical Review Research 2, 013056 (2020).
  * Babbush  _et al._ [2021] R. Babbush, J. R. McClean, M. Newman, C. Gidney, S. Boixo, and H. Neven, Focus beyond quadratic speedups for error-corrected quantum advantage, PRX Quantum 2, 010103 (2021).
  * Campbell  _et al._ [2019] E. Campbell, A. Khurana, and A. Montanaro, Applying quantum algorithms to constraint satisfaction problems, Quantum 3, 167 (2019).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Cerezo  _et al._ [2021a] 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 3, 625 (2021a).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, _et al._ , Noisy intermediate-scale quantum algorithms, Reviews of Modern Physics 94, 015004 (2022).
  * Cerezo  _et al._ [2023] M. Cerezo, M. Larocca, D. García-Martín, N. Diaz, P. Braccia, E. Fontana, M. S. Rudolph, P. Bermejo, A. Ijaz, S. Thanasilp, _et al._ , Does provable absence of barren plateaus imply classical simulability? or, why we need to rethink variational quantum computing (2023), arXiv:2312.09121 [quant-ph] .
  * Stilck França and Garcia-Patron [2021] D. Stilck França and R. Garcia-Patron, Limitations of optimization algorithms on noisy quantum devices, Nature Physics 17, 1221 (2021).
  * Bittel and Kliesch [2021] L. Bittel and M. Kliesch, Training variational quantum algorithms is np-hard, Physical review letters 127, 120502 (2021).
  * Anschuetz and Kiani [2022] E. R. Anschuetz and B. T. Kiani, Quantum variational algorithms are swamped with traps, Nature Communications 13, 7760 (2022).
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 4812 (2018).
  * Wang  _et al._ [2021a] 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 12, 6961 (2021a).
  * García-Martín  _et al._ [2023a] D. García-Martín, M. Larocca, and M. Cerezo, Effects of noise on the overparametrization of quantum neural networks, arXiv preprint arXiv:2302.05059 (2023a).
  * Fuller  _et al._ [2021] B. Fuller, C. Hadfield, J. R. Glick, T. Imamichi, T. Itoko, R. J. Thompson, Y. Jiao, M. M. Kagele, A. W. Blom-Schieber, R. Raymond, and A. Mezzacapo, Approximate solutions of combinatorial problems via quantum relaxations (2021), arXiv:2111.03167 [quant-ph] .
  * Patti  _et al._ [2022] T. L. Patti, J. Kossaifi, A. Anandkumar, and S. F. Yelin, Variational quantum optimization with multibasis encodings, Physical Review Research 4, 033142 (2022).
  * Tan  _et al._ [2021] B. Tan, M.-A. Lemonde, S. Thanasilp, J. Tangpanitanon, and D. G. Angelakis, Qubit-efficient encoding schemes for binary optimisation problems, Quantum 5, 454 (2021).
  * Huber  _et al._ [2023] E. X. Huber, B. Y. L. Tan, P. R. Griffin, and D. G. Angelakis, Exponential qubit reduction in optimization for financial transaction settlement (2023), arXiv:2307.07193 [quant-ph] .
  * Leonidas  _et al._ [2023] I. D. Leonidas, A. Dukakis, B. Tan, and D. G. Angelakis, Qubit efficient quantum algorithms for the vehicle routing problem on nisq processors (2023), arXiv:2306.08507 [quant-ph] .
  * Tene-Cohen  _et al._ [2023] Y. Tene-Cohen, T. Kelman, O. Lev, and A. Makmal, A variational qubit-efficient maxcut heuristic algorithm (2023), arXiv:2308.10383 [quant-ph] .
  * Perelshtein  _et al._ [2023] M. R. Perelshtein, A. I. Pakhomchik, A. A. Melnikov, M. Podobrii, A. Termanova, I. Kreidich, B. Nuriev, S. Iudin, C. W. Mansell, and V. M. Vinokur, Nisq-compatible approximate quantum algorithm for unconstrained and constrained discrete optimization, Quantum 7, 1186 (2023).
  * Abbas  _et al._ [2023] A. Abbas, A. Ambainis, B. Augustino, A. Bärtschi, H. Buhrman, C. Coffrin, G. Cortiana, V. Dunjko, D. J. Egger, B. G. Elmegreen, N. Franco, 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, V. Markov, 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 (2023), arXiv:2312.02279 [quant-ph] .
  * Choi and Ye [2000] C. Choi and Y. Ye, Solving sparse semidefinite programs using the dual scaling algorithm with an iterative solver, https://web.stanford.edu/~yyye/yyye/cgsdp1.pdf (2000).
  * Burer and Monteiro [2003] S. Burer and R. D. C. Monteiro, A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization, Mathematical Programming 10.1007/s10107-002-0352-8 (2003).
  * Dunning  _et al._ [2018] I. Dunning, S. Gupta, and J. Silberholz, What works best when? a systematic evaluation of heuristics for Max-Cut and QUBO, INFORMS Journal on Computing 30, 608 (2018).
  * Ye [shed] Y. Ye, Gset test problems, https://web.stanford.edu/~yyye/yyye/Gset/ (unpublished).
  * Cerezo  _et al._ [2021b] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nature communications 12, 1791 (2021b).
  * Wang  _et al._ [2021b] 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 12, 10.1038/s41467-021-27045-6 (2021b).
  * Fontana  _et al._ [2023] 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 (2023), arXiv:2309.07902 [quant-ph] .
  * Ragone  _et al._ [2023] M. Ragone, B. N. Bakalov, F. Sauvage, A. F. Kemper, C. O. Marrero, M. Larocca, and M. Cerezo, A unified theory of barren plateaus for deep parametrized quantum circuits (2023), arXiv:2309.09342 [quant-ph] .
  * Arrasmith  _et al._ [2022] A. Arrasmith, Z. Holmes, M. Cerezo, and P. J. Coles, Equivalence of quantum barren plateaus to cost concentration and narrow gorges, Quantum Science and Technology 7, 045015 (2022).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Phys. Rev. X 7, 021050 (2017).
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error mitigation for short-depth quantum circuits, Phys. Rev. Lett. 119, 180509 (2017).
  * Cai  _et al._ [2023] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O’Brien, Quantum error mitigation (2023), arXiv:2210.00921 [quant-ph] .
  * Chermoshentsev  _et al._ [2022] D. A. Chermoshentsev, A. O. Malyshev, M. Esencan, E. S. Tiunov, D. Mendoza, A. Aspuru-Guzik, A. K. Fedorov, and A. I. Lvovsky, Polynomial unconstrained binary optimisation inspired by optical simulation (2022), arXiv:2106.13167 [quant-ph] .
  * Biamonte [2008] J. D. Biamonte, Nonperturbative k-body to two-body commuting conversion hamiltonians and embedding problem instances into ising spins, Physical Review A 77, 10.1103/physreva.77.052331 (2008).
  * Karp [1972] R. M. Karp, Reducibility among combinatorial problems, in  _Complexity of Computer Computations: Proceedings of a symposium on the Complexity of Computer Computations_, edited by R. E. Miller, J. W. Thatcher, and J. D. Bohlinger (Springer US, Boston, MA, 1972) pp. 85–103.
  * Håstad [2001] J. Håstad, Some optimal inapproximability results, J. ACM 48, 798–859 (2001).
  * Trevisan  _et al._ [2000] L. Trevisan, G. B. Sorkin, M. Sudan, and D. P. Williamson, Gadgets, approximation, and linear programming, SIAM Journal on Computing 29, 2074 (2000), https://doi.org/10.1137/S0097539797328847 .
  * Goemans and Williamson [1995] M. X. Goemans and D. P. Williamson, Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming, J. ACM 42, 1115–1145 (1995).
  * Poljak and Turzík [1986] S. Poljak and D. Turzík, A polynomial time heuristic for certain subgraph optimization problems with guaranteed worst case bound, Discrete Mathematics 58, 99 (1986).
  * Bylka  _et al._ [1999] S. Bylka, A. Idzik, and Z. Tuza, Maximum cuts: improvementsand local algorithmic analogues of the edwards-erdos inequality, Discrete Mathematics 194, 39 (1999).
  * Rinaldi [shed] G. Rinaldi, Rudy graph generator, http://www-user.tu-chemnitz.de/~helmberg/rudy.tar.gz (unpublished).
  * Pataki and Schmieta [shed] G. Pataki and S. H. Schmieta, The DIMACS library of mixed semidefinite-quadratic-linear programs, http://archive.dimacs.rutgers.edu/Challenges/Seventh/Instances/ (unpublished).
  * Benlic and Hao [2013] U. Benlic and J.-K. Hao, Breakout local search for the Max-Cut problem, Engineering Applications of Artificial Intelligence 26, 1162 (2013).
  * Efthymiou  _et al._ [2021] S. Efthymiou, S. Ramos-Calderer, C. Bravo-Prieto, A. Pérez-Salinas, D. García-Martín, A. Garcia-Saez, J. I. Latorre, and S. Carrazza, Qibo: a framework for quantum simulation with hardware acceleration, Quantum Science and Technology 7, 015018 (2021).
  * Efthymiou  _et al._ [2022] S. Efthymiou, M. Lazzarin, A. Pasquale, and S. Carrazza, Quantum simulation with just-in-time compilation, Quantum 6, 814 (2022).
  * Patti  _et al._ [2021] T. L. Patti, J. Kossaifi, S. F. Yelin, and A. Anandkumar, Tensorly-quantum: Quantum machine learning with tensor methods (2021), arXiv:2112.10239 [quant-ph] .
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Quantum supremacy using a programmable superconducting processor, Nature (London) 574, 505 (2019), arXiv:1910.11333 [quant-ph] .
  * García-Martín  _et al._ [2023b] D. García-Martín, M. Larocca, and M. Cerezo, Deep quantum neural networks form gaussian processes (2023b), arXiv:2305.09957 [quant-ph] .

## Supplementary Information

Figure 5: Pauli correlations under different loss functions. Histograms of
expectation values $\langle\Pi_{i}\rangle$ for the quadratic Pauli-correlation
encoding $\Pi^{(2)}$ after training under different loss functions, for 250
random graph instances of $m=108$ vertices and 5 random parameter
initializations per instance. Top left: quadratic loss function
$\mathcal{L}^{(\text{qua})}=\sum_{(i,j)\in
E}W_{ij}\,\langle\Pi_{i}\rangle\,\langle\Pi_{j}\rangle$. Top right: quadratic
loss function with regularization,
$\mathcal{L}^{(\text{qua})}+\mathcal{L}^{(\text{reg})}$. Bottom left: non-
linear loss function from Eq. (2) with regularization removed,
$\mathcal{L}-\mathcal{L}^{(\text{reg})}$. Bottom right: Complete loss function
$\mathcal{L}$ given by Eq. (2). The panels show also the average approximation
ratios $\overline{r}$ obtained in each case.

##  Appendix A Choice of loss function

Here we motivate the specific loss function chosen in Eq. (2). $\mathcal{L}$
leverages two main features: the non-linearities from the hyperbolic tangents
and the regularization term $\mathcal{L}_{\text{reg}}$, given by (5), which
forces all the correlators to have small values. In Fig. 5, we show a
comparison of the distribution of expectation values $\langle\Pi_{i}\rangle$
of Pauli string correlators at the end of a training process based on four
different loss functions. Namely, we compare Eq. (2) with similar loss
functions obtained by removing the hyperbolic tangent factors (i.e., a
quadratic function of the $\langle\Pi_{i}\rangle$) and/or the regularization
term $\mathcal{L}^{(\text{reg})}$.

We see that, for a quadratic loss function without the regularization term
(top left), the distribution of expectation values is approximately flat with
a peak around the origin and small peaks around $\pm 1$. The introduction of
the non-linear function $\tanh(\cdot)$ (bottom left) causes the expectation
values to cluster in heavy-tailed distributions around two symmetric points.
This alters the optimization landscape, thereby discouraging extremal values,
which is a particularly important feature for our encoding due to its
sensitivity to frustration constraints. We emphasize that the specific choice
of the hyperbolic tangents is not particularly important: we observed that any
sigmoid-like non-linear function leads to the same concentration phenomenon;
this is a consequence of their vanishing gradients close to the extrema. The
addition of the regularization term to each of those loss functions (top and
bottom right plots, respectively) further incentivizes the
$\langle\Pi_{i}\rangle$ to stay close to zero, reducing the tails of the
distributions and shifting their mean magnitude closer to zero. The strength
of this shift is modulated by the hyperparameter $\beta$ appearing in Eq. (5),
whose value we also fine-tune by extensive numerical exploration on random
graph instances following the same procedure of Sec. C for $\alpha$. The
result is shown in Fig. 6.

Figure 6: Tuning of $\beta$. The plots shows the average approximation ratio
versus $\beta$ over 250 randomly generated graph instances (over 5
initialization) for increasing number of qubits, at fixed depth $l=11$ . Each
plot display an order of polynomial compression: $k=2$ and $k=3$.

##  Appendix B Frustration bounds

Here we prove (for $k=2$) frustration constraints on the magnitude of Pauli
string correlators $\langle\Pi_{i}\rangle$ used in our Pauli-correlation
encodings $\Pi^{(k)}$ and verify numerically without rigorous proof the
scaling for general $k$.

Figure 7: Tuning of $\alpha$. Average approximation ratios versus $\alpha$
over 250 randomly generated graph instances (with 5 initializations each) for
increasing number of qubits, at fixed depth $l=5$ . Each plot displays an
order of polynomial compression, from $k=2$ (leftmost) to $k=6$ (rightmost).
The precise value of $\alpha$ used by our solver is fine-tuned to maximize
$r$. For $k=1$ (not shown), we observe a constant behaviour $\alpha\sim 1.5$.
The optimal values observe are thus compatible with the scaling
$\alpha=\mathcal{O}(n^{\lfloor k/2\rfloor})$, which is equivalent to the
scaling $\big{|}\langle\Pi_{i}\rangle\big{|}=\mathcal{O}(n^{-\lfloor
k/2\rfloor})$ for the Pauli correlators. For $k=2$, this agrees with the
analytical upper bound derived in App. B.

The $k=2$ Pauli-correlation encoding is characterized by two-body correlators
of the form $\langle A_{\ell_{1}}A_{\ell_{2}}\rangle$ where $A$ is one of the
traceless Pauli matrices (i.e., $X,Y,$ or $Z$) and $A_{\ell}$ indicates its
action on qubit $\ell\in[n]$. Assume that all the non-trivial correlators have
the same magnitude, i.e., $|\langle A_{\ell_{1}}A_{\ell_{2}}\rangle|=c\geq 0$
for every $\ell_{1}\neq\ell_{2}$ – say, $n_{+}$ of them are equal to $+c$ and
$n_{-}$ are equal to $-c$, with $n_{+}+n_{-}=3n(n-1)$. Our goal here is to
establish an upper bound on the strength $c$. For that, let us introduce
$A_{\text{tot}}:=\sum_{\ell\in[n]}A_{\ell}$ for each $A\in\\{X,Y,Z\\}$ and the
positive-semidefinite operator $\sum_{A\in\\{X,Y,Z\\}}A_{\text{tot}}^{2}$. It
follows that $0\leq\sum_{A\in\\{X,Y,Z\\}}\langle
A_{\text{tot}}^{2}\rangle=3n+(n_{+}-n_{-})c$. For $n_{+}>n_{-}$ this is
trivially satisfied for all $c$, so without loss of generality one can assume
$n_{-}>n_{+}$. In this case, the inequality contrains $c$ to $c\leq
3n/(n_{-}-n_{+})=1/(n-1-2n_{+}/3n)$. A sufficient condition to attain any sign
configuration corresponds to the worst case situation where all signs are
negative ($n_{+}=0$), namely $c\leq 1/(n-1)=\mathcal{O}(n^{-1})$.

In the case of arbitrary compression degree $k$, an analytic proof has been
elusive to us. However, in App. C we numerically show (see Fig. 7) that
$\absolutevalue{\langle A_{\ell_{1}}\ldots
A_{\ell_{k}}\rangle}=\order{n^{-\lfloor k/2\rfloor}}$ is a sufficient
condition.

##  Appendix C Choice of $\alpha$

We studied the behavior of the rescaling parameter $\alpha$ by looking at the
average approximation ratios achieved at the end of the optimization for 250
random graph instances of increasing size, generated in the standard way
detailed in Sec. III.3. The value of $\alpha$ was increased until a plateau in
solution quality was reached (see Fig. 7). Based on this analysis, we fine-
tune $\alpha$ for the solver at each compression rate $k$ so as to maximize
$r$. We observe that its optimal value scales with $n$ and $k$ as
$\alpha=\mathcal{O}(n^{\lfloor k/2\rfloor})$, which generalizes the analytical
bound obtained for $k=2$ in App. B.

##  Appendix D Possible classical analogues of our algorithm

There is a natural approach to classically mimic the algorithmic pipeline of
our solver. This consists of substituting the quantum circuit on $n$ qubits by
a classical generative neural network that generates a probability
distribution $P$ over, for instance, $3n$ bits, $n$ bits for each of the three
mutually-commutting sets in $\Pi^{(k)}$. With this, one can encode the binary
variables into classical correlations across $k$ bits, described by $k$-body
marginal distributions of $P$. If the network samples from $P$, one can Monte-
Carlo estimate all $m$ marginals efficiently in the same fashion as we do with
measurements on the quantum circuit. Then, one can train the network so that
its output classical correlations minimize a loss function analogous to that
in Eq. (2).
Benchmarking our scheme against the plethora of classical generative neural
models is beyond the scope of the current work. Nevertheless, our algorithm
has promising prospects coming from the fact that brickwork quantum circuits
are explicitly known [60] to produce local expectation values that cannot be
efficiently simulated classically, as well as to sample from classically
intractable output probability distributions.

##  Appendix E Training complexity

Here we provide further numerical details of the number of optimization
parameters and training epochs for compressions of degree $k=2$ and $k=3$. In
Fig. 8 we show the observed scaling with $m$ of these figures of merit over
random MaxCut instances at two different target solution qualities:
$\overline{r}=0.941$ excluding the final local search step, as in Fig. 2
(upper left and right panels); and $P_{95}(r)=1$, i.e. the $95^{\text{th}}$
percentile of $r$ was equal to one, now including the final local search step
and increasing circuits depths accordingly (lower left and right panels). The
latter was done to give an idea of the resources required to observe the exact
solution with high probability in the practical average case. In terms of gate
complexity, we observed a linear scaling in the first scenario (upper right),
and a quadratic one in the second (lower right). The number of optimization
epochs, on the other hand, was observed to scale linearly in both cases (upper
left and upper right, respectively).

Figure 8: Training complexity. Upper left: number of epochs needed to reach
optimization convergence with an average approximation ratio $\overline{r}\geq
16/17\approx 0.941$ (over 250 random MaxCut instances and 5 random
initializations per instance) without the local bit-swap search (quantum-
circuit’s output $\bm{x}$ alone) versus problem size $m$, both for quadratic
and cubic compressions. A linear scaling was observed in both cases; Upper
right: scaling of the number of parameters in the same set up of upper left
panel. Lower left: observed number of epochs under the requirement that the
$95^{\text{th}}$ percentile of $r$ was equal to one ($P_{95}(r)=1$) and
including the final local search step (full solver’s output $\bm{x}^{*}$).
Lower right: scaling of the number of parameters in the same set up of lower
left panel.

##  Appendix F Sample complexity

Here we upper-bound the minimum number of measurements needed to estimate
$\mathcal{L}(\bm{\theta})$ at some arbitrary point $\bm{\theta}$. More
precisely, given $\varepsilon$, $\delta>0$, and a vector $\bm{\theta}$ of
variational parameters, consider the problem of estimating
$\mathcal{L}(\bm{\theta})$ up to additive precision $\varepsilon$ and with
statistical confidence $1-\delta$.

For each Pauli correlator $\expectationvalue{\Pi_{i}}$, $i\in[m]$, assume
that, with confidence $1-\delta$, one has an unbiased estimator $\Pi_{i}^{*}$
with statistical error at most $\eta>0$, that is,

$\Delta\langle\Pi_{i}\rangle:=\langle\Pi_{i}\rangle-\Pi_{i}^{*}\quad\text{is
such that}\quad\big{|}\Delta\langle\Pi_{i}\rangle\big{|}\leq\eta\,.$ (6)

The corresponding error in the loss function is
$\Delta\mathcal{L}:=\mathcal{L}-\mathcal{L}^{*}$, with $\mathcal{L}^{*}$ given
by (2) computed using $\Pi_{i}^{*}$ instead of $\langle\Pi_{i}\rangle$. The
multivariate Taylor theorem ensures that there is a $\xi\in[-1,1]^{m}$ such
that

$\displaystyle\Delta\mathcal{L}=\sum_{i\in[m]}\frac{\partial\mathcal{L}}{\partial{\langle\Pi_{i}\rangle}}\bigg{|}_{\xi}\,\Delta\langle\Pi_{i}\rangle\,.$
(7)

We next restrict to MaxCut, for which the expressions take simple forms in
terms of the number of vertices $m$ and edges $|E|$, but the extension to
weighted MaxCut is straightforward. For the loss function (2), using the basic
inequalities $\tanh(x)\leq 1$ and $\frac{d}{dx}\tanh(x)=\sech^{2}(x)\leq 1$,
one can show that
$\big{|}\big{(}\frac{\partial\mathcal{L}}{\partial\langle\Pi_{i}\rangle}\big{|}_{\xi}\big{)}\big{|}\leq
2\alpha\big{[}d(i)+{2\beta\nu}/{m}\big{]}$ , where
$d(i):=\sum_{j\in[m]}\\!|W_{ij}|$ is the degree of vertex $i$. As a result,

$\displaystyle|\Delta\mathcal{L}|$ $\displaystyle\leq
2\,\eta\,\alpha\sum_{i\in[m]}\left[d(i)+\frac{2\beta\nu}{m}\right]\leq\eta\,\alpha\,(6|E|+m),$
(8)

where the first step follows from the triangle inequality together with (6),
while the second uses the identity $\sum_{i\in[m]}d(i)=2|E|$,
$\nu=(2|E|+m-1)/4$ (see Sec. III.2), and $\beta<1$.

To ensure $|\Delta\mathcal{L}|\leq\varepsilon$ we require that

$\displaystyle\eta\leq\frac{\varepsilon}{\alpha\,(6|E|+m)}\,.$ (9)

The minimum number $S$ of samples needed to achieve such precision can be
upper-bounded by standard arguments using the union bound and Hoeffding’s
inequality, which gives $S\leq(4/\eta^{2})\log(2m/\delta)$ Then, by virtue of
Eq. (9), it suffices to take

$\displaystyle S$
$\displaystyle=\left\lfloor\frac{4\alpha^{2}}{\varepsilon^{2}}\,(6|E|+m)^{2}\log(\frac{2m}{\delta})\right\rfloor.$
(10)

This is the general form of our upper bound. However, for the particular cases
$k=2$ and $k=3$, $\alpha=\mathcal{O}(n^{\lfloor
k/2\rfloor})=\mathcal{O}(m^{1/2})$ (see App. C), hence
$S=\mathcal{O}\big{(}m\,(6|E|+m)^{2}\log(\frac{2m}{\delta})/\varepsilon^{2}\big{)}$.
Moreover, in practice, it is often the case (as in all the instances in Table
1) that $|E|$ is linear in $m$, making the statistical overhead
$\tilde{\mathcal{O}}({m^{3}})$.

Graph | $m$ | Alg. | $k$ | $n$ | 2-q | Par. | Runs |  $r$ |
---|---|---|---|---|---|---|---|---|---
Mean | Max
torus | 512 | BM | - | - | - | 512 | 100 | 0.939 | 0.969
MBE | 1 | 256 | 768 | 1792 | 30 | 0.948 | 0.978
Our | 4 | 10 | 125 | 510 | 30 | 0.961 | 0.987
G14 | 800 | BM | - | - | - | 800 | 1 | 0.984 | -
Our | 5 | 11 | 200 | 811 | 10 | 0.985 | 0.991
G23 | 2000 | BM | - | - | - | 2000 | 1 | 0.989 | -
Our | 6 | 12 | 498 | 2004 | 10 | 0.992 | 0.995
G60 | 7000 | BM | - | - | - | 7000 | 1 | 0.970 | -
Our | 5 | 15 | 1750 | 7014 | 5 | 0.975 | 0.978
Table 3:  Single-shot resources and quality of solutions. For each instance we
display, for each heuristic, $m$, $k$, $n$, the 2-qubit gate counts, the
number of optimization parameters used during training, and the number of
random initializations. The last two columns report the average approximation
ratios and the best one observed at the end of training. In the case of
pm3-8-80 (torus), for which bigger statistics are available, we observed
significant improvements in the amount of required resources compared to the
single-qubit multi-basis encoding of [27], which is equivalent to our PCE with
$k=1$. Additionally, both the average and peak approximation ratios showed
improvement over both MBE and BM. For the remaining three instances, where
data from only a single initialization is available for the BM algorithm, we
noted the average performance of our method to be higher.

##  Appendix G Details on the comparison with Burer-Monteiro

Here we provide more detailed information on the results of our method on the
benchmark instances reported in Tab. 1. In their paper, Burer and Monteiro
provide an effective method to carry out their non-convex optimization (see
Algorithm-1 in Ref.[35]). After reaching a local minimum, and extensive local
search is performed (one- and two-bit swap search). Then, the obtained
parameters are perturbed, and a new minimization is carried out until
convergence is reached. If a local search on the new solution leads to a
better value of the cut, the parameters are updated, and the procedure
repeated. If, after $N$ perturbations, no better cut is found, the
optimization is halted. For all the benchmarked instances, they provide the
results obtained with different choices of number of initializations and of
$N$. On the other hand, in our method, we execute a single (quantum)
optimization followed by a final local search, which effectively places it on
equal footing to the BM algorithm with $N=0$. Given that, in the table (3) we
provide a comparison of the approximation ratios with the single-shot version
of BM.

##  Appendix H Comparison of the experimental and naive solutions

Here we compare the solutions found in our experimental demonstrations to
“naive” solutions obtained by randomly picking a graph partition and
performing a local search over it. Fig. 9 shows the cut value distributions
over 1000 naive solutions together with the cut value of our experimental
solutions. The distribution appears Gaussian, as indicated by the Gaussian fit
(pink curve). The vertical lines locate the $3\sigma$ right tail, the hardness
threshold $0.941$, and our experimental solutions for the $k=2$ and $k=3$
encodings. The results clearly indicate the non-trivial character of our
solutions, which lie beyond $3\sigma$ for the first two instances and near
$3\sigma$ for the last one. We recall that the hardness bound is not shown for
pm3-8-50 (left plot) since this is a weighted MaxCut instance.

Figure 9: Experimental versus naive solutions. Cut value distributions of
“naive” solutions obtained by performing local search over a randomly picked
bitstring (1000 random cut assignments were used) for the three benchmarking
instances used in the experiment (see Table 1). Each plot shows the $3\sigma$
right tail of the distribution based on a Gaussian fit (pink curve), the
hardness bound $0.941$, as well as the observed experimental results with
$k=2$ and $k=3$. Clearly, the experimental solution is non-trivial with
respect the naive ones.

##  Appendix I Analytical barren plateau characterization

In this section, we analytically compute the variance of our loss function for
deep circuits. That is, we compute the value to which the variance converges
as the circuit depth increases. The nonlinear hyperbolic tangent appearing in
the loss renders an exact computation involved. Hence, we begin by computing
the variance for a simplified loss function $\mathcal{L}_{\rm sim}$. This
simplified calculation will be instrumental to obtain the variance of the
actual loss function. More precisely, we first show that, under the assumption
that the circuit ensemble under random parameter initializations is a 4-design
over the special unitary group, the variance of $\mathcal{L}_{\rm sim}$ is
equal to $\frac{1}{d^{2}}\sum_{(i,j)\in
E}w_{ij}^{2}+\mathcal{O}\left(\frac{1}{d^{4}}\right)$, where $d=2^{n}$ is the
dimension of the Hilbert space. Then, we show that, under the assumption that
the circuit is fully Haar random, the variance of $\mathcal{L}$ is given by
$\frac{\alpha^{4}}{d^{2}}\sum_{(i.j)\in
E}w_{ij}^{2}+\mathcal{O}\left(\frac{\alpha^{6}}{d^{3}}\right)$.

The simplified loss function takes the form

$\mathcal{L}_{\rm sim}=\sum_{(i,j)\in
E}w_{ij}\Tr\left[U({\bm{\theta}})\,\varrho\,U({\bm{\theta}})^{\dagger}\Pi_{i}\right]\Tr\left[U({\bm{\theta}})\,\varrho\,U({\bm{\theta}})^{\dagger}\Pi_{j}\right]\,,$
(11)

with $\varrho$ a pure state. For arbitrary depths, one would be interested in
the variance of $\mathcal{L}_{\rm sim}$ over a uniform sampling of parameter
values in the interval $[0,2\pi]$, i.e. in computing ${\rm
Var}_{\bm{\theta}}\left[\mathcal{L}_{\rm sim}\right]$. However, such
computation is complicated. Instead, we will resort to representation-
theoretic techniques and compute the variance assuming that the quantum
circuit is a design over the special unitary group, that we denote
$\mathbb{SU}(d)$. In practice, if the circuit is deep enough it will always
form a design over the dynamical Lie group associated to the circuit’s
generators [41], thus justifying the utility of the computation. When the
circuit’s generators are traceless and universal (as it happens in our case),
the circuit’s dynamical Lie group is $\mathbb{SU}(d)$. Since we will work at
the special unitary group level, we will henceforth drop the explicit
dependence of the unitaries on the variational parameters ${\bm{\theta}}$.

We start by computing

$\mathbb{E}_{\mathbb{SU}(d)}\left[\mathcal{L}^{2}_{\rm
sim}\right]=\sum_{(i,j)\in E}\sum_{(k,l)\in E}w_{ij}w_{kl}\int
d\mu(U)\Tr\left[U^{\otimes 4}\varrho^{\otimes 4}(U^{\dagger})^{\otimes
4}\,\Pi_{i}\otimes\Pi_{j}\otimes\Pi_{k}\otimes\Pi_{l}\right]\,,$ (12)

where $d\mu(U)$ is the volume element from the Haar measure, and we used the
property that $\Tr[A\otimes B]=\Tr[A]\Tr[B]$. (In fact, for the simplified
loss functions it suffices that $d\mu(U)$ defines a 4-design, the fully-random
Haar measure will be needed only for the actual loss function below.) Using
standard Weingarten calculus techniques [61], we have that

$\begin{split}&\mathbb{E}_{\mathbb{SU}(d)}\left[\Tr\left[U^{\otimes
4}\varrho^{\otimes 4}(U^{\dagger})^{\otimes
4}\,\Pi_{i}\otimes\Pi_{j}\otimes\Pi_{k}\otimes\Pi_{l}\right]\right]=\\\
&\frac{1}{d^{4}}\sum_{\sigma\in S_{4}}\Tr[\varrho^{\otimes
4}P_{d}(\sigma)]\Tr[P_{d}(\sigma^{-1})\,\Pi_{i}\otimes\Pi_{j}\otimes\Pi_{k}\otimes\Pi_{l}]+\frac{1}{d^{4}}\sum_{\sigma,\pi\in
S_{4}}c_{\sigma,\pi}\Tr[\varrho^{\otimes
4}P_{d}(\sigma)]\Tr[P_{d}(\pi)\Pi_{i}\otimes\Pi_{j}\otimes\Pi_{k}\otimes\Pi_{l}]=\\\
&\frac{1}{d^{4}}\sum_{\sigma\in
S_{4}}\Tr[P_{d}(\sigma^{-1})\,\Pi_{i}\otimes\Pi_{j}\otimes\Pi_{k}\otimes\Pi_{l}]+\frac{1}{d^{4}}\sum_{\sigma,\pi\in
S_{4}}c_{\sigma,\pi}\Tr[P_{d}(\pi)\Pi_{i}\otimes\Pi_{j}\otimes\Pi_{k}\otimes\Pi_{l}]\,,\end{split}$
(13)

where $d=2^{n}$ is the dimension of the Hilbert space,
$c_{\sigma,\pi}\in\mathcal{O}(1/d)$, and $P_{d}$ is the representation of the
Symmetric group $\mathbb{S}_{t}$ that permutes the $d$-dimensional subsystems
in the $t$-fold tensor product Hilbert space, $\mathcal{H}^{\otimes t}$ (i.e.
$P_{d}(\sigma)=\sum_{i_{1},\dots,i_{t}=0}^{d-1}|i_{\sigma^{-1}(1)},\dots,i_{\sigma^{-1}(t)}\rangle\langle
i_{1},\dots,i_{t}|$, for a permutation $\sigma\in\mathbb{S}_{t}$).
Furthermore, we used that since $\varrho$ is pure, $\Tr[\varrho^{\otimes
4}P_{d}(\sigma)]=1$ $\forall\sigma\in\mathbb{S}_{4}$.

Let us now take a look at the permutations in $\mathbb{S}_{4}$. We have
$4!=24$ permutations, namely: the identity, six transpositions, three double-
transpositions, eight 3-cycles and six 4-cycles. We now note that for any
$\sigma$ containing an odd-length cycle, the term
$\Tr[P_{d}(\sigma)\Pi_{i}\otimes\Pi_{j}\otimes\Pi_{k}\otimes\Pi_{l}]$
vanishes, since $\Pi_{i},\Pi_{j},\Pi_{k},\Pi_{l}$ are all traceless. Hence, we
are left with the double transpositions and the 4-cycles. Moreover, for the
double transpositions $(ik)(jl)$ and $(il)(jk)$ we find that
$\Tr[P_{d}(\sigma)\Pi_{i}\otimes\Pi_{j}\otimes\Pi_{k}\otimes\Pi_{l}]$ is equal
to $d^{2}\delta_{\Pi_{i}\Pi_{k}}\delta_{\Pi_{j}\Pi_{l}}$ and
$d^{2}\delta_{\Pi_{i}\Pi_{l}}\delta_{\Pi_{j}\Pi_{k}}$, respectively. For the
double transposition $(ij)(kl)$,
$\Tr[P_{d}(\sigma)\Pi_{i}\otimes\Pi_{j}\otimes\Pi_{k}\otimes\Pi_{l}]=0$ since
$\Pi_{i}\neq\Pi_{j}$. Noting that $\delta_{\Pi_{i}\Pi_{k}}=\delta_{ik}$, we
thus arrive at

$\mathbb{E}_{\mathbb{SU}(d)}\left[\mathcal{L}^{2}_{\rm
sim}\right]=\frac{1}{d^{2}}\sum_{(i,j)\in E}\sum_{(k,l)\in
E}w_{ij}w_{kl}\,\delta_{ik}\delta_{jl}+\mathcal{O}\left(\frac{1}{d^{4}}\right)=\frac{1}{d^{2}}\sum_{(i.j)\in
E}w_{ij}^{2}+\mathcal{O}\left(\frac{1}{d^{4}}\right)\,,$ (14)

where the term in $\mathcal{O}\left(\frac{1}{d^{4}}\right)$ comes from the
4-cycles contributions. On the other hand, we have

$\mathbb{E}_{\mathbb{SU}(d)}\left[\mathcal{L}_{\rm sim}\right]=\sum_{(i,j)\in
E}w_{ij}\left(\frac{1}{d^{2}}\sum_{\sigma\in S_{2}}\Tr[\varrho^{\otimes
2}P_{d}(\sigma)]\Tr[P_{d}(\sigma^{-1})\,\Pi_{i}\otimes\Pi_{j}]+\frac{1}{d^{2}}\sum_{\sigma,\pi\in
S_{2}}c_{\sigma,\pi}\Tr[\varrho^{\otimes
2}P_{d}(\sigma)]\Tr[P_{d}(\pi)\Pi_{i}\otimes\Pi_{j}]\right)=0\,,$ (15)

where we used that $\Tr[\Pi_{i}\otimes\Pi_{j}]=\Tr\left[{\rm
SWAP}\,\Pi_{i}\otimes\Pi_{j}\right]=0$. We finally find that the variance is

${\rm Var}_{\mathbb{SU}(d)}\left(\mathcal{L}_{\rm
sim}\right)=\frac{1}{d^{2}}\sum_{(i,j)\in
E}w_{ij}^{2}+\mathcal{O}\left(\frac{1}{d^{4}}\right)\,.$ (16)

In particular, we notice that for the unweighted version of the MaxCut
problem, the variance is given by ${\rm
Var}_{\mathbb{U}(d)}\left(\mathcal{L}_{\rm sim}\right)=\frac{|E|}{d^{2}}$.
Equation (16) implies that if the quantum circuit is a $4$-design over the
special unitary group, then the variance of the loss function is exponentially
suppressed as $\mathcal{O}(1/2^{2n})$.

We now compute the variance for the actual loss function in the main text,

$\mathcal{L}=\sum_{(i,j)\in
E}w_{ij}\tanh\left(\alpha\Tr\left[U({\bm{\theta}})\,\varrho\,U({\bm{\theta}})^{\dagger}\Pi_{i}\right]\right)\tanh\left(\alpha\Tr\left[U({\bm{\theta}})\,\varrho\,U({\bm{\theta}})^{\dagger}\Pi_{j}\right]\right)+\mathcal{L}^{({\rm
reg})}\equiv\mathcal{L}^{({\tanh})}+\mathcal{L}^{({\rm reg})}\,.$ (17)

We proceed by using the Taylor-series expansion of the hyperbolic tangent,
namely $\tanh(x)=\sum_{s=1}^{\infty}C_{s}\,\,x^{2s-1}$, where
$C_{s}\equiv\frac{2^{2s}(2^{2s}-1)B_{2s}}{(2s)!}$ and $B_{2s}$ are the
Bernoulli numbers. We are interested in computing

$\begin{split}\mathbb{E}_{\mathbb{SU}(d)}\left[\left(\mathcal{L}^{({\tanh})}\right)^{2}\right]=\sum_{(i,j)\in
E}\sum_{(k,l)\in
E}w_{ij}w_{kl}\sum_{s_{1},s_{2},s_{3},s_{4}=1}^{\infty}C_{s_{1}}C_{s_{2}}C_{s_{3}}C_{s_{4}}\,\int
d\mu(U)&\left(\alpha\Tr\left[U\,\varrho\,U^{\dagger}\Pi_{i}\right]\right)^{2s_{1}-1}\left(\alpha\Tr\left[U\,\varrho\,U^{\dagger}\Pi_{j}\right]\right)^{2s_{2}-1}\\\
&\left(\alpha\Tr\left[U\,\varrho\,U^{\dagger}\Pi_{k}\right]\right)^{2s_{3}-1}\left(\alpha\Tr\left[U\,\varrho\,U^{\dagger}\Pi_{l}\right]\right)^{2s_{4}-1}\,.\end{split}$
(18)

Here, we will be dealing with quantities of the form

$\alpha^{t}\int d\mu(U)\Tr\left[U({\bm{\theta}})^{\otimes t}\varrho^{\otimes
t}\left(U({\bm{\theta}})^{\dagger}\right)^{\otimes t}\,\Pi_{i}^{\otimes
t_{1}}\otimes\Pi_{j}^{\otimes t_{2}}\otimes\Pi_{k}^{\otimes
t_{3}}\otimes\Pi_{l}^{\otimes t_{4}}\right]\,,$ (19)

where $t_{\gamma}=2s_{\gamma}-1$ and $t=t_{1}+t_{2}+t_{3}+t_{4}$. Using
asymptotic Weingarten calculus (see e.g. Ref. [61]), it can be shown that

$\begin{split}\int d\mu(U)\Tr\left[U^{\otimes t}\varrho^{\otimes
t}\left(U^{\dagger}\right)^{\otimes t}\,\Pi_{i}^{\otimes
t_{1}}\otimes\Pi_{j}^{\otimes t_{2}}\otimes\Pi_{k}^{\otimes
t_{3}}\otimes\Pi_{l}^{\otimes
t_{4}}\right]&=\frac{\left|\mathcal{T}_{t_{1}+t_{3}}\right|\delta_{ik}+\left|\mathcal{T}_{t_{2}+t_{4}}\right|\delta_{jl}}{d^{t/2}}\\\
&+\frac{(1-\delta_{ik})(\left|\mathcal{T}_{t_{1}}\right|+\left|\mathcal{T}_{t_{3}}\right|)+(1-\delta_{jl})(\left|\mathcal{T}_{t_{2}}\right|+\left|\mathcal{T}_{t_{4}}\right|)}{d^{t/2}}+\mathcal{O}\left(\frac{1}{d^{t/2+1}}\right)\,,\end{split}$
(20)

where e.g. $\left|\mathcal{T}_{t_{1}+t_{3}}\right|$ is the number of
permutations that consist exactly of $t_{1}+t_{3}$ disjoint transpositions,
i.e.

$\left|\mathcal{T}_{t_{1}+t_{3}}\right|=\frac{(t_{1}+t_{3})!}{2^{(t_{1}+t_{3})/2}\left(\frac{t_{1}+t_{3}}{2}\right)!}\,.$
(21)

The dimensional dependence obtained in Eq. (20) implies that, for large-enough
$d$, Eq. (18) is well-approximated by its first-order terms, i.e. the terms
where $t=4$. These first-order terms lead to a total contribution
$\mathcal{O}\left(1/d^{2}\right)$, while the rest contribute with
$\mathcal{O}\left(1/d^{3}\right)$. More precisely, we find that

$\mathbb{E}_{\mathbb{SU}(d)}\left[\left(\mathcal{L}^{({\tanh})}\right)^{2}\right]=\frac{\alpha^{4}\,C_{1}^{4}}{d^{2}}\sum_{(i.j)\in
E}w_{ij}^{2}+\mathcal{O}\left(\frac{\alpha^{6}}{d^{3}}\right)\,.$ (22)

Furthermore, we have that $\Tr[P_{d}(\sigma)\,\Pi_{i}^{\otimes
t_{1}}\otimes\Pi_{j}^{\otimes t_{2}}]=0$
$\forall\sigma\in\mathbb{S}_{t_{1}+t_{2}}$ (this is true because the
hyperbolic tangent is an odd function, which implies that $t_{1}$ and $t_{2}$
are odd). It then follows that
$\mathbb{E}_{\mathbb{U}(d)}\left[\left(\mathcal{L}^{({\tanh})}\right)\right]=0\,.$
Let us now see what is the effect on the variance of the regularization term.
This term is given by

$\mathcal{L}^{({\rm reg})}=\beta\,\nu\left[\frac{1}{m}\sum_{i\in
V}\tanh\big(\alpha\,\langle\Pi_{i}\rangle\big{missing})^{2}\right]^{2}.$ (23)

Here, it suffices to notice that
$\mathbb{E}_{\mathbb{U}(d)}\left[\left(\mathcal{L}^{({\rm
reg})}\right)^{2}\right]\in\mathcal{O}\left(\frac{1}{d^{4}}\right)\,,$
$\mathbb{E}_{\mathbb{U}(d)}\left[\mathcal{L}^{(\rm
reg)}\mathcal{L}^{({\tanh})}\right]\in\mathcal{O}\left(\frac{1}{d^{3}}\right)\,,$
and $\mathbb{E}_{\mathbb{U}(d)}\left[\mathcal{L}^{(\rm
reg)}\right]^{2}\in\mathcal{O}\left(\frac{1}{d^{4}}\right)\,,$ which follow
from applying Eq. (20). Combining these observations with Eq. (22) (and noting
that $C_{1}=1$), the final result is

${\rm
Var}_{\mathbb{SU}(d)}\left(\mathcal{L}\right)=\frac{\alpha^{4}}{d^{2}}\sum_{(i.j)\in
E}w_{ij}^{2}+\mathcal{O}\left(\frac{\alpha^{6}}{d^{3}}\right)\,.$ (24)

We remark here that since the hyperbolic tangent is expanded as an infinite
Taylor series, we require the circuit to be fully Haar random under the
parameter initialization in order for Eq. (24) to hold exactly. However, if
the circuit ensemble is already a $4$-design, we do not expect higher order
terms to differ significantly but rather to become negligible as
$n\rightarrow\infty$.