論文 / 著書情報
Article / Book Information
Title

Dynamics of reverse annealing for the fully connected p-spin model

AuthorsYu Yamashiro, Masaki Ohkuwa, Hidetoshi Nishimori, Daniel A. Lidar

Citation

Physical Review A, Vol. 100,

Pub. date

2019, 11

DOI

http://dx.doi.org/10.1103/PhysRevA.100.052321

Copyright(C) 2019 American Physical Society

Powered by T2R2 (Tokyo Institute Research Repository)

PHYSICAL REVIEW A 100, 052321 (2019)
Editors’ Suggestion

Dynamics of reverse annealing for the fully connected p-spin model
Yu Yamashiro
Department of Physics, Tokyo Institute of Technology, Nagatsuta-cho, Midori-ku, Yokohama 226-8503, Japanand Jij, Inc., High Tech Hongo Building 1F, 5-25-18 Hongo, Bunkyo, Tokyo 113-0033, Japan

Masaki Ohkuwa
NTT DATA Mathematical Systems, Inc., Shinanomachi, Shinjuku-ku, Tokyo 160-0016, Japan

Hidetoshi Nishimori
Institute of Innovative Research, Tokyo Institute of Technology, Nagatsuta-cho, Midori-ku, Yokohama 226-8503, Japan
and Graduate School of Information Sciences, Tohoku University, Sendai 980-8579, Japan

Daniel A. Lidar
Departments of Electrical Engineering, Chemistry, and Physics & Astronomy, University of Southern California,
Los Angeles, California 90089, USA
and Center for Quantum Information Science & Technology, University of Southern California, Los Angeles, California 90089, USA(Received 12 July 2019; published 18 November 2019)
Reverse annealing is a relatively new variant of quantum annealing in which one starts from a classical stateand increases and then decreases the amplitude of the transverse field, in the hope of finding a better classicalstate than the initial state for a given optimization problem. We numerically study the unitary quantum dynamics
of reverse annealing for the mean-field-type p-spin model and show that the results are consistent with thepredictions of equilibrium statistical mechanics. In particular, we corroborate the equilibrium analysis predictionthat reverse annealing provides an exponential speedup over conventional quantum annealing in terms of solving
the p-spin model. This lends support to the expectation that equilibrium analyses are effective at revealing
essential aspects of the dynamics of quantum annealing. We also compare the results of quantum dynamics
with the corresponding classical dynamics to reveal their similarities and differences. We distinguish between
two reverse annealing protocols we call adiabatic and iterated reverse annealing. We further show that iterated
reverse annealing, as has been realized in the D-Wave device, is ineffective in the case of the p-spin model butnote that a recently introduced protocol (“h-gain”), which implements adiabatic reverse annealing, may lead to
improved performance.
DOI: 10.1103/PhysRevA.100.052321

I. INTRODUCTION

Quantum annealing (QA) is a quantum-mechanical metaheuristic for combinatorial optimization problems [1–6], andits strict adiabatic realization is known as adiabatic quantum
computation [7–9]. In the conventional formulation of QA,
the initial condition is chosen to be the ground state of
the transverse field, which is a uniform superposition of all
possible classical states. Perdomo-Ortiz et al. proposed sometime ago [10] to instead start from an appropriately chosenclassical state and gradually increase and then decrease theamplitude of the transverse field to find a better solution
than the initial classical state. They called this protocol“sombrero adiabatic quantum computation,” based on the
shape of the amplitude of the transverse field, and showedby means of a few numerical examples that the method is
indeed effective if the initial condition is properly chosen.
Further refinements were proposed in Ref. [11], and theprotocol whereby the initial state is classical rather than aquantum superposition is now called reverse annealing. The
current generation of the D-Wave quantum annealer, D-Wave
2469-9926/2019/100(5)/052321(11)2000Q, implements reverse annealing, and several recent
studies have been reported on using this feature in different
contexts, e.g., quantum simulation [12], matrix factorization[13], portfolio optimization [14], and mid-anneal pausing[15].
A theoretical study of reverse annealing was initiatedrecently by some of us in Ref. [16], in which equilibrium
statistical mechanics was employed to study when and how
reverse annealing is effective. It was shown explicitly thatreverse annealing can turn first-order quantum phase transitions into second-order transitions in a simple problem, thep-spin model, by choosing an appropriate annealing path. It is
generally the case that a first-order quantum phase transitionis associated with an exponentially closing energy gap 
between the ground state and the first excited state as afunction of system size, whereas a second-order transition isassociated with a polynomially closing gap (although exceptions are known [17–19]). On the other hand, according to
the adiabatic theorem of quantum mechanics, the computation
time τ is inversely proportional to a small power of the052321-1©2019 American Physical Society

YAMASHIRO, OHKUWA, NISHIMORI, AND LIDARPHYSICAL REVIEW A 100, 052321 (2019)energy gap  [20–22]. Taken together, these facts suggest thatreverse annealing provides an exponential speedup relative tothe conventional QA protocol for the p-spin model. In thiswork we study this intriguing connection from the dynamical perspective by numerically solving the time-dependent
Schrödinger equation of the reverse annealing protocol forthe p-spin model. Our goal is to test whether the exponential
speedup suggested by the static, equilibrium statistical mechanics analysis does indeed hold.
This is nontrivial, since, while a quantum annealer at very
long annealing times is likely to experience a quasistatic evolution, returning a final population that is close to a Boltzmann
distribution of the Hamiltonian at a single (freeze-out) point
during the annealing [23], results from equilibrium statisticalmechanics clearly do not reveal all aspects of the dynamical behavior of quantum systems. In particular, quantum
annealing on relatively short timescales compared to the quantum gap can give rise to distinctly nonequilibrium features
such as coherent oscillations due to quantum interference[24,25].
In the present work we evaluate the unitary, closed-system
dynamics of reverse annealing by direct numerical integrationof the time-dependent Schrödinger equation for large systemsizes. This is possible due to a special symmetry of thep-spin model. To this end, it is useful to distinguish two
types of reverse annealing. The first is the one adopted inRef. [16] to facilitate analytical treatment by equilibrium
statistical mechanics. In this formulation, an additional term isintroduced into the Hamiltonian, which is usually composedof just two terms representing the classical cost function
and the transverse field. The additional, third term serves toenforce a given classical ground state as an initial condition.
The system is then expected to follow an adiabatic path to thefinal state. We call this protocol adiabatic reverse annealing(ARA). In the second formulation the initial state is a classicalstate, an excited state of the initial (classical) Hamiltonian (thecost function), as formulated in Refs. [10,11]. The annealing
process therefore follows a complicated combination of instantaneous eigenstates of the conventional two-term Hamiltonian, without any additional term. It is clear that dynamical
studies are indispensable for this case. This latter protocol isreverse annealing, as is widely used in practice on the D-Wave
device [12–15] and is often used iteratively in practice, thatis, the final state of a single cycle is fed into the system as
its initial state of the next cycle under the expectation thatthe result improves iteratively. We therefore call this protocol
iterated reverse annealing.
This paper is organized as follows. In Sec. II, we analyze
adiabatic reverse annealing. After reviewing the formulation
and results on static properties of the p-spin model, we
solve the time-dependent Schrödinger equation numericallyfor large systems to test the consistency of dynamics and static
properties. We also study two classical versions of the model
to clarify the similarities and differences between quantum
and classical dynamics. Section III deals with iterated reverse annealing by numerical solution of the time-dependent
Schrödinger equation, as well as by spectral analysis of theinstantaneous Hamiltonian, with an emphasis on the special
structure of the p-spin model. The last section is devoted toconclusions.

II. ADIABATIC REVERSE ANNEALING

In adiabatic reverse annealing (ARA), we add a term Ĥinit
to the Hamiltonian to fix the initial state to the ground state of
the initial Hamiltonian [16],Ĥ (t ) = s(t )Ĥ0 + [1 − s(t )][1 − λ(t )]Ĥinit+ [1 − s(t )]λ(t )V̂TF ,(1)

where Ĥ0 is the cost function to be minimized and is a function
only of the Pauli matrices {σ̂iz }. The last term represents thetransverse field,
NV̂TF = −σ̂ix ,(2)i=1

where N is the total number of sites (qubits). The parameter controls the strength of the transverse field relative to theother terms. The initialization term is written asNĤinit = −i σ̂iz ,(3)
i=1

where {i (= ±1)} denotes the given classical initial state. Thisinitial Hamiltonian fixes σ̂iz to i as its ground state. The timedependent parameters s(t ) and λ(t ) both change from 0 to 1
as time t proceeds from 0 to τ , where τ is the computation
time. Thus the initial Hamiltonian is Ĥinit and the final one isĤ0 . The quantum term V̂TF is effective only in intermediatetimes. The existence of the term Ĥinit allows us to analyze thesystem properties by equilibrium statistical mechanics under
the assumption of adiabatic evolution.1
We study the p-spin model in this paper,

p
N
1  z
Ĥ0 = −Nσ̂,(4)
N i=1 iwhere p  3 is an integer. The ground state of this model istrivial, σiz = 1 ∀i (and σiz = −1 ∀i for p even). This simple
problem is known to be hard for conventional quantum annealing due to a first-order phase transition [27]. Since Ĥ0 andV̂TF are symmetric with respect to an arbitrary permutation of
the site index i, we do not lose generality by the following
assignment of i :

+1 for i  Nc
(5)
i =−1 for i > Nc,
where the probability c (1/2  c  1) for i to be +1 ischosen such that Nc is an integer. This means that the initial
magnetization is 2c − 1. Thus, the larger c is, the closer
the initial state is to the correct answer to the optimization
problem, σiz = 1, ∀i.2
Permutation symmetry allows us to write the Hamiltonian
in terms of only a set of total spin operators of two subsystems,
1  x,zσ̂ ,2 i=1 iNc

Ŝ1x,z ≡Ŝ2x,z ≡N
1  x,zσ̂ ,2 i=Nc+1 i(6)

1A recent update of the D-Wave quantum annealer introduced afunction ‘‘h-gain,” which can be used to realize ARA [26].
2
For p even, another degenerate ground state σiz = −1 (∀i) exists.
But we do not lose generality by restricting ourselves to the subspacewith non-negative magnetization.

052321-2

DYNAMICS OF REVERSE ANNEALING FOR THE FULLY …

as

PHYSICAL REVIEW A 100, 052321 (2019)



 p
2 z
Ŝ1 + Ŝ2z
Ĥ0 = −N
,
N




Ĥinit = −2 Ŝ1z − Ŝ2z , V̂TF = − Ŝ1x + Ŝ2x .
2

(7)
(8)

2

Since (Ŝ1 ) and (Ŝ2 ) commute with the total Hamiltonian,
the system stays in the subspace with the largest values of(Ŝ1 )2 and (Ŝ2 )2 because the initial state lies in this subspace.
This greatly reduces the size of the Hilbert space to be studied,from exponential to polynomial in N, O(N 2 ).
A. Static properties

We first review the results presented in Ref. [16] withsome generalizations. A statistical-mechanical analysis of thep-spin model for the case of ARA of Eq. (1) with s and λ taken
as fixed parameters, leads to the following expression for thefree energy at zero temperature (the ground state):
f = s(p − 1)m p−1− c [spm p−1 + (1 − s)(1 − λ)]2 +  2 (1 − s)2 λ2− (1 − c) [spm p−1 − (1−s)(1−λ)]2 +  2 (1−s)2 λ2 ,(9)where m satisfies the following self-consistency condition tominimize the free energy:m=c

spm p−1 + (1 − s)(1 − λ)[spm p−1 + (1 − s)(1 − λ)]2 +  2 (1 − s)2 λ2+ (1 − c)spm p−1 − (1−s)(1−λ).
[spm p−1 − (1−s)(1−λ)]2 +  2 (1−s)2 λ2(10)

Derivations are detailed in Appendix A of Ref. [16], where 
was fixed to 1, but a generalization to arbitrary  as given
here is straightforward. We focus our analysis on the caseof p = 3, but other values of p(4) show qualitatively thesame behavior. Figure 1 is the λ-s phase diagram for  = 1.
The axis λ = 1 corresponds to conventional QA, where s ischanged from 0 to 1. Along this axis, a first-order phasetransition at s  0.4 marks where QA becomes inefficient,since the energy gap closes exponentially there. As λ isreduced from 1 toward the inner part of the phase diagram,the effects of the initial Hamiltonian Ĥinit come into play,
and the line of first-order phase transitions is split into twoparts, one for larger λ and the other for smaller λ, if c isabove a threshold, i.e., if the initial condition is reasonably
close to the final answer. This means that a path exists thatconnects the initial point λ = s = 0 and the goal λ = s = 1
without crossing a first-order transition. In particular, when wechoose the path s = λ, the diagonal of the phase diagram, thesystem encounters a first-order phase transition when c = 0.7
and 0.8 but not for c = 0.9. Since a larger c means that theinitial condition is closer to the correct answer, the phase
diagram shows that a good choice of the initial classical state
exponentially accelerates the computation by removing an
exponentially closing energy gap.
The relation between the order of a phase transition predicted by statistical mechanics and the scaling of the energy
gap for large but finite N is not proven in general and should

FIG. 1. Static phase diagram for  = 1 and p = 3. The axis λ =1 corresponds to conventional QA. ARA starts from λ = s = 0 andends at λ = s = 1. Lines for c = 0.7, 0.8, and 0.9 represent firstorder phase transitions.

be checked by other methods. Indeed, while usually a first(second) order transition is accompanied by an exponentially(polynomially) closing gap, as mentioned in the Introduction,
there exist examples in which the gap closes nonexponentially at a first-order phase transition [17–19]. We therefore
evaluated the energy gap by direct numerical diagonalization,
which is possible for relatively large system sizes due to thehigh symmetry of the problem. As seen in Fig. 2, which isa log-log plot of the gap as a function of system size fors = λ, only the data for c = 0.9 behave nonexponentially.
This conforms with the standard expectation from the phase
diagram of Fig. 1, where the path s = λ avoids the line of

FIG. 2. System size dependence of the energy gap for conventional QA (λ = 1), shown blue dashed line, and ARA (λ = s) for
several values of c. The linear behavior in the c = 0.9 case isconsistent with inverse polynomial scaling, while the curvature seen
for c = 0.7 and 0.8 is consistent with inverse exponential scaling.
Here, as in Fig. 1, we set  = 1 and p = 3.

052321-3

YAMASHIRO, OHKUWA, NISHIMORI, AND LIDARPHYSICAL REVIEW A 100, 052321 (2019)

(a)

(b)

(c)

FIG. 3. Phase diagrams in the s-λ plane for p = 3 for different values of the amplitude of transverse field for (a)  = 2, (b)  = 5, and(c)  = 10. Curves indicate first-order phase transitions.

first-order transitions only when c = 0.9. It is interesting thatthe data for c = 0.8 behaves almost polynomially up to about
N = 50. This should reflect the fact that the line s = λ crosses
the first-order transition line close to its termination point in
Fig. 1: close to the termination point, which is a critical point
or a second-order transition point, the width and height of theenergy barrier between two coexisting states should be small,and the system’s behavior is close to that at a second-ordertransition, as long as the size is not too large.
Larger values of  change the phase diagram quantitatively, if not qualitatively. As seen in Figs. 3(a) ( = 2), 3(b)( = 5), and 3(c) ( = 10), the line of first-order transitions
for c = 0.8 and 0.9 on the right part of the phase diagram
extends toward the left part as  increases. This is to be
expected, because for large  the intermediate Hamiltonian ismostly dominated by the transverse-field term V̂TF and is close
to the initial Hamiltonian of conventional QA. It then follows
that the first-order transition of conventional QA persists evenfor relatively large values of c.

transverse field  = 1, ARA does not necessarily performbetter than conventional QA. For a larger  = 2, ARA has ashorter TTS than QA for the parameter values we have tested.
(a)B. Schrödinger dynamics

We next report our results on the closed system dynamics.
Let |φ denote the ground state of the cost function Ĥ0
and |ψ (τ ) the actual state reached after time τ , subject to
Schrödinger dynamics. Thenpe (τ ) = 1 − |φ|ψ (τ )|2

(b)

(11)

is the error probability, i.e., the probability of not reaching the
ground state. This quantity is plotted in Fig. 4 as a functionof system size. For  = 2, the error probability increases
polynomially for large c, which is consistent with the phase
diagram in Fig. 3.
A standard measure of computation time in quantum annealing is the time to solution (TTS), as defined by [28]TTS(τ, pd ) = τlog(1 − pd ).
log pe (τ )(12)

The TTS is the effective time it takes using “runs” lasting timeτ to find the correct solution at least once with probability pd ,which we set to 0.99 or higher.
Results for the TTS are depicted in Fig. 5 as a function ofτ for  = 1 and  = 2, and fixed system size. For a weaker

FIG. 4. Error probability as a function of system size for p =3 and τ = 100. In panel (a) we set  = 1 and in (b)  = 2. Bluedashed lines are for conventional QA. In panel (a), the final error
probabilities of ARA are close to those of conventional QA. In panel(b), in contrast, ARA with c = 0.8 and 0.9 has much smaller errors
than QA.

052321-4

DYNAMICS OF REVERSE ANNEALING FOR THE FULLY …

PHYSICAL REVIEW A 100, 052321 (2019)(a)

FIG. 6. Size dependence of TTS of conventional QA (bluedashed line), ARA with  = 1, c = 0.8 (green dash-dotted line), and = 2, c = 0.8 (green solid line). As seen in the phase diagramsin Figs. 1 and 3(a), ARA with  = 1, c = 0.8 encounters a firstorder transition whereas  = 2, c = 0.8 does not. This difference instatics is reflected in the dynamics as the exponential and polynomial
dependence of the optimal TTS.

(b)

whereas the static phase diagram represents the long-time andthermodynamic limits, τ → ∞, N → ∞.
C. Classical spin-vector dynamics

FIG. 5. TTS for N = 45 and p = 3 as a function of τ for (a)  =1 and (b)  = 2. In the case of larger , the TTS for ARA is shorter
than that for QA. Note the existence of minima in the TTS for  = 2.

Interestingly, minima exist in Fig. 5 for  = 2. This in factholds for a broad range of system sizes, which allows us to extract the optimal annealing time τ , τopt = argminτ TTS(τ, pd )as a function of N and hence the optimal TTS scaling withsize [28,29]. The optimal TTS, estimated from the data for
c = 0.8 as in Fig. 5, is plotted in Fig. 6 as a function of systemsize. It is clearly seen that ARA with  = 2, c = 0.8 exhibits
polynomial scaling, whereas ARA with  = 1, c = 0.8 andconventional QA are exponential.3 This is consistent with thestatic phase diagram of Figs. 1 and 3(a), where the path s = λdoes not cross a first-order phase transition if  = 2, c = 0.8
in Fig. 3(a), but conventional QA along λ = 1 and ARA with
 = 1, c = 0.8 in Fig. 1 do. This is quite nontrivial because
the TTS is a purely dynamical measure for finite-size systems,

3

Since this scaling is derived without the existence of an optimal
annealing time, the true scaling can only be worse, as shown in
Ref. [30]. Essentially, the reason is that the optimal TTS at small
sizes N is smaller than the TTS shown.

It is instructive to compare the results of the previous sections with those of a counterpart classical model. In the spinvector dynamics (SVD) model [31–34], the system evolves
according to the classical Hamilton dynamics under the semiclassical potential VSC corresponding to the Hamiltonian ofEq. (1),
VSC ≡  (s)|Ĥ (s)| (s)= −sN (n1 sin θ1 cos φ1 + n2 sin θ2 cos φ2 ) p− (1 − s)(1 − λ)N (n1 sin θ1 cos φ1 − n2 sin θ2 cos φ2 )− (1 − s)λN (n1 cos θ1 + n2 cos θ2 ),

(13)

where
Nc
N − Nc, n2 =,
NN
and | (s) is the spin-coherent state,
n1 =

Nc

| (s) =cos
i=1

θ1 (s)
θ1 (s) iφ1|−i + sin
e |+i
2
2

N

×

cos
i=Nc+1(14)

θ2 (s) iφ2
θ2 (s)|−i + sin
e |+i .
2
2(15)

The initial conditions are φ1 = φ2 = 0 and θ1 = −θ2 = π /2.
We used the equation of motion, which is derived from thepath-integral formulation and permutation symmetry [34].
Thetimeevolutionofthemagnetization
 (s)| N2 (Ŝ1z + Ŝ2z )| (s) under the SVD model is compared052321-5

YAMASHIRO, OHKUWA, NISHIMORI, AND LIDARPHYSICAL REVIEW A 100, 052321 (2019)with its ARA quantum counterpart in Fig. 7. No essential
differences are observed between classical and quantum
dynamics in Figs. 7(a) and 7(b). However, in Fig. 7(c) thequantum dynamics succeeds in coming close to the rightanswer m = 1 at the end, but the classical dynamics fails.
To explain this behavior, we note that for the parametervalues of Fig. 7(c) [which lies between what is shown in
Figs. 3(a) and 3(b)], the evolution path crosses the first-orderphase transition line near its termination in the middle of the
λ-s phase diagram, where the energy barrier across the firstorder transition is thin and low. Quantum dynamics apparently
tunnels through the barrier, but the classical SVD algorithm
gets stuck because there is no classical mechanism to go
through or over the barrier, however thin or low the barrier is.
Figure 8 supports this viewpoint by illustrating the appearanceof an energy barrier in the semiclassical potential. Figure 9
conspicuously demonstrates the marked difference between
classical and quantum dynamics, which is not clear in thestatic phase diagram. It shows the final value of the magnetization as a function of . Quantum dynamics exhibits a gradual
deterioration as  increases beyond 3.4, corresponding to alower tunneling rate, whereas the classical SVD suddenly fails
beyond a threshold close to  = 3.4, where the energy barrier
seen in Fig. 8 appears.

(a)

(b)

D. Classical spin-vector Monte Carlo

Since the energy barrier is thin and low for Fig. 7(c), theclassical system can hop over the barrier if we introduce
thermal fluctuations. Figure 10 shows the results of spinvector Monte Carlo [33], a finite-temperature version of SVD,
in which the angle variables in the semiclassical potential are
updated stochastically according to the usual Monte Carlorule. As seen in this figure, the system succeeds in reachinga good answer, m close to 1, if the temperature is chosen
appropriately, but if the temperature is too high, the final value
is far from 1. Thus, SVMC requires some fine-tuning of thetemperature, which can be viewed as a limitation.

(c)III. ITERATED REVERSE ANNEALING

In iterated reverse annealing (IRA), one measures the statein the computational (classical) basis after a single cycle of
RA and starts the next cycle from the classical state thus
obtained. IRA, which is a typical reverse annealing protocol as
implemented on D-Wave devices, is fundamentally different
from ARA because it expects to obtain a ground state by
nonadiabatic transitions. Therefore, to investigate whether
the IRA protocol works like ARA under unitary evolution.
Let Û denote the time evolution operator generated by the
Hamiltonian of reverse annealing without the initialization
term present in ARA,
(16)

FIG. 7. Time evolution of the magnetization under ARA (orange
full line) and SVD (blue dashed line) as a function of the normalized
annealing time s for (a)  = 1, (b)  = 2, and (c)  = 4 with N =50, c = 0.8, p = 3, and τ = 40.

where A(t ) = s(t ), B(t ) = 1 − s(t ) with an appropriate choiceof the function s(t ). The key difference from standard QA is
that now s(t ) is nonmonotonic and has a minimum smin , asshown in Fig. 11.
The value of this minimum relative to the critical value sc
plays an important role, since it determines whether or not

the system crosses a phase transition: it does when smin < sc .
Time changes from t = 0 to t = τ , where τ is the annealing
time of a single cycle, and s(0) = s(τ ) = 1. The total annealing time with r cycles is rτ . The initial state is classical, i.e.,a computational basis state |i (where i ∈ [0, . . . , 2N − 1] isĤ (t ) = A(t )Ĥ0 + B(t )V̂TF ,052321-6

DYNAMICS OF REVERSE ANNEALING FOR THE FULLY …
(a)

(b)

(c)

(d)

PHYSICAL REVIEW A 100, 052321 (2019)(a)

FIG. 8. Contour plot of the semiclassical potential VSC as a function of the angles sin θ1 and sin θ2 . Panels (a) s = λ = 0.2 and (b) s =λ = 0.3 for  = 2 show that the minimum of the potential evolves
continuously as a function of time s(= λ) without encountering aphase transition, in agreement with the phase diagram Fig. 3(a). This
leads to a large value of the magnetization at the end of the anneal
as in Fig. 9. By contrast, if  is larger than a threshold value close to
3 as in panels (c) and (d) where  = 3.4, the minimum at an earlier
time s = λ = 0.3 in (c) splits into two (almost) degenerate minima
as in (d) s = λ = 0.49, and the system has to transition from a local
minimum near the center of the figure to the global minimum at theright-top. This transition is possible by quantum tunneling, since thebarrier width and height are not very large at these parameter values,
but is impossible for the classical SVD algorithm, resulting in the
difference seen in Fig. 9.

the state index), and the probabilityof the final state being | jafter a cycle is p ji (τ ) = k∈S(e j ) |k|Û (τ )|i|2 , where S(e j )is the index set of the degenerate states with energy e j . This
process is realized in the current version of the D-Wave 2000Q
device, except that Û is replaced by open-quantum-systemevolution, which means that thermal effects play a role as wellduring each cycle in the real device.

FIG. 9. Final magnetization as a function of the strength oftransverse field  for N = 50, τ = 40, c = 0.8, and p = 3. Theorange solid line represents the case of ARA and blue-dashed is for
SVD.

(b)

(c)

FIG. 10. Time evolution of the magnetization under spin-vector
Monte Carlo at finite temperature β −1 with N = 50, c = 0.8,  = 4,
and 500 sweeps (a sweep is a complete Metropolis update of all thespins at a given s) where β = (a) 10, (b) 5, and (c) 1. The width ofthe light-blue hue is the standard deviation calculated from 100 runs.
(a) Since β = 10 is too large (temperature too low) for the system
to escape the basin of the initial state, the magnetization just wiggles
around the initial value. (b) At a higher temperature β = 5, the state
can hop over the energy barrier to reach the correct ferromagnetic
state in the end. (c) If the temperature is too high β = 1, the state
becomes completely random.

052321-7

YAMASHIRO, OHKUWA, NISHIMORI, AND LIDARPHYSICAL REVIEW A 100, 052321 (2019)

(a)

(b)
FIG. 11. Annealing schedule of a single cycle of IRA for τ = 10.
The schedule in IRA is nonmonotonic as a function of time, startingand ending at s = 1.

IRA is an algorithm that repeats the cycle of RA andmeasurement. Therefore, IRA can be considered as a series ofclassical Markov transitions, and a single cycle is a stochastictransition to the final classical state from the initial classicalstate. Using the overlap p ji (τ ) above, the classical transition
probability matrix P has elements Pji = p ji (τ ). Let π denotea classical probability vector whose ith component πi denotes
the probability of obtaining the ith computational basis state.
After r cycles, the probability becomes
π j(r) =(Pr ) ji πi(0) ,(17)

FIG. 12. Distribution of magnetization of the final state withvarious initial conditions (indicated by the values of c) with two
different minimum values of s [above (a) smin = 0.5 and below(b) smin = 0.3 the critical s value], and N = 50 and τ = 10.

iwhere πi(0) is the given initial state. The probability of obtaining the ground state is π0(r) .

A sufficient condition for IRA to succeed is that the
ground-state probability increases after a single cycle. Correspondingly, excitations should be suppressed in a single cycle.
Unfortunately, this condition does not seem to be satisfied
for the p-spin model as shown below, at least for the specificschedule choice we tested.
We choose a quadratic function of s(t ) with a minimum
smin as in Fig. 11. Figures 12(a) and 12(b) show the distribution of the magnetization of the final states after a single
cycle, starting from various values of c with N = 50, τ =10, and p = 3. In Fig. 12(a), since smin = 0.5 is above thecritical value sc ,4 and hence the system does not undergo aphase transition, the magnetization stays close to the initialvalue (given by c), though it is shifted somewhat towards
m = 0. If smin is below the transition point as in Fig. 12(b),the system does cross the phase transition and, as expected,
the final magnetization deviates more from the initial value.
Specifically, the final magnetization moves toward m = 0,
especially for c ∈ [0.3, 0.7]. However, if the initial state is theground state, the system seems to stay near the initial state.
To better understand what happens in Fig. 12, we havecalculated the energy spectrum and occupation probability.
Results are shown in Fig. 13 for smin = 0.3.
4

See Fig. 3 (left panel) of Ref. [35].

Because of the symmetry of the annealing schedule of IRA
with respect to s = 21 , avoided crossings in the first half appear
also in the second half. Therefore, when annealing starts fromthe ground state as shown in Fig. 13(a), the state becomes
excited once and then returns to the ground state at the second
avoided level crossing, resulting in coherent oscillations of theground-state probability as a function of the total time τ [24].
It can also happen that due to the large number of avoided
crossings in the first and second halves, an excited state canbe further excited to higher energy levels, as in Fig. 13(b).
The situation is far from simple, and the general tendency isthat the distribution among states becomes broader as cycles
are iterated except for the limiting case of adiabatic annealing,where the total evolution time τ is large on the scale set by
the inverse of the minimum gap between the ground and firstexcited state.
The density plots in Figs. 14 and 15 show the transition
probabilities (Pr ) ji from the initial state i (horizontal axis) to
the final state j (vertical axis) in single cycles (r = 1) andafter multiple cycles (r = 3 and 5) with smin = 0.5 (Fig. 14)
and smin = 0.3 (Fig. 15). For example, in Fig. 14(a), if westart from the ground state i = 0, the system reaches the sameground state ( j = 0) with a large probability, as indicatedby the yellow box, since the computation time τ = 30 isrelatively long and close to adiabatic. Other states along thecolumn i = 0 have lower probabilities (dark colors). Similarly, if the initial condition is an excited state i > 0, the samestate j = i is reached with a high probability, resulting in a

052321-8

DYNAMICS OF REVERSE ANNEALING FOR THE FULLY …

(a)

PHYSICAL REVIEW A 100, 052321 (2019)
(a)

(b)

(c)

(d)

(b)

FIG. 13. Energy spectrum of IRA during a single cycle andoccupation probability (represented by the red thickness) for N = 10,
p = 3, τ = 10, and smin = 0.3. The black dashed lines represent theenergy spectrum. In (a) the system starts in the ground state (c = 1)
and in (b) it starts in an excited state (c = 0.8).

series of bright colors along the diagonal. Repetition of IRA[panels (b) and (d)] or shorter computation times [panels (c)and (d)] results in scattering of the probabilities away fromthe diagonal. The situation is similar for smin = 0.3 < sc in
Fig. 15 but with a slightly different structure. In any case, theprobability to reach the ground state at the end of the process,shown along the top row in each panel, remains low unless
we start from the ground state itself, implying that IRA undercoherent dynamics is unsuccessful in the present problem.
To what extent this is a special feature of the p-spin
model, and whether IRA may work better in other cases, is an
important topic for future work. Moreover, it is to be expectedthat thermal relaxation will reduce excited-state occupation
probabilities, so that an open-system study of IRA is likely tofind better performance for the finite- (but low) temperaturep-spin model than reported here.
IV. CONCLUSIONS

Reverse quantum annealing protocols have generated muchinterest recently, as they may potentially overcome some ofthe obstacles that cause conventional forward quantum an-

FIG. 14. Transition probabilities (Pr ) ji from state i (horizontalaxis) to state j (vertical axis) with i = 0 the ground state, i = 1 thefirst excited state, etc., for p = 3, N = 10, smin = 0.5, where (a) τ =30, r = 1; (b) τ = 30, r = 5; (c) τ = 8, r = 1; and (d) τ = 8, r = 3.
Panels (b) and (d) are for the results after r repetitions. In this caseof smin = 0.5, the phase transition is not crossed so the probabilityof staying in the ground state is high, as seen in the left-top elementin each panel. The general tendency is clearly that transitions to theground state, shown in the top row, are small except for the top-leftcorner, which represents the probability to stay in the initially given
ground state. Panel (a) has a rather large τ and the system stays close
to adiabatic, which is reflected in the high probabilities along the
diagonal. In panel (c), with a smaller τ , some diabatic transitions
take place as seen by the increased brightness of the off-diagonal
elements. Since smin = 0.5 > sc , the energy gap 01 between the
ground state and the first excited state remains relatively large, buttransitions do occur between excited states because the energy gaps
between higher energy states are smaller than 01 .

nealing to fail. In this work we presented a numerical study ofthe closed-system quantum dynamics of two types of reverse
annealing, adiabatic reverse annealing (ARA) and iterated
reverse annealing (IRA), for the p-spin model. For ARA, wedemonstrated that the dynamical behavior of the system isconsistent with its static phase diagram. The latter suggested
that ARA may be able to exhibit an exponential speedup
over QA for properly chosen initial conditions by finding apath that avoids a first-order quantum phase transition thatis encountered by QA. Our study confirms this expectation
using the time-to-solution (TTS) measure: by establishing an
optimal annealing time for each problem size we were able
to extract the optimal TTS scaling and show that it scales
polynomially for ARA. At the same time we found thatthe lower bound for the TTS of QA is exponential. This isremarkable not only because ARA is thus a demonstrably
better protocol than QA for at least one class of (admittedly
trivial) optimization problems, but also since it lends credence
to the predictions of the static analysis, which is often much
simpler to perform than solving the dynamics.

052321-9

YAMASHIRO, OHKUWA, NISHIMORI, AND LIDAR
(a)

(b)

(c)

(d)

PHYSICAL REVIEW A 100, 052321 (2019)

FIG. 15. Transition probabilities (Pr ) ji from state i to state jfor p = 3, N = 10, smin = 0.3, after r repetitions, where (a) τ =30, r = 1; (b) τ = 30, r = 5; (c) τ = 10, r = 1; and (d) τ = 10, r =5. Since smin = 0.3 < sc , a phase transition occurs, accompanied
by a relatively small energy gap between the ground state and thefirst excited state, which causes a diabatic transition to the firstexcited state, after which the state returns to the ground state at thesecond avoided crossing, which is represented as bright yellow inthe top-left corner of panel (a). When started from an excited statei > 0, the system does not reach the ground state. The reason isthat the existence of many avoided crossings causes both upwardand downward transitions and does not unilaterally favor downward
transitions toward the ground state.

to work and provide an enhanced probability of finding the
ground state after multiple cycles, the probability distributionof the final state after a single cycle should shift toward
lower energy states than the initial state. Unfortunately, thiscondition is not satisfied in the p-spin model, as is clear
from Figs. 14 and 15. The generality of this conclusion iscurrently unclear, and it may be associated with the particular
structure of the energy spectrum of the p-spin model, depicted
in Fig. 13. In particular, the gap structure is such that thep-spin model does not appear to lend itself to a diabatic cascade, wherein upward diabatic transitions are accompanied by
downward transitions for properly tuned annealing schedules[34,36,37]. However, our p-spin model result suggests theimportant lesson that there exist examples where IRA does
not yield a ground-state probability enhancement.
The present study concentrates on purely unitary quantum
dynamics at zero temperature. In practice, real devices are
exposed to an environment and are open quantum systems
with state transitions affected by thermal effects. In the DWave device, for example, thermal relaxation seems to lead
to better performance, as exemplified by the protocol of
mid-anneal pausing [15]. It is an interesting future directionof research to test whether and how environmental effects
change our conclusions, e.g., using methods such as those inRefs. [38,39].
ACKNOWLEDGMENTS

We also showed that ARA uses tunneling in the semiclassical potential to avoid a trap that the corresponding classical (zero-temperature) spin-vector dynamics cannot avoid.
However, this trap is not an obstacle for the classical spinvector Monte Carlo algorithm, which uses thermal activation
to hop over the corresponding potential barrier as long as thetemperature is appropriately chosen.
In contrast to ARA, we found a negative result for IRA in
the context of the p-spin model. In order for the IRA protocol

This research is based upon work supported by the Officeof the Director of National Intelligence (ODNI), Intelligence
Advanced Research Projects Activity (IARPA), via the U.S.
Army Research Office through Contract No. W911NF-17-C0050. The views and conclusions contained herein are those of
the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressedor implied, of the ODNI, IARPA, or the U.S. Government.
The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding
any copyright annotation thereon. D.L. is additionally partially supported by a DOE/HEP QuantISED program grant
and the Quantum Machine Learning and Quantum Computation Frameworks for HEP (QMLQCF) through Award No.
DE-SC0019227. The work of H.N. is also supported partiallyby the JSPS through KAKENHI Grant No. 26287086.

[1] T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355 (1998).
[2] J. Brooke, D. Bitko, T. F. Rosenbaum, and G. Aeppli, Science
284, 779 (1999).
[3] G. E. Santoro, R. Martoňák, E. Tosatti, and R. Car, Science 295,
2427 (2002).
[4] G. E. Santoro and E. Tosatti, J. Phys. A 39, R393 (2006).
[5] A. Das and B. K. Chakrabarti, Rev. Mod. Phys. 80, 1061 (2008).
[6] S. Morita and H. Nishimori, J. Math. Phys. 49, 125210 (2008).
[7] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, arXiv:quantph/0001106.
[8] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and
D. Preda, Science 292, 472 (2001).

[9] T. Albash and D. A. Lidar, Rev. Mod. Phys. 90, 015002 (2018).
[10] A. Perdomo-Ortiz, S. E. Venegas-Andraca, and A. AspuruGuzik, Quant. Info. Proc. 10, 33 (2011).
[11] N. Chancellor, New J. Phys. 19, 023024 (2017).
[12] A. D. King, J. Carrasquilla, J. Raymond, I. Ozfidan, E.
Andriyash, A. Berkley, M. Reis, T. Lanting, R. Harris, F.
Altomare, K. Boothby, P. I. Bunyk, C. Enderud, A. Fréchette, E.
Hoskinson, N. Ladizinsky, T. Oh, G. Poulin-Lamarre, C. Rich,
Y. Sato et al., Nature (London) 560, 456 (2018).
[13] D. Ottaviani and A. Amendola, arXiv:1808.08721.
[14] D. Venturelli and A. Kondratyev, Quantum Mach. Intell. 1, 17(2019).

052321-10

DYNAMICS OF REVERSE ANNEALING FOR THE FULLY …

PHYSICAL REVIEW A 100, 052321 (2019)[15] J. Marshall, D. Venturelli, I. Hen, and E. G. Rieffel, Phys. Rev.
Appl. 11, 044083 (2019).
[16] M. Ohkuwa, H. Nishimori, and D. A. Lidar, Phys. Rev. A 98,
022314 (2018).
[17] P. Pfeuty, Ann. Phys. 57, 79 (1970).
[18] C. R. Laumann, R. Moessner, A. Scardicchio, and S. L. Sondhi,
Phys. Rev. Lett. 109, 030502 (2012).
[19] J. Tsuda, Y. Yamanaka, and H. Nishimori, J. Phys. Soc. Jpn. 82,
114004 (2013).
[20] T. Kato, J. Phys. Soc. Jpn. 5, 435 (1950).
[21] S. Jansen, M.-B. Ruskai, and R. Seiler, J. Math. Phys. 48,
102111 (2007).
[22] D. A. Lidar, A. T. Rezakhani, and A. Hamma, J. Math. Phys.
50, 102106 (2009).
[23] M. H. Amin, Phys. Rev. A 92, 052323 (2015).
[24] H. Munoz-Bauza, H. Chen, and D. Lidar, npj Quantum Inf. 5, 2
(2019).
[25] V. Karanikolas and S. Kawabata, arXiv:1806.08517.
[26] https://docs.dwavesys.com/docs/latest/c_fd_hg.html.
[27] T. Jörg, F. Krzakala, J. Kurchan, A. C. Maggs, and J. Pujos,
Europhys. Lett. 89, 40004 (2010).

[28] T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D.
Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer, Science
345, 420 (2014).
[29] T. Albash and D. A. Lidar, Phys. Rev. X 8, 031016 (2018).
[30] I. Hen, J. Job, T. Albash, T. F. Ronnow, M. Troyer, and D. A.
Lidar, Phys. Rev. A 92, 042325 (2015).
[31] J. A. Smolin and G. Smith, Front. Phys. 2, 52 (2014).
[32] T. Albash, T. Rønnow, M. Troyer, and D. Lidar, Eur. Phys. J.:
Spec. Top. 224, 111 (2015).
[33] S. W. Shin, G. Smith, J. A. Smolin, and U. Vazirani,arXiv:1401.7087.
[34] S. Muthukrishnan, T. Albash, and D. A. Lidar, Phys. Rev. X 6,031010 (2016).
[35] Y. Seki and H. Nishimori, Phys. Rev. E 85, 051112 (2012).
[36] R. D. Somma, D. Nagaj, and M. Kieferová, Phys. Rev. Lett.
109, 050501 (2012).
[37] L. T. Brady and W. van Dam, Phys. Rev. A 95, 032335 (2017).
[38] G. Passarelli, G. De Filippis, V. Cataudella, and P. Lucignano,Phys. Rev. A 97, 022319 (2018).
[39] G. Passarelli, V. Cataudella, and P. Lucignano, Phys. Rev. B
100, 024302 (2019).

052321-11