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

TL;DR

This paper investigates the dynamics of reverse annealing in the fully-connected $p$-spin model, demonstrating exponential speedup over traditional quantum annealing and comparing classical and quantum approaches.

Contribution

It provides a numerical analysis of reverse annealing dynamics, confirms equilibrium predictions, and compares different reverse annealing protocols for the $p$-spin model.

Findings

Reverse annealing offers exponential speedup over conventional quantum annealing.

Iterated reverse annealing is ineffective for the $p$-spin model.

Adiabatic reverse annealing with $h$-gain may improve performance.

Abstract

Reverse annealing is a relatively new variant of quantum annealing, in which one starts from a classical state and increases and then decreases the amplitude of the transverse field, in the hope of finding a better classical state 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 the predictions of equilibrium statistical mechanics. In particular, we corroborate the equilibrium analysis prediction that 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, but note that a recently-introduced protocol ("$h$-gain"), which implements adiabatic reverse annealing, may lead to improved performance.

Equations

$$\begin{array}{rl}\hat{H}(t)=& s(t){\hat{H}}_0+\text{quantity}(1-s(t))\text{quantity}(1-\lambda (t)){\hat{H}}_{\text{init}}\\ & +\Gamma \text{quantity}(1-s(t))\lambda (t){\hat{V}}_{\text{TF}},\end{array}$$

$${\hat{V}}_{\text{TF}}=-\sum _{i=1}^N{\hat{\sigma }}_i^x,$$

$${\hat{H}}_{\text{init}}$$

$${\hat{H}}_0=-N\text{quantity}(\frac{1}{N}\sum _{i=1}^N{\hat{\sigma }}_i^z{)}^p,$$

$${ϵ}_i=\{\begin{array}{ll}+1& \text{for}i\le Nc\\ -1& \text{for}i>Nc,\end{array}$$

$${\hat{S}}_1^{x,z}\equiv \frac{1}{2}\sum _{i=1}^{Nc}{\hat{\sigma }}_i^{x,z},{\hat{S}}_2^{x,z}\equiv \frac{1}{2}\sum _{i=Nc+1}^N{\hat{\sigma }}_i^{x,z}$$

$${\hat{H}}_0$$

$${\hat{H}}_{\text{init}}$$

$$\begin{array}{rl}& f=s(p-1)m^{p-1}\\ & -c\sqrt{{\left[sp{m}^{p-1}+(1-s)(1-\lambda )\right]}^2+{\Gamma }^2(1-s{)}^2{\lambda }^2}\\ & -(1-c)\sqrt{{\left[sp{m}^{p-1}-(1-s)(1-\lambda )\right]}^2+{\Gamma }^2(1-s{)}^2{\lambda }^2},\end{array}$$

$$\begin{array}{rl}& m=c\frac{sp{m}^{p-1}+(1-s)(1-\lambda )}{\sqrt{{\left[sp{m}^{p-1}+(1-s)(1-\lambda )\right]}^2+{\Gamma }^2(1-s{)}^2{\lambda }^2}}\\ & +(1-c)\frac{sp{m}^{p-1}-(1-s)(1-\lambda )}{\sqrt{{\left[sp{m}^{p-1}-(1-s)(1-\lambda )\right]}^2+{\Gamma }^2(1-s{)}^2{\lambda }^2}}.\end{array}$$

$$p_{\text{e}}(\tau )=1-\text{quantity}|⟨\varphi ||\psi (\tau )⟩{|}^2$$

$$\text{TTS}(\tau ,p_d)=\tau \frac{\log (1-p_d)}{\log p_e(\tau )}.$$

$$\begin{array}{rl}{V}_{\text{SC}}\equiv & ⟨\Omega (s)|\hat{H}(s)|\Omega (s)⟩\\ =& -sN{\left(n_1\sin {\theta }_1\cos {\varphi }_1+n_2\sin {\theta }_2\cos {\varphi }_2\right)}^p\\ & -(1-s)(1-\lambda )N(n_1\sin {\theta }_1\cos {\varphi }_1-n_2\sin {\theta }_2\cos {\varphi }_2)\\ & -\Gamma (1-s)\lambda N(n_1\cos {\theta }_1+n_2\cos {\theta }_2),\end{array}$$

$$n_1=\frac{Nc}{N},n_2=\frac{N-Nc}{N},$$

$$\begin{array}{rl}& |\Omega (s)⟩\\ & =\underset{i=1}{\overset{Nc}{⨂}}\text{quantity}[\cos \frac{{\theta }_1(s)}{2}{|-⟩}_i+\sin \frac{{\theta }_1(s)}{2}{e}^{i{\varphi }_1}{|+⟩}_i]\\ & \times \underset{i=Nc+1}{\overset{N}{⨂}}\text{quantity}[\cos \frac{{\theta }_2(s)}{2}{|-⟩}_i+\sin \frac{{\theta }_2(s)}{2}{e}^{i{\varphi }_2}{|+⟩}_i].\end{array}$$

$$\hat{H}(t)=A(t){\hat{H}}_0+B(t){\hat{V}}_{\text{TF}}$$

$${\pi }_j^{(r)}=\sum _i(P^r{)}_{ji}{\pi }_i^{(0)},$$

Full text

I Introduction

Quantum annealing (QA) is a quantum-mechanical metaheuristic for combinatorial optimization problems Kadowaki and Nishimori (1998); Brooke et al. (1999); Santoro et al. (2002); Santoro and Tosatti (2006); Das and Chakrabarti (2008); Morita and Nishimori (2008), and its strict adiabatic realization is known as adiabatic quantum computation Farhi et al. (2000, 2001); Albash and Lidar (2018a). 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 some time ago Perdomo-Ortiz et al. (2011) to instead start from an appropriately chosen classical state and gradually increase and then decrease the amplitude 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 showed by 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. Chancellor (2017), and the protocol whereby the initial state is classical rather than a quantum superposition is now called reverse annealing. The current generation of the D-Wave quantum annealer, D-Wave 2000Q, implements reverse annealing, and several recent studies have been reported on using this feature in different contexts, e.g., quantum simulation King et al. (2018), matrix factorization Ottaviani and Amendola (2018), portfolio optimization Venturelli and Kondratyev (2018), and mid-anneal pausing Marshall et al. (2019).

A theoretical study of reverse annealing was initiated recently by some of us in Ref. Ohkuwa et al. (2018), in which equilibrium statistical mechanics was employed to study when and how reverse annealing is effective. It was shown explicitly that reverse annealing can turn first order quantum phase transitions into second order transitions in a simple problem, the $p$-spin model, by choosing an appropriate annealing path. It is generally the case that a first-order quantum phase transition is associated with an exponentially closing energy gap $\Delta$ between the ground state and the first excited state as a function of system size, whereas a second order transition is associated with a polynomially closing gap (though exceptions are known Pfeuty (1970); Laumann et al. (2012); Tsuda et al. (2013)). On the other hand, according to the adiabatic theorem of quantum mechanics, the computation time $\tau$ is inversely proportional to a small power of the energy gap $\Delta$ Kato (1950); Jansen et al. (2007); Lidar et al. (2009). Taken together, these facts suggest that reverse annealing provides an exponential speedup relative to the conventional QA protocol, for the $p$-spin model. In this work we study this intriguing connection from the dynamical perspective, by numerically solving the time-dependent Schrödinger equation of the reverse annealing protocol for the $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 non-trivial 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 Amin (2015), results from equilibrium statistical mechanics 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 non-equilibrium features such as coherent oscillations due to quantum interference Munoz-Bauza et al. (2019); Karanikolas and Kawabata (2018).

In the present work we evaluate the unitary, closed-system dynamics of reverse annealing by direct numerical integration of the time-dependent Schrödinger equation for large system sizes. This is possible due to a special symmetry of the $p$-spin model. To this end, it is useful to distinguish two types of reverse annealing. The first is the one adopted in Ref. Ohkuwa et al. (2018) to facilitate analytical treatment by equilibrium statistical mechanics. In this formulation, an additional term is introduced into the Hamiltonian, which is usually composed of just two terms representing the classical cost function and the transverse field. The additional, third term serves to enforce a given classical ground state as an initial condition. The system is then expected to follow an adiabatic path to the final state. We call this protocol adiabatic reverse annealing (ARA). In the second formulation the initial state is a classical state, an excited state of the initial (classical) Hamiltonian (the cost function), as formulated in Refs. Perdomo-Ortiz et al. (2011); Chancellor (2017). 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 is reverse annealing as is widely used in practice on the D-Wave device King et al. (2018); Ottaviani and Amendola (2018); Venturelli and Kondratyev (2018); Marshall et al. (2019), and is often used iteratively in practice, that is, the final state of a single cycle is fed into the system as its initial state of the next cycle under the expectation that the result improves iteratively. We therefore call this protocol iterated reverse annealing (IRA).

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 numerically for 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 the instantaneous Hamiltonian with an emphasis on the special structure of the $p$-spin model. The last section is devoted to conclusions.

II Adiabatic reverse annealing

In adiabatic reverse annealing (ARA), we add a term $\hat{H}_{\text{init}}$ to the Hamiltonian to fix the initial state to the ground state of the initial Hamiltonian Ohkuwa et al. (2018),

[TABLE]

where $\hat{H}_{0}$ is the cost function to be minimized and is a function only of the Pauli matrices $\{\hat{\sigma}_{i}^{z}\}$. The last term represents the transverse field,

[TABLE]

where $N$ is the total number of sites (qubits). The parameter $\Gamma$ controls the strength of the transverse field relative to the other terms. The initialization term is written as

[TABLE]

where $\{\epsilon_{i}(=\pm 1)\}$ denotes the given classical initial state. This initial Hamiltonian fixes $\hat{\sigma}_{i}^{z}$ to $\epsilon_{i}$ as its ground state. The time-dependent parameters $s(t)$ and $\lambda(t)$ both change from 0 to 1 as time $t$ proceeds from 0 to $\tau$, where $\tau$ is the computation time. Thus the initial Hamiltonian is $\hat{H}_{\text{init}}$ and the final one is $\hat{H}_{0}$. The quantum term $\hat{V}_{\text{TF}}$ is effective only in intermediate times. The existence of the term $\hat{H}_{\text{init}}$ allows us to analyze the system properties by equilibrium statistical mechanics under the assumption of adiabatic evolution.111 A recent update of the D-Wave quantum annealer introduced a function ‘$h$-gain’, which can be used to realize ARA hga .

We study the $p$-spin model in this paper,

[TABLE]

where $p\geq 3$ is an integer. The ground state of this model is trivial, $\sigma_{i}^{z}=1~{}\forall i$ (and $\sigma_{i}^{z}=-1~{}\forall i$ for $p$ even). This simple problem is known to be hard for conventional quantum annealing due to a first order phase transition Jörg et al. (2010). Since $\hat{H}_{0}$ and $\hat{V}_{\text{TF}}$ are symmetric with respect to an arbitrary permutation of the site index $i$, we do not lose generality by the following assignment of $\epsilon_{i}$,

[TABLE]

where the probability $c~{}(1/2\leq c\leq 1)$ for $\epsilon_{i}$ to be $+1$ is chosen 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, $\sigma_{i}^{z}=1,~{}\forall i$.222For $p$ even, another degenerate ground state $\sigma_{i}^{z}=-1~{}(\forall i)$ exists. But we do not lose generality by restricting ourselves to the subspace with non-negative magnetization.

Permutation symmetry allows us to write the Hamiltonian in terms of only a set of total spin operators of two subsystems,

[TABLE]

as

[TABLE]

Since $(\hat{\bf{S}}_{1})^{2}$ and $(\hat{\bf{S}}_{2})^{2}$ commute with the total Hamiltonian, the system stays in the subspace with the largest values of $(\hat{\bf{S}}_{1})^{2}$ and $(\hat{\bf{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$, $\mathcal{O}(N^{2})$.

II.1 Static properties

We first review the results presented in Ref. Ohkuwa et al. (2018) with some generalizations. A statistical-mechanical analysis of the $p$-spin model for the case of ARA of Eq. (1) with $s$ and $\lambda$ taken as fixed parameters leads to the following expression for the free energy at zero temperature (the ground state):

[TABLE]

where $m$ satisfies the following self-consistency condition to minimize the free energy:

[TABLE]

Derivations are detailed in Appendix A of Ref. Ohkuwa et al. (2018), where $\Gamma$ was fixed to 1, but a generalization to arbitrary $\Gamma$ as given here is straightforward. We focus our analysis on the case of $p=3$, but other values of $p(\geq 4)$ show qualitatively the same behavior. Figure 1 is the $\lambda$-$s$ phase diagram for $\Gamma=1$.

The axis $\lambda=1$ corresponds to conventional QA, where $s$ is changed from 0 to 1. Along this axis, a first order phase transition at $s\simeq 0.4$ marks where QA becomes inefficient, since the energy gap closes exponentially there. As $\lambda$ is reduced from 1 toward the inner part of the phase diagram, the effects of the initial Hamiltonian $\hat{H}_{\text{init}}$ come into play, and the line of first order phase transitions is split into two parts, one for larger $\lambda$ and the other for smaller $\lambda$, if $c$ is above a threshold, i.e., if the initial condition is reasonably close to the final answer. This means that a path exists that connects the initial point $\lambda=s=0$ and the goal $\lambda=s=1$ without crossing a first-order transition. In particular, when we choose the path $s=\lambda$, the diagonal of the phase diagram, the system 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 the initial 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 be checked by other methods. Indeed, while usually a first (second) order transition is accompanied by an exponentially (polynomially) closing gap, as mentioned in Introduction there exist examples in which the gap closes non-exponentially at a first order phase transition Pfeuty (1970); Laumann et al. (2012); Tsuda et al. (2013). We therefore evaluated the energy gap by direct numerical diagonalization, which is possible for relatively large system sizes due to the high symmetry of the problem. As seen in Fig. 2, which is a log-log plot of the gap as a function of system size for $s=\lambda$, only the data for $c=0.9$ behave non-exponentially. This conforms with the standard expectation from the phase diagram of Fig. 1, where the path $s=\lambda$ avoids the line of first order transitions only when $c=0.9$. It is interesting that the data for $c=0.8$ behaves almost polynomially up to about $N=50$. This should reflect the fact that the line $s=\lambda$ 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 the energy barrier between two coexisting states should be small, and the system’s behavior is close to that at a second-order transition, as long as the size is not too large.

Larger values of $\Gamma$ change the phase diagram quantitatively, if not qualitatively. As seen in Figs. 3(a) ($\Gamma=2$), 3(b) ($\Gamma=5$), and 3(c) ($\Gamma=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 $\Gamma$ increases. This is to be expected, because for large $\Gamma$ the intermediate Hamiltonian is mostly dominated by the transverse-field term $\hat{V}_{\text{TF}}$, and is close to the initial Hamiltonian of conventional QA. It then follows that the first order transition of conventional QA persists even for relatively large values of $c$.

II.2 Schrödinger dynamics

We next report our results on the closed system dynamics. Let $\ket{\phi}$ denote the ground state of the cost function $\hat{H}_{0}$ and $\ket{\psi(\tau)}$ the actual state reached after time $\tau$, subject to Schrödinger dynamics. Then

[TABLE]

is the error probability, i.e., the probability of not reaching the ground state. This quantity is plotted in Fig. 4 as a function of system size. For $\Gamma=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 Rønnow et al. (2014)

[TABLE]

The TTS is the effective time it takes using “runs” lasting time $\tau$ to find the correct solution at least once with probability $p_{d}$, which we set to $0.99$ or higher.

Results for the TTS are depicted in Fig. 5 as a function of $\tau$ for $\Gamma=1$ and $\Gamma=2$, and fixed system size. For a weaker transverse field $\Gamma=1$, ARA does not necessarily perform better than conventional QA. For a larger $\Gamma=2$, ARA has a shorter TTS than QA for the parameter values we have tested.

Interestingly, minima exist in Fig. 5 for $\Gamma=2$. This in fact holds for a broad range of system sizes, which allows us to extract the optimal annealing time $\tau$, $\tau_{\text{opt}}=\text{argmin}_{\tau}\text{TTS}(\tau,p_{d})$ as a function of $N$, and hence the optimal TTS scaling with size Rønnow et al. (2014); Albash and Lidar (2018b). The optimal TTS, estimated from the data for $c=0.8$ as in Fig. 5, is plotted in Fig. 6 as a function of system size. It is clearly seen that ARA with $\Gamma=2,~{}c=0.8$ exhibits polynomial scaling, whereas ARA with $\Gamma=1,~{}c=0.8$ and conventional QA are exponential  333Since this scaling is derived without the existence of an optimal annealing time, the true scaling can only be worse, as shown in Ref. Hen et al. (2015). Essentially, the reason is that the optimal TTS at small sizes $N$ is smaller than the TTS shown.. This is consistent with the static phase diagram of Figs. 1 and 3(a), where the path $s=\lambda$ does not cross a first order phase transition if $\Gamma=2,~{}c=0.8$ in Fig. 3(a), but conventional QA along $\lambda=1$ and ARA with $\Gamma=1,~{}c=0.8$ in Fig. 1 do. This is quite non-trivial because the TTS is a purely dynamical measure for finite-size systems, whereas the static phase diagram represents the long-time and thermodynamic limits, $\tau\to\infty,~{}N\to\infty$.

II.3 Classical spin vector dynamics

It is instructive to compare the results of the previous subsections with those of a counterpart classical model. In the spin-vector dynamics (SVD) model Smolin and Smith (2014); Albash et al. (2015); Shin et al. (2014); Muthukrishnan et al. (2016), the system evolves according to the classical Hamilton dynamics under the semi-classical potential $V_{\rm SC}$ corresponding to the Hamiltonian of Eq. (1),

[TABLE]

where

[TABLE]

and $\ket{\Omega(s)}$ is the spin-coherent state,

[TABLE]

The initial conditions are $\phi_{1}=\phi_{2}=0$ and $\theta_{1}=-\theta_{2}=\pi/2$. We used the equation of motion, which is derived from the path integral formulation and permutation symmetry Muthukrishnan et al. (2016).

The time evolution of the magnetization $\bra{\Omega(s)}{\frac{2}{N}(\hat{S}_{1}^{z}+\hat{S}_{2}^{z})}\ket{\Omega(s)}$ under the SVD model is compared 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) the quantum dynamics succeeds in coming close to the right answer $m=1$ at the end, but the classical dynamics fails.

To explain this behavior, we note that for the parameter values of Fig. 7(c) [which lies between what is shown in Figs. 3(a) and 3(b)], the evolution path crosses the first order phase transition line near its termination in the middle of the $\lambda$-$s$ phase diagram, where the energy barrier across the first order 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 appearance of an energy barrier in the semiclassical potential. Figure 9 conspicuously demonstrates the marked difference between classical and quantum dynamics, which is not clear in the static phase diagram. It shows the final value of the magnetization as a function of $\Gamma$. Quantum dynamics exhibits a gradual deterioration as $\Gamma$ increases beyond $3.4$, corresponding to a lower tunneling rate, whereas the classical SVD suddenly fails beyond a threshold close to $\Gamma=3.4$, where the energy barrier seen in Fig. 8 appears.

II.4 Classical spin vector Monte Carlo

Since the energy barrier is thin and low for Fig. 7(c), the classical system can hop over the barrier if we introduce thermal fluctuations. Figure 10 shows the results of spin vector Monte Carlo Shin et al. (2014), a finite-temperature version of SVD, in which the angle variables in the semiclassical potential are updated stochastically according to the usual Monte Carlo rule.

As seen in this figure, the system succeeds in reaching a 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. Clearly, the temperature needs to be tuned in order to reach the correct solution. This is simple to do in the present case, but in general one does not know the right value in advance. This is in contrast with the quantum dynamics, where there is no such adjustable parameter.

III Iterated Reverse annealing

In iterated reverse annealing (IRA), one measures the state in the computational (classical) basis after a single cycle of RA and starts the next cycle from the classical state thus obtained. Let $\hat{U}$ denote the time evolution operator under the Hamiltonian of reverse annealing without the initialization term present in ARA,

[TABLE]

where $A(t)=s(t),B(t)=1-s(t)$ with an appropriate choice of the function $s(t)$. The key difference from standard QA is that now $s(t)$ is non-monotonic and has a minimum, $s_{\min}$, as shown in Fig. 11.

The value of this minimum relative to the critical value $s_{c}$, plays an important role, since it determines whether or not the system crosses a phase transition: it does when $s_{\min}<s_{c}$.

Time changes from $t=0$ to $t=\tau$, where $\tau$ is the annealing time of a single cycle, and $s(0)=s(\tau)=1$. The total annealing time with $r$ cycles is $r\tau$. The initial state is classical, i.e., a computational basis state $\ket{i}$ (where $i\in[0,\dots,2^{N}-1]$ is the state index) and the probability of the final state being $\ket{j}$ after a cycle is $p_{ji}(\tau)=\sum_{k\in S(e_{j})}|\bra{k}\hat{U}(\tau)|i\rangle|^{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 $\hat{U}$ is replaced by open quantum system evolution, which means that thermal effects play a role as well during each cycle in the real device.

IRA is an algorithm that repeats the cycle of RA and measurement. Therefore, IRA can be considered as a series of classical Markov transitions, and a single cycle is a stochastic transition to the final classical state from the initial classical state. Using the overlap $p_{ji}(\tau)$ above, the classical transition probability matrix $P$ has elements $P_{ji}=p_{ji}(\tau)$. Let $\pi$ denote a classical probability vector whose $i$th component $\pi_{i}$ denotes the probability of obtaining the $i$th computational basis state. After $r$ cycles, the probability becomes

[TABLE]

where $\pi_{i}^{(0)}$ is the given initial state. The probability of obtaining the ground state is $\pi_{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 specific schedule choice we tested.

We choose a quadratic function of $s(t)$ with a minimum $s_{\min}$ 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$, $\tau=10$ and $p=3$. In Fig. 12(a), since $s_{\min}=0.5$ is above the critical value $s_{c}$,444 See Fig. 3 (left panel) of Ref. Seki and Nishimori (2012). and hence the system does not undergo a phase transition, the magnetization stays close to the initial value (given by $c$), though it is shifted somewhat towards $m=0$. If $s_{\min}$ 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\in[0.3,0.7]$. However, if the initial state is the ground state, the system seems to stay near the initial state.

To better understand what happens in Fig. 12, we have calculated the energy spectrum and occupation probability. Results are shown in Fig. 13 for $s_{\min}=0.3$.

Because of the symmetry of the annealing schedule of IRA with respect to $s=\frac{1}{2}$, avoided crossings in the first half appear also in the second half. Therefore, when annealing starts from the 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. It can also happen that due to the large number of avoided crossings in the first and second halves, an excited state can be further excited to higher energy levels, as in Fig. 13(b). The situation is far from simple, and the general tendency is that the distribution among states becomes broader as cycles are iterated except for the limiting case of adiabatic annealing, where the total evolution time $\tau$ is large on the scale set by the inverse of the minimum gap between the ground and first excited state.

The density plots in Figs. 14 and 15 show the transition probabilities $(P^{r})_{ji}$ from the initial state $i$ (horizontal axis) to the final state $j$ (vertical axis) in single cycles ($r=1$) and after multiple cycles ($r=3$ and $5$) with $s_{\min}=0.5$ (Fig. 14) and $s_{\min}=0.3$ (Fig. 15). For example, in Fig. 14(a), if we start from the ground state $i=0$, the system reaches the same ground state ($j=0$) with a large probability, as indicated by the yellow box, since the computation time $\tau=30$ is relatively long and close to adiabatic. Other states along the column $i=0$ have lower probabilities (dark colors). Similarly, if the initial condition is an excited state $i>0$, the same state $j=i$ is reached with a high probability, resulting in a series of bright colors along the diagonal. Repetition of IRA (panels (b) and (d)) or shorter computation time (panels (c) and (d)) results in scattering of the probabilities away from the diagonal. The situation is similar for $s_{\min}=0.3<s_{c}$ in Fig. 15, but with a slightly different structure. In any case, the probability 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 under coherent 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 expected that thermal relaxation will reduce excited state occupation probabilities, so that an open system study of IRA is likely to find better performance for the finite (but low) temperature $p$-spin model than reported here.

IV Conclusions

Reverse quantum annealing protocols have generated much interest recently, as they may potentially overcome some of the obstacles that cause conventional forward quantum annealing to fail. In this work we presented a numerical study of the 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, we demonstrated that the dynamical behavior of the system is consistent 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 a path that avoids a first-order quantum phase transition that is 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 that the lower bound for the TTS of QA is exponential. This is remarkable 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.

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 spin vector Monte Carlo algorithm, which uses thermal activation to hop over the corresponding potential barrier as long as the temperature 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 to work and provide an enhanced probability of finding the ground state after multiple cycles, the probability distribution of the final state after a single cycle should shift toward lower energy states than the initial state. Unfortunately, this condition is not satisfied in the $p$-spin model as is clear from Figs. 14 and 15. The generality of this conclusion is currently 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 the $p$-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 Somma et al. (2012); Muthukrishnan et al. (2016); Brady and van Dam (2017). However, our $p$-spin model result suggests the important 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 system with state transitions affected by thermal effects. In the D-Wave device, for example, thermal relaxation seems to lead to better performance as exemplified by the protocol of mid-anneal pausing Marshall et al. (2019). It is an interesting future direction of research to test whether and how environmental effects change our conclusions, e.g., using methods such as those in Refs. Passarelli et al. (2018, 2019).

acknowledgement

This research is based upon work supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office contract W911NF-17-C-0050. 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 expressed or 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 DOE/HEP QuantISED program grant, QMLQCF (Quantum Machine Learning and Quantum Computation Frameworks for HEP), award number DE-SC0019227.