# Effectiveness of quantum annealing for continuous-variable optimization

Shunta Arai1<EMAIL_ADDRESS>Hiroki Oshiyama2 Hidetoshi
Nishimori1,2,3 1International Research Frontiers Initiative, Tokyo Institute
of Technology, Shibaura, Minato-ku, Tokyo 108-0023, Japan
2 Graduate School of Information Sciences, Tohoku University, Sendai 980-8579,
Japan
3 RIKEN, Interdisciplinary Theoretical and Mathematical Sciences (iTHEMS),
Wako, Saitama 351-0198, Japan

###### Abstract

The application of quantum annealing to the optimization of continuous-
variable functions is a relatively unexplored area of research. We test the
performance of quantum annealing applied to a one-dimensional continuous-
variable function with a rugged energy landscape. After domain-wall encoding
to map a continuous variable to discrete Ising variables, we first benchmark
the performance of the real hardware, the D-Wave 2000Q, against several state-
of-the-art classical optimization algorithms designed for continuous-variable
problems to find that the D-Wave 2000Q matches the classical algorithms in a
limited domain of computation time. Beyond this domain, classical global
optimization algorithms outperform the quantum device. Next, we examine
several optimization algorithms that are applicable to the Ising formulation
of the problem, such as the TEBD (time-evolving block decimation) to simulate
ideal coherent quantum annealing, simulated annealing, simulated quantum
annealing, and spin-vector Monte Carlo. The data show that TEBD’s coherent
quantum annealing achieves far better results than the other approaches,
demonstrating the effectiveness of coherent tunneling. From these two types of
benchmarks, we conclude that the hardware realization of quantum annealing has
the potential to significantly outperform the best classical algorithms if
thermal noise and other imperfections are sufficiently suppressed and the
device operates coherently, as demonstrated in recent short-time quantum
simulations.

##  I Introduction

Quantum annealing (QA) is a metaheuristic designed to solve combinatorial
optimization problems by controlling quantum fluctuations [1, 2, 3, 4, 5, 6,
7, 8, 9, 10, 11]. The hardware implementation of QA is provided by D-Wave
Quantum Inc. [12, 13, 14, 15, 16] with more than 5000 qubits in the latest
version [17].

Combinatorial optimization problems with discrete variables are the main
target of QA, but the basic idea of using quantum fluctuations to search for
the lowest-energy state in the rugged energy landscape is also applicable to
problems with continuous variables. An example is the early paper by Finnilla
et al. [18] although they used an imaginary-time version of the Schrödinger
equation without quantum coherence. Stella et al. [19] studied a few one-
dimensional problems with and without rugged energy landscapes and compared QA
with its classical counterpart of simulated annealing (SA) by directly solving
the Schrödinger equation for QA and the Fokker-Planck equation for SA.

They found that the relative efficiency of QA and SA depends on the shape of
the energy landscape. Simulated annealing is influenced by the height of the
energy barrier, whereas QA reflects the structure of the energy spectrum,
which affects the probability of quantum tunneling. An apparent difference in
performance between SA and QA was clearly seen in a parabolic washboard-like
potential with a single global minimum and many local minima as proposed by
Shinomoto and Kabashima [20]. These latter authors used a coarse-grained
phenomenological version of the problem to show that the best possible
performance of SA can be achieved when the rate of temperature decrease is
proportional to the inverse of logarithm, which results in a very slow
logarithmic decrease of the residual energy as a function of the annealing
time. Stella et al. discussed also QA phenomenologically for the same
washboard-like potential by considering a coarse-grained tight-binding
Hamiltonian under a semi-classical approximation, which consists of a harmonic
potential and a Laplacian term with a coefficient representing a semi-
classical WKB tunneling effect. They found that the residual energy shows a
power law decrease, which is faster than the logarithmic decrease in SA as
found by Shinomoto and Kabashima. Inack and Pilati studied the performance of
QA by the projective quantum Monte Carlo method [21] and found that this
imaginary-time formulation behaves better than SA in the long-time regime. Koh
and Nishimori recently analyzed a problem with a Shinomoto-Kabashima-like
potential with many local minima [22], which is known as the Rastrigin
function in the field of nonconvex continuous optimization problems [23]. They
solved the time-dependent Schrödinger equation directly numerically for QA and
the master equation for SA and confirmed that the residual energy under QA
decreases in a power law with an arbitrarily large power achievable by
adjusting the annealing schedule, whereas SA has a smaller power. This is a
simple but nontrivial example of a continuous optimization problem, in which
QA has an advantage over SA.

The studies mentioned above were based on theoretical or numerical approaches.
Several experimental investigations on real hardware have recently been
reported using domain wall encoding [24] to represent continuous variables in
terms of Ising spins. Examples include Abel et al. [25, 26] who observed
quantum tunneling behavior from a local minimum to the global minimum in a
double-well potential. In another paper, Abel et al. [27] have applied the
same method to two-dimensional continuous optimization problems with many
local minima using reverse annealing (RA) [28, 29, 30, 31], which is a variant
of QA starting from a classical spin configuration 111In Ref. [27], RA is
referred to as QA.. They graphically compared the output of QA on the D-Wave
hardware with data from classical continuous-variable optimization algorithms
and found that the former reached the global minimum much more reliably than
the latter. This is a very interesting development, but detailed analysis is
still needed to quantitatively demonstrate a clear advantage of QA over state-
of-the-art classical algorithms tailored for continuous-variable optimization.

In the present paper, we systematically investigate whether QA has an
advantage over classical algorithms for optimization of the Rastrigin function
in one dimension, adopting domain-wall encoding to map a continuous variable
to Ising spins. The first part of our test compares the hardware D-Wave 2000Q
with several well-established classical algorithms for continuous-variable
optimization, Nelder-Mead (NM) [33], conjugate gradient descent (CGD) [34],
basin-hopping (BH) [35], and differential evolution (DE) [36]. We have used
the D-Wave 2000Q, not the latest D-Wave Advantage, because the former produced
better results than the latter as will be described later, in particular in
Appendix B.1. The second part also runs the hardware, but we additionally test
a classical algorithm to faithfully simulate quantum annealing, TEBD (time-
evolving block decimation) [37, 38] based on a matrix-product representation
of the wave function, in order to clarify to what extent the hardware runs
coherently. We also study several classical algorithms for combinatorial
optimization, simulated annealing (SA) [39], simulated quantum annealing (SQA)
[40], and spin-vector Monte Carlo (SVMC) [41], to investigate how classical
algorithms are affected by the barrier height compared to the quantum case.

This paper is organized as follows. In Sec. II, we explain the Rastrigin
function and how to implement continuous-variable optimization with domain
wall encoding. In Sec. III, we perform benchmark tests for the Rastrigin
function on the D-Wave 2000Q and several classical algorithms developed for
continuous-variable optimization. In Sec. IV, we next compare D-Wave 2000Q,
TEBD, SA, SQA, and SVMC, all of which are for discrete variable problems, to
see how thermal effects influence the performance. Finally, Sec. V discusses
the results and concludes the paper. A few technical details are described in
the Appendixes.

##  II The Problem and Methods

In this section, we explain the mathematical definition of the Rastrigin
function and how to optimize the continuous function with the domain wall
encoding method on the D-Wave 2000Q.

Figure 1: Illustration of the Rastrigin function with $k=0.5$ and $w_{0}=0.2$
for different values of $h_{0}$.

###  II.1  Rastrigin function

The Rastrigin function is one of the standard benchmark functions for
nonconvex optimization problems [23]. We consider the one-dimensional
Rastrigin function which consists of a harmonic potential term with a unique
global minimum at $x=0$ and a cosine term that gives many local minima as

$\displaystyle
V_{\mathrm{RF}}(x)=\frac{1}{2}kx^{2}+\frac{h_{0}}{2}\left[1-\cos\left(\frac{2\pi
x}{w_{0}}\right)\right],$ (1)

where $k$ is the spring constant, and $h_{0}$ and $w_{0}$ are the control
parameters of the height and width of each local minimum. The above
mathematical formulation of the function was introduced by Koh and Nishimori
[22] and they called it the Shinomoto-Kabashima-like potential. The original
formulation of the Rastrigin function has no control parameters and
corresponds to $k=2$, $w_{0}=1$, and $h_{0}=20$ in Eq. (1). In this study, we
do not adopt the original formulation, but Koh and Nishimori’s formulation to
examine the dependence of algorithm performance on the shape of the energy
landscape. We illustrate the one-dimensional Rastrigin function with $k=0.5$
and $w_{0}=0.2$ for different values of $h_{0}$ in Fig. 1. Due to the second
cosine term, the simple gradient descent algorithm stops at a local minimum
and fails to find the global minimum efficiently.

###  II.2 Domain wall encoding for continuous function optimization

We briefly formulate QA and then describe the way to perform continuous
function optimization by quantum annealing with domain wall encoding.

The generic Hamiltonian to be implemented by QA is

$\displaystyle H(s)$
$\displaystyle=\frac{B(s)}{2}H_{0}-\frac{A(s)}{2}H_{\mathrm{TF}},$ (2)
$\displaystyle H_{0}$ $\displaystyle=\sum_{(i,j)\in
E}J_{ij}\hat{\sigma}_{i}^{z}\hat{\sigma}_{j}^{z}+\sum_{i=1}^{N}h_{i}\hat{\sigma}_{i}^{z},$
(3) $\displaystyle H_{\mathrm{TF}}$
$\displaystyle=\sum_{i=1}^{N}\hat{\sigma}_{i}^{x},$ (4)

where $\hat{\sigma}_{i}^{z}$ and $\hat{\sigma}_{i}^{x}$ are the Pauli
operators acting at site $i$, $J_{ij}$ is the coupling constant on the edge
$(i,j)$ in the graph $E$ describing the connectivity between qubits, $h_{i}$
is the local field at site $i$, and $N$ represents the number of spins
(qubits). The coefficients of the transverse field and the problem are denoted
as $A(s)$ and $B(s)$, respectively, controlled by the time-dependent annealing
parameter $s=t/t_{a}\in[0,1]$. The initial state is the ground state of Eq.
(4) with the coefficient $-A(0)/2$, and the system evolves from $t=0$ to
$t=t_{a}$ following the Schrödinger equation with the coefficients evolving as
$A(0)>0\rightarrow A(1)=0$ and $B(0)\approx 0\rightarrow B(1)>0$. The precise
forms of the scheduling functions of $A(s)$ and $B(s)$ depend on the physical
devices [42]. If the system evolves adiabatically, the final state of Eq. (2)
approaches the ground state of Eq. (3).

Following Refs. [24, 25, 26], we use the domain-wall encoding method to
rewrite a continuous variable optimization problem to a corresponding problem
with Ising spins. In this method, we map the value of a continuous variable to
the position of a domain wall within a spin chain. Concretely, we discretize a
continuous variable $x\in[x_{\mathrm{min}},x_{\mathrm{max}}]$ into $N-1$
values as

$\displaystyle x_{j}$ $\displaystyle~{}=x_{\mathrm{min}}+(j-1)\Delta
x~{}~{}(j=1,\dots,N-1)$ (5)

with the step size

$\displaystyle\Delta x$
$\displaystyle~{}=\frac{x_{\mathrm{max}}-x_{\mathrm{min}}}{N-1},$ (6)

and map each value $x_{j}$ to an $N$-spin state $\ket{x_{j}}$ which has a
domain wall at the $j$th bond: All spins at sites 1 to $j$ are down, and the
spins from $j+1$ to $N$ are up. Using the fact that
$\bra{x_{j}}(\hat{\sigma}_{i+1}-\hat{\sigma}_{i})/2\ket{x_{j}}=\delta_{ij}$
under this spin configuration, where $\delta_{ij}$ is Kronecker delta, we can
construct an $N$-spin Hamiltonian $H_{\mathrm{RF}}$ such that its diagonal
element in the basis $\ket{x_{j}}$ corresponds to $V_{\mathrm{RF}}(x_{j})$ as

$\displaystyle
H_{\mathrm{RF}}=\frac{1}{2}\sum_{i=1}^{N-1}V_{\mathrm{RF}}(x_{i})(\hat{\sigma}_{i+1}^{z}-\hat{\sigma}_{i}^{z}).$
(7)

To enforce that the ground state is a single domain wall state, we use a
penalty Hamiltonian that punishes multi-kink and no-kink states as

$\displaystyle
H_{\mathrm{DW}}=-J\sum_{i=1}^{N-1}\hat{\sigma}_{i}^{z}\hat{\sigma}_{i+1}^{z}+h\left(\hat{\sigma}_{1}^{z}-\hat{\sigma}_{N}^{z}\right),$
(8)

where the first term with $J>0$ is the ferromagnetic coupling to align
neighboring spins and the second term with $h>J$ is the magnetic field at both
ends of the chain to fix the spin at the first site in the down state and the
spin at the $N$th site in the up state. The degenerate ground states of
$H_{\mathrm{DW}}$ are $\ket{x_{j}}$ for $j=1,\dots,~{}N-1$ and the penalty
energy of each additional kink is $2J$.

Combining these two Hamiltonians, the domain wall encoding Hamiltonian for
optimizing $V_{\mathrm{RF}}(x)$ is given as

$\displaystyle H_{0}$ $\displaystyle=\lambda H_{\mathrm{RF}}+H_{\mathrm{DW}},$
(9)

where $\lambda$ is a parameter that controls the relative strength of the
penalty term $H_{\mathrm{DW}}$. For sufficiently small $\Delta x$ and
$\lambda$, the ground state of $H_{0}$ corresponds to the global minimum of
$V_{\mathrm{RF}}(x)$. If $\lambda$ is too large, the ground state of $H_{0}$
is no longer a single domain wall state.

Assuming that the potential function is differentiable and $\Delta x\ll 1$,
Eq. (7) can be reduced to

$\displaystyle H_{\mathrm{RF}}$ $\displaystyle=-\frac{1}{2}\Bigg{(}\Delta
x\sum_{i=2}^{N-1}\frac{V_{\mathrm{RF}}(x_{i})-V_{\mathrm{RF}}(x_{i-1})}{\Delta
x}\hat{\sigma}_{i}^{z}$
$\displaystyle\quad+V_{\mathrm{RF}}(x_{1})\hat{\sigma}_{1}^{z}-V_{\mathrm{RF}}(x_{N-1})\hat{\sigma}_{N}^{z}\Bigg{)}$
$\displaystyle\approx-\frac{1}{2}\Bigg{(}\Delta
x\sum_{i=2}^{N-1}\left.\frac{dV_{\mathrm{RF}}(x)}{dx}\right|_{x=x_{i}}\hat{\sigma}_{i}^{z}$
$\displaystyle+V_{\mathrm{RF}}(x_{1})\hat{\sigma}_{1}^{z}-V_{\mathrm{RF}}(x_{N-1})\hat{\sigma}_{N}^{z}\Bigg{)}$
$\displaystyle=\sum_{i=1}^{N}h_{i}\hat{\sigma}_{i}^{z},$ (10) $\displaystyle
h_{i}$
$\displaystyle=\begin{cases}-\frac{V_{\mathrm{RF}}(x_{1})}{2}\qquad(i=1)\\\
-\frac{\Delta
x}{2}\left.\frac{dV_{\mathrm{RF}}(x)}{dx}\right|_{x=x_{i}}\,(i=2,...,N-1)\\\
\frac{V_{\mathrm{RF}}(x_{N-1})}{2}\qquad(i=N)\end{cases}.$ (11)

The above construction of the Ising Hamiltonian for continuous function
optimization can be easily generalized to higher-dimensional problems [26].

##  III Comparison of QA with classical algorithms for continuous optimization

In this section, we show the results of benchmarks for QA and several
classical algorithms designed for the optimization of continuous variable
functions to investigate how QA applied to continuous variable optimization
compares with dedicated classical algorithms. In the present section, QA means
the hardware implementation, the D-Wave 2000Q.

###  III.1 Optimization algorithms

We run QA and the four classical algorithms for the optimization of
continuous-variable functions: Nelder-Mead (NM) [33], conjugate gradient
descent (CGD) [34], basin hopping (BH) [35], and differential evolution (DE)
[36], which are easily performed on the classical computer and are commonly
used for the optimization of continuous functions with a complex potential
landscape. While the NM and CGD methods are for local optimization, BH and DE
are for global optimization. The NM method searches for minimum by updating
the polyhedron with the observed functional value. The CGD method utilizes the
gradient of the potential function with the adjustable step size to find the
minimum. Basin hopping performs two steps iteratively: local optimization
around the proposed point and perturbation to search for the neighborhood
domain of the current solution. Differential evolution is a kind of
evolutionary computation and updates a proposed point by following the
evolutionary procedure. Details of the classical algorithms are explained in
Appendix A. Domain-wall encoding is applied only to QA on the D-Wave 2000Q,
and classical algorithms use the continuous variable directly.

We apply these algorithms to Eq. (1) using the Scipy package [43]. All
numerical simulations for the classical optimization algorithms are performed
on a laptop with Apple M1 Max chip with 10 CPU cores, which is one of the
fastest chips 222Extensive CPU benchmarks are performed in Geekbench [75]. The
single-core score of the Apple M1 Max chip in Geekbench ver. 6 is higher than
those of the Intel Xeon E5-1650 used in an old benchmark [76] and the Intel
Xeon E5-2695 v4 in a recent benchmark [77]. . For QA, we optimize Eq. (9) with
the D-Wave 2000Q, more precisely the “DW_2000Q_6” device. We have also run the
same experiments on two types of the more recent D-Wave Advantage. Since
“DW_2000Q_6” yielded the best results, we utilize it in the main text. The
comparison of performance among a few models of D-Wave devices is described in
the Appendix B.1.

###  III.2 Measured quantities

We compute two observables. The first is the ground state probability,

$\displaystyle P_{\mathrm{GS}}=\frac{\rm\\#GS}{n_{\mathrm{init}}},$ (12)

where $\rm\\#GS$ is the number of cases where the global minimum was obtained
by each algorithm, and $n_{\mathrm{init}}$ stands for the number of trials
(the number of different initial conditions). The second is the absolute error
in the energy,

$\displaystyle E_{\mathrm{abs}}=\big{|}\langle V_{\mathrm{RF}}(x)\rangle-
V_{0}\big{|},$ (13)

where the brackets represent the average over the different initial conditions
and $V_{0}=0$ is the minimum value of Eq. (1).

###  III.3 Parameter settings

For the classical optimization algorithms, we set $n_{\mathrm{init}}=10^{3}$
and fix the search space as $x\in[-3,3]$ 333In the original formulation of the
Rastrigin function, the search space is defined in $x\in[-5.12,5.12]$ [23].
Taking into account the computational precision of QA, we reduce the search
space . For a fair comparison with QA, we regard the output $x$ obtained by
the classical optimization algorithms as the solution to Eq. (1) if the
condition $|x|<\Delta x$ is satisfied. We control the execution time
$t_{\mathrm{ex}}$ by changing the maximum number of iterations of classical
algorithms as $t_{\mathrm{max}}=2^{i}$ with the integer $i$ running from 0 to
7. The execution time can be computed as
$t_{\mathrm{ex}}=\sum_{k=1}^{t_{\mathrm{max}}}\tau_{k}$ where $\tau_{k}$
represents the $k$th iteration time.

For the NM and CGD methods, the default parameter settings in Scipy are used.
For BH, we apply the CGD method in the local optimization process. For DE, the
number of populations is set to 10 and 120. Each algorithm stops on the way if
it satisfies the stopping criterion as
$|x^{t}-x^{t-1}|<10^{-7}\quad(t=1,\dots,t_{\mathrm{max}})$.

Our experimental settings are carefully designed for a fair comparison of QA
and the classical algorithms in execution time and performance. Although there
may be room for further improvements in execution time and performance, our
studies are expected to lead to a general understanding of the features of QA
for continuous function optimization against classical optimization
algorithms.

For QA, we set the system size $N=211$ and fix the coefficients as $J=1,h=2$,
and $\lambda=1$. The value $N=211$, which corresponds to $\Delta x=0.0286$
with $-x_{\rm min}=x_{\rm max}=3$, is large enough to represent the Rastrigin
function with sufficient precision in its hills and valleys. The gradient of
the potential in Eq. (1) is equal to $0$ at the global minimum $x=0$ and
yields the vanishing of the magnetic field $h_{j^{\prime}}=0$ where
$j^{\prime}$ represents the place of the zero magnetic field. The vanishing of
the magnetic field according to Eq. (11) gives two degenerate ground states of
Eq. (9) as
$\bm{\hat{\sigma}^{z}}_{\mathrm{\pm}}=(-1,\dots,-1,\hat{\sigma}_{j^{\prime}}^{z}=\pm
1,1,\dots,1)$. The ground state $\bm{\hat{\sigma}^{z}}_{\mathrm{-}}$ yields
$x_{j^{\prime}}=0$ with Eq. (5). On the other hand, the global minimum of Eq.
(1) can not be recovered from the ground state
$\bm{\hat{\sigma}^{z}}_{\mathrm{+}}$ 444By substituting $x_{j}=0$ into the
right-hand of Eq. (5), the location of the zero magnetic field can be computed
as $j^{\prime}=106$ with the present numerical values, $N=211$ and $x_{\rm
max}=-x_{\rm min}=3$. For $\hat{\sigma}_{106}^{z}=-1$, the ground state $x=0$
can be recovered as $x_{106}=0$. For $\hat{\sigma}_{106}^{z}=+1$, the domain
wall is located between $j=105$ and $j=106$, and $x_{105}=-\Delta x(\neq 0)$
results . Since the problem is negligible when $N$ is large enough, we regard
the two degenerate ground states of Eq. (9) as the equally qualified solution
of Eq. (1).

We operate the execution time of QA by varying the annealing time $1(\mu
s)\leq t_{a}\leq 500(\mu s)$. The total execution time of QA is defined as
$t_{\mathrm{ex}}=t_{p}+t_{s}$ where $t_{p}$ is programming time and $t_{s}$ is
sampling time which includes annealing time and readout time. The
$t_{\mathrm{ex}}$ can be obtained from ”qpu_access_time” in the D-Wave Ocean
SDK [42]. To compute the observables, we repeat the annealing process 20 times
independently and sample $n_{\mathrm{reads}}=10^{3}$ spin configurations for
each run. To reduce control errors, we apply flux bias calibration and gauge
transformation [47, 48]. The flux bias is estimated before each run. A gauge
is generated from a uniform distribution and fixed in each independent run.
The ground-state probability is calculated from $n_{\mathrm{reads}}=10^{3}$
spin configurations. The absolute error in the energy is computed from samples
that meet the single-domain-wall condition.

###  III.4 Result for high energy barrier

We first consider the problem with a high energy barrier $h_{0}=1.0$. We show
the dependence of the ground state probability and the absolute error in the
energy on the execution time in Fig. 2. Data points denoted as “D-Wave (mean)”
are the averaged result for $20$ independent runs. Data with the maximum
ground state probability among all runs are denoted as “D-Wave (best)”. The
error bars show the standard error computed by bootstrapping. Execution time
$t_{\mathrm{ex}}$ is always expressed in seconds.

Figure 2 (a) shows that around the region $0.7({\rm sec})\leq
t_{\mathrm{ex}}\leq 0.9({\rm sec})$, the ground state probability obtained by
“D-Wave (best)” can be larger than the other algorithms except for “DE
(popsize=120)”. The data for the NM and CGD methods are not shown in the left-
side range beyond each left-most point because they failed to find the
solution in this short-time region. Also, the data for the NM and CGD methods
as well as for QA are not shown to the right of their respective rightmost
points since they show saturation and no improvements in these long-time
region 555The NM and CGD methods get stuck in local minima before the number
of iterations of the algorithms reaches the maximum value. Therefore, an
increase in the maximum number of iterations does not affect the execution
time, and their rightmost points have the same values of the execution time.
On the other hand, BH and DE find the solution in the longer time region
$10({\rm sec})\leq t_{\mathrm{ex}}$ and their performance improves as the
execution time increases. In the shorter time region, the minimum execution
time for BH and DE cannot be reduced further than shown in Fig. 2 with our
computational resources.

Figure 2 (b) shows that the absolute error in the energy obtained by QA is
smaller than that of the NM and CGD methods and BH around the region
$0.4(\rm{sec})\leq t_{\mathrm{ex}}\leq 0.9(\rm{sec})$. Therefore, QA can find
the local minima closer to the global minima than the NM and CGD methods and
BH in this region. The outputs of BH and DE with large populations concentrate
around the global minima as we increase the execution time as we see in the
right half of Fig. 2 (b). If the number of populations is small, DE falls into
local minima even in the long execution time region. Similarly to the ground-
state probability, the absolute error in the energy of QA saturates at a
finite value.

For the higher energy barrier $h_{0}=3.0$, the qualitative behaviors of the
classical algorithms are similar to those of $h_{0}=1.0$. On the other hand,
the performance of QA deteriorates significantly, as will be described in Sec.
IV.1.

###  III.5 Result for low energy barrier

Next, we consider the problem with a lower energy barrier $h_{0}=0.2$. Figure
3 (a) shows an overall improvement of classical algorithms compared to the
case of a higher energy barrier in Fig. 2 (a). Nevertheless, the performance
of QA deteriorates slightly. The best result from QA does not reach the level
obtained by BH and DE, even with a minimum execution time of about 0.8
seconds. The saturation of the ground-state probability of QA for a longer
execution time is observed beyond about 0.9 seconds as in Fig. 2 (a) and is
thus not drawn explicitly in Fig. 3 (a).

Figure 3 (b) demonstrates that the absolute error in the energy generally
decreases compared to the case of a high-energy barrier in Fig. 2 (b). The
lowest value of the absolute error in the energy for the NM method is seen to
be slightly lower than that of the CGD method in contrast to the case of Fig.
2 (b). This result implies that the one-dimensional structure with unique
global minima fits the search strategy of the NM method. Except for this
observation, the behaviors of the absolute error in the energy remain similar
to those in Fig. 2 (b).

Figure 2:  Dependence of observables on the execution time for $k=0.5$,
$w_{0}=0.2$, and $h_{0}=1.0$. The vertical axes denote the observables: (a)
the probability of ground state and (b) the absolute error in energy. The
insets enlarge the region $0.25(\mathrm{sec})\leq t_{\mathrm{ex}}\leq
1(\mathrm{sec})$.  Figure 3: Dependence of observables on the execution time
for $k=0.5$ and $w_{0}=h_{0}=0.2$. Both axes are the same as those in Fig. 2.
The insets enlarge the region $0.25(\mathrm{sec})\leq t_{\mathrm{ex}}\leq
1(\mathrm{sec})$.

In summary, the performance of QA is slightly better than those of local
optimization algorithms and is comparable to that from global algorithms in a
limited region of execution time for the potential function with a higher
energy barrier. On the other hand, if the energy barrier is low, we did not
observe the clear benefits of QA over classical algorithms for the present
problem of continuous function optimization. Whereas increasing the execution
time enhances the performance of the global optimization algorithms, QA shows
saturation beyond intermediate time scales. In the next section, we
investigate whether this saturation originates in the problem in hardware
implementation or in other reasons including the mapping of the continuous
variable to Ising spins by the domain-wall encoding.

Figure 4: Dependence of the observables obtained by the D-Wave on the
annealing time for three values of $h_{0}$. The vertical axes are the
observables: (a) the kink density, (b) the probability satisfying the single
domain wall condition, (c) the residual energy, and (d) the ground state
probability. The horizontal axes represent the annealing time.

##  IV Comparison of QA with algorithms for discrete-variable optimization

In this section, we compare the performance of the D-Wave 2000Q with several
classical discrete-variable optimization algorithms, as well as with TEBD, a
direct implementation of noise-free coherent QA on a classical computer, to
examine how the D-Wave 2000Q compares with these other methods on the same
footing of domain-wall coding. The results may help us estimate how the
hardware would perform if the noise level were significantly reduced. The
classical algorithms used in this section are SA, SQA, and SVMC. We carefully
distinguish between coherent QA simulated by TEBD and hardware implementation
of QA, the latter being represented by the product name. For the sake of
simplicity, the D-Wave 2000Q is referred to as the D-Wave hereafter.

Table 1: The fitting exponents of the observables in the short annealing time
region: (a) the kink density and (b) the residual energy.

$h_{0}$ | D-Wave | Coherent-QA | SA | SQA
---|---|---|---|---
0.2 | $0.458\pm 0.014$ | $0.517\pm 0.003$ | $0.616\pm 0.006$ | $0.631\pm 0.014$
1.0 | $0.359\pm 0.011$ | $0.529\pm 0.003$ | $0.474\pm 0.017$ | $0.571\pm 0.01$
3.0 | - | $0.596\pm 0.008$ | $0.164\pm 0.011$ | $0.148\pm 0.026$

$h_{0}$ | D-Wave | Coherent-QA | SA | SQA
---|---|---|---|---
0.2 | $0.974\pm 0.022$ | $0.543\pm 0.010$ | $0.632\pm 0.012$ | $0.651\pm 0.015$
1.0 | $0.553\pm 0.025$ | $0.601\pm 0.006$ | $0.483\pm 0.014$ | $0.591\pm 0.011$
3.0 | - | $0.657\pm 0.002$ | $0.182\pm 0.01$ | $0.244\pm 0.015$

Table 2: The fitting exponents of the residual energy in the long annealing
time region.

$h_{0}$ | D-Wave | SA | SQA | SVMC
---|---|---|---|---
0.2 | $0.124\pm 0.03$ | $0.307\pm 0.025$ | $0.66\pm 0.044$ | $1.043\pm 0.059$
1.0 | $0.201\pm 0.02$ | $0.291\pm 0.031$ | $0.786\pm 0.032$ | $1.722\pm 0.169$
3.0 | - | - | $0.751\pm 0.009$ | -

Figure 5: Number of single kink states encoded in the lowest energy states of
the classical Hamiltonian $H_{0}$ as a function of the height of energy
barrier $h_{0}$. The parameters of the Rastrigin function are set to $k=0.5$
and $w_{0}=0.2$, and the weight of $H_{\mathrm{RF}}$ is set to $\lambda=1.0$.

We compute four quantities. The first one is the kink density

$\displaystyle\rho=\frac{1}{2N}\sum_{i=1}^{N-1}\left(1-\langle\hat{\sigma}_{i}^{z}\hat{\sigma}_{i+1}^{z}\rangle\right).$
(14)

The second is the probability $P_{\mathrm{const}}$ that the result satisfies
the condition of single domain wall, and the third one is the residual energy

$\displaystyle E_{\mathrm{res}}=\frac{\langle H_{0}\rangle-E_{0}}{N}.$ (15)

Here, the brackets $\langle\cdots\rangle$ denote the sample average and
$E_{0}$ is the ground state energy of Eq. (9). The last quantity is the ground
state probability $P_{\mathrm{GS}}$ that the ground state is successfully
found. The kink density $\rho$ represents the number of defects in spin
chains. When a single domain wall exists, $\rho=1/N$ holds. Unlike the
absolute error $E_{\mathrm{abs}}$ used in Figs. 2 (b) and 3 (b), the residual
energy $E_{\mathrm{res}}$ is calculated from all samples, including the
samples which do not satisfy the single domain wall condition. To understand
the behavior of the entire system, we do not adopt $E_{\mathrm{abs}}$ but use
$E_{\mathrm{res}}$ in this section.

Since the execution time of the D-Wave contains the measurement,
initialization and other auxiliary-processing times, we mainly use the pure
annealing time as the $x$ axis of each figure for fair comparison. The
experimental setting of the D-Wave is the same as those in Figs. 2 and 3.

###  IV.1 D-Wave

Figure 4 shows the dependence of the observables obtained from the D-Wave on
the annealing time for $h_{0}=0.2,1.0$ and 3.0. The error bars represent the
standard error obtained by bootstrapping and computed in $2\times 10^{4}$
samples for $\rho$ and $E_{\rm res}$ and in $20$ independent runs for $P_{\rm
const}$ and $P_{\rm GS}$. All error bars in figures of this section are
computed in this way. The data for $h_{0}=0.2$ and $h_{0}=1.0$ are the same as
those used in Figs. 2 and 3.

####  IV.1.1 Kink density

Let us first look at the kink density $\rho$ in Fig. 4 (a). We see that it
needs some time that the data reaches $\rho=1/N$, where the condition of
single domain wall is satisfied. In other words, this is the minimum amount of
time it takes to find the correct solution on average. If the kink density
behaves as $\rho\propto t_{a}^{-b}$ below this threshold time, the annealing
time required for the kink density to reach $\rho=1/N$ scales as $t_{a}\propto
N^{1/b}$. Therefore, the value of the exponent $b$ is important to understand
how the performance of the domain wall encoding depends on $N$, and a larger
$b$ is preferable. Except for the case of $h_{0}=3.0$, Fig. 4 (a) demonstrates
that the kink density gradually decreases as we increase the annealing time.

Each line is fitted with the polynomial function $\rho\propto t_{a}^{-b}$ in
the range $1(\mu s)\leq t_{a}\leq 10(\mu s)$. Extracted exponents are
summarized in Table 1. A smaller value of the exponent reflects the effect of
the higher energy barrier as seen in Table 1 (a). In addition, the higher
energy barrier gives a larger absolute value of $\rho$, not just a smaller
exponent. For $h_{0}=3.0$, the kink density behaves non-monotonically. Such
behavior has also been observed in several experiments on the D-Wave devices
and is believed to be caused by the open-system nature of the hardware [50,
48].

Kink densities for $h_{0}=1.0$ and $3.0$ do not reach $\rho=1/N$ and saturate
at finite values for long annealing time. Since the behaviors of all
observables in Fig. 4 do not change even though we increase the annealing time
beyond $t_{a}=500(\mu s)$, we only show the data in $1(\mu s)\leq t_{a}\leq
500(\mu s)$.

####  IV.1.2 Single domain wall condition

Figure 4 (b) illustrates that almost all samples satisfy the single domain
wall condition for $h_{0}=0.2$ beyond $t_{a}=10(\mu s)$. On the other hand,
the probability of satisfying the single domain wall condition saturates
around $P_{\mathrm{const}}\approx 0.8$ for $h_{0}=1.0$ and
$P_{\mathrm{const}}\approx 0.1$ for $h_{0}=3.0$. For those cases with higher
energy barriers, the number of samples that satisfy the single domain wall
condition deteriorates, and an excessive number of domain walls appear.

####  IV.1.3 Residual energy

Figure 4 (c) indicates that the residual energy polynomially decreases as
$E_{\mathrm{res}}\propto t_{a}^{-c}$ in the short annealing time region $1(\mu
s)\leq t_{a}\leq 10(\mu s)$ except for $h_{0}=3.0$. As with the kink density,
the exponent of $E_{\mathrm{res}}$ becomes smaller for higher energy barriers,
see Table 1 (b). Higher energy barriers yield larger absolute values of the
residual energy in addition to smaller exponents.

Next, we investigate the behaviors of $E_{\mathrm{res}}$ in the large
annealing time region $t_{a}\geq 100(\mu s)$. The scaling analysis in this
region gives us a hint on the performance of the D-Wave as an optimization
solver. We omit the data for $h_{0}=3.0$ because few samples satisfy the
single domain wall condition. We fit the polynomial function and show the
exponents in Table 2. For the problem with $h_{0}=1.0$, the slope is slightly
larger than in the case of $h_{0}=0.2$. This result indicates interestingly
that the process is more efficient in the problem with the higher energy
barrier.

Although the residual energy has a polynomial decay as found in Koh and
Nishimori [22], the exponent of $E_{\mathrm{res}}$ is not close their result
$E_{\mathrm{res}}\sim t_{a}^{-2}$. This discrepancy may come from the effects
of noise on the hardware. In Fig. 4, we omit the data for $t_{a}\geq 500(\mu
s)$ because the behavior of $E_{\mathrm{res}}$ hardly changes, and the effect
of noise increases in the longer annealing time. This result corresponds to
the fact in Figs. 2 and 3 that the observables hardly show improvements if we
increase the execution time beyond a certain value and the data are omitted
over $t_{\mathrm{ex}}\geq 0.9(\mathrm{sec})$. Therefore, the execution time
region $t_{\mathrm{ex}}\geq 0.9(\mathrm{sec})$ corresponds to the annealing
time region $t_{a}\geq 500(\mu s)$.

####  IV.1.4 Ground state probability

The ground state probability for $h_{0}=1.0$ is slightly larger than for
$h_{0}=0.2$ but remains small for $h_{0}=3.0$ as seen in Fig. 4 (d). For
$h_{0}=0.2$ and $1.0$, the residual energy shows slightly decreasing behavior
even in the large annealing time region $t_{a}\geq 100(\mu s)$, but the ground
state probability saturates around $P_{\mathrm{GS}}\approx 0.1$. This result
indicates that the transition to the ground states hardly appears, whereas the
transition to the neighborhoods of the ground states occurs in this time
scale. The ground state probability for $h_{0}=3.0$ is suppressed due to the
existence of many defects with few samples satisfying the single domain wall
condition as shown in Figs. 4 (a) and (b).

Figure 6: Results of TEBD simulations for different values of $h_{0}$ as
functions of annealing time. The vertical axes are the same as those in Fig.
4. The annealing time $t_{a}$ is in the unit of $J$ and does not correspond to
that of Fig. 4.

####  IV.1.5 Kinks in low-energy states

To understand the reason for the significant performance deterioration
observed for $h_{0}=3.0$, we examine the behavior of the low energy states of
the classical Hamiltonian $H_{0}$ in Eq. (9). To evaluate how many single-kink
states are properly encoded as low energy states of $H_{0}$, we introduce
$n_{\mathrm{enc}}$ defined as

$\displaystyle n_{\mathrm{enc}}=\min\\{j\in 0,1,\dots\mid
n_{\mathrm{kink}}(\bm{S}_{j})\neq 1\\},$ (16)

where $n_{\mathrm{kink}}(\bm{S})$ denotes the number of kinks in the spin
state $\bm{S}$, and $\bm{S}_{j}$ is the spin state corresponding to the $j$th
energy level. The value of $n_{\mathrm{enc}}$ represents the index of the
lowest energy level for which the number of kinks in the corresponding spin
state is not equal to 1. If $n_{\mathrm{enc}}>0$, the ground state is a single
kink state. If $n_{\mathrm{enc}}=0$ on the other hand, the ground state is a
multi-kink state, indicating a failure of the encoding. The value of
$n_{\mathrm{enc}}$ depends on the relative magnitude of $\lambda H_{\rm RF}$
and $H_{\mathrm{DW}}$. If we set $\lambda$ to be sufficiently small relative
to $J$, we can make $n_{\mathrm{enc}}>0$.

Figure 5 shows the $h_{0}$ dependence of $n_{\mathrm{enc}}$ obtained by
dynamic programming [51]. As $h_{0}$ increases, the values of the magnetic
field $h_{i}$ determined by Eq. (17) become larger, resulting in a smaller
$n_{\mathrm{enc}}$. For $h_{0}=0.2$ and $1.0$, we have $n_{\mathrm{enc}}=196$
and $160$, respectively, indicating that most of the single kink states are
encoded as low energy states of $H_{0}$. For $h_{0}=3.0$, we have
$n_{\mathrm{enc}}=26$, indicating that the lowest energy states of $H_{0}$ are
contaminated with multi-kink states. The reason for the extremely low
constraint satisfaction probability $P_{\mathrm{const}}$ for $h_{0}=3.0$ in
Fig. 4 (b) is considered to be that the low energy multi-kink states trap the
annealing process. Reducing $\lambda$ increases $n_{\mathrm{enc}}$, which
leads to an increase in $P_{\mathrm{const}}$. However, reducing $\lambda$ also
increases the thermal effect, and therefore, the optimization performance may
not necessarily improve. The dependence of the observables on $\lambda$ is
shown in Appendix. B.2.

####  IV.1.6 Summary

The values of observables obtained by the D-Wave depend on the height of the
energy barrier. From the viewpoint as an optimization solver, achievement of
the condition $P_{\mathrm{const}}=1$ is important. For the problem with the
higher energy barrier, the probability of satisfying the single-domain wall
condition does not reach 1 in our experiments even if we increase the
annealing time. The ground state probability also saturates at a finite value
and does improve any further. Since the decoherence effects by thermal
fluctuations may be dominant in the long annealing time region, these effects
hamper the transition to the ground states and degrade the performance.
Although the residual energy $E_{\mathrm{res}}$ shows a polynomial decay in
this region, the exponents are different from the direct simulation of the
Schrödinger dynamics [22] due to the effects of thermal noise.

Figure 7: Results of TEBD simulations. (a) Relation between the residual
energy and the deviation of the ground state probability from unity. The
dashed line denotes $(1-P_{\mathrm{GS}})(E_{1}-E_{0})$ where $E_{1}$ is the
first excited energy of $H_{0}$. (b) Deviation of the ground state probability
from unity as a function of the annealing time. (c) Residual energy and its
fitting curve as a function of the annealing time. The solid lines represent
the curves of $E_{\mathrm{res}}=C^{\prime}\exp(-Ct_{a})$ fitted to the data in
the range of $3\times 10^{3}\leq t_{a}\leq 10^{4}$. Figure 8: Local number of
kinks at $t_{a}=5\times 10^{4}$ for $\lambda=0.01$ as a function of spatial
coordinate $x_{i}=x_{\mathrm{min}}+i\Delta x$, where $i$ is the site index.
The solid line represents the potential function of Eq. (1) for $h_{0}=0.2$.
The dashed lines are a guide to the eye.

###  IV.2 TEBD for coherent QA

In order to understand the origin of the suboptimal performance of the D-Wave
described in the previous subsection, we now present the results of coherent
QA based on TEBD to elucidate the differences between the D-Wave and the ideal
QA that follows the Schrödinger dynamics without thermal noise. The TEBD
algorithm is based on the Suzuki-Trotter decomposition [40] of the time
evolution operator and the matrix product state (MPS) representation to
approximately follow the time evolution of the wave function [52, 53]. Time
and energy are measured in the unit of $J$ in Eq. (8) and we adopt the natural
unit where $\hbar=1$. In this unit, if we choose $B(1)/2$ for the D-Wave as
$J$, then $t_{a}=1$ corresponds to approximately 0.026 ns. The time step size
of the Suzuki-Trotter decomposition is set to $\Delta t=0.025/J$. The maximum
bond dimension of the MPS, which controls the accuracy of the approximation,
is $32$. We confirmed that the results converged well with respect to these
parameters. The annealing schedules are the linear function, as $A(s)/2=1-s$
and $B(s)/2=s$.

Figure 6 illustrates the results obtained by TEBD. A significant difference
from the experimental results by the D-Wave in Fig. 4 is that the behavior of
the kink density $\rho$ exhibits only very slight dependence on the height of
the energy barrier $h_{0}$, which may originate in the weak dependence of
coherent quantum tunneling effects on the barrier height. Moreover, in the
case of coherent QA, for $t_{a}\geq 10^{3}$, the single domain wall condition
is almost completely satisfied, and the residual energy begins to decay quite
rapidly, while in the case of the D-Wave in Fig. 4, the single domain wall
condition is not satisfied even in the long time regime, and the residual
energy does not decay. These results suggest that deviations from the ideal
Schrödinger dynamics in the D-Wave lead to low performance, particularly for
high $h_{0}$.

In Figs. 6 (a) and (c), the kink density and the residual energy have the
polynomial decay in the region $t_{a}\leq 100$, and the fitting exponents are
summarized in Table 1. In the limit where $N$ is large and $\lambda$ is small,
the theory of the Kibble-Zurek mechanism [54, 55, 56] predicts that
$\rho\propto t_{a}^{-b}$ and $E_{\mathrm{res}}\propto t_{a}^{-c}$ with
$b=c=1/2$ for the one-dimensional transverse-field Ising model without
longitudinal fields. In the simulation for Fig. 6, $\lambda$ is set to 1, and
therefore the magnitude of the longitudinal magnetic field is not negligible,
see Eqs. (7) and (9), but the exponents show relatively mild dependence on
$h_{0}$.

As shown in Figs. 6 (b) and 6 (d), for $t_{a}>100$, the constraint
satisfaction probability $P_{\mathrm{const}}$ starts to increase from 0 and
the ground state probability $P_{\mathrm{GS}}$ also begins to increase. In
this region, the residual energy deviates from the polynomial behavior and
decreases more rapidly.

In a longer annealing time region ($t_{a}>2\times 10^{3}$), where the single
domain wall condition is completely satisfied, $P_{\mathrm{GS}}$ approaches 1.
In this regime, the energy expectation value can be expressed as

$\displaystyle\langle
H_{0}\rangle=E_{0}P_{\mathrm{GS}}+\mathcal{O}(1-P_{\mathrm{GS}})\approx
E_{0}+\mathcal{O}(1-P_{\mathrm{GS}}),$ (17)

and therefore the residual energy is proportional to $1-P_{\mathrm{GS}}$ as

$\displaystyle E_{\mathrm{res}}=\frac{\langle H_{0}\rangle-
E_{0}}{N}=\mathcal{O}(1-P_{\mathrm{GS}}).$ (18)

We can confirm this behavior by plotting the relation between
$E_{\mathrm{res}}$ and $1-P_{\mathrm{GS}}$ in Fig. 7 (a), as the slope is
approximately 1 in the region where $1-P_{\mathrm{GS}}<0.1$. The annealing
time dependence of $1-P_{\mathrm{GS}}$ is shown in Fig. 7 (b), and in the
region where $1-P_{\mathrm{GS}}<0.1$, it can be seen that $1-P_{\mathrm{GS}}$
follows an exponential decay as $1-P_{\mathrm{GS}}\propto\exp(-Ct_{a})$, which
is reminiscent of the Landau-Zener formula. Consequently, the annealing time
dependence of the residual energy can be expressed as
$E_{\mathrm{res}}=C’\exp(-Ct_{a})$, where $C^{\prime}$ and $C$ are positive
constants. The solid lines in Fig. 7 (c) represent the result of fitting this
expression to the data in the range of $3\times 10^{3}\leq t_{a}\leq 10^{4}$.
In Ref. [22], an analysis of QA by direct numerical solution of the time-
dependent Schrödinger equation for a single particle system in a continuous
space showed that for a linear annealing schedule, the residual energy decays
as the inverse square of the annealing time beyond the Landau-Zener regime as
suggested in Ref. [57]. The deviation of the residual energy for $h_{0}=0.2$
from the exponential behavior in the longer-time region $t_{a}>10^{4}$ may be
an indication of a transition to such a power law behavior, although the
exponent is much larger than two. To confirm the existence of this transition
and to determine the power law exponents, simulations for much longer
annealing times are necessary, which is challenging due to the accumulation of
numerical errors in TEBD.

Figure 9: Dependence of the observables obtained from SA on the annealing time
for different values of $h_{0}$. The vertical axes are the same as those in
Fig. 4. The horizontal axes represent the Monte Carlo step.

To understand the effect of the height of the energy barrier on the coherent
QA, we analyze the spatial distribution of kinks in the wave function after
the process of QA has terminated. Figure 8 shows the local number of kinks
$(1-\langle\hat{\sigma}_{i}^{z}\hat{\sigma}_{i+1}^{z}\rangle)/2$ at
$t_{a}=5\times 10^{4}$ as a function of spatial coordinate, as well as the
values of the ground state probability $P_{\mathrm{GS}}$ and residual energy
$E_{\mathrm{res}}$. For a small energy barrier $h_{0}=0.2$, kinks are
distributed around the global minimum of the potential function at $x=0$
without significantly avoiding the maxima. This result suggests that for the
low energy barrier, the kink is more likely to move around, making it
difficult to remain at the bottom of the potential minima. On the other hand,
for a large energy barrier $h_{0}=3.0$, kinks localize in the vicinity of the
minima and hardly exists in the vicinity of the maxima due to the high energy
barrier. The increase in the population at the local (but not global) minima
for $h_{0}=3.0$ leads to larger residual energy than that of $h_{0}=0.2$.
However, the probability of obtaining the ground state $P_{\mathrm{GS}}$ is
higher when $h_{0}=3.0$. This behavior suggests that a lower energy barrier
case is not necessarily easier to optimize, contrary to naive intuition. The
characteristic of kink spatial distribution discussed here becomes indistinct
when $\lambda$ is larger, but it can still be observed. Indeed, for
$\lambda=1$, the optimization performance does not decrease even in the case
of high energy barrier as shown in Fig. 6.

In summary, we have found that the optimization performance of coherent QA
realized by TEBD is almost independent of the height of energy barrier
$h_{0}$, and the system converges to the ground state with high probability
for sufficiently long annealing times. Therefore, the significant dependence
of performance on $h_{0}$ and the premature saturation of the ground state
probability in the data of the D-Wave can be attributed to imperfect
realization of QA in the hardware, possibly mainly by thermal noise. As
discussed in Appendix C, in order to fully understand the thermal effects on
the hardware, it is necessary to analyze an open quantum system, not just an
isolated system as done in this study. The performance of QA with a more
coherent device would approach the results of TEBD simulations presented in
this subsection and can be expected to be equivalent to or even better than
that of global optimization algorithms such as BH and DE. Additional evidence
will be provided in the following subsections to show that the $h_{0}$
dependence of data originates in thermal effects.

###  IV.3 Simulated annealing

Next, we consider SA to investigate the effects of thermal fluctuations. In
SA, we use a single-spin update with the Metropolis rule. The temperature is
linearly decreased as

$\displaystyle
T(t)=T_{0}+\frac{T_{1}-T_{0}}{t_{\mathrm{MCS}}-1}\,t\quad(t=0,\dots,t_{\mathrm{MCS}}-1),$
(19)

with $T_{0}=1$ and $T_{1}=10^{-5}$. The maximum number of Monte Carlo steps
(MCS) is denoted as $t_{\mathrm{MCS}}$. The number of spin updates is
$t_{\mathrm{MCS}}N$. We repeat SA $20$ times independently and sample $10^{3}$
spin configurations for each independent run.

####  IV.3.1 Kink density

Figure 9 shows the results obtained by SA. As in Fig. 4 (a), Fig. 9 (a)
indicates the power law behaviors of the kink density in the short annealing
time region $2\leq t_{\mathrm{MCS}}\leq 100$. Unlike the results obtained from
the D-Wave in Fig. 4 (a), the power law behaviors can be seen even in the case
of high energy barrier, $h_{0}=3.0$. For $h_{0}=0.2$, the exponents are larger
than those of the results obtained by the D-Wave and the coherent QA, see
Table 1. On the other hand, for larger $h_{0}$ with higher energy barrier, the
exponents are larger than those of the D-Wave and smaller than those of the
coherent QA for the problem. In this time region, the values of the kink
density become larger than those of the D-Wave and take similar values to the
case of coherent QA except for $h_{0}=3.0$. For $h_{0}=3.0$, the kink density
does not reach $\rho=1/N$ even at $t_{\mathrm{MCS}}=10^{5}$. We need more
annealing time to reduce the defects for the problem with the larger value of
$h_{0}$.

####  IV.3.2 Single domain wall condition

Figure 9 (b) illustrates that few samples satisfy the single domain wall
condition in the short-time region where the kink density has the polynomial
decay. While this result is similar to the result obtained by the coherent QA
in the short-time region as seen in Fig. 6 (b), a clear difference from the
D-Wave in Fig. 4 (b) is observed. Since the system is not sufficiently
equilibrated in the short annealing time region, many defects remain. For the
D-Wave result in Fig. 4 (b), $P_{\mathrm{const}}$ has a finite value except
for $h_{0}=3.0$ in this region. Although this D-Wave result is similar to the
behavior of SA in the intermediate annealing time region $10^{2}\leq
t_{\mathrm{MCS}}\leq 10^{3}$ in Fig. 9 (b), it is impossible to reproduce all
D-Wave data quantitatively by SA, for example, the residual energy as shown
below.

In the long annealing time region $t_{\mathrm{MCS}}\geq 10^{4}$, almost all
samples satisfy the single domain wall condition except for $h_{0}=3.0$. For
$h_{0}=3.0$, $P_{\mathrm{const}}$ is suppressed. Similarly to the D-Wave data
in Fig. 4 (b), SA is also affected by the height of the energy barrier.

####  IV.3.3 Residual energy

As with the kink density in Fig. 9 (a), Fig. 9 (c) shows that the residual
energy scales polynomially in the short annealing time region $2\leq
t_{\mathrm{MCS}}\leq 100$. The exponents are smaller than the D-Wave results
as seen in Table 1 (b). The exponents are larger than those by the coherent QA
for $h_{0}=0.2$ and smaller for $h_{0}=1.0$ and $3.0$. These results imply
that the dynamics based on thermal fluctuations are hampered by the higher
energy barrier, and quantum fluctuations are more efficient than thermal
fluctuations for higher barriers.

The residual energy of SA decreases polynomially in the long annealing time
regions $2\times 10^{3}\leq t_{\mathrm{MCS}}\leq 10^{5}$ for $h_{0}=0.2$ and
$10^{4}\leq t_{\mathrm{MCS}}\leq 10^{5}$ for $h_{0}=1.0$ where
$P_{\mathrm{const}}\approx 1$ holds. According to Table 2, the slope of the
fitting line becomes smaller for the problem with the higher energy barrier.
The higher energy barrier impedes the transition to the ground states by
thermal fluctuations. On the other hand, the exponents obtained by the D-Wave
show the opposite trend. Note that $P_{\mathrm{const}}=1$ does not hold for
$h_{0}=1.0$ in Fig. 4 (b) for the D-Wave. Whether or not this is related to
the quantumness of the D-Wave is not clear at this point. Although the
exponents of SA are better than those of the D-Wave, the inverse square
scaling of $E_{\mathrm{res}}$ is not observed.

####  IV.3.4 Ground state probability

Figure 9 (d) shows that the ground state probability starts to increase
monotonically from $t_{\mathrm{MCS}}=10^{2}$. The performance depends on the
height of the energy barrier. The higher the energy barrier is, the smaller
the success probability becomes. The behavior is different from that in Figs.
4 (d) and 6 (d). Since the ground state probability is not close to 1 even in
the long annealing time region, the exponential decay of $E_{\mathrm{res}}$
observed in the coherent QA in Fig. 7 (c) does not appear in SA as seen in
Fig. 9 (c).

####  IV.3.5 Summary

The observables of SA clearly depend on $h_{0}$ while the observables of the
coherent QA hardly depend on $h_{0}$. This is a significant difference between
SA and coherent QA. In coherent QA, quantum tunneling effects play an
important role in the state transition between local minima. On the other
hand, SA shows deteriorated performance under higher energy barriers. This
behavior is consistent with the physical picture that SA goes over the energy
barrier with thermal fluctuations. Also in the case of the D-Wave, dependence
of the observables on $h_{0}$ can be clearly seen. Premature saturation of
$P_{\mathrm{const}}$ and $P_{\mathrm{GS}}$ at lower values appear in the long
annealing time region for the D-Wave. Such saturation of these observables
does not exist in SA. This suggests that the degraded performance of the
D-Wave may not be due to thermal noise alone, but to other factors as well
[58, 59].

Figure 10: Dependence of the observables obtained by SQA on the annealing time
for different values of $h_{0}$. Both axes are the same as those in Fig. 4.

###  IV.4 Simulated quantum annealing

Next, we study SQA and see its difference from SA and other protocols in some
detail. Simulated quantum annealing is the algorithm to simulate some aspects
QA on a classical computer. In SQA, we apply the Suzuki-Trotter decomposition
[40] to the partition function of Eq. (2). The quantum system at a finite
temperature is then mapped to the effective classical Ising model as

$\displaystyle H_{\mathrm{eff}}=$
$\displaystyle\frac{B(s)}{2M}\sum_{k=1}^{M}H_{0}(k)+\frac{M}{2\beta}\ln\tanh\left(\frac{\beta
A(s)}{2M}\right)\sum_{k=1}^{M}\sum_{i=1}^{N}\sigma_{ik}\sigma_{ik+1},$ (20)
$\displaystyle
H_{0}(k)=\lambda\sum_{i=1}^{N}h_{i}\sigma_{ik}-J\sum_{i=1}^{N-1}\sigma_{ik}\sigma_{i+1k}+h\left(\sigma_{1k}-\sigma_{Nk}\right),$
(21)

where $M$ is the Trotter number, $\beta=1/T$ is the inverse temperature, and
$\sigma_{ik}\in\\{\pm 1\\}(i=1,\dots,N,k=1,\dots,M)$ is the classical Ising
configuration. The coefficients of the problem and transverse field are
linearly controlled as

$\displaystyle\frac{A(s)}{2}$
$\displaystyle=1-s=1-\frac{t}{t_{\mathrm{MCS}}-1},$ (22)
$\displaystyle\frac{B(s)}{2}$ $\displaystyle=s=\frac{t}{t_{\mathrm{MCS}}-1},$
(23)

where $t=0,\dots,t_{\mathrm{MCS}}-1$. We set the effective temperature
$\beta/M=1$ and $M=10^{3}$. We utilize a single-spin update with the
Metropolis rule. We sweep the spin configurations $MN$ times for each MCS. We
perform SQA 20 times independently. As shown in Ref. [60], the final spin
configurations in Trotter replicas are utilized for computing the observables.

####  IV.4.1 Kink density

Figure 10 shows the results obtained from SQA. In Fig. 10 (a), the kink
density does not depend on the value of $h_{0}$ in the shortest annealing time
region $2\leq t_{\mathrm{MCS}}\leq 3$. This behavior is similar to that of the
coherent QA in Fig. 6 (a). In the short annealing time region $4\leq
t_{\mathrm{MCS}}\leq 100$, the polynomial decay of $\rho$ can be seen. Except
for $h_{0}=3.0$, SQA has the largest exponent in annealing protocols for
discrete-value optimization, as shown in Table 1. The exponents of $\rho$
become smaller as the height of the energy barrier increases. The tendency is
shared by other annealing protocols except for the coherent QA. For
$h_{0}=1.0$ and $3.0$, SQA requires less annealing time to satisfy the single
domain wall condition as compared to SA, and thus classically-simulated
quantum fluctuations are more effective than thermal fluctuations in the
higher barrier case. This result reflects the difference in the dynamics of
both methods and can also be understood from the behaviors of
$P_{\mathrm{const}}$ in Fig. 10 (b). For $h_{0}=0.2$, we cannot observe a
clear difference between SA and SQA. Compared with the coherent QA, the
dependence of $\rho$ on $h_{0}$ appears. The dynamics of SQA is fundamentally
different from the dynamics of coherent QA [61] and may be regarded to be
rather closer to the D-Wave behavior and SA.

####  IV.4.2 Single domain wall condition

Figure 10 (b) shows that few samples satisfy the single domain wall condition
in the short-time region similarly to the cases of SA and coherent QA in Figs.
9 (b) and 6 (b). The condition$P_{\mathrm{const}}\approx 1$ holds in the long
annealing time region $t_{\mathrm{MCS}}\geq 10^{4}$ including the case of
$h_{0}=3.0$. This behavior is different from the D-Wave and SA in Figs. 4 and
9. No suppression of $P_{\mathrm{const}}$ for $h_{0}=3.0$ is observed for SQA,
which is similar to coherent QA.

####  IV.4.3 Residual energy

Figure 10 (c) illustrates that the residual energy decreases polynomially in
the short annealing time region. The range of fitting is the same as that in
Fig. 10 (a) to extract the values in Table 1. As seen in this Table, as the
height of the energy barrier increases, the exponents of SQA decrease. The
exponents of SQA are larger than those of SA. The result indicates that SQA
can solve the problem more efficiently than SA. For $h_{0}=3.0$, the exponent
of SQA is smaller than that of the coherent QA. Although SQA can simulate some
aspects of quantum effects, the higher energy barrier degrades the
performance, unlike the coherent QA. In the long annealing time region
$10^{3}\leq t_{\mathrm{MCS}}\leq 10^{4}$, the polynomial decay of
$E_{\mathrm{res}}$ can be observed. As seen in Table 2, SQA can solve the
problem more efficiently than the D-Wave and SA since SQA gives larger values
of the exponent than the D-Wave and SA. In our experiments, the exponential
decay of $E_{\mathrm{res}}$ as in coherent quantum annealing is not observed
in SQA. To see the behavior, we may need to investigate longer annealing time
regions where $P_{\mathrm{GS}}\approx 1$ holds.

Figure 11: Dependence of the observables attained by SVMC on the annealing
time for different values of $h_{0}$. Both axes are the same as those in Fig.
4.

####  IV.4.4 Ground state probability

Figure 10 (d) shows that the success probability becomes smaller as the height
of the energy barrier increases in the large annealing time region
$t_{\mathrm{MCS}}\geq 10^{4}$. Since the values of $P_{\mathrm{GS}}$ obtained
by SQA are larger than those by SA in the region $t_{\mathrm{MCS}}\geq 2\times
10^{3}$, the effect of the energy barrier is weakened by classically-simulated
quantum fluctuations.

####  IV.4.5 Summary

The observables obtained by SQA depend on the value of $h_{0}$, which is
different from the case of coherent QA. This attributes to the fact that the
dynamics of SQA does not follow the Schrödinger dynamics but is based on
classical Monte Carlo updates under thermal fluctuations. Compared with SA,
SQA suffers less performance degradation under the higher energy barrier,
possibly by partly simulating tunneling effects [62].

###  IV.5 Spin-vector Monte Carlo

Spin-vector Monte Carlo is a classical algorithm proposed to explain the
outputs from the D-Wave device [41]. In SVMC, a quantum spin is replaced by a
classical rotor of unit length with a continuous angle variable
$\theta_{i}\in[0,\pi]$ as

$\displaystyle\hat{\sigma}_{i}^{z}\rightarrow\cos\theta_{i},\quad\hat{\sigma}_{i}^{x}\rightarrow\sin\theta_{i}.$
(24)

The total Hamiltonian can then be described as :

$\displaystyle H$
$\displaystyle=\frac{B(s)}{2}H_{0}-\frac{A(s)}{2}\sum_{i=1}^{N}\sin\theta_{i},$
(25) $\displaystyle H_{0}$
$\displaystyle=\lambda\sum_{i=1}^{N}h_{i}\cos\theta_{i}-J\sum_{i=1}^{N-1}\cos\theta_{i}\cos\theta_{i+1}$
$\displaystyle+h\left(\cos\theta_{1}-\cos\theta_{N}\right).$ (26)

The annealing schedules of $A(s)/2$ and $B(s)/2$ are the same as those in SQA.
The temperature in Monte Carlo process is set in our simulation to $T=10^{-5}$
and the number of all spin updates is written as $Nt_{\mathrm{MCS}}$. Each
continuous angle is initialized to $\theta_{i}=\pi/2$, and spins are randomly
chosen and updated with the Metropolis rule as
$\theta_{i}\rightarrow\theta_{i}^{\prime}\in[0,\pi]$. In addition, we apply
the transverse-field-dependent (TFD) update [63, 64], which is proposed to
capture the freezing phenomena in the D-Wave device. In the TFD update, the
current angle is updated as $\theta_{i}\rightarrow\theta_{i}+\epsilon_{i}(s)$
where $\epsilon_{i}(s)$ is the random variable dependent on the annealing
parameter and the energy coefficients,

$\displaystyle\epsilon_{i}(s)\in\left[-\mathrm{min}\left(1,\frac{A(s)}{B(s)}\right)\pi,\mathrm{min}\left(1,\frac{A(s)}{B(s)}\right)\pi\right].$
(27)

The range of updating the angle gradually decreases in the late annealing
process $A(s)<B(s)$. We execute SVMC $10$ times independently and sample
$10^{3}$ spin configurations for each run.

####  IV.5.1 Kink density

Figure 11 shows the results of SVMC. In Fig. 11 (a), nonlinear behaviors of
$\rho$ are observed. In the region of very short annealing time $2\leq
t_{\mathrm{MCS}}\leq 3$, the kink density does not depend on the values of
$h_{0}$. Dependence of $\rho$ on the values of $h_{0}$ appears for
$t_{\mathrm{MCS}}\geq 3$. For $h_{0}=3.0$, the kink density saturates at a
larger value than that of SA in the region $t_{\mathrm{MCS}}\geq 50$. This
result represents that SVMC is more affected by the high energy barrier than
SA.

####  IV.5.2 Single domain wall condition

Figure 11 (b) shows that most samples satisfy the single domain wall condition
$t_{\mathrm{MCS}}\geq 10^{5}$ for $h_{0}=0.2$ and $1.0$. Compared with other
annealing protocols, SVMC needs more annealing time to satisfy
$P_{\mathrm{const}}\approx 1$. The result indicates that the dynamics of SVMC
is slower than that of other annealing protocols.

####  IV.5.3 Residual energy

As in $\rho$ in Fig. 11 (a), Fig. 11 (c) illustrates that nonlinear behavior
of $E_{\mathrm{res}}$ emerges and the saturation of $E_{\mathrm{res}}$ for
$h_{0}=3.0$ appears. In the long annealing time region, the polynomial decay
of $E_{\mathrm{res}}$ can be observed. The fitting regions are $5\times
10^{4}\leq t_{\mathrm{MCS}}\leq 5\times 10^{5}$ for $h_{0}=0.2$ and
$10^{5}\leq t_{\mathrm{MCS}}\leq 5\times 10^{5}$ for $h_{0}=1.0$. In these
regions, $P_{\mathrm{const}}\approx 1$ holds as seen in Fig. 11 (b). From
Table 2, the exponents of SVMC are the largest among all annealing protocols.
Therefore, SVMC can reach the ground states faster than other annealing
protocols in the long time region. Compared with SQA in the region where the
single domain wall condition holds, the values of $E_{\mathrm{res}}$ obtained
by SVMC are smaller for $h_{0}=0.2$ and $1.0$.

####  IV.5.4 Ground state probability

Figure 11 (d) demonstrates that the ground state probability starts to
increase more slowly than SA and SQA. The higher energy barrier suppresses the
increase of $P_{\mathrm{GS}}$. At a later stage, SVMC produces a larger value
of $P_{\mathrm{GS}}$ than SQA does for $h_{0}=0.2$ and $1.0$ in the region
where $P_{\mathrm{const}}\approx 1$ holds.

Figure 12: Dependence of the observables on the execution time for different
annealing protocols with $k=0.5$, $w_{0}=0.2$, and $h_{0}=1.0$. Both axes are
the same as those in Fig. 4.

####  IV.5.5 Summary

Values of observables in SVMC depend on the value of $h_{0}$ as in SA and SQA.
The dynamics of SVMC is more affected by the higher energy barrier. This is a
consequence of the fact that the dynamics of SVMC is activated by thermal
fluctuations in the semi-classical potential function. Although the behaviors
of the observables are unique to SVMC and are different from those obtained by
the D-Wave, the performance of SVMC is the best in the region where
$P_{\mathrm{const}}\approx 1$ holds except for the coherent QA. This result
suggests that the dynamics of SVMC is more effective than other classical
annealing protocols for the Rastrigin function in the space where only a
single domain wall exists.

Considering SVMC as a practical solver for continuous function optimization, a
smaller absolute value of the execution time is required, not just larger
scaling exponents. In the next section, we analyze the dependence of
observables on the absolute execution time.

###  IV.6 Actual execution time

Finally, we show the dependence of observables on the actual execution time in
seconds for different annealing protocols in Fig. 12. We use the same data as
those presented in this Sec. IV for $k=0.5$, $w_{0}=0.2$, and $h_{0}=1.0$. The
definition of $t_{\mathrm{ex}}$ is the same as that defined in Sec. III.2.
Since the TEBD is the deterministic algorithm and different from the other
annealing protocols based on the Monte Carlo method, we avoid to show the TEBD
data for direct comparison in this section.

Figure 12 (a) illustrates that SQA reaches $\rho=1/N$ most quickly among
different annealing protocols 666 Notice that the D-Wave becomes stuck at a
value larger than $\rho=1/N$ as shown in Fig. 4 (a) though this behavior is
not clearly seen in the scale of Fig. 12 (a).. Although the D-Wave reaches the
lowest value of $\rho$ around $0.4({\rm sec})\leq t_{\mathrm{ex}}\leq 0.9({\rm
sec})$, its value of $\rho$ remains suboptimal in the longer time region due
to the noise and other imperfections as observed more clearly in Fig. 4 (a).

Around the moderately-short time region $0.4({\rm sec})\leq
t_{\mathrm{ex}}\leq 0.9({\rm sec})$, $P_{\mathrm{const}}$ reaches close to 1
by the D-Wave as shown in Fig. 12 (b). If the goal is to satisfy the condition
of single domain wall as quickly as possible in the actual time of execution,
the D-Wave is the choice.

Figure 12 (c) shows that the D-Wave has the smallest value of
$E_{\mathrm{res}}$ around $0.4({\rm sec})\leq t_{\mathrm{ex}}\leq 0.9({\rm
sec})$. Around $2({\rm sec})\leq t_{\mathrm{ex}}\leq 370({\rm sec})$, SQA
performs best. From the analysis of the exponents in Table 2, SVMC surpasses
SQA in the long time region, and the difference of SQA and SVMC widens with
time. To satisfy the single domain wall condition, SVMC needs more time, see
Figs. 12 (a) and (b). However, once the condition is satisfied, the residual
energy decreases most rapidly by SVMC.

Figure 12 (d) indicates that the success probability of the D-Wave has the
largest value around $0.4({\rm sec})\leq t_{\mathrm{ex}}\leq 0.9({\rm sec})$,
$P_{\mathrm{GS}}\approx 0.124$ as indicated by the horizontal line. The D-Wave
reaches this value most quickly among all annealing protocols. However, it
shows saturation at around this finite value, and the performance hardly
improves with time. In Fig. 2, we see that the ground state probability
obtained by BH and DE are large $P_{\mathrm{GS}}>0.9$ even in the longest time
region $10({\rm sec})\leq t_{\mathrm{ex}}\leq 100({\rm sec})$. As far as our
implementation is concerned, other annealing protocols, SA, SQA, SVMC, are
similarly less effective than the global optimization algorithms for
continuous-variable optimization since the former do not achieve such high
values in this long-time region.

Overall, in the short time window of $0.4({\rm sec})\leq t_{\mathrm{ex}}\leq
0.9({\rm sec})$, the D-Wave is the best choice. At around the intermediate
time $t_{\mathrm{ex}}\approx 2({\rm sec})$, SQA is competitive. Beyond this
time region, dedicated global optimization algorithms for continuous
variables, DE and BH, keep improving with time, whereas classical
metaheuristics for discrete optimization remain relatively mediocre. The
D-Wave stays stagnant. It is nevertheless expected from the data of TEBD that
the ideal QA without noise is likely to be most competitive, if realized at
the hardware level, also in the intermediate and long-time regions.

##  V Conclusion

We have performed extensive benchmarks of QA and classical optimization
algorithms for the continuous-variable Rastrigin function in one dimension
with a rugged energy landscape. By using the domain wall encoding, the
continuous variable problem was represented by a one-dimensional Ising chain.

We have first compared the D-Wave with classical algorithms for continuous-
variable optimization with local or global updates and have demonstrated
quantitatively that, in the case of higher energy barrier, the performance of
the D-Wave is better or comparable to that of the local optimization
algorithms in a limited region of execution time and is comparable to the
global optimization algorithms in the same time region. For longer execution
times, the global algorithms show outstanding performance whereas the D-Wave
remains stagnated. For the potential function with lower energy barrier, no
clear advantage of the D-Wave has been observed. While the performance of the
global optimization algorithms keeps improving with increasing execution time,
no improvements were identified in the D-Wave.

For better understanding of the performance of QA, in its hardware realization
of the D-Wave as well as in an ideal setting without noise, we next compared
QA with protocols for discrete-variable optimization, TEBD for coherent QA,
SA, SQA, and SVMC. The observed values of physical quantities show a clear
dependence on the height of the energy barrier except for TEBD for coherent
QA. This is natural for SA, SQA, and SVMC because these classical algorithms
make explicit use of thermal noise to go over the barrier. In the case of the
D-Wave, its barrier-height dependence implies that the process realized in the
D-Wave is explicitly affected by thermal noise. Particularly in the case of
the highest barrier with $h_{0}=3.0$, we found no sign of the successful
solution by the D-Wave, whereas coherent QA by TEBD reached the ground state
with high probability for relatively long annealing time with little-to-no
dependence on the barrier height. This result suggests that the performance
degradation of the D-Wave at long annealing time may be significantly reduced
by improvements in the hardware to reduce noise and possibly other
imperfections as demonstrated in recent experiments of quantum simulation [48,
66]. Stated otherwise, quantum annealing has the potential to become quite
competitive even for continuous-variable problems under the domain wall
encoding. This is achieved by quantum tunneling effects, as implied by the
independence of the barrier height observed in the TEBD data. Classical
algorithms, including those for continuous-variable optimization, would
struggle to find a solution if the barrier become even higher than we tested,
while QA, if realized coherently, would be much less affected by the height.
This conclusion has become possible by the analysis of the continuous-variable
problem with an explicit and easy control of the barrier height, which is
generally difficult to achieve in problems expressed directly by the Ising
model. In other words, we have demonstrated the effectiveness of coherent
quantum tunneling effects explicitly and quantitatively in an optimization
problem. Simulated quantum annealing may be regarded as a quantum-inspired
classical algorithm that partly takes into account quantum effects.
Nevertheless, comparison of data for TEBD in Fig. 6 and SQA in Fig. 10
suggests that a complete realization of quantum effects in QA is necessary for
better results.

This paper represents a first step toward systematic and quantitative studies
of continuous-variable optimization by QA in comparison with a series of well-
established classical algorithms. In particular, the present one-dimensional
problem with a periodic structure of the barrier is one of the simplest
examples of continuous-variable optimization. The global optimization
algorithm DE based on the genetic algorithm may find it more difficult to
reach the correct solution for aperiodic potentials than for the periodic
Rastrigin function. Also, higher-dimensional problems with the complex energy
landscape as studied qualitatively in Refs. [25, 26, 27], which are generally
difficult to analyze by TEBD, may be the interesting next target of systematic
studies.

## Acknowledgment

The work of H.N. is based on a project JPNP16007 commissioned by the New
Energy and Industrial Technology Development Organization (NEDO).

##  Appendix A Classical algorithms for continuous-variable optimization

In this Appendix, we explain in some detail the classical optimization
algorithms used for benchmarking. We consider the optimization problem as

$\displaystyle\min_{\bm{x}\in\mathbb{R}^{n}}f(\bm{x}),$ (28)

where the objective function is defined by
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$. The implementation of all classical
algorithms for benchmarking utilizes Scipy [43].

###  A.1 Nelder-Mead

The Nelder-Mead (NM) method is commonly used for unconstrained continuous
optimization problems [33, 67]. The NM method does not need any information on
derivatives and utilizes only the observed functional value. In the NM method,
we hold a simplex constructed from the convex hull of $n+1$ points in
$\mathbb{R}^{n}$. We iteratively update the simplex by removing the worst
functional point and substituting another better point for it. The candidate
is selected by the following process: reflecting, expanding, or contracting
the simplex along the line connecting the worst point with the centroid of the
other points [68]. If the solution cannot be improved anymore, the simplex is
shrunk as all other points get closer to the current best solution. The
detailed algorithm of the NM method is explained in Refs. [68, 67].

###  A.2 Conjugate gradient descent

The conjugate gradient decent (CGD) method is an iterative algorithm with the
gradient of the objective function [68]. The original algorithm is designed
for solving a linear system of equations [34]. The extended algorithm for
nonlinear functions was proposed by Fletcher and Reeves [69]. In the CGD
method, we set the initial point $\bm{x}^{(1)}\in\mathbb{R}^{n}$ and its
direction $\bm{d}^{(1)}=-\nabla f(\bm{x}^{(1)})$, and update the current point
$\bm{x}^{(t)}$ as follows:
$\bm{x}^{(t+1)}\leftarrow\bm{x}^{(t)}+\alpha^{(t)}\bm{d}^{(t)}$ where
$\alpha^{(t)}$ is the step length obtained via a line search for minimizing
$f(\bm{x}^{(t)})$ along the direction $\bm{d}^{(t)}$. The current direction is
calculated by using the gradient as $\bm{d}^{(t+1)}\leftarrow-\nabla
f(\bm{x}^{(t+1)})+\beta^{(t+1)}\bm{d}^{(t)}$. In the Fletcher and Reeves
formulations, the hyperparameter $\beta^{(t+1)}$ is determined by
$\beta^{(t+1)}\leftarrow\nabla f(\bm{x}^{(t+1)})^{T}\nabla
f(\bm{x}^{(t+1)})/\nabla f(\bm{x}^{(t)})^{T}\nabla f(\bm{x}^{(t)})$. We
iterate the above procedure until the algorithm converges to a local minimum
or the number of iterations reaches the maximum iteration.

###  A.3 Basin-Hopping

The Basin-hopping (BH) method is a widely used global optimization algorithm
for nonlinear multimodal optimization problems, where the local optimization
algorithm often fails to find the global minimum due to the rugged energy
landscape [35]. The BH method is combined with the gradient-based local
optimization and the perturbation based on Monte Carlo random walks with the
Metropolis rule [70, 71]. Starting from an initial point $\bm{x}^{(1)}$, we
apply the local optimization to map a current point $\bm{x}^{(t)}$ to the
nearest minima $\bm{y}^{(t)}$ in the local optimization process. Next, we add
a perturbation to the current solution as
$\bm{x}_{\mathrm{new}}=\bm{y}^{(t)}+\bm{\epsilon}$ where $\bm{\epsilon}$ is a
random coordinate generated from the uniform distribution $U(a,b)$. We perform
the local optimization once again and attain another local minimum
$\bm{y}_{\mathrm{new}}$. Following the Metropolis rule, we accept the proposed
point as $\bm{y}^{(t+1)}=\bm{y}_{\mathrm{new}}$ or reject it as
$\bm{y}^{(t+1)}=\bm{y}^{(t)}$ with the probability,
$\mathrm{min}\left[1,\exp(-(f(\bm{y}_{\mathrm{new}})-f(\bm{y}^{(t)}))/T)\right]$
where $T$ is the temperature to adjust the scale of the energy. The above
iteration continues until $t=t_{\mathrm{max}}$. In our experiments, we utilize
the default setting $a=-0.5$, $b=0.5$, and $T=1$, and adopt the CGD method in
the local optimization process.

###  A.4 Differential evolution

The differential evolution (DE) method is a stochastic search algorithm for
minimizing real-valued global optimization problems based on the evolutionary
strategy [72, 36]. In DE, we initialize a population vector for each
generation $G$ as $\bm{x}_{i}^{G}\in\mathbb{R}^{n}\hskip
2.0pt(i=1,\dots,p,G=0,\dots,G_{\mathrm{max}})$, where $p$ is the number of
populations. The population vector is generated from the uniform distribution
to cover the entire parameter space. Next, we create a new population vector
with different ones as
$\tilde{\bm{x}}_{i}^{G+1}={\bm{x}}_{r_{1}}^{G}+F\left({\bm{x}}_{r_{2}}^{G}-{\bm{x}}_{r_{3}}^{G}\right)$,
where $r_{1},r_{2},$ and $r_{3}$ are indices of the population determined by
the “crossover” rule and $F>0$ is the amplification factor. The new trial
vector $\bm{u}_{i}^{G+1}$ can be obtained by mixing the target vector
$\bm{x}_{i}$ and the new population vector $\tilde{\bm{x}}_{i}^{G+1}$ with the
“crossover” rule [73]. If the new trial vector has a lower objective function
value than the target vector, the new trial vector supersedes the target
vector and is inherited by the next generation. The above process is called
the “selection”. The algorithm stops at $G=G_{\mathrm{max}}$ or the mean
objective function value over all populations becomes lower than the threshold
$\epsilon$. In our experiments, we adopt the “best/2/bin” strategy [72, 36]
and set $\epsilon=10^{-7}$. The number of populations $G_{\mathrm{max}}$
corresponds to the maximum number of iterations $t_{\mathrm{max}}$ in the main
text. The other hyperparameters are used in the default settings.

##  Appendix B Additional data for QA

In this Appendix, we show some additional data from the D-Wave and consider
the problem with $k=0.5$, $w_{0}=0.2$, and $h_{0}=1.0$. We compute the
observables as illustrated in Sec. III.2. The experimental setting of the
D-Wave is almost the same as those in Fig. 4.

Figure 13: Dependence of the observables on the annealing time for the different quantum annealers. Both axes are the same as those in Fig. 4.  Table 3: The fitting exponents of the observables for the different quantum annealers. The ”$E_{\mathrm{res}}$ (short)” represents the residual energy in the short annealing time region $1\leq t_{a}\leq 10$. The ”$E_{\mathrm{res}}$ (long)” is the residual energy in the long annealing time region $100\leq t_{a}\leq 500$. machine | 2000Q | Advantage4.1 | Advantage6.1 | Prototype1.1
---|---|---|---|---
$\rho$ | $0.359\pm 0.011$ | $0.337\pm 0.008$ | $0.301\pm 0.005$ | $0.239\pm 0.015$
$E_{\mathrm{res}}$ (short) | $0.553\pm 0.025$ | $0.383\pm 0.012$ | $0.338\pm 0.006$ | $0.319\pm 0.021$
$E_{\mathrm{res}}$ (long) | $0.201\pm 0.02$ | $0.341\pm 0.036$ | $-0.218\pm 0.012$ | $0.201\pm 0.031$

###  B.1 Comparison of D-Wave 2000Q and D-Wave Advantage

Figure 13 illustrates the dependence of observables on the annealing time for
a few different quantum annealers. In addition to the “DW_2000Q_6” device, we
utilize the “Advantage_system 4.1”, “Advantage_system 6.1”, and
“Advantage2_prototype 1.1” devices and compare their performance.

Regardless of the type of annealers, the polynomial decay of $\rho$ can be
observed in the short annealing time region $1(\mu s)\leq t_{a}\leq 10(\mu s)$
as shown in Fig. 13 (a). The fitting exponents are summarized in Table 3. The
largest exponent of $\rho$ is obtained from the D-Wave 2000Q. Except for the
D-Wave Advantage 6.1, the kink density gradually decreases as we increase the
annealing time. The kink density obtained by the D-Wave Advantage 6.1 shows
similar behavior to the D-Wave Advantage 4.1 in the region $1(\mu s)\leq
t_{a}\leq 10(\mu s)$ and starts to increase at $t_{a}=30$ $(\mu s)$. This
behavior may be due to decoherence. In the long annealing time region
$t_{a}\geq 100(\mu s)$, no device reaches the target $\rho=1/N$. In this
region, the D-Wave 2000Q takes similar values as the D-Wave Advantage2
Prototype 1.1 does.

Figure 13 (b) illustrates that the D-Wave 2000Q has the largest
$P_{\mathrm{const}}$ in the region $1(\mu s)\leq t_{a}\leq 100(\mu s)$. The
D-Wave Advantage2 Prototype 1.1 takes similar values of $P_{\mathrm{const}}$
to those obtained by the D-Wave 2000Q in the region $100(\mu s)\leq t_{a}\leq
500(\mu s)$. In the long annealing time region, most samples obtained by the
above two devices satisfy the single domain wall condition. Both devices yield
the larger value of $P_{\mathrm{const}}$ than the D-Wave Advantage 4.1 and
6.1.

Figure 13 (c) shows that the residual energy decays polynomially irrespective
of the types of the quantum annealer in the short annealing time region as
also seen in Fig. 13 (a). From Table 3, the largest exponent of
$E_{\mathrm{res}}$ can be obtained by the D-Wave 2000Q. In our experiments,
the D-Wave 2000Q gives the lowest value of $E_{\mathrm{res}}$ among all
devices. In the long annealing time region $t_{a}\geq 100(\mu s)$, we can see
a polynomial decay of $E_{\mathrm{res}}$. The fitting exponents obtained by
the D-Wave 2000Q and the D-Wave Advantage2 prototype 1.1 take similar values
to those in Table 3. The result indicates that qualitative behaviors of the
two devices show little differences. The D-Wave Advantage 4.1 has the largest
exponent among all devices, from which one may expect Advantage 4.1 to yield
smaller values of $E_{\mathrm{res}}$ than other models if the device would
operate with minimal noise in the longer annealing time region.

Figure 13 (d) demonstrates that the D-Wave 2000Q gives the largest value of
$P_{\mathrm{GS}}$. Although the D-Wave 2000Q and the D-Wave Advantage2
Prototype 1.1 give similar exponents in $E_{\mathrm{res}}$ in the long
annealing time region from Table 3, the performance gap appears in
$P_{\mathrm{GS}}$. This result implies that the D-Wave Advantage2 Prototype
1.1 may reach low energy states, but not necessarily the ground states.

The observable values obtained by different quantum annealers show different
behaviors. We generally conclude that the D-Wave 2000Q is suitable for
optimizing the Rastrigin function under the present problem setting.

Figure 14: Dependence of the observables attained by the D-Wave on the annealing time for different values of $\lambda$. Both axes are the same as those in Fig. 4.  Table 4: The fitting exponents of the observables for different values of $\lambda$. The fitting region is the same as used in Table 3. $\lambda$ | $0.1$ | $0.5$ | $1.0$
---|---|---|---
$\rho$ | $0.359\pm 0.011$ | $0.372\pm 0.011$ | $0.359\pm 0.011$
$E_{\mathrm{res}}$ (short) | $0.536\pm 0.026$ | $0.553\pm 0.024$ | $0.553\pm 0.025$
$E_{\mathrm{res}}$ (long) | $0.545\pm 0.073$ | $0.392\pm 0.019$ | $0.201\pm 0.02$

###  B.2 Dependence of the observables for different values of $\lambda$

Figure 14 shows the dependence of the observables on the annealing time for
different values of $\lambda$. All data are obtained from the D-Wave 2000Q. As
in the main text, we refer to the D-Wave 2000Q as the D-Wave. The data for
$\lambda=1.0$ is the same as used in Fig. 4. In Fig. 14 (a), the polynomial
decay of $\rho$ can be seen in the short annealing time region $1(\mu s)\leq
t_{a}\leq 10(\mu s)$ for different values of $\lambda$. The fitting exponents
are summarized in Table 4. In this region, the coefficient $\lambda$ hardly
affects the behavior of the system. On the other hand, in the long annealing
time region $t_{a}\geq 100(\mu s)$, the saturation seen in $\rho$ for
$\lambda=1.0$ does not occur in $\lambda=0.1$ and $0.5$. The small $\lambda$
weakens the effect of the problem term Eq. (10) and enhances the effect of the
constraint term Eq. (8).

Figure 14 (b) demonstrates the value of $\lambda$ hardly affects
$P_{\mathrm{const}}$ in the short annealing time region. Similarly to the kink
density in Fig. 14 (a), the effect of $\lambda$ can be seen in the long
annealing time region $t_{a}\geq 100(\mu s)$. Premature finite saturation of
$P_{\mathrm{const}}$ appears for $\lambda=1.0$. The number of samples that
satisfy the single domain wall condition increases by adjusting $\lambda$. In
the long annealing time region, the effect of quantum fluctuation becomes
weaker due to the decoherence. In this region, thermal fluctuations affect the
system. Since small values of $\lambda$ largely reduce the height of the
energy barrier and the lower energy barrier is favorable for thermal
fluctuations as shown in Fig. 9, the number of samples that satisfy the single
domain wall condition increases with decreasing $\lambda$ in this region.

Figure 14 (c) shows that the polynomial decay of $E_{\mathrm{res}}$ can be
observed in the short annealing time region $1(\mu s)\leq t_{a}\leq 10(\mu
s)$. From Table 4, the values of $\lambda$ hardly affects the exponents of
$E_{\mathrm{res}}$ as in the kink density. In the long annealing time region
$t_{a}\geq 100(\mu s)$, the residual energy decreases polynomially and the
fitting exponents are written in Table 4 (b). A smaller $\lambda$ yields a
larger exponent. This result implies that a smaller $\lambda$ leads to lower
energy states and reduces the hardness of optimization due to the higher
energy barrier.

Figure 14 (d) illustrates that larger values of $\lambda$ give larger
$P_{\mathrm{GS}}$. This result is in contrast to the behavior of
$E_{\mathrm{res}}$ in Fig. 14 (c). Smaller values of $\lambda$ decrease the
values of the magnetic field physically encoded into the D-Wave. Since the bit
width to represent the problem Hamiltonian is limited in the D-Wave, the
actual implemented value in the D-Wave may be different. Due to this
discretization error, the actual energy landscape slightly changes. Because
the success probability is more affected by variation of the energy landscape
than the residual energy, the success probability for the smaller $\lambda$ is
degraded.

##  Appendix C Comparison of the data of the D-Wave 2000Q and thermal
equilibrium state

The D-Wave is a device to reproduce the system behavior under the time
dependent Hamiltonian $H(s)$ in Eq. (2). However, as the hardware is an open
system, its dynamics is affected by thermal noise. To describe such dynamics,
the hypothesis of quasi-static time evolution and freeze-out phenomenon has
been proposed [74]. In this hypothesis, during the annealing process, the
density matrix of the system at a given time is approximately equal to the
density matrix of the thermal equilibrium state of the Hamiltonian at that
time. The dynamics freezes at a time denoted as $s^{*}$, when the relaxation
time to equilibrium becomes much larger than the time scale of annealing, and
the density matrix of the final result becomes
$\exp(-H(s^{*})/T_{\mathrm{phys}})$, where $T_{\mathrm{phys}}$ denotes the
physical device temperature that is about 13.5 mK. If the strength of the
transverse field at the freeze-out point, $A(s^{*})$, is much smaller than the
classical part of the Hamiltonian, $B(s^{*})$, then the density matrix of the
system can be approximated as

$\displaystyle\rho_{\mathrm{sys}}(s^{*})\approx\frac{e^{-(B(s^{*})/2)H_{0}/T_{\mathrm{phys}}}}{\mathrm{Tr}[e^{-(B(s^{*})/2)H_{0}/T_{\mathrm{phys}}}]}.$
(29)

This can be interpreted as the density matrix at the thermal equilibrium state
of $H_{0}$ with a temperature of $2T_{\mathrm{phys}}/B(s^{*})$. By determining
this rescaled temperature from the output samples of the D-Wave, it is
possible to estimate the freeze-out point $s^{*}$. It should be noted that if
the obtained effective temperature value is too high, the estimation of
$s^{*}$ may become inaccurate and it would be appropriate to estimate the
temperature using $H(s^{*})$ instead of $H_{0}$.

We now discuss the consistency between the experimental results of the D-Wave
in this study and the hypothesis of quasi-static time evolution and freeze-
out. The effective temperature $T_{\mathrm{eff}}$ is defined as the value
obtained by solving the equation $E(T)=\langle H_{0}\rangle$ for $T$, where
$E(T)$ is the internal energy of the thermal equilibrium state of $H_{0}$ at
temperature $T$ and $\langle H_{0}\rangle$ is the energy expectation value of
the data of the D-Wave. The annealing time dependence of the effective
temperature from the experimental results of the D-Wave is shown in Fig. 15.
For $h_{0}=0.2$, the effective temperature decreases toward
$T_{\mathrm{phys}}$ while for $h_{0}=1.0$ and $3.0$, the effective temperature
saturates at a constant value several times higher than
$2T_{\mathrm{phys}}/B(1)$.

Figure 15: Effective temperature obtained from the data of the D-Wave as a
function of annealing time. The dashed horizontal line represents the physical
temperature of the D-Wave. Figure 16: Relation between the inverse effective
temperature and freeze-out point. The lines represent the annealing schedule
functions used to determine $s^{*}$ from $T_{\mathrm{eff}}$.

Regarding $T_{\mathrm{eff}}$ as $2T_{\mathrm{phys}}/B(s^{*})$, the freeze-out
point $s^{*}$ is estimated as the value obtained by solving the equation
$2T_{\mathrm{phys}}/B(s)=T_{\mathrm{eff}}$ for $s$. The relation between
$s^{*}$ and $1/T_{\mathrm{eff}}$ is shown in Fig. 16 and the annealing
schedule function $B(s)$ used to determine $s^{*}$ and the schedule function
of the transverse field $A(s)$ are also shown. For $h_{0}=0.2$, at
$t_{a}=500(\mu s)$ the freeze-out point $s^{*}$ reaches approximately 0.6,
where $B(s^{*})\gg A(s^{*})$, indicating that the classical term is dominant
in $H(s^{*})$. For $h_{0}=1.0$ and $h_{0}=3.0$, $s^{*}$ reaches approximately
0.35 and 0.4, respectively. For these $s^{*}$ values, $B(s^{*})\approx
A(s^{*})$ and $B(s^{*})<A(s^{*})$, respectively, which may lead to the
deviation of the data of the D-Wave from the thermal equilibrium of $H_{0}$.

Figure 17: Comparison between data of the D-Wave (shown by symbols with lines)
and the thermal equilibrium of $H_{0}$ at the corresponding effective
temperature (shown by filled triangles). (a) Constraint satisfaction
probability and (b) ground state probability as functions of annealing time.
Figure 18: Comparison between data of coherent QA obtained by TEBD (shown by
symbols with lines) and the thermal equilibrium of $H_{0}$ at the
corresponding effective temperature (shown by filled triangles). (a)
Constraint satisfaction probability and (b) ground state probability as
functions of annealing time. The annealing time $t_{a}$ is in the unit of $J$.

To examine how closely the D-Wave data matches a thermal equilibrium value, we
compare the data of the D-Wave and the observables of the thermal equilibrium
of $H_{0}$ at the obtained effective temperature as shown in Fig. 17. For
$h_{0}=0.2$ and for $t_{a}>10$, where $s^{*}>0.5$, the ground state
probability $P_{\mathrm{GS}}$ and the constraint satisfaction probability
$P_{\mathrm{const}}$ are close to the value of the thermal equilibrium state,
implying that the hypothesis of quasi-static time evolution is valid. On the
other hand, in the region where $s^{*}<0.5$ for all $h_{0}$ values,
$P_{\mathrm{GS}}$ of the D-Wave is larger than that of the thermal
equilibrium. Furthermore, $P_{\mathrm{const}}$ of the D-Wave is smaller than
that of thermal equilibrium. This result suggests that in the case of
$s^{*}<0.5$, the hypothesis of quasi-static time evolution may not hold, or
that the final state of the annealing is affected by the transverse field.

A similar analysis of the effective temperature is performed for the results
of TEBD, as shown in Fig. 18. As in the case of the D-Wave, $P_{\mathrm{GS}}$
is larger, and $P_{\mathrm{const}}$ is smaller than the value of the
corresponding thermal equilibrium state. This result suggests that it may be
generally observed that the single kink state probability is larger, and the
ground state probability is smaller, in non-equilibrium states compared to the
thermal equilibrium states at the corresponding effective temperature.

In summary, the annealing time dependence of the effective temperature varies
significantly depending on $h_{0}$. For $h_{0}=0.2$, the effective temperature
decreases toward the hardware temperature while it saturates at higher
temperatures for $h_{0}=1.0$ and $h_{0}=3.0$. Moreover, the value of effective
temperature changes non-monotonically with $h_{0}$. The mechanism of this
nontrivial dependence of the effective temperature on $h_{0}$ remains to be
clarified. In the region where the freeze-out point is small and not in
thermal equilibrium of $H_{0}$, the behaviors of $P_{\mathrm{GS}}$ and
$P_{\mathrm{const}}$ exhibit similar characteristics to coherent QA, although
the underlying mechanism remains unclear. For a deeper understanding of the
results, analysis taking the effects of thermal noise into account is
necessary.

## References

  * Kadowaki and Nishimori [1998] T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355 (1998).
  * Farhi  _et al._ [2001] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science 292, 472 (2001).
  * Santoro  _et al._ [2002] G. E. Santoro, R. Martoňák, E. Tosatti, and R. Car, Science 295, 2427 (2002).
  * Santoro and Tosatti [2006] G. E. Santoro and E. Tosatti, J. Phys. A: Math. Gen. 39, R393 (2006).
  * Das and Chakrabarti [2008] A. Das and B. K. Chakrabarti, Rev. Mod. Phys. 80, 1061 (2008).
  * Morita and Nishimori [2008] S. Morita and H. Nishimori, J. Math. Phys. 49, 125210 (2008).
  * Tanaka  _et al._ [2017] S. Tanaka, R. Tamura, and B. K. Chakrabarti,  _Quantum Spin Glasses, Annealing and Computation_ (Cambridge University Press, 2017).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Rev. Mod. Phys. 90, 015002 (2018).
  * Hauke  _et al._ [2020] P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, and W. D. Oliver, Rep. Prog. Phys. 83, 054401 (2020).
  * Grant and Humble [2020] E. K. Grant and T. S. Humble, in  _Oxford Research Encyclopedia of Physics_ (Oxford University Press, 2020) pp. 1–26.
  * Crosson and Lidar [2021] E. J. Crosson and D. A. Lidar, Nat. Rev. Phys. 3, 466 (2021).
  * Johnson  _et al._ [2010] M. W. Johnson, P. Bunyk, F. Maibaum, E. Tolkacheva, A. J. Berkley, E. M. Chapple, R. Harris, J. Johansson, T. Lanting, I. Perminov, E. Ladizinsky, T. Oh, and G. Rose, Supercond. Sci. Technol. 23, 065004 (2010).
  * Berkley  _et al._ [2010] A. J. Berkley, M. W. Johnson, P. Bunyk, R. Harris, J. Johansson, T. Lanting, E. Ladizinsky, E. Tolkacheva, M. H. S. Amin, and G. Rose, Supercond. Sci. Technol. 23, 105014 (2010).
  * Harris  _et al._ [2010] R. Harris, M. W. Johnson, T. Lanting, A. J. Berkley, J. Johansson, P. Bunyk, E. Tolkacheva, E. Ladizinsky, N. Ladizinsky, T. Oh, F. Cioata, I. Perminov, P. Spear, C. Enderud, C. Rich, S. Uchaikin, M. C. Thom, E. M. Chapple, J. Wang, B. Wilson, M. H. S. Amin, N. Dickson, K. Karimi, B. Macready, C. J. S. Truncik, and G. Rose, Phys. Rev. B 82, 024511 (2010).
  * Dickson  _et al._ [2013] N. G. Dickson, M. W. Johnson, M. H. Amin, R. Harris, F. Altomare, A. J. Berkley, P. Bunyk, J. Cai, E. M. Chapple, P. Chavez, F. Cioata, T. Cirip, P. deBuen, M. Drew-Brook, C. Enderud, S. Gildert, F. Hamze, J. P. Hilton, E. Hoskinson, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Lanting, T. Mahon, R. Neufeld, T. Oh, I. Perminov, C. Petroff, A. Przybysz, C. Rich, P. Spear, A. Tcaciuc, M. C. Thom, E. Tolkacheva, S. Uchaikin, J. Wang, A. B. Wilson, Z. Merali, and G. Rose, Nat. Commun. 4, 1903 (2013).
  * Lanting  _et al._ [2014] T. Lanting, A. J. Przybysz, A. Y. Smirnov, F. M. Spedalieri, M. H. Amin, A. J. Berkley, R. Harris, F. Altomare, S. Boixo, P. Bunyk, N. Dickson, C. Enderud, J. P. Hilton, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, R. Neufeld, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, S. Uchaikin, A. B. Wilson, and G. Rose, Phys. Rev. X 4, 021041 (2014).
  * McGeoch and Farré [2022] C. McGeoch and P. Farré,  _Advantage Processor Overview: Technical Report_, Tech. Rep. (2022).
  * Finnila  _et al._ [1994] A. Finnila, M. Gomez, C. Sebenik, C. Stenson, and J. Doll, Chem. Phys. Lett. 219, 343 (1994).
  * Stella  _et al._ [2005] L. Stella, G. E. Santoro, and E. Tosatti, Phys. Rev. B 72, 014303 (2005).
  * Shinomoto and Kabashima [1991] S. Shinomoto and Y. Kabashima, J. Phys. A: Math. Gen. 24, L141 (1991).
  * Inack and Pilati [2015] E. M. Inack and S. Pilati, Phys. Rev. E 92, 053304 (2015).
  * Koh and Nishimori [2022] Y. W. Koh and H. Nishimori, Phys. Rev. A 105, 062435 (2022).
  * Dieterich and Hartke [2012] J. M. Dieterich and B. Hartke, arXiv:1207.4318 [cs.NE] (2012).
  * Chancellor [2019] N. Chancellor, Quantum Sci. Technol. 4, 045004 (2019).
  * Abel  _et al._ [2021] S. Abel, N. Chancellor, and M. Spannowsky, Phys. Rev. D 103, 016008 (2021).
  * Abel and Spannowsky [2021] S. Abel and M. Spannowsky, PRX Quantum 2, 010349 (2021).
  * Abel  _et al._ [2022] S. Abel, A. Blance, and M. Spannowsky, Phys. Rev. A 106, 042607 (2022).
  * Perdomo-Ortiz  _et al._ [2011] A. Perdomo-Ortiz, S. E. Venegas-Andraca, and A. Aspuru-Guzik, Quantum Inf. Process 10, 33 (2011).
  * Ohkuwa  _et al._ [2018] M. Ohkuwa, H. Nishimori, and D. A. Lidar, Phys. Rev. A 98, 022314 (2018).
  * Yamashiro  _et al._ [2019] Y. Yamashiro, M. Ohkuwa, H. Nishimori, and D. A. Lidar, Phys. Rev. A 100, 052321 (2019), 1906.10889 .
  * Arai  _et al._ [2021] S. Arai, M. Ohzeki, and K. Tanaka, Phys. Rev. Research 3, 033006 (2021).
  * Note [1] In Ref. [27], RA is referred to as QA.
  * Nelder and Mead [1965] J. A. Nelder and R. Mead, Comput. J. 7, 308 (1965).
  * Hestenes and Stiefel [1952] M. R. Hestenes and E. Stiefel, J. Res. Natl. Bur. Stand. 49, 409 (1952).
  * Leary [2000] R. H. Leary, J. Global Optim. 18, 367 (2000).
  * Storn and Price [1997] R. Storn and K. Price, J. Global Optim. 11, 341 (1997).
  * Vidal [2004] G. Vidal, Phys. Rev. Lett. 93, 040502 (2004).
  * White and Feiguin [2004] S. R. White and A. E. Feiguin, Phys. Rev. Lett. 93, 076401 (2004).
  * Kirkpatrick  _et al._ [1983] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Science 220, 671 (1983).
  * Suzuki [1976] M. Suzuki, Prog. Theor. Phys. 56, 1454 (1976).
  * Shin  _et al._ [2014] S. W. Shin, G. Smith, J. A. Smolin, and U. Vazirani, arXiv:1401.7087 [quant-ph] (2014).
  * [42] D-wave ocean software documentation, https://docs.ocean.dwavesys.com/en/stable/index.html, accessed: 2023-03-08.
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, A. Vijaykumar, A. P. Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C. N. Woods, C. Fulton, C. Masson, C. Häggström, C. Fitzgerald, D. A. Nicholson, D. R. Hagen, D. V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G. A. Price, G.-L. Ingold, G. E. Allen, G. R. Lee, H. Audren, I. Probst, J. P. Dietrich, J. Silterra, J. T. Webber, J. Slavič, J. Nothman, J. Buchner, J. Kulick, J. L. Schönberger, J. de Miranda Cardoso, J. Reimer, J. Harrington, J. L. C. Rodríguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. Thoma, M. Newville, M. Kümmerer, M. Bolingbroke, M. Tartre, M. Pak, N. J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P. A. Brodtkorb, P. Lee, R. T. McGibbon, R. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T. J. Pingel, T. P. Robitaille, T. Spura, T. R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upadhyay, Y. O. Halchenko, Y. Vázquez-Baeza, and S. . . Contributors, Nat. Methods 17, 261 (2020).
  * Note [2] Extensive CPU benchmarks are performed in Geekbench [75]. The single-core score of the Apple M1 Max chip in Geekbench ver. 6 is higher than those of the Intel Xeon E5-1650 used in an old benchmark [76] and the Intel Xeon E5-2695 v4 in a recent benchmark [77].
  * Note [3] In the original formulation of the Rastrigin function, the search space is defined in $x\in[-5.12,5.12]$ [23]. Taking into account the computational precision of QA, we reduce the search space.
  * Note [4] By substituting $x_{j}=0$ into the right-hand of Eq. (5), the location of the zero magnetic field can be computed as $j^{\prime}=106$ with the present numerical values, $N=211$ and $x_{\rm max}=-x_{\rm min}=3$. For $\hat{\sigma}_{106}^{z}=-1$, the ground state $x=0$ can be recovered as $x_{106}=0$. For $\hat{\sigma}_{106}^{z}=+1$, the domain wall is located between $j=105$ and $j=106$, and $x_{105}=-\Delta x(\neq 0)$ results.
  * Nishimura  _et al._ [2020] K. Nishimura, H. Nishimori, and H. G. Katzgraber, Phys. Rev. A 102, 042403 (2020).
  * King  _et al._ [2022] A. D. King, S. Suzuki, J. Raymond, A. Zucca, T. Lanting, F. Altomare, A. J. Berkley, S. Ejtemaee, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, J. D. Whittaker, J. Yao, R. Harris, D. A. Lidar, H. Nishimori, and M. H. Amin, Nat. Phys. 18, 1324 (2022).
  * Note [5] The NM and CGD methods get stuck in local minima before the number of iterations of the algorithms reaches the maximum value. Therefore, an increase in the maximum number of iterations does not affect the execution time, and their rightmost points have the same values of the execution time.
  * Weinberg  _et al._ [2020] P. Weinberg, M. Tylutki, J. M. Rönkkö, J. Westerholm, J. A. Åström, P. Manninen, P. Törmä, and A. W. Sandvik, Phys. Rev. Lett. 124, 090502 (2020).
  * Hamasaki and Nishimori [2004] T. Hamasaki and H. Nishimori, J. Phys. Soc. Jpn. 73, 1490 (2004).
  * Suzuki and Okada [2007] S. Suzuki and M. Okada, Interdisciplinary Information Sciences 13, 49 (2007).
  * Schollwöck [2011] U. Schollwöck, Ann. Phys. 326, 96 (2011), january 2011 Special Issue.
  * Kibble [1976] T. W. B. Kibble, J. Phys. A: Math. Gen. 9, 1387 (1976).
  * Zurek [1985] W. H. Zurek, Nature (London) 317, 505 (1985).
  * Dziarmaga [2005] J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005).
  * Suzuki and Okada [2005] S. Suzuki and M. Okada, J. Phys. Soc. Jpn. 74, 1649 (2005).
  * Zaborniak and de Sousa [2021] T. Zaborniak and R. de Sousa, IEEE Trans. Quantum Eng. 2, 1 (2021).
  * Chancellor  _et al._ [2022] N. Chancellor, P. J. D. Crowley, T. Durić, W. Vinci, M. H. Amin, A. G. Green, P. A. Warburton, and G. Aeppli, npj Quantum Inf. 8, 73 (2022).
  * Könz  _et al._ [2021] M. S. Könz, W. Lechner, H. G. Katzgraber, and M. Troyer, PRX Quantum 2, 040322 (2021).
  * Bando and Nishimori [2021] Y. Bando and H. Nishimori, Phys. Rev. A 104, 022607 (2021).
  * Jiang  _et al._ [2017] Z. Jiang, V. N. Smelyanskiy, S. V. Isakov, S. Boixo, G. Mazzola, M. Troyer, and H. Neven, Phys. Rev. A 95, 012322 (2017).
  * Albash and Marshall [2021] T. Albash and J. Marshall, Phys. Rev. Appl. 15, 014029 (2021).
  * Bando  _et al._ [2022] Y. Bando, K.-W. Yip, H. Chen, D. A. Lidar, and H. Nishimori, Phys. Rev. Appl. 17, 054033 (2022).
  * Note [6] Notice that the D-Wave becomes stuck at a value larger than $\rho=1/N$ as shown in Fig. 4 (a) though this behavior is not clearly seen in the scale of Fig. 12 (a).
  * King  _et al._ [2023] A. D. King, J. Raymond, T. Lanting, R. Harris, A. Zucca, F. Altomare, A. J. Berkley, K. Boothby, S. Ejtemaee, C. Enderud, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, R. Molavi, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, N. Tsai, M. Volkmann, J. D. Whittaker, J. Yao, A. W. Sandvik, and M. H. Amin, Nature 617, 61 (2023).
  * Gao and Han [2012] F. Gao and L. Han, Comput. Optim. Applic. 51, 259 (2012).
  * Nocedal and Wright [2006] J. Nocedal and S. J. Wright,  _Numerical Optimization_ , 2nd ed. (Springer, New York, NY, USA, 2006).
  * Fletcher and Reeves [1964] R. Fletcher and C. M. Reeves, Comput. J. 7, 149 (1964).
  * Iwamatsu and Okabe [2004] M. Iwamatsu and Y. Okabe, Chem. Phys. Lett. 399, 396 (2004).
  * Luo  _et al._ [2012] Z. Luo, B. Olson, I. Hashmi, K. Molloy, and A. Shehu, Adv. Artif. Intell 2012, 674832 (2012).
  * Storn [1996] R. Storn, in  _Proc. North Am. Fuzzy. Inf. Process_ (1996) pp. 519–523.
  * Das and Suganthan [2011] S. Das and P. N. Suganthan, IEEE Trans. Evol. Comput. 15, 4 (2011).
  * Amin [2015] M. H. Amin, Phys. Rev. A 92, 052323 (2015).
  * [75] Geekbench browser, https://browser.geekbench.com/, accessed: 2023-03-08.
  * Denchev  _et al._ [2016] V. S. Denchev, S. Boixo, S. V. Isakov, N. Ding, R. Babbush, V. Smelyanskiy, J. Martinis, and H. Neven, Phys. Rev. X 6, 031015 (2016).
  * Tasseff  _et al._ [2022] B. Tasseff, T. Albash, Z. Morrell, M. Vuffray, A. Y. Lokhov, S. Misra, and C. Coffrin, arXiv:2210.04291 [math.OC] (2022).