Off-Diagonal Expansion Quantum Monte Carlo

1.2 MB·

by  Tameem Albash

Off-Diagonal Expansion Quantum Monte Carlo
Tameem Albash,1,2 Gene Wagenbreth,3 and Itay Hen1,2,∗
1Information Sciences Institute, University of Southern California, Marina del Rey, California 90292, USA
2Department of Physics and Astronomy and Center for Quantum Information Science & Technology,
University of Southern California, Los Angeles, California 90089, USA
3Cray, Seattle, Washington 98164, USA
(Dated: January 9, 2017)
We propose a Monte Carlo algorithm designed to simulate quantum as well as classical systems
7
at equilibrium, bridging the algorithmic gap between quantum and classical thermal simulation
1
algorithms. The method is based on a novel decomposition of the quantum partition function that
0
2 canbeviewedasaseriesexpansionaboutitsclassicalpart. Wearguethatthealgorithmisoptimally
suitedtotacklequantummany-bodysystemsthatexhibitarangeofbehaviorsfrom‘fully-quantum’
n
to ‘fully-classical’, in contrast to many existing methods. We demonstrate the advantages of the
a
techniquebycomparingitagainstexistingschemes. Wealsoillustratehowourmethodallowsforthe
J
unification of quantum and classical thermal parallel tempering techniques into a single algorithm
5 and discuss its practical significance.

]
h
c I. INTRODUCTION type will have wide-ranging applicability in diverse
e areasrangingfromstatisticalphysicsthroughquan-
m Quantum Monte Carlo (QMC) algorithms are tumchemistrytoquantumcomputing,tomentiona
- known to be notoriously inefficient in ‘almost clas- few areas.
t
a sical’ parameter regimes where updates resulting Here, we propose an algorithm that has the algo-
t from thermal fluctuations are expected to be far rithmicflexibilitytosimulateinteractingmany-body
s
. more dominant than those resulting from quantum systemsrangingfromthefully-quantumtothefully-
t
a fluctuations. This is particularly true in models classical. We present a Monte Carlo scheme that
m that can be parametrically tuned from quantum to is based on a novel decomposition of the canonical
- classical regimes, such as the transverse-field Ising quantum partition function into a sum of ‘general-
d model [1–3], the XXZ model [4, 5] or the Bose- ized’ Boltzmann weights and that converges to the
n
Hubbard model [6–8]1. Since QMC methods evolve usual decomposition of the classical partition func-
o
c via configuration updates that are based on quan- tioninthelimitwheretheHamiltonianofthesystem
[ tum fluctuations, the acceptance rates of quantum becomes classical. Based on this unique decomposi-

updates,e.g.,singlespin(orcluster)flipsintheIsing tion,ouralgorithmaimstoimprovetheconvergence
1
system, decrease dramatically in classical regimes, rates of simulated systems for which existing tech-
v
9 often causing QMC algorithms to dramatically slow niques are often inefficient.
9 down or ‘freeze’ (see, e.g., Ref [9]). Within our approach the quantum imaginary-
4 Efficient classical thermal updates are typically time dimension of the algorithm is ‘elastic’, i.e., it
1 hardtoimplementwithintheframeworkofQMCal- can stretch or shrink dynamically depending on the
0
gorithms because these algorithms do not normally strength of the quantum part of the system — the
.
1 converge to classical Monte Carlo algorithms in the off-diagonalportionoftheHamiltonian. Inaddition,
0 limit where the model becomes classical. For this the proposed method does not introduce Trotter-
7 reason, there are almost no algorithms that effi- type errors as in path-integral QMC (PIQMC), a
1
ciently simulate systems that exhibit the full range sourceoferrorsthatnormallyoccursfromaninsuffi-
:
v of behavior from being ‘fully-quantum’ to ‘fully- cientdiscretizationoftheimaginary-timedimension
i classical’. For the successful simulation of systems (over-discretization tends to sharply reduce the ac-
X
exhibiting the above characteristics, it is therefore ceptance rates of the QMC updates). Moreover, in
r
a important to devise Monte Carlo schemes that can the classical limit where off-diagonal terms vanish,
function both as quantum as well as classical algo- our algorithm naturally reduces to a classical ther-
rithms when necessary. Efficient methods of this mal simulation. As we illustrate, the above proper-
ties allow our method to naturally overcome certain
inefficiencies typically encountered by other QMC
algorithms.
∗ [email protected]
The structure of the paper is as follows. In
1 TheBose-Hubbardmodelsexhibitsaquantumphasetran-
sition from a highly delocalized superfluid at one extreme Sec.II,wedescribethedecompositionofthecanoni-
ofitsparameterspacetoaclassicallocalizedMottinsulator calquantumpartitionfunctionintowhatwereferto
attheotherextreme. asgeneralizedBoltzmannweights(GBWs). Wethen
2
proceed in Sec. III to present the basic steps, up- replace the trace operation Tr[·] with the explicit
(cid:80)
dates and measurements of our off-diagonal expan- sum (cid:104)z| ·|z(cid:105) and then expand the exponent in
z
sion (ODE) quantum Monte Carlo algorithm that the partition function in a Taylor series:
buildsontheabovedecomposition. Toillustratethe
practicality of our algorithm, we examine in Sec. IV
simulationsofthetransversefieldIsingmodel,espe- (cid:88)(cid:88)∞ βn
cially inside the spin-glass phase, where the model Z = (cid:104)z|(−H)n|z(cid:105) (2a)
n!
is known to be hard to simulate. We also discuss in
z n=0
thiscontexttheunificationofquantumandclassical (cid:88)(cid:88)∞ βn (cid:88)
parallel tempering. We present some conclusions in = (cid:104)z|(−H +Γ V )n|z(cid:105) (2b)
n! c j
Sec. V. z n=0 j
(cid:88)(cid:88)∞ (cid:88) βn
= (cid:104)z|S |z(cid:105), (2c)
n! n
II. GENERALIZED BOLTZMANN
z n=0{Sn}
WEIGHTS
A. Decomposition of the partition function
whereinthelaststepwehavealsoexpanded(−H)n,
and {S } denotes the set of all sequences of length
n
Themaininsightattheheartofourapproachisa
n composed of products of basic operators H and
c
novel decomposition of the canonical quantum par-
V .
j
tition function which, as we argue, allows for the
development of a QMC algorithm that has certain We proceed by removing all the diagonal Hamil-
advantages over existing methods. Our work builds tonian terms from the sequence (cid:104)z|S |z(cid:105). We do
n
in part on the Stochastic Series Expansion (SSE) so by evaluating their action on the relevant basis
algorithm, a well-known and successful QMC algo- states, leaving only the off-diagonal operators un-
rithm pioneered by Sandvik [10, 11], which—unlike evaluated inside the sequence. At this point, the
thetraditional‘slicing’ofthepartitionfunctioninto partition function can be written as:
Trottersegments—involvesaTaylorseriesexpansion
in the inverse-temperature β = 1/T (in our units
k = 1) of the partition function as was originally
B
suggested by Handscomb [12, 13]. (cid:88)(cid:88)∞ (cid:88) (cid:32)(cid:88)∞ βn(−1)n−q
The canonical quantum partition function of a Z = Γq(cid:104)z|Sq|z(cid:105) n!
system described by a Hamiltonian H is given by z q=0{Sq} n=q
Z = Tr (cid:2)e−βH(cid:3). Our decomposition begins by first 
writing the Hamiltonian in the form × (cid:88) Ek0(z0)·...·Ekq(zq) , (3)
H =H −(cid:88)Γ V . (1) (cid:80)ki=n−q
c j j
j
Here, H is a ‘classical’ Hamiltonian, i.e., a diag-
c
onal operator in some known basis, which we re- where E(zi) = (cid:104)zi|Hc|zi(cid:105) and {Sq} denotes the set
fer to as the computational basis, and whose basis ofallsequencesoflengthq of‘bare’off-diagonal op-
states will be denoted by {|z(cid:105)}. The {Γj} are gen- erators Vj. The term in parenthesis sums over the
erally complex-valued parameters, and {Vj} are off- diagonal contribution of all (cid:104)z|Sn|z(cid:105) terms that cor-
diagonal operators satisfying [Vj,Hc] (cid:54)= 0 that give respondtoasingle(cid:104)z|Sq|z(cid:105)term. Thevarious{|zi(cid:105)}
the system its ‘quantum dimension’. In an analo- states are the states obtained from the action of the
gous way to standard SSE, in order for the decom- orderedVj operatorsinthesequenceSq on|z0(cid:105),then
position of the partition function to be feasible, we on |z1(cid:105), and so forth2. Figure 1 gives a schematic
require the off-diagonal operators to be chosen such representation of (cid:104)z|Sq|z(cid:105). After a change of vari-
that they obey V |z(cid:105)=|z(cid:48)(cid:105) for every basis state |z(cid:105),
j
where |z(cid:48)(cid:105) (cid:54)= |z(cid:105) is also a basis state. For simplic-
itywehenceforthassumethatalltheΓ parameters
j
are identical, namely that Γ = Γ,∀j, however as
j
will become evident shortly this restriction is by no
means necessary.
We now present the main steps for the decompo- 2 For example, for Sq = Viq...Vi2Vi1, we obtain |z0(cid:105) =
sition of the quantum partition function. We first |z(cid:105),Vi1|z0(cid:105)=|z1(cid:105),Vi2|z1(cid:105)=|z2(cid:105),etc.
3
input variables [E ,...,E ]. In our case, F[·] is the
0 q
⟨𝑧|𝑉 𝑉 … 𝑉 |𝑧⟩ function
% % %
& ’ )
F[E ,...,E ]=e−β[E0,...,Eq]. (8)
0 q
A feature of divided differences is that they are in-
𝑧 = 𝑧, 𝑧- 𝑧/ 𝑧.0- 𝑧. = 𝑧 variant under rearrangement of the input values, so
the input sequence forms a multiset, i.e., a general-
ization of the mathematical set which allows repe-
Figure 1. A diagrammatic representation
titions but where order does not play a role. The
of the term (cid:104)z|S |z(cid:105). The sequence of operators
q partition function in its close-to-final form is thus
S =V ·V ···V issandwichedbetweenclassicalbra
q i1 i2 iq given by:
(cid:104)z| and ket |z(cid:105), inducing a sequence of classical states
(|z (cid:105),...,|z (cid:105)). The multiset of classical energies of the
0 q ∞
states|z (cid:105),namely,E =E(z )=(cid:104)z |H |z (cid:105)generatethe (cid:88)(cid:88)(cid:88)
i i i i c i Z = (cid:104)z|S |z(cid:105)(−Γ)qe−β[E0,...,Eq]. (9)
generalized Boltzmann weight. q
z q=0{Sq}
ables, n→n+q, we arrive at:
We note that a single divided-difference term is a
(cid:88)(cid:88)∞ (cid:88) (cid:32) (cid:88)∞ (−β)n sum of an infinite number of terms in the usual
Z = (cid:104)z|S |z(cid:105) (βΓ)q breakdown of SSE. This can be immediately seen in
q (n+q)!
Eq. (3), which relates the standard SSE weight, in-
z q=0{Sq} n=0
 volving sequences of diagonal as well as off-diagonal
× (cid:88) Ek0(z0)···Ekq(zq) . (4) ablognodristh(dmentohtaetdobnylySnin)v,otlovethoeffw-deiiagghotnsaolfbthoendcus3r.rent
(cid:80)ki=n Furthermore, the mean value theorem for divided
differences [14, 15] together with the monotonicity
AbbreviatingE ≡E(z )(notethatthevarious{E }
i i i of the exponential function allows us to write
are functions of the |z (cid:105) states created by the opera-
i
torsequenceSq),thepartitionfunctionisnowgiven dq(cid:0)e−βE(cid:1)(cid:12)(cid:12)
by: e−β[E0,...,Eq] = (cid:12) (10)
dEq (cid:12)
(cid:12)
(cid:88)∞ (cid:88) E=E(0,...,q)
Z = (−Γ)q (cid:104)z|Sq|z(cid:105) (5) (−β)qe−βE(0,...,q)
= ,
q=0 z,{Sq} q!
 
× (∞(cid:88),...,∞) (q+(−(cid:80)β)qk )! (cid:89)q (−βEj)kj . fEor ∈a(min[Esin,g.l.e.,E ],mreaaxl-[vEalu,e.d..,E ])energy
i (0,...,q) 0 q 0 q
{ki}=(0,...,0) j=0 calculated from the multiset [E ,...,E ] and β.
0 q
A feature of the above infinite sum is that the term Specifically, E(0,...,q) lies within the spectrum of
in parentheses can be further simplified to give the the classical Hamiltonian. This allows us to write
exponent of divided differences of the E ’s (we give the partition function in terms of classical ‘effective
i
a short description of divided differences and an energies’ as
accompanying proof of the above assertion in Ap-
pendixA),namelyitcanbesuccinctlyrewrittenas: Z =(cid:88)(cid:88)(cid:88) (βqΓ!)q(cid:104)z|Sq|z(cid:105)e−βE(0,...,q). (11)
(cid:88) (−β)q (cid:89)q {z} q {Sq}
(q+(cid:80)k )! (−βEj)kj =e−β[E0,...,Eq],(6)
i
{ki} j=0
where [E0,...,Eq] is a multiset of energies and 3 We note here however that the computational cost of cal-
where a function F[·] of a multiset of input values is culatingtheweightsofconfigurationshereishigherthanin
defined by standardSSE,duetotheneedtoevaluatethedivideddif-
ferencesoftheexponentialfunctionovertheenergiesofthe
q configuration. Asweshowbelow,thecostofeachsucheval-
F[E ,...,E ]≡(cid:88) F(Ej) (7) uationcanbeshowntobeproportionalintheworstcaseto
0 q j=0 (cid:81)k(cid:54)=j(Ej −Ek) tdhireecstqudaivreidoefd-tdhieffenruemncbeesrcoaflctuelramtisonin,steheeAspeqpuenendcixeASqa(nbdy,
e.g., Ref. [14]), which scales linearly with the inverse tem-
and is called the divided differences [14, 15] of the peratureβ andwiththenumberofparticlesinthesystem
function F[·] with respect to the list of real-valued N (thisisdiscussedindetaillateron).
4
To calculate E , one may use the divided dif-
(0,...,q)
ferences recursion relation (see Appendix A)
F[E ,...,E ] (12)
i i+j
F[E ,...,E ]−F[E ,...,E ]
= i+1 i+j i i+j−1 ,
E −E
i+j i
whichintermsofeffectiveclassicalenergiesbecomes
(−β)q
e−βE(0,...,q) =
q!
(−β)q−1(cid:0)e−βE(0,...,q−1) −e−βE(1,...,q)(cid:1)
.(13)
(q−1)! E −E
0 q
Isolating E , we arrive at
(0,...,q)
1 2qsinhβ∆E
E =E¯− log , (14)
(0,...,q) β β(E −E )
q 0
where
2E¯ =E +E and
(1,...,q) (0,...,q−1)
2∆E =E −E .
(1,...,q) (0,...,q−1)
In the limiting case where all energies in the se-
quence are equal, the above relation neatly becomes Figure 2. Calculating the effective classical en-
E = E = E . The initial condition for the ergy of the generalized Boltzmann weight using
(0,...,q) (0) 0
a ‘pyramid’ structure. The evaluation of the divided
above recursion is simply e−βE(i) =e−βEi. An illus-
differencesoftheexponentialfunctionofq+1inputen-
tration of how the recursion relation is used to cal-
ergies consists of calculating each level of the pyramid
culate E is depicted in Fig. 2 and is discussed
(0,...,q) starting at its base. The values at the base of the pyra-
in more detail in Appendix B. midE aresimply theenergyinputsE (shown asthe
(j) j
Finally, since by construction the term (cid:104)z|Sq|z(cid:105) bluelineatthebottomofthepyramid),withallidentical
evaluates to either 0 or to 1 (the operation Sq|z(cid:105) energiesplacedtogetherasagroup(theexactorderingof
returns a basis state |z(cid:48)(cid:105) and therefore (cid:104)z|S |z(cid:105) = theenergiesisnotimportant). Tocalculatetheelements
q
(cid:104)z|z(cid:48)(cid:105)=δ ), the partition function can be rewrit- at the next level of the pyramid, we use the relation in
z,z(cid:48)
teninitsfinalformasasumoveronlynon-vanishing Eq.(14). Thisprocedureiscontinueduntilthefinallevel
ofthepyramidisevaluated,whichcorrespondstotheef-
terms:
fective classical energy in the GBW, namely, E .
(0,...,q)
(cid:88) (βΓ)q
Z = e−βE(0,...,q). (15)
q!
{Sq:(cid:104)z|Sq|z(cid:105)=1} positive, i.e., if the off-diagonal elements are non-
positive, which is the case for the so-called stoquas-
We interpret the terms in the sum in Eq. (15) as
(cid:80) ticHamiltonians[17,18]. Asisalsoevidentfromthe
weights, i.e., Z = W , where the set of config-
{C} C aboveexpression,evenvaluesofq alsoyieldpositive
urations {C} is all the distinct pairs {|z(cid:105),S }. Be-
q weightsregardlessofthesignofΓ. Thiscorresponds
cause of the form of W ,
C to a scenario where off-diagonal operators must be
(βΓ)q injectedalongtheimaginarytimedimensioninpairs
WC = q! e−βE(0,...,q), (16) in order to ensure nonzero weights. One such exam-
ple is the transverse-field Ising Hamiltonian
we refer to it as a ‘generalized Boltzmann weight’
(cid:88) (cid:88) (cid:88)
(or,GBW).WeshallrefertoE asthe‘effective H = J σzσz+ h σz−Γ σx, (17)
(0,...,q) ij i j j j j
classicalenergy’oftheconfigurationC anddenoteit (cid:104)i,j(cid:105) j j
at times for brevity simply by E .
In order to interpret the WCC terms as actual where Hc =(cid:80)(cid:104)i,j(cid:105)Jijσizσjz+(cid:80)jhjσjz and the off-
weights, they must be non-negative for any simu- diagonal operators are the spin-flip terms V = σx.
j j
lated system that is not plagued by the sign prob- In order for the (cid:104)z|S |z(cid:105) terms to evaluate to one
q
lem [16]. It is therefore interesting to note that rather than to zero, off-diagonal operators must al-
the above weights are automatically positive if Γ is waysbeproducedandannihilatedinpairs,implying
5
that the total sign of the weight, Eq. (16), is posi- Second, for any arbitrary energy shift ∆E of
tive. We have thus established a decomposition of the diagonal part of the Hamiltonian, the following
thecanonicalquantumpartitionfunctionintoasum holds:
of positive-valued weights.
e−β[E0+∆E,...,Eq+∆E] =e−β∆Ee−β[E0,...,Eq]. (20)
This identity reflects the fact that the addition of
B. Properties of the GBWs constants to the simulated Hamiltonian has a triv-
ial effect on the various weights. Specifically, ratios
One property of the above decomposition of the of weights, which in turn determine the acceptance
canonical quantum partition function is that it may rates of the QMC updates, are invariant under the
beviewedasunifyingtheclassicalandquantumpar- above addition of a constant, as they should be.
tition functions. Specifically, it contains as a sub On a more academic note, it is interesting to ob-
sum the classical partition function decomposition serve that the proposed algorithm also has close
of its diagonal part H . Writing the quantum par- relations to continuous-time QMC methods (e.g.,
c
tition function as a series in the ‘quantum strength’ Ref. [22]), via the Hermite-Genocchi formula [15]:
parameterΓ,oneobtainstheclassicalpartitionfunc- (cid:90)
tion as the zeroth term, namely, e−β[E0,...,Eq] = dt0···dtqe−β(E0t0+E1t1+...+Eqtq),
Ω
Z = (cid:88) e−βE(z) (21)
{S0:(cid:104)z|z(cid:105)=1} where ti ≥ 0 and the area of integration Ω is
(βΓ)2 (cid:88) bounded by t0+t1+...+tq from above.
+ e−βE(0,1,2) +... (18)
2
{S2:(cid:104)z|S2|z(cid:105)=1}
C. A simple analytical example
Furthermore, in classical regimes where Γ is zero or
very small, the dominant configurations, i.e., those
As a first illustration, let us consider as a sim-
with highest weights, have no off-diagonal terms,
ple example the transverse-field Ising Hamiltonian
and only the q = 0 terms survives. In this case
where the classical Ising part vanishes, namely,
the typical weights are
where H = 0. In this case the model becomes the
c
trivial system H = −Γ(cid:80) σx. The partition func-
(βΓ)q (cid:12) i i
e−βE(0,...,q)(cid:12) =e−βE(0) =e−βE(z), (19) tion in this special case is decomposed as:
q! (cid:12)q=0
(cid:88) (βΓ)q
Z = , (22)
where E(z) is the classical energy of the spin con-
q!
figuration z. Our decomposition thus automatically {Sq:(cid:104)z|Sq|z(cid:105)=1}
reduces to the usual sum over Boltzmann weight of where the classical energies are E = ··· = E =
0 q
classical Hamiltonians4. 0, corresponding to E = 0. In this case, we
(0,...,q)
The GBW, Eq. (16), also has several attractive have (cid:104)z|S |z(cid:105) = 1 if and only if all σx off-diagonal
q i
propertiesthatmakeitusefulforMonteCarlosimu- operators in S appear an even number of times.
q
lations. First,aswasalreadymentioned,itisstrictly Denoting by N (q) the number of nonzero weights
p
positive for stoquastic systems. This feature auto- foreachvalueof(even)qandevery|z(cid:105),thepartition
matically resolves the ‘diagonal sign problem’ that function can be simplified to
sometimesappearsinotherschemes[19],wherecon-
stantsmustbeaddedtothediagonalbondstorectify Z =2N (cid:88) (βΓ)qN (q). (23)
the problem. Moreover, since the addition of such q! p
q even
constants considerably affects the convergence rate
ofthealgorithm[10,11,19–21],theseconstantsusu- A simple calculation (see Appendix C) reveals
ally have to be optimized for faster convergence. A (cid:18) (cid:19)
1 (cid:88) N
QMC algorithm based on the GBW decomposition Np(q)= 2N k (N −2k)q, (24)
is in this respect parameter-free, a property that is k=0
expected to facilitate computations.
which yields
(cid:88)(βΓ)q (cid:88)(cid:18)N(cid:19)
Z = (N −2k)q (25)
q! k
q k=0
4 This is to be contrasted with other decompositions of the (cid:18) (cid:19)
partitionfunctionwheretheclassicallimitiseitherunnat- =(cid:88) N (cid:88) [βΓ(N −2k)]q .
uralorill-defined. k
k=0 q≥0,even
6
Carrying out the sum over q, we end up with: (a)
⟨𝑧|𝑉 …𝑉 |𝑧⟩ ⟨𝑧′|𝑉 …𝑉 |𝑧′⟩
(cid:18) (cid:19) %& %( %& %(
Z =(cid:88) N cosh[βΓ(N −2k)]=(2coshβΓ)N , (b) 𝑧′ 𝑧′′
k
k=0 ⟨𝑧|𝑉 …𝑉 𝑉 …𝑉 |𝑧⟩ ⟨𝑧|𝑉 …𝑉 𝑉 …𝑉 |𝑧⟩
(26) %& %+ %+,& %( %& %+,& %+ %(
which is the correct result for the partition function (c) 𝑧′ 𝑧
for the non-interacting system H =−Γ(cid:80) σx.
i i ⟨𝑧|𝑉%&…𝑉%+𝑉%+,&…𝑉%(|𝑧⟩ ⟨𝑧′|𝑉%+,&…𝑉%(𝑉%&…𝑉%+|𝑧′⟩
𝑧′ 𝑧′′
(d)
III. OFF-DIAGONAL EXPANSION QMC ⟨𝑧|𝑉%&…𝑉%+𝑉%+,&…𝑉%(|𝑧⟩ ⟨𝑧|𝑉%&…𝑉%+𝑉𝑉𝑉%+,&…𝑉%(|𝑧⟩
ALGORITHM
Figure 3. Basic update moves of the ODE al-
We now describe the basic ingredients of our Off-
gorithm. (a) Classical moves (e.g., a single bit flip),
Diagonal Expansion (ODE) algorithm that is based whereby only the initial state z is changed to z(cid:48) leaving
on the above partition function decomposition. For
S unchanged. (b) Local swap, whereby two adjacent
q
concreteness we discuss the algorithm as it applies operators V V are interchanged changing the state
ik ik+1
tothetransverse-fieldIsingmodel,Eq.(17),however between them from z(cid:48) to z(cid:48)(cid:48). (c) Block swap, whereby
we note that generalization to other systems should two partitions of the sequence are interchanged. This
be straightforward. We first establish the compu- changes the initial state from z to z(cid:48) as well as the or-
tational complexity associated with implementing dering of the sequence. (d) Pair creation/annihilation,
this new algorithm, discussing in detail generic up- whereby a new pair of operators is inserted or deleted.
dates as well as measurements. We then present
some results that allow us to fully characterize and
to some extent quantify the advantages of the algo- i.e.,theclassicalBoltzmannweightoftheinitialran-
rithm over generic QMC methods, specifically path- dom state |z(cid:105). Here the effective classical energy
integral QMC. ECinit is the classical energy of |z(cid:105).
B. Updates
A. General description of the algorithm
We next describe the basic update moves for the
An ODE configuration is a pair C = {|z(cid:105),S }
q
algorithm. We consider here only generic local up-
where |z(cid:105) corresponds to a classical bit configura-
dates but note that updates of the global type can
tionandS =V V ···V isasequenceof(possibly
q i1 i2 iq be specifically tailored to the system in question.
repeated) off-diagonal operators. As was discussed
An update is considered local if it changes the mul-
above, each configuration C induces a list of states
tiset M by a finite (i.e., by a system-size indepen-
Z = {|z (cid:105) = |z(cid:105),|z (cid:105),...,|z (cid:105) = |z(cid:105)} (see Fig. 1), C
0 1 q
dent) number of terms, e.g., M →M +{E(z )}−
which in turn also generates a corresponding mul- C C i
{E(z )}. The basic updates are summarized in
tiset of diagonal energies M = {E ,E ,...,E } j
C 0 1 q
Fig. 3 and are discussed in detail below.
of not-necessarily-distinct values (recall that E =
i
(cid:104)z |H |z (cid:105)). Forsystemswithdiscretizedenergyval-
i c i
ues,themultisetcanbestoredefficientlyina‘multi-
plicity table’ {m ,m ,...,m ,...}, where m is the 1. Classical moves
0 1 j j
multiplicity of the energy E in the multiset. Given
j
MC, the evaluation of the effective classical energy Classical moves are any moves that involve a ma-
EC and the GBW WC follow from the definition of nipulation of the classical state |z(cid:105) while leaving Sq
the GBW, Eq. (16). The actual evaluation of the unchanged [see Fig. 3(a)]. In a single bit-flip clas-
effective classical energy is schematically given in sical move, a spin from the classical bit-string state
Fig. 2 and is discussed in more technical detail in |z(cid:105) of C is picked randomly and is flipped, gener-
Appendix B. ating a state |z(cid:48)(cid:105) and hence a new configuration C(cid:48).
The initial configuration of the ODE algorithm is Performingthischangerequiresrecalculatingtheen-
a random classical configuration |z(cid:105) and the empty ergies associated with the sequence S leading to a
q
sequenceS =1. Theweightofthisinitialconfig- new multiset M and can become computationally
q=0 C(cid:48)
uration is intensive if q is large. Classical moves should there-
fore be attempted with low probabilities if q large.
W =e−βE(z), (27) Simplyenough,theacceptanceprobabilityforaclas-
Cinit
7
sical move is 4. Creation/annihilation of off-diagonal operators
(cid:18) (cid:19)
p=min 1,WC(cid:48) =min(cid:0)1,e−β∆E(cid:1) , (28) ThemovespresentedsofarhaveleftthesizeofS
q
W
C unchanged. The creation/annihilation move shown
in Fig. 3(d) has the effect of changing the value of
where ∆E =E −E is the difference between the
C(cid:48) C q by 2, i.e., q → q ± 2, which in the transverse-
effective classical energies of the proposed configu-
fieldIsingmodelcorrespondstocreatingordestroy-
ration C(cid:48) and current configuration C. ing off-diagonal operators σx in pairs. We imple-
j
In the absence of a quantum part to the Hamil- ment this via the insertion or deletion of two adja-
tonian (Γ = 0), not only are classical moves the cent,identicaloperators. Withprobabilityp (e.g.,
del
only moves necessary, but they are also the only p = 1/2) we try to annihilate an adjacent pair,
del
moveswithanonzeroacceptanceprobability. Inthis and with probability 1−p we try to insert a pair.
del
case, theODEalgorithmautomaticallyreducestoa
For pair insertion, we randomly pick an internal
classical thermal algorithm keeping the size of the
insertion point in the sequence (we denote this in-
imaginary-time dimension at zero (q = 0) for the ternal state by |z(cid:48)(cid:105)) and a random V to insert. This
duration of the simulation. addstwonewenergiesE(z(cid:48))andE(z(cid:48)(cid:48))tothemulti-
set, where |z(cid:48)(cid:48)(cid:105)=V|z(cid:48)(cid:105). The acceptance probability
for pair creation is given by
2. Local swap
(cid:18) p Nβ2Γ2 (cid:19)
p=min 1, del e−β∆E (30)
A local swap is the swapping of neighboring off- 1−pdel(q+2)(q+3)
diagonal operators in the sequence S . A random
q
pair of adjacent off-diagonal operators in the se- where as before ∆E = EC(cid:48) − EC is the difference
quenceispickedandswapped[asshowninFig.3(b)]. between the effective classical energies of the pro-
If the state between V and V is |z(cid:105) and is |z(cid:48)(cid:105) posed configuration C(cid:48) and current configuration C
aftertheswap,thenthiekswapinikv+o1lvesaddinganen- and MC(cid:48) = MC +{E(z(cid:48)),E(z(cid:48)(cid:48))}. For deletion, we
ergy E(z(cid:48)) and removing an energy E(z) from the randomly pick an internal point in the sequence. If
energy multiset [note that E(z) and E(z(cid:48)) may be the two operators to the side of the insertion point
the same]. The acceptance probability for the move are not identical, no deletion is performed, and the
isasinEq.(28)withM =M +{E(z(cid:48))}−{E(z)}. move is rejected. If the two operators are identical,
C(cid:48) C
theyaredeletedandtherelevantenergiesE(z(cid:48))and
E(z(cid:48)(cid:48)) are removed from the multiplicity table. The
probability of acceptance for the deletion move is
3. Block-swap
(cid:18) (cid:19)
1−p q(q+1)
p=min 1, del e−β∆E , (31)
A block swap [Fig. 3(c)] is a local update that p β2Γ2N
del
involves a change of the classical state z. Here,
a random position k in the sequence S is picked
q where as before ∆E = EC(cid:48) −EC and MC(cid:48) = MC −
suchthatthesequenceissplitintotwo(non-empty) {E(z(cid:48)),E(z(cid:48)(cid:48))}.
parts, S = S S , with S = V ···V and S =
q 1 2 1 i1 ik 2 Thesizeoftheimaginarytimedimensionq comes
V ···V . The classical state |z(cid:48)(cid:105) at position k in
ik+1 iq strictlyfromoff-diagonaltermsandshrinksorgrows
the sequence is given by
depending on the strength of the ‘quantum compo-
nent’ of the model. This property is expected to
(cid:104)z(cid:48)|=(cid:104)z|S =(cid:104)z|V ...V , (29)
1 i1 ik be heavily utilized in order to overcome the freezing
of QMC algorithms in almost classical regimes. In
where|z(cid:105)istheclassicalstateofthecurrentconfigu- these regimes, q is small, and the algorithm reduces
ration. Thestate|z(cid:48)(cid:105)hasenergyE(z(cid:48)),andthestate to being a classical thermal algorithm5.
|z(cid:105) has energy E(z). We consider a new configura-
tion defined by (cid:104)z(cid:48)|S S |z(cid:48)(cid:105). The multiplicity table
2 1
of this configuration differs from that of the current
configuration by having one fewer E(z) state and
one additional E(z(cid:48)) state. The weight of the new 5 This is to be contrasted with the standard SSE formalism
where one normally introduces an additional parameter L
configuration is then proportional to e−βMC(cid:48) where in order to fix the size of imaginary time dimension for
themultisetMC(cid:48) =MC+{E(z(cid:48))}−{E(z)}. Theac- moreefficientweightcalculations. Thefixingofthesizeof
ceptanceprobabilityisasinEq.(28)withtheafore- imaginarytimemayadverselyaffecttheconvergenceofthe
mentioned M . algorithm. Here,thisparametertooisspurious.
C(cid:48)
8
C. Measurements where ∆E = EC(cid:48) −EC and C(cid:48) is the configuration
associatedwiththemultisetM =M −{E(z)}. In
C(cid:48) C
An integral part of any QMC algorithm is the ac- theaboveform,wecanreinterprettheweightWC as
quisition of various properties of the model such as contributing
averageenergy,magnetization,specificheatandcor-
qe−β∆E
relation functions. In the ODE algorithm (as in v =δ , (36)
SSE), diagonal (classical) measurements are mea- k k,iq βΓ
sured differently than off-diagonal ones.
to (cid:104)V (cid:105).
k
As in the case of the diagonal measurements, one
can take advantage of the periodicity in the imag-
1. Diagonal measurements
inary time direction to improve statistics by rotat-
ing the sequence such that any of the elements of
Diagonal operators D obey D|z(cid:105) = d(z)|z(cid:105) where
S becomes the last element of the sequence (see
d(z) is a number that depends both on the oper- q
Sec.IIIB3),weightedaccordinglybytheblock-swap
ator and the state it acts on. Since (cid:104)z|DS |z(cid:105) =
q probability. By doing so, v becomes
d(z)(cid:104)z|S |z(cid:105), for any given configuration C = k
q
(as|tgzao(cid:105)tn,isSatlqic)o,sp,tewhreeartecoarinstahalescroomncatorlnibasuvideteirorangreodt(cid:104)a=Dti(cid:105)od.n(zsT)ionto(imtthhpeeropdveie-- vk =(cid:88)i βqΓ(cid:80)jqe=−−01βeE−CβiECj ee−−ββEECCi(cid:48) = ZqNβkΓe−βEC(cid:48)
(37)
riodic)imaginarytime. Todothat,wemayconsider
(cid:80)
where M = M +{E(z )}−{E(z)}, the sum
‘virtual’ block-swap moves (see Sec. IIIB3) that ro- Ci C i i
is over all rotated configurations C(cid:48) whose S ends
tateS andasaresultalsochangetheclassicalcon- q
q
with V , and N is the number of times V appears
figuration from |z(cid:105) to |z (cid:105). The contribution to the k k k
i
in the sequence S .
expectation value of a diagonal operator D thus be- q
comes:
1 (cid:88)q−1 3. Products of off-diagonal measurements
d= Z d(zi)e−βECi . (32)
i=0 Thesamplingoftheexpectationvaluesoftheform
whereE istheeffectiveclassicalenergyassociated (cid:104)V V (cid:105)proceedsverysimilarlytothesingleopera-
Ci k1 k2
withconfigurationC whosemultisetisM =M + torcaseexceptthatnowbothoperatorsmustappear
i Ci C
{E(z )}−{E(z)}(recallthatz ≡z,soM =M ). at the end of the sequence. The argument proceeds
i 0 C0 C
The normalization factor Z above is the sum similarlytothesingleoff-diagonalmeasurement,and
we have that the contribution to the expectation
q−1
Z =(cid:88)e−βECj =(cid:88)mje−βECj (33) value of (cid:104)Vk1Vk2(cid:105) is
j=0 j q(q−1)e−βEC(cid:48)
v =δ δ , (38)
overallnonzeromultiplicitiesmj. Inthecasewhere k1,k2 iq−1,k1 iq,k2 β2Γ2 e−βEC
D =H the above expression simplifies to:
c
with M = M −{E(z),E(z )}. As in the sin-
C(cid:48) C q−1
q−1
1 (cid:88) 1 (cid:88) gle off-diagonal operator case, we can use the block-
d= Z E(zi)e−βECi = Z mjE(zj)e−βECj . swap move to alter the elements at the end of the
i=0 j sequence, and for each pair of adjacent operators in
(34)
the sequence obtain an improved contribution
q(q−1)(cid:88) e−βECi e−βECi(cid:48)
2. Off-diagonal measurements vk1,k2 = β2Γ2 (cid:80)q−1e−βECj e−βECi
i j=0
Wenextconsiderthecaseofmeasuringtheexpec- = q(q−1)(cid:88)e−βECi(cid:48) , (39)
tation value of an off-diagonal operator V , namely, Zβ2Γ2
k
i
(cid:104)V (cid:105). Todothis,weinterprettheinstantaneouscon-
k
figuration as follows where MCk = MC + {E(zk)} − {E(z)}, MCi(cid:48) =
M −{E(z),E(z(cid:48)(cid:48))} with |z(cid:48)(cid:48)(cid:105) = V |z(cid:48)(cid:105) and |z(cid:48)(cid:105) is
(βΓ)qe−βEC C k2
W = (cid:104)z|S |z(cid:105) (35) the classical state after the block swap. Similar to
C q! q (cid:80)
the single off-diagonal operator case, the sum is
=(cid:18) βΓ (cid:19)(cid:20)(βΓ)q−1e−βEC(cid:48)(cid:104)z|S V |z(cid:105)(cid:21) , overallrotatedconfigurationsC(cid:48)whoseSq endswiith
qe−β∆E (q−1)! q−1 iq Vk1Vk2.
9
Measurements of thermal averages of products of to exact diagonalization is feasible. An example is
more than two off-diagonal operators can also be given in Fig. 5 illustrating the excellent agreement
derived in a straightforward manner. ofODEwiththeexact-numericalvalues,eveninthe
high-β but low-Γ regime where PIQMC begins to
show deviations from the exact results. Increasing
IV. RESULTS the number of measurements for PIQMC rectifies
this discrepancy, but the deviation already suggests
Having laid the groundwork for the ODE QMC thattheODEalgorithmmayrequirefewermeasure-
algorithm, we present in this section some results mentsoverPIQMCinthelow-Γbutlarge-β regime.
thathighlightsomeofitspropertiesandadvantages
overexistingQMCtechniques,specificallyacluster-
updatesPIQMCalgorithm6. Forbenchmarkingpur-
ODE
poses, we study random 3-regular MAX2SAT in-
PIQMC
stancesaugmentedwithatransversefield. Thisclass
Exact diag
ofinstancescorrespondstoaparticularchoiceofthe
Ising Hamiltonian given in Eq. (17), whereby each
spin is coupled antiferromagnetically (with strength
J =1)withexactlythreeotherspinspickedatran-
ij
dom (see Fig. 4 for an illustration). We study this
classofinstancesasitisknowntoexhibitaquantum
spin glass phase transition and is notoriously diffi-
cult to simulate by standard QMC techniques (see,
e.g., Refs. [9, 24]).
(a)
ODE
PIQMC
Exact diag
Figure 4. Connectivity of a randomly generated
N = 36-spin 3-regular MAX2SAT instance. Here
thediamondsdenotespinsandtheedgesdenoteantifer-
romagnetic couplings with strength J = 1. Each spin
ij
is connected to three other randomly chosen spins.
(b)
Figure 5. Agreement between ODE, PIQMC
A. Correctness of algorithm and elastic and exact diagonalization for small systems.
imaginary time The thermal expectation value of the internal en-
ergy per spin (cid:104)H(cid:105)/N and specific heat per spin C =
β2(cid:0)(cid:104)H2(cid:105)−(cid:104)H(cid:105)2(cid:1)/N fora3-regularMAX2SATinstance
As a preliminary test, we verify that we are able
ofsizeN =12forarangeofβwithΓ=(10β)−1/2ascal-
to reproduce the correct thermal expectation val-
culated using ODE, PIQMC (with 5120 Trotter slices),
ues for sufficiently small systems where comparison
and exact diagonalization. Error bars correspond to 2σ
generated by performing 1000 bootstraps over the mea-
surements.
6 WeuseWolffclusterupdates[23]alongtheimaginarytime
direction. WenextstudyinFig.6thedependenceoftheav-
10
erage size of the imaginary time dimension, namely, show next, this allows us to naturally unify the
q on system size N, inverse-temperature β, and classical Parallel Tempering (CPT) algorithm (also
transverse field strength Γ7. As was discussed ear- known as ‘exchange Monte Carlo’) [25, 26] and its
lier,theODEQMCdoesnotpresumea-prioriasize quantumcounterpart(QPT,seee.g.,Ref[27]). CPT
for the imaginary time dimension but rather allows is a refinement of the simulated annealing algo-
it to be set dynamically during the simulation. As rithm [28], whereby N replicas of an N-spin sys-
T
is shown in Fig. 6(a), as the simulation advances, tem at inverse-temperatures β <β <...<β
1 2 NT
the instantaneous q which starts at q =0 gradually undergo Metropolis spin-flip updates independently
grows and eventually fluctuates around an average of one another and in addition, replicas with neigh-
value indicating the size of the imaginary time di- boringtemperaturesregularlyattempttoswaptheir
mension. Asweexpect,theaveragevalueofq,which temperatureswithprobabilitiesthatsatisfydetailed
we denote (cid:104)q(cid:105), scales linearly with N and β with balance [29]. In this way, each replica performs a
√ √
fluctuations on the order of N and β [Figs. 6(b) random-walk on the temperature axis, which gen-
and (c), respectively]. Moreover, we find that (cid:104)q(cid:105) erally allows for quicker equilibration of the system
does indeed grow with the quantum strength of the in comparison to other techniques. Analogously in
model. Specifically, we find it to scale quadratically QPT, temperature is replaced by a parameter Γ of
with Γ as indicated in Fig. 6(d). the (quantum) Hamiltonian, e.g., the strength of
the transverse magnetic field in the transverse Ising
model, and each replica performs a random-walk on
B. ODE vs PIQMC the Γ axis.
Both CPT and QPT are two widely used vari-
Since the value of q determines the cost of cal- ations on Monte Carlo schemes but have so far
culating the GBWs, our results in Fig. 6 indicate been considered as separate algorithms. The cur-
thattheODEalgorithmcanhavesignificantadvan- rent formulation allows to unify the two tempering
tages in the low-Γ but large-β regime. For the 3- algorithms in a straightforward manner. A natural
regular MAX2SAT class, this would be in the spin- generalization is to consider a tempering algorithm
glass phase, where we can expect QMC algorithms that traces an arbitrary curve in the this classical-
to become less efficient. We quantify this possible quantum β-Γ plane. This opens up the opportu-
advantage by comparing the performance of our al- nity to study, e.g., certain properties of experimen-
gorithmagainstPIQMCinthisregime. InFig.7,we tal quantum annealers (see for example Ref. [30])
compare the warm-up time required to reach close which trace such quantum-classical curves as well
to the thermal state for the two algorithms. We as to study classical-quantum optimization tech-
observethatinorderforthe(discrete-time)PIQMC niques and equilibration methods, by, e.g., looking
algorithmtoachievethis,weneedasufficientlylarge forcurvesthatwouldallowonetobypassfirstorder
Trotterslicing(>1024), whichinturn increasesthe phase transitions.
time cost of performing a sweep in the simulations. If we consider replicas along a curve in the β-Γ
In this regard, the ODE algorithm reaches the ther- plane at points {(β1,Γ1),...,(βNT,ΓNT)}, then a
mal state in less computational time, with even a paralleltemperingswapprobabilitybetweenthei-th
factor of 10 advantage when compared to PIQMC and (i+1)-th replica is given by:
with 2048 Trotter slices.
(cid:18) W (β ,Γ )W (β ,Γ )(cid:19)
P =min 1, Ci i+1 i+1 Ci+1 i i (,40)
W (β ,Γ )W (β ,Γ )
Ci i i Ci+1 i+1 i+1
C. Quantum-classical parallel tempering
where the above weight ratio is conveniently simpli-
As we demonstrated in Sec. II, the ODE parti- fied to:
tionfunctiondecompositionnaturallyreducestothe
W (β ,Γ )W (β ,Γ )
classical one when the strength of the off-diagonal Ci i+1 i+1 Ci+1 i i = (41)
W (β ,Γ )W (β ,Γ )
terms in the Hamiltonian are sent to zero. As we Ci i i Ci+1 i+1 i+1
(cid:18) βiΓi (cid:19)qi+1−qi e−βi(EC(cid:48)i+1−ECi) ,
βi+1Γi+1 e−βi+1(ECi+1−EC(cid:48)i)
7 Thewarm-upofthesimulationsinvolvedalinearannealin where E and E are the effective classical ener-
βfromaninitialvaluethatisafactor103smallerthanthe Ci Ci+1
gies of configurations C and C , respectively and
targetβtothetargetβ. 106 sweepsareperformedintotal i i+1
duringthewarm-up. Afterthewarm-up,104measurements EC(cid:48)i and EC(cid:48)i+1 are the effective classical energies of
are performed, with 102 sweeps between measurements to these configurations when calculated with switched
ensurethesubsequentmeasurementsareuncorrelated. β and Γ.