Quantum Monte Carlo Annealing with Multi-Spin Dynamics

TL;DR

This paper presents a new Simulated Quantum Annealing algorithm using multi-spin quantum fluctuations, which improves performance over traditional transverse field methods in large disordered systems like spin glasses.

Contribution

Introduction of a multi-spin driver Hamiltonian in SQA and development of an efficient Quantum Monte Carlo algorithm for disordered systems.

Findings

Multi-spin driver improves annealing performance.

Efficient Quantum Monte Carlo simulation for large disordered systems.

Potential for broader application of multi-spin quantum fluctuations.

Abstract

We introduce a novel Simulated Quantum Annealing (SQA) algorithm which employs a multispin quantum fluctuation operator. At variance with the usual transverse field, short-range two-spin flip interactions are included in the driver Hamiltonian. A Quantum Monte Carlo algorithm, capable of efficiently simulating large disordered systems, is described and tested. A first application to SQA, on a random square lattice Ising spin glass reveals that the multi-spin driver Hamiltonian improves upon the usual transverse field. This work paves the way for more systematic investigations using multi-spin quantum fluctuations on a broader range of problems.

Equations

$$H(t)=H_P+(t_{\mathrm{final}}-t)\text{leavevmode}{H}_Q\text{leavevmode}.$$

$$H_P=\sum _{i,j}{J}_{ij}{\sigma }_i^z{\sigma }_j^z\text{leavevmode},$$

$$H_Q^{TF}=-\Gamma \sum _i{\sigma }_i^x\text{leavevmode},$$

$$H_Q^{FI}=-\Gamma \sum _i{\sigma }_i^x-\Lambda \sum _{i,j}{\sigma }_i^x{\sigma }_j^x\text{leavevmode},$$

$$H=-\left(J\sum _{i,j}{\sigma }_i^z{\sigma }_j^z+\Lambda \sum _{i,j}{\sigma }_i^x{\sigma }_j^x\right)\text{leavevmode}.$$

$$h_k=\sum _{b\in \{b{\}}_k}{H}_{i_b,j_b}=\sum _{b\in \{b{\}}_k}{H}_b\text{leavevmode},$$

$$Z=tr\text{leavevmode}{e}^{-\beta H}=\underset{M\to \infty }{\lim }tr\text{leavevmode}{\left(\prod _{k=1}^K{e}^{-\frac{\beta }{M}{h}_k}\right)}^M.$$

$$Z=\sum _{\{s_{il}^z\}}W(\{s_{il}^z\})=\sum _{\{s_{il}^z\}}\prod _p{W}_p(\{s_p\})\text{leavevmode},$$

$$W_p(s_p)=⟨s(i_b,l)s(j_b,l)|e^{-\Delta H}|s(i_b,l+1)s(j_b,l+1)⟩\text{leavevmode},$$

$$W(1)$$

$$W(2)$$

$$W(3)$$

$$W(4)$$

$$\sum _j{w}_{ij}$$

$$w_{ij}$$

$$w_{12}$$

$$w_{13}$$

$$w_{34}$$

$$w_{12}$$

$$w_{13}$$

$$w_{11}$$

$$H=\sum _{i,j}{J}_{ij}{\sigma }_i^z{\sigma }_j^z-\Gamma \sum _i{\sigma }_i^x-\Lambda \sum _{i,j}{\sigma }_i^x{\sigma }_j^x\text{leavevmode},$$

$$\mathcal{C}(t_{final})=t_{final}\text{leavevmode}\bar{n},$$

$$\Gamma (t)$$

$$\Lambda (t)$$

Full text

I Introduction

Discrete combinatorial optimization problems can be encoded into minimizing the energy of classical Ising-type HamiltoniansFu and Anderson (1986) and it has been proposed that adding quantum mechanics could help in finding the ground state of these frustrated Ising system faster than any classical technique. One can define a suitable time-dependent quantum Hamiltonian that connects an initial Hamiltonian, whose ground state is easy to prepare, to the final classical Ising Hamiltonian. At sufficient low temperatures a quantum system undergoing such adiabatic relaxation reaches its ground state and solves the optimization problem. This technique is called quantum annealing (QA)Farhi et al. (2000); Kadowaki and Nishimori (1998); Farhi et al. (2001); Das and Chakrabarti (2008).

Early numerical studies predicted QA Kadowaki and Nishimori (1998); Farhi et al. (2001); Santoro et al. (2002) to be a competitive computational resource for solving spin glass problems, compared to its closely related classical counterpart, the simulated annealing (SA) algorithm Kirkpatrick et al. (1983). However, so far no quantum speedup has been observed in experiments Boixo et al. (2014); Rønnow et al. (2014). The essential difference between these two heuristic methods relies on the type of fluctuations which drive the system away from the multiple local minima, occuring in rugged energy landscapes.

Quantum fluctuations are expected to give an advantage to quantum annealing in particular when the free energy landscape displays tall but narrow barriers. These are easier to tunnel through quantum-mechanically, compared to climbing over them by means of thermally activated escape events.

In QA the system closely follows the ground state of a time-dependent Hamiltonian $H(t)$ which at $t=0$ is dominated by a pure quantum fluctuation part $H_{Q}$ whereas the final Hamiltonian $H(t_{\rm final})$ encodes only the cost function $H_{P}$ of the combinatorial optimization problem. The Hamiltonian as a function of the time $t$ may read, in the case of simple linear annealing schedules,

[TABLE]

Since random ensembles of hard problems are closely connected to spin glass models, we choose, for our tests, the problem Hamiltonian as an Ising spin glass

[TABLE]

where $\sigma_{i}^{z}$ are Pauli matrices, acting on spins $i$, and $J_{ij}$ is the coupling between spins $i$ and $j$, which are randomly uniformly distributed in a range $[-1,1]$ .

Most QA studies so far employed a single type of the quantum term $H_{Q}$, a transverse-field (TF)

[TABLE]

where the $\sigma_{i}^{x}$ operator acts locally on the spin index $i$ inducing quantum fluctuations. This is simplest to implement both in physical devices, such as the D-Wave devicesJohnson et al. (2011); Harris et al. (2010); Bunyk et al. (2014); Boixo et al. (2014) and in simulated quantum annealing (SQA) by means of quantum Monte Carlo (QMC) methodsSantoro et al. (2002); Heim et al. (2015); Isakov et al. (2016); Denchev et al. (2015).

However, as pointed out by Ref. Suzuki et al., 2007, it is possible to choose also different types of quantum fluctuations operators $H_{Q}$. Following Ref. Suzuki et al., 2007 we define a two body fluctuation operator with ferromagnetic interactions (FI) as

[TABLE]

where i,j are nearest neighbours, so that the interaction is short-range.

The adiabatic theorem states that the system follows the instantaneous ground state as long as the total annealing time $t_{final}\gg 1/\Delta^{2}$ where $\Delta$ is the minimum energy gap from the ground state, which is encountered along the annealing run of Eq. (1). It has been conjectured and shown for simple models, that employing quantum fluctuations beyond the TF term could be beneficial to avoid these small gap events Suzuki et al. (2007); Seoane and Nishimori (2012), which are the bottlenecks of QASantoro et al. (2002); Altshuler et al. (2010); Farhi et al. (2012); Knysh (2016). Another argument in support of this possibility is that, in Eq. (4), we are effectively adding a new schedule parameter $\Lambda$ which can be also optimized to increase the QA performance.

Moreover, it has been found that quantum annealing with transverse fails to identify all degenerate ground-state configurations, preventing a fair sampling of equally probable statesMatsuda et al. (2009); Mandrà et al. (2016). Multi-spin quantum fluctuations could alleviate or possibly remove this issue.

In this paper we devise a QMC algorithm to perform SQA with this two-spin transverse ferromagnetic interaction. We note that the minus sign in front of $\Lambda$ makes this Hamiltonian stoquastic for $\Lambda>0$. This means that it is sign-problem free and can thus be simulated by QMC methods. On the other hand, the existence of a sign problem precludes the possibility to simulate the two spin transverse antiferromagnetic ($\Lambda<0$) interactionSeki and Nishimori (2012); Hormozi et al. (2016); Nishimori (2016) within the same approach unless the graph of couplings is bipartite.

The paper is organized as follows, in Sect. II we describe the QMC algorithm, which is tested against exact results for small systems in Sect. III. In Sect. IV we report a first illustrative application to the relevant problem of Ising glass SQA. We discuss the results as well as possible future developments in Sect. VI.

II Quantum Monte Carlo Methods

II.1 Overview

Quantum Monte Carlo (QMC) methods are the only classical approaches for simulating quantum annealing on systems as large as the ones realized experimentally, which consider more than $N=1000$ qubits. Other numerical methods such as unitary evolutionKadowaki and Nishimori (1998); Farhi et al. (2001) scale exponentially with $N$ and are therefore limited to much smaller numbers of spins ($N\approx 20$).

The most widely used QMC technique in the context of QA is Path Integral Monte Carlo (PIMC). PIMC relies on the path integral formalism of quantum mechanics and samples the density matrix corresponding to the quantum Hamiltonian $H$ by means of a classical Hamiltonian $H_{cl}$ on an extended system having an additional dimension, the imaginary time directionSantoro et al. (2002); Heim et al. (2015); Isakov et al. (2016). The original transverse field quantum spin system is then mapped into a classical one, which can be simulated by standard Metropolis Monte Carlo.

It has recently been shown Isakov et al. (2016); Jiang et al. (2017) that for tunneling through a barrier the PIMC efficiency scales as a function of the system size as a physical QA would . This connection exists due to the fact that the tunneling rate in QMC and in the exact real-time quantum evolution scales in the same way, to leading order, as a function of the relevant parameters.

This connection has two implications. The first is that a QMC annealing simulation is relevant for investigating the behaviour of a real quantum annealer, as long as it is stoquastic and the QMC simulations thus do not suffer from a sign problem. The second is that it makes QMC-SQA a competitive tool compared to QA in real devices, by precluding a scaling advantage of stoquastic QA over SQA (cfn. also Ref. Denchev et al., 2015).

Starting from this mapping, other quantum inspired algorithms can be developed such as PIGS, where open boundary conditions in imaginary time are usedIsakov et al. (2016) for further acceleration. In this case QMC represents only a classical optimization algorithm and the physical meaning of the simulation is lost. Similarly, it has been recently shown that performing PIMC away from the converged physical limit, provides a more efficient algorithm compared to the continuos time oneHeim et al. (2015). This can be also considered a quantum inspired algorithm. For this reason we focus here on implementing a discrete-time PIMC algorithm rather than a continuos time PIMC.

II.2 Local versus Cluster Updates

The simplest Metropolis algorithm performs local updates: one generates trial configurations by simply flipping one spin at a time. While this algorithm is ergodic, it can be very inefficient. For non-frustrated systems considerable improvements is obtained by cluster algorithmsSwendsen and Wang (1987); Wolff (1989). They allow simultaneous flips of clusters of spins and are especially advantegeous when large scale fluctuations are important, usually around phase transitions, where local update algorithms slow down.

Standard cluster algorithmsSwendsen and Wang (1987); Wolff (1989) allow the clusters to grow without restrictions. While efficient for non-frustrated systems this approach breaks down for frustrated spin systems, as spin glasses. Due to the frustration the clusters grow to fill a very large fraction of the lattice, if not the entire lattice. The clusters are then larger than the physical domains that should be flipped and the algorithm becomes inefficientKandel et al. (1992); Kandel and Domany (1991); Coddington and Han (1994); Zhu et al. (2015). Essentially the whole system freezes out and one flips almost all spins.

Therefore, in the context of SQA of spin glasses with transverse field an hybrid approach is usually employed, in which restricted clusters are builtRieger and Kawashima (1999); Heim et al. (2015). Here, local clusters can grow on a single site only in the imaginary time dimension, where there is no frustration. In this paper we employ the same strategy with an important modification. Due to the non local nature of the two-spin interaction as cannot restrict ourselves to strictly local clusters but allow a small extent in the spatial direction.

II.3 Cluster QMC for Pure Two-Spin Couplings

We start with the case $\Gamma=0$, i.e. considering only the two-spin term in Eq. (4). The Hamiltonian $H=H_{P}+(t_{\rm final}-t)\leavevmode\nobreak\ H_{Q}$, at any time $t$ is as a special case of the anisotropic quantum XYZ model Takahashi and Suzuki (1972). In this Section we present a modified version of the loop algorithm Evertz (2003) used to simulate the thermodynamics of this model. Readers unfamiliar with the loop algorithm arte referred to the review Evertz (2003) for details.

For sake of simplicity we consider a ferromagnetic Hamiltonian with homogeneous couplings $J>0$ in this discussion, defined on an arbitrary graph of $N$ sites connected by $B$ bonds (the disordered case $J_{ij}$ will be considered at the end):

[TABLE]

The generalization of the algorithm to spatially varying couplings is straightforward.

We first identify the non-commuting pieces of this Hamiltonian to perform a Trotter breakup. This procedure is less straightforward compared to the transverse field case where the Hamiltonian splitting was a trivial splitting $H=H_{P}+H_{Q}^{TF}$. Here we need to split the Hamiltonian $H=\sum_{k=1}^{K}h_{k}$ into a sum of commuting bond Hamiltonians $h_{k}$

[TABLE]

where $i_{b},j_{b}$ are two sites connected by the bond $b$ and $H_{b}=-(J\sigma_{i_{b}}^{z}\sigma_{j_{b}}^{z}+\Lambda\sigma_{i_{b}}^{x}\sigma_{j_{b}}^{x}).$ Each set $\{b\}_{k}$ is defined such that its elements don’t share any site, i.e. no site $i$ appears twice (see Fig. 1.a) . The number $K$ of non-commuting terms $h_{k}$ depends on the graph’s connectivity. For example, in a linear chain we have $K=2$ and the standard checkerboard decompositionEvertz (2003), whereas for a square lattice $K=4$.

The partition function then reads

[TABLE]

Inserting complete sets of $\sigma^{z}$ eigenstates, and following Ref. Evertz, 2003 we write

[TABLE]

where the outer summation is carried over all possible spin configuration on the extended lattice. The index $i$ runs over the original graph sites $i=1,\cdots,N$, whereas $l$ label the imaginary time coordinate $l=1,\cdots,MK$. The index $p$ extends over all the shaded plaquettes of the extended lattice (cfn. Ref. Evertz, 2003). We define a plaquette as the 4-spin configuration $p=p(b,l)=\{(i_{b},l),(j_{b},l),(i_{b},l+1),(j_{b},l+1)$. The statistical weight of each plaquette $p=p(b,l)$ is

[TABLE]

where $\Delta=\beta/M$. $W_{p}(s_{p})$ is different from $1$ only when the piece of Hamiltonian $h_{k}$ acts on the bond $b$ at the correct time $l=(\ell-1)K+k$, with $\ell=1,\cdots,M$. This condition defines the shaded plaquettes. For each bond $b$ we have exactly $M$ shaded plaquettes, one for each Trotter time step $\ell$ (see Fig. 1.b) . The non-trivial Hamiltonian evolution in imaginary time occurs only on these special 4-spins configurations. For the simple case of the spin chain, the shaded plaquettes construction resembles a checkerboard lattice. There are only eight non vanishing matrix elements on the shaded plaquettes, which can be divided in four types:

[TABLE]

expanding the exponential for small $\Delta$, and taking into account that, for example $W(1)=\langle++|e^{-\Delta H_{b}}|++\rangle=\langle--|e^{-\Delta H_{b}}|--\rangle$ since the one body magnetic field operator $\sigma_{z}$ is absent. Since $H_{b}$ preserves the parity of the magnetization, only 8 out of the 16 possible states give a non-zero contribution. Notice that this property would not hold in the transverse field case.

The simplest Monte Carlo algorithm usually employs local updates, i.e. by proposing to flip one spin at a time. This is not a good choice for XX couplings as unconstrained single spin flips lead to forbidden plaquette configurations with zero weight. Our strategy is therefore to implement directly a cluster MC algorithm which automatically avoids sampling of forbidden configurations and minimizes the autocorrelation times. The key feature of the loop algorithmEvertz (2003) is that it allows nonlocal changes of spins configurations, by making only local stochastic decisions on the shaded plaquettes. We refer the reader to the original Ref. Evertz, 2003 for an exhaustive and rigorous justification of the algorithm. Here we simply describe the outline of our optimized version for XX couplings, which, as it usually happens in most MC algorithms, can be divided in two steps: cluster formation and cluster flip.

II.3.1 Outline of the Algorithm

We first restrict the growth of the cluster to a subset of bonds $\mathcal{B}_{m}$ such that $\bigcup_{m}\mathcal{B}_{m}=\mathcal{B}$ contains all the $B$ bonds of the graph (see Fig. 2). For example, a possible choice could be these set of bonds which constitute all the possible smallest loops in the graph, as in Sect. IV.2. This is necessary to keep the cluster size small in the $N\rightarrow\infty$ limit (see Sect. II.2).

Next we explicitely show how to combine to standard Loop algorithmEvertz (2003) with the restricted update scheme, therefore in the following we implicitely refer to the Loop algorithm terminology and we suggest that the reader should become familiar first with the standard version described in Ref. Evertz, 2003.

The algorithm proceeds as follows:

i. We randomly select the set $\mathcal{B}_{m}$ of spins on which we want to build a cluster. Then we choose one bond $b_{O}\in\mathcal{B}_{m}$ and a Trotter time-slice $\ell=1,\cdots,M$. This identifies the corresponding shaded plaquette $p=p(b_{O},l_{O})$. Finally we randomly select one of the two sites, i.e. $(i_{O},l_{O})$, connected by the bond $b_{O}$, at the imaginary time $l_{O}$. This will be the first site added to the cluster

ii. We follow the Loop Algorithm rules to add spins to the cluster. In this case we add only spins which are connected by a bond in the set $\mathcal{B}_{m}$.

iii. We propose to flip this cluster with suitable probability $p_{flip}$, i.e. we reverse each spin $s\rightarrow-s$ belonging to the cluster. $p_{flip}=\mathrm{min}[1,\mathrm{exp}(\delta E\leavevmode\nobreak\ K/\Delta)]$ is given by the change in energy $\delta E$ due to the cluster flip. Notice that only the bonds outside the set $\mathcal{B}_{m}$ will be taken into account in computing $\delta E$ since the cluster building rules take care of detailed balance for all bonds within $\mathcal{B}_{m}$. In this way we evaluate the energy gain of the cluster flip restricted to $\mathcal{B}_{m}$, which is embedded in the graph. This is the main difference compared to Ref. Evertz, 2003.

Notice also that if $\mathcal{B}_{m}=\mathcal{B}$ already contains all the possible edges of the graph, then step iii. simplifies to flipping the cluster with probability $1$, and we recover the single cluster formulation of the Loop Algorithm. We label this type of move as unrestricted or global update. If instead the cluster is restricted to be on a localized set of bonds, we are performing a semi-local update.

II.3.2 Breakup Probabilities

The essential ingredients of the Loop Algorithm are the so-called “breakups weights”. In short, this set of weights $w_{ij}$ ($i,j=1,2,3,4$) determines the shape and the size of the clusters. For example, the weight $w_{ij}$ sets the probability that a shaded plaquette of type $i$ is changed into a plaquette type $j$, after the accepted MC update [cfn. Eq. (10)]. In particular, this probability is given by $p_{i\rightarrow j}=w_{ij}/W(i)$. For example, if we connect and flip the spins at the bottom of a type-1 plaquette, we obtain a type-4 plaquette. The flip of two spins inside a shaded plaquette always produce a different type of configuration, which is in turn a valid plaquette configuration, by construction. This is the main idea of the Loop algorithm. One needs to define every possible $i\rightarrow j$ transition weight. Some of these weights can also be zero, this means that some particular transitions are not-allowed.

We refer the interested reader to the original Loop AlgorithmEvertz (2003) for details, derivations and discussions. Here we simply provide our choice for the breakups weights for future reproducibility.

While the plaquette weights $W(i)$, with $i=1,2,3,4$, are fixed by the Hamiltonian [cfn. Eq. (10)], there is freedom in choosing the breakup weights $w_{ij}$, provided that the following constraints are met

[TABLE]

Experience showed that one should optimize these weights for an efficient algorithm. Indeed, the autocorrelation time dramatically increases with non-optimal choices of $w_{ii}$, which gives the probability to remain in the same plaquette configuration. We thus aim to minimize this weight. The optimal weight’s set varies as a function of the annealing schedule. Let us assume that $\Delta$ is small so that $W(1),W(2)>W(3),W(4)>0$.

Let us consider first the case $\Lambda>|J|$, i.e. at the beginning of the annealing. Under this condition it is possible to set all the freezings to zero. The non-zero weights are:

[TABLE]

with $w_{ij}=w_{ji},\leavevmode\nobreak\ i,j=1,2,3,4$. In the later stage of the annealing we have $\Lambda<|J|$ instead, so we find

[TABLE]

In this case is not possible to have $w_{ii}=0$, but this is the optimal choice as pointed out also in Ref. Kawashima, 1996. Notice that these two sets always satisfy Eq. (14) and work both for ferromagnetic and anti-ferromagnetic couplings.

II.4 Adding the Transverse Field

To obtain an algorithm for the full Hamiltonian

[TABLE]

we need to include the transverse field term $\sigma_{i}^{x}$ to the Loop Algorithm ($\Gamma=0$) described in Sect. II.3. This extension is hinted in Ref. Evertz, 2003 and is not trivial as the transverse field term $\sigma_{i}^{x}$ breaks the parity of the magnetization within each interaction plaquettes. Instead of enlarging the total number of allowed plaquette configurations, taking into account all these possible states, we treat this operator stochastically by adding additional single-site breakups. The transverse field operator can end a loop cluster at any spacetime point $(i,l)$, whereever $\sigma_{i}^{x}$ acts. In our algorithm we add the possibility to stop the cluster when entering into a shaded plaquette at site $i$, with probability given by $p_{x}=\sinh{(\Delta\Gamma/K_{i})}$, where $K_{i}$ denotes the number of physical neighbours of the site $i$ in the graph. This is the generalization to arbitrary graphs of the procedure sketched in Ref. Evertz, 2003 for a linear chain In the following we give details on the actual implementation of this idea.

i. Suppose we start the cluster from position $(i_{0},l_{0})$.

ii. When we jump into a shaded plaquette, say at position $(i,l)$, if the plaquette type is not of type 3 or 4111For simplicity we don’t break clusters on these types of plaquettes. This approximation is justified in the $M\rightarrow\infty$ limit., we stop the cluster with probability $p_{x}$.

iii. Each time we stop we keep track of this position by inserting a $\sigma_{i}^{x}$ operator label. We put this label on the bond above (below) site $i$, if the direction was upward (downward) in imaginary time direction. Otherwise we continue to build the cluster following the procedure described in Sect.II.3.1.

iv. If we stop the cluster then we restart from $(i_{0},l_{0})$ and proceed in the opposite direction until we stop again, notice that now we can stop either by inserting a new $\sigma_{i}^{x}$ operator or by touching an existing one already in place.

v. The cluster flipping decision remains unchanged.

vi. Finally, remove the $\sigma_{i}^{x}$ operator labels between segments of equal orientation. This completes one update.

The last issue concerns the existence of plaquettes having an odd number of spin up (down), which fall outside the breakup selection rule of the XX loop algorithm. These plaquettes certainly occur if $\Gamma>0$. Since in our approximation we can encounter the $\sigma_{i}^{x}$ operator only going vertically along the imaginary time direction, the decision rule can be adapted from the standard transverse field cluster algorithmRieger and Kawashima (1999) in this case: suppose we are at site $(i,l)$ and we are proceeding upward, if there is not a $\sigma_{i}^{x}$ operator already in place above $(i,l)$, then we check whether the spin at $(i,l+1)$ is parallel to $(i,l)$, or not . In the latter case we stop the cluster.

Finally we notice that, in the $\Lambda=0$ case, this algorithm does not reduce to the common TF algorithm for disordered systemsRieger and Kawashima (1999), in which clusters are built only along the imaginary time direction, if a non-empty bond set $\mathcal{B}_{m}$ is considered. Indeed, the cluster can still span the entire bond region $\mathcal{B}_{m}$, due to the occurence of the freezing plaquettes, which make the cluster non-local. We use therefore the same cluster algorithm, within the same choice of $\mathcal{B}_{m}$’s, to compare FI and TF Hamiltonians, also considering the limiting $\Gamma=0$ case.

III Phase Diagram of the XX Model with a Transverse Field

We test the accuracy of the algorithm against exact diagonalization (ED) resultsBauer et al. (2011), for a small uniform ferromagnetic chain and square lattices with spacial periodic boundary conditions, having Hamiltonian given by Eq. (22), i.e. without bond disorder. In Fig. 3 we plot the $\langle\sigma^{z}_{i}\sigma^{z}_{j}\rangle$ correlation function, for nearest neighbours, as a function of the inverse temperature $\beta$, for several choices of the parameters $\Lambda$ and $\Gamma$. The agreement of QMC with ED is satisfactory in the full temperature range.

Next, in Fig. 4, we focus on the low-temperature regime relevant for SQA and benchmark simulations obtained with different combinations of $\Lambda$ and $\Gamma$ values, which can be realized along a generic SQA run.

Finally in Fig. 5, we show, as an example, the low-temperature phase diagram of the Hamiltonian in Eq. (22), with $J_{ij}=1$, on a 2D square lattice. Generalizing the breakup weights in Sect. (II.3.2) for an arbitrary XYZ Hamiltonian, it would be possible to study the phase diagram of sign-problem free XYZ models with transverse field on arbitrary lattice.

IV Simulated Quantum Annealing

IV.1 Annealing schedule

We now report a first application of the algorithm to SQA, to assess the annealing sensitivity to the $H^{FI}_{Q}$ driver Hamiltonian. We consider a spin glass problem Hamiltonian, defined on a $10\times 10$ square lattice, with periodic boundary conditions, and uniformly randomly distributed couplings in the range $[-1,1]$. Following previous SQA studiesHeim et al. (2015); Santoro et al. (2002), we plot the residual energy $E_{res}(t_{final})=E(t_{final})-E_{0}$ as a function of different annealing times $t_{final}$, where $E_{0}$ is the exact solutionspi and $E(t_{final})$ is the ground state of the Hamiltonian found at the end of the annealing.

In order to compare the performance of different annealing strategies in the $\Gamma-\Lambda$ parameters space, we need to rigorously define the computational effort for a QMC algorithm having unrestricted and semi-local updates. Different annealing schedules lead to different average cluster sizes, and, the larger is the cluster built, the heavier is the computational cost of each update. This cost is proportional to the cluster size $\bar{n}$ and to the number of bonds one has to check for evaluate $\delta E$ in the acceptance/rejection step. We notice that, in the algorithm, the computationally expensive operations are made on the shaded plaquettes. Therefore, the total effort is still proportional to the number of sites $N$ (or generically to the bond number $B$) and the number of Trotter slices $M$, although the total number of slices along the imaginary time direction is $M\times K$.

Here we thus define the computational cost associated we each annealing run as

[TABLE]

where $t_{final}$ is the total number of Monte Carlo updates.

We notice that the QMC algorithm is efficient, for low $\Lambda$ and in the thermodynamic limit, only if we restrict the cluster growth, at each update, to a small bond subset with $\#\mathcal{B}_{m}\ll B$. In this way $\bar{n}\ll N\times MK$

Usually, with TF-SQA, a linear schedule with starting field $\Gamma_{0}=\Gamma(t=0)>\mathrm{max}\leavevmode\nobreak\ |J_{ij}|$ is employedHeim et al. (2015), in our notations, this would correspond to decrease linearly the transverse field parameter as $\Gamma(t)=\Gamma_{0}(1-t)/t_{final}$ and set $\Lambda=0$, $\forall t$. In the following we will compare to the simplest annealing path for SQA with ferromagnetic interactions (FI-SQA), given by

[TABLE]

namely, decreasing linearly both the control parameters. We will also compare to the choice $\Lambda(t)=\Lambda_{0}(1-t)/t_{final}$ and $\Gamma=0$, $\forall t$, i.e. employing only the two pure two-body operator. The optimization of the annealing path to non-trivial $\Lambda(t),\Gamma(t)$ time-dependences is left for future studies.

IV.2 2D Spin Glass Results

In Fig. 6 we show the residual energy as function of the annealing time, for $10\times 10$ lattice, and using the semi-local loop update, an approach that remains efficient for large disordered systems. In this case, each possible bond subset $\mathcal{B}_{m}$, with $m=1,\cdots,N$, is defined to be the smallest 4-bonds loop that can be constructed on the square lattice (see Fig. 2). The four lattice sites which belong to these sets are given by $(i,j),(i+1,j),(i+1,j+1),(i,j+1)$, with $i,j=1,\cdots,L$.

With this choice the QMC algorithm is ergodic as the original loop algorithm (with global updates) and maintains efficiency in the disordered case, in the $N\rightarrow\infty$ limit. In the general case, each bond subset $\mathcal{B}_{m}$ has to be provided as an input, and varies with the graph under consideration.

The bond subsets defining each non-commuting Hamiltonian $h_{k},\leavevmode\nobreak\ k=1,\cdots,K$ is also an input of the algorithm. Generally, performing this decomposition can be also an hard optimization task, which is related to the edge coloring problem, but it can always be solved by using a sub-optimal $K$, one larger than the minimal value (which depends on the graph complexity)Misra and Gries (1992). For regular graphs, such as the square lattice, these sets arte always easy to construct.

In this study we compare TF and FI-SQA as classical optimization algorithms. Therefore we renormalized the annealing time in such a way to fairly compare the theoretical computational effort of each run, as discussed above. We also define the residual energy as the minimum energy $E_{min}$ among all the possible $M$ Trotter slices. We perform the simulations at different inverse temperatures $\beta=20$, 40, and 80. We keep the Trotter time-step $\beta/M=0.3125$, constant and close to convergence to the physical limit (cfn. also Fig. 8).

In Fig. 6 we observe that, for fixed $M$ and $\beta$, TF and FI-SQA display different behaviour, as the two curves are not related by a trivial shift along the time axis. We notice that the FI driver always outperforms the standard TF, for sufficiently large annealing times, at each temperature, despite the larger complexity of the FI driver. Interestingly, the lower is temperature, the larger is the difference between the FI and TF residual energies.

Surprisingly, performing the annealing at zero transverse field, $\Gamma=0$, results in a drastic decrease of the performance. These observations can be explained in the following way: i. For short annealing time, a large part of the residual energy can be recovered by simple one-flip moves, then the TF driver is more effective in eliminating the defects in the extended lattice. ii. The pure two body driver Hamiltonian $\sigma_{i}^{x}\sigma_{j}^{x}$ is way less efficient in this regards, as the clusters are usually bigger and always closed in the extended lattice. Therefore, the transverse-field operator is still an important ingredient for a QA driver Hamiltonian.

Finally, we also check the performance using the global update algorithm. In this case the cluster can freely traverse the extended lattice. We consider the same system and perform the same set of simulations and show the results in Fig. 7. We notice that the trend is preserved, although the annealing profiles are slightly different. In particular, the pure two-body driver performs significantly worse than the mixed FI-SQA, having both $\Lambda$ and $\Gamma$ non-zero. However, FI-SQA still outperforms TF-SQA for sufficiently long annealing times, as in the previous semi-local algorithm. Interestingly, this global updates version of the algorithm displays a better efficiency compared to the previous approach. Indeed such system size can be considered still small and the global update scheme does not suffer much from a critical slowing down.

Note though the non-monotonicity as a function of the annealing time compared to the TF (cfn. also Ref. Heim et al., 2015), which hints at inefficiencies of the global updates. This feature disappears when the average energy is used (see Fig. 8), instead of the lowest one among the Trotter slices.

V Efficiency of Quantum Annealing

In this last section we perform simulations indicative of the relative efficiency between the FI and TF operators in a real QA device instead of as a classical optimization algorithm. QMC is expected to reproduce the performance of physical QA for typical tunneling problems as discussed in Refs. Isakov et al., 2016; Jiang et al., 2017. This time we take the continuous time limit by using a large enough number of time slices $M$, and we measure the annealing time simply using the Monte Carlo steps, as the computational effort related to the different update scheme is not relevant for the study of the physical machine. Moreover, following Ref. Heim et al., 2015 we average the final energy over all the $M$ Trotter slices. In Fig. 8 we observe the same trend as in SQA, where the FI operator eventually outperforms the TF at later stages of the annealing.

VI Conclusions

We introduced a novel QMC algorithm to perform SQA with a transverse ferromagnetic two-spin driver Hamiltonian Eq. (4). This type of fluctuation extends the standard TF-SQA, adding also a two-body transverse operator of the form $\sigma_{i}^{x}\sigma_{j}^{x}$. Though this possibility has been introduced already some years agoSuzuki et al. (2007), it has never been used before with QMC simulations, and therefore never applied so far on large optimization problems. We notice that, since QMC techniques are limited by the so-called sign problem, we can simulate only the ferromagnetic version of two body transverse interaction.

Our discrete-time path-integral Monte Carlo algorithm is an extension of the well established Loop AlgorithmEvertz (2003), with the inclusion of the transverse field operator, and implements restricted cluster updates to simulate efficiently large disordered systems.

A first application to quantum annealing, on a random square lattice, reveals that the new driver Hamiltonian improves upon the usual transverse field, though more systematic studies are required to address conclusively this question. Indeed it will be interesting to optimize the interplay of the two parameters $\Gamma(t)$ and $\Lambda(t)$ along the annealing schedule and to perform a size scaling analysis. Morover it will be important to find other classes of problem Hamiltonians which can benefit more from this technique.

We notice that the range of applicability of the present algorithm can go beyond SQA, since it would be also possible to explore phase diagrams of XYZ sign-problem free Hamiltonians with transverse field, defined on arbitrary lattices, by performing equilibrium simulations with a converged number of Trotter time steps.

Acknowledgements.

This work has been supported by the Swiss National Science Foundation through the National Competence Center in Research QSIT, by the European Research Council through ERC Advanced Grant SIMCOFE and by ODNI, IARPA via MIT Lincoln Laboratory Air Force Contract No. FA8721-05-C-0002. We acknowledge useful discussions with G. E. Santoro, H. Katzgraber, D. Herr and M. Dolfi.