# On the approximability of random-hypergraph MAX-3-XORSAT problems with
quantum algorithms

Eliot Kapit1,3<EMAIL_ADDRESS>Brandon A. Barton2,3 Sean Feeney3 George
Grattan3 Pratik Patnaik1,3 Jacob Sagal3 Lincoln D. Carr1,2,3 Vadim
Oganesyan4,5,6 1 Department of Physics, Colorado School of Mines, 1523
Illinois St, Golden CO 80401  2 Department of Applied Mathematics and
Statistics, Colorado School of Mines, 1500 Illinois St, Golden CO 80401  3
Quantum Engineering Program, Colorado School of Mines, 1523 Illinois St,
Golden CO 80401  4 Department of Physics and Astronomy, College of Staten
Island, CUNY, Staten Island, NY 10314, USA  5 Physics program and Initiative
for the Theoretical Sciences, The Graduate Center, CUNY, New York, NY 10016,
USA  6 Center for Computational Quantum Physics, Flatiron Institute, 162 5th
Avenue, New York, NY 10010, USA

###### Abstract

A huge range of important problems in computer science–including task
optimization, formal logic, encryption, and machine learning–can be solved by
finding the sequence of binary variables that optimizes a cost function
defined by a series of few-variable constraint relationships. These problems
define the complexity class NP, and are in the worst case, and often the
typical case, exponentially hard in the number of variables for all known
methods. This hardness applies both to exact and approximate optimization,
e.g. finding configurations with a value within a defined fraction of the
global optimum. Fundamentally, the lack of any guided local minimum escape
method ensures the hardness of both exact and approximate optimization
classically, but the intuitive mechanism for approximation hardness in quantum
algorithms based on Hamiltonian time evolution is not well understood. In this
work, using the prototypically hard MAX-3-XORSAT problem class, we explore
this question. We conclude that the mechanisms for quantum exact and
approximation hardness are fundamentally distinct. We qualitatively identify
why traditional methods such as high depth quantum adiabatic optimization
algorithms are not reliably good approximation algorithms. We propose a new
spectral folding optimization method that does not suffer from these issues
and study it analytically and numerically. We consider random rank-3
hypergraphs including extremal planted solution instances, where the ground
state satisfies an anomalously high fraction of constraints compared to truly
random problems. We show that, if we define the energy to be
$E=N_{\mathrm{unsat}}-N_{\mathrm{sat}}$, then spectrally folded quantum
optimization will return states with energy $E\leq AE_{\mathrm{GS}}$ (where
$E_{\rm GS}$ is the ground state energy) in polynomial time, where
conservatively, $A\simeq 0.6$. We thoroughly benchmark variations of
spectrally folded quantum optimization for random classically approximation-
hard (planted solution) instances in simulation, and find performance
consistent with this prediction. We do not claim that this approximation
guarantee holds for all possible hypergraphs, though our algorithm’s mechanism
can likely generalize widely. These results suggest that quantum computers are
more powerful for approximate optimization than had been previously assumed.

###### Contents

  1. I Introduction
  2. II Hardness mechanisms, problem definitions and previous approaches
    1. II.1 Inability of previous methods to find or even approach the ground state at large $N$
    2. II.2 Mechanisms for quantum solution hardness: exponentially small gaps and transverse field chaos
    3. II.3 Approximation hardness and conventional AQC/QAOA
    4. II.4 The MAX-3-XORSAT problem and approximation-hard instance construction
    5. II.5 Absence of angle fine tuning in this work
  3. III Spectrally folded quantum optimization algorithms
  4. IV Analytical performance predictions
    1. IV.1 Preliminaries
    2. IV.2 Paramagnet to spin glass transition scaling for MAX-3-XORSAT
    3. IV.3 Scaling of inter-valley tunneling in a $p=3$ quantum spin glass
    4. IV.4 Achievable approximation ratio with spectrally folded trial minimum annealing
    5. IV.5 Further Comments and Caveats
  5. V Numerical tests of approximation hardness for traditional methods
    1. V.1 Setup and summary of results
    2. V.2 Performance of greedy classical algorithms
    3. V.3 Performance of high-depth QAOA for this problem
  6. VI Numerical tests of spectral folding variations
    1. VI.1 Spectrally folded adiabatic interpolation performance
    2. VI.2 Trial minimum annealing performance
    3. VI.3 Discussion of our spectrally folded optimization results
    4. VI.4 Extensions of these results
  7. VII Conclusions and Outlook
  8. VIII Acknowledgments
  9. A Tunneling between two $p$-spin wells
  10. B Band structure considerations for linear spectral folding

##  I Introduction

Combinatorial optimization of constraint satisfaction problems (CSPs) is an
enormously important–and often, enormously difficult–area of modern computer
science Korte et al. (2011). Specifically, a huge array of problems in
optimization, cybersecurity, machine learning, and more amount to finding low-
energy configurations of large collections of $M$ few-body constraints over
$N$ binary variables, see Fig. 1. Generically, the energy landscape of these
cost functions is extremely rough (see Fig. 3) with exponentially many local
minima, making it very difficult to find the true ground state or even a
sufficiently low energy configuration.

From the point of view of statistical physics these cost functions are often
equivalent to the Hamiltonian of a _disordered spin glass_ , and the core
hardness mechanism comes from the inability of the system to efficiently
escape local minima by flipping small numbers of spins at each step Garey et
al. (1974); Ullman (1975); Grest et al. (1986); Crescenzi et al. (1995);
Hochba (1997); Monasson (1998); Woeginger (2003); Bapst et al. (2013); Hillar
and Lim (2013); Venturelli et al. (2015); Bellitti et al. (2021); Jones et al.
(2022). In the hardest problems–and often, even, in the typical case–the time
to find the solution grows exponentially for all known classical and quantum
methods. Remarkably, these problems are not only hard to solve, but also hard
to approximate, where approximation is defined as finding any configuration
within a defined fraction of the global optimum Arora and Barak (2009). And
while a host of clever algorithms have been proposed over the years to attack
these problems, supplemented by rapid growth in computing power since the
invention of integrated circuits, the exponential worst-case difficulty
scaling remains, and is believed (in the as yet unproven statement that ${\rm
P\neq NP}$) by most computer scientists to be a fundamental and insurmountable
fact of our reality, at least for classical computers.

To make further progress, a host of heuristic quantum algorithms have been
developed, such as analog quantum annealing Finnila et al. (1994); Kadowaki
and Nishimori (1998); Das and Chakrabarti (2008); Johnson et al. (2011); Boixo
et al. (2014); Hauke et al. (2019); King et al. (2023), and its closed system
and digital cousins, adiabatic quantum computing (AQC) Farhi et al. (2000);
Albash and Lidar (2017) and quantum approximate optimization algorithms
(QAOAs) Farhi et al. (2014). In all of these methods the spin glass problem
Hamiltonian (diagonal in the computational $z$ basis) is combined with a
transverse field term (typically, a uniform field along $x$), which allows the
system to escape from local minima by multiqubit tunneling, a collective
process where large clusters of qubits all change configuration simultaneously
to tunnel from one minimum to another Grattan et al. (2023). In some
artificial problem instances this process has been shown to produce
exponential speedups compared to classical simulated annealing, with
polynomial speedups observed in experiment for short-ranged graphs Albash and
Lidar (2018); King et al. (2018); Ebadi et al. (2022); King et al. (2023).
However, generally applicable beyond-quadratic speedups for NP optimization
problems have not been realized, and given the present noisy state of quantum
hardware and the projected severe overhead of error correction Babbush et al.
(2021), this degree of speedup is insufficient for practical advantage over
classical machines.

In this work, we approach the question of approximation hardness with the goal
of identifying both core intuitive mechanisms ensuring it, and opportunities
for exponential quantum advantage. In particular, for MAX-3-XORSAT (a
particularly difficult CSP class, defined in detail in section II.4), we argue
the following:

  * •

The hardness mechanism(s) for directly finding the ground state of a given
problem Hamiltonian $H_{P}$ with heuristic quantum methods can be readily
identified, and likely cannot be circumvented in the worst cases.

  * •

Unlike classical algorithms based on local updates, the mechanisms which
ensure that it is hard to find the ground state do not generalize to ensure
approximation hardness for quantum algorithms.

  * •

However, there are good reasons to believe that traditional quantum approaches
such as AQC or QAOA are not effective approximators (e.g. reliably returning
low energy states in polynomial time) in the worst case. We present relatively
large scale numerical simulations that support this expectation.

  * •

Understanding why this is the case suggests a novel _spectral folding_ quantum
strategy, that transforms $H_{P}$ and then solves the transformed Hamiltonian
through more traditional means. For some variations the performance of
spectrally folded quantum optimization can be predicted analytically and
promises an efficient approximation guarantee for an extremely large fraction
of instances, well into the classically hard regime. For random hypergraphs in
this regime, spectrally folded quantum optimization provides an exponential
speedup for returning low energy states.

We back up all of these claims with extensive theoretical analysis and
numerical tests of all core predictions, to the largest feasible system sizes
for simulating quantum algorithms, and to the largest sizes needed to ensure
asymptotic scaling has been reached for classical algorithms. In doing so, we
define instance construction rules that ensure classical approximation
hardness, at least for algorithms based on local updates, and numerically
establish that high depth QAOA does not exhibit meaningful quantum advantage
in finding either exact or approximate solutions to these hard instances. We
analytically and numerically show that spectrally folded quantum optimization
can efficiently approximate these problems, in practice in a linearly growing
number of cost function evaluations.

This paper is structured as follows. In section II, we first provide an
overview of classical and quantum approximation hardness and the MAX-3-XORSAT
problem, and define the construction rules for the problems studied in this
work. In section III, we define spectrally folded quantum optimization and
propose two variations of our algorithm. Then, in section IV, we develop an
extensive-order resummed quantum perturbation theory capable of predicting the
performance of these algorithms, and establish an approximation guarantee for
random hypergraphs. To verify all of these claims, in section V, we present
extensive numerical tests and simulations for a variety of algorithms and
problem parametrizations. A summary of the key simulation results is presented
in table 1. We finally offer concluding remarks, and include additional
technical details in the appendices.

##  II Hardness mechanisms, problem definitions and previous approaches

###  II.1 Inability of previous methods to find or even approach the ground
state at large $N$

To motivate our novel methods, it is important to first review the qualitative
reasons classical and established quantum approaches are unable to efficiently
solve or approximate these problems. The classical failure mechanism is
straightforward: hard problems display a high density of poor quality local
minima, e.g. high energy as compared to the true ground state. Once a local
minimum is reached, as no general mechanism for guided local minimum escape
exists it is impossible to know in general how close one is to a ground state,
either in energy or Hamming distance (number of bit flips separating two
states), at least unless the found minimum happens to satisfy a large fraction
of constraints.

We emphasize that we are concerned here with an approximation guarantee, not
just a method that works well in practice. For completely random problems, one
can often predict the average ground state energy using statistical physics
arguments, such as the famous Parisi solution to the Sherrington-Kirkpatrick
model Parisi (1979). If one is able to find configurations close to this
energy in a given instance, that is not sufficient to rule out the existence
of some other, much deeper minimum far away in configuration space, even if
randomly drawing problems with such deep minima is exponentially unlikely. In
other words, there is no efficient classical algorithm to know if a given
instance is extremal in this way unless P=NP Arora and Barak (2009).

Since there are in many cases exponentially many more poor quality minima than
low-lying ones, the basin of attraction of the true ground state is generally
an exponentially small fraction of total configuration space and thus very
hard to find, a phenomenon that has been referred to as an “entropic barrier”
to problem solving Bellitti et al. (2021). In particular, the authors of
Bellitti et al. (2021) showed that for the 3-XORSAT problem they considered,
once a local minimum was found it was more efficient to simply restart the
algorithm from a random state instead of attempting to climb out of the found
minimum through penalized local operations, as in simulated annealing or
parallel tempering Earl and Deem (2005). We expect this may be a generic
feature in some of the hardest CSP classes. And interestingly, these arguments
apply equally well to approximation hardness, not just finding the optimal
solution. While for a given class and system size approximation is nominally
an easier problem, given that there are many more valid approximate solutions,
both tasks scale exponentially in the worst case, for fundamentally the same
reason. We now turn to quantum algorithms, for which the situation is
considerably more complex.

###  II.2 Mechanisms for quantum solution hardness: exponentially small gaps
and transverse field chaos

We consider a broad class of heuristic quantum algorithms, which derive from
quantum annealing, AQC and QAOA. These algorithms can be supplemented with
additional gate model techniques, such as amplitude amplification Shaydulin et
al. (2023), which improve performance but do not circumvent exponential
runtimes. In these algorithms the system is initialized in the ground state of
a trivial Hamiltonian, which is then slowly interpolated into $H_{P}$; the
system is then measured. Other variations based on energy matching Baldwin and
Laumann (2018); Smelyanskiy et al. (2018); Kechedzhi et al. (2018);
Smelyanskiy et al. (2019) initialize a known or planted low energy state of
$H_{P}$ and then use collective quantum tunneling to try to find other low
energy states.

These algorithms are fundamentally distinct from classical approaches in two
ways. The first is the presumed mechanism for quantum advantage: multiqubit
quantum tunneling (or collective quantum phase transitions in general), where
local minima can be efficiently escaped through many-body quantum effects that
have no classical analog. Second and more subtly, where classical local update
algorithms all start from a random high energy state and attempt to cool to
low energy states, quantum algorithms start below the energy of the problem
ground state. They then attempt to tunnel into the global minimum and into
other low energy states as the energy of the initial state crosses the target
state as the total Hamiltonian changes.

Unfortunately, when compared to other applications such as quantum simulation
or factoring large numbers with Shor’s algorithm, the realistic performance
advantage of these algorithms is generally much more modest. In the worst
case–and for many problem classes, the typical case–the macroscopic quantum
tunneling rate into the ground state decreases exponentially in $N$. This is a
fairly generic expectation, as in many cases, including the MAX-3-XORSAT
problem discussed below, the tunneling rate at the crossing point can be
computed using $N$th order perturbation theory and the convergence of such a
method implies exponential decay. For some specific problems this can be
circumvented by introducing additional terms into the evolution, such as the
well-known result of ramping transverse field terms down one-by-one in mean-
field $p$-spin ferromagnets Susa et al. (2018), though generalizations of this
and other methods Hauke et al. (2019); Graß (2019); Sels and Polkovnikov
(2017); Bao et al. (2018); Yang et al. (2017); Čepaitė et al. (2023) to
realistic disordered problems are not expected to show similar advantages at
large scales.

Further, for a given class, even if one can somehow ensure that the
paramagnet-to-spin-glass transition decays polynomially (as it does in the
Sherrington-Kirkpatrick problem Montanari (2021); Farhi et al. (2022) and,
likely, for MAXCUT Boulebnane and Montanaro (2021); Basso et al. (2021)), that
is not sufficient to ensure that the solution can be found in polynomial time.
This is because of a phenomenon known as transverse field chaos (TFC)
Altshuler et al. (2010); Knysh (2016), where energetic corrections from the
transverse field can change the energy hierarchy of classical minima in the
quantum spin glass phase, potentially pushing local minima below the energy of
the true ground state of $H_{P}$ (when all transverse terms are turned off).
Consequently, optimization methods will steer the system toward these false
ground states first, with additional phase transitions that occur as the
transverse field is further weakened. As these transitions occur at weak field
values, from the analysis in section IV they are generically exponentially
slow, and indeed, engineering this effect intentionally is an elegant way to
craft hard benchmark problems for quantum algorithms Tang and Kapit (2021).
This effect can be avoided by restricting the algorithm to unstructured driver
Hamiltonians Farhi et al. (2010) or very weak transverse fields, but in either
case performance is very poor. The combination of exponentially small gaps and
TFC make it extremely unlikely that any quantum algorithm of this type can
reliably and directly find the ground state of NP-complete problems.

###  II.3 Approximation hardness and conventional AQC/QAOA

Approximation hardness, however, is another story. For many NP-hard problems
guaranteeing an approximation better than random guessing by a constant
fraction is also NP-hard Sahni and Gonzalez (1976); Håstad (2001); Arora and
Barak (2009). As TFC involves crossings between states that were close in
energy to begin with, it cannot by itself lead to quantum approximation
hardness at this level. So for example, if the random state energy of a given
class is chosen to be zero and the ground state $-N$, a sufficiently general
convention, then any algorithm which could return states in polynomial time
with energy $\leq-cN$ for constant $c>0$ for all instances would promise a
potentially exponential speedup. Thus, a hypothetical algorithm which always
returned states with energy $-N/3$ or below in polynomial time assuming
spectral continuity would still promise an exponential speedup for the hardest
problem classes even if TFC reduced that guarantee to $-N/4$. In other words,
while TFC can prove ruinous for finding an exact solution when the problem
exhibits a clustering phase Mézard et al. (2002, 2005); Hartmann and Weigt
(2006); Krzakala and Kurchan (2007); Altarelli et al. (2008); Ibrahimi et al.
(2012); Jones et al. (2022) and there are exponentially many states very close
to the ground state in energy but well separated from it in Hamming distance,
it is not going to push zero energy states into competition with
$O\left(1\right)$ fractions of the ground state energy.111It’s important to
note that for some problems where the classical approximation hardness
threshold is relatively close to 1, such as MAXCUT, TFC might well become a
serious obstacle to achieving quantum advantage for approximation.

Exponentially small gaps are a more serious problem, but those too cannot so
easily be assumed to ensure approximation hardness. This is because the
empirical “difficulty exponents” of phase transitions in low-order CSPs are
often quite small; for random hypergraph MAX-3-XORSAT instances, the minimum
gap at the paramagnet-spin glass transition scales as approximately
$\Omega_{0}\left(N\right)\sim 2^{-cN}$ with $c\simeq 0.14$ (see section IV.2).
And while the fraction of states $p_{E}\left(N\right)$ below the approximation
threshold is typically exponentially small, the total number of such states is
exponentially large. If we could assume that the mixing rate between the
paramagnetic initial state and the dressed excited states of $H_{P}$ is
equivalent to that of the ground state, then the naive product of
$\Omega_{0}^{2}\left(N\right)\times 2^{N}\times p_{E}\left(N\right)$ is often
exponentially large as well, which would ensure fast approximation. But this
naive analysis is inaccurate, because as we argue momentarily, mixing matrix
elements with excited states decay more quickly, and the decay exponents grow
with excitation energy. The exponential number of target states in
approximation problems allow an algorithm suffering from both TFC and
exponentially decaying matrix elements to still guarantee an efficient
approximation, provided, at least, that those decay exponents are not too
small. In some sense, our novel spectrally folded quantum algorithm achieves
fast approximation of classically hard instances by ensuring that the
exponential decay rate of tunneling into the problem’s excited states matches
or exceeds that of the ground state in more traditional methods.

Given all this, should we expect that AQC or QAOA will efficiently approximate
NP-hard problems? Formally, it has not been proven that these algorithms are
not efficient approximators when the circuit depth is allowed to grow as a
low-order polynomial in system size, though there are good reasons to be
doubtful, supported by a number of recent works Basso et al. (2022);
Boulebnane and Montanaro (2022); Anshu and Metger (2023); Benchasattabuse et
al. (2023). Let us consider, qualitatively, how AQC/QAOA solves such a
problem, assuming the limit of quasi-continuous time; we present a quantiative
analysis below in Sec. IV.2. The system is initialized in the paramagnetic
ground state of a uniform transverse field Hamiltonian $H_{D}$, and the system
evolves in time interpolating between $H_{D}$ and $H_{P}$ by lowering the
coefficient of one and raising the coefficient of the other, raising the
energy of the paramagnetic state until it crosses the problem ground state
from below. We again emphasize the fundamental distinction of crossing from
below, rather than cooling from above, in quantum and classical algorithms.
The minimum gap at the transition is expected to decay exponentially in $N$
and is approximately given by the overlap of the dressed problem ground state
$\left|G_{D}\right>$ with the paramagnet state $\left|S\right>$; we perform
this calculation in Sec. IV.2. If the dressings are weak enough that
$\left|G_{D}\right>=\left|G\right>$ (the Grover limit Roland and Cerf (2002)),
then $\left<S|G_{D}\right>=2^{-N/2}$; however for random MAX-3-XORSAT problems
the perturbative corrections spread $\left|G_{D}\right>$ over a more
significant (if still exponentially small) fraction of Hilbert space and
reduce the decay exponent to around a quarter of that in the Grover case.

Assuming that we evolve time too quickly (e.g. not exponentially long) and
miss the primary phase transition, we can ask how efficiently $\left|S\right>$
will mix with the problem’s excited states, as these rates ultimately
determine the algorithm’s efficacy as an approximator. Calculating them
directly is very difficult, but we can qualitatively predict that they should
be much smaller for two reasons. First, these crossings occur as the
transverse field strength $\kappa$ is reduced toward zero, and since the
dressings that reduce the decay exponent all scale with extensive powers of
$\kappa$, they will be much reduced by any reduction in $\kappa$ itself.
Second, when considering a dressed excited state $\left|E_{D}\right>$, the
perturbative dressings that come from mixing with states with lower energy now
have opposite sign and destructively interfere with other corrections in the
overlap with $\left|S\right>$, in contrast to $\left|G_{D}\right>$ where all
higher order terms are positive definite. Both of these effects can, and do,
considerably worsen the mixing rates with excited states and make AQC/QAOA a
poor approximator in the worst cases, and often in practice.

###  II.4 The MAX-3-XORSAT problem and approximation-hard instance
construction

Given the severe overhead of fault tolerance Babbush et al. (2021), quantum
hardware is expected to exhibit enormous prefactor disadvantages as compared
to parallel silicon, particularly when the comparison is made to hardware with
equivalent financial value (e.g. millions of USD). Quantum algorithms thus
have the most promise when the problem is hard or outright impossible for
classical machines. NP-hard constraint satisfaction problems are no exception,
so when benchmarking a proposed quantum algorithm it is important to ensure
that the problem classes we consider are sufficiently hard for classical
machines, and in the present NISQ era, exhibit their exponential difficulty at
small enough $N$ that numerical simulations of quantum algorithms can
demonstrate meaningful improvements. MAX-3-XORSAT problems are ideal for
benchmarking quantum algorithms because their exponential difficulty scaling
is obvious at small $N$ for both classical and prior quantum approaches, in
contrast to other problems where the asymptotic exponential scaling often does
not set in until system sizes that are prohibitively large for simulation.

Figure 1: Graphical representation of a random 3-uniform hypergraph
$G_{\mathcal{H}}=(\mathcal{V},\mathcal{E})$ used in the MAX-3-XORSAT problem.
The set of vertices $\mathcal{V}$ labelled $\\{v_{1},\dots,v_{6}\\}$ are
connected to a hyperedge (i.e., constraint) from the set $\mathcal{E}$ on the
right labeled $\\{c_{1},\dots,c_{8}\\}$ if the vertex is in constraint. Each
constraint necessarily contains a random set of three unique vertices. The
legend specifies the value of the variable $V_{ijk}=\pm 1$ where the grey and
white filled $c_{j}$ boxes denote a $(+1)$, and $(-1)$ valued constraint
respectively. This example has a constraint density of $N_{C}=\frac{4}{3}N$
for 6 vertices and 8 constraints, although we use other values of $N_{C}$
throughout this work.

To prototype our quantum algorithms, we thus consider MAX-3-XORSAT Arora and
Barak (2009), a well-studied problem that consists in the Pauli basis of a
hypergraph of $N_{C}$ three-body constraint terms, as sketched in figure 1:

$\displaystyle H_{P}=-\sum_{ijk}^{N_{C}}V_{ijk}Z_{i}Z_{j}Z_{k},\;\;V_{ijk}=\pm
1.$ (1)

A constraint is said to be satisfied if, for a given bitstring,
$V_{ijk}Z_{i}Z_{j}Z_{k}=1$, and unsatisfied otherwise. This is called a
hypergraph because $V_{ijk}$ has three indices rather than the usual two found
in graph theory. Thanks to the linearity of the problem, one can use Gaussian
elimination to check if a solution exists that satisfies all the constraints
in $O\left(N^{3}\right)$ time, but if the problem is not fully satisfiable,
finding the lowest energy state(s) is NP hard. A random state satisfies half
the constraints on average and is thus energy zero. Further, it was shown by
Håstad Håstad (2001) that if the true ground state satisfies a fraction
$\left(1-\epsilon\right)$ of the constraints, then finding any configuration
that satisfies more than $\left(1/2+\epsilon\right)$ of them is also NP-hard.
The hardest instances are thus those with small but finite $\epsilon$, e.g.
almost satisfiable problems, as both finding the true ground state, and even
finding an approximate solution, is exponentially difficult. Note that if the
problem graph is sparse (e.g. $N_{C}/N$ is on the order of 1) finding
approximate solutions can still be easy, since one can randomly select a
fraction of the constraints $<cN$ (for some $O\left(1\right)$ constant $c$),
and solve that new, much easier problem; solutions to this sub-problem will
satisfy half the remaining constraints, on average.

To ensure that we are studying problems that are both hard to solve and hard
to approximate for all known methods, we consider a family of instances we
call planted partial solution problems (PPSPs). To construct a PPSP, we choose
a small unsatisfied fraction $\epsilon$ and pick a random hypergraph of
$N_{C}\gg N$ unique triplets; we use $\epsilon=0.1$ in all simulations here.
We then pick a random bitstring $G$ and randomly select
$\left(1-\epsilon\right)N_{C}$ of the constraints to be satisfied in $G$, by
picking the sign of $V_{ijk}$ appropriately, with the rest unsatisfied. If
$\epsilon$ is small and $N_{C}/N\gg 1$, $G$ will be the problem ground state
with very high probability, as the SAT/UNSAT transition for this problem is at
$N_{C}/N\sim 0.92$ Dubois and Mandler (2002) and at densities much higher than
the SAT/UNSAT threshold ground states for random graphs satisfy
$N_{C}/2+O\left(\sqrt{N_{C}}\right)$ constraints. This property also ensures
that $G$ is a unique ground state with high probability at large $N$. When we
refer to random hypergraph problems throughout this work, we refer to this
construction rule: a random, potentially fairly dense, hypergraph where one
can optionally randomly chose an anomalously large fraction of constraints to
be satisfied by matching the signs to a randomly chosen ground state
bitstring.

Our PPSP construction is necessary because truly approximation hard
problems–where the practical polynomial time approximation difficulty
approaches the random guessing limit of the complexity class separation–are
rare in the space of all possible instances. Sufficiently sparse problems are
approximation-easy, and for denser random problems one can always find strings
that satisfy $N_{C}/2+O\left(\sqrt{N_{C}}\right)$ constraints in polynomial
time Håstad and Venkatesh (2002), with a smaller prefactor in front of the
$\sqrt{N_{C}}$ than the prefactor in the average satisfied in the ground
state. We formulated our PPSP construction to ensure our algorithm was being
benchmarked on instances with a plausible claim to true classical
approximation hardness. We note that commonly studied 3-regular problems
Bellitti et al. (2021); Kowalsky et al. (2021) do not display strong
approximation hardness, as they are sparse, and can be solved efficiently if
satisfiable. And intriguingly, we show that, by some metrics, the performance
of our novel algorithm progressively improves with increasing $N_{C}/N$ in
this regime, at constant $\epsilon$.

###  II.5 Absence of angle fine tuning in this work

Before presenting our algorithm and main results, we want to make one last
point about what we mean by QAOA when talking about quantum algorithms.
Specifically, following the original proposal Farhi et al. (2014), most
studies of QAOA allow the individual angles governing the magnitudes of the
driver and problem Hamiltonians, $H_{D}$ and $H_{P}$, to vary independently at
each timestep. For finding ground states with high depth circuits, this
technique is analogous to the schedule fine tuning concept that dates back to
the adiabatic formulation of Grover’s algorithm Roland and Cerf (2002), and
yields a quadratic speedup. And for many problems, these algorithms exhibit
_concentration_ , where the set of angles that is optimal for one randomly
generated instance is close to optimal for another with high probability Farhi
et al. (2022); Anshu and Metger (2023). More general formulations such as
ADAPT-QAOA Zhu et al. (2022) can potentially offer more significant speedups,
albeit with the challenge of a much larger search space for optimizing control
parameters.

We avoid this approach in our work, and instead focus only on simple
heuristics where the schedule is determined from smooth functions based on
intuitive guessing and a small amount of trial and error. We do this for two
reasons. First, as argued in Kapit and Oganesyan (2021), quadratic speedups
from schedule fine tuning are extremely fragile and unlikely to be viable at
large $N$ except in very narrow circumstances. Second, as the truly
approximation hard instances of MAX-3-XORSAT (and we suspect, many other CSPs)
are extremal in the space of random problems, it is less obvious that
concentration arguments would apply to the cases we consider. So even if, for
example, these instances exhibit concentration and some set of angles is near-
optimal for a specific $\epsilon$ and $N_{C}/N$ scaling, we do not think it
can be easily assumed that those angles would generalize to other extremal
parametrizations.

The fundamental challenge of approximation hardness is to find, or rule out
the existence of, one or more very deep minima in an exponentially large
search space. Consequently, optimizing angles based on average energies
returned is not expected to improve algorithm performance. If the algorithm
finds a state in the basin of attraction of the deep minimum, we can likely
halt as a successful approximation has been found, but if it does not, sets of
angles that return lower energies for random, uncorrelated shallower minima
are not expected to improve the probability of reaching states near the deep
minimum. And again, even if a given subset of extremal problems does exhibit
concentration for angle optimization, those angles may not generalize to any
other extremal part of parameter space.

We therefore restrict all of our simulations in this work to linearly
increasing total runtime with $N$, and smooth, simple functions to control the
relative magnitudes of $H_{D}$ and $H_{P}$. Unsurprisingly, the best scaling
and prefactor choices differ somewhat for QAOA and variations of spectrally
folded optimization, and the results here are not optimal for any specific
algorithm variation or PPSP subclass, but rather represent a decent choice,
found by intuition and trial and error, for a broad set of parameters.

##  III Spectrally folded quantum optimization algorithms

The core idea of spectral folding (and related, more general spectral
deformations) is to modify how the Hamiltonian is applied to the quantum state
through the introduction of a filter function. Specifically, the algorithms we
consider solve problems through simulating the time evolution of a quantum
state, as $\left|\psi\right>\to e^{-iH\left(t\right)dt}\left|\psi\right>$,
with the exponentiated Hamiltonian discretized as a series of layers e.g.
$e^{iaH_{D}}e^{ibH_{P}}$. The driver Hamiltonian, and any other additional
Hamiltonian terms, are not changed by spectral folding, so we will ignore them
for now and focus on the problem Hamiltonian itself. Specifically, we write
$\left|\psi\right>$ in the computational basis as
$\left|\psi\right>=\sum_{m=0}^{2^{N}-1}c_{m}\left|m\right>$, where $m$ is the
decimal integer representation of a given bitstring. Then, for (arbitrary)
control angle $\gamma$:

$\displaystyle e^{i\gamma
H_{P}}\left|\psi\right>=\sum_{m=0}^{2^{N}-1}e^{i\gamma
E\left(m\right)}c_{m}\left|m\right>,$ (2) $\displaystyle
E\left(m\right)=\left<m\right|H_{P}\left|m\right>=-\sum_{ijk}^{N_{C}}\left<m\right|V_{ijk}Z_{i}Z_{j}Z_{k}\left|m\right>.$
(3)

In other words, the phase of each component state advances proportionally to
its energy under the problem Hamiltonian, and that energy is computed at each
step by applying a sequence of gates to implement each constraint.

In spectral folding, the phase of each component instead advances proportional
to an arbitrary function $f$ of the diagonal $H_{P}$,

$\displaystyle\left|\psi\right>\to e^{i\gamma
f\left(H_{P}\right)}\left|\psi\right>=\sum_{m=0}^{2^{N}-1}e^{i\gamma
f\left(E\left(m\right)\right)}c_{m}\left|m\right>.$ (4)

This can be accomplished by introducing a register of auxiliary qubits,
applying a gate sequence that maps the sum of the constraint terms to a
fraction of that register to store $E$, using a second fraction of that
register to compute $f\left(E\right)$, applying a sequence of controlled-phase
gates to advance the phase by $f\left(E\right)$, and then uncomputing the
previous steps to return the register to its initial state. The entire process
is sketched in figure 2. Provided $f$ is a relatively simple function, this
adds a multiplicative overhead which is polylogarithmic in $N$, since
$E\left(m\right)$ is bounded by a polynomial in $N$ and each arithmetic
operation takes $O\left(\log{N}\right)$ steps.222We note that for any choice
of $f$ more complex than multiplying $E$ by a constant (something that does
not require auxiliary qubits to begin with), any spatial locality the graph
might have is lost in this step, since $E\left(m\right)$ is a global quantity
which we are deforming with $f$.

Figure 2: A schematic for implementing spectrally deformed time evolution
$U=\exp(if(H_{P})dt)$. A set of gates maps the value of each of the $N_{C}$
constraint to adding or subtracting 1 to a register of ancilla qubits, using
binary signed adder circuits controlled by the value of
$V_{ijk}Z_{i}Z_{j}Z_{k}$. A sequence of additional gates computes
$f\left(E\right)$ from $E$ (likely using more ancilla qubits), and then a set
of local $Z$ rotations is applied to the register storing $f\left(E\right)$ to
advance the phase of each component of $\left|\psi\right>$ proportionally. The
computation of $f\left(E\right)$ and entangling with “constraint-controlled”
gates are then uncomputed, returning the ancillas to their initial state and
disentangling them from the $N$ primary qubits over which the problem is
defined. The net result of this entire process is to enact the operation in
Eq. (4), evolving time under an arbitrary function of the diagonal Hamiltonian
(transverse field layers and other operations on the primary qubits are not
shown). For relatively simple functions, the net overhead of this entire
process (compared to enacting $\exp\left(iH_{P}dt\right)$ directly) is
polylogarithmic in $N$.

We define _spectral warping_ as any $f$ which applies a nonlinear rescaling of
$E$, such as $f\left(E\right)=cE^{2}$, and _spectral folding_ as a choice of
$f$ that mirrors $E$ about a specific value, e.g.
$f\left(E\right)=\left|E-E_{t}\right|$. The core idea is sketched in figure 3,
and these two methods can of course be composed. Incorporating this operation
enormously expands the space of quantum optimization algorithms we can define;
in this work we focus on two choices, _linear_ and _quadratic_ spectral
folding. Specifically, if we choose our problem normalization via including a
multiplicative constant333Formally, this choice assumes that we know the
fraction of the $N_{C}$ constraints which are satisfied in the ground state,
something that we cannot know in advance of running our quantum optimization
algorithm! However, we can simply repeatedly run the algorithm with different
normalization choices to guess its value, a prefactor overhead of at most
$O\left(N_{C}\right)$. so that the ground state energy is
$E_{\mathrm{GS}}=-N$, and let $E_{t}=AN$ for a constant $A$, then we define
linear and quadratic spectral folding as

$\displaystyle
f_{\mathrm{lin}}\left(E\right)=\frac{\left|E+AN\right|}{A},\;\;\;f_{\mathrm{quad}}\left(E\right)=\frac{\left(E+AN\right)^{2}}{A^{2}N}$
(5)

Here, $A<1$ defines the approximation target. And critically, it is defined
using the conventions that random states have energy zero, so returning a
state with energy $AE_{\mathrm{GS}}$ approximates, by a factor of $A$, the
degree to which the true ground state itself improves on random states. These
normalization choices ensure that the energy difference between the new ground
states of the folded problem, and random states, is $N$ as in the original
renormalized problem. Making this choice simplifies the analysis
significantly.

A good choice of $A$ is important for spectrally folded quantum optimization
to succeed; if $A$ is chosen to be too close to 1, then we risk failing to
well-approximate $H_{P}$ due to the interference effects mentioned in section
II.3. A choice of $A$ which is too small will return a suboptimal
approximation ratio, and if $A$ is too close to zero, cause instabilities from
having a poorly defined problem to solve. Fortunately, for random hypergraph
MAX-3-XORSAT instances–and here random refers to the graph itself and not,
critically, on how many constraints are satisfied in the ground state–we can
predict the threshold $A$ for which we expect a polynomial depth circuit to
return states with $E\simeq AE_{\mathrm{GS}}$ from first principles. The ideal
value of $A$ depends both on the problem class and on the variation of
spectral folding employed; for MAX-3-XORSAT, $A\geq 0.6$ is achievable as
derived below in section IV.4. This is a significant leap over the best known
classical approximation algorithm for this problem Barak et al. (2015), which
offers a weaker guarantee with much more restricted viability, in the worst
case, equivalent to $A\to 0$ in our notation. It is also a significant leap
over recent quantum approaches to this problem Anshu et al. (2021); Marwaha
and Hadfield (2022). We note that minimizing $\left(H-E\right)^{2}$ is not
itself a novel idea and has been used in classical and quantum algorithms for
finding states close to specific energies in chemical and many-body systems
Wang and Zunger (1994a, b); McClean et al. (2016); Santagati et al. (2018);
Zhang et al. (2021); Tilly et al. (2022); Tazi and Thom (2023). To our
knowledge, however, the use of spectral folding for approximate optimization
of CSPs is novel, both in concept and in the analysis we present below to
choose $A$ and understand at a deeper level why it presents significant
advantages over optimizing the problem $H_{P}$ directly.

From hereon, we let $H_{\mathrm{fold}}$ be the spectrally folded problem
Hamiltonian. Having defined it, there are a number of ways we can attempt to
find its ground states. The simplest choice, and one that performs well, is
AQC/QAOA-inspired state preparation, where we interpolate in Trotterized
evolution between the transverse field driver and the folded problem over a
time $t_{f}$:

$\displaystyle
H\left(t\right)=f\left(t\right)H_{D}+g\left(t\right)H_{\mathrm{fold}},H_{D}=-\sum_{j}X_{j},$
(6) $\displaystyle
f\left(0\right)=g\left(t_{f}\right)=1,\;\;f\left(t_{f}\right)=g\left(0\right)=0.$
(7)

This prescription, with the quadratic folding choice in Eq. (5), is the most
straightforward to benchmark using standard quantum simulation packages as
each call requires $O\left(N_{C}^{2}\right)$ multiqubit
$Z_{i}Z_{j}Z_{k}Z_{l}...$ rotations. Classically simulating the linear folding
prescription requires auxiliary qubits to implement the absolute value
operation or, much more practically, saving the phase oracle as a pre-computed
diagonal operator.

We can also consider trial minimum annealing (TMA), originally proposed in
Kapit (2021a). We explore the TMA formulation in depth because we can predict
its scaling analytically. In this scheme, a simple classical algorithm is used
to find an initial local minimum of $H_{P}$; the quality of the minimum does
not particularly matter and for approximation-hard instances we assume it is
far above the true ground state energy, in the worst case asymptotically
approaching random guessing. Let this classical minimum state be
$\left|L\right>$. We will use the linear folding prescription in Eq. 5 for
$H_{P}$ itself, and add to it a new, diagonal _lowering Hamiltonian_ $H_{L}$
which has $\left|L\right>$ as its ground state, and assign to it a time-
dependent coefficient $C\left(t\right)$. Our total cost function Hamiltonian
is

$\displaystyle
H_{\mathrm{cost}}\left(t\right)=\frac{\left|H_{P}+AN\right|}{A}+C\left(t\right)H_{L}.$
(8)

Recall that $H_{P}$ is normalized so that its ground state energy is $-N$. To
go further, we need to specify a form for $H_{L}$. For this analysis will
choose a new random hypergraph of $N_{C}$ triples which, critically, has no
correlation to the hypergraph of $H_{P}$; we choose the same $N_{C}$ as the
problem for convenience here but any $O\left(N\right)$ quantity should be
fine. We choose the signs of all constraints so that $\left|L\right>$
satisfies all of them. $H_{L}$ is not included in the folding procedure so
applied separately in time evolution. We then choose $C\left(t=0\right)$ such
that the initial energy of $\left|L\right>$ (defined by
$H_{\mathrm{fold}}+C\left(t\right)H_{L}$) is well below $-N$ but remains
$O\left(N\right)$. Our algorithm simulates appropriately discretized time
evolution in the following sequence:

  * •

Initialize $\left|L\right>$, with $H_{\mathrm{cost}}$ always on, and evolve
time smoothly ramping up the transverse field from 0 to $\kappa$ in time
$t_{r}$. We assume $t_{r}$ increases linearly with $N$ and choose
$\kappa\leq\kappa_{c}$; $\kappa_{c}\simeq 1.3$ for MAX-3-XORSAT but can vary
for other problems, and we expect weak variation from one instance to the
next. We want to choose $\kappa$ at or just below this value, so we remain in
the dressed problem phase (DPP, defined in section IV) at all times. Leaving
and then re-entering the DPP does not mean the algorithm will fail but makes
predictions harder. This smoothly evolves the state to $\left|L_{D}\right>$,
the dressed version of $\left|L\right>$.

  * •

Evolve time for a total time $T$, likely also $O\left(N\right)$, where
$C\left(t\right)$ is smoothly ramped down to zero, ensuring that
$\left|L\right>$ crosses the hyperspherical shell of ground states of
$H_{\mathrm{fold}}=\left|H_{P}+AN\right|/A$. Note that this crossing occurs
when the ground state energy of $C\left(t\right)H_{L}$ is $O\left(-N\right)$,
and if we assume the initial minimum was uncorrelated with the true ground
state $\left|G\right>$, the mean Hamming distance between $\left|L\right>$ and
any of the ground states of the folded Hamiltonian is $N/2$ flips.
Consequently, $H_{L}$ adds an $O\left(\sqrt{N}\right)$ energy uncertainty to
these states that has no meaningful impact on the approximation ratio.

  * •

Finally, ramp the transverse down to zero smoothly over $t_{r}$ and measure
the system in the $z$ basis. For an appropriate $O\left(1\right)$ choice of
$A$ and a random problem hypergraph, this algorithm will return states with
energies close to $AE_{\mathrm{GS}}$ with constant probability. We can
optionally repeat the algorithm many times, starting from different choices of
$\left|L\right>$, to ensure a fairer sampling of states in that energy range.

The total gate count of this algorithm is as follows. We have a factor of
$O\left(N+N_{C}\;{\rm polylog}\left(N\right)\right)$ per timestep for the
layers of transverse field, $H_{\mathrm{fold}}$ and $H_{L}$ terms, which we
simplify to $N_{C}\;{\rm polylog\left(N\right)}$. We obtain, in the worst
case, a factor of $O\left(N_{C}\right)$ for the number of guesses one needs to
make to correctly set the normalization for a chosen $A$. We assume, on
empirical grounds, that the total quantum evolution time is
$O\left(N\right)$.This choice works well in practice in our simulations, and
more intuitively, the very simplest classical optimization routine, steepest
descent, requires $O\left(N\right)$ Hamiltonian calls to halt. We do not think
it reasonable, ultimately, that a quantum algorithm should perform well with
fewer steps per shot. Finally, in the worst case we expect a timestep $dt\sim
1/N$, for graphs where a small number of variables connect to significant
fractions of the $N_{C}$ constraints, but $dt$ constant or increasing
logarithmically is empirically and intuitively fine in the typical case. Taken
together, and we emphasize assuming that the algorithm is capable of returning
states with $E<AE_{\mathrm{GS}}$ in constant probability, we estimate a total
runtime between $O\left(N_{C}^{2}N^{2}\;{\rm polylog}\left(N\right)\right)$ in
the worst case and $O\left(N_{C}N\;{\rm polylog}\left(N\right)\right)$ in more
typical cases.

All that said, justifying the assumption of constant success probability is
the central task of the paper. We now provide a theoretical analysis of the
performance of this formulation of spectral folding on random hypergraph
PPSPs, and a more qualitative analysis of the expected performance of other
variations.

Figure 3: Illustration of the spectral folding procedure. (Top) Sketch of the
rough energy landscape of an approximation-hard CSP, with a single deep
minimum whose basin of attraction is an exponentially small fraction of the
configuration space. Directly optimizing this cost function through quantum
approaches often misses the deep minimum entirely, for reasons explained in
the text. (Bottom) Spectrally folded energy landscape in Eq. (5), where the
problem Hamiltonian energies are mirrored around an approximation target
$E=AE_{\mathrm{GS}}$. This can be implemented in a gate model algorithm with
modest overhead, as shown in the text. Doing so promotes the states near the
fold to an exponentially large ground state band while eliminating an
interference effect that reduces tunneling into them from trivial initial
states; for a wide range of problem instances (and likely, low-order problem
classes), this works out to an approximation guarantee. Detailed performance
predictions, and numerical benchmarking, are shown in the text.

##  IV Analytical performance predictions

###  IV.1 Preliminaries

In this section, we will predict from first principles the average macroscopic
quantum tunneling rate–and thus, achievable approximation ratio–for the TMA
variation of spectrally folded quantum optimization. This version may not be
the optimal choice, and there are intuitive reasons to believe that other
variations could offer better performance, but being able to make direct
analytical predictions enormously strengthens our argument and bolsters the
scaling expectations one can infer from our numerical results. To do so, we
have developed a somewhat novel resummed extensive order perturbation theory
based on previous _forward approximation_ results Bravyi et al. (2011);
Pietracaprina et al. (2016); Baldwin et al. (2016, 2017); Scardicchio and
Thiery (2017); Baldwin and Laumann (2018); Grattan et al. (2023).

For the random hypergraph problems we study here, the two key factors in
determining the macroscopic quantum tunneling rate are the transverse field
strength $\kappa_{c}$ where a phase transition occurs between the paramagnet
and the quantum spin glass, which we call the dressed problem phase (DPP) in
this work, and the energy cost $E\left(x\right)$ for $x$ random flips away
from the ground state. We first derive the energy cost, and remarkably, for
random flip sequences it turns out to be graph independent. Specifically, for
MAX-3-XORSAT, our problem is defined as a hypergraph of $N_{C}$ $p$-body
constraints (e.g. $V_{ijk}Z_{i}Z_{j}Z_{k}$) over $N$ variables, where $p=3$
here, and each constraint returns $\pm 1$ and flips to the opposite value when
any one of the spins flips. Let us say the system is in some classical
configuration $s$; the energy is then given by
$E\left(s\right)=N_{C}\left(n_{\mathrm{unsat}}-n_{\mathrm{sat}}\right)$, where
a sat constraint returns 1 in this notation, and $n$ implies a density.

Now we flip one spin at random. Each spin participates in, on average,
$pN_{C}/N$ constraints, and consequently, the average energy change for a
single spin flip is

$\displaystyle\Delta
E_{\mathrm{avg}}=+2p\frac{N_{C}}{N}\left(N_{\mathrm{sat}}-N_{\mathrm{unsat}}\right)=-2pE$
(9)

Now imagine we have flipped $y$ spins from our initial configuration. If we
flip one more spin at random, once again $\Delta
E_{\mathrm{avg}}\left(y\right)=-2pE\left(y\right)\Delta y$. However, we have
already flipped $y$ spins, so when we flip one more at random, with
probability $1-y/N$ we have flipped a spin back and are computing the energy
change associated with reducing $y$ by 1. Consequently

$\displaystyle\left(1-\frac{2y}{N}\right)\frac{\Delta
E_{\mathrm{avg}}\left(y\right)}{\Delta y}=-2pE\left(y\right)$ (10)

If we interpret this as a differential equation, it has a straightforward
solution: starting from the ground state $G$, for $x$ unique random flips away
the average energy is

$\displaystyle
E_{\mathrm{avg}}\left(x\right)=E_{\mathrm{GS}}\left(1-\frac{2x}{N}\right)^{p}.$
(11)

Note that this statement is graph independent, and is only an average;
individual trajectories will of course display substantial variations. It is a
rederivation of a familiar result for dense hypergraphs with Gaussian
distributed constraint energies Derrida (1980); Baldwin and Laumann (2018),
but is applicable in much broader contexts and is easy to confirm numerically.
Throughout this work, we rescale all problems by a multiplicative constant so
that the ground state energy is $-N$.

###  IV.2 Paramagnet to spin glass transition scaling for MAX-3-XORSAT

Before proceeding to our main calculation, we note that, from this result, we
can predict the typical case difficulty scaling of finding the ground state
with AQC or QAOA, for random hypergraphs as defined in section II.4. To
predict the phase transition rate, we will use high order perturbation theory
to compute the perturbative dressings to the problem ground state; call this
state $\left|G_{D}\right>$. Using fourth order perturbation theory, for
uniform transverse field strength $\kappa$, the energy of the dressed ground
state is, again on average444In this step the graph structure of $H_{P}$ can
potentially be important. This is because the energy denominator $2p$, the
average cost per flip just a few flips away from $G$, is more sensitive to the
details of the graph than the energy cost many flips away, which can in turn
shift $\kappa_{c}$. We thus do not claim this result is valid for all graphs.
However, we find that for our random PPSPs the transition is consistently near
the $\kappa_{c}$ we predict here and the scaling of the uniform field result
closely matches our prediction.

$\displaystyle
E_{\mathrm{GS}}\simeq-N\left(1+\frac{\kappa^{2}}{2p}+\frac{\kappa^{4}}{8p^{3}}+...\right)\,.$
(12)

This crosses the energy of the paramagnetic state, which is $-N\kappa$ and
minimally perturbed by the problem Hamiltonian, at $\kappa_{c}\simeq 1.29$ for
our family of $p=3$ PPSPs. The splitting $2\Omega_{0}$ at the phase transition
is expected to be proportional to the overlap of the uniform superposition
state $\left|S\right>$, the ground state of $H_{D}$, with the dressed state
$\left|G_{D}\right>$:

$\displaystyle\Omega_{0}\propto\left<S|G_{D}\right>.$ (13)

To compute $\Omega_{0}$, we must thus compute the perturbative corrections to
$\left|G_{D}\right>$, to high orders. Again using Eq. (11), if we let
$\left|G^{(i,j,k...)}\right>\equiv X_{i}X_{j}X_{k}...\left|G\right>$, and
$\tilde{E}_{\mathrm{avg}}\left(k\right)\equiv
E_{\mathrm{avg}}\left(k\right)-E_{G}$, at the transition point we have

$\displaystyle\left|G_{D}\right>\simeq\left|G\right>+\frac{\kappa_{c}}{\tilde{E}_{\mathrm{avg}}\left(1\right)}\sum_{j}\left|G^{(j)}\right>$
$\displaystyle+2!\frac{\kappa_{c}^{2}}{\tilde{E}_{\mathrm{avg}}\left(1\right)\tilde{E}_{\mathrm{avg}}\left(2\right)}\sum_{i\neq
j}\left|G^{(i,j)}\right>$
$\displaystyle+3!\kappa_{c}^{3}\prod_{m=1}^{3}\frac{1}{\tilde{E}_{\mathrm{avg}}\left(m\right)}\sum_{i\neq
j\neq k}\left|G^{(i,j,k)}\right>+...$ (14)

The factorials come from the combinatorics of ordering the $m$ spin flips to
reach each term. Now, since all states are present in $\left|S\right>$ with
equal amplitude $2^{-N/2}$ and all terms in $\left|G_{D}\right>$ are positive
definite, we can immediately conclude

$\displaystyle\Omega_{0}\left(N\right)\simeq
2^{-N/2}\left(1+\sum_{m=1}^{N}\kappa_{c}^{m}\binom{N}{m}m!\prod_{n=1}^{m}\frac{1}{\tilde{E}_{\mathrm{avg}}\left(n\right)}\right).$
(15)

For $\kappa_{c}=1.29$, this function is well fit by
$\Omega_{0}=a\sqrt{N}2^{-bN}$, where $b\simeq 0.14$, in decent agreement with
the result for a mean-field $p=3$-spin ferromagnet derived in Jörg et al.
(2010). We note also the hardness equivalence between dense and random sparse
graphs drawn in Basso et al. (2022), though we emphasize that the results we
derive here are not limited to completely random graphs. If we are in the
diabatic regime for $t_{f}\ll 1/\Omega_{0}\left(N\right)$, the probability of
finding the ground state after a sweep is
$P_{\mathrm{GS}}\left(t_{f}\right)\simeq\Omega_{0}^{2}t_{f}/W$, where $W\sim
O\left(N\right)$ is the energy range swept over in the evolution; for
$t_{f}\propto N$ we thus have $P_{\mathrm{GS}}\left(t_{f},N\right)\simeq
N2^{-2bN}$. Our numerical simulations for PPSPs, shown below, are in good
agreement with this result. Note that this assumes a single avoided crossing
and does not consider transverse field chaos, e.g. crossings late in evolution
between states close in energy; such transitions can in principle be
substantially more difficult Schiffer et al. (2023), but are not relevant to
approximation hardness in this problem as argued earlier.

###  IV.3 Scaling of inter-valley tunneling in a $p=3$ quantum spin glass

This calculation is less easy to generalize to a spectrally folded
Hamiltonian, however, where we have exponentially many competing ground
states, clustered in a thin hyperspherical shell around the true minimum, at
least in approximation-hard problems. So instead we will consider tunneling
between two $p$-spin wells in an $N$ spin system, spaced $N/2$ flips apart,
and show that the average per-state tunneling rate in the TMA formulation of
spectrally folded optimization should have near-identical scaling. Our total
$N$-spin Hamiltonian, consequently, is

$\displaystyle
H=-\frac{1}{N^{p-1}}\left[\left(\sum_{j}Z_{j}a_{j}\right)^{p}+\left(\sum_{j}Z_{j}b_{j}\right)^{p}\right]$
$\displaystyle-\kappa\sum_{j}X_{j},$ (16)

where the $a_{j}$ and $b_{j}$ are all equal to $\pm 1$ and specify the minimum
position for each of the two terms. We assume for the remainder of this
writeup that $a$ and $b$ correspond to bitstrings $M=N/2$ flips apart; we let
these states be $\left|0\right>$ and $\left|1\right>$. We will also define the
bare classical energy, with no corrections from transverse fields. We start
from one of the two minima, and consider a state which is $m+n$ random flips
away from it. We let $m$ of these flips be ones which move toward the other
minimum (e.g. reduce the Hamming distance to it), and the $n$ flips be flips
that move away from it. Then the bare energy $E_{m,n}^{(0)}$ is given by:

$\displaystyle
E_{m,n}^{(0)}=-N\left[\left(1-2\frac{m+n}{N}\right)^{p}+\left(2\frac{m-n}{N}\right)^{p}\right].$
(17)

We assume that the transverse field strength $\kappa$ is below $\kappa_{c}$,
the $p$-dependent critical point where a transition to the paramagnetic state
occurs. The ground states are thus the symmetric and antisymmetric
combinations of the two dressed classical minima, with splitting
$2\Omega_{0}$, where $\Omega_{0}$ decays exponentially in $N$ and our goal in
this section is to predict its decay rate.

Computing $\Omega_{0}$ proceeds through the following steps:

  * •

We compute the renormalized cost per flip away from either minimum,
incorporating transverse field corrections, which we will then use in the
energy denominators of our $M$th order perturbation theory. This step is
analogous to commonly used resummation schemes in diagrammatic quantum field
theory, where self-energy corrections are incorporated into the propagators
used to compute higher order processes.

  * •

We divide the system between primary spins, which flip between the classical
minima, and secondary spins, which do not. We then compute the dressed states
$\left|0_{D}\right>$ and $\left|1_{D}\right>$ that comprise all the primary
flip sequences up to order $M/2$ away from each minimum. It is at this order
that the two states have nonzero overlap.

  * •

These dressed states are then normalized; incorporating this normalization,
their overlap gives the primary spin contribution to the tunneling rate,
$\Omega_{0}^{(p)}$.

  * •

We then compute the secondary spin contributions to tunneling, which take two
forms: an increase of the tunneling rate from the constructive contribution of
many additional tunneling sequences in which secondary spins participate, and
a decrease from normalization corrections and the spread of the classical
minima away from the core classical configurations from which the tunneling
calculation begins.

  * •

Incorporating both sets of secondary spin contributions gives us a closed form
expression for $\Omega_{0}$ which can then be evaluated numerically and
compared to exact diagonalization.

We first want to compute the energy shifts, in second order perturbation
theory, to these states. These corrections arise from a single spin being
flipped and flipped back, and are opposite in sign to the cost of the local
flip. Let $u_{m,n}$ be the difference between energies of a state $m,n$ and a
ground state, incorporating these corrections. Then

$\displaystyle
u_{m,n}=E_{m,n}^{(0)}+N\left(1+\frac{\kappa^{2}}{2p}\right)-\left(\frac{N}{2}-2m\right)\frac{\kappa^{2}}{\partial_{m}E_{m,n}^{(0)}}$
$\displaystyle-\left(\frac{N}{2}-2n\right)\frac{\kappa^{2}}{\partial_{n}E_{m,n}^{(0)}}+O\left(\kappa^{4}\right).$
(18)

Note that $\partial_{m}E_{m,n}|_{m,n=0}=\partial_{n}E_{m,n}|_{m,n=0}=2p$, so
$u_{0,0}=0$. One can observe that if $p=2$,
$\partial_{m}E_{m,0}=4\left(1-\frac{4m}{N}\right)$, and the transverse field
corrections to state energies are $m$-independent, so that the energy barrier
between the two competing ground states is not renormalized by the transverse
field. But for $p=3$ and higher these corrections are nontrivial and act to
reduce the effective energy barrier between the states, increasing the
tunneling matrix element. This process is effectively a resummation of higher
order corrections and is necessary to obtain quantitatively accurate results.

We start by computing the dressed states, summing over primary spin
corrections only. They take the form

$\displaystyle\left|0_{D}\right>\equiv\left|0\right>+\sum_{j}\frac{\kappa}{u_{1,0}}X_{j}\left|0\right>+\sum_{j,k\left(j\neq
k\right)}\frac{2\kappa^{2}}{u_{1,0}u_{2,0}}X_{j}X_{k}\left|0\right>$
$\displaystyle+\sum_{j,k,l\left(j\neq k\neq
l\right)}\frac{3!\kappa^{3}}{u_{1,0}u_{2,0}u_{3,0}}X_{j}X_{k}X_{l}\left|0\right>+...$
(19)

The expression for $\left|1_{D}\right>$ is identical. We stop our expansion at
order $M/2$, which is the lowest nontrivial order needed to connect the
states. For simplicity we assume $M$ is even though the argument is easy to
generalize to odd $M$ as well. Note that this state is not normalized, and in
fact the norm of the state written above is exponentially large, so we will
need to incorporate normalization corrections into the definition of the
states. Thanks to the dressing of the states, we obtain a primary-spin energy
splitting

$\displaystyle\frac{1}{2}\left(\left<0_{D}\right|+\left<1_{D}\right|\right)H_{P}\left(\left|0_{D}\right>+\left|1_{D}\right>\right)=2\Omega_{0}^{(p)},$
$\displaystyle\frac{1}{2}\left(\left<0_{D}\right|-\left<1_{D}\right|\right)H_{P}\left(\left|0_{D}\right>-\left|1_{D}\right>\right)=0.$
(20)

Evaluating these expressions, the degeneracy splitting from only considering
primary spins is:

$\displaystyle\Omega_{0}^{(p)}$ $\displaystyle=$
$\displaystyle\kappa^{M}\binom{M}{M/2}\frac{1}{u_{M/2,0}}\left(\frac{M}{2}!\prod_{k=1}^{M/2-1}\frac{1}{u_{k,0}}\right)^{2}\times\mathcal{N}^{(p)},$
$\displaystyle\mathcal{N}^{(p)}$ $\displaystyle=$
$\displaystyle\left(1+\sum_{k=1}^{M/2}\binom{M}{k}\left(\kappa^{k}k!\prod_{j=1}^{k}\frac{1}{u_{j,0}}\right)^{2}\right)^{-1}.$
(21)

Here, $\mathcal{N}_{p}$ is the normalization correction. This covers the
primary spin portion of the macroscopic quantum tunneling rate.

We now turn to the secondary spins. To introduce secondary spin corrections,
consider a single secondary spin $j$ out of the $\left(N-M\right)\sim N/2$
total, whose bit value is the same in both classical minima. Since the same
transverse field is acting on it as all other spins, when we consider the sum
of all processes that connect the two minima, we can now divide them between
those where $M/2$ primary spins flip from each minima to meet in the middle,
with secondary spin $j$ unchanged, and a new set of processes where spin $j$
flips starting from each minimum and the two wavefunctions overlap at the set
of states where $M/2$ primary spins have flipped along with $j$. The first set
of processes is what was considered in Eq. IV.3; the second is new, and we
want to calculate its matrix element. It is most useful to express these
matrix elements as a ratio of the new term to the original, primary-spin-only
process, since both decay exponentially in $M$. Let the primary spin
perturbative matrix element to reach $M/2$ flips be $\xi_{M/2}$, so that

$\displaystyle\xi_{M/2}=\kappa^{M/2}\left(\frac{M}{2}!\prod_{k=1}^{M/2}\frac{1}{u_{k,0}}\right).$
(22)

To define the analogous process where $j$ flips, we need to sum over all the
points during the perturbative sequence when that can happen. We thus have:

$\displaystyle\xi_{M/2}^{s}=\kappa^{M/2+1}\frac{M}{2}!\sum_{n=1}^{M/2}\prod_{k=1}^{n}\frac{1}{u_{k,0}}\prod_{k=n}^{M/2}\frac{1}{u_{k,1}}.$
(23)

And noting that we have to make this insertion in the matrix elements from
both minima, the total tunneling term is increased by

$\displaystyle\Omega_{0}^{(p)}\to\left(1+\left(\frac{\xi_{M/2}^{s}}{\xi_{M/2}}\right)^{2}\right)\Omega_{0}^{(p)},$
(24)
$\displaystyle\left(1+\left(\frac{\xi_{M/2}^{s}}{\xi_{M/2}}\right)^{2}\right)\equiv\gamma_{T}.$
(25)

We now need to consider the rest of the secondary spins. Formally of course,
the additional energy cost of each secondary spin flip changes as more
secondary spins flip, but since these corrections are fairly weak (though they
are appreciable and necessary for an accurate prediction of the scaling
exponent) the total tunneling rate is going to be dominated by the set of
processes where a comparatively small fraction of secondary spins have
flipped, and we can thus approximate them as independent contributions. In
this limit, since there are $N-M$ secondary spins,

$\displaystyle\Omega_{0}^{(p)}\to\Omega_{0}^{(p)}\gamma_{T}^{N-M}.$ (26)

Alongside this, the secondary spin corrections also spread the competing
ground state wavefunctions out over Hilbert space, which exponentially reduces
the weight of the core classical configurations from which the tunneling
calculation begins. To be consistent with the independence approximation made
above, we simply compute all the corrections to the ground state from each
secondary spin independently and multiply them. Noting that we must apply this
calculation to both competing ground states, this reduces the tunneling rate
by

$\displaystyle\Omega_{0}^{(p)}\to\Omega_{0}^{(p)}\left(\frac{\gamma_{T}}{\gamma_{R}}\right)^{N-M},\;\;\gamma_{R}\equiv
1+\left(\frac{\kappa}{u_{0,1}}\right)^{2}.$ (27)

Note that, if we set the cost per flip of a given secondary spin to some
constant $U$, independent of the configuration of the other spins, that would
imply it is disconnected from the primary spins as there are no couplings to
shift the energy. In this limit a direct evaluation of the two functions shows
that $\gamma_{T}=\gamma_{R}$ (for any choice of $U$) and this now disconnected
spin plays no role in tunneling at all. This factorization of disconnected
spins is reassuring, and lends support to the correctness of this approach.
Taking into account all these effects, our total tunneling rate is

$\displaystyle\Omega_{0}=\kappa^{M}\binom{M}{M/2}$
$\displaystyle\frac{1}{u_{M/2,0}}\left(\frac{M}{2}!\prod_{k=1}^{M/2-1}\frac{1}{u_{k,0}}\right)^{2}$
$\displaystyle\times\mathcal{N}^{(p)}\times\left(\frac{\gamma_{T}}{\gamma_{R}}\right)^{N-M}.$
(28)

Taking all of these effects into account yields a highly accurate prediction
of the minimum gap scaling for a wide range of values for $p$ and $\kappa$,
with only an $O\left(1\right)$ discrepancy in the prefactor and few percent
discrepancies in the scaling exponent (empirically, Eq. IV.3 tends to slightly
overestimate the decay compared to the exponent extracted from numerical
diagonalization). We refer to the Appendix A for more details.

###  IV.4 Achievable approximation ratio with spectrally folded trial minimum
annealing

With this result in hand, we will now predict the macroscopic quantum
tunneling rate–and thus, achievable approximation ratio–for the spectrally
folded Hamiltonian. From this, we can calculate our target value of $A$,
assuming that there is a single deep minimum far below the energy of any local
minima. Relaxing this assumption will improve the performance of the algorithm
by virtue of there being many more target states, as confirmed in our
simulations for lower constraint densities. We consider the protocol in
section III, with an initial state $\left|L\right>$. We assume that the
process of ramping the transverse field up and down is itself at least roughly
adiabatic, i.e., we can assume approximate spectral continuity with respect to
the folding and lowering Hamiltonians, noting that the lowering Hamiltonian
will itself create $O\left(\sqrt{N}\right)$ shifts to the energies of states
near the fold. It follows from our assumption that the ramping process itself
does not meaningfully heat the system. We then consider the set $\mathcal{T}$
of all states within $O\left(1\right)$ shifts of $-AE_{\mathrm{GS}}$ in
$H_{P}$, the states closest to the fold, and compute, as a function of all our
various algorithm parameters, the total probability of tunneling into any one
of them.

Since the tunneling rate into any individual state is exponentially small, and
the time over which we slowly turn off the lowering Hamiltonian is $T\propto
O\left(N\right)$, we can assume that tunneling will be diabatic with respect
to any individual state. A Fermi’s Golden rule analysis as in Kapit and
Oganesyan (2021) suggests that

$\displaystyle
P_{\mathrm{tot}}\propto\frac{T}{W}\sum_{j\in\mathcal{T}}\Omega_{0,Lj}^{2},$
(29)

where $W\sim O\left(N\right)$ is the energy range swept over by reducing
$C\left(t\right)$ to 0, and $\Omega_{0,Lj}$ is the tunnel splitting at
degeneracy between $\left|L_{D}\right>$ and the target state
$\left|j_{D}\right>$, which we assume are an average of $\sim N/2$ flips
apart. To go further, we need to compute the average value of
$\Omega_{0,Lj}^{2}$, noting that while of course there will be substantial
state-to-state variations, given that there are exponentially many states in
$\mathcal{T}$ the average value is going to dominate Eq. (29). As in the
previous calculation, the most important quantity here is the average cost per
flip away from the typical ground state in the folded Hamiltonian, which
remarkably turns out to be $A$-independent.

To see this, we start from the average cost per flip away from
$\left|G\right>$, given by Eq. (11), and note that we can invert that equation
to find the mean number of flips $x_{A}N$ for which
$\left<E\left(x_{A}N\right)\right>=-AN$. To be specific,

$\displaystyle x_{A}=\frac{1-A^{1/3}}{2}.$ (30)

We can therefore assume that the typical state in $\mathcal{T}$ is $x_{A}N$
flips away from $\left|G\right>$. If we consider the sequences of primary spin
flips connecting $\left|j_{D}\right>$ and $\left|L_{D}\right>$, the typical
flip sequence starts $x_{A}N$ flips away from $\left|G\right>$, and notice
that with probability $1-x_{A}$ an additional random flip towards
$\left|L\right>$ will also move closer to $\left|G\right>$. Taking all these
effects and the division by $A$ into account, so that the bare unperturbed
energy of both $\left|j\right>$ and $\left|L\right>$ when they cross are both
$\sim-N$, a bit of algebra shows that for $y$ flips away from a ground state
of the folded Hamiltonian, not only is the total average cost $\Delta
E\left(y\right)$ $A$-independent, it is precisely equal to the cost given by
Eq. (11).

This is again only an average, but noting that since it appears the
denominators of equations like (IV.3), variations about it are more likely to
increase the tunneling rate than decrease it. And likewise, since $H_{L}$ is a
random 3-XORSAT problem itself the mean cost per flip away from
$\left|L\right>$ is going to be given by Eq. 11 as well, so Eqns. IV.3 through
IV.3 can faithfully predict the average tunneling rate between
$\left|L\right>$ and a randomly chosen ground state of the folded
Hamiltonian.555We expect that using the average cost per flip in Eq. 17 will
if anything underestimate the per-state tunneling rate in real disordered
problems. This is because all of these energy costs appear in denominators,
which leads to likely small asymmetries in how much the deviations from the
average in any individual flip sequence contribute to the total matrix
element, giving lower energy sequences proportionally higher weight. We do not
really expect this effect to be significant but rather highlight it as another
point where our prediction is conservative by design. For more discussion of
this approximation, see the appendix.

Of course, this rate decays exponentially; assuming the two states are $M\sim
N/2$ flips away for $\kappa=1.29$, $\Omega_{0}\propto 2^{-bN}$ where $b\simeq
0.2$. But this is balanced by the fact that there are on the order of
$\binom{N}{x_{A}N}$ target states. We can further note that out of these
states, while the mean distance to $\left|L\right>$ is $M\sim N/2$, ones which
are $k$ flips closer have tunneling rates which are larger by a factor of
$2^{2bk}$ on average, and though those states are proportionally rare their
increased weight is enough to meaningfully impact our choice of $A$. Since our
total runtime is linear, simple diabatic scaling predicts that the probability
of tunneling into the typical state $M-k$ flips away is proportional to
$\Omega_{0}^{2}$, e.g. $2^{-4b\left(M-k\right)}$.

We now take this result and plug it into Eq. (29), so that we can determine
the choice of $A$ where the returned $P_{\mathrm{tot}}$ provides an
approximation guarantee. If we use Stirling’s approximation to write the
binomial coefficients as exponentials, and ignore slowly varying polynomial
factors, the total number of states in $\mathcal{T}$ scales as:

$\displaystyle N_{\mathcal{T}}\propto\exp\left(-\left[x_{A}\ln
x_{A}+\left(1-x_{A}\right)\ln\left(1-x_{A}\right)\right]N\right).$ (31)

Likewise, if the average probability of tunneling into a target state $k$
primary spin flips closer to $\left|L\right>$ is increased by a factor of at
least $2^{4bk}$, the weighted per-state average of the diabatic tunneling rate
into states in $\mathcal{T}$ can be approximated as

$\displaystyle\frac{\log\left<\Omega_{0}^{2}\right>}{-N}\approx 2b\ln
2-x_{A}\left(\ln\left(1+2^{4b}\right)-\ln 2-2b\ln 2\right).$ (32)

Note that this average comes from considering only $x_{A}N$ flips away from
$\left|G\right>$ but varying Hamming distance from $\left|L\right>$ and thus
neglects the influence of comparatively rarer states larger distances from
$\left|G\right>$. Taking all these terms into account, the probability of
returning a state with $E\simeq AE_{\mathrm{GS}}$ as measured relative to the
original $H_{P}$ becomes constant, or at least stops decaying exponentially,
when

$\displaystyle-2b\ln 2-\left[x_{A}\ln
x_{A}+\left(1-x_{A}\right)\ln\left(1-x_{A}\right)\right]+$ $\displaystyle
x_{A}\left(-\ln 2-2b\ln 2+\ln\left(1+2^{4b}\right)\right)=0.$ (33)

The achievable approximation ratio is thus determined by the per-state decay
exponent $b$, computed in section IV.3 as a function of $N$ and $\kappa$ by
fitting Eq. (IV.3) to $\Omega_{0}\left(N\right)\propto\sqrt{N}2^{-bN}$, and
then choosing $x_{A}$ using Eq. (30) to solve Eq. (IV.4). This analysis only
counts states within $O\left(1\right)$ shifts of $AE_{\mathrm{GS}}$ (recall
that $E_{\mathrm{GS}}=-N$ in our normalization) and ignores low-order
polynomial prefactors; for $b=0.2$, which again depends on $\kappa=1.29$ in
this calculation, this is solved when $x_{A}\simeq 0.08$, or $A\simeq 0.59$.

This means that if the true ground state satisfies a constraint fraction $F$
beyond random guessing—e.g. a total fraction $1/2+F$, so $F$ is at most 1/2
here—our algorithm will return states which satisfy a fraction $1/2+AF$ with
high probability. If we choose $A$ to be too large compared to the target
value established by expressions like Eq. (IV.4), we risk failing to well-
approximate the problem; conversely, choosing $A$ below it will reduce the
returned approximation ratio to $A$ and thus perform suboptimally. And, we
emphasize again, this prediction assumes a random, potentially dense
hypergraph but is fundamentally independent of the fraction satisfied in
$E_{\mathrm{GS}}$ itself and so applies to the planted partial solution
instances we use for numerical benchmarking below.

###  IV.5 Further Comments and Caveats

We expect that this analysis underestimates the choice of $A$ that will return
states with $E\leq AE_{\mathrm{GS}}$ with constant probability. This is
because our counting here only counts states very close to the fold, when in
reality the probability of tunneling into states a small extensive fraction
larger than $AE_{\mathrm{GS}}$ is still going to be appreciable due to the
continued exponential growth of the number of targets, even if the per-state
tunneling rate does tend to decrease with increasing $E$ due to the
interference effect mentioned earlier, in which perturbative corrections that
mix with states of lower energy have opposite sign. In addition to this
consideration, because the target ground states and low lying excitations of
the folded Hamiltonian in $\mathcal{T}$ very roughly form a hyperspherical
shell, any individual target state will have other states in $\mathcal{T}$
that are relatively close to it in Hamming distance. Consequently we expect
these states to have a band dispersion, centered around the mean energy given
by the corrections in Eq. (IV.3). Since we are already assuming off-resonant
tunneling, i.e. a per-state tunneling rate $\propto\Omega_{0}^{2}$, and so
summing over squared matrix elements, if we consider states near the band
center where the density is highest this will alter the average tunneling rate
by at most a prefactor. However, there are good reasons to suspect that
tunneling into extremal states near the bottom of the band can be
substantially enhanced, enough to increase the optimal value of $A$. This
calculation is difficult to do quantitatively so we do not attempt it here; we
instead see the approximation of considering only states near the band center
of $\mathcal{T}$ as another choice that likely underestimates the achievable
approximation ratio. We discuss this point in more detail in Appendix B.

We also want to emphasize that this variation is not necessarily the optimal
spectrally folded optimization algorithm, but instead merely the one where we
were able to analytically compute the threshold $A$. For example, one can
perform trial minimum annealing with a simple local $Z$ bias lowering
Hamiltonian (e.g. $H_{L}=\sum_{j}h_{j}Z_{j}$), or standard AQC interpolation
using the quadratic folding procedure in Eq. (5) as the cost function. The
linear lowering Hamiltonian is expected to have equal or better tunneling
rates to a 3-XORSAT-based minimum as the overall cost-per-flip curve is
shallower, though the local energy shifts to the ground states of
$H_{\mathrm{fold}}$ from $H_{L}$ are expected to be larger. The total gate
count at each time step is lower. Empirical performance in testing up through
$N=25$ showed fairly similar performance to 3XOR-based $H_{L}$ for all other
parameters equal, but with more significant non-monotonic behaviors that made
fitting difficult; see Sec. V for details. That the two schemes could
asymptotically converge to the same achievable approximation ratios seems
plausible to us but we cannot simulate large enough systems to be sure.

For quadratic folding AQC as in Eq. (6), if we choose $A=1$ the gap is
efficiently computable using the methods in Jörg et al. (2010) and decays as
$\Omega_{0}\sim 2^{-0.16N}$. Given that, like linear folding, the cost per
flip curve is $A$-independent, if we assume that the tunneling rate per state
for $A<1$ is basically equal to this, then the total decay exponent vanishes
if $x_{A}\simeq 0.06$ and $A\simeq 0.68$ using the arguments of the previous
few paragraphs. We do not think that can be simply assumed as easily as with
tunneling between semiclassical minima and a linearly folded problem
Hamiltonian, in the DPP, and more theoretical work is needed here to
analytically determine the optimal choice of $A$. Interestingly however, our
simulation data in Sec. VI supports this conclusion, with a worst case
polynomial time approximation ratio of 0.7 found in our simulations. These
simulations show that this method performs similarly too, or slightly worse
than, the 3XORSAT-TMA algorithm, which is better able to outperform the
approximation guarantee of $\sim 0.6$ derived here.666For smaller systems we
also tested linear folding AQC and quadratic folding TMA. In very preliminary
studies we found that quadratic AQC modestly outperformed linear AQC, and
linear TMA more significantly outperformed quadratic TMA. So we chose not to
pursue those methods for larger simulations and do not present those results
here, but they may be viable or even superior for other problem classes. We
also expect that the average tunneling rate–and thus, achievable approximation
ratio–can likely be further increased by using other, potentially many-
frequency, AC methods such as RFQA Kapit (2021b); Kapit and Oganesyan (2021);
Tang and Kapit (2021); Mossi et al. (2023); Grattan et al. (2023). For
simplicity, we do not incorporate these methods in this work, but they could
be a novel way to further improve the performance of this algorithm and are
worth exploring in future research.

In summary, through a relatively novel resummed extensive order perturbation
theory, we have shown that random hypergraph MAX-3-XORSAT instances, including
extremal ones with planted partial solutions, are efficiently approximable to
a fairly large constant fraction through the spectrally folded quantum
optimization algorithms. We do not expect this to be the case for QAOA, an
expectation supported by the numerical evidence we present below. We similarly
do not expect such guarantees to be possible for directly finding global
optima, for the reasons set forth in the introduction. Evidently, one cannot
so easily summit a mountain in hyperspace, but one can reach the rim of a
crater. We now present a series of numerical simulations to further support
these claims.

##  V Numerical tests of approximation hardness for traditional methods

###  V.1 Setup and summary of results

To confirm our predictions–or at least, verify that any serious issues with
our calculations and interpretation of the problem are subtle and not apparent
at system sizes within reach of present or near-future classical
simulations–we performed a series of numerical simulations of various
classical and quantum algorithms applied to our PPSPs. For all quantum
simulation tasks we used the Qulacs package Suzuki et al. (2021). For smaller
systems and algorithm prototyping we ran our simulations on local
workstations; this includes all the spectral folding TMA simulations. For all
QAOA and spectral folding AQC simulations presented, we used the Fujitsu
Quantum Simulator, a classical HPC system. This allowed us to probe larger
system sizes while still averaging over enough instances to have reliable
statistics for the problem class. Unless otherwise stated, each datapoint
represents the average over 960 or 1000 (for the QAOA and spectral folding AQC
results) random problem instances constructed with the prescriptions outlined
in section II.4. In all cases in this work we used an unsatisfied fraction
$\epsilon=0.1$ for our partial planted solutions; we expect similar
phenomenological behavior for other constant $\epsilon$ values, though we
expect the precise thresholds we measure will vary with $\epsilon$ for
classical optimization algorithms and QAOA (but not folded optimization, below
its predicted performance floor).

The results of our simulations are summarized in Table 1, which lists the best
per-shot polynomial time approximation ratio achievable through each algorithm
studied, both classical and quantum. To estimate these values, we measured the
average per-shot probability $P_{q}\left(N\right)=P\left(E\leq
qE_{\mathrm{GS}}\right)$ of returning states with energies at or below
$qE_{\mathrm{GS}}$, where $q\leq 1$ is the approximation ratio and
$E_{\mathrm{GS}}$ is negative in our conventions. These probabilities were
computed by choosing bins of size 0.05 $E_{\rm GS}$.

To determined the polynomial time hardness threshold, we applied a very simple
rule where we fitted $P_{q}\left(N\right)$ to a simple exponential function,
assuming any observed decay faster than $2^{-0.005N}$ corresponded
asymptotically to exponential decay, and any positive exponents, e.g.
exponential growth, represent small-$N$ growth toward some constant saturation
value. This threshold of $2^{-0.005N}$ was chosen to reduce the influence of
uncertainty from fitting in a small handful of cases. For the greedy algorithm
with $N_{C}/N=1.5\sqrt{N}$ we fitted decay to an exponential in $\sqrt{N}$ as
discussed below. The polynomial time approximation hardness threshold for a
given algorithm, unsatisfied fraction $\epsilon$ and $N_{C}/N$ scaling choice
is defined to be $q_{a}$, the largest values of $q$ for which we do not
observe exponential decay. For $N_{C}/N=3\sqrt{N}/2$ this threshold is
decaying with system size, and the asterisk next to the result for QAOA is to
highlight the fact that we only ran these simulations out to $N=28$ so are
likely not capturing the asymptotic threshold. As discussed elsewhere, the
number of cost function calls is $O\left(N\right)$ for all approaches studied.
Cases for which a range is quoted are where we felt there was some ambiguity
to the fitting, and all values are the result of extrapolating fits to
numerical simulations and are naturally somewhat approximate.

$N_{C}/N$ | Classical | QAOA | AQC/0.75 | AQC/0.85 | TMA-3/0.75 | TMA-3/0.85 | TMA-L/0.75
---|---|---|---|---|---|---|---
2 | 0.75-0.8 | 0.75 | 0.75 | 0.75 | 0.75 | 0.75 | 0.7
4 | 0.55 | 0.55 | 0.7 | 0.7 | 0.75 | 0.8 | 0.75
6 | 0.45 | 0.45 | 0.75 | 0.7 | 0.8 | 0.8 | 0.75
3 $\sqrt{N}/2$ | decaying | 0.25/decaying* | 0.75 | 0.75 | 0.8 | 0.8 | 0.75-0.8
Table 1: Approximation hardness thresholds for the classical greedy search,
high-depth QAOA and spectral folding variations. This table lists $q_{a}$, the
largest value of the approximation ratio $q$ before exponential decay is
reported, drawn from the numerical experiments in figures 5 through 8 and 14.
Spectral folding results are labeled as protocol/$A$ (where $A$ is the
approximation target); AQC is quadratic spectral folding in the AQC
formulation, TMA-3 is trial minimum annealing with a linear folded Hamiltonian
and 3-XORSAT lowering Hamiltonian, and TMA-L is the same with local $Z$ biases
for the lowering Hamiltonian. These results support the predictions in section
IV that random hypergraph problems are efficiently approximable through
spectrally folded quantum optimization.

###  V.2 Performance of greedy classical algorithms

To explore the classical difficulty of our PPSPs, we applied a greedy local
search algorithm adapted from Bellitti et al. (2021). This algorithm is
straightforward; we start with a random bitstring. Then, beginning at each
step, we calculate $k=$ ”number of unsatisfied constraints minus the number of
satisfied constraints” associated with each bit. We then calculate the
fraction of bits, $f_{k}$, belonging to each $k$ value. Note that we only care
when $k>0$ because these are the cases where flipping a particular bit will
lower the energy. Using some weight function, $w(k)$, we select a $k$ value
with normalized probabilities $\propto w(k)*f_{k}$ and flip a bit with that
$k$ value. If there are multiple bits belonging to the $k$ value chosen, we
choose a bit in this set with uniform probability to flip. This is repeated
until the configuration finds itself in a minimum, and the algorithm halts.

We found the algorithm performed best with a weight profile quadratic with k,
however this can be experimented with for different results. Notably, Bellitti
et al. (2021) found that when applying a highly optimized version of this
search to 3-regular 3-XORSAT problems it performed well even when compared to
more sophisticated algorithms such as simulated annealing and parallel
tempering. The intuitive reason for this can be inferred from the typical
energy landscapes of these problems, which are rough and contain exponentially
many high energy local minima. Once one is found, it is more efficient to
simply restart the algorithm from a new random configuration instead of
attempting to “climb out” using penalized operations in simulated annealing or
parallel tempering. As the locations of these minima are uncorrelated with the
true ground state, finding one provides no useful information in a ground
state search.

Figure 4: Per-shot probability of finding the true ground state with the
quasi-greedy classical algorithm defined in section V.2 as adapted from
Bellitti et al. (2021), for the four constraint densities studied in this work
and planted solution unsatisfied fraction $\epsilon=0.1$. The probability
decays superpolynomially with $N$; fits are to $a2^{-bN}$ for
$N_{C}=\left\\{2,4,6\right\\}N$ and $a2^{-b\sqrt{N}}$ for $N_{C}=1.5N^{3/2}$;
the scaling fairly well matches the empirically observed Eq. (34). Each
datapoint is the average of $10^{5}$ shots. As the constraint density
increases, the basin of attraction widens, and the problem becomes
easier–though still exponentially scaling–for local update classical routines.
Each shot consists of $O\left(N\right)$ local updates so the time to solution
scales essentially as the inverse of this probability.

This expected inefficiency of simulated annealing/parallel tempering for this
problem can easily be inferred from the results plotted in Figs. 4 and 5.
Namely, in all cases the algorithm will find a single relatively deep minimum
with high probability at each shot, leading to a super-polynomial cost to
escape from it in algorithms simulating a thermal bath. Interestingly, as
$N_{C}/N$ increases for fixed unsatisfied ground state fraction $\epsilon$, we
find that the decay exponent of the per-shot probability of finding the ground
state, $P_{\mathrm{GS}}\left(N\right)$, monotonically decreases, suggesting
that the basin of attraction of the true ground state is widening as the
problem becomes more extremal. In fact, for $N_{C}\geq 2N$, we empirically
observe that the per-shot probability of finding the planted ground state has
the approximate scaling

$\displaystyle\log\left(P_{\mathrm{GS}}\left(N\right)\right)\simeq-
c_{g}\frac{N^{2}}{N_{C}}\,.\;\;\;({\rm PPSPs,}\;N_{C}\geq 2N)$ (34)

However, the approximation ratio $q_{a}$–defined as the minimum energy for
which the probability of finding states at or below it stays constant as $N$
increase–steadily worsens. We attribute this to there being a high density of
local minimia with energies $\geq q_{a}E_{\mathrm{GS}}$ (recall
$E_{\mathrm{GS}}$ is negative in our conventions), but below that threshold
the number of minima quickly decreases and the probability of finding one
decreases exponentially. This results in the scaling collapse seen in Fig.
5–for sufficiently low energy there are no minima aside from the ground state,
so the approximation probability scales nearly identically to
$P_{\mathrm{GS}}\left(N\right)$. This _high energy clustering phase_ is a
feature of our PPSP construction, and is responsible for its classical
approximation hardness. We again contrast this to problems near the
statistical SAT/UNSAT threshold such as three-regular instances, where the
_clustering energy_ , the lowest energy where there are still exponentially
many local minima and they are thus easy to find, is close to
$E_{\mathrm{GS}}$ and they are not approximation-hard in practice as a result.
We conjecture that high energy clustering behavior is a generic feature of
low-degree constraint problem classes that are approximation-hard for local
update algorithms.

As shown in Fig. 5, these construction rules yield in a set of instances which
are hard to approximate in practice. If we let $N_{C}/N$ grow slowly with $N$,
e.g. as $\ln\left(N\right)$ or $\sqrt{N}$, then as $N\to\infty$, the
probability of finding any states with energies any $O\left(1\right)$ fraction
better than random guessing, decays superpolynomially in $N$. And since the
unsatisfied fraction $\epsilon$ in the ground state is small but nonzero,
Gaussian elimination cannot be used to efficiently find the solution, forcing
classical computers to rely on local update algorithms stymied by entropic
barriers. It is of course possible that some clever algorithm could be written
to exploit our PPSP structure to efficiently solve or approximate these
instances classically; we merely claim hardness for generic methods based on
local updates. Our PPSP construction rules can easily be generalized to other
CSPs, and we suggest that they could prove to be a useful tool for exploring
practical approximation hardness in other contexts.

Figure 5: Classical approximation hardness of PPSPs using the quasi-greedy
classical algorithm, for constraint densities (clockwise from top left)
$N_{C}=\left\\{2N,4N,6N,1.5N^{3/2}\right\\}$. In each figure we plot the
probability that a given shot returns an energy below $qE_{\mathrm{GS}}$ for
various choices of $q$. For all the fixed $N_{C}/N$ problem classes we find an
empirical approximation threshold $q=A_{g}$ below which finding states becomes
superpolynomially hard, and that this value decreases as $N_{C}/N$ increases.
For the case $N_{C}=1.5N^{3/2}$ (bottom left), this value steadily drops, as
discussed in the text.

###  V.3 Performance of high-depth QAOA for this problem

Figure 6: Performance of QAOA as an approximator for
$N_{C}=\left\\{2,4,6,3\sqrt{N}/2\right\\}N$ (clockwise from top left), with
problem and algorithm parameters described in-text. Plotted are the
probabilities $P_{q}\left(N\right)=P\left(E\leq qE_{\mathrm{GS}}\right)$ for
$q$ running from 0.5 to 0.85 in steps of 0.05 (top panels) and 0.25 to 0.6
(bottom panels). The results for the larger constraint densities thus plot a
weaker approximation range. Thick straight lines correspond to simple
exponential fits where $P_{q}\left(N\right)$ is not decaying, dashed straight
lines correspond to exponential decay, and thin lines between points are
included for visual clarity. As summarized in table 1, these results do not
represent a meaningful improvement over the classical greedy result (Fig. 5),
though in some cases where $P_{q}\left(N\right)$ decays exponentially for both
approaches, the exponent for QAOA may be better. These results were obtained
using the Fujitsu Quantum Simulator, a classical HPC system.

To make firm points of comparison, alongside the simulations of spectrally
folded optimization itself we extensively benchmarked high-depth QAOA and
quasi-greedy Bellitti et al. (2021); Kowalsky et al. (2021) classical
algorithms on these instances. For our high-depth QAOA simulations, we
formulated our algorithm to mimic trotterized time evolution over a total time
$t_{f}$, with:

$\displaystyle\left|\psi\left(t+dt\right)\right>=e^{-2\pi
if\left(t\right)dtH_{D}}e^{-2\pi
ig\left(t\right)dtH_{P}}\left|\psi\left(t\right)\right>,$ (35) $\displaystyle
f\left(t\right)=\sqrt{1-t/t_{f}},\;g\left(t\right)=\sqrt{t/t_{f}}.$

In all the presented data we used $t_{f}=N/32$ and $dt\simeq 0.05$. Individual
shots use a random evolution time between $2t_{f}/3$ and $4t_{f}/3$; this
runtime averaging produces substantially more reliable scaling, particularly
when probabilities are small. These parameters were chosen by trial and error
for smaller systems; we observed that the probabilities of finding the ground
state and other low-energy states increased sublinearly with $t_{f}$ beyond
this point. The relative improvements of the probabilities of returning low-
lying states were similarly sublinear. No sophisticated numerical or iterative
optimization methods were used, for the reasons discussed in section II.5.

The results of our QAOA simulations, of 1000 random PPSPs for each choice of
$N$ running from 8 to 30, are shown in Fig. 6. The probability of finding the
ground state, shown in Fig. 11, decays exponentially with an exponent very
close to that in Eq. (15), which we find remarkable given the simplifying
assumptions in that derivation, and that it does not use the more
sophisticated techniques used to compute tunneling rates between semiclassical
minima. Further, this exponent displays only small variations with constraint
density and is nearly identical in all four cases. Smaller system studies for
other constraint densities all yielded very similar results for
$P_{\mathrm{GS}}$, as predicted by Eq. (15).

Turning to approximation hardness, being relatively sparse, the $N_{C}=2N$
problems are fairly well-approximated by QAOA, with the algorithm returning
strings within $q=0.75$ with constant or saturating probability; we attribute
this to the presence of many competing minima with energies not far from
$E_{\mathrm{GS}}$. In contrast, for $N_{C}=4N$ the algorithm’s performance for
approximation degrades, with clear exponential decay for approximation ratios
better than $q=0.55$. For higher constraint densities approximation becomes
even more difficult, decaying exponentially below $q=0.45$ for $N_{C}=6N$ and
0.25 for $N_{C}=3N^{3/2}/2$. We expect decay at sufficiently large $N$ for any
constant fraction in that case, but cannot simulate larger system sizes.
Crucially, the thresholds $q_{a}$ we measure are nearly identical to those
found by the classical greedy algorithm (figure 5), and no signatures of an
exponential quantum advantage in these instances can be seen.

Interestingly, as $N_{C}/N$ increases, the probabilities of finding states
comparatively close to $\left|G\right>$ in Hamming distance improve (see Fig.
11), but the probabilities of finding states close in energy worsen. We
attribute this behavior to the high energy clustering phase conjectured in
Sec. V.2. Empirically for our PPSPs there is a high density of local minima
with energies $E\geq q_{a}E_{\mathrm{GS}}$, and if $q_{a}$ is relatively close
to 1 it becomes harder for high-depth QAOA to find local excitations near the
planted ground state, as the probability amplitudes will be spread over
increasingly many competing minima and their own basins of attraction.
Conversely, as $q_{a}$ decreases with increasing $N_{C}/N$, the probability of
finding local excitations relatively near $\left|G\right>$ increases though
still decays exponentially, as the low energy minima far from $\left|G\right>$
are at proportionally higher energies and thus do not compete directly with
few-flip states. That the thresholds $q_{a}$ for QAOA match those of the
greedy classical algorithm further supports the interpretation of an energy
threshold above which local minima become common.

Comparing the classical and established quantum methods, we find that, for
these approximation hard instances, QAOA performs very poorly for finding the
ground state but less poorly for approximation below the classical hardness
threshold $q_{a}$ (see Table 1), with decay exponents that are much closer to
the classical result. In some higher constraint density cases our fits
produced favorable exponents for approximation with QAOA but our range of $N$
here is smaller than we would prefer to claim any relative quantum advantage
absent theoretical justification. Nonetheless, both methods show clear super-
polynomial decay per-shot for approximate optimization below the energy range
where local minima are dense. With these results in hand, we now turn to
folded quantum optimization, which maintains an approximation guarantee
regardless of the problem’s constraint density.

##  VI Numerical tests of spectral folding variations

Having numerically confirmed the expectation that our PPSPs are
superpolynomially hard for classical and prior quantum approaches, for both
exact and approximate optimization, we now present the results of our folded
quantum optimization simulations. We simulated both the trial minimum
annealing and interpolation (e.g. AQC) variations. In all cases we chose
runtimes increased linearly in $N$, albeit with larger prefactors than in the
QAOA simulations (which used $t_{f}=N/32$). A longer runtime further helps
reduce the potential influence of diabatic local heating as the Hamiltonian
parameters are varied. Run and ramp times that are too short can lead to
artificially poor scaling, arising from the formation of local excitations as
the transverse field is turned on or off too quickly. This is fundamentally a
different, and much more prosaic, issue, than decaying collective tunneling
rates, but can be difficult to distinguish when our only measures are energy
and Hamming distance from $\left|G\right>$.

###  VI.1 Spectrally folded adiabatic interpolation performance

We first present simulations of the AQC variation of spectral folding, in
figures 7 and 12. For the AQC variation, we followed the procedure in Eq.
(35), with a quadratic $H_{\mathrm{fold}}$ (Eq. 5) in place of $H_{P}$ and
$f\left(t\right)=\left(1-t/t_{f}\right)^{1/4}$ instead of $\sqrt{1-t/t_{f}}$,
with an average runtime $t_{f}=N/24$. This schedule modification was found to
improve scaling at higher approximation ratios. In all cases individual shots
are runtime averaged between $2t_{f}/3$ and $4t_{f}/3$ as in our QAOA
simulations.

All these parameter choices are the result of intuition, trial and error, and
a desire for simulations to complete in reasonable amounts of time; none are
particularly optimal. Nonetheless, as shown in figure 7, in all but one case
we were able to meet the approximation target $A=0.75$, but not exceed it; for
$N_{C}=4N$. Our spectral folding methods are also much better at reliably
returning states close to $G$ in Hamming distance, consistent with the
approximation guarantee (see Fig. 12). We likewise tested increasing $A$ to
0.85 in simulations up to $N=24$ with this variation, and found improved
prefactors but no improvement in scaling (see Fig. 14). This suggests that we
have found the performance ceiling for this approach. Interestingly, the worst
case $q_{a}=0.7$ observed for this method is very close to the approximation
ratio of $\sim 0.68$ predicted in section IV.5, using a more simplified
analysis than was employed for the TMA variation.

Figure 7: Performance of the quadratic AQC formulation of spectral folding,
for $N_{C}=\left\\{2,4,6,3\sqrt{N}/2\right\\}N$ (clockwise from top left),
with $N$ running from 10 to 27 (top row) or 26 (bottom row), $A=0.75$,
$dt=0.0325$, $t_{f}=N/24$ and other parameters as stated in text. In each plot
the 8 curves plot $P_{q}\left(N\right)=P\left(E\leq qE_{\mathrm{GS}}\right)$
for $q$ running from 0.5 to 0.85 (top to bottom) in steps of 0.05. Thick
straight lines correspond to simple exponential fits where
$P_{q}\left(N\right)$ is not decaying, dashed straight lines correspond to
exponential decay, and thin lines between points are included for visual
clarity. In all four cases spectrally folded optimization is able to meet its
approximation target of $A=0.75$, returning states at or below this energy
with constant probablity in a linearly growing number of cost function calls.
This is in stark contrast to our classical greedy algorithm (FIG. 5) and QAOA
(FIG. 6) results, where the achievable polynomial time approximation ratio
steadily worsens with increasing $N_{C}/N$, and supports the theoretical
analysis of section IV. These results were obtained using the Fujitsu Quantum
Simulator, a classical HPC system.

###  VI.2 Trial minimum annealing performance

We also tested the TMA formulation of spectral folding–the formulation for
which we can make the most reliable analytical predictions–as plotted in Fig.
8. For these variations we used runtimes $t_{f}=N/12$ and $dt=0.025$; note
that this choice of $t_{f}$ is a factor of $8/3$ larger than in our QAOA
simulations but with the same scaling. In smaller system studies similar
qualitative performance was observed for shorter $t_{f}$ (such as $N/24$ or
$N/32$). For this variation we used a 3-XORSAT $H_{L}$–the formulation for
which we could predict performance in Sec. IV–with the minimum energy set to
$-2N$ via $C\left(t\right)$, which was linearly ramped down to zero by
$t_{f}$, and simple sinusoidal ramp profiles with $t_{r}=N/24$; the transverse
field strength $\kappa$ during the main evolution was 1.3. The careful reader
may note that choosing $C\left(t\right)$ to set the minimum energy to $-2N$
instead of $-N$ is naively suboptimal, as all other things being equal
sweeping over a larger energy range increases $W$ in Eq. (29), and should
reduce the returned probabilities $P_{q}\left(N\right)$ by an appropriate
prefactor. However in our simulations this choice consistently improved both
the prefactors and scaling, e.g. the value of $q_{a}$, as compared to choosing
a minimum energy $-N$ for $H_{L}$. We suspect this has to do with the band
structure considerations described in Appendix B but due to the complexity of
the problem, are unable to make a quantitative prediction.

The performance of the two approaches is qualitatively similar with subtle
differences as we vary the returned approximation ratio $q$. At lower
approximation ratios the AQC formulation returns higher probabilities, at
lower total gate count since there are no ramping steps and no additional
gates associated with adding $H_{L}$. However, at higher approximation ratios
the TMA formulation appears to be better able to approach the approximation
target $A=0.85$. As discussed in the algorithm definition, folded optimization
will definitionally fail to consistently return energies significantly below
$AE_{\mathrm{GS}}$, and we expect it to break down as $A$ gets too close to 1
given that QAOA and similar methods fail to reliably approximate these
problems. Choosing the best value for $A$ is thus a subtle issue that depends
on the problem class; for extensions of this method to hard CSPs it will
necessarily change from one problem class to the next.

Likewise, as mentioned in section IV.5, one can replace the random 3-XORSAT
lowering Hamiltonian in TMA for a simpler set of linear $Z$ biases, reducing
the gate count per timestep and, potentially, increasing the per-state
tunneling rate by implementing a shallower cost-per-flip curve. In comparison
to the 3-XORSAT variation discussed in the previous paragraph, to achieve good
performance we needed to double the ramp time. This protocol seemed to be more
sensitive to performance degradation from heating during ramps. The algorithm
also benefitted from adjusting $C\left(t\right)$ so that the minimum energy of
$H_{L}$ was $-3N$, as compared to $-2N$ for the 3-XORSAT $H_{L}$. Relative
performance for $A=0.75$ is comparable to the other variations, as illustrated
in Fig. 9, though the individual $P_{q}\left(N\right)$ show more significant
non-monotonicity that makes reliable curve fitting challenging. This issue is
even more pronounced for $A=0.85$ (data not shown), to the point that we did
not quote $q_{a}$ values for that variation in the Table 1.

Figure 8: Performance of the trial minimum annealing formulation of spectral
folding with a 3-XORSAT lowering Hamiltonian $H_{L}$, with $A=0.85$, plotted
for constraint densities (clockwise from top left)
$N_{C}=\left\\{2N,4N,6N,1.5\times N^{3/2}\right\\}$ and approximation ratios
between $q=0.5$ and $0.85$ with $N$ running from 10 to 24, for the parameters
detailed in the text. All data is derived from averaging over 960 random
instances and choices of $t_{f}$. Compared to the AQC formulation shown in
figure 7, the achievable approximation ratio $q_{a}$ is often slightly higher,
though the total gate count in this formulation is larger by a constant
prefactor. In all cases $q_{a}$ well exceeds the value of 0.6 conservatively
predicted for this formulation.

Figure 9: Performance of the trial minimum annealing formulation of spectral
folding with a local $Z$ bias $H_{L}$, with $A=0.75$, for $N$ running from 8
to 25. $q_{a}$ in each case is comparable to other variations, though greater
non-monotonicity in the individual $P_{q}\left(N\right)$ curves makes fitting
more difficult.

###  VI.3 Discussion of our spectrally folded optimization results

The reader may note that all of these simulations targeted energies in the
range $A=0.75$ to $A=0.85$, well above the theoretical prediction of $\sim
0.6$. We also tested $A=0.65$ at smaller scales (data not shown) for both
formulations but found no region of the parameter space where it showed
meaningful improvements over choosing $A=0.75$ or higher. When increasing $A$
to 0.85 we were not able to meet this target with constant probability, though
for some cases we did see improvements in the returned achievable
approximation ratio $q_{a}$ when compared to $A=0.75$. We did not run tests
with $A>0.85$ due to system size constraints, though comparing performance
between $A=0.75$ and $A=0.85$ in table 1 suggests we have found the
performance ceiling for these methods. Given the comparatively small system
sizes available to classical simulation our ability to draw meaningful scaling
distinctions with very small changes in $A$ is limited. Considering the
results of all our numerical simulations, the minimum achievable approximation
ratio for our PPSP problem classes is at least 0.7 (often 0.75) for the
quadratic AQC variation, and at least 0.75 (often 0.8) for the 3-XORSAT TMA
variation.

As summarized in table 1, the contrast between the clear super-polynomial
decay of higher approximation ratios with classical methods and QAOA (Figs.
4-6), and the constant probabilities returned by folded quantum optimization
(Figs. 7-9) is stark. In particular, for the established methods, the
polynomial time approximation threshold is set by the problem structure,
specifically the relative energy of the high energy clustering phase, and as
our PPSPs become more extremal this threshold steadily worsens as a fraction
of the energy of the planted ground state. We found no evidence for an
exponential separation in approximation power between local classical searches
and QAOA/AQC. At lower constraint densities, the clustering energy is not far
above the ground state and spectrally folded quantum optimization provides no
benefits for approximation, though it can still show scaling advantages for
returning states close in Hamming distance to the ground state (see
appendices). At higher constraint densities however, spectral folding is able
to return states close to $G$, in both energy and Hamming distance,
consistently and without degradation as $N_{C}/N$ increases. The performance
of all spectral folding variations tested was broadly similar, with the
3-XORSAT TMA variation returning the highest approximation ratio but at the
highest prefactor cost in gate count.

This illustrates the fundamentally different structure of collective tunneling
in this approach. Where high depth QAOA, AQC and other direct methods attempt
to find the ground state, asymptotically fail (given exponential scaling), and
return approximate states passed after this missed opportunity, spectrally
folded optimization deforms the cost function to search for states in an
exponentially large hyperspherical shell, avoiding the interference and
weakening field issues that are qualitatively responsible for QAOA’s lack of
obvious quantum advantage. In the interest of fair comparisons the total
runtime per-shot of all routines is $O\left(N\right)$ Hamiltonian calls. With
this simple linear scaling we are able to ensure that any exponential scaling
of the time to solution not visible from our figures or fits must involve very
small exponents, though of course we cannot rule out such behavior on numerics
alone. We thus do not see our numerical simulations as a precision scaling
benchmark of folded quantum optimization, but rather a test of the basic
veracity of our conservative analytical predictions in section IV. The results
of our simulations, for a wide range of parameters, all suggest that our
theory is sound.

###  VI.4 Extensions of these results

Having thoroughly explored MAX-3-XORSAT in this work, it would be interesting
to test spectral folding methods on other hard CSPs or problems that can be
straightforwardly formulated as such. We expect that the derivation of the
achievable approximation ratio through the overlap of dressed states could
generalize with some modifications to many other problem classes. We focused
on MAX-3-XORSAT due to the simplicity of its structure, classical
approximation hardness, and the fact that its exponential difficulty scaling
is obvious at small $N$ for more standard classical and quantum approaches.

It also strikes us as noteworthy that while spectral folding can be
implemented in traditional classical heuristics, there is no benefit to doing
so. As we discussed at length in Sec. II, classical algorithms generally start
from high energy states and attempt to cool towards the ground state through
local updates. In the case of linear spectral folding, the folding procedure
makes no difference whatsoever to the returned energy until the fold energy
has been reached, and in the worst case reaching that energy is an
exponentially difficult task for classical computers unless P=NP, due to
entropic barriers as discussed in Sec. II and confirmed in our simulations.
Changes to the high energy spectrum from quadratic folding are similarly not
expected to make approximation easier. So while spectral folding is not in and
of itself a quantum operation, we expect that it is only valuable in quantum
algorithms and we consider it an irreducibly quantum method as a result.
Further, because of its nonlocal nature, even in the locally gapped dressed
problem phase we expect volume-scaling entanglement in low-lying states of the
folded quantum spin glass, and significant classical simulation difficulty.

One interesting potential exception to this argument is quantum Monte Carlo.
Being a stoquastic problem, a folded quantum spin glass can in principle be
efficiently simulated using QMC Isakov et al. (2016); Andriyash and Amin
(2017); Jiang et al. (2017a, b); King et al. (2019), which for a uniform field
has in some cases proven to be an effective quantum-inspired classical solver
Heim et al. (2015); Albash and Lidar (2018). Incorporating spectral folding
into these algorithms is possible, though the loss of locality makes
evaluating each update much more expensive. As QMC is, fundamentally, an
energy-based classical local update rule combined with many replicas, we do
not expect it to overcome the entropic barrier in MAX-3-XORSAT the way that
true quantum evolution can.777As remarked earlier, we found the scaling of
linear spectral folding in the AQC framework, which has no influence on the
high energy spectrum at all, to be close to that of quadratic folding for the
same schedule. As finding the deep minimum is an entropic barrier challenge
for classical algorithms, we do not expect simulated quantum annealing through
QMC to be able to replicate the true quantum result. All that said, there may
be some narrow cases where incorporating spectral folding in a QMC calculation
could prove useful as a classical solver; this would be an interesting avenue
for future research. If it turns out that many or all of these instances are
efficiently approximable classically through QMC with a folded spectrum, that
would be a very significant discovery in its own right.

##  VII Conclusions and Outlook

Using the NP-hard MAX-3-XORSAT problem class, we explored the question of
classical and quantum approximation hardness from a practical, mechanism-
focused point of view. Guided by theoretical intuition, we proposed a class of
instances called planted partial solution problems (PPSPs) which we showed are
empirically hard for both exact and approximate optimization for classical
searches and established quantum methods such as the quantum approximate
optimization algorithm (QAOA) and adiabatic quantum computing (AQC). Through
extrapolation and qualitative analysis of a rigorous calculation of the
typical-case minimum gap at the first order paramagnet to spin glass
transition bottlenecking MAX-3-XORSAT, we were able to identify two effects
that significantly impede high-depth QAOA’s ability to approximate the hardest
instances, namely weakening transverse field and destructive interference. We
then proposed a novel algorithmic update, called _spectral folding_ , that
does not suffer from these two issues. Spectral folding is conceptually
simple; the analytical analysis of its performance is perhaps less so. And
rather intriguingly, it works by applying classical modifications to the
classical cost function being optimized, which only provide benefits for
quantum optimization.

Using a resummed extensive-order perturbation theory, we were able to predict
a constant fraction approximation guarantee for our difficult random
hypergraph PPSPs, and consequently, an exponential quantum speedup in the
classically hard regime. To further support our claims, we performed a series
of numerical simulations of high-depth QAOA, classical optimization routines,
and spectral folding methods, out to the largest sizes that we could feasibly
reach while still being able to gather good statistics. These numerical
results support our claims and we did not discover any meaningful
discrepancies or red flags in them; every major prediction here has been
numerically checked to the extent we found reasonably possible with current
supercomputing resources.

The implications of an effective fast approximation guarantee through spectral
folding–or any quantum algorithm, for that matter–are profound, and even
restricting ourselves to MAX-3-XORSAT it is important to ask what types of
hypergraph might cause it to fail. To be clear, “failure” in this case means
that, for a class of hypergraphs, some property invalidates the analytical
predictions we made and restricts polynomial-time approximation to values of
$A$ near zero. If a subsequent, more sophisticated analysis shows, for
example, that due to some subtle effect missed in our resummed extensive order
perturbation theory, the variations we propose are asymptotically limited to
$A=0.4$ for random hypergraphs instead of our prediction of $\geq 0.6$, we
would not consider that a general failure of the algorithm. Given the
exponential density of target states, for folded quantum optimization to fail
the per-state tunneling rate needs to be much worse than what our calculations
and simulations return; even doubling the per-state tunneling decay exponent
in Eq. (IV.4) yields $A\simeq 0.18$, a much worse approximation ratio than we
predict and observe in simulations, but still a constant fraction better than
random guessing.

We are not so hubristic as to claim it is impossible that there is some hidden
effect that substantially worsens scaling at large $N$ for random hypergraphs,
which is not captured by our theory and is invisible in our numerics. But we
see no evidence for it, have no idea what such an issue could be, and
formulated our theory such that the simplifying approximations we made were
more likely to underestimate the achievable approximation ratio $A$ than
overestimate it. Excellent empirical performance in simulation supports this
interpretation. We see it as much more likely that one can structure PPSP
hypergraphs in some non-random way as to violate the core assumptions of our
calculations and become inapproximable. The conditions on such hypergraphs
are, however, fairly strict; besides needing to more than double the
exponential decay rate of per-state tunneling as compared to random
hypergraphs, whatever property is responsible for the slowdown must be
resilient both to the addition of a random, uncorrelated problem with similar
constraint density and ground state energy, and to spectral deformations such
as quadratic folding. Further, such graphs must be relatively dense, as sparse
problems are easy to approximate by solving random sub-problems, and their
construction rules must not inadvertently render them amenable to classical
optimization.

Identifying hypergraph properties that cause this entire method to break down
could lead to valuable new discoveries about macroscopic quantum tunneling
physics and problem hardness, and further algorithm innovations in the steps
needed to mitigate their effects. As we stated in the introduction, the
underlying reasons for classical approximation hardness are generic and
intuitive if often only applicable to extremal problem instances, but their
quantum equivalents are not, and much more opaque. The correctness of our
predictions here would imply that quantum approximation hardness may not be
generic at all, and instead specific to as-yet undiscovered sets of problem
properties. Let’s go exploring.

##  VIII Acknowledgments

We would like to thank Matthew Jones, Chris Laumann, Gianni Mossi, Eleanor
Rieffel, Antonello Scardicchio and Davide Venturelli for valuable discussions
of the issues in this work. We would like to thank Caleb Rotello for detailed
feedback on the manuscript. We would also like to thank Takuto Komatsuki and
Joey Liu for support with HPC calculations. This work was supported by the
DARPA Reversible Quantum Machine Learning and Simulation program under
contract HR00112190068, as well as by National Science Foundation grants
PHY-1653820, PHY-2210566, DGE-2125899, and by the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers,
Superconducting Quantum Materials and Systems Center (SQMS) under contract
number DE-AC02-07CH11359. The SQMS Center supported EK’s advisory role in this
project, as well as time improving and fine tuning the algorithms and writing
the paper. Many of the numerical simulations in this work were performed with
a generous grant of HPC access from the Fujitsu Corporation. Part of this
research was performed while the one of the authors (BAB) was visiting the
Institute for Pure and Applied Mathematics (IPAM), which is supported by the
National Science Foundation (Grant No. DMS-1925919). The Flatiron Institute is
a division of the Simons Foundation.

##  Appendix A Tunneling between two $p$-spin wells

To benchmark our prediction in Eq. IV.3, we performed a series of numerical
simulations, shown in FIG. 10. $\left(p,\kappa\right)$ values used were
$\left\\{\left(2,1.25\right),\left(3,1.25\right),\left(4,1.2\right),\left(5,1.1\right),\left(7,1.06\right),\left(9,1.03\right)\right\\}$,
and the resulting exponents found from numerical fitting of
$\Omega_{0}\left(N\right)=a\sqrt{N}2^{-bN}$, as compared to Eq. IV.3, are
(listed as $\left(p,b_{num},b_{th}\right)$)
$\\{\left(2,0.175,0.180\right),\left(3,0.235,0.239\right),\left(4,0.275,0.287\right)$,
$\left(5,0.482,0.512\right),\left(7,0.591,0.623\right),\left(9,0.676,0.722\right)\\}$,
demonstrating excellent overall agreement, even with the prefactor, as shown
in FIG 10.

While fully sufficient for our purposes here (where $p=3$), we want to
highlight two issues with the formulation in Eq. IV.3 that suggest a more
refined treatment should be developed to tackle $p>3$ and/or asymmetric
minima. First, $\partial_{m}E_{m,n}$ is guaranteed to vanish at some point,
since the energy curve at the peak of the barrier is smooth, and since it
enters the calculation in the denominator, it predicts a diverging scale for
energy corrections. In symmetric cases, this divergence is canceled by the
factor of $N/2-2n$ in the numerator and the result is finite, but if the
energy curve is not symmetrical between the two minima those two factors won’t
generally coincide (as the barrier peak won’t occur at $M/2$ flips for
asymmetric minima) and the divergence will not be canceled.

Second, even we consider the symmetric case and the divergence _is_ canceled,
this formulation of perturbation theory can still lead to unphysical results.
For $p\geq 4$, while the energy corrections remain finite throughout, for
$\kappa$ relatively close to $\kappa_{c}$, the renormalized energy barrier
peaks at some $m=m_{c}<M/2$, and decreases from there. This is unphysical, and
we can mitigate it with the ad hoc prescription that $u_{m,0}$ simply stops
increasing beyond $m=m_{c}$ and maintains that value until $m=M/2$. The sum in
Eq. 23 is similarly modified so that the insertion of an additional spin flip
cannot lower the energy. Using this prescription produces quantitatively good
results for $p$ running from 4 to 9, as shown in FIG. 10.

Figure 10: Comparison of Eq. IV.3 to the tunnel splitting found from exact
diagonalization of the Hamiltonian in Eq. IV.3, for
$p=\left\\{2,3,4,5,7,9\right\\}$ (blue, gold, green, purple, red, brown) and
$\kappa$ stated in the text. There are no free parameters in this; exponents
found by extrapolating Eq. IV.3 out to $N=80$ agree with the numerical results
to within few percent. To make sensible predictions for $p\geq 4$ an
additional modification was needed to regularize the transverse field
corrections to state energies– see text for details. Only $p=3$ is directly
relevant to this work; simulation of other $p$ values is to demonstrate the
accuracy and flexibility of our theory.

One can also work through this analysis for the splitting $p=2$ mean field
all-to-all ferromagnet, where the transition out of the ferromagnet phase
occurs when $\kappa=2$ for this normalization. In that case all $N$ spins must
flip and there are no secondary spin corrections; further, due to a symmetry
cancellation the self energy corrections (e.g. Eq. IV.3) are energy-
independent and thus do not renormalize the tunneling barrier. Repeating the
same steps for this somewhat simpler calculation yields

$\Omega_{0}\left(N\kappa\right)\simeq\sqrt{\frac{N}{2\pi}}\frac{\kappa^{N}}{w\left(N,\kappa\right)^{2}}\left(\frac{e}{4}\right)^{N}.$
(36)

We of course still need to evaluate the normalization factor
$w\left(N,\kappa\right)$. Explicitly, it is given by

$\displaystyle
w\left(N,\kappa\right)=\sqrt{1+\sum_{k=1}^{N/2}\binom{N}{k}\kappa^{2k}\left(k!\prod_{j=1}^{k}\epsilon_{k}^{-1}\right)^{2}}$
(37)