quantum physics, quantum information, computer science

Tadashi Kadowaki

# Greedy parameter optimization for diabatic quantum annealing

Tadashi Kadowaki1 and Hidetoshi Nishimori2,3,4 1DENSO CORPORATION, Kounan,
Minato-ku, Tokyo 108-0075, Japan
2Institute of Innovative Research, Tokyo Institute of Technology, Yokohama,
Kanagawa 226-8503, Japan
3Graduate School of Information Sciences, Tohoku University, Sendai, Miyagi
980-8579, Japan
4RIKEN Interdisciplinary Theoretical and Mathematical Sciences Program
(iTHEMS), Wako, Saitama 351-0198, Japan<EMAIL_ADDRESS>
###### Abstract

A shorter processing time is desirable for quantum computation to minimize the
effects of noise. We propose a simple procedure to variationally determine a
set of parameters in the transverse-field Ising model for quantum annealing
appended with a field along the $y$ axis. The method consists of greedy
optimization of the signs of coefficients of the $y$-field term based on the
outputs of short annealing processes. We test the idea in the ferromagnetic
system with all-to-all couplings and spin-glass problems, and find that the
method outperforms the traditional form of quantum annealing and simulated
annealing in terms of the success probability and the time to solution, in
particular in the case of shorter annealing times, achieving the goal of
improved performance while avoiding noise. The non-stoquastic $\sigma^{y}$
term can be eliminated by a rotation in the spin space, resulting in a non-
trivial diabatic control of the coefficients in the stoquastic transverse-
field Ising model, which may be feasible for experimental realization.

###### keywords:

quantum annealing, counterdiabatic driving, greedy optimization

##  1 Introduction

Quantum optimization algorithm is an active field of research, including
quantum annealing (QA) [1, 2, 3, 4, 5, 6, 7, 8], adiabatic quantum computing
[9], and quantum approximate optimization algorithm (QAOA) [10]. If QA
operates in the adiabatic regime, it often encounters a difficulty of
exponential computation time coming from exponential closing of the energy gap
due to a first-order phase transition, examples of which are found in Refs.
[11, 12, 13]. A number of approaches have been tried to mitigate this problem,
including non-stoquastic Hamiltonians [14, 15, 16, 17], inhomogeneous field
driving [18, 19, 20, 21, 22, 23], reverse annealing and pausing [24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34], and counterdiabatic driving [35, 36, 37, 38,
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49].

The last method of counterdiabatic driving is particularly appealing because
it allows us to reduce the computation time as compared to the adiabatic
process, yet keeping the success probability large. The reduction of
computation time is important not just theoretically but also experimentally
from the viewpoint of avoiding the effects of noise in quantum devices. One of
the drawbacks of the idea of counterdiabatic driving is the excessive
complexity of the additional counterdiabatic term, see e.g. Ref. [40], which
generally makes it impossible to be implemented experimentally. An interesting
development to avoid this problem is the variational approach to optimize the
coefficients of approximate counterdiabatic terms composed of manageable
operators [45, 44, 46, 50, 49, 51]. In the present paper, we further simplify
this strategy and propose to determine just the signs of coefficients of
counterdiabatic terms composed of $y$ components of the Pauli matrix. Since
the determination of signs is processed in a greedy way variationally, we call
the method quantum greedy optimization (QGO). This approach is appealing
because of its simplicity and yet a significantly improved success probability
in comparison with the original QA and the classical method of simulated
annealing.

This paper is organized as follows. Section 2 defines the model system and
explains our strategy to solve optimization problems. The performance of the
method is analyzed numerically in section 3 for the mean-field, ferromagnetic,
and spin glass problems, and the results are compared with simulated annealing
(SA) [52] and the original QA. Also discussed are a further simplification of
the algorithm and a possible improvement based on an oracle setting. The final
section 4 summarizes and discusses the results.

##  2 Formulation

Let us define the following time-dependent Hamiltonian

$\displaystyle\mathcal{H}=A(t)\mathcal{H}^{z}+B(t)\mathcal{H}^{x}+\sum_{i=1}^{N}C_{i}(t)\mathcal{H}_{i}^{y},$
(1)

where $\mathcal{H}^{z}$, $\mathcal{H}^{x}$, and $\mathcal{H}_{i}^{y}$ are the
problem (Ising) Hamiltonian, the transverse field term (to be called the
$x$-field), and the $y$-field term, respectively,

$\displaystyle\mathcal{H}^{z}=-\sum_{i<j}J_{ij}\sigma_{i}^{z}\sigma_{j}^{z},~{}\mathcal{H}^{x}=-\sum_{i=1}^{N}\sigma_{i}^{x},~{}\mathcal{H}_{i}^{y}=-\sigma_{i}^{y}.$
(2)

The symbols $\sigma_{i}^{x},\sigma_{i}^{y}$, and $\sigma_{i}^{z}$ denote the
components of the Pauli matrix at site (qubit) $i$, and $J_{ij}$ is for the
interaction between sites $i$ and $j$. The problem size is $N$. The term
$\mathcal{H}_{i}^{y}$ renders the system non-stoquastic [53] and was
introduced in Ref. [45] to approximately realize counterdiabatic driving. See
also Refs. [44, 48, 50, 54]. The time dependent coefficients $A(t),B(t)$, and
$C_{i}(t)$ should satisfy the initial ($t=0$) and final $(t=\tau)$ conditions,

$\displaystyle A(0)=C_{i}(0)=0,B(0)\neq 0,~{}A(\tau)>0,B(\tau)=C_{i}(\tau)=0$
(3)

such that the system starts with the simple transverse field term
$\mathcal{H}^{x}$ and ends with the problem Hamiltonian $\mathcal{H}^{z}$. The
$y$-field is applied only in the middle of annealing. We arbitrarily choose
the following functions for these coefficients for simplicity though other
(similar) functions can be considered [45],

$\displaystyle
A(t)=\frac{at}{\tau},~{}B(t)=b\Big{(}1-\frac{t}{\tau}\Big{)},~{}C_{i}(t)=c_{i}\,\sin^{2}\Big{(}\frac{\pi
t}{\tau}\Big{)}$ (4)

with $a=1$ to set the overall energy scale. The annealing time $\tau$ can be
chosen arbitrarily and we test the cases of $\tau=1$ and $\tau=5$ in the
following.

Our strategy is to optimize $b$ and $c_{i}$ under a variational principle to
optimize appropriate measures. We simplify the process by fixing the amplitude
of $c_{i}$ to a value found optimal in the mean-field case (to be described
below) and choosing only the sign of $c_{i}$ in a greedy way. This facilitates
the process considerably and yet will turn out to lead to significant
performance improvements.

Our variational optimization is carried out by minimization of the expectation
value of the final energy
$E=\braket{\psi(t=\tau)}{\mathcal{H}^{z}}{\psi(t=\tau)}$, where $\psi(t)$ is
the wave function at time $t$, as well as by maximization of the ground-state
probability, or the fidelity, $P_{\rm gs}=|\braket{\psi(t=\tau)}{\psi_{\rm
gs}}|^{2}$, where $|\psi_{\rm gs}\rangle$ is the true ground state of the
Ising Hamiltonian. Since we do not know the true ground state in a generic
problem, we refer to the latter measure as an oracle.

The additional $y$-field term may not be straightforward to be implemented
experimentally. Nevertheless, by a simple rotation of the spin axes [45, 50],

$\displaystyle
U_{g}(t)=\exp\left(\frac{i}{2}\sum_{i}\theta_{i}(t)\sigma_{i}^{z}\right)$ (5)

with $\theta_{i}(t)=\arctan(C_{i}(t)/B(t))$, we can transform the Hamiltonian
into a transverse-field Ising model with non-trivial time dependence of
coefficients,

$\displaystyle\mathcal{H}_{\rm
eff}=A(t)\mathcal{H}^{z}-\frac{1}{2}\sum_{i}\frac{d\theta_{i}(t)}{dt}\,\sigma_{i}^{z}-\sum_{i}\sqrt{B(t)^{2}+C_{i}(t)^{2}}\,\sigma_{i}^{x}.$
(6)

which may be easier to realize experimentally than the original Hamiltonian of
Eq. (1). An example of the behavior of the coefficients in Eq. (6) is shown in
Appendix A.

##  3 Numerical results

In this section, we present numerical results for the mean-field,
ferromagnetic, and random (spin glass) systems. In addition, we discuss an
additional algorithm which shortcuts the process of iterative determination of
the signs of coefficients.

###  3.1 Mean-field theory

As a preliminary to the next section of a ferromagnetic system, we first
analyze the properties of the mean-field theory with the Hamiltonian (see Ref.
[44] for a related idea)

$\mathcal{H}=-t\braket{\sigma^{z}}\sigma^{z}-b(1-t)\sigma^{x}-c\sin^{2}(\pi
t)\sigma^{y},$ (7)

where $\braket{\sigma^{z}}$ is the magnetization of the system
$\braket{\psi(t)}{\sigma^{z}}{\psi(t)}$. We set $\tau=1$ in this section. The
state vector can be expressed as a qubit
$\ket{\psi(t)}=\alpha\ket{0}+\beta\ket{1}$, where $\ket{0}=(1,0)^{T}$ and
$\ket{1}=(0,1)^{T}$, and the Schrödinger equation reads

$\frac{d}{dt}\left(\begin{array}[]{c}\alpha\\\
\beta\end{array}\right)=\left(\begin{array}[]{cc}it\braket{\sigma^{z}}&ib(1-t)+c\sin^{2}(\pi
t)\\\ ib(1-t)-c\sin^{2}(\pi
t)&-it\braket{\sigma^{z}}\end{array}\right)\left(\begin{array}[]{c}\alpha\\\
\beta\end{array}\right),$ (8)

where $\braket{\sigma^{z}}$ is $|\alpha|^{2}-|\beta|^{2}$. As the ground state
of the final Hamiltonian is doubly degenerate, we choose the spin-up state
$\braket{\sigma^{z}}>0$ as the desired ground state. Two measures of success
as defined in the previous section, the fidelity $P_{\rm
gs}=(\braket{\sigma^{z}}+1)/2$ and the energy $E=-\braket{\sigma^{z}}^{2}$,
can be expressed by the magnetization $\braket{\sigma^{z}}$ in the mean-field
theory, which is not the case generally.

Figure 1 shows the dependence of magnetization $\braket{\sigma^{z}(t=1)}$ on
$b$ and $c$. In this landscape of the parameter space, we find the optimal
values to maximize $\braket{\sigma^{z}}$ by the classical optimizer, the
Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [55, 56, 57, 58], as $b_{\rm
opt}^{\rm mf}=0.539$ and $c_{\rm opt}^{\rm mf}=1.565$ with
$\braket{\sigma^{z}}=1.000$. Notice that no greedy search is applied here
since the problem is very simple in the mean-field case.

Figure 1:  Parameter dependence of the magnetization $\braket{\sigma^{z}}$ for
the mean-field theory. The optimal point $(b_{\rm opt}^{\rm mf}=0.539,c_{\rm
opt}^{\rm mf}=1.565)$ is marked in an open asterisk.

Although finding optimal values of the coefficients is sufficient for our
purpose in the following analysis, we illustrate in Appendix B the time
developments of the wave function, the coefficients, and the magnetization in
the mean-field context, comparing mean-field results with the exact solution
of the counterdiabatic term, for better understanding of the physics of the
mean-field theory.

###  3.2 Ferromagnetic system

####  3.2.1 Parameter dependence of the energy and fidelity

We move a step forward and analyze the ferromagnetic system with all-to-all
interactions,

$\displaystyle\mathcal{H}^{z}=-\frac{J}{N-1}\sum_{i<j}\sigma_{i}^{z}\sigma_{j}^{z}.$
(9)

The mean-field theory of the previous section is expected to give the correct
solution to this model in the thermodynamic limit [59] though finite-size
effects may reveal differences.

As the interactions are homogeneous, we drop the $i$ dependence of $c_{i}$ and
set all of them to $c$. To break the $\mathbb{Z}_{2}$ symmetry, we choose
positive values of $b$ and $c$ and focus our attention on the region around
the optimal value obtained in the mean-field case, $b\in[0,1]$ and
$c\in[1,2]$. Dependence of the fidelity and energy on $b$ and $c$ is shown in
Fig. 2 (a) and (b), respectively, for $N=20$.

\begin{overpic}[height=91.04881pt]{fig2a.pdf} \put(0.0,75.0){(a)}
\end{overpic}\begin{overpic}[height=91.04881pt]{fig2b.pdf} \put(0.0,78.5){(b)}
\end{overpic}

\begin{overpic}[height=91.04881pt]{fig2c.pdf} \put(0.0,75.0){(c)}
\end{overpic}\begin{overpic}[height=91.04881pt]{fig2d.pdf} \put(0.0,78.5){(d)}
\end{overpic}

Figure 2:  Parameter ($b,c$) dependence of (a) the fidelity and (b) the energy
in the ferromagnetic model with $N=20$. Dependence on the annealing time and
the parameter $c$ of (c) the fidelity and (d) the energy when annealing is
terminated before the designated annealing time ($\tau_{\rm actual}<\tau$) and
is carried through beyond ($\tau_{\rm actual}>\tau$). Color scales are located
on the right side of each figure.

We also conducted the same analysis for the case where annealing is terminated
before the designated annealing time ($\tau_{\rm actual}<\tau$) or carried
through beyond ($\tau_{\rm actual}>\tau$), and the results are displayed in
Fig. 2 (c) and (d). These plots indicate that the optimization is robust in
$b$ and sensitive to $c$ and the actual annealing time. Similarities of Fig. 2
(a) and (b) and of Fig. 2 (c) and (d) suggest that the two measures of
success, the fidelity and the energy, give essentially the same result, as
anticipated from the mean-field case.

System size dependence of the optimal parameter values for each $N$, $b_{\rm
opt}^{N}$ and $c_{\rm opt}^{N}$, is shown in Fig. 3(a) and (b). They show
convergence beyond $N\approx 12$ toward the values at $N=20$ ($b_{\rm
opt}^{N=20}=0.539,c_{\rm opt}^{N=20}=1.564$). These values are in good
agreement with the results from the mean-field theory, $b_{\rm opt}^{\rm
mf}=0.539$ and $c_{\rm opt}^{\rm mf}=1.565$, as expected. The error in
magnetization $1-\braket{\sigma^{z}}$ decreases with the system size, and we
do not see a visible difference between the two measures in Fig. 3 (c) and
(d). The residual error decreases polynomially as a function of the system
size $N$. We will use the size-dependent parameters $b_{\rm opt}^{N}$ and
$c_{\rm opt}^{N}$ obtained from fidelity optimization for further analysis.

\begin{overpic}[height=113.81102pt]{fig3a.pdf} \put(0.0,70.0){(a)}
\end{overpic}\begin{overpic}[height=113.81102pt]{fig3b.pdf}
\put(0.0,70.0){(b)} \end{overpic}

\begin{overpic}[height=113.81102pt]{fig3c.pdf} \put(0.0,70.0){(c)}
\end{overpic}\begin{overpic}[height=113.81102pt]{fig3d.pdf}
\put(0.0,70.0){(d)} \end{overpic}

Figure 3:  (a) Size dependence of the parameter $b_{\rm opt}^{N}$ optimized
with respect to the fidelity and energy measures in the ferromagnetic model.
(b) Size dependence of the optimized parameter $c_{\rm opt}^{N}$. (c) Time
development of the magnetization, and (d) size dependence of the
magnetization.

####  3.2.2 Greedy optimization of the parameters

In preparation of the random case, we next try a step-wise greedy optimization
of parameters assuming site dependence of the sign of $c_{i}$ as in the well-
established methods like quadratic pseudo-boolean optimization (QPBO) and the
roof duality algorithm solutions [60, 61, 62, 63] and a sampling-based
algorithm [64]. More precisely, we choose the optimal values of $c_{i}$ one by
one with other coefficients fixed at their constant values. This approach
reduces the search space considerably and also makes the energy (or fidelity)
landscape simple.

The first step is to pick up an arbitrary $i$ (e.g., $i=1$) and optimize
$c_{1}$ by setting other $c_{j}~{}(j\neq 1)$ to zero. As shown in the rows of
Fig. 4 (a) and (b) for $N=8$, we choose optimal $c_{1}$ by minimizing
$1-P_{\rm gs}$ (in (a)) or $E$ (in (b)). See also Fig. C.1, where the scale of
the vertical axis is changed step by step for a better resolution. In
practice, when we use $E$ as the function to be minimized, this function is
symmetric with respect to the change of sign of $c_{1}$ due to the double
degeneracy of states. We therefore arbitrarily set $c_{1}$ to a positive value
$c_{1}^{\rm BFGS}$ by BFGS. When we use $1-P_{\rm gs}$ for minimization, we
choose the all-up state as the target ground state, so no problem of
degeneracy exists. We next choose another $i$ (e.g., $i=2$) and optimize
$c_{2}$ by fixing $c_{1}$ to the already-optimized value and other
$c_{j}~{}(j\neq 1,2)$ to zero. The $c_{2}$ dependence of $1-P_{\rm gs}$ and
$E$ is shown in the second row of Fig. 4. This process is repeated until all
parameters are fixed.

\begin{overpic}[width=199.16928pt]{fig4.pdf} \put(1.0,82.0){(a)}
\put(53.0,82.0){(b)} \end{overpic} Figure 4:  The parameter $c_{i}$ dependence
of (a) the fidelity and (b) the energy in each iteration (from top to bottom)
for the ferromagnetic model with $N=8$. Notice that the site index $i$ changes
from step to step as explained in the text. See also Fig. C.1 in which the
same data are shown but with different vertical scales for different iteration
steps for better resolution at each step.

It is seen in Fig. 4 that $1-P_{\rm gs}$ and $E$ become smaller as the
iteration proceeds.

We have found that the optimal values $c_{i}^{\rm BFGS}$ are always close to
$c_{\rm opt}^{N}$: The resulting optimal values vary from site to site,
ranging from $1.507$ to $1.566$, while $c_{\rm opt}^{N=8}$ is $1.563$
according to the analysis of the previous section. This means that our
sequential greedy optimization strategy finds the ground state with good
precision for the ferromagnetic system with all-to-all couplings.

###  3.3 Spin glass problem

We apply essentially the same procedure as developed above to the prototypical
hard optimization problem of spin glass.

The algorithm proceeds as follows. The value of $b$ is set to $b_{\rm
opt}^{N}$ obtained for the ferromagnetic case, and the absolute value of
$c_{i}$ is fixed to $c_{\rm opt}^{N}$ as demonstrated in the ferromagnetic
case. We only choose the sign of each $c_{i}$ step by step, site to site, in a
greedy way. It will be seen that this simple process leads to significant
improvements in performance as compared to the traditional QA and classical
SA.

The sign of $c_{i}$ is chosen by the average gradient of the optimization
measure ($E$ or $1-P_{\rm gs}$) near the origin $0\leq c_{i}\leq\Delta=0.1$
with the values of other $c_{j}$’s $(j\neq i)$ fixed. The reason is that the
optimization measure has only a single minimum, whose position can be detected
near the origin, as a function of a single $c_{i}$, as will be illustrated
later.

More precisely, we first choose an arbitrary $i$ (e.g. $i=8$) and fix all
other $c_{j}$ to zero. The optimization measure $E$ is symmetric with respect
to the inversion $c_{8}\to-c_{8}$, and we therefore look instead at the
curvature of $E$ in the asymmetric range $0\leq c_{8}\leq\Delta$. This process
is repeated for all $i$ and we choose the site $i$ with the largest absolute
average curvature. If this turns out to be $i=8$ with a positive derivative,
we set $c_{8}=-c_{\rm opt}^{N}$, which breaks the symmetry. We next choose
another $i$ (e.g. $i=6$) and see the sign of the average derivative of $E$
near the origin with $c_{8}$ fixed to $-c_{\rm opt}^{N}$ and all other
$c_{i}$’s to zero. We choose another $i$ and the same process is repeated for
all $i(=1,2,3,4,5,7)$ and we select the site with the largest value of the
absolute average derivative (e.g. $i=6$ with a negative derivative) and assign
the sign as $c_{6}=+c_{\rm opt}^{N}$. Iterating this process, all $c_{i}$’s
are assigned fixed values. See Algorithm 1.

Algorithm 1 Sequential QGO

0: QA measure $f(b,{\bm{c}})$, $b_{\rm opt}^{N}$, $c_{\rm opt}^{N}$

0: a solution of the cost function

1: $b\leftarrow b_{\rm opt}^{N}$

2: ${\bm{c}}\leftarrow(0,\cdots,0)$

3: repeat

4: ${\bm{g}}\leftarrow\nabla
f(b,{\bm{c}})=(\frac{f(b,c_{1}+\Delta,\cdots,c_{N})-f(b,{\bm{c}})}{\Delta},\cdots,\frac{f(b,c_{1},\cdots,c_{N}+\Delta)-f(b,{\bm{c}})}{\Delta})$

5: $i\leftarrow\operatorname*{arg\,max}_{j\in\\{j|c_{j}=0\\}}|g_{j}|$

6: $c_{i}\leftarrow-c_{\rm opt}^{N}\operatorname{sgn}g_{i}$

7: until $c_{i}\neq 0$ for all $i$

8: return $\operatorname{sgn}{\bm{c}}$

An example of this process is illustrated in Fig. 5 for the all-to-all
coupling Sherrington-Kirkpatrick model of spin glass [65] with the
distribution function of interaction,

$P(J_{ij})\sim\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{J_{ij}^{2}}{2\sigma^{2}}\right),$
(10)

with $\sigma^{2}=1/(N-1)$.

Figure 5:  Parameter dependence of the energy measure in each iteration of
sequential QGO in the random model ($N=8$) from the first step (top row) to
the last step (bottom row). Each column represents a site $i$, 1 (leftmost
column) to 8 (rightmost column). Blue-dashed and red-bold curves are
parameters for sites that are under evaluation in each iteration. In each
iteration, the site with the largest absolute average gradient marked in red-
bold is selected to fix the parameter. Gray-dotted curves located under the
red curves are additional plots for the site which we do not evaluate for
further iterations. See also Fig. C.2 in which the same data are shown with
different vertical scales for different iteration steps for a better
resolution.

Red curves show $E$ as a function of a given $c_{i}$ which has the largest
absolute value of the average gradient. Note that the procedure works for the
top row for the initial step as the average gradients are calculated in the
asymmetric range $[0,\Delta]$ and resulted in all positive even though the
curve is symmetric at $c_{i}=0$. Blue dotted curves are for other sites with
smaller gradient, and light gray curves show the behavior of $E$ when the
already-fixed $c_{i}$ is tentatively changed with other $c_{i}$’s fixed. One
sees that the minimum value of $E$ on red curves decreases as the greedy
optimization proceeds from the top row to the bottom row.

The final values of $c_{i}$ thus obtained are compared with those from direct
brute-force optimization by the BFGS algorithm,

$\displaystyle{\bm{c}}_{\rm QGO}$
$\displaystyle=(1.563,-1.563,-1.563,-1.563,1.563,-1.563,1.563,-1.563)$
$\displaystyle{\bm{c}}_{\rm BFGS}$
$\displaystyle=(1.570,-1.564,-1.562,-1.567,1.569,-1.561,1.561,-1.563).$

We observe that the difference is minimal, at most 0.4%, and the signs are
correctly reproduced by the greedy algorithm.

For a more systematic analysis, we applied the algorithm to the spin glass
with various sizes from $N=4$ to 16 with 100 random instances for each $N$.
The resulting success probability is plotted as a function of size in Fig.
6(a) for $\tau=1$ and 5. QGO and QA obey the Schrödinger dynamics, and the
calculation was conducted by QuTiP [66].

\begin{overpic}[height=128.0374pt]{fig6a.pdf} \put(0.0,70.0){(a)}
\end{overpic}\begin{overpic}[height=128.0374pt]{fig6b.pdf} \put(0.0,72.0){(b)}
\end{overpic}

Figure 6:  (a) Size dependence of the average success probability of QGO
(blue), SA (orange), and QA (green) for $\tau=1$ (circle) and $5$ (square) in
100 instances of the random system. (b) Size dependence of time to solution of
QGO (blue), SA (orange), and QA (green) for $\tau=1$ (circle) and $5$
(square). Dotted lines with open markers are plotted by taking into account
the factor $N(N+3)/2$ multiplying the original QGO to compensate for multiple
runs of QA. The shaded region represents 95% confidence intervals for each
line.

Also plotted are the success probabilities by the traditional QA and SA. The
latter classical algorithm has been tested by the following schedule of
temperature decrease (inverse temperature increase)

$\displaystyle\beta(t)=\frac{(t/\tau)}{1.1-(t/\tau)}$ (11)

where $\beta(0)=0$ and $\beta(\tau)=10$ and by solving the classical master
equation. QGO returns not a state vector in QA or a probability distribution
in SA but a specific solution of spin configuration, which is the output from
the Algorithm 1. Thus the success probability of QGO refers to a portion of
instances with success. Confidence intervals are calculated by bootstrapping
from 10,000 resamples. It is clear that the present algorithm QGO far
outperforms the other two.

Another viewpoint is provided by the time-to-solution (TTS), a standard
measure of computation time for heuristic algorithms [67]. It is defined as

$\displaystyle{\rm TTS}=\tau\frac{\log(1-P)}{\log(1-p_{\rm success})},$ (12)

where $p_{\rm success}$ is the empirical probability of success for the
computation time $\tau$, and $P$ is the target success probability often set
to 0.99. The result is plotted in Fig. 6 (b), which apparently shows a low
computational cost of the present algorithm QGO drawn in blue. We should
notice, however, that QGO carries an overhead of repeated computation of the
parameter dependence of $E$ at each iteration.

QGO calls QA as a subroutine $N(N+3)/2$ times in total to determine the
gradients. This number is derived as follows: In the $n$-th iteration
($n=1,\dots,N$), we run QA at $(c_{1},\cdots,c_{N})$ once and at
$(c_{1},\cdots,c_{i}+\Delta,\cdots,c_{N})$ for $i$-th site $N-(n-1)$ times as
we already fixed $n-1$ sites. Thus, the total number is
$\sum_{n=1}^{N}(1+N-n+1)=N(N+3)/2$. After we take into account this overhead,
the TTS of QGO becomes comparable with the other two as seen in Fig. 6(b) in
blank blue symbols. It is nevertheless concluded that QGO achieves much better
success probability with similar computation times.

In order to further understand how the algorithm leads to better results, we
have analyzed the time dependence of overlaps of the wave function with the
instantaneous ground states of three kinds of the Hamiltonian, (i) the final
Ising model $\mathcal{H}^{z}$, (ii) the transverse-field Ising model
$A(s)\mathcal{H}^{z}+B(s)\mathcal{H}^{x}$, and (iii) the full Hamiltonian
$A(s)\mathcal{H}^{z}+B(s)\mathcal{H}^{x}+C(s)\mathcal{H}^{y}$, for QGO and QA.
Results are plotted in Fig. 7 (a) for the ferromagnetic system and (b) for a
spin glass problem.

\begin{overpic}[height=128.0374pt]{fig10a.pdf} \put(0.0,73.0){(a)}
\end{overpic}\begin{overpic}[height=128.0374pt]{fig10b.pdf}
\put(0.0,73.0){(b)} \end{overpic}

Figure 7:  Time dependence of overlaps of the wave function with the
instantaneous ground states of three kinds of the Hamiltonian, (blue) the
final Ising model $\mathcal{H}^{z}$, (orange) the transverse-field Ising model
$A(s)\mathcal{H}^{z}+B(s)\mathcal{H}^{x}$, and (green) the full Hamiltonian
$A(s)\mathcal{H}^{z}+B(s)\mathcal{H}^{x}+C(s)\mathcal{H}^{y}$, for (solid) QGO
and (dashed) QA. The problems are (a) the ferromagnetic system, and (b) a spin
glass problem, both with eight spins. Curves in the range [0.97, 100] beyond
the black vertical line are magnified for better resolution.

It is observed that QGO does not directly let the system follow the
instantaneous ground state of the full Hamiltonian but it succeeds in
enhancing the overlap with the final Ising Hamiltonian at a relatively early
stage of computation.

###  3.4 Further simplifications of the algorithm

To see if we can further simplify the algorithm without much compromising the
performance, we test two possibilities in this section.

####  3.4.1 $y$-field optimization

The first example is to drive the state only by $\mathcal{H}^{y}$, i.e. by a
rotation around the $y$ axis,

$\displaystyle
U(\theta)=\exp\left(\frac{i}{2}\sum_{i}\theta_{i}(t)\sigma_{i}^{y}\right)\equiv\prod_{i}R_{y}(\theta_{i})$
(13)

starting from the ground state of $\mathcal{H}^{x}$ and to variationally
optimize the parameters $\theta$ to minimize $E$ as illustrated symbolically
in Fig. 8 (a).

\begin{overpic}[width=184.9429pt,clip,trim=80.29999pt 311.16249pt 70.2625pt
321.2pt]{fig7a.pdf} \put(-5.0,32.0){(a)}
\end{overpic}\begin{overpic}[height=128.0374pt]{fig7b.pdf} \put(0.0,75.0){(b)}
\end{overpic}

Figure 8:  (a) Quantum circuit of the $y$-field optimization. (b) Size
dependence of the success probability of QGO (blue) for $\tau=1$ (circle) and
$5$ (square) and the $y$-field optimization (purple diamond). The shaded
region represents 95% confidence intervals for each line.

Results for 100 spin glass instances are shown in Fig. 8 (b). It is seen that
this $y$-field optimization performs slightly worse than QGO for $\tau=1$. It
is also noted by comparison with Fig. 6 (a) that the $y$-field optimization
leads to better results than QA and SA with $\tau=1$. Thus the present method
may be useful for some purposes due to its simplicity and its easiness of
implementation possibly on a gate-based hardware.

####  3.4.2 Single-shot QGO

Let us next study what happens if we fix the signs of coefficients of $c_{i}$
in a single shot without iteration. We first fix one of the parameters, e.g.
$c_{1}=c_{\rm opt}^{N}$ (spin-up), to break the $\mathbb{Z}_{2}$ symmetry.
Signs of other sites are next fixed according to the gradient of the
optimization function, $1-P_{\rm gs}$ or $E$, near the origin as shown in
Algorithm 2.

Algorithm 2 Single-shot QGO

0: QA measure $f(b,{\bm{c}})$, $b_{\rm opt}^{N}$, $c_{\rm opt}^{N}$

0: a solution of the cost function

1: $b\leftarrow b_{\rm opt}^{N}$

2: ${\bm{c}}\leftarrow(c_{\rm opt}^{N},0,\cdots,0)$

3: ${\bm{g}}\leftarrow\nabla
f(b,{\bm{c}})=(\frac{f(b,c_{1}+\Delta,\cdots,c_{N})-f(b,{\bm{c}})}{\Delta},\cdots,\frac{f(b,c_{1},\cdots,c_{N}+\Delta)-f(b,{\bm{c}})}{\Delta})$

4: $c_{i}\leftarrow-c_{\rm opt}^{N}\operatorname{sgn}g_{i}$ for $i\geq 2$

5: return $\operatorname{sgn}{\bm{c}}$

Figure 9 (a) shows the $c_{i}$ dependence of $1-P_{\rm gs}$ for each spin to
confirm our assumption that the optimal $c_{i}$ is close to $\pm c_{\rm
opt}^{N}$. Blue-dashed and solid curves represent that the optimal values are
negative and positive, respectively. Gray-dotted curves indicate re-evaluation
of the fixed-parameter, similarly to Fig. 5. The correct solution has been
obtained by this single-shot QGO using $1-P_{\rm gs}$ for this instance. When
the energy measure $E$ is employed for optimization as shown in Fig. 9 (b),
two spins (the sixth and seventh) are fixed incorrectly. As illustrated in
Fig. 9 (c), the ground state is retrieved if the two spins are re-evaluated
after the other six spins are fixed. This suggests that we need iterations of
the process for better results if we use $E$ for optimization.

\begin{overpic}[width=327.20668pt]{fig8a.pdf} \put(3.0,13.0){(a)}
\end{overpic}\begin{overpic}[width=327.20668pt]{fig8b.pdf} \put(3.0,14.0){(b)}
\end{overpic}\begin{overpic}[width=327.20668pt]{fig8c.pdf} \put(3.0,14.0){(c)}
\end{overpic}

Figure 9:  $c_{i}$ dependence of the single-shot QGO with (a) $1-P_{\rm gs}$
and (b) $E$ for optimization in an eight-spin spin glass problem. (c) The
sixth and seventh parameters are re-evaluated by fixing other parameters to
$\pm c_{\rm opt}$ according to the results in (b). Blue-dashed and solid
curves represent that the optimal values are negative and positive,
respectively. Gray-dotted curves indicate re-evaluation of the already fixed
parameters similarly to Fig. 5. Vertical scales are normalized for each
estimation in (b) and (c).

It is also noticed that minimum values in the sixth spin in Fig. 9 (b) are not
around $\pm c_{\rm opt}^{N=8}=\pm 1.563$ but are around $0.7$ (red arrow).
This implies that our assumption, optimal values are around $\pm c_{\rm
opt}^{N}$ is not valid in this case, which would lead to the difference in the
performance between (a) and (b).

Figure 10 shows the performances of the single-shot QGO for two measures of
optimization for 100 spin glass instances. The single-shot QGO with the
fidelity measure always succeeds in identifying the ground state, while the
success probability quickly dumps under the energy measure. The fidelity
measure is the overlap between the final state and the ground state and thus
includes information of the ground state, leading to better results.

Figure 10:  Size dependence of the success probability of the single-shot QGO
with fidelity (blue circle) and energy (orange x) measures averaged out the
100 instances of spin glass problem. The shaded region represents 95%
confidence intervals for each line.

##  4 Summary and discussion

We have formulated a variational algorithm of quantum greedy optimization QGO
aiming to achieve improved performance of quantum annealing through a simple
process of adjustment of signs of $y$-field coefficients. The result showed
notable improvements in the success probability over the original method as
well as over classical simulated annealing under comparable computational
costs. Since we adjust the signs of coefficients sequentially, the energy
landscape of each step is very simple with only a single minimum at almost the
same absolute values of coefficients, the adjustment or optimization process
does not encounter the problems of barren plateau in deep variational quantum
curcuits [68] or highly complicated landscapes which plague well-known
variational algorithms such as VQE and QAOA [69]. Another advantage of the
present method is the simplicity of the additional term in the Hamiltonian,
the $y$-field, which can be rewritten in terms of the transverse-field Ising
model without $y$ field by a rotation in the spin space. The latter can be
implemented experimentally if the hardware can be designed to allow for the
non-monotonic time dependence of coefficients of the $x$ field and the Ising
part of the Hamiltonian. It is in principle possible to apply the same idea to
a final Hamiltonian with $x$ and $y$ components of the Pauli matrix, typically
for problems of quantum state preparation in chemistry, although it is non-
trivial if the present protocol would lead to satisfactory results in such
cases. In the real settings of quantum processors, we can not measure the
fidelity because the target state is unknown. The fidelity analysis provides
the algorithm’s upper limit, revealing that the success probability can be
100% and suggesting that the goodness of the cost function defines the
algorithm’s performance. This finding is informative for the further
development of algorithms.

The code and data are available in the supplementary materials.

##  Appendix A Time dependence of the coefficients in the rotated frame

An example of the behavior of the coefficients in the rotated frame, Eq. (6),
is given in this Appendix. Let us write the Hamiltonian as

$\mathcal{H}_{\rm
eff}=A(t)\mathcal{H}^{z}-\sum_{i}B_{i}^{\prime}(t)\sigma_{i}^{x}+\sum_{i}C_{i}^{\prime}(t)\sigma_{i}^{z},$
(A.1)

where

$\displaystyle B_{i}^{\prime}(t)$ $\displaystyle=$
$\displaystyle\sqrt{B(t)^{2}+C_{i}(t)^{2}}$ (A.2) $\displaystyle
C_{i}^{\prime}(t)$ $\displaystyle=$
$\displaystyle-\frac{bc_{i}[\pi(1-t/\tau)\sin(2\pi t/\tau)+\sin^{2}(\pi
t/\tau)]}{2\tau[b^{2}(1-t/\tau)^{2}+c_{i}^{2}\sin^{4}(\pi t/\tau)]}.$ (A.3)

Figure A.1 shows the time dependence of $A(t)$, $B_{i}^{\prime}(t)$ and
$C_{i}^{\prime}(t)$ for $a=1$, $b=0.5$, $c_{i}=1.5$ and $\tau=1$. It is
observed that $B^{\prime}(t)$ and $C^{\prime}(t)$ are non-monotonic, which
suppresses system’s transition to excited states in the final state.

Figure A.1:  Time development of the coefficients $A(t)$, $B_{i}^{\prime}(t)$
and $C_{i}^{\prime}(t)$ for $a=1$, $b=0.5$, $c_{i}=1.5$, and $\tau=1$.

##  Appendix B System behavior of the mean-field theory

To visualize the system behavior under the mean-field theory discussed in part
3.1 of section 3, we plot the trajectory of the system on the Bloch sphere in
Fig. B.1 for a number of parameter values with $b=0.5$ (green) and $c=1.5$
(long dashed) being closest to the optima.

Figure B.1:  Trajectories of the system under the mean-field theory for $b=$
(0.0, 0.5, 1.0, 1.5, 2.0) and $c=$ (0.1, 0.5, 1.0, 1.5, 2.0).

The system starts along the $x$ axis and is supposed to end ideally at the
north pole $|0\rangle$. The optimal parameter values are observed to
approximately follow a straightforward path connecting these points.

Figure B.2 compares the exact counterdiabatic driving function for the mean-
field theory [40] and our approximate result. The parameterized Ising
Hamiltonian is

$\mathcal{H}^{z}=-g\braket{\sigma^{z}}\sigma^{z}-h\sigma^{z}.$ (B.1)

Although QGO tested in the main text is $g=1$ and $h=0$, the parameter $h$
should be non-zero in the exact calculation to avoid the double degeneracy of
states. Exact solutions are obtained for $(g,h)=$ $(0,1)$, $(1,0.01)$,
$(1,0.1)$, $(1,0.3)$ and $(1,1)$. The system follows the instantaneous ground
state under these exact solutions whereas it does not under our approximate
scheme, to which $(g,h)=(1,0.3)$ shows the closest behavior.

\begin{overpic}[height=128.0374pt]{figS3a.pdf} \put(0.0,78.0){(a)}
\end{overpic}\begin{overpic}[height=128.0374pt]{figS3b.pdf}
\put(3.0,75.0){(b)} \end{overpic}

Figure B.2:  Exact solution of (a) counter-diabatic Y-field, and (b)
magnetization for QGO and $(g,h)=$ $(0,1)$, $(1,0.01)$, $(1,0.1)$, $(1,0.3)$
and $(1,1)$.

##  Appendix C Enhanced scales of performance measure

For better resolution of the vertical axes of Figs. 4 and 5, we show in Figs.
C.1 and C.2 the same data with enhanced scales of vertical axes.

\begin{overpic}[width=199.16928pt]{figS4.pdf} \put(1.0,82.0){(a)}
\put(53.0,82.0){(b)} \end{overpic} Figure C.1:  The same plot as in Fig. 4
with different vertical scales.  Figure C.2:  Same plot as in Fig. 5 with
different scales of vertical axis.

##  Appendix D $y$-field optimization

The simple $y$-field optimization algorithm can be executed the circuit model
in Fig. 8(a), and described as follows:

Algorithm 3 Simple Y-field optimization

0: energy measure $f({\bm{\theta}})$

0: a solution of the cost function

1: ${\bm{\theta}}\leftarrow(0,\cdots,0)$

2: repeat

3: ${\bm{g}}\leftarrow\nabla
f({\bm{\theta}})=(\frac{f(\theta_{1}+\Delta,\cdots,\theta_{N})-f({\bm{\theta}})}{\Delta},\cdots,\frac{f(\theta_{1},\cdots,\theta_{N}+\Delta)-f({\bm{\theta}})}{\Delta})$

4: $i\leftarrow\operatorname*{arg\,max}_{j\in\\{j|\theta_{j}=0\\}}|g_{j}|$

5: $\theta_{i}\leftarrow-\frac{\pi}{2}\operatorname{sgn}g_{i}$

6: until $\theta_{i}\neq 0$ for all $i$

7: return $\operatorname{sgn}{\bm{\theta}}$

## References

  * Kadowaki and Nishimori [1998] Tadashi Kadowaki and Hidetoshi Nishimori.  Quantum annealing in the transverse Ising model.  _Phys. Rev. E_ , 58(5):5355–5363, 1998.  10.1103/PhysRevE.58.5355.
  * Kadowaki [1998] Tadashi Kadowaki.  _Study of Optimization Problems by Quantum Annealing_.  PhD thesis, Tokyo Institute of Technology, 1998.
  * Brooke et al. [1999] J. Brooke, D Bitko, T Rosenbaum, and G Aeppli.  Quantum Annealing of a Disordered Magnet.  _Science_ , 284(5415):779–781, 1999.  10.1126/science.284.5415.779.
  * Santoro [2002] Giuseppe E. Santoro.  Theory of Quantum Annealing of an Ising Spin Glass.  _Science_ , 295(5564):2427–2430, 2002.  10.1126/science.1068774.
  * Santoro and Tosatti [2006] Giuseppe E Santoro and Erio Tosatti.  Optimization using quantum mechanics: quantum annealing through adiabatic evolution.  _J. Phys. A_ , 39(36):R393–R431, 2006.  10.1088/0305-4470/39/36/R01.
  * Das and Chakrabarti [2008] Arnab Das and Bikas K. Chakrabarti.  Colloquium : Quantum annealing and analog quantum computation.  _Rev. Mod. Phys._ , 80(3):1061–1081, 2008.  10.1103/RevModPhys.80.1061.
  * Morita and Nishimori [2008] Satoshi Morita and Hidetoshi Nishimori.  Mathematical foundation of quantum annealing.  _J. Math. Phys_ , 49(12):1–47, 2008.  10.1063/1.2995837.
  * Tanaka et al. [2017] Shu Tanaka, Ryo Tamura, and B. K. Chakrabarti.  _Quantum spin glasses, annealing and computation_.  Cambridge University Press, 2017.  ISBN 1107113199.
  * Farhi et al. [2001] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, Joshua Lapan, Andrew Lundgren, and Daniel Preda.  A Quantum Adiabatic Evolution Algorithm Applied to Random Instances of an NP-Complete Problem.  _Science_ , 292(5516):472–475, 2001.  10.1126/science.1057726.
  * Farhi et al. [2014] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A Quantum Approximate Optimization Algorithm.  2014\.  arXiv:1411.4028.
  * Jörg et al. [2008] Thomas Jörg, Florent Krzakala, Jorge Kurchan, and Anthony C Maggs.  Simple Glass Models and Their Quantum Annealing.  _Phys. Rev. Lett._ , 101(14):147204, 2008.  10.1103/PhysRevLett.101.147204.
  * Young et al. [2010] A. P. Young, S. Knysh, and V. N. Smelyanskiy.  First-order phase transition in the quantum adiabatic algorithm.  _Phys. Rev. Lett._ , 104(2):020502, 2010.  10.1103/PhysRevLett.104.020502.
  * Jörg et al. [2010] Thomas Jörg, Florent Krzakala, Guilhem Semerjian, and Francesco Zamponi.  First-order transitions and the performance of quantum algorithms in random optimization problems.  _Phys. Rev. Lett._ , 104(20):207206, 2010.  10.1103/PhysRevLett.104.207206.
  * Seki and Nishimori [2012] Yuya Seki and Hidetoshi Nishimori.  Quantum annealing with antiferromagnetic fluctuations.  _Phys. Rev. E_ , 85(5):051112, 2012.  10.1103/PhysRevE.85.051112.
  * Seki and Nishimori [2015] Yuya Seki and Hidetoshi Nishimori.  Quantum annealing with antiferromagnetic transverse interactions for the Hopfield model.  _J. Phys. A_ , 48(33):335301, 2015.  10.1088/1751-8113/48/33/335301.
  * Seoane and Nishimori [2012] Beatriz Seoane and Hidetoshi Nishimori.  Many-body transverse interactions in the quantum annealing of the p-spin ferromagnet.  _J. Phys. A_ , 45(43):435301, 2012.  10.1088/1751-8113/45/43/435301.
  * Nishimori and Takada [2017] Hidetoshi Nishimori and Kabuki Takada.  Exponential enhancement of the efficiency of quantum annealing by non-stoquastic Hamiltonians.  _Frontiers in ICT_ , 4(FEB):2, 2017.  10.3389/fict.2017.00002.
  * Susa et al. [2018a] Yuki Susa, Yu Yamashiro, Masayuki Yamamoto, and Hidetoshi Nishimori.  Exponential Speedup of Quantum Annealing by Inhomogeneous Driving of the Transverse Field.  _J. Phys. Soc. Jpn._ , 87(2):023002, 2018a.  10.7566/JPSJ.87.023002.
  * Susa et al. [2018b] Yuki Susa, Yu Yamashiro, Masayuki Yamamoto, Itay Hen, Daniel A Lidar, and Hidetoshi Nishimori.  Quantum annealing of the p-spin model under inhomogeneous transverse field driving.  _Phys. Rev. A_ , 98(4):042326, 2018b.  10.1103/PhysRevA.98.042326.
  * Mohseni et al. [2018] Masoud Mohseni, Johan Strumpfer, and Marek M Rams.  Engineering non-equilibrium quantum phase transitions via causally gapped Hamiltonians.  _New J. Phys._ , 20(10):105002, 2018.  10.1088/1367-2630/aae3ed.
  * Adame and McMahon [2020] Juan Ignacio Adame and Peter McMahon.  Inhomogeneous driving in quantum annealers can result in orders-of-magnitude improvements in performance.  _Quantum Sci. Tech._ , 5(3):035011, 2020.  10.1088/2058-9565/ab935a.
  * Rams et al. [2016] Marek M. Rams, Masoud Mohseni, and Adolfo del Campo.  Inhomogeneous quasi-adiabatic driving of quantum critical dynamics in weakly disordered spin chains.  _New J. Phys._ , 18(12):123034, 2016.  10.1088/1367-2630/aa5079.
  * Gómez-Ruiz and Del Campo [2019] F. J. Gómez-Ruiz and A. Del Campo.  Universal Dynamics of Inhomogeneous Quantum Phase Transitions: Suppressing Defect Formation.  _Phys. Rev. Lett._ , 122(8):080604, 2019.  10.1103/PhysRevLett.122.080604.
  * Ohkuwa et al. [2018] Masaki Ohkuwa, Hidetoshi Nishimori, and Daniel A Lidar.  Reverse annealing for the fully connected p-spin model.  _Phys. Rev. A_ , 98(2):022314, 2018.  10.1103/PhysRevA.98.022314.
  * Yamashiro et al. [2019] Yu Yamashiro, Masaki Ohkuwa, Hidetoshi Nishimori, and Daniel A Lidar.  Dynamics of reverse annealing for the fully-connected $p$-spin model.  _Phys. Rev. A_ , 100:052321, 2019.  10.1103/PhysRevA.100.052321.
  * Passarelli et al. [2020a] Gianluca Passarelli, Ka-Wa Yip, Daniel A. Lidar, Hidetoshi Nishimori, and Procolo Lucignano.  Reverse quantum annealing of the p -spin model with relaxation.  _Phys. Rev. A_ , 101:022331, 2020a.  doi.org/10.1103/PhysRevA.101.022331.
  * Venturelli and Kondratyev [2019] Davide Venturelli and Alexei Kondratyev.  Reverse quantum annealing approach to portfolio optimization problems.  _Quantum Machine Intelligence_ , 1:17–30, 2019.  10.1007/s42484-019-00001-w.
  * Chancellor [2020] Nicholas Chancellor.  Fluctuation-guided search in quantum annealing.  _Phys. Rev. A_ , 102(6):062606, 2020.  10.1103/PhysRevA.102.062606.
  * Marshall et al. [2019] Jeffrey Marshall, Davide Venturelli, Itay Hen, and Eleanor G. Rieffel.  Power of Pausing: Advancing Understanding of Thermalization in Experimental Quantum Annealers.  _Phys. Rev. Appl._ , 11(4):044083, 2019.  10.1103/PhysRevApplied.11.044083.
  * Ikeda et al. [2019] Kazuki Ikeda, Yuma Nakamura, and Travis S Humble.  Application of Quantum Annealing to Nurse Scheduling Problem.  _Sci. Rep._ , 9(1):12837, 2019.  10.1038/s41598-019-49172-3.
  * Grant et al. [2020] Erica Grant, Travis Humble, and Benjamin Stump.  Benchmarking Quantum Annealing Controls with Portfolio Optimization.  _Phys. Rev. Appl._ , 15(1):014012, 2020.  10.1103/PhysRevApplied.15.014012.
  * Rocutto et al. [2021] Lorenzo Rocutto, Claudio Destri, and Enrico Prati.  Quantum Semantic Learning by Reverse Annealing of an Adiabatic Quantum Computer.  _Adv. Quant. Tech._ , 4(2):2000133, 2021.  10.1002/qute.202000133.
  * Arai et al. [2021] Shunta Arai, Masayuki Ohzeki, and Kazuyuki Tanaka.  Mean field analysis of reverse annealing for code-division multiple-access multiuser detection.  _Phys. Rev. Res._ , 3(3):033006, 2021.  10.1103/PhysRevResearch.3.033006.
  * Golden and O’Malley [2021] John K. Golden and Daniel O’Malley.  Reverse Annealing for Nonnegative/Binary Matrix Factorization.  _PLoS One_ , 16:e0244026, 2021.  10.1371/journal.pone.0244026.
  * Demirplak and Rice [2003] Mustafa Demirplak and Stuart A. Rice.  Adiabatic population transfer with control fields.  _J. Phys. Chem. A_ , 107(46):9937–9945, 2003\.  10.1021/jp030708a.
  * Demirplak and Rice [2005] Mustafa Demirplak and Stuart A. Rice.  Assisted adiabatic passage revisited.  _J. Phys. Chem. B_ , 109(14):6838–6844, 2005\.  10.1021/jp040647w.
  * Berry [2009] M. V. Berry.  Transitionless quantum driving.  _J. Phys. A_ , 42(36):365303, 2009.  10.1088/1751-8113/42/36/365303.
  * Chen et al. [2010] Xi Chen, A. Ruschhaupt, S. Schmidt, A. Del Campo, D. Guéry-Odelin, and J. G. Muga.  Fast optimal frictionless atom cooling in harmonic traps: Shortcut to adiabaticity.  _Phys. Rev. Lett._ , 104(6):063002, 2010.  10.1103/PhysRevLett.104.063002.
  * Chen et al. [2011] Xi Chen, E. Torrontegui, and J. G. Muga.  Lewis-Riesenfeld invariants and transitionless quantum driving.  _Phys. Rev. A_ , 83(6):062116, 2011.  10.1103/PhysRevA.83.062116.
  * Takahashi [2013] Kazutaka Takahashi.  Transitionless quantum driving for spin systems.  _Phys. Rev. E_ , 87(6):062117, 2013.  10.1103/PhysRevE.87.062117.
  * Jarzynski [2013] Christopher Jarzynski.  Generating shortcuts to adiabaticity in quantum and classical dynamics.  _Phys. Rev. A_ , 88(4):040101, 2013.  10.1103/PhysRevA.88.040101.
  * Torrontegui et al. [2012] E. Torrontegui, S. Ibáñez, S. Martínez-Garaot, M. Modugno, A. del Campo, D. Guéry-Odelin, A. Ruschhaupt, Xi Chen, and J. G. Muga.  Shortcuts to adiabaticity.  _Adv. Atomic, Mol., Opt. Phys._ , 62:117–169, 2012.  10.1016/B978-0-12-408090-4.00002-5.
  * del Campo [2013] Adolfo del Campo.  Shortcuts to adiabaticity by counter-diabatic driving.  _Phys. Rev. Lett._ , 111(10):100502, 2013.  10.1103/PhysRevLett.111.100502.
  * Takahashi [2017] Kazutaka Takahashi.  Shortcuts to adiabaticity for quantum annealing.  _Phys. Rev. A_ , 95(1):012309, 2017.  10.1103/PhysRevA.95.012309.
  * Sels and Polkovnikov [2017] Dries Sels and Anatoli Polkovnikov.  Minimizing irreversible losses in quantum systems by local counterdiabatic driving.  _Proc. Nat. Acad. Sci._ , 114(20):E3909–E3916, 2017.  10.1073/pnas.1619826114.
  * Özgüler et al. [2018] A. Barış Özgüler, Robert Joynt, and Maxim G. Vavilov.  Steering random spin systems to speed up the quantum adiabatic algorithm.  _Phys. Rev. A_ , 98(6):062311, 2018.  10.1103/PhysRevA.98.062311.
  * Claeys et al. [2019] Pieter W. Claeys, Mohit Pandey, Dries Sels, and Anatoli Polkovnikov.  Floquet-Engineering Counterdiabatic Protocols in Quantum Many-Body Systems.  _Phys. Rev. Lett._ , 123(9):090602, 2019.  10.1103/PhysRevLett.123.090602.
  * Hartmann and Lechner [2019] Andreas Hartmann and Wolfgang Lechner.  Rapid counter-diabatic sweeps in lattice gauge adiabatic quantum computing.  _New J. Phys._ , 21(4):043025, 2019.  10.1088/1367-2630/ab14a0.
  * Hatomura and Takahashi [2021] Takuya Hatomura and Kazutaka Takahashi.  Controlling and exploring quantum systems by algebraic expression of adiabatic gauge potential.  _Phys. Rev. A_ , 103(1):012220, 2021.  10.1103/PhysRevA.103.012220.
  * Prielinger et al. [2021] Luise Prielinger, Andreas Hartmann, Yu Yamashiro, Kohji Nishimura, Wolfgang Lechner, and Hidetoshi Nishimori.  Two-parameter counter-diabatic driving in quantum annealing.  _Phys. Rev. Res._ , 3(1):13227, 2021.  10.1103/physrevresearch.3.013227.
  * Imoto et al. [2021] Takashi Imoto, Yuya Seki, Yuichiro Matsuzaki And, and Shiro Kawabata.  Quantum annealing with twisted fields.  2021\.  arXiv:2111.15283.
  * Kirkpatrick et al. [1983] S Kirkpatrick, ; C D Gelatt, and ; M P Vecchi.  Optimization by Simulated Annealing.  _Science_ , 220:671–680, 1983.  10.1126/science.220.4598.671.
  * Bravyi and Terhal [2009] Sergey Bravyi and Barbara Terhal.  Complexity of stoquastic frustration-free Hamiltonians.  _SIAM J. Comput._ , 39(4):1462, 2009.  10.1137/090750688.
  * Passarelli et al. [2020b] Gianluca Passarelli, Vittorio Cataudella, Rosario Fazio, and Procolo Lucignano.  Counterdiabatic driving in the quantum annealing of the p-spin model: a variational approach.  _Phys. Rev. Res._ , 2:013283, 2020b.  10.1103/PhysRevResearch.2.013283.
  * Broyden [1970] C. G. Broyden.  The convergence of a class of double-rank minimization algorithms 1. General considerations.  _IMA J. Appl. Math._ , 6(1):76–90, 1970.  10.1093/imamat/6.1.76.
  * R. Fletcher [1970] R. Fletcher.  New approach to variable metric algorithms.  _Comp. J._ , 13(3):317–322, 1970.  10.1093/comjnl/13.3.317.
  * Goldfarb [1970] Donald Goldfarb.  A family of variable-metric methods derived by variational means.  _Mathematics of Computation_ , 24(109):23–23, 1970.  10.1090/s0025-5718-1970-0258249-6.
  * Shanno [1970] D. F. Shanno.  Conditioning of quasi-Newton methods for function minimization.  _Math. Comp._ , 24(111):647–647, 1970.  10.1090/s0025-5718-1970-0274029-x.
  * Nishimori and Ortiz [2011] H. Nishimori and G. Ortiz.  _Elements of Phase Transitions and Critical Phenomena_.  Oxford, United Kingdom: Oxford University Press, 2011.
  * Kolmogorov and Zabih [2004] Vladimir Kolmogorov and Ramin Zabih.  What energy functions can be minimized via graph cuts?  _IEEE Trans. Patt. Anal. Mach. Intell._ , 26(2):147–159, 2004.  10.1109/TPAMI.2004.1262177.
  * Boros and Hammer [2002] Endre Boros and Peter L. Hammer.  Pseudo-Boolean optimization.  _Discrete Appl. Math._ , 123(1-3):155–225, 2002\.  10.1016/S0166-218X(01)00341-9.
  * Boros et al. [2006] Endre Boros, Peter L Hammer, and Gabriel Tavares.  Preprocessing of Unconstrained Quadratic Binary Optimization.  _Rutgers Research Report_ , 188(1):99–104, 2006\.
  * Rother et al. [2007] Carsten Rother, Vladimir Kolmogorov, Victor Lempitsky, and Martin Szummer.  Optimizing binary MRFs via extended roof duality.  2007\.  10.1109/CVPR.2007.383203.
  * Karimi et al. [2017] Hamed Karimi, Gili Rosenberg, and Helmut G. Katzgraber.  Effective optimization using sample persistence: A case study on quantum annealers and various Monte Carlo optimization methods.  _Phys. Rev. E_ , 96(4):043312, 2017.  10.1103/PhysRevE.96.043312.
  * Sherrington and Kirkpatrick [1975] David Sherrington and Scott Kirkpatrick.  Solvable Model of a Spin-Glass.  _Phys. Rev. Lett._ , 35(26):1792–1796, 1975.  10.1103/PhysRevLett.35.1792.
  * Johansson et al. [2013] J. R. Johansson, P. D. Nation, and Franco Nori.  QuTiP 2: A Python framework for the dynamics of open quantum systems.  _Comp. Phys. Commun._ , 184(4):1234–1240, 2013\.  10.1016/J.CPC.2012.11.019.
  * Albash and Lidar [2018] Tameem Albash and Daniel A. Lidar.  Adiabatic quantum computation.  _Rev. Mod. Phys._ , 90(1):015002, 2018.  10.1103/RevModPhys.90.015002.
  * McClean et al. [2018] Jarrod R. McClean, Sergio Boixo, Vadim N. Smelyanskiy, Ryan Babbush, and Hartmut Neven.  Barren plateaus in quantum neural network training landscapes.  _Nature Communications 2018 9:1_ , 9(1):1–6, 2018\.  10.1038/s41467-018-07090-4.
  * Harrigan et al. [2020] Matthew P. Harrigan, Kevin J. Sung, Matthew Neeley, Kevin J. Satzinger, Frank Arute, Kunal Arya, Juan Atalaya, Joseph C. Bardin, Rami Barends, Sergio Boixo, Michael Broughton, Bob B. Buckley, David A. Buell, Brian Burkett, Nicholas Bushnell, Yu Chen, Zijun Chen, Ben Chiaro, Roberto Collins, William Courtney, Sean Demura, Andrew Dunsworth, Daniel Eppens, Austin Fowler, Brooks Foxen, Craig Gidney, Marissa Giustina, Rob Graff, Steve Habegger, Alan Ho, Sabrina Hong, Trent Huang, L. B. Ioffe, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Cody Jones, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Seon Kim, Paul V. Klimov, Alexander N. Korotkov, Fedor Kostritsa, David Landhuis, Pavel Laptev, Mike Lindmark, Martin Leib, Orion Martin, John M. Martinis, Jarrod R. McClean, Matt McEwen, Anthony Megrant, Xiao Mi, Masoud Mohseni, Wojciech Mruczkiewicz, Josh Mutus, Ofer Naaman, Charles Neill, Florian Neukart, Murphy Yuezhen Niu, Thomas E. O’Brien, Bryan O’Gorman, Eric Ostby, Andre Petukhov, Harald Putterman, Chris Quintana, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Andrea Skolik, Vadim Smelyanskiy, Doug Strain, Michael Streif, Marco Szalay, Amit Vainsencher, Theodore White, Z. Jamie Yao, Ping Yeh, Adam Zalcman, Leo Zhou, Hartmut Neven, Dave Bacon, Erik Lucero, Edward Farhi, and Ryan Babbush.  Quantum Approximate Optimization of Non-Planar Graph Problems on a Planar Superconducting Processor.  _Nature Physics_ , 17(3):332–336, 2020.  10.1038/s41567-020-01105-y.