Quantum and classical annealing in a continuous space with multiple local minima

Let us consider the following potential inspired by Shinomoto and Kabashima Shinomoto and Kabashima (1991),

The first is a harmonic term forming the envelope of the overall potential, and it sets a unique global minimum at $x=0$. Without loss of generality, we shall let the spring constant $k$ be unity from now on. The cosine term introduces local minima onto the energy landscape, and the structure of these minima can be tuned using two parameters $h_{0}$ and $w_{0}$. As an example, Fig. 1 shows the graph of $V_{\rm SK}(x)$ when $h_{0}=w_{0}=0.2$. The parameter $h_{0}$ corresponds approximately to the height of the energy barrier between two adjacent minima, and $w_{0}$ to the horizontal distance between them.

Our proposed potential $V_{\rm SK}(x)$ tries to capture the essential physics underlying Shinomoto and Kabashima’s idea. Nevertheless, there are some differences because in their original formulation, it is assumed that the local minima are all uniform such that the barrier height and distance between neighboring minima remain constant ad infinitum along the $x$-axis. On the other hand, in Eq. (1), due to the presence of the quadratic term, the local minima of $V_{\rm SK}(x)$ are smoothed out far away from the origin and so the total number of local minima is in fact finite. This is illustrated by Fig. 1, where the inset shows the global view of the same graph. Notwithstanding this difference, Eq. (1) is expected to lead to essentially the same asymptotic behavior of the system because the probability of the system being far from the origin is very small in the long-time limit. We adopt Eq. (1) because of its simpler numerical implementation and yet with essentially the same asymptotic properties to be expected.

Let us define $N_{\rm min}$ as the number of local minima in the positive $x$-direction. It is seen that $N_{\rm min}=15$ when the parameters are $h_{0}=w_{0}=0.2$. Notice that although the local minima are relatively uniform in the vicinity of the origin, they become shallower further away, so our studies should be focused on the former region in order to stay faithful to Shinomoto and Kabashima’s proposal.

We calculated $N_{\rm min}$ numerically as a function of the parameters $h_{0}$ and $w_{0}$, and the results are shown in Fig. 2 in the form of a phase diagram. Individual regions of $N_{\rm min}$ from 0 up to 7 are shown. In the region $N_{\rm min}\geq 8$, the boundaries between successive regions become increasingly closely spaced, with $N_{\rm min}$ approaching infinity as one approaches the vertical axis. The graph of Fig. 1 is indicated by the circle, located inside the region $N_{\rm min}\geq 8$. Roughly speaking, one can interpret the phase diagram of $N_{\rm min}$ as a map of the difficulty of optimizing the potential $V_{\rm SK}(x)$. Regions where $N_{\rm min}$ is large are, intuitively, more difficult to optimize than regions where it is smaller. Within a region of constant $N_{\rm min}$, it is generally more difficult to optimize if either (or both) of the parameters $h_{0}$ and $w_{0}$ increases.

It is not feasible to perform a thorough computational study of the entire phase space. Instead, we shall focus on just one point, $(w_{0},h_{0})=(0.2,0.2)$, which has been introduced in Figs. 1 and 2. Our choice is due to the need to balance two competing factors. On the one hand, the energy landscape of the potential should be difficult enough to be non-trivial for our annealing algorithms to optimize. On the other hand, we are limited by the computational resources needed to solve the time-dependent Schrödinger equation. For instance, the top left corner of the phase space where $N_{\rm min}$ and $h_{0}$ are large would pose a good challenge to both quantum and simulated annealing; however, a very large number of grid points is required to accurately represent the wavefunction in this region of phase space, making it prohibitively expensive to compute energy eigenvalues and perform long annealing simulations. The energy landscape at the phase point $(0.2,0.2)$ is complex and non-trivial, but at the same time also computationally feasible to simulate, allowing us to focus on the physics of annealing rather than on the optimization of a potential with true technical difficulties. Hence, in the rest of the paper, unless otherwise stated, it is implicit that all numerical calculations are performed with the parameters $h_{0}=w_{0}=0.2$.

Quantum annealing of the potential $V_{\rm SK}(x)$ is performed by solving the time-dependent Schrödinger equation

where $\psi(x,t)$ is the wavefunction, and the time dependence of mass $m(t)$ serves as the annealing schedule. Conceptually speaking, one starts annealing from zero mass where the kinetic energy is dominating, and slowly increase it to infinity such that the potential $V_{\rm SK}(x)$ becomes the dominating term. If the change of $m(t)$ with time is infinitesimally slow, then according to the adiabatic theorem Albash and Lidar (2018), a wavefunction initially at the ground state will evolve into the global optimized solution of $V_{\rm SK}(x)$ upon the completion of annealing.

Let us restrict the time evolution under Eq. (6) to within the time interval $0\leq t\leq T$ where $T$ is the total annealing time. Consider the following linear form for $m(t)$

where $m_{i}$ and $m_{f}$ are the masses at the initial $(t=0)$ and final $(t=T)$ times, respectively. When $m_{i}$ and $m_{f}$ are fixed, $T$ determines the annealing speed, with a large $T$ giving slow annealing. The superscript (1) denotes ‘first-order’, and is to differentiate Eq. (7) from higher-order, nonlinear schedules which we shall introduced later. As suggested by Morita Morita (2007), only the initial and final time derivatives of annealing schedule determine the final excitation probability for sufficiently large annealing time. It is thus expected that an explicit form of the annealing function is unimportant, i.e. it is irrelevant whether we change the mass as in Eq. (7) or change the coefficient $\Gamma(t)\equiv\hbar/2m$ of the kinetic energy term as a whole as was done in Ref. Stella et al. (2005), as long as the initial and final time derivatives are shared among different functions. The linear schedule of Eq. (7) represents the case with finite derivatives at both ends of annealing.

One operational constraint faced when quantum annealing is implemented numerically via the solution of Eq. (6) concerns the choice of the mass parameters $m_{i}$ and $m_{f}$. In principle, one should be able to use arbitrary small $m_{i}$ and arbitrary large $m_{f}$ in the schedule Eq. (7). For instance, one would like to start from the high energy state shown in Fig. 3(a) and anneal towards the low energy state shown in Fig. 4(c). In practice, however, it is not feasible to solve Eq. (6) numerically with such a large energy difference because the number of grid points needed to represent the wavefunction is very large. In particular, one sees from Figs. 3(a) and 4(c) that the wavefunction evolves from a delocalized to a localized form, and so the grid must be wide enough to represent the initial wavefunction, while at the same time its points must also be dense enough to resolve the final wavefunction. Hence, in practice, the choice of $m_{i}$ and $m_{f}$ is partly constrained by the computational resources, both in terms of memory and processing speed, one has at one’s disposal. In all our quantum annealing simulations, we restricted ourselves to a plane wave grid with 2048 grid points, and then choose $m_{i}$ and $m_{f}$ that are computationally feasible with this grid size.

Despite the above issue, our ultimate goal is still to attain some understanding of annealing when the initial mass is small and the final mass is large. Consider the masses $m_{a}$ and $m_{d}$ listed in Table 1, where the difference $m_{d}-m_{a}$ spans six orders of magnitude. The ground state at $m_{a}$ was discussed in Fig. 3(a). To get an intuitive sense of the ground state at $m_{d}$, note that $m_{d}$ is larger than $m_{c}$, whose ground state was discussed in Fig. 4(c). Our objective is to study quantum annealing from the high energy, delocalized state at $m_{a}$ to the low energy, localized state at $m_{d}$.

To traverse this large mass difference, we perform annealing in stages. Referring to Table 1, let us define the annealing from $m_{a}$ to $m_{b}$ as Stage 1, that from $m_{b}$ to $m_{c}$ as Stage 2, and that from $m_{c}$ to $m_{d}$ as Stage 3. The properties of the ground states at $m_{a}$, $m_{b}$, and $m_{c}$ have been discussed in Sec. III. Apart from keeping the computations manageable, the stages are also set up based on energy considerations. In Stage 1, the energy decreases from a high level to half the barrier height of $V_{\rm SK}(x)$, as indicated by the blue lines in Figs. 3(a) and 4(a). In Stage 2, the energy continues its descent to the potential level of the two lowest local minima, as seen in Fig. 4(c). Finally, in Stage 3 the energy descends to close to the potential of the global minimum.

A graphical perspective of the annealing stages is also given in Figs. 5(a) and 6 where they are superposed upon the $E_{0}$ and $\Delta$ curves. In these figures, the vertical lines indicate the initial and final masses of each of the stages.

In numerical calculations, the grid in each stage is determined by making sure that it is wide enough to accommodate the ground state wavefunction at $m_{i}$ and its grid points dense enough to resolve the ground state wavefunction at $m_{f}$. This is done quantitatively by making sure that the eigenenergy $E_{0}$ is converged with respect to the number of grid points at both $m_{i}$ and $m_{f}$ ¹¹1More concretely, we first determine the limit $x_{\rm max}$ such that the ground state wavefunction at $m_{i}$ lies comfortably within the interval $[-x_{\rm max},x_{\rm max}]$. We then compute the energy eigenvalue $E_{0}$ using a grid spanning the interval with 2048 grid points. The computed $E_{0}$ must remain unchanged when the number of grid points is decreased to 1024 and increased to 4096. With the converged grid at $m_{i}$, we repeat the procedure at $m_{f}$ to ensure that the same grid must guarantee the convergence of $E_{0}$ at the final mass as well. . A grid that passed the convergence test at both ends of a stage is used for quantum annealing of that stage.

Quantum annealing of a stage is performed by inserting the particular $m_{i}$ and $m_{f}$ into the schedule Eq. (7), which is then substituted into Eq. (6). The resulting Schrödinger equation is solved using the WavePacket program package Schmidt and Lorenz (2017); Schmidt and Hartmann (2018); Lorenz and Schmidt (2020). Atomic units are used, and $m(t)$ is expressed in units of the electron rest mass. Time propagation is performed using the Runge-Kutta algorithm with a time step of $dt=0.1$ ²²2Convergence of trajectories with respect to $dt$ is checked by halving it (to 0.05) and then seeing that recalculated results remain invariant.. The ground state wavefunction at $m_{i}$ is used as the initial wavefunction at the start of annealing.

To assess the performance of annealing, we consider the residual energy, defined as

The first term is the expectation value of the Hamiltonian [given in Eq. (6)] taken with respect to the wavefunction obtained when annealing finishes at $t=T$. The second term is the ground state energy at the final mass, and is the target energy which annealing seeks to attain. In the limit $T\rightarrow\infty$, the time evolution becomes adiabatic and the residual energy approaches zero. The rate at which $R_{\rm q.a.}(T)$ decays with the total annealing time $T$ is a measure of how quickly the target energy is being attained by the annealing process.

Figure 7 shows $R_{\rm q.a.}(T)$ as a function of $T$ obtained under the linear schedule. The three curves show the results for the three stages discussed earlier. The straight lines (black) are proportional to $T^{-2}$ and indicate the asymptotic behavior of the curves when $T$ is large ³³3Due to the periodic oscillatory behavior of the data points, instead of peforming a least square fit, we have chosen instead to manually adjust the vertical displacement of the line $T^{-2}$ to fit those data points at the crests of the oscillatory curve. We find that this reflects more accurately the observed gradients of the curves.. We also performed some additional simulations at other $h_{0}$-$w_{0}$ phase points and found that the decay rate $T^{-2}$ holds quite generally in phase space. This decay rate for the linear schedule is in agreement with previous studies on spin systems Suzuki and Okada (2005); Morita (2007).

The curve of Stage 2 (blue squares) is slightly different from the other two stages, first undergoing a transient phase before crossing over into the asymptotic $T^{-2}$ regime. Similar two-phased behavior in the residual energy has also been reported in a two-level system, with the earlier phase being attributed to diabatic transitions and explained by the Landau-Zener theorem Morita (2007). The transient phase here is probably also due to transitions into excited states. Figure 8 shows an example of the final wavefunction that is obtained in Stage 2, with the corresponding residual energy indicated in Fig. 7 (arrowed $T=3250$). For comparison, the target ground state is shown in Fig. 4(c). One sees that the annealed wavefunction exhibits two spurious peaks at the neighboring local minima, which are absent in the target wavefunction. Furthermore, we found that the first excited level in the ‘flat gap’ region (c.f. Fig. 6) is doubly-degenerate, spanned by two compact gaussians each centered at one of the two local minima. Hence, the spurious peaks in the annealed wavefunction in Fig. 8 are most likely caused by diabatic transition into the first excited level when the annealing traverses the ‘flat gap’ region.

It is interesting to see that in quantum annealing, parts of the wavefunction can become trapped near local minima if the annealing time is not long enough. Normally, this may be what one would expect more in simulated annealing. This similarity between quantum and simulated annealing is inspiring from a physics point of view, especially since two annealing protocols obey very different rules. The mechansim underlying the trapping in simulated annealing is fairly well-understood. On the other hand, we might try to explain the trapping observed for quantum annealing using Landau-Zener transition. In the case of $V_{\rm SK}(x)$, we have seen that the system encounters a region of ‘flat gap’ during the course of annealing. This may pose a problem for the Landau-Zener theory as there is no well-defined avoided crossing, which is one of the basic tenets of the Landau-Zener theory. Hence, although the physical process underlying the trapping observed here might be a diabatic transition, the detailed mechanism could differ from the traditional Landau-Zener one.

We now discuss the time evolution of the wavefunction during the course of annealing. Figure 9(a) shows a snapshot of the probability density at time $t=44$ in Stage 1, with $T=560$. The single-modal form is quite representative of the shape of the wavefunction in Stage 1, so we can express its temporal behavior more succinctly using the width $\langle\psi(t)|x^{2}|\psi(t)\rangle$, which is illustrated in Fig. 9(a). Figure 10(a) shows the time evolution of the width in Stage 1 for several values of $T$. Overall, one sees that all three graphs oscillate, with the amplitude of oscillation decaying with the passage of time. Among the graphs themselves, as $T$ increases, the oscillation becomes suppressed to yield a smoother decay.

Although informative, looking at the width alone gives us rather limited understanding of the behavior of the wavefunction. More insight can be gained by looking at the so-called probability current, defined as

where $\psi^{*}(x,t)$ is the complex conjugate of $\psi(x,t)$, $i=\sqrt{-1}$, and $m(t)$ is again the time-dependent mass. The current $j(x,t)$ is a real-valued quantity, and indicates the direction in which the probability is flowing at a point in space ⁴⁴4For the plane wave $e^{ikx}$, one has $j(x,t)=\frac{\hbar k}{m}=\frac{p}{m}=v$, which is simply the velocity of the particle. Hence, the probability current can be interpreted as a generalization of the velocity field, indicating both the speed and direction of motion of the particle. . Figure 9(b) shows the probability current corresponding to the probability density shown in panel (a). It is seen that the current is positive (negative) when $x>0$ ($x<0$), meaning that probability is flowing outwards away from the global minimum at the origin. With respect to the purpose of optimization, such an instantaneous behavior of the wavefunction should be considered unfavorable because the particle is diffusing away from the optimal solution at $x=0$.

As with the width of the probability density, let us also introduce a parameter associated with $j(x,t)$ to monitor the time evolution of the current. Define the amplitude as the largest displacement of $j(x,t)$ from zero in the region $x<0$. This is illustrated in Fig. 9(b), where the amplitude is $j(x_{0},t)$ at the point $x_{0}$. The amplitude can take on negative values, with the sign indicating the direction of probability flow.

Figure 10(b) shows the time evolutions of the current amplitude $j(x_{0},t)$ associated with the time evolutions of the width shown in panel (a). We see that the amplitude oscillates in time between positive and negative values, meaning that the probability current flow towards and away from the global minimum in a periodic manner. This is consistent with the pulsating behavior exhibited by the width. Another feature is that as $T$ increases, the amplitude tends to zero, meaning that the current itself approaches zero both spatially, as well as temporally throughout the entire duration of annealing. As $j(x,t)$ is identically zero for stationary states, this is a sign of the annealing process approaching adiabaticity. Viewed differently, even for the moderate value of $T=560$, where the residual energy approaches the asymptotic $T^{-2}$ law for Stage 1 as seen in Fig. 7 (at $\log_{10}T=2.748$), the current keeps oscillating, see Fig. 10(b), implying a diabatic evolution.

So far, we have focused our discussion on Stage 1. The probability current in Stage 2 is slightly more complicated. An example is shown in Fig. 21, where a snapshot of $j(x,t)$ at $t=10$ with $T=320$ is shown. One sees that the current direction varies spatially, unlike the spatially-coherent flow we have seen in Fig. 9(b). Despite such differences, the overall behaviors in Stages 2 and 3 are generally quite similar to Stage 1. For completeness, these results are discussed in Appendix B.

In the previous section, we have seen that asymptotically the residual energy decays with the total annealing time as $T^{-2}$. Analysis related to the adiabatic theorem shows that this scaling is largely a result of annealing with the linear schedule Eq. (7) Suzuki and Okada (2005); Morita (2007). Nonlinear schedules have been proposed in which the residual energy decays much faster with the annealing time, leading to better efficiency Morita (2007). In this section, we apply them to the potential $V_{\rm SK}(x)$.

Nonlinear schedules can be obtained by generalizing Eq. (7) as

where $n$ denotes the order of the schedule, and the linear function $t/T$ is replaced by higher-order functions $f_{n}(t/T)$. Details concerning the derivation of $f_{n}$ can be found in Ref. Morita (2007). Physically speaking, these functions are designed such that annealing is constrained to proceed slowly both at the beginning as well as at the end of the process. Earlier studies on spin systems revealed that the residual energy decays with annealing time as $T^{-2n}$, so higher-order schedules generally lead to better annealing performances. In the following, we shall investigate the performances of the functions $f_{1}(s)$ to $f_{4}(s)$ given in Ref. Morita (2007),

Figure 11 compares the residual energies obtained under the schedules given by Eq. (10), with $n=1$ to 3. The results from Stages 1 to 3 are shown in panels (a) to (c) of Fig. 11, respectively. The general behaviors in Stages 1 and 3 are qualitatively quite similar. It is seen that for an $n$th-order schedule the residual energy decays as $T^{-2n}$, verifying the theoretical prediction of Ref. Morita (2007). To emphasize the importance of varying the schedule slowly both at the start as well as towards the end of annealing, we also considered a straightforward quadratic schedule $m^{(q)}(t)$, partly motivated by Ref. Stella et al. (2005), which is obtained by simply substituting $f_{n}$ with $(t/T)^{2}$ in Eq. (10). This schedule proceeds slowly only at the start but not at the end of annealing. In panels (a) and (c), we see that the residual energy of $m^{(q)}(t)$ (solid blue squares) decays only as $T^{-2}$. That both linear and quadratic schedules should show the same asymptotic decay rate can also be explained based on the theoretical analysis of Ref. Morita (2007).

The situation in Stage 2 is different from the other two in that the residual energy decays in the manner of the transient phase discussed in Sec. IV.3. Note that this is in the domain of strongly diabatic transitions and so is not the regime in which the function $f_{n}$ were originally designed for. Nevertheless, it is seen from panel (b) of Fig. 11 that the residual energies of $m^{(2)}(t)$ and $m^{(3)}(t)$ decrease faster than that of $m^{(1)}(t)$. However, we also found that the curve for $m^{(4)}(t)$ (not shown) is almost exactly the same as that of $m^{(3)}(t)$, so we should be cautious about generalizing the improvement in performance seen for $n=2$ and 3 to higher orders. We were unable to pursue the curves of $m^{(2)}(t)$ and $m^{(3)}(t)$ to larger $T$ to examine the asymptotic crossover because we have exhausted the numerical precision offered by the simulation package. Within the time span we could study as shown in Fig. 11(b), diabatic processes lead to significantly faster convergence than $T^{-2}$. Interestingly, the quadratic schedule gives the best performance here, unlike in the other two stages. Overall, nonlinear schedules generally show improvement over the linear one in this stage. Interpretation of the results, however, is not as clear-cut as in Stages 1 and 3.

Earlier in Sec. IV.4, we have seen that under the linear schedule, increasing the annealing time suppresses of the probability current and leads to a more adiabatic process. Conventionally, adiabaticity is considered to be closely related with how the global minimum is ultimately to be attained in quantum annealing. With nonlinear schedules, however, we are now allowed to explore the schedule’s order $n$ as a new parameter dimension while keeping the annealing time constant. We shall see that this leads to tunneling, as opposed to adiabaticity, as a dominant mechanism for quantum annealing to arrive at the optimal solution.

Let us focus our discussion on Stage 1. The probability currents under nonlinear schedules also have the spatially-coherent form of Fig. 9(b) (no node in each of the region $x>0$ and $x<0$), so once again we can use the amplitude $j(x_{0},t)$ to monitor their temporal behaviors. Figure 12 shows the time evolutions of $j(x_{0},t)$ obtained under the various schedules, with annealing time $T=560$. The salient feature that is common to all the nonlinear schedules and which distinguishes them from the linear one is that the amplitude is consistently positive throughout the entire duration of annealing. The probability current is therefore always flowing towards the origin, and no current ever flows away. This is a highly diabatic process in which the optimization is targeted at the global minimum. In stark contrast, the current of the linear schedule oscillates in an unguided manner, resulting in inefficient annealing.

Apart from being persistent, the current of nonlinear schedules also has the additional feature of exhibiting tunneling. To explain, let us first note from Fig. 9 that the spatial extent of both the probability density and the current in the domain $x>0$ can be approximated by the width $\sqrt{\langle\psi(t)|x^{2}|\psi(t)\rangle}$. The red curve in Fig. 13 shows the time evolution of $\sqrt{\langle\psi(t)|x^{2}|\psi(t)\rangle}$ for the schedule $m^{(3)}(t)$ in Stage 1 with $T=560$. The (green) shaded areas show the classically forbidden regions where the total energy $\langle\psi(t)|H(t)|\psi(t)\rangle$ is less than the potential energy $V_{\rm SK}(x)$. Shaded areas lying below the red curve are regions where the current resides in classically forbidden regimes, and these currents are tunneling in nature. Hence, nonlinear schedules induce persistent currents that tunnel across potential barriers to greatly improve the efficiency of annealing. This is in contrast to the linear schedule which, on the contrary, seeks to suppress the flow of current in order to bring about a more adiabatic annealing process.

So far, our discussion has focused on Stage 1. For Stage 2, we mentioned earlier that under the linear schedule the current is not spatially-coherent (c.f. Fig. 21). Under nonlinear schedules, however, the current resumes a spatially-coherent form and so we can again use the amplitude $j(x_{0},t)$ as a monitoring parameter to study its temporal evolution. Despite the different behavior of their residual energies, it was found that the observations reported for Stage 1 are valid for Stage 2 as well. For completeness, these results are summarized in Appendix B.

The situation in Stage 3 is quite similar to that of Stage 1, with the current showing persistence under the nonlinear schedules. Tunneling, however, is not a prominent feature in Stage 3 as the wavefunction is already well-localized within the global minimum [c.f. Fig. 4(c)]. The optimization is therefore taking place within an effective quadratic potential, which is a relatively simple problem for quantum annealing.

We perform simulated annealing numerically using the Metropolis algorithm. Let the position of the particle be $x_{t}$ at time $t$. A new position $x^{\prime}$ is proposed as

where $\Delta x$ is a random number drawn from the interval $[-\frac{s}{2},\frac{s}{2}]$ of width $s$. The quantity $s$ is called the Monte Carlo step size and it determines how far the particle is allowed to jump with each move. The proposed position $x^{\prime}$ is then accepted with probability

where the temperature is provided by the annealing schedule $\beta(t)$, which will be specified later. After $x_{t}$ is updated, the time $t$ advances by a time step $dt$ and the Metropolis algorithm is repeated ⁵⁵5The simulated annealing calculations in Secs. VII and VI are performed with $dt=0.1$, to be consistent with the time step used in quantum annealing. .

To compute averages, we work with an ensemble of $N$ particles where each of them evolves independently in time according to the Metropolis algorithm. The average $\langle V_{\rm SK}(x)\rangle_{t}$ at time $t$ is calculated via the ensemble average as

where $x^{\nu}_{t}$ is the position of the $\nu$th particle at time $t$. In practice, one works with a large but finite $N$, typically $N=10^{7}$ and $10^{9}$. The ensemble average then becomes an approximation, subjected to errors due to finite-size fluctuations.

Let us consider the following linear form for $\beta(t)$

As in quantum annealing, we restrict ourselves to the time interval $0\leq t\leq T$, where $T$ again denotes the total annealing time. The parameters $\beta_{i}$ and $\beta_{f}$ are the inverse temperatures at the initial ($t=0$) and final ($t=T$) times, respectively. The schedule $\beta^{(1)}(t)$ is the simulated annealing counterpart to the mass schedule $m^{(1)}(t)$ that we have used in quantum annealing.

We would like to perform the simulated annealing calculations in such a way that they can be compared to the results of quantum annealing obtained in Sec. IV. To enable a fair comparison, we use energy as the criterion. The boundary parameters $\beta_{i}$ and $\beta_{f}$ in Eq. (16) are chosen in such a way that the initial and final energies of a stage in simulated annealing are approximately the same as the initial and final energies of the corresponding stage in quantum annealing. These boundary parameters $\beta_{l}$ and their energies $U(\beta_{l})$ are summarized in Table 1, alongside their corresponding quantum counterparts $m_{l}$ and $E_{0}(m_{l})$. To facilitate our subsequent discussions, let us call annealing from $\beta_{a}$ to $\beta_{b}$ Stage A, that from $\beta_{b}$ to $\beta_{c}$ Stage B, and that from $\beta_{c}$ to $\beta_{d}$ Stage C. The equilibrium state properties at the stage boundaries have been discussed in Sec. III [c.f. Figs. 3(b), 4(b), and 4(d)]. A graphical perspective of the three stages is also given in Fig. 5(b), where they are superposed upon the curve of the internal energy $U(\beta)$. It is seen that the stages are positioned in the non-trivial regime of the system exhibiting ‘excess energy’.

The annealing performance is assessed by considering the residual energy

The first term is the average potential energy obtained when annealing ends at $t=T$. The second term is the average potential energy at equilibrium at the final temperature [c.f. Eq. (5)], and is the target energy that annealing seeks to attain. The quantity $R_{\rm s.a.}(T)$ is the simulated annealing counterpart to the residual energy $R_{\rm q.a.}(T)$ considered in quantum annealing.

Figure 7 shows the behavior of $R_{\rm s.a.}(T)$ as a function of $T$ in Stage A under the linear schedule. We experimented with a wide range of step sizes and some examples of the residual energy curves we obtained are shown in the figure. There exists an optimum step size that gives the lowest-lying $R_{\rm s.a.}(T)$ curve, and in Stage A this step size was found to be $s=1$ (red circles). The solid line (black) is obtained by fitting a straight line to the data points of $s=1$. We see that the residual energy decreases as $T^{-1}$. This scaling of $R_{\rm s.a.}(T)$ with $T$ is quite general. We identified the optimum step sizes in Stages B and C as well, and the lowest-lying $R_{\rm s.a.}(T)$ curve in all three stages are summarized in Fig. 15, shown by the three slanting curves. The two solid lines (black) are obtained by fitting a straight line to the data points of Stages B and C. It is seen that the residual energy also decreases as $T^{-1}$ in these two stages.

The scaling of $T^{-1}$ is fairly robust in the sense that it is attainable over quite a wide range of step sizes. However, in Stages B and C, we also found that if the step size falls below the inter-minima distance of $V_{\rm SK}(x)$ (i.e. $\frac{s}{2}<w_{0}$), then the $T^{-1}$ scaling will no longer be achieved. In Fig. 15, the two uppermost curves (magenta squares) show the residual energies obtained in these two stages with the annealing performed using $s=0.1$. It is seen that they decrease much slower than the previous curves. The reason for the slow decay is trapping by local minima, and can be understood by looking at the initial Boltzmann distributions of these two stages, shown in Figs. 4(b) and 4(d). The probability channels connecting the global minimum to the two adjacent local ones are suppressed, so if the step size is smaller than the inter-minima distance $w_{0}$, particles at the local minima will need a long time to hill-climb the potential barrier in order to reach the global minimum, possibly leading to a logarithmic decrease of the residual energy as suggested by Shinomoto and Kabashima Shinomoto and Kabashima (1991). See also discussions in the next section. Viewed differently, larger step sizes satisfying $\frac{s}{2}>w_{0}$ may be regarded to realize quasi-global searches across energy barriers without hill-climbing processes, leading to a faster polynomial decrease of the residual energy. It is noteworthy that, even with such a quasi-global search, the decrease rate $T^{-1}$ is slower than the quantum counterpart $T^{-2}$.

One might inquire if it is possible to obtain a better scaling than $T^{-1}$ by making some adjustments to the annealing algorithm. We tried a simple quadratic schedule where the linear term $t/T$ in Eq. (16) is squared. We also experimented with a variant algorithm in which the step size decreases adaptively with temperature during the course of annealing. Overall, we did not manage to achieve a faster decrease than $T^{-1}$ with these attempts.

In Ref. Shinomoto and Kabashima (1991), Shinomoto and Kabashima studied simulated annealing of the original form of $V_{\rm SK}(x)$, where the barrier height and distance between neighboring minima are kept constant for all $x$, under a logarithmic annealing schedule with time running from 0 to $\infty$ in contrast to the finite-time protocols in the previous sections. To test their prediction, let us consider the following schedule for the inverse temperature

where $\beta_{i}$ is the initial inverse temperature at time $t=0$. Time $t$ is supposed to run from 0 to $\infty$. The superscript $(l)$ denotes ‘logarithmic’ and is to differentiate Eq. (18) from the other schedules. Shinomoto and Kabashima showed, by assuming that particle are located only at local minima and taking the continuum limit ignoring discreteness of positions, that the average of the potential energy during simulated annealing under Eq. (18) should decrease with time asymptotically as

where $\langle V_{\rm SK}(x)\rangle_{t}$ is given by the right side of Eq. (5), with a small difference that $\beta$ is now further dependent on time via the annealing schedule. We indicate this time dependence using the subscript $t$ in $\langle\cdot\rangle_{t}$. In the following, we shall study the dependence of $\langle V_{\rm SK}(x)\rangle_{t}$ on time from the perspective of Monte Carlo simulations to see if their prediction is observed numerically.

Figure 16 shows the decay of $\langle V_{\rm SK}(x)\rangle_{t}$ with time under Eq. (18), with $\beta_{i}=10^{0.29}$. The initial ensemble at the start of annealing is the one drawn from the Boltzmann distribution of $V_{\rm SK}(x)$ with temperature $\beta_{i}$. Each curve in the figure represents the result obtained using a particular step size $s$. It is seen that when $t$ is large, the curves collapse onto a single curve. We fitted a straight line to the curve of $s=4$ using the data points from $\log_{10}\log_{10}t>0.5$, and the result is shown by the solid line (black). We see that $\langle V_{\rm SK}(x)\rangle_{t}$ decays asymptotically as $(\log_{10}t)^{-0.74}$. The decay exponent $-0.74$ is different from the $-1$ proposed by Shinomoto and Kabashima.

To account for our results, we propose generalizing the scaling Eq. (19) slightly as

where the decay exponent $\alpha_{\rm s.a.}$ is obtained from numerical simulations. We repeated the simulated annealing calculations with different values of $\beta_{i}$ and the results are summarized in Table 2. These results capture the non-trivial aspects of the potential $V_{\rm SK}(x)$. The column $\beta_{\rm final}$ indicates the final inverse temperature attained by the schedule for that annealing simulation. By referring $\beta_{\rm final}$ to Fig. 5(b), we see that our annealing was performed in the regime where the system exhibits ‘excess energy’.

The exponents $\alpha_{\rm s.a.}$ may be understood on the basis of the equilibrium properties of the system. Let us define

where $\langle V_{\rm SK}(x)\rangle_{\beta}$ is given by Eq. (5) and is the average potential energy at equilibrium. The quantity $\alpha_{\rm eq.}(\beta^{\prime})$ is the gradient of the curve $\log_{10}\langle V_{\rm SK}(x)\rangle_{\beta}$ at $\log_{10}\beta^{\prime}$. In Table 2, the column $\alpha_{\rm eq.}(\beta_{\rm final})$ shows the gradient evaluated at the final attained temperature $\beta_{\rm final}$, and we see that the values agree quite well with the decay exponent $\alpha_{\rm s.a.}$. This agreement between the two quantities can be explained by noting that when $t$ is large the schedule Eq. (18) varies very slowly. The system is therefore very close to equilibrium during the late stages of annealing, and it traces out a narrow segment of the $\langle V_{\rm SK}(x)\rangle_{\beta}$ curve around $\beta_{\rm final}$. The relation Eq. (20) can therefore be approximated by $\langle V_{\rm SK}(x)\rangle_{\beta_{\rm final}}\sim(\beta_{\rm final})^{\alpha_{\rm s.a.}}$, from which the interpretation of $\alpha_{\rm s.a.}$ as the gradient $\alpha_{\rm eq.}(\beta_{\rm final})$ follows immediately via Eq. (21).

We close this section with some comments on the decay exponent $-1$ obtained by Shinomoto and Kabashima and the $\alpha_{\rm s.a.}$ we obtained. One possible reason for the discrepancy is that the potential $V_{\rm SK}(x)$ we used, although inspired by Shinomoto and Kabashima’s work, is not exactly the same as the one they treated. Another cause may be the coarse-graining assumption $(w_{0})^{2}\beta_{\rm final}\ll 1$ which the authors made during their derivation. Substituting $w_{0}=0.2$ and the values of $\beta_{\rm final}$ in Table 2, we see that the inequality is not very well-satisfied for the numerical simulations which we performed here. A third, and possibly more important, reason may be as follows. If the system is close to equilibrium and the particle stays almost exclusively in the valley around the global minimum at $x=0$ for $t$ large and thus $\beta$ large, then the equipartition law suggests that the average of the potential energy is proportional to $\beta^{-1}$, which would imply $\langle V_{\rm SK}(x)\rangle_{\beta}\sim(\log_{10}t)^{-1}$ under the annealing schedule of Eq. (18). Also the kinetic energy term leads to the same behavior due to equipartition. Deviations from this inverse-log law in numerical simulations would imply a few possibilities including a spreading of probability distribution away from the central valley to nearby local minima as well as deviations from equilibrium. It is anyway the case that the contribution to the equilibrium internal energy from kinetic energy would yield the inverse-log law under the quasi-equilibrium context, which is more dominant than the potential term with larger powers as listed in Table 2.