Proposal of a new quantum annealing schedule for studying quantum annealing of transverse field Ising models

TL;DR

This paper introduces a new quantum annealing schedule based on the smallest effective transverse field, demonstrating improved performance over classical annealing in Ising spin glass models through Monte Carlo simulations.

Contribution

A novel QA schedule utilizing the smallest effective transverse field is proposed, enabling better comparison and performance evaluation against classical annealing.

Findings

QA outperforms CA with sufficient annealing time

The new schedule provides a more effective comparison framework

Monte Carlo simulations confirm the schedule's advantages

Abstract

Recently, Heim, Ronnow, Isakov and Troyer [Science 348 (2015) 215] have reported that Monte Carlo simulations for the Ising spin glass model on the square lattice in the physically relevant continuous-imaginary-time limit do not show superiority of quantum annealing (QA) using transverse field against classical annealing (CA). Although the QA schedule that they had used has been using conventionally, however the QA schedule mathematically has no guarantee that the used schedule is the best QA schedule for performance of optimization. We propose a new QA schedule for studying transverse-field-based quantum versus classical annealing of the Ising model. The present QA schedule utilizes a smallest effective transverse field derived in this article. This QA schedule is made for the comparison between the system with no transverse field and the system with the smallest effective transverse field. As a case study, we study QA of the Ising spin glass model on the square lattice at low but finite temperature. A Monte Carlo algorithm using the physically relevant continuous-imaginary-time limit is performed. As the simulation results, we show superiority of QA against CA when the annealing time is sufficiently spent.

Equations

$$\mathcal{H}(t)=J(t){\mathcal{H}}_P+\gamma (t){\mathcal{H}}_Q,$$

$${\mathcal{H}}_P=-\sum _{<ij>}{\tau }_{ij}{S}_i{S}_j,$$

$$P({\tau }_{ij})=\frac{1}{2}{\delta }_{{\tau }_{ij},-1}+\frac{1}{2}{\delta }_{{\tau }_{ij},1}$$

$${\mathcal{H}}_Q=-\sum _{i=1}^N{\sigma }_i^x,$$

$$\displaystyle\sigma^{x}=\left(\begin{array}[]{cc}0&1\\ 1&0\end{array}\right),\quad\sigma^{y}=\left(\begin{array}[]{cc}0&-i\\ i&0\end{array}\right)\quad{\rm and}\quad\sigma^{z}=\left(\begin{array}[]{cc}1&0\\ 0&-1\end{array}\right)\,.$$

$$\mathcal{H}(t)=-J(t)\sum _{<ij>}{\tau }_{ij}{\sigma }_i^z{\sigma }_j^z-\gamma (t)\sum _{i=1}^N{\sigma }_i^x.$$

$$e^{-\beta ({\mathcal{H}}_A+{\mathcal{H}}_B+{\mathcal{H}}_C)}=\underset{\Delta \tilde{t}\to 0}{\lim }(e^{-{\mathcal{H}}_A\Delta \tilde{t}}{e}^{-{\mathcal{H}}_B\Delta \tilde{t}}{e}^{-{\mathcal{H}}_C\Delta \tilde{t}}{)}^{\frac{\beta }{\Delta \tilde{t}}},$$

$$⟨{\sigma }_x^z|e^{\beta \gamma (t){\sigma }^x}|{\sigma }_y^z⟩=\frac{1}{2}(e^{\beta \gamma (t)}+{\sigma }_x^z{\sigma }_y^z{e}^{-\beta \gamma (t)}),$$

$$\displaystyle Z(t)=\lim_{\Delta\tilde{t}\to 0}\sum_{\{\sigma^{z}_{ik}\}}\exp\biggl{\{}J(t)\Delta\tilde{t}\sum_{k\in\rho}\sum_{ij\in\rho}\tau_{ij}\sigma^{z}_{ik}\sigma^{z}_{jk}$$

$$\displaystyle+\sum_{kl\in\rho}\sum_{i\in\rho}\log\frac{\exp(\gamma(t)\Delta\tilde{t})+\sigma^{z}_{ik}\sigma^{z}_{il}\exp(-\gamma(t)\Delta\tilde{t})}{2}\biggr{\}}\,,$$

$$Z(t)$$

$$w(\sigma^{z}_{ik},\sigma^{z}_{il})=\exp\biggr{[}\log(\gamma(t)\Delta\tilde{t})\,\frac{1-\sigma^{z}_{ik}\sigma^{z}_{il}}{2}\biggr{]}\,.$$

$$w_{\mathrm{para}}=1.$$

$$w_{\mathrm{anti}}=\gamma (t)\Delta \tilde{t}.$$

$$w_{\mathrm{para}}=w(g_{\mathrm{conn}})+w(g_{\mathrm{cut}}),$$

$$w_{\mathrm{anti}}=w(g_{\mathrm{cut}}).$$

$$w(g_{\mathrm{conn}})=1-\gamma (t)\Delta \tilde{t},$$

$$w(g_{\mathrm{cut}})=\gamma (t)\Delta \tilde{t}.$$

$$P_{\mathrm{para}}(g_{\mathrm{cut}})=\frac{w(g_{\mathrm{cut}})}{w(g_{\mathrm{conn}})+w(g_{\mathrm{cut}})}=\gamma (t)\Delta \tilde{t},$$

$$P_{\mathrm{anti}}(g_{\mathrm{cut}})=\frac{w(g_{\mathrm{cut}})}{w(g_{\mathrm{cut}})}=1.$$

$$P(n_{\mathrm{add}})$$

$$=$$

$${\tilde{t}}_{\mathrm{add}}=\beta R,$$

$$n_{\mathrm{cut}}=n_{\mathrm{add}}+n_{\mathrm{kink}}.$$

$$P(E\to E^{\prime })=\frac{e^{-\beta E^{\prime }}}{e^{-\beta E}+e^{-\beta E^{\prime }}}.$$

$$[⟨E{⟩}_T{]}_J=-\frac{\partial }{\partial \beta }[\log Z{]}_J$$

$$E=-\frac{J(t)}{\beta }\sum _{<ij>}{\int }_0^{\beta }d\tilde{t}{\tau }_{ij}{\sigma }_i^z(\tilde{t}){\sigma }_j^z(\tilde{t})-\frac{1}{\beta }\sum _{i=1}^N{n}_{\mathrm{kink}}(i)$$

$$S_{i}^{\rm result}=\frac{\int_{0}^{\beta}\!d\tilde{t}\,\sigma^{z}_{i}(\tilde{t})}{\big{|}\int_{0}^{\beta}\!d\tilde{t}\,\sigma^{z}_{i}(\tilde{t})\big{|}}\,.$$

$$E_{P}=-\sum_{<ij>}\tau_{ij}\,\frac{\int_{0}^{\beta}\!d\tilde{t}\,\sigma^{z}_{i}(\tilde{t})}{\big{|}\int_{0}^{\beta}\!d\tilde{t}\,\sigma^{z}_{i}(\tilde{t})\big{|}}\,\frac{\int_{0}^{\beta}\!d\tilde{t}\,\sigma^{z}_{j}(\tilde{t})}{\big{|}\int_{0}^{\beta}\!d\tilde{t}\,\sigma^{z}_{j}(\tilde{t})\big{|}}$$

$$E_{\mathrm{res}}=[⟨E_P{⟩}_T-E_G{]}_J.$$

$$P_{\mathrm{exact}}=[⟨{\delta }_{E_P,E_G}{⟩}_T{]}_J.$$

Full text

1 Introduction

The quantum annealing (QA) is an optimization method utilizing quantum effect[1, 2, 3, 4, 5, 6]. By using QA, the solutions of given optimization problems are approximately or exactly obtained. QA is also related to a quantum computation by adiabatic evolution[7, 8].

QA is proposed as an alternative of the classical annealing (CA)[9]. In CA, the temperature is operated from high to low in order to utilize thermal fluctuations. On the other hand, in QA, external fields (and/or exchange interactions) generating quantum effect are operated from large to zero in order to utilize quantum fluctuations. These annealings help to make the systems escape from local minima of the free energy, and the systems reach low-energy states effectively.

The understandings of QA are recently becoming significant in relation to the D-Wave chip[10]. This chip is designed to perform QA for solving optimization problems. The occurrence of quantum properties for this chip is shown[10].

We simulate QA of real-physical systems in this article, and we do not simulate QA as a classical optimization method in this article. Recently, it is reported in Ref. [11] that Monte Carlo simulations for the Ising spin glass model on the square lattice in the physically relevant continuous-imaginary-time limit do not show superiority of QA using transverse field against CA. In order to study physically relevant systems, we use a continuous-imaginary-time quantum Monte Carlo algorithm in this article. In addition, we investigate the Ising spin glass model on the square lattice.

Although the QA schedule used in Ref. [11] has been using conventionally, however the QA schedule mathematically has no guarantee that the used schedule is the best QA schedule for performance of optimization. We propose a new QA schedule for studying transverse-field-based quantum versus classical annealing of the Ising model. This QA schedule utilizes a smallest effective transverse field derived in this article. We derive this field based on a discussion of percolation of spin correlation per spin along the continuous-imaginary-time direction. The present QA schedule is made for the comparison between the system with no transverse field and the system with the smallest effective transverse field. By this comparison, the superiority of QA (or CA) can be made clear.

We apply a quantum Monte Carlo algorithm [14] proposed by Nakamura and Ito. This algorithm uses the physically relevant continuous-imaginary-time limit. This algorithm is a Glauber-dynamics-like algorithm (a heat-bath-like algorithm) which performs single-spin flips in the real space and cluster flips in the imaginary-time space. The Glauber dynamics is suited to the study of the dynamical features of physical systems[12, 13]. By applying this algorithm, observing a true relaxation of the original system is expected[14]. The mathematical form of this algorithm is directly related to the derivation of a smallest effective transverse field derived in this article, thus we use this algorithm in this article

We compare three annealings: the present one, a conventional quantum one (used for a quantum computation by adiabatic evolution) and a classical one. The simulation results for the annealings are shown in this article.

This article is organized as follows. The quantum annealing and the spin glass model are explained in §2. The derivation of a quantum Monte Carlo algorithm by Nakamura and Ito [14] is described in §3. The way of applying this algorithm is written in §4. The differences between this algorithm and the related algorithms [15, 16] are mentioned in §4. Quantities that we treat in this article are described in §5. A smallest effective transverse field is derived in §6. A new QA schedule is proposed in §6. Simulation results are shown in §7. Concluding remarks are in §8.

2 The quantum annealing and the spin glass model

In terms of theoretical physics, the quantum annealing can be to investigate a time dependent Hamiltonian ${\cal H}(t)$ which has quantum effect terms. Usually, the Hamiltonian ${\cal H}(t)$ can be written as

[TABLE]

where $t$ is the time from $t=0$ to $t=T_{I}$. ${\cal H}_{P}$ is the problem Hamiltonian that an optimization problem is written. ${\cal H}_{Q}$ is a quantum effect Hamiltonian for generating quantum effect. When $T_{I}$ is a small number, this annealing corresponds to a fast annealing, and, when $T_{I}$ is a large number, this annealing corresponds to a slow annealing. $J(t)$ is an increasing function for $t$, or $J(t)$ is a constant. $J(t)$ increases from [math] to $1$ for example. $\gamma(t)$ is a decreasing function for $t$. $\gamma(t)$ decreases from a large value to zero for example, or $\gamma(t)$ decreases from $1$ to [math] for example. A method using a pulse of the transverse field is also proposed[16].

When $\gamma(t)$ gives zero from $t=0$ to $t=T_{I}$ while the value of $J(t)$ (or the temperature) changes, this annealing is a classical annealing. When there is a time that $\gamma(t)$ gives non-zero, this annealing is a quantum annealing.

The problem Hamiltonian ${\cal H}_{P}$ is written as the form of the Ising spin glass model called the Edwards-Anderson model. The Hamiltonian ${\cal H}_{P}$ is given by [17, 18, 19, 20]

[TABLE]

where $<ij>$ denotes nearest-neighbor pairs, $S_{i}$ is a state of a spin at site $i$, $S_{i}=\pm 1$, and this spin is called the Ising spin. The value of $\tau_{ij}$ is given by the problem that is asked to solve. In the spin glass model, the probability of giving the value of $\tau_{ij}$ is

[TABLE]

for example. Here, $\delta$ is the Kronecker delta. We use the value of $\tau_{ij}$ obtained by using Eq. (3) in this article. The Ising spin glass model using $P(\tau_{ij})$ is called the $\pm J$ Ising spin glass model. We investigate the $\pm J$ Ising spin glass model on the square lattice with periodic boundary conditions in this article.

The quantum effect Hamiltonian ${\cal H}_{Q}$, which is investigated in this article, is given by [1, 3, 4, 5]

[TABLE]

where $\sigma^{k}_{i}$ is the k component of the Pauli matrix at site i. $N$ is the number of sites (spins) in the whole system. The Pauli matrices are

[TABLE]

When the measurement axis is chosen to be $z$, the external field by $\sigma^{x}$ is called the transverse field. An Ising spin under the transverse field corresponds to a quantum bit in this study.

By using Eqs.(1), (2) and (4), the Hamiltonian ${\cal H}(t)$ is also written as

[TABLE]

We study this Hamiltonian ${\cal H}(t)$.

3 The derivation of a continuous-imaginary-time quantum Monte Carlo algorithm

The derivation of a quantum Monte Carlo algorithm by Nakamura and Ito [14] is described in this section, since the mathematical form of this algorithm is directly related to the derivation of a smallest effective transverse field derived in this article. The way of applying this algorithm is described in §4.

The Glauber dynamics is suited to the study of the dynamical features of physical systems[12, 13]. By applying the Glauber dynamics, discussions for dynamical features as a physical system are possible[12, 13]. The Glauber dynamics is a heat-bath algorithm for the Ising model. We use a Glauber-dynamics-like algorithm (a heat-bath-like algorithm) which performs single-spin flips in the real space and cluster flips in the imaginary-time space. For the real-space direction, this algorithm performs single-spin flips like as the Glauber dynamics does. For the imaginary-time direction, the continuous-imaginary-time limit is applied.

We use the Suzuki-Trotter decomposition given by [21]

[TABLE]

where $\Delta\tilde{t}\equiv\frac{\beta}{M}$, $M$ is called the Trotter number, and $\tilde{t}$ is called the imaginary time. $\beta$ is the inverse temperature, $\beta=1/(k_{B}T_{E})$, $T_{E}$ is the temperature, and $k_{B}$ is the Boltzmann constant. The physically relevant continuous-imaginary-time limit mentioned in Ref. [11] corresponds to the limit of $\Delta\tilde{t}\to 0$ in this article. By this decomposition, the $d$-dimensional quantum system is treated as the $d+1$-dimensional classical system. A continuous-imaginary-time quantum Monte Carlo algorithm is used in this study. This means that the system for $\Delta\tilde{t}\to 0$ is directly investigated. This study is not affected by the error of the Suzuki-Trotter decomposition. There is a relation[21]:

[TABLE]

where $\sigma^{z}_{x},\sigma^{z}_{y}=\pm 1$. Thus, the partition function $Z(t)$ is written as [21]

[TABLE]

where $\rho$ is the set of the coordinates of the weights for the decomposed partition function in the $d+1$-dimensional classical system. The coordinates of $\rho$ is somewhat complicated[21], but there is no need to understand the coordinates, because of taking the continuous-imaginary-time limit. By performing $\Delta\tilde{t}\to 0$, the imaginary time becomes a continuous space. In Eq. (15), the term for $\sigma^{z}_{ik}\sigma^{z}_{jk}$ is a term for the interaction of real-space direction, and the term for $\sigma^{z}_{ik}\sigma^{z}_{il}$ is a term for the interaction of imaginary-time direction. The imaginary-time direction has periodic boundary conditions. The size of the imaginary-time direction is the same with the inverse temperature $\beta$. Eq. (15) is also written for $\Delta\tilde{t}\to 0$ as

[TABLE]

In order to derive this algorithm, we derive weights for graph representation. The framework of deriving weights for graph representation is described in Ref. [22]. We define the weight of two spins along imaginary-time direction as $w(\sigma^{z}_{ik},\sigma^{z}_{il})$. $w(\sigma^{z}_{ik},\sigma^{z}_{il})$ is given by

[TABLE]

We define the weight for $\sigma^{z}_{ik}\sigma^{z}_{il}=1$ as $w_{\rm para}$. By using Eq. (17), we obtain

[TABLE]

We define the weight for $\sigma^{z}_{ik}\sigma^{z}_{il}=-1$ as $w_{\rm anti}$. By using Eq. (17), we obtain

[TABLE]

We define the weight of graph for connecting two spins as $w(g_{\rm conn})$. We define the weight of graph for cutting two spins as $w(g_{\rm cut})$. We are able to write

[TABLE]

and

[TABLE]

By using Eqs. (18) - (21), we obtain

[TABLE]

and

[TABLE]

We define the probability of cutting two spins for $\sigma^{z}_{ik}\sigma^{z}_{il}=1$ as $P_{\rm para}(g_{\rm cut})$. We define the probability of cutting two spins for $\sigma^{z}_{ik}\sigma^{z}_{il}=-1$ as $P_{\rm anti}(g_{\rm cut})$. We are able to write

[TABLE]

and

[TABLE]

We calculate the probability of giving the number of cuts per spin for $P_{\rm para}(g_{\rm cut})$ of Eq. (24) in the continuous-imaginary-time limit. We define $n_{\rm add}$ as the number of cuts per spin for parallel spins $\sigma\,\sigma^{\prime}=1$ along the imaginary-time direction. By using Eq. (24), the probability $P(n_{\rm add})$ of giving $n_{\rm add}$ is obtained as

[TABLE]

where \left(\begin{array}[]{c}x\\ y\end{array}\right)\equiv\frac{x!}{y!(x-y)!}, and $n_{\rm kink}$ is the number of kinks per spin. The kink means the position of antiparallel spins $\sigma\,\sigma^{\prime}=-1$ along the imaginary-time direction. The position $\tilde{t}_{\rm add}$ of the cut for parallel spins $\sigma\,\sigma^{\prime}=1$ along the imaginary-time direction is obtained by

[TABLE]

where $R$ is a pseudorandom number for $0\leq R<1$. By using Eq. (25), cuts are added for antiparallel spins $\sigma\,\sigma^{\prime}=-1$ along the imaginary-time direction with probability one. Thus the number of cuts per spin for antiparallel spins along the imaginary-time direction is $n_{\rm kink}$. Therefore, the number of all cuts per spin along the imaginary-time direction, $n_{\rm cut}$, is given by

[TABLE]

This algorithm performs single-spin flips in the real space. The transition probability $P(E\to E^{\prime})$ for the energy transition $E\to E^{\prime}$ is given by

[TABLE]

Cuts are added along the imaginary-time direction, and the spin correlations along the imaginary-time direction are cut at the positions of the cuts. Therefore, $P(E\to E^{\prime})$ is only calculated for real spaces. When there is no cut, this algorithm becomes the Glauber dynamics.

4 The way of applying this Monte Carlo algorithm

We describe the way of applying this algorithm. The derivation of this algorithm is in §3. When the transverse field is not imposed, this algorithm becomes the Glauber dynamics.

One Monte Carlo step (sweep) is as follows:

The number of cuts for parallel spins along the imaginary-time direction, $n_{\rm add}$, is calculated by using Eq. (29) for each spin, and the positions of cuts for parallel spins along the imaginary-time direction, $\tilde{t}_{\rm add}$, are calculated $n_{\rm add}$ times by using Eq. (30) for each spin. Moreover, cuts are placed on all kinks with probability one (Eq. (25)). Clusters are detected by all the cuts. The number of clusters, $N_{C}$, is counted up. When there is no cut, $N_{C}$ gives $N$ (the number of spins). 2. 2.

One of the detected clusters is randomly picked up, and the cluster is flipped with the probability of Eq. (32). The trials of flipping clusters are performed $N_{C}$ times. 3. 3.

All the cuts are deleted. 4. 4.

Quantities are sampled if needing to sample the quantities.

One Monte Carlo step is a unit of time used in this study. This time is called the Monte Carlo time. We investigate physical systems under this time evolution.

Fig. 1 shows a sketch of one Monte Carlo step. Full lines correspond to worldline segments with spin up, and broken lines correspond to worldline segments with spin down. The horizontal axis indicates the real space, and the longitudinal axis indicates the imaginary time. A discrete space is assumed for the real-space direction, and a continuous space is assumed for the imaginary-time direction. The imaginary-time direction has periodic boundary conditions. (a) Spin states are shown. (b) Cuts are added. Cuts for kinks are represented as bars, and cuts for parallel spins along the imaginary-time direction are represented as crosses. Cuts for kinks are added with the probability one according to Eq. (25). Cuts for parallel spins along the imaginary-time direction are added according to Eqs.(29) and (30). (c) Cluster flips are performed. One of clusters is randomly picked up, and the cluster is flipped with the probability of Eq. (32). Cluster flips are tried for the number of clusters. (d) Cuts are deleted, and quantities are sampled.

Fig. 2 shows a sketch of one cluster flip. The horizontal axis indicates the real space, and the longitudinal axis indicates the imaginary time. (a) Spin states are shown before a cluster $C_{1}$ is flipped. (b) Spin states are shown after the cluster $C_{1}$ was flipped. If the transition probability $P(E\to E^{\prime})$ of Fig. 2 is calculated, $P(E\to E^{\prime})$ is obtained as $P(E\to E^{\prime})=1/\{1+\exp[-2J(t)(\tilde{t}_{1}-\tilde{t}_{3})\tau_{12}+2J(t)(\tilde{t}_{1}-2\tilde{t}_{2}+\tilde{t}_{3})\tau_{23}]\}$.

One simulation is as follows:

A problem (a set of $\tau_{ij}$) is prepared. 2. 2.

Initial states of spins at $t=0$ are made by using some Monte Carlo steps, and quantities are sampled for $t=0$. 3. 3.

One Monte Carlo step is performed. $t\leftarrow t+1$ is executed. Quantities are sampled for $t$. 4. 4.

When $t=T_{I}$, the simulation is closed.

By repeating the above procedures, many simulations are performed, and the average values of quantities are calculated.

We describe the differences between related algorithms. The algorithm by Rieger and Kawashima is that for each site one generates new cuts in addition to the old ones from the already existing segments via a Poisson process with decay time $1/\gamma$ along the imaginary-time direction[15]. This algorithm directly detects the correlation length along the imaginary-time direction. On the other hand, this algorithm directly treats the number of cuts for detecting the correlation lengths along the imaginary-time direction. This algorithm generates the number of cuts for each site, and this algorithm adds the cuts randomly. The algorithm by Morita, Suzuki and Nakamura is that the initial cluster state is made by giving the position of cuts by the Poisson process with the mean value $\beta\gamma$ [16]. If the number of generated domain walls is odd, one of the domain walls is removed to obey the periodic boundary conditions along the imaginary-time direction[16]. This algorithm does not judge whether the number of generated domain walls is odd or even. If this process for removing one of the domain walls was combined with this algorithm, this combined algorithm gives incorrect results. Therefore, this algorithm and the algorithm by Morita, Suzuki and Nakamura should be considered to be separated.

The simulation results by these algorithms can be the same for all situations. The orders of calculation costs for this algorithm and the related algorithms can also be the same for the system sizes, thus we do not compare the algorithms in detail in this article. The aim of this article is to propose a new quantum annealing schedule and to investigate the efficiency of the schedule. The mathematical form of this Monte Carlo algorithm is directly related to the derivation of a smallest effective transverse field, and the schedule that we propose in this article utilizes the field. We derive the field in §6.

5 Quantities

We describe the quantities that we investigate in this article.

By using Eq. (16), the expectation value $[\langle E\rangle_{T}]_{J}$ of the energy $E$ of the Hamiltonian ${\cal H}(t)$ (Eq. (12)) is written as

[TABLE]

where $\langle\,\rangle_{T}$ is the thermal average, and $[\,]_{J}$ is the average for $\tau_{ij}$. $[\,]_{J}$ is the average for the problems that are asked to solve, and, if the word for the spin glass model is used, $[\,]_{J}$ is the random configuration average. Therefore, in order to calculate the expectation value of the energy by the Monte Carlo simulation,

[TABLE]

is sampled, where $n_{\rm kink}(i)$ is the number of kinks at site $i$. By the way, the solution of a given optimization problem at site $i$ can be obtained as

[TABLE]

The aim of the quantum annealing is to solve optimization problems. Therefore, by using Eqs.(2) and (35), investigating an energy $E_{P}$ given by

[TABLE]

can be suited for this study, since, if a lot of kinks remains in the system after annealing, there is a possibility that $E$ of Eq. (34) gives an unexpected low value. In addition, there is a possibility that this unexpected low value gives a false impression for success of the quantum annealing. Thus we investigate the energy $E_{P}$ instead of the energy $E$ in this article.

We define the exact ground-state energy of Eq. (2) as $E_{G}$, where the exact ground-state is the exact solution of the given problem. $E_{G}$ depends on $\tau_{ij}$. The residual energy $E_{\rm res}$ is given by

[TABLE]

We investigate the residual energy $E_{\rm res}$ in this article.

The accuracy $P_{\rm exact}$ is given by

[TABLE]

We investigate the accuracy $P_{\rm exact}$ in this article.

6 A new quantum annealing schedule

The mathematical form of this Monte Carlo algorithm derived in §3 is directly related to the derivation of a smallest effective transverse field, and the quantum annealing schedule that we propose in this article utilizes the field. Firstly, we derive the field based on a discussion of percolation of spin correlation per spin along the continuous-imaginary-time direction.

By cuts used in Eqs.(25) - (30), spin correlations are cut along the imaginary-time direction. We consider the fewest number of cuts per spin. Since the imaginary-time direction has periodic boundary conditions, by considering a percolation of spin correlations, the fewest number of cuts per spin, $n_{\rm fewest\,cut}$, is obtained as two,

[TABLE]

By using Eqs. (31) and (39), the fewest number of cuts per spin for parallel spins along the imaginary-time direction, $n_{\rm fewest\,add}$, is obtained as

[TABLE]

By using Eq.(29), the average number of cuts per spin for parallel spins along the imaginary-time direction is obtained as

[TABLE]

Therefore, when there is no kink ($n_{\rm kink}=0$), by using Eqs. (40) and (41), we obtain a smallest effective transverse field $\gamma_{S}$ as

[TABLE]

When there are kinks, by using Eq. (40), $n_{\rm fewest\,add}$ is obtained as

[TABLE]

since $n_{\rm kink}>0$. Therefore, when there are kinks, a smallest effective transverse field $\tilde{\gamma}_{S}$ is obtained as

[TABLE]

If there is no kink, the transverse field $\gamma$ for $\gamma<\gamma_{S}$ is not effective for generating quantum effect. $\gamma$ for $\gamma=\gamma_{S}$ is effective for generating quantum effect regardless of the number of kinks. A smallest effect transverse field that does not depend on the number of kinks is desired in this study. Therefore, we estimate that $\gamma_{S}$ is the smallest effective field that should be used in this study.

We propose a QA schedule for studying transverse-field-based quantum versus classical annealing of the Ising model by utilizing the field $\gamma_{S}$ as

[TABLE]

where $H_{E}(x)$ is the Heaviside step function, using the half-maximum convention, which gives $1$ if $x>0$, gives $0.5$ if $x=0$, and gives 0 if $x<0$. $\theta$ is an adjusting value for $0\leq\theta<1$. The value $\theta$ is set to $0.5$ for example. Here, in this study, $t$ is the Monte Carlo time, and $T_{I}$ is the ending Monte Carlo time. The present QA schedule is made for the comparison between the system with no transverse field and the system with the smallest effective transverse field. The inverse temperature $\beta$ is set to a low but finite temperature. For $J(t)$ in Eq. (1), using $J(t)=t/T_{I}$ for $0\leq t\leq T_{I}$ may be the simplest form as used in Ref. [7]. Therefore, we investigate

[TABLE]

as the present-quantum-annealing schedule in this article. The Hamiltonian ${\cal H}_{2}(t)$ for a conventional-quantum-annealing schedule is given by [7]

[TABLE]

We also investigate ${\cal H}_{2}(t)$ as a conventional-quantum-annealing schedule in this article for comparison. As a classical-annealing schedule, we investigate

[TABLE]

in this article for comparison. As for calculating annealing processes, this classical annealing schedule is mathematically the same with a classical annealing for $J=1$, $\beta(t)=\beta t/T_{I}$ and $\gamma=0$ in Eq. (12), thus we investigate ${\cal H}_{3}(t)$.

In Refs. [3, 11], an another quantum-annealing schedule has also been applied. This QA schedule is that the transverse field $\gamma(t)$ is initially much larger than the couplings, $\gamma(0)\gg J(0)|\tau_{ij}|$, $J(t)$ is a constant for $t$, and, during QA, $\gamma(t)$ is slowly reduced to zero. For this annealing schedule, it is reported in Ref. [11] that Monte Carlo simulations for the Ising spin glass model on the square lattice in the physically relevant continuous-imaginary-time limit do not show superiority of QA against CA. We do not investigate this schedule in this article because the results are already shown in Ref. [11]. This schedule is the conventional schedule of the quantum annealing, and the schedule of ${\cal H}_{2}(t)$ is the conventional schedule used in a quantum computation by adiabatic evolution[7, 8]. The schedule of ${\cal H}_{2}(t)$ has not been tested in the framework of the quantum annealing, thus we investigate the schedule of ${\cal H}_{2}(t)$ in this article

A method using a pulse of the transverse field is also proposed in Ref. [16], however the annealing schedule for this method is unclear in Ref. [16]. Thus we do not compare the present annealing schedule with this method in this article.

We investigate the effectivities for ${\cal H}_{1}(t)$, ${\cal H}_{2}(t)$ and ${\cal H}_{3}(t)$ in this article.

In this study, the exchange interaction is tuned instead of the temperature. Even if the temperature is tuned, the temperature is tuned from infinity to a low temperature since our aim is to simulate real-physical systems. It can be considered that the temperature does not reach absolute zero in real-physical systems. It can be reasonable that the final temperature of annealing is a low temperature. In this study, the exchange interaction changes instead of the temperature, thus, if the transverse field is imposed, the annealing becomes QA, and, if the transverse field is not imposed, the annealing becomes CA.

By using Eqs. (31) and (41) on condition of no kink and a fixed number for adding cuts, the transverse field is quantized as

[TABLE]

in the Ising model for generating quantum effect. Even if a transverse field $\gamma_{Q}(n)$ is imposed, the quantum effect increases with increasing the number of kinks. However, to know that the transverse field is quantized as in Eq. (49) can be convenient for adjustment of the quantum effect. $\gamma_{Q}(2)$ is the smallest effective field used in this study ($\gamma_{Q}(2)=\gamma_{S}$). If a transverse field $\gamma_{Q}$ for $n=3$ is imposed, the generation of quantum effect for $n=3$ is guaranteed at least. When a transverse field $\gamma_{\rm ex}$ is imposed for example, the quantum number $n_{\rm ex}$ is obtained as $n_{\rm ex}=\beta\gamma_{\rm ex}$. If $n_{\rm ex}<2$, $\gamma_{\rm ex}$ is estimated as no effective for generating quantum effect. If $n_{\rm ex}\geq 2$, the generation of quantum effect for $n_{\rm ex}$ is guaranteed at least. Since very large $n$ disorders spin orders too much, searching proper $n$ for time $t$ is, of course, required for good performance of optimization. This quantization is a quantization for the number of cutting spin correlation per spin along the continuous-imaginary-time direction. Because the quantum effect by the transverse field depends on the number of the cuts for the correlation lengths along the continuous-imaginary-time direction, we focus on this quantization, and we utilize the fewest quantum number in this article.

7 Simulation results

We firstly tested our program source code for the two-spin-model case in order to fix software bugs, because static properties in the two-spin model are exactly calculated. We secondly tested our program source code for the one-dimensional ferromagnetic Ising spin chain with periodic boundary conditions under the transverse field ($\tau_{ij}=1$ and $\gamma(t)=\gamma$) in order to fix software bugs. The critical properties at $T_{E}=0$ are identical to those of the classical two-dimensional Ising model and are well known exactly (the critical exponents are $\gamma_{C}=1$, $\nu=1$, $\beta=1/8$ and $z=1$) [23]. We confirmed the values of the critical exponents $\nu$ and $\beta$ by using this algorithm and the finite-size-scaling method. We do not show the results of the two-spin-model case and the one-dimensional-ferromagnetic-Ising-spin-chain case, because the simulations were performed for fixing software bugs. Software bugs are fixed as above.

Next, we performed simulations for the $\pm J$ Ising spin glass model on the square lattice with periodic boundary conditions. For each system size and each ending Monte Carlo time, $100$ realizations of $\tau_{ij}$ were investigated. For each realization, $10$ simulations were performed. Thus, $1000$ simulations were performed for each system size and each ending Monte Carlo time. $100$ Monte Carlo steps were used for preparing the initial states at each simulation. The inverse temperature $\beta$ was set to $\beta=10$. The CA schedule, the conventional-quantum-annealing (CQA) schedule and the present-quantum-annealing (PQA) schedule were performed. All the instances used for CA, CQA and PQA are the same.

Fig. 3 shows the relation between the ending Monte Carlo time $T_{I}$ and the residual energy $E_{\rm res}$ given in Eq. (37). The number of spins, $N$, is $100$. The solid square represents the result for the CA schedule, the solid circle represents the result for the CQA schedule, and the open circle represents the result for the PQA schedule. At $T_{I}=200$, the result of CA is better than the results of CQA and PQA. At $T_{I}=1600$, the results of CQA and PQA are better than the result of CA. Therefore, by fast annealing, CA is seen to outperform QA, but, by slow annealing, QA is seen to outperform CA. At $T_{I}=200$ and $T_{I}=400$, the results of PQA are better than the results of CQA.

Fig. 4 shows the relation between the ending Monte Carlo time $T_{I}$ and the accuracy $P_{\rm exact}$ given in Eq. (38). The number of spins, $N$, is $100$. The solid square represents the result for the CA schedule, the solid circle represents the result for the CQA schedule, and the open circle represents the result for the PQA schedule. At $T_{I}=200$, the result of CA is better than the results of CQA and PQA. At $T_{I}=1600$, the results of CQA and PQA are better than the result of CA. Therefore, by fast annealing, CA is seen to outperform QA, but, by slow annealing, QA is seen to outperform CA. At $T_{I}=400$ and $T_{I}=1600$, the results of PQA are better than the results of CQA.

From Figs. 3 and 4, by fast annealing, CA is seen to outperform QA at a number of spins, but, by slow annealing, QA is seen to outperform CA at the number of spins. From Figs. 3 and 4, PQA is seen to outperform CQA.

Fig. 5 shows the relation between the number of spins, $N$, and the residual energy $E_{\rm res}$ given in Eq. (37). The ending Monte Carlo time $T_{I}$ is $1600$. The solid square represents the result for the CA schedule, the solid circle represents the result for the CQA schedule, and the open circle represents the result for the PQA schedule. At $N=100$, the result of QA is better than the result of CA. QA is seen to outperform CA at large system sizes. At $N=36$, the result of PQA is better than the result of CQA. We can see that the results of CQA and PQA are better than the result of CA for $N$ dependence. This shows the superiority of QA against CA. We could not see the superiority of PQA against CQA for $N$ dependence.

Fig. 6 shows the relation between the number of spins, $N$, and the accuracy $P_{\rm exact}$ given in Eq. (38). The ending Monte Carlo time $T_{I}$ is $1600$. The solid square represents the result for the CA schedule, the solid circle represents the result for the CQA schedule, and the open circle represents the result for the PQA schedule. At $N=100$, the result of QA is better than the result of CA. QA is seen to outperform CA at large system sizes. At $N=36$ and $N=100$, the results of PQA are better than the results of CQA. We can see that the results of CQA and PQA are better than the result of CA for $N$ dependence. This shows the superiority of QA against CA. We could not see the superiority of PQA against CQA for $N$ dependence.

From Figs. 5 and 6, QA is seen to outperform CA for the dependences of the number of spins. From Figs. 5 and 6, PQA is seen to outperform CQA.

8 Concluding remarks

We performed Monte Carlo simulations for the spin glass model on the square lattice. As the simulation results, QA was seen to outperform CA at a number of spins when the annealing time was sufficiently spent. As the simulation results, QA was seen to outperform CA for the dependences of the number of spins when the annealing time was sufficiently spent. Therefore, our conclusion is that there is a superiority of QA against CA for the spin glass model by using Monte Carlo simulations.

QA was seen to outperform CA for the system-size dependences. However, we were not able to numerically estimate whether the dependences are power dependences or exponential dependences, because the dependences were not clear enough. Although the numerical estimations for the system-size dependences would be hard tasks, the estimations are tasks for the future.

Conventionally, operating the temperature $T_{E}$ with a fixed $J(t)(=J)$ and no $\gamma(t)$ term of Eq. (1) is called CA, however, in this article, we called that operating $J(t)$ with a fixed $T_{E}$ and no $\gamma(t)$ term of Eq. (1) is CA. There are two reasons. One is that, as for calculating annealing processes, operating $T_{E}$ with a fixed $J$ and no $\gamma(t)$ term is mathematically the same with operating $J(t)$ with a fixed $T_{E}$ and no $\gamma(t)$ term. One is that operating $J(t)$ with a fixed $T_{E}$ and no $\gamma(t)$ term is more practical in relation to the D-Wave chip[10]. These may be trivial, but we consider that these are also significant for understanding quantum annealing computing.

In this article, results at a low temperature are only shown. There is a problem of whether a temperature, that changes the superiority QA over CA, exists or not. However, investigating each optimized temperatures for each annealings can be more important. These are tasks for the future

The size of the transverse field used in the conventional QA schedule is larger than the size of that used in the present QA schedule, however, as the simulation results, the present QA schedule was more effective than the conventional QA schedule. The reason may be that the relaxation time for disorders by large transverse field is required, although the tunneling effect by the transverse field helps to make the systems escape from local minima of the free energy. Therefore, this means that a thorough investigation for setting of QA schedule is needed. This can be an important task for the future.

Very recently, it is reported in Ref. [24] that, by fast annealing, QA is seen to outperform CA for their benchmark spin-glass problems. On the other hand, our results showed that, by slow annealing, QA is seen to outperform CA for our benchmark spin-glass problems. Since the used annealing schedules and the used models are different, we can not say anything for the differences between the results, but they have also proposed a way of comparing QA with CA. Their schedule uses observables, on the other hand, our schedule uses a smallest effective field. Their schedule can be a better schedule, however their schedule depends on models. On the other hand, our schedule does not depend on models. For this difference between the QA schedules, we believe that the present study is also significant. The detailed comparison is a task for the future.

Other QA schedules using the smallest effective field derived in this article may also be worth considering. The present QA schedule is made based on the concept of perturbation theory. Other QA schedules using the smallest effective field based on the same concept can also be considered, and some schedules among them may be better for performance of optimization.

A set of the computer program source codes, used in this study, is available on payment basis [25].

Acknowledgments

The spin glass server [26] was used to obtain the ground-state energies for realizations. We thank the persons concerned with the spin glass server.