Simulated Quantum Annealing Can Be Exponentially Faster than
Classical Simulated Annealing
Elizabeth Crosson∗

Aram W. Harrow†

arXiv:1601.03030v2 [quant-ph] 23 Jun 2016

June 24, 2016

Abstract
Can quantum computers solve optimization problems much more quickly than classical computers? One major piece of evidence for this proposition has been the fact that Quantum
Annealing (QA, also known as adiabatic optimization) finds the minimum of some cost functions exponentially more quickly than Simulated Annealing (SA), which is arguably a classical
analogue of QA.
One such cost function is the simple “Hamming weight with a spike” function in which
the input is an n-bit string and the objective function is simply the Hamming weight, plus
a tall thin barrier centered around Hamming weight n/4. While this problem can be solved
by inspection, it is also a plausible toy model of the sort of local minima that arise in realworld optimization problems. It was shown by Farhi, Goldstone and Gutmann [16] that for this
example SA takes exponential time and QA takes polynomial time, and the same result was
generalized by Reichardt [28] to include barriers with width nζ and height nα for ζ + α ≤ 1/2.
This advantage could be explained in terms of quantum-mechanical “tunneling.”
Our work considers a classical algorithm known as Simulated Quantum Annealing (SQA)
which relates certain quantum systems to classical Markov chains. By proving that these chains
mix rapidly, we show that SQA runs in polynomial time on the Hamming weight with spike
problem in much of the parameter regime where QA achieves exponential advantage over SA.
While our analysis only covers this toy model, it can be seen as evidence against the prospect
of exponential quantum speedup using tunneling.
Our technical contributions include extending the canonical path method for analyzing
Markov chains to cover the case when not all vertices can be connected by low-congestion
paths. We also develop methods for taking advantage of warm starts and for relating the quantum state in QA to the probability distribution in SQA. These techniques may be of use in
future studies of SQA or of rapidly mixing Markov chains in general.

1

Introduction

Classical algorithms are often useful but not provably so, with justifications for their success
coming from a combination of empirical and heuristic evidence. For example, the simplex
algorithm for linear programming was successful for decades before being proven to run in
polynomial time, and for a long time was the most practical LP solver even while the ellipsoid
algorithm was the only provably poly-time solver. Another example is MCMC (Markov chain
∗
†

Caltech IQIM. crosson@caltech.edu
MIT CTP. aram@mit.edu

1

Monte Carlo) which is used for applications in statistics, simulation, optimization and elsewhere,
but almost never in regimes that are covered by formal proofs of correctness.
With quantum algorithms, there has been necessarily a greater emphasis on provable correctness. The present state of quantum computing technology does not yet allow us to test
large-scale quantum algorithms empirically, nor can we usually empirically determine whether
a proposed quantum algorithm outperforms all classical algorithms on worst-case inputs. Nevertheless, heuristic quantum algorithms are likely to be important for practical problems, just
as they have been throughout the history of classical computing.
A particularly compelling heuristic proposal for optimization problems is quantum annealing
(QA), also known as quantum adiabatic optimization [22, 17]. (In this work we use the term
“quantum annealing” to mean adiabatic optimization in thermal equilibrium at a low but nonzero temperature, though in some other contexts QA may be taken to include non-equilibrium
thermal effects.) The
P idea of QA is to interpolate between a static problem-independent Hamiltonian such as − i σxi for which we can efficiently prepare the ground state, and a final
Hamiltonian whose ground state yields the desired answer. If we want to minimize a function f : {0, 1}n → R then we can take this final Hamiltonian to be proportional to diag(f ). This
can be thought of as a quantum version of classical simulated annealing (SA) with the diagonal
terms playing the role of bias and the off-diagonal terms causing hopping. Like SA its performance is hard to make provable general statements about, but it is a promising general-purpose
heuristic, and rigorous statements about its performance are known for many illustrative cases.
Intriguingly, QA has been shown to have an exponential asymptotic advantage over simulated annealing for certain cost functions [16]. Two examples are given in [16]: one in which the
quantum algorithm could be said to be taking advantage of symmetry (“the bush of implications”), and another which models tunneling (“the spike”) that will be our primary focus here.
The cost function for the spike is,

|z| + nα : n/4 − nζ /2 < |z| < n/4 + nζ /2
,
(1)
f (z) :=
|z|
: o.w.
where |z| is the Hamming weight of the string z, and α > 0, ζ ≥ 0 are independent of n. The
global minimum of f is the string with |z| = 0, but the spike term creates a local minimum at
|z| = n/4 + nζ /2. The spike presents a problem for a simulated annealing algorithm which only
proposes moves that flip k bits for k ≤ nζ . First recall the definition of SA: starting with a
point x ∈ {0, 1}n , it repeatedly chooses a random nearby point y (say within the Hamming ball
of radius k) and moves to y with probability min(1, e(f (x)−f (y))/T ), where T is a temperature
parameter that is gradually lowered. Following [16], consider first a SA algorithm that flips
one bit at a time, i.e. with k = 1. Since SA begins at high temperature the initial state is
overwhelmingly likely to have Hamming weight near n/2, and as the temperature of the system
is lowered the random walk will move to strings of lower Hamming weight until reaching the
local minimum at n/4 + nζ /2. This will happen for T = O(1), so at this point the probability of
ζ
accepting a move onto the spike with probability e−Ω(n ) , so classical SA requires exponential
time to find the global minimum with high probability. This argument applies to flipping any
k < nζ bits at once. Now suppose that the SA algorithm flips nc bits at a time for some c ≤ 1.
Once the Hamming weight is ≈ n/4, flipping nc random bits will change the√Hamming weight
by a random variable with expectation ≈ 21 nc and standard deviation ≈ 23 nc/2 . This has
c
probability e−Ω(n ) probability of being negative, so with high probability the SA algorithm will
not even attempt to move past the spike. In contrast, for any α + ζ < 1/2 it can be shown
that QA finds the global minimum with high probability in time O(n) [28], showing that an
exponential separation in the perfomance of SA and QA is possible.
While the spike is clearly a toy problem and can be solved efficiently by classical algorithms
that exploit its structure, an important aspect of both QA and SA is that a single, general
implementation of these algorithms is meant to be useful for solving a large variety of different

2

problems without knowledge of their structure. Moreover, the spike arguably demonstrates
a general advantage of QA over SA in tunneling through thin, high barriers in the energy
landscape.
On the other hand, the standard formulation of QA uses a stoquastic Hamiltonian (i.e. a
local Hamiltonian with non-positive off-diagonal matrix elements in the computational basis),
and computational models based on ground states or thermal states of such systems are believed
to be less powerful than universal quantum computation. In addition to complexity theoretic
evidence [7, 6], suggestive evidence for this belief is also provided by the quantum-to-classical
mapping of Suzuki et al. [31, 30], which allows for properties of low-energy states of a stoquastic Hamiltonian to be estimated using classical Markov chain Monte Carlo methods. These
algorithms are known as Quantum Monte Carlo (QMC) methods, and despite the name, are
algorithms for classical computers. While QMC for stoquastic Hamiltonians is always a welldefined algorithm, its performance depends on the rate at which a Markov chain converges to its
stationary distribution. This can range from polynomial to exponential time, and few general
conditions are known in which it is provably polynomial-time. A few cases where the simulation
can be made provably efficient are adiabatic evolution with frustration-free stoquastic Hamiltonians with a unique ground state [8] and ferromagnetic transverse Ising models in a large range
of temperatures [5], but while these have some physics significance, they do not translate into
nontrivial cost functions for QA.
When QMC is applied to QA Hamiltonians the result is an algorithm called simulated
quantum annealing (SQA). Although there are examples for which standard versions of SQA
take exponentially longer than the quantum evolution being simulated [18], the general challenge
from SQA to QA remains: for any purported speedup of QA we should see whether it can also be
achieved by SQA. Moreover, since SQA is a Markov chain based algorithm on a domain that can
be interpreted as a classical spin system, and since SQA is designed to sample from the output
of a quantum optimization procedure, SQA can be considered as yet-another physics-inspired
classical optimization method in its own right, which can naturally be compared to SA.
The main result of this paper is that the standard version of SQA, which does not use any
structure of the problem, finds the minimum of the cost function (1) in polynomial-time when
α + ζ < 1/2.
Theorem 1. Simulated quantum annealing based on the path-integral Monte Carlo method
efficiently samples the output distribution of QA for the spike cost function (1) when α+ζ < 1/2.
The running time using single qubit worldline updates is Õ(n7 ), and the running time using
single-site spin flips is Õ(n17 ). (Worldline and single-site spin flips are defined in Section 2.2.)
Thus SQA obtains an exponential speedup over SA for this particular problem. This result
suggests that the benefit of adiabatic evolution in tunneling through barriers should not be
thought of as an exclusively quantum advantage, since it can also be achieved by a generalpurpose classical optimization algorithm. We remark that the Õ(n17 ) scaling is almost certainly
an artifact of our analysis and that numerical evidence suggests that O(n4 ) single-site updates
should be sufficient [10].

Previous Work
There have been many past studies comparing the performance of SA and SQA using numerics [26, 1, 19] and more recently using analytical methods of physics such as the instanton
approximation to tunneling [20]. Studies comparing QA to SQA have also begun to emerge
since [2] found the success probabilities of SQA are highly correlated with the results of QA
performed on D-Wave quantum hardware with hundreds of qubits, while the distribution of
success probabilities for SA on the same set of instances bears little resemblance to that of QA
and SQA. More recently, the performance of QA, SQA, and SA was empirically compared on
an ensemble of spin glass instances with were designed to have tall, thin barriers [12], as a step

3

acronym

definition

notes

SA

simulated
annealing

Classical MCMC algorithm that simulates thermal annealing [23]

QA

quantum
annealing

Adiabatic optimization [17] generalized to low temperatures

SQA

simulated
quantum
annealing

Classical MCMC that simulates quantum annealing; cf. Section 2.2

Table 1: Our paper repeatedly uses the similar sounding acronyms SA, QA and SQA, which are
defined here. QA is a quantum algorithm while SA and SQA run on classical computers.

towards understanding the kinds of instances for which QA has an advantage over SA. In that
work QA and SQA were found to have roughly the same scaling with system size for that particular ensemble of instances, though it was also pointed out that the large constant overhead in
SQA made it less competitive in the sense of wall-clock times using modern classical hardware.
Without access to quantum hardware, comparison of SQA and QA is either limited to small
system sizes where QA Hamiltonians can be exactly diagonalized (. 50 qubits), or to models
for which analytical solutions of the quantum system are known (such as the spike problem
we study here). We remark that the spike and related objective functions have the subject of
recent analytic work [24, 4], and that there have also been numerical studies of SQA [10, 3],
with findings that are consistent with our main result.

Proof Outline
Our proof of the efficient convergence of SQA on the spike problem involves bounding the
mixing time of the underlying Markov chain, and there is an interesting parallel between a
method which was used to lower bound the QA spectral gap when α + ζ < 1/2 [28]. There,
a lower bound on the quantum gap can be found using a variational method with a trial wave
function equal to the ground state of the system when no spike term is present (i.e. QA for
the spikeless Hamming weight cost function f˜(z) = |z|). Similarly, we compare the spectral
gap λ of the SQA Markov chain for the spike system with the spectral gap λ̃ of the spikeless
system (throughout the subsequent sections we use tildes to distinguish quantities belonging to
the spikeless system). Without a spike term, the quantum Hamiltonian H̃ is a tensor product
operator with no interactions between the qubits. This trivial system translates in SQA to a
collection of n non-interacting 1D classical ferromagnetic Ising models in a uniform magnetic
field (which will become clear when the SQA Markov chain is described in detail in Section 2.2),
and upper bounding the mixing time for this system is relatively straightforward.
Let π and π̃ be the stationary distribution of the SQA Markov chain with and without the
spike. These stationary distributions are close in a sense, kπ − π̃k1 < poly(n−1 ), but on the other
hand there are exponentially many points x ∈ Ω for which the ratio π(x)/π̃(x) is exponentially
small. A review of existing comparison techniques in Section 3.1 concludes that none is quite
suited to the present problem; indeed the review [13] states that there have been “relatively few
successes in comparing chains with very different stationary distributions”. To overcome this
we introduce a comparison method which involves partitioning the state space into “good” and
“bad” sets of vertices, Ω = ΩG ∪ ΩB . In Section 3.2 we begin with a set of canonical paths
yielding a bound ρ̃ on the congestion of the easy-to-analyze chain, and show that the paths

4

which lie entirely within ΩG can be used to construct an upper bound on the congestion ρ of
the difficult-to-analyze chain, albeit within the set ΩG of measure less than 1. This “most-paths
comparison” method may be of independent interest, and so the exposition in Section 3.2 is
given without dependence on the specific details of SQA.
There are two main ingredients to the comparison method. First we show that if two chains
have stationary distributions and transition probabilities that are similar on most of the points,
then we can convert a known gap for one chain into a set of canonical paths for most of the
state space of the other.
Theorem 2 (Most-paths comparison). Let (π, P ) and (π̃, P̃ ) be reversible Markov chains with
the same state space graph (Ω, E). Let a = maxx∈Ω π(x)/π̃(x) and define Ωθ := {x ∈ Ω : π(x) <
θπ̃(x)}. If there is a set of canonical paths for (π̃, P̃ ) achieving congestion ρ̃ and satisfying
3a2 ρ̃π̃(Ωθ ) < 1, then there is a subset ΩG ⊂ Ω with π(ΩG ) ≥ 1 − 3a2 ρ̃ θπ(Ωθ ), and a canonical
flow for (π, P ) that connects every x, y ∈ ΩG with paths contained in ΩG for which the congestion
ρ of any edge in ΩG satisfies
#
"
P̃ (x, y) 2
a ρ̃.
(2)
ρ ≤ 16 θ max
x,y∈ΩG P (x, y)
There is a caveat here, which is that this bound on the congestion applies to the transitions
of the Markov chain P on the subset ΩG . If we assume that walkers leaving ΩG are deleted
then, since P is not restricted to ΩG , these transitions form a substochastic “leaky” random
walk on the set ΩG , with a quasi-stationary distribution equal to π within this subset. (The
term “quasi-stationary” refers to the fact that in the infinite time limit repeated applications of
P |ΩG will converge to zero, but there may be a long intermediate time when we are close to π.)
Thus it is necessary to show that the chain mixes before it leaves the good set ΩG . One
way to guarantee this is to use a “warm start,” meaning a starting sample from a distribution
that is close to the quasi-stationary distribution. Our analysis of SQA will rely on following the
adiabatic path used by the quantum algorithm to guarantee the warm-start condition.
Theorem 3. Let (π, Ω, P ) be a reversible Markov chain and suppose Ω = ΩG ∪ΩB is a partition.
Let PG be the substochastic transition matrix PG (x, y) := P (x, y)1x∈ΩG 1y∈ΩG . Suppose there is
a set of canonical paths connecting every pair of points x, y ∈ ΩG , and the congestion of the
walk P on this set of paths is ρ. If µ is a warm start with µ(x) ≤ M π(x) for all x ∈ ΩG then
the distribution obtained by starting from µ and applying t steps of the random walk satisfies
−1 −t/ρ
kµPGt − πk1 ≤ M tπ(ΩB ) + πmin
e

(3)

The SQA state space can be interpreted as a path (worldline) representation of the original
quantum system, and the bad states which constitute ΩB will be those for which the paths
spend too much “time” on the location of the spike (i.e. on strings with Hamming weight
between n/4 − nζ /2 and n/4 + nζ /2). States that spend too much time on the spike are those
for which π(x)/π̃(x) is exponentially small, and naturally those are the ones we will need to
exclude. In Section 4 we show that the mean spike time is proportional to the square of the
ground state amplitude on the spike, while the m-th moment of the spike time distribution can
also be bounded using the properties of the corresponding quantum system. Finally, we use
the derived upper bound on the m-th moment of the spike time distribution to upper bound
the probability of large deviations from the mean spike time, which yields an upper bound on
π(ΩB ) that suffices to complete the proof.

Discussion
Our proof does not bound the convergence time for SQA α + ζ > 1/2, although QA does work
for some values of (α, ζ) in this range, such as when ζ = 0, α = O(1) [24] or α + 2ζ < 1 [4]. We

5

conjecture that SQA will be efficient for these values as well (which is supported by numerical
evidence), though this will require extensions of the present techniques. The approach used in
this work is also suggestive of a more general connection between the quantum spectral gap of
Hamming symmetric barrier problems (including barriers of various shapes and widths) and the
corresponding performance of SQA, which we sketch here.
Assume we are near the critical value s∗ of the adiabatic parameter at which the system
tunnels through the barrier i.e. for s < s∗ the ground state probability mass is concentrated on
one side of the barrier, while for s > s∗ the opposite occurs, and for s = s∗ the probability mass
on both sides of the barrier is O(1)). While a general understanding of what barriers admit
such a tunneling description has not yet been found, there are many examples for which this is
known to occur[16, 24, 27]. Assume that the QA spectral gap at s∗ is ∆ = 1/ poly(n), so that
QA will be able to pass through the barrier efficiently.
We expect the first excited state to have a node at some location |ki inside the barrier region,
and some properties of the ground state wave function near |ki can be inferred from the spectral
gap. Since the spectral gap is at least 1/ poly(n), we know that the ground state wave function
cannot be too small inside of the barrier, or else there would be a balanced-weight cut in the
ground state that could be used to construct an orthogonal state with low energy. This benefits
the comparison approach because the ground state amplitudes inside the barrier need to be at
least 1/ poly(n) in order for a 1/ poly(n) fraction of the canonical paths to transfer from the
spikeless system to the system with a barrier.
On the other hand, the amplitudes inside the barrier cannot become too large because the
barrier is a classically forbidden region, and because an upper bound on ∆ will also upper
bound the amplitudes near |ki (by the same argument using a balanced-weight cut together
with the assumption that the first excited state has a node at |ki). The fact that the total
probability mass inside the barrier is not too large could be useful for deriving the necessary
upper bounds on π(Ωθ ) that are used in the comparison approach. The primary ingredient that
would be needed to make this sketch rigorous is a better understanding of the wave functions
and spectral gap of the quantum system with a barrier, and depending on the outcome of this
understanding the comparison approach for analyzing SQA may require modifications as well.
More generally, we believe that our most-paths comparison methods should have wider applicability. For example, consider a collection of classical particles with weak repulsive interactions.
If the particles were non-interacting the thermal distribution of the particles would be easy to
sample from, and if the interactions are weak enough, then they do not shift the typical probabilities (or energies) by very much. While some configurations will have exponentially lower
probability in the interacting case (if many particles are very close to each other), these configurations should be overall very unlikely. In this setting our framework should imply that the
repulsive interactions do not significantly worsen the mixing time.
Finally, while relatively few rigorous facts are known about the performance of SQA or
Quantum Monte Carlo more generally, it remains in practice a successful and widely used class
of algorithms. This strikes us as an area where theorists should work to catch up with current
practice.

Overview of remaining sections
In Section 2 we review the QA and SQA algorithms in the forms that our paper will use
them. Section 3 fleshes out our Theorems 2 and 3 on incomplete sets of canonical paths and
leaky random walks. This section presents its results for a general Markov chains and can be
understood without reference to the specific SQA Markov chains discussed elsewhere in the
paper. Our main result is proved in Section 4; this entails showing that the Markov chains from
Section 2 meet the mixing conditions laid out in Section 3.

6

2
2.1

Background
Quantum annealing

Quantum annealing associates a cost function f : {0, 1}n → R with a Hamiltonian that is
diagonal in the computational basis,
X
Hf :=
f (z)|zihz| ,
(4)
z∈{0,1}n

so that the ground state of Hf is a computational basis state corresponding to the bit string
that minimizes f . To prepare the ground state of Hf the system is initialized in the ground
state of a uniform transverse field, which can be easily prepared,
H0 := −

n
X

σix

,

i=1

1
|ψinit i := √
2n

X

|zi,

(5)

z∈{0,1}

and then linearly interpolates between H0 and Hf ,
H := H(s) = (1 − s) H0 + sHf ,

(6)

where the adiabatic parameter s sweeps through the interval 0 ≤ s ≤ 1. The total run time tmax
of the algorithm depends on how quickly the adiabatic parameter is adjusted, which defines a
time-dependent Hamiltonian H(t) := H(s = t/T ). At zero temperature the system evolves acd
cording to the Schrödinger equation, dt
|ψ(t)i = −iH(t)ψ(t), and the adiabatic theorem ensures
that the state ψ(T ) at the end of the evolution has a high overlap with the ground state of Hf
as long as T ≥ poly(n, ∆−1 ), where ∆ = mins E1 (s) − E0 (s) is the minimum gap between the
two lowest eigenvalues of H(s) during the evolution.
More generally (and realistically) we can take the state of the system to be not the ground
state but a thermal state with inverse temperature β < ∞. The equilibrium thermal state of
the system evolves with the adiabatic parameter,
σ(s) :=

e−βH(s)
Z(s)

Z(s) := tr e−βH(s) .

,

(7)

We can see that if β = ∞ the system will be in the ground state, and we will assume that β
is sufficiently large so that the system will remain close to the ground state. Just as adiabatic
theorem has an error term corresponding to transitions out of the ground state, in thermal
annealing the system will in general not be exactly in the equilibrium state (7). Nevertheless it
is a useful idealization, and one that our simulation algorithms will aim to reproduce. Specifically
the simulated quantum annealing algorithm described in the next section will produce samples
from the distribution Πs (z) := hz|σ(s)|zi. The minimum gap ∆ of the Hamiltonian (6) with the
cost function (1) is constant when α + ζ < 1/2 [28], and the density of states is such that taking

β = n for a constant  > 0 suffices to make k ρ(s) − |ψ0 (s)ihψ0 (s)| k1 < O(ne−n ). Thus this
low-temperature thermal-equilibrium version of QA produces a final state which can be sampled
to obtain the minimum of f .

2.2

Simulated quantum annealing

In principle stoquastic Hamiltonians such as (6) are amenable to a variety of classical Markov
chain based simulation algorithms, which are collectively known as quantum Monte Carlo
(QMC) methods (the term “stoquastic” is a combination of “quantum” + “stochastic” in the
sense of stochastic matrices [7]). Any QMC method applied to the QA Hamiltonian (6) defines

7

a version of SQA. The version we consider here is based on the path-integral representation of
the thermal state (7).
The starting point of the method is to express Z(s) as a sum over an exponential number
of nonnegative terms, or in physics language to write the quantum partition function as an
imaginary-time path integral over trajectories basis states. Since the Hamiltonian is stoquastic
in the standard basis, these trajectories will have the form (x1 , ..., xL ), where xi ∈ {0, 1}n ,
Z = tr e−βH = tr

L
Y

e−

βH
L

=

L
X Y

hxi |e−

βH
L

|xi+1 i,

(8)

x1 ,...,xL i=1

i=1

where xL+1 := x1 , L is to be chosen in the next step and we neglect the dependence on s to
keep the notation simple. The individual bit strings xi in (x1 , ..., xL ) are sometimes called“time
slices.” When βkHk/L is sufficiently small the Suzuki-Trotter approximation provides a way to
split up the non-commuting terms in e−βH/L while only incurring a small error in the partition
function. Define A = −βsHf and B = −β(1 − s)H0 , so that the partition function is Z =
tr eA+B . Define the Suzuki-Trotter approximation to the partition function,
h A B iL
Z := tr e L e L .

(9)

According to Lemma 3 and the surrounding discussion of [5], taking L = Θ
for sufficiently small δ > 0 achieves

p

(βkHk)3 δ −1

(1 − δ)Z ≤ Z ≤ Z(1 + δ).



(10)

We will take δ = n−1 , which implies L = O(n2 β 3/2 ) and subsequently ignore the error in (10)
because it does not affect the convergence time of SQA, and creates only a negligible error in
the distribution which SQA will sample from for the spike cost function (1). Expanding (9) as
was done for (8),
Z=

L
X Y

A

B

hxi |e L e L |xi+1 i

(11)

x1 ,...,xL i=1

=

βs

X

e− L

PL

i=1

f (xi )

x1 ,...,xL

=

B

hxk |e L |xk+1 i

(12)

k=1
βs

X

L
Y

e− L

PL

i=1

f (xi )

x1 ,...,xL

n Y
L
Y

x

hxj,k |eωσj |xj,k+1 i,

(13)

j=1 k=1

where ω := β(1 − s)/L. Using the identity exp (ωσ x ) = cosh(ω)I + sinh(ω)σ x , the individual
factors of the product in (13) become


x
hxj,k |eωσj |xj,k+1 i = cosh(ω) 1xj,k =xj,k+1 + tanh(ω)1xj,k 6=xj,k+1 ,
(14)
and after suppressing the uninteresting multiplicative factor of cosh(ω)nL , the partition function
is expressed as
Z=

X
x1 ,...,xL

βs

e− L

PL

i=1

f (xi )

n Y
L
Y


1xj,k =xj,k+1 + tanh(ω)1xj,k 6=xj,k+1 ,
j=1 k=1

8

(15)

and so Z can be viewed as the normalizing constant of a probability distribution,
π(x1 , ..., xL ) =

=

n,L

1 − βs PLi=1 f (xi ) Y 
e L
1xj,k =xj,k+1 + tanh (ω) 1xj,k 6=xj,k+1
Z

(16)

1 − βs PLi=1 f (xi )
e L
φ(x̄j )
Z
j=1

(17)

j,k=1
n
Y

where x̄j := (xj,1 , ..., xj,L ) is called “the worldline of the j-th qubit”, and φ(x̄j ) := tanh(ω)|{k:xj,k 6=xj,k+1 }|
counts the number of consecutive bits which disagree in that worldline.
Performing a similar calculation for Π(x) = hx|σ|xi (where σ = e−βH /Z) one finds it can be
expressed as the marginal of π on the first time slice,
X
Π(x) =
π(x, x2 , ..., xL )
(18)
x2 ,...,xL

From this point SQA proceeds by discretizing the adiabatic path and using the Markov
chain Monte Carlo method to sample from π at various values of the adiabatic parameter
s1 ≈ 0, . . . , smax ≈ 1. We will analyze two discrete-time Markov chains which both have
stationary distribution π on the state space Ω = {0, 1}n×L . The first chain consists of singlesite Metropolis updates. If x, x0 ∈ Ω differ by a single bit, the transition probability from x
to x0 is


π(x0 )
1
0
min 1,
,
(19)
PM (x, x ) =
2nL
π(x)
and otherwise the transition probability is zero. We interpret this as follows: with probability
1/2 propose to flip one of the nL bits in x to create x0 , and accept this move with probability
min(1, π 0 (x)/π(x)). Otherwise the configuration x remains unchanged.
Note that π is supported on all of Ω when s < 1, while at s = 1 the state space becomes
disconnected under the local move transitions described above. This is not a major limitation
for our application, however, since sampling from Π when smax = 1 − n1 suffices to find the true
minimum of f . In practical implementations it is important to avoid this “critical slowing down”
as s → 1 by using cluster updates that flip many spins at once. Since local moves in classical
SA correspond to flipping a single bit {1, . . . , n}, one could argue that a move which arbitrarily
updates the worldline of a single qubit should be considered a local move in SQA as well.
More importantly, such moves can be performed efficiently for any easily computable objective
function f . Therefore in addition to the single-site Metropolis updates defined above we will
analyze the single qubit heat-bath worldline updates, which are a form of generalized
heat-bath updates [14].
Definition. The heat-bath worldline update (x̄1 , . . . x̄n ) → (x̄01 , . . . , x̄0n ) proceeds as follows:
1. Select a site i ∈ {1, . . . , n} uniformly at random.
2. Set x̄0j = x̄j for all j 6= i.
3. Choose x̄0i from the conditional distribution π(x̄0i |x̄01 , . . . , x̄0i−1 , x̄0i+1 , . . . , x̄0n ).
As with all generalized heat-bath updates these transitions define a Markov chain which is
reversible with respect to π. An algorithm that efficiently implements these transitions is given
in [15] and without using any structure of the cost function (1) it runs in Õ(n2 β 2 ) elementary
steps with high probability. One practical advantage of these updates is that L can be taken to
be effectively ∞ (or more precisely, the runtime can be made to scale as poly log(L)). This is
because the typical number of bit flips per worldline is independent of L, and so we can represent
x̄i succinctly by specifying only the locations of the bit flips. This is a major improvement in

9

practice but for the purposes of establishing poly-time mixing, we will not explore this further
here.
At each value of the adiabatic parameter, the run-time of SQA will be determined by the
mixing time of either of the Markov chains described above. This quantity can be defined in
terms of the total variation distance from π to the distribution P t (x, ·) obtained by running the
chain for t steps starting from x,
1 X t
|P (x, x0 ) − π(x0 )|,
(20)
dx (t) := max |P t (x, A) − π(A)| =
A⊆Ω
2 0
x ∈Ω

with the mixing time τ () being the worst-case time needed to be within variation distance  of
the stationary distribution,
tmix () := max min{t : dx (t0 ) ≤  ∀t ≥ t0 }.
t

x∈Ω

(21)

A standard way to bound the mixing time is to relate it to the spectral gap λ of the transition
matrix P [25]. For all x ∈ Ω,
−1 −λt
kP t (x, ·) − πk1 ≤ πmin
e .
(22)


which implies tmix () ≤ λ−1 log π1min . These bounds are worst-case in the sense that they
handle any starting vertex x ∈ Ω. In the next section of our paper, we will develop slightly
different methods to deal with the fact that our Markov chain may not mix efficiently from some
starting vertices.

3

Incomplete sets of canonical paths

In this section we discuss how to show rapid mixing even in the presence of a small number
of bad vertices. While we will freely make assumptions specific to our particular problem we
will introduce some techniques that apply more generally to the analysis of Markov chains, and
may be of use elsewhere. The notion of “good” and “bad” sets in Markov chains has been used
before [21, 9], but always (to our knowledge) in a setting where a separate argument shows the
bad set can always be quickly escaped. By contrast we model the bad set quite pessimistically
and can assume that the walker gets absorbed (or equivalently, trapped for an exponential
amount of time) upon hitting a bad vertex. Despite this we will show that the overall algorithm
works with high probability.

3.1

Markov chains and the path comparison method

This subsection reviews some standard facts from section 13.5 of the book by Levin, Peres and
Wilmer [25]. Let (π, P ),(π̃, P̃ ) be reversible Markov chains
and define Q(x, y) := π(x)P (x, y)
P
and likewise for Q̃. Define the inner product hf, giπ := x π(x)f (x)g(x) and the Dirichlet form
E(f, g) := hEf, giπ := h(I − P )f, giπ .

(23)

Lemma 13.11 of [25] states that
E(f ) := E(f, f ) =

1 X
[f (x) − f (y)]2 Q(x, y).
2

(24)

x,y∈Ω

This can be used to define the gap (cf. Remark 13.13 of [25]),
λ :=

min

Ω

E(f ) =

f ∈R
f ⊥π 1,kf k2 =1

min

Ω

f ∈R
Varπ (f )6=0

10

E(f )
.
Varπ (f )

(25)

To estimate the gap we will use various comparison methods. Given x, y ∈ Ω let Pxy be the
set of simple paths from x to y, and suppose that νxy is a measure over this set. Then we can
define a congestion ratio:
ρ(e) :=

1 X
Q̃(x, y)
Q(e)
x,y∈Ω

X

νxy (Γ)|Γ|

(26)

Γ:e∈Γ∈Pxy

ρ := max ρ(e)

(27)

e∈E

This last sum ranges over all paths Γ ∈ Pxy that contain the edge e and thus is a measure of
the total load on edge e. Then Corollary 13.26 of [25] states that


π(x)
ρλ.
(28)
λ̃ ≤ max
x∈Ω π̃(x)
If νxy has all its weight on a single flow γxy then we can slightly simplify (26) to
ρ(e) =

1 X
Q̃(x, y)|γxy |.
Q(e)

(29)

x,y∈Ω
e∈γxy

Another variant is when we compare a chain Q̃ with the complete graph for which the x → y
transition probability is simply π̃(y). The paths used for this purpose are labelled {γ̃xy }. The
corresponding congestion is
ρ̃(e) :=

1 X
π̃(x)π̃(y)|γ̃xy |.
Q̃(e) x,y∈Ω

(30)

e∈γ̃xy

Since the complete graph has gap 1 and the same stationary distribution, we have
λ̃ ≥

3.2

1
.
ρ̃

(31)

Most-paths comparison

In this section we will describe a comparison method for two reversible Markov chains (π̃, P̃ ),
(π, P ) defined on the same state space graph (Ω, E). The method is designed for chains which
are comparable on a subset ΩG containing most of the stationary probability of both chains,
but which may have π and π̃ differ greatly in other regions of the state space. Assuming there
is a set of canonical paths connecting all x, y ∈ Ω with congestion ρ̃ for (π̃, P̃ ), then we will
show it is possible to construct a canonical flow for (π, P ) connecting all x, y ∈ ΩG which has
comparable congestion ρ within ΩG .
In our application of this method the easy-to-analyze chain (π̃, P̃ ) will be the SQA Markov
chain for the spikeless distribution, and (π, P ) will be the SQA chain for the spike system.
Therefore π(x) will never be much larger than π̃(x), and the quantity
a := max
x∈Ω

π(x)
π̃(x)

(32)

will be O(1) for our application. However, there will be exponentially many configurations x ∈ Ω
with π(x)/π̃(x) = O(e−n ) and so we define a subset parameterized by the value of this ratio,


π(x)
Ωθ := x ∈ Ω :
<θ .
(33)
π̃(x)

11

Let {γ̃xy } be the set of paths for (π̃, P̃ ) for which the congestion (defined in (30)) is ρ̃. These
paths may lead to far larger congestion for (π, P ) if they involve routing flow over edges e where
Q(e)  Q̃(e). However, our first claim is that most of these paths should avoid the bad set Ωθ .
First observe that for any e the definition of ρ̃ implies that
X
π̃(x)π̃(y) ≤ ρ̃Q̃(e).
(34)
x,y:γxy 3e

Define Eθ to be the set of edges incident upon Ωθ . Reversibility implies that
1
π̃(Ωθ ) ≤ Q̃(Eθ ) ≤ π̃(Ωθ ).
2

(35)

(The inequalities are because an edge may be counted once or twice depending on whether both
end points are in Ωθ .) Additionally for e = (v, w) 6∈ Eθ we have
Q̃(e) = π̃(v)P̃ (v, w) ≤ θπ(v)P̃ (v, w) ≤ θRπ(v)P (v, w) = θRQ(e),

(36)

where R := max(v,w)∈E
/ θ P̃ (v, w)/P (v, w). Note that in our application R will be O(1) for all
the types of transitions we consider.
Next we Let C := Ω2 , CB := {(x, y) ∈ C : γxy ∩ Eθ 6= ∅} and CG := C − CB . Now we sum
(34) over all e ∈ Eθ to obtain
X
π̃(x)π̃(y)|γ̃xy | ≤ ρ̃Q̃(Eθ ) = ρ̃π̃(Ωθ ).
(37)
(x,y)∈CB

We conclude that not many of the γ̃xy go through any edges in Eθ .
Now we define a second partition of Ω into good and bad vertices. Let CB (x) := {y : (x, y) ∈
CB } and define CG (x) similarly. Then define
ΩB := {x : π̃(CB (x)) ≥ 1/3a},

(38)

and ΩG = Ω − ΩB . From (37) we have π̃(ΩB ) ≤ 3aρ̃π̃(Ωθ ).
Using first (32) and the assumption that π(Ωθ ) is sufficiently small we have
π(ΩG ) = 1 − π(ΩB ) ≥ 1 − 3a2 ρ̃π̃(Ωθ ) ≥

11
.
12

(39)

Thus ΩG is a large-measure set that is mostly well connected. Indeed ∀x ∈ ΩG , π(CG (x)) ≥ 2/3.
We can now define canonical flows between all pairs x, y ∈ ΩG . Observe that
π(CG (x) ∩ CG (y) ∩ ΩG ) ≥ 1/4.

(40)

This means that even if x, y are not directly connected to each other by a path that avoids Eθ ,
they are still indirectly connected via pairs of paths that route through a large-measure (≥ 1/4)
set. We will see that this allows us to construct canonical flows that are not too much worse
than our original canonical paths.
Observe next that the conditional distribution πxy := π|CG (x)∩CG (y)∩ΩG satisfies πxy ≤ 4π.
We will construct our flow νxy by choosing a random r ∼ πxy and concatenating the paths γ̃xr

12

and γ̃ry . The load on edge e from these flows is 0 if e 6∈ Eθ , or if e ∈ Eθ can be bounded as
X
X
π(x)π(y)
πxy (r)(1e∈γ̃xr + 1e∈γ̃ry )(|γ̃xr | + |γ̃ry |)
(41)
x,y∈ΩG

r

X

≤4

π(x)π(y)

X

≤ 16

X
X

by symmetry (42)

X

π(x)

X

π(r)1e∈γ̃xr |γ̃xr |

since πxy ≤ 4π

(43)

r∈CG (x)∩ΩG

x∈ΩG

≤ 16

X

π(x)π(y)

x,y∈ΩG

= 16

πxy (r)1e∈γ̃xr |γ̃xr |

r

x,y∈ΩG

π(r)1e∈γ̃xr |γ̃xr |

since π is normalized (44)

r∈CG (x)∩ΩG

π(x)π(r)1e∈γ̃xr |γ̃xr |

passing to a superset (45)

(x,r)∈CG

≤ 16a2

X

π̃(x)π̃(r)1e∈γ̃xr |γ̃xr |

from (32) (46)

(x,r)∈CG

≤ 16a2 ρ̃Q̃(e)

from (34) (47)

2

≤ 16θRa ρ̃Q(e)

from (36) (48)

We conclude that our set of canonical flows achieves congestion
ρ ≤ 16θRa2 ρ̃,

(49)

albeit on a set of measure less than one, namely ΩG . In the next section we will discuss the
implications of this for mixing.

3.3

Leaky random walks

In this section we will prove Theorem 3 by analyzing the general case of a substochastic random
walk PG which is defined by the (unnormalized) restriction of the transition probabilities of a
Markov chain P to a subset ΩG of its full state space Ω. Let Ω = ΩG ∪ ΩB be a partition of Ω,
and define
PG (x, y) := P (x, y)1x∈ΩG 1y∈ΩG .
(50)
If P is ergodic then the stationary distribution π has support on ΩB , so limt→∞ PGt (x, ·) = 0
for any starting point x ∈ Ω. However there may be an intermediate range of times, tmix ≤ t ≤
tfail , for which kPGt (x, ·) − πk1 < δ for some desired δ, if the initial point x ∈ ΩG is a warm start
taken from x ∼ µ with kπ − µk1 sufficiently small. Define an M -warm distribution to be one
satisfying µ(x) ≤ M π(x) for all x. Let ΠG to be the projector onto RΩG , so that PG can be
expressed as
PG = ΠG P ΠG .
(51)
To bound the gap of PG observe that the path comparison method (Thm 13.23 of [25]) can be
applied even to substochastic matrices. Here it implies that
ΠG EΠG π ρ−1 ΠG .

(52)

with ρ from (49), and with π denoting the semidefinite ordering with respect to the h, iπ inner
product (i.e. A π 0 means hf, Af iπ ≥ 0). Eq. (52) does not directly give bounds on the gap
of P . Indeed we could in principle have P (x, x) = 1 for all x ∈ ΩB in which case P would have
many eigenvalues equal to 1. Thus without additional information we cannot say anything more
about P .

13

Instead we will analyze a slightly different Markov chain. We will define a “filled-in” chain
PF by adding nonnegative numbers to the entries of PG in such a way that PF is stochastic and
has π as a stationary distribution. This does not uniquely specify PF and indeed P also satisfies
these conditions. WePwill instead choose PF in a way that guarantees fast mixing. Specifically
with probability 1 − y∈Ω PG (x, y) we will forget our current location x and jump according to
some “fill-in” measure ϕ. For this to result in π being the stationary distribution we must have
ϕ = π(I − PG ). Defining the column vector 1 := (1, 1, . . . , 1)T , we have
PF = PG + 1ϕ

(53)

πPF = π(PG + 1ϕ) = π(PG + 1π(I − PG )) = π.

(54)

Now that PF is a proper Markov chainP
we will bound its Dirichlet form. Recall that in (25)
we can assume that f ⊥π 1 meaning that x π(x)f (x) = 0. Note as well that
EF = I − PF = I − PG − 1π(I − PG ) = (I − 1π)(I − PG ).

(55)

Define Dπ to have π along the diagonal and zeros elsewhere. Then
hEF f, f i = h(I − 1π)(I − PG )f, f i
T

= f (I −
T

= f (I −

PGT )(I −
PGT )Dπ f

T

(56a)

T

π 1 )Dπ f
since f ⊥π 1

T

T

T

−1 T

= f (I − Π)Dπ f + f (Π −
≥ f (I − Π)Dπ f + ρ
≥ρ

−1 T

=ρ

−1

(56b)

PGT )Dπ f

f ΠDπ f

f Dπ f
Varπ [f ]

(56c)
(56d)

using (52)

(56e)

since ρ ≥ 1

(56f)

since f ⊥π 1

(56g)

We conclude that PF mixes rapidly, while differing from PG only in the events where leakage
occurs. More precisely we can define a coupling between two processes: the first a walk evolving
according to PF and the second a walk in
P which x moves to y with probability PG (x, y) and
to a special state * with probability 1 − y PG (x, y). The state * is absorbing, meaning that
when the second walker is at state * it stays there. Conditioned on the second walker not being
in state * the two walkers will be at the same location. Additionally the probability of the
second walker ending in * equals the probability that the first walker ever passed through ΩB .
If we start with a probability distribution µ and take t steps then the probability of ending in
* is µ(PFt − PGt )1. This quantity also equals the variational distance between the two walkers’
probability distributions.
Note the warm start property is strictly preserved by PF and PG since
µPGt ≤ µPFt ≤ M πPFt = M π.

(57)

In this case the probability that a path x1 , . . . , xt (with x1 ∼ µ and xi ∼ PF (xi−1 ) for i > 1)
passes through ΩB is
t
X
≤
Pr[xi ∈ ΩB ] ≤ M tπ(ΩB ),
(58)
i=1

where we have used first the union bound and then (57). By the above arguments we have
kµ(PFt − PGt )k1 ≤ M tπ(ΩB ).

(59)

Using the spectral gap of PF (from (56)) and the mixing bound in (22) we conclude that
−1 −t/ρ
kπ − µPGt k1 ≤ M tπ(ΩB ) + πmin
e

(60)

We see that the leakage probability increases linearly with t while the usual distance to
stationarity decreases exponentially with t. In many cases this leaves a wide range of t in which
the RHS of (60) can be small.

14

4

Efficient convergence of SQA for the spike cost function

In this section we apply the method from Section 3.2 with the easy-to-analyze chain (π̃, P̃ ) taken
to be the SQA Markov chain with heat-bath worldline updates for the system without the spike,
and (π, P ) equal to the corresponding chain for the spike system.
The subset ΩB will be shown to satisfy π̃(ΩB ) ≤ O(n−c ) for a constant c that we will choose
so that we can prove a walker beginning in ΩG is likely to mix before it hits a point in ΩB .

Congestion of the spikeless chain. Recall that ΩB is defined in terms of a set of
canonical paths {γ} on Ω with congestion ρ̃ for the spikeless chain, together with a subset Ωθ
of points which are excluded from paths in {γ} to obtain a new set of paths with congestion
ρ ≤ O(θρ̃) for the chain with the spike, within the subset ΩG . The spikeless distribution π̃
corresponds to a collection of n non-interacting 1D ferromagnetic Ising models of length L in
the presence of 1-local fields that bias the distribution towards configurations of lower Hamming
weight. The spin-spin coupling is such that each broken bond in x lowers π̃(x) by a factor of
Θ(tanh(ω)).
First we will bound the congestion of heat-bath worldline updates (defined in Section 2.2).
Here it is convenient to represent states x ∈ Ω by their worldlines x = (x̄1 , . . . , x̄n ), where
x̄i := (xi,1 , . . . , xi,L ). For the spikeless system (π̃, P̃ ) spins in different worldlines do not interact
and so the conditional distribution of the i-th worldline is equal to the marginal of the stationary
distribution on that worldline,
X
π̃(x̄0i |x̄01 , . . . , x̄0i−1 , x̄0i+1 , . . . , x̄0n ) =
π̃(x̄1 , . . . , x̄0i , . . . , x̄n )
(61)
x̄j :j6=i

Including a 1/2 probability of the chain not moving at each step in order to make it irreducible,
the probability of a transition P ((z̄1 , . . . , z̄i , . . . , z̄n ), (z̄10 , . . . , z̄i0 , . . . , z̄n0 )) that updates the i-th
worldline (z̄j = z̄j0 for all j 6= i) is
P ((z̄1 , . . . , z̄i , . . . , z̄n ), (z̄10 , . . . , z̄i0 , . . . , z̄n0 )) =

1 X
π̃(z̄100 , . . . , z̄i0 , . . . , z̄n00 )
2n

(62)

z̄j :j6=i

The path γxy from x = (x̄1 , . . . , x̄n ) to y = (ȳ1 , . . . , ȳn ) proceeds by updating the worldlines
in order {1, . . . , n}. The paths have length |γxy | = n. The k-th step of the path γxy will go
through the edge (z, z 0 ) with z = (ȳ1 , . . . , ȳk−1 , x̄k , . . . , x̄n ) and z 0 = (ȳ1 , . . . , ȳk , x̄k+1 , . . . , x̄n ).
To evaluate the sum (30) we apply the standard encoding trick [29]. Define an injective function
η(z,z0 ) which maps the paths through (z, z 0 ) into the state space Ω,
η(z,z0 ) (γxy ) = (x̄1 , . . . , x̄k−1 , ȳk , . . . , ȳn ),
which is injective since z and η(z,z0 ) (γxy ) provide sufficient data to uniquely determine γxy .
Notice that π̃(x)π̃(y) = π̃(z)π̃ η(z,z0 ) (γxy ) , and so the congestion ρ̃(z, z 0 ) is
X
1
π̃(x)π̃(y)|γxy |
π̃(z)P (z, z 0 )
γxy 3(z,z 0 )
X

n
=
π̃ η(z,z0 ) (γxy )
0
P̃ (z, z ) y 3(z,z0 )
xy
X
n
=
π̃(x̄1 , . . . , x̄k−1 , ȳk , . . . , ȳn )
P̃ (z, z 0 ) x̄1 ,...,x̄k−1

ρ̃(z, z 0 ) =

(63)
(64)
(65)

ȳk+1 ,...,ȳn
2

= 2n ,

(66)

15

where in going from (64) to (65) we do not sum over ȳk because it is fixed by z 0 . Finally, since
the edge (z, z 0 ) was arbitrary we have
ρ̃ = O(n2 ).

(67)

Now we will analyze the single-site Metropolis chain for the spikeless system (π̃, P̃M ). The
path from x = ((x1,1 , . . . , x1,L ), . . . , (xn,1 , . . . , xn,L )) to y = ((y1,1 , . . . , y1,L ), . . . , (yn,1 , . . . , yn,L ))
proceeds as follows: for each i in order from {1, . . . , n} update spin (i, j) for j going in order
from {1, . . . , L}. The paths have length |γxy | = nL. Any edge (z, z 0 ) along such a path will
create at most two new broken bonds (i.e. pairs of spins which disagree) along the direction
{1, . . . , L}, and so



π̃(z 0 )
1
min 1,
= Ω (nL)−1 tanh(ω)
(68)
P̃M (z, z 0 ) =
2nL
π̃(z)
Just as for the heat-bath worldline updates, we apply an encoding function to injectively map
the paths through an arbitrary edge into the state space. The only difference from the version
above is that z and η(z,z0 ) (γxy ) may in the worst-case each have two broken imaginary-time
bonds which are not present in either x or y, and so

π̃(x)π̃(y) = O coth(ω)2 π̃(z)π̃ η(z,z0 ) (γxy ) .
(69)
Applying the same calculations in (63)-(65) but with (68) and (69), and using the fact that
−1
coth(ω) = O(ωmin
) = O(nLβ −1 ), the congestion is
ρ̃M = O(L5 n5 β −3 ) = O(n15 β −9/2 ).

(70)

Comparing (70) with (67) shows that there is a large polynomial overhead resulting from singlesite updates, which arises from the strong interactions that occur in the imaginary-time direction
when the transverse field term of the QA Hamiltonian is small. Nevertheless, the bound (70)
is still efficient in the sense of polynomial time, and will suffice to obtain an overall polynomial
run time for single-site Metropolis chain of the spike system.

Ωθ and the spike time distribution. The states in Ωθ which will be excluded from the
set of paths {γ} are those which have |xi | ∈ IS := (n/4 − nζ /2, n/4 + nζ /2) for too many i.
Define 1S : {0, 1}n → {0, 1} to be the indicator function for the spike i.e. 1S (z) = 1 if z ∈ IS ,
and 1S (z) is zero otherwise. The spike time for x ∈ Ω is defined to be
ST(x) =

L
X

1S (xi ).

(71)

i=1

Let  :=

1
2

− α, and define

Ωθ =

x ∈ Ω : ST(x) ≥

L
1

n 2 (1−)−ζ


.

Set β = n/2 so that every x ∈
/ Ωθ satisfies
!



π(x)
Z̃
βnα
= O(1),
≥
exp −
· min ST(x)
π̃(x)
Z
L
x∈Ω
/ θ

(72)

(73)

which shows θ = O(1) in the congestion bound (49). The remainder of the section will be
devoted to computing the m-th moment of the random variable ST ∼ π̃, with m = c/, in order
to show,


L
Pr ST ≥ 1 (1−)−ζ
≤ O(n−c ),
(74)
n2
π̃

16

which is equivalent to the statement π̃(Ωθ ) ≤ O(n−c ).
To calculate the moments hSTm iπ̃ we will relate them to expectation values of the spikeless
quantum system, and use the fact that the latter is exactly solvable because the qubits are
non-interacting. Recall that the quantum expectation value of an operator A can be expressed
as a derivative of the partition function,
hAi =



1 ∂
1 ∂Z(λ)
1  −βH 
=
tr Ae
tr e−βH+λA |λ=0 =
|λ=0 .
Z
Z ∂λ
Z(0) ∂λ

(75)

Passing from Z(λ) to the corresponding Suzuki-Trotter approximation Z(λ) in the above expression changes the value by at most a multiplicative factor of (1 ± O(L−1 ) (a proof of this fact
will appear in a future work [11]). Let {|ki : k = 0, . . . , n} be P
a basis of states for the symmetric
subspace which are labeled by Hamming weight, and let S = k∈IS |kihk|. Since the observable
S is diagonal in the computational basis we can include the term λS into the diagonal part of
the Hamiltonian for the quantum-to-classical mapping and compute,


! n
L
X
1 ∂  X
β f˜(xi ) λhxi |S|xi i Y
hSiσ̃ =
+
φ(x̄j ) |λ=0
(76)
exp
−
Z ∂λ x ,...,x
L
L
j=1
i=1
1
L
#
"
L
n
Y
β PL
1 X
1X
−L
f˜(xi )
i=1
1S (xi ) e
φ(x̄j )
(77)
=
Z x ,...,x L i=1
j=1
1

=

X

L

L−1 ST(x)π̃(x) = L−1 hSTiπ̃

(78)

x∈Ω

Since the low-temperature thermal state and the ground state have a large overlap, the expectation value (78) can be computed within the ground state. To compute the error caused
by this replacement we will examine the density of states of the spikeless system. First add a
constant shift to the Hamiltonian
pso that H̃|ψ̃0 i = 0. Next, let |ψ̃1 i, . . . , |ψ̃n i denote the excited
eigenstates of H̃. Define ∆ := 2 (1 − s)2 + s2 and observe that |ψ̃k i is an eigenstate of H̃ with
eigenvalue k∆, and that the degeneracy of the k-th energy level is nk so
kσ̃ − |ψ̃0 ihψ̃0 |k1 ≤

n
X

e−β∆k

k=1

 

n
= O(ne−n ),
k

(79)

and since  > 0 is a constant this error will be sub-leading.
Since the spikeless system is non-interacting the ground state can be written explicitly, and
the ground state probability distribution on the Hamming weights is a binomial distribution [28].
It follows that hSiσ̃ is asymptotically never larger than the central binomial coefficient times
the width nζ of the spike term, therefore


hSTiπ̃ = O Lnζ−1/2 .
(80)
To obtain (74) we will use the moment inequality,
Pr [ST ≥ b]π̃ ≤

17

hSTm iπ̃
,
bm

(81)

1

with b = Ln− 2 (1−) . By definition we have,
m

hST iπ̃ =

X

L
X

π̃(z1 , . . . , zL )

z1 ,...,zL

=

X
L
X

1S (zt )

(82)

t=1
L
X

π̃(z1 , . . . , zL )

z1 ,...,zL

=

!m

1S (zt1 ) . . . 1S (ztm )

(83)

t1 ,...,tm

h1S (zt1 ) . . . 1S (ztm )iπ̃ .

(84)

t1 ,...,tm

To compute these m-point correlation functions we return to the quantum description (this
generates a multiplicative error of size 1 ± O(L−1 ), which will make a sub-leading contribution
to the m-th moment and thus will be ignored),
h1S (zt1 ) . . . 1S (ztm )iπ̃ = he−τ1 H Se−(τ2 −τ1 )H S . . . e−(τm −τm−1 )H Se−(β−τm )H iσ̃

(85)

where τi := βti /L. Once again we replace the low-temperature thermal state with the ground
state and incur a sub-leading error as in (79). Therefore,
hSTm iπ̃ =

L
X

hψ̃0 |e−τ1 H Se−(τ2 −τ1 )H S . . . e−(τm −τm−1 )H Se−(β−τm )H |ψ̃0 i

(86)

t1 ,...,tm

Since the ground state, the Hamiltonian, and the operator S are all bit-symmetric, the expectation can be evaluated in the symmetric subspace. Expanding each of the terms in (86) in the
basis of symmetric energy eigenstates {|ψ̃k i},
X
hSTm iπ̃ =
e−(τ2 −τ1 )∆k1 · · · e−(τm −τm−1 )∆k1 hψ̃0 |S|ψ̃k1 ihψ̃k2 |S|ψ̃k3 i · · · hψ̃km |S|ψ̃0 i (87)
k1 ,...,km
t1 ,...,tm

States with higher energy will contribute less to the sum over all times t1 , . . . , tm in (86), because
the exponentials in (87) decay more quickly. For ki > 0, the sum over ti can be truncated
whenever τi − τi−1  1/ki ∆.
Since the ground state wave function is a binomial distribution the mean spike time will
only be large when the peak of the ground state is near the support of the spike IS . In the
range of the adiabatic parameter in which this occurs the excited spikeless eigenstates satisfy
hψ̃i |ki ≤ |hψ̃0 |ki| ≤ O(n−1/4 ) for all i = 1, . . . , n and k ∈ IS , because the ground state wave
function is centered on the spike and the excited state wave functions have a greater spread,
which can be seen from the explicit form of the spikeless eigenfunctions given in [24]. Now we
define gi := ti − ti−1 and relabel the sum of t1 , . . . , tm = 0, . . . , L by a sum over the P
gi . For the
purpose of obtaining an upper bound on the m-th moment we relax the constraint i gi = L,
and instead sum over the full range gi = 1, . . . , L for each i. Using these facts we can upper
bound (87) by
X
hSTm iπ̃ ≤ nm(ζ−1/2)
e−g1 ∆k1 · · · e−gm ∆km
(88)
k1 ,...,km
g1 ,...,gm

We will now organize the terms
of (88) according to the number ` of excited energies Ẽki > 0

they contain. There are m
terms
of (88) that contain ` eigenstates above the ground state,
`
and for each ` we must sum over the ga1 , . . . , ga` for which the corresponding ka1 , . . . , ka` are

18

non-zero. Now (88) becomes
hSTm iπ̃ ≤ nm(ζ−1/2)

m  
`
X
X Y
m Lm−`
e−gai ∆kai
` (m − `)! ga ,...ga i=1
`=1

(89)

1
`
ka1 ...ka`

where the factor of Lm−` /(m − `)! results from performing the sum over the m − ` of the gi
PL
which have ki = 0. Now we sum over ga1 , . . . , ga` using the fact that g=1 e−gk ≤ k −1 ,
hSTm iπ̃ ≤ Lm nm(ζ−1/2)

m  
X
m
`=0

Using

m
`



≤ m` and

Pn

`

(β∆)−`

n
X
ka1 ...ka`

`
Y
1
.
k
i=1 i

≤ log(n) + 1, at last this becomes
`
m 
X
log(n) + 1
hSTm iπ̃ ≤ Lm nm(ζ−1/2)
mβ∆

(90)

k=1

(91)

`=0

Since m and ∆ are constant and β = n/2 with fixed  > 0 the terms with inverse powers of β
are sub-leading and so hSTm iπ̃ ≤ O(Lm nm(ζ−1/2) ). Finally, applying (81) with m = c/ yields
the desired result (74).

Adiabatic schedule. Here we show that a discretization of the adiabatic path with 1/ poly(n)
step size is sufficient to fulfill the statement we need for the warm starts in Section 3.3.
We will take the largest value of the adiabatic parameter to be smax = 1 − n−1 so that
||ψ0 (1)i − |ψ0 (smax )i|1 ≤ poly(n−1 ), and the global minimum of the cost function can be obtained by sampling from Π at s = smax with essentially the same probability at it would be
obtained by sampling from the ground state probability distribution.
We will sample from π at several values of the adiabatic parameter s1 , . . . , smax . Define
si := s0 − i∆s, ωi := β(1 − si )/L, ∆ω := β∆s/L, and let π i be the stationary distribution
(17) when the adiabatic parameter is si . At each stage we simulate the Markov chain (19)
for sufficiently many steps to achieve a variational distance to the stationary distribution of
exp(−nΩ(1) ). These errors then add up to a negligible amount. To choose a step size ∆s
satisfying the warm start condition πi+1 ≤ 2πi we’ll use a claim which is inspired by Lemma
5.1 of [8] but is a bit simpler in the classical case.
P
Lemma 4. Let E1 , E2 : A → R be energy functions on a domain A and define Zi := x∈A e−Ei (x)
and pi (x) = e−Ei (x) /Zi for i = 1, 2. If maxx |E1 (x) − E2 (x)| ≤ δ, then | log(Z1 /Z2 )| ≤ δ and
maxx | log(p1 (x)/p2 (x))| ≤ 2δ.
Proof.
Applying the uniform bound |E2 (x) − E1 (x)| ≤ δ for all x ∈ A to the sum Z2 =
P
−E1 (x)
e
leads to e−δ Z1 ≤ Z2 ≤ eδ Z1 , therefore | log(Z1 /Z2 )| ≤ δ. Since | log(p1 (x)/p2 (x))| =
x∈A
| log(p2 (x)/p1 (x))|, we treat the case p1 (x) ≥ p2 (x) (which implies E2 (x) ≥ E1 (x)), while the
other case follows similarly.
max | log(p1 (x)/p2 (x))| = max | log(Z2 /Z1 ) + log(e−(E1 (x)−E2 (x)) )|
x

x

= | log(Z2 /Z1 )| + max | log(eE2 (x)−E1 (x) )|
x

≤ 2δ
Applying Lemma 4 with the form of the stationary distribution (17) we may take


tanh(ωi − ∆ω)
β
δ = ∆sfmax + nL log
= O (β∆s n log(n)) ,
L
tanh(ωi )

therefore taking ∆s = O (βn log(n))−1 fulfills the warm start condition.

19

Quasi-stationary mixing. In this section we will bound the overall run time of SQA
applied to the spike cost function. First we need to show that taking c in (74) to be a sufficiently
large constant will allow the leaky walk for the spike system sample from ΩG according to the
quasi-stationary distribution π for an expected time of nq , for any desired constant q, before it
is eventually likely to escape into ΩB .
−1
Inserting ρ̃ from (67) into (49) yields ρ = O(n2 ). To apply (60) we next must consider πmin
.
−1
From (17) we have log πmin ≤ O(nL log(n)) because there can be L pairs of bits which disagree
in each of the n worldlines. However, according to (17) these disagreements follow a binomial
distribution with mean β(1 − s), and so we may abort the algorithm with exponentially small
probability if it ever encounters a configuration with Ω(β log n) jumps in any worldline, which
−1
allows us to take log πmin
= O(nβ log(n). This implies a mixing time of tmix = O(n3 β log n)
within ΩG for the SQA spike chain with single qubit worldline updates at each value si of the
adiabatic path.
Meanwhile, from (39) together with (74) we have π(ΩB ) ≤ Θ(ρ̃π̃(Ωθ )) = Θ(n2−c ). At each
step si of the adiabatic path, after time t ≥ tmix the leaky random walk mixes to within a
distance O(tπ(ΩB )) so by (60) it suffices to take c = 2 + q + log(1/δ) in order for the leaky walk
to be with distance δ to the stationary distribution π for times tmix ≤ t ≤ Θ(nq ).
Finally, since there are Õ(nβ) steps of the adiabatic path, and each worldline update takes
Õ(n2 β 2 ) time to implement, we obtain a total run time of Õ(n6 β 4 ), which implies the Õ(n7 )
stated in Theorem 1. Repeating the analysis above for single-site Metropolis updates, we com−1
bine the congestion (70) with log πmin
= Õ(nβ) and Õ(nβ) steps of the adiabatic path to arrive
17
bound the total run time by Õ(n ) as stated in Theorem 1.

Acknowledgments
We thank Dave Bacon, Wim van Dam and Alistair Sinclair for helpful conversations. Elizabeth
Crosson gratefully acknowledges funding provided by the Institute for Quantum Information and
Matter, an NSF Physics Frontiers Center (NSF Grant PHY-1125565) with support of the Gordon
and Betty Moore Foundation (GBMF-12500028), and is also grateful for support received while
completing a portion of this work at the MIT Center for Theoretical Physics with funding from
NSF grant number CCF-1111382. Aram Harrow was funded by NSF grants CCF-1111382 and
CCF-1452616 and ARO contract W911NF-12-1-0486.

References
[1] D. Battaglia, G. E. Santoro, and E. Tosatti. Optimization by quantum annealing: Lessons
from hard 3-SAT cases. Phys. Rev. E, 2005, arXiv:cond-mat/0502468.
[2] S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis,
and M. Troyer. Quantum annealing with more than one hundred qubits. Nature Phys.,
page 218, 2014, arXiv:1304.4595.
[3] L. T. Brady and W. van Dam. Quantum Monte Carlo simulations of tunneling in quantum
adiabatic optimization, 2015, arXiv:1509.02562.
[4] L. T. Brady and W. van Dam. Spectral gap analysis for efficient tunneling in quantum
adiabatic optimization, 2016, arXiv:1601.01720.
[5] S. Bravyi. Monte Carlo simulation of stoquastic Hamiltonians, 2014, arXiv:1402.2295.
[6] S. Bravyi, A. J. Bessen, and B. M. Terhal. Merlin-Arthur games and stoquastic complexity,
2006, arXiv:quant-ph/0611021.

20

[7] S. Bravyi, D. P. DiVincenzo, R. I. Oliveira, and B. M. Terhal. The complexity
of stoquastic local Hamiltonian problems. Quant. Inf. Comp., 8(5):0361–0385, 2006,
arXiv:quant-ph/0606140.
[8] S. Bravyi and B. M. Terhal. Complexity of stoquastic frustration-free Hamiltonians. SIAM
J. Comput., 39(4):1462–1485, 2009, arXiv:0806.1746.
[9] P. Caputo, E. Lubetzky, F. Martinelli, A. Sly, and F. L. Toninelli. Dynamics of (2 + 1)dimensional SOS surfaces above a wall: Slow mixing induced by entropic repulsion. Ann.
Probab., 42(4):1516–1589, 07 2014, arXiv:1205.6884.
[10] E. Crosson and M. Deng. Tunneling through high energy barriers in simulated quantum
annealing, 2014, arXiv:1410.8484.
[11] E. Crosson and A. Harrow. Rapidly mixing Monte Carlo for 1D stoquastic Hamiltonians,
2016. In preparation.
[12] V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush, V. Smelyanskiy, J. Martinis, and H. Neven. What is the computational value of finite range tunneling?, 2015,
arXiv:1512.02206.
[13] M. Dyer, L. A. Goldberg, M. Jerrum, and R. Martin. Markov chain comparison. Probability
Surveys, 3:89–111, 2006, arXiv:math/0410331.
[14] M. Dyer, A. Sinclair, E. Vigoda, and D. Weitz. Mixing in time and space for lattice
spin systems: a combinatorial view. In Randomization and approximation techniques in
computer science, pages 149–163. Springer, 2002.
[15] E. Farhi, J. Goldstone, D. Gosset, S. Gutmann, H. B. Meyer, and P. Shor. Quantum
adiabatic algorithms, small gaps, and different paths. arXiv preprint arXiv:0909.4766,
2009.
[16] E. Farhi, J. Goldstone, and S. Gutmann. Quantum adiabatic evolution algorithms versus
simulated annealing, 2002, arXiv:quant-ph/0201031.
[17] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser. Quantum computation by adiabatic
evolution, 2000, arXiv:quant-ph/0001106.
[18] M. B. Hastings. Obstructions to classically simulating the quantum adiabatic algorithm.
Quantum Information & Computation, 13(11-12):1038–1076, 2013, arXiv:1302.5733.
[19] E. Inack and S. Pilati. Simulated quantum annealing of double-well and multiwell potentials. Physical Review E, 92(5):053304, 2015, arXiv:1510.04650.
[20] S. V. Isakov, G. Mazzola, V. N. Smelyanskiy, Z. Jiang, S. Boixo, H. Neven, and M. Troyer.
Understanding quantum tunneling through quantum Monte Carlo simulations, 2015,
arXiv:1510.08057.
[21] M. Jerrum, A. Sinclair, and E. Vigoda. A polynomial-time approximation algorithm for
the permanent of a matrix with nonnegative entries. J. ACM, 51(4):671–697, 2004.
[22] T. Kadowaki and H. Nishimori. Quantum annealing in the transverse Ising model. Phys.
Rev. E, 58:5355–5363, Nov 1998, arXiv:cond-mat/9804280.
[23] S. Kirkpatrick, M. P. Vecchi, et al. Optimization by simmulated annealing. science,
220(4598):671–680, 1983.
[24] L. Kong and E. Crosson. The performance of the quantum adiabatic algorithm on spike
Hamiltonians, 2015, arXiv:1511.06991.
[25] D. Levin, Y. Peres, and E. Wilmer. Markov Chains and Mixing Times. American Mathematical Soc., 2008.
[26] R. Martoňák, G. E. Santoro, and E. Tosatti. Quantum annealing by the path-integral
Monte Carlo method: The two-dimensional random Ising model. Phys. Rev. B, 66:094203,
Sep 2002.

21

[27] S. Muthukrishnan, T. Albash, and D. A. Lidar. Tunneling and speedup in quantum optimization for permutation-symmetric problems, 2015, arXiv:1511.03910.
[28] B. W. Reichardt. The quantum adiabatic optimization algorithm and local minima. In
Proceedings of the thirty-sixth annual ACM symposium on Theory of computing, pages
502–510. ACM, 2004.
[29] A. Sinclair. Improved bounds for mixing rates of Markov chains and multicommodity flow.
Combinatorics, Probability and Computing, 1(4):351–370, 1992.
[30] M. Suzuki. Quantum statistical Monte Carlo methods and applications to spin systems.
Journal of Statistical Physics, 43(5-6):883–909, 1986.
[31] M. Suzuki, S. Miyashita, and A. Kuroda. Monte Carlo simulation of quantum spin systems.
Prog. Theor. Phys., 58(5):1377–1387, 1977.

22

