Simulating Quantum Annealing

with Quantum Monte Carlo

Guglielmo Mazzola (ETH)
!
M. Troyer (ETH)
S. Isakov, V. Smelyanskiy, S. Boixo, H. Neven (Google)
Z. Jiang (NASA)
1
| |
1. Introduction: QA vs SQA (the story so far)!
2. Introduction: Quantum Monte Carlo (QMC) for SQA
3. Tunneling with QMC and with QA
4. Conclusions

Outline

Guglielmo Mazzola | QuAASI 2016 | 2

 SQA
~2000

Simulations

QA SA
2011
1983

Hardwares

CA
~10.000 b.C. Guglielmo Mazzola | QuAASI 2016 | 3
Implement “Quantum Mechanics”

SQA

Implement classical
optimization

QA SA

CA Guglielmo Mazzola | QuAASI 2016 | 4

 SQA

Run on Classical Hardwares

QA SA

Quantum Hardware

Guglielmo Mazzola | QuAASI 2016 | 5

 SQA

SA
Santoro et. al. Science 295, 2427 (2001)

2001: early evidence that SQA

outperform SA for 2D random Ising
glasses 6
Guglielmo Mazzola | QuAASI 2016 |
2014: D-Wave
correlates with SQA SQA

QA

Boixo et. al. Nat. Phys. 10, 218 (2014)

Guglielmo Mazzola | QuAASI 2016 | 7
2014: D-Wave

cl ee
as d
does not display

sp

si up
du m

ca
p
sp antu
quantum speed-

l
ee
up compared to

qu
SA.

Ronnow et. al. Science, 345, 420 (2014)

QA SA

Guglielmo Mazzola | QuAASI 2016 | 8

 SQA
2014: D-Wave
correlates with SQA 2001: SQA with QMC
outperform SA for 2D
random Ising glasses

QA SA
2014: D-Wave
does not display
N.B. stoquastic drivers such
as transverse field quantum speed-
hamiltonians up compared to
always assumed
SA.
Guglielmo Mazzola | QuAASI 2016 | 9

V. Denchev et. al. arXiv:1510.08057

 SQA
2014: D-Wave
correlates with SQA 2001: SQA with QMC
outperform SA for 2D
random Ising glasses

QA SA
2014: D-Wave
does not display
N.B. stoquastic drivers such
as transverse field quantum speed-
hamiltonians up compared to
always assumed
SA.
Guglielmo Mazzola | QuAASI 2016 | 10

V. Denchev et. al. arXiv:1510.08057

1. Introduction: QA vs SQA (the story so far)
2. Introduction: Quantum Monte Carlo (QMC) for SQA!
3. Tunneling with QMC and with QA
4. Conclusions

Outline

Guglielmo Mazzola | QuAASI 2016 | 11

 Simulating Quantum Annealing

Ideally we would like to solve time- d

dependent Schroedinger equation:
| (t)i = iH(t)| (t)i
dt
Hilbert space scales exponentially with the size! System size: ~40 spins

| |
Quantum tunneling in the smallest droplet of water

Richardson et al. "Concerted hydrogen-

bond breaking by quantum tunneling in
the water hexamer prism."

Science 351, 1310 (2016)

Guglielmo Mazzola | QuAASI 2016 | 13

 Quantum Monte Carlo simulations

Ultracold bose gas

S. Trotzky, L. Pollet, F. Gerbier, U.
Schnorrberger, I. Bloch, N. V. Prokof’ev, B.
Svistunov & M. Troyer

Nature Physics 6, 998–1004 (2010)
Guglielmo Mazzola | QuAASI 2016 | 14
 Simulating Quantum Annealing

Ideally we would like to solve time- d

dependent Schroedinger equation:
| (t)i = iH(t)| (t)i
dt
Hilbert space scales exponentially with the size! System size: ~40 spins

Quantum Monte Carlo (PIMC)

Computational effort scales linearly with size!

System size: ~10.000 spins

Guglielmo Mazzola | | 15
Santoro et. al. Science 295, 2427 (2001) AQC 2016
 Path Integral Monte Carlo

H H
Density matrix is an operator: ⇢=e Z = Tr e
Matrix elements are difficult to H 0 H = H1 + H 2
hq|e |q i because [H1 , H2 ] 6= 0
compute:
h iM
⌧H ⌧ H1 ⌧ H2 H H1 H2
instead e ⇡e e e = lim e M e M
M !1

Z
hq|e H
|q 0 i = dq1 dq2 · · · dqM hq|e ⌧H
|q1 i q

hq1 |e ⌧H
|q2 i · · · hqM |e ⌧H 0
|q i
··· q0

0 ⌧ 2⌧ 3⌧

Guglielmo Mazzola | QuAASI 2016 | 16

 Path Integral Monte Carlo

H H
Density matrix is an operator: ⇢=e Z = Tr e
Matrix elements are difficult to H 0 H = H1 + H 2
hq|e |q i because [H1 , H2 ] 6= 0
compute:
h iM
⌧H ⌧ H1 ⌧ H2 H H1 H2
instead e ⇡e e e = lim e M e M
M !1

Z
hq|e H
|q 0 i = dq1 dq2 · · · dqM hq|e ⌧H
|q1 i q
⌧H ⌧H 0 q0
hq1 |e |q2 i · · · hqM |e |q i
q↵

Guglielmo Mazzola | QuAASI 2016 | 17

 Path Integral Monte Carlo

H H
Density matrix is an operator: ⇢=e Z = Tr e
Matrix elements are difficult to H 0 H = H1 + H 2
hq|e |q i because [H1 , H2 ] 6= 0
compute:
h iM
⌧H ⌧ H1 ⌧ H2 H H1 H2
instead e ⇡e e e = lim e M e M
M !1

Z
H ⌧H
hq|e |qi = dq1 dq2 · · · dqM hq|e |q1 i q
⌧H ⌧H
hq1 |e |q2 i · · · hqM |e |qi q↵

Guglielmo Mazzola | QuAASI 2016 | 18

 Path Integral Monte Carlo

Z
H ⌧H
hq|e |qi = dq1 dq2 · · · dqM hq|e |q1 i q
⌧H ⌧H
hq1 |e |q2 i · · · hqM |e |qi

Guglielmo Mazzola | QuAASI 2016 | 19

 Path Integral Monte Carlo
Z Z
H ⌧H
Z= dq hq|e |qi = dq dq1 dq2 · · · dqM hq|e |q1 i
⌧H ⌧H
hq1 |e |q2 i · · · hqM |e |qi

X
S[path]
Z= e
NB: only possible for stoquastic H
paths

Sum of all possible paths or

trajectories in imaginary time q(⌧ )
q
S[q(⌧ )]
Each path contributes with e
@S[q(⌧ )]
Dominant contributions come from paths =0
@q(⌧ )
Form of S[q(⌧ )] is system dependent. Guglielmo Mazzola | QuAASI 2016 | 20
 Path Integral Monte Carlo: spin systems

X X
z z
H= Jij i j + sxi
ij i

?
J

S[q(⌧ )] H/M Jij

e =e 0 1
M
X X X
H= @ k k
Jij si sj + J ? k k+1 A
Jij si si
k=1 ij i

Guglielmo Mazzola | QuAASI 2016 | 21

 Path Integral Monte Carlo: physical limit
h iM PIMC samples the correct partition
H M H1 M H2
e = lim e e function in the physical limit M ! 1
M !1

M =4
M =8 M = 16 M =1
Continuous time limit (CT)
One could also optimise M at given

Guglielmo Mazzola | QuAASI 2016 |

Heim et. al. Science 348, 215 (2015) 22
1. Introduction: QA vs SQA (the story so far)
2. Introduction: Quantum Monte Carlo (QMC) for SQA
3. Tunneling with QMC and with QA!
4. Conclusions

Outline

Guglielmo Mazzola | QuAASI 2016 | 23

Experiment: compare runtimes of QA device vs QMC

The runtime of the (ideal) Quantum How does the runtime of a QMC
device is dictated by the simulated annealing algorithm scale?
Quantum Adiabatic theorem:
2 TQM C / TQA ?
TQA /
Guglielmo Mazzola | QuAASI 2016 | 24
 Path Integral Monte Carlo pseudodynamics
Doing the integral with Monte Carlo by sampling ring-
S[q(⌧ )]
polymer configurations (paths) with Metropolis weight e
y
q
Evolution of the classical
path as a function of the
simulation time t

q(⌧, t)
given by the Metropolis
pseudo-dynamics
(updates).

@q(⌧, t) S[q(⌧, t)]

= + ⌘(⌧, t)
@t q(⌧, t)

qx Guglielmo Mazzola | QuAASI 2016 | 25

 Path Integral Monte Carlo pseudodynamics
Doing the integral with Monte Carlo by sampling ring-
S[q(⌧ )]
polymer configurations (paths) with Metropolis weight e
y
q Evolution of the classical
path as a function of the
0⌧ < simulation time t

0  t < TQM C q(⌧, t)

given by the Metropolis
pseudo-dynamics

⌧ (updates).

d
| (t)i = iH(t)| (t)i
dt

qx Guglielmo Mazzola | QuAASI 2016 | 26

 Transition states of the dynamics: instantons
Consider double well problems @q(⌧, t) S[q(⌧, t)]
= + ⌘(⌧, t)
@t q(⌧, t)
S[q(⌧, t)] Transition state (TS)
Transition states defined by: =0 is the saddle point
q(⌧, t) with non-trivial
boundary conditions
x x
⇤⇤
q(⌧, t = TQM C ) for double well is known: the instanton x (⌧ )

S[x⇤⇤ (⌧ )] 2
e ⇠
q ⇤⇤ (⌧ )

q(⌧, t = 0) 0
y Guglielmo Mazzola | AQC 2016 |
⌧
27
 Transition states of the dynamics: instantons
Consider double well problems @q(⌧, t) S[q(⌧, t)]
= + ⌘(⌧, t)
@t q(⌧, t)
S[q(⌧, t)] Transition state (TS)
Transition states defined by: =0 is the saddle point
q(⌧, t) with non-trivial
boundary conditions
x
⇤⇤
for double well is known: the instanton x (⌧ )

S[x⇤⇤ (⌧ )] 2
e ⇠

0
Guglielmo Mazzola | AQC 2016 |
⌧
28
 Transition states of the dynamics: instantons
The classical field, in a PIMC simulation, evolves through a Langevin equation,

@q(⌧, t) S[q(⌧, t)]

= + ⌘(⌧, t)
@t q(⌧, t)
x
In a double well model, we know the transition state
⇤⇤
(transition path or trajectory in imaginary time) q (⌧ )
⌧

The escape rate of this classical thermally activated event is given by

Kramers theory (Boltzmann weight at the TS)

S[q ⇤⇤ (⌧ )] 2
k/e ⇠

Therefore we expect that the QMC tunneling rate must scale as ⇠ 2

Guglielmo Mazzola | QuAASI 2016 | 29

 Continous space model: 1D double well
@2
H= 2
+ V (x)
@x Z
S= d⌧ ẋ2 (⌧ ) + V (x(⌧ ))
0

< Ttunn >QM C

2
↵
q ⇤⇤ (⌧, t = Ttunn )
1/

q(⌧, t = 0)
Ttunn : Number of QMC updates required to generate q ⇤⇤ (⌧ ) x0
Guglielmo Mazzola | AQC 2016 | 30
 Ferromagnetic Ising system
X X
z z x
Consider now a spin system. H=J i j + i
i,j i

L R

z
1 m =< > 1

Path integral construction lead to

1
0 =p ( L + R)
an extended lattice (formally 2
similar to the previous ring
polymer) Guglielmo Mazzola | QuAASI 2016 | 31
 QMC tunneling rate in ferromagnetic Ising system

Let’s measure QMC tunneling time as a function

of the system size L.

2
< TP BC >⇠ 1/

“Understanding Quantum Tunneling through Quantum

Monte Carlo Simulations”
S.V. Isakov, G. Mazzola, V.N. Smelyanskiy, Z. Jiang,
S. Boixo, H. Neven, and M. Troyer
arXiv:1510.08057 Guglielmo Mazzola | AQC 2016 | 32
 Transition states of the dynamics: instantons
Computing the thermal density matrix requires closed paths in imag. time
Z
H S[x⇤⇤ (⌧ )]
Z= dq hq|e |qi e ⇠ 2

x x
S[x⇤ (⌧ )]
e ⇠

⇤
x (⌧ )

⇤⇤
x (⌧ )
y ⌧
0 Guglielmo Mazzola | QuAASI 2016 | 33
 QMC tunneling rate in ferromagnetic Ising system

Let’s measure QMC tunneling time as a function

of the system size L.

< TOBC >⇠ 1/

bL
p Le
1 ⇡a
2
< TP BC >⇠ 1/

“Understanding Quantum Tunneling through Quantum

Monte Carlo Simulations”
S.V. Isakov, G. Mazzola, V.N. Smelyanskiy, Z. Jiang,
S. Boixo, H. Neven, and M. Troyer
arXiv:1510.08057 Guglielmo Mazzola | AQC 2016 | 34
 Conclusions /1

We show that, for a double well model,

we can identify the transition state of the PIMC
⇤⇤
pseudodynamics, the instanton q (⌧ ),
2
and the autocorrelation time scales as

We propose a simple modification of the PIMC

algorithm, with open boundary conditions in
imaginary time, which implies a quadratic
speedup (quantum inspired classical algorithm).
V. Denchev et. al. arXiv:1510.08057

We are currently working on finding possible exceptions to this matching, exploring models
(borrowed from chemistry) which feature multidimensional tunnelling.

Clearly, also Hamiltonians with sign-problem (non stoquastic), beyond the Transverse
field, may show scaling advantage compared to QMC.
Guglielmo Mazzola | QuAASI 2016 | 35
Thank you!

36
| |