Quantum Annealing with Markov Chain Monte Carlo Simulations

37 Pages·2016·1.73 MB·English

StatisticalScience
2016,Vol.31,No.3,362–398
DOI:10.1214/16-STS560
©InstituteofMathematicalStatistics,2016
Quantum Annealing with Markov Chain
Monte Carlo Simulations and D-Wave
Quantum Computers
Yazhen Wang, Shang Wu and Jian Zou
Abstract. Quantum computation performs calculations by using quantum
devicesinsteadofelectronicdevicesfollowingclassicalphysicsandusedby
classical computers. Although general purpose quantum computers of prac-
tical scale may be many years away, special purpose quantum computers
are being built with capabilities exceeding classical computers. One promi-
nent case is theso-called D-Wavequantum computer, which is acomputing
hardware device built to implement quantum annealing for solving combi-
natorial optimization problems. Whether D-Wave computing hardware de-
vices display a quantum behavior or can be described by a classical model
hasattractedtremendousattention,anditremainscontroversialtodetermine
whetherquantumorclassicaleffectsplayacrucialroleinexhibitingthecom-
putationalinput–outputbehaviorsoftheD-Wavedevices.Thispaperconsists
of two parts where the first part provides a review of quantum annealing
anditsimplementations,andthesecondpartproposesstatisticalmethodolo-
gies to analyze data generated from annealing experiments. Specifically, we
introduce quantum annealing to solve optimization problems and describe
D-Wave computing devices to implement quantum annealing. We illustrate
implementations of quantum annealing using Markov chain Monte Carlo
(MCMC) simulations carried out by classical computers. Computing exper-
iments have been conducted to generate data and compare quantum anneal-
ingwithclassicalannealing.Weproposestatisticalmethodologiestoanalyze
computingexperimentaldatafromaD-Wavedeviceandsimulateddatafrom
the MCMC based annealing methods, and establish asymptotic theory and
checkfinitesampleperformancesfortheproposedstatisticalmethodologies.
Our findings confirm bimodal histogram patterns displayed in input–output
datafromtheD-WavedeviceandbothU-shapeandunimodalhistogrampat-
ternsexhibitedininput–outputdatafromtheMCMCbasedannealingmeth-
ods. Further statistical explorations reveal possible sources for the U-shape
patterns. On the other hand, our statistical analysis produces statistical evi-
dencetoindicatethatinput–outputdatafromtheD-Wavedevicearenotcon-
sistent with the stochastic behaviors of any MCMC based annealing models
under the study. We present a list of statistical research topics for the future
studyonquantumannealingandMCMCsimulations.
Keywordsandphrases: Quantumannealing,quantumcomputing,Markov
chain Monte Carlo, Ising model, ground state success probability, Hamilto-
nian,quantumbit(qubit).
YazhenWangisChairandProfessor,Departmentof GraduateStudent,DepartmentofStatistics,Universityof
Statistics,UniversityofWisconsin-Madison,Madison,WI Wisconsin-Madison,Madison,WI53706,USA(e-mail:
53706,USA(e-mail:[email protected]).ShangWuis [email protected]).JianZouisAssistantProfessor,
362
QUANTUMANNEALINGANDMCMCSIMULATIONS 363
1. INTRODUCTION certifytheirbehaviorsaccordingtoclaimsand/orspec-
ifications. One prominent example is the commercial-
Quantum computation is based on the idea of us-
ization of a special purpose quantum computer man-
ing quantum devices to process information and per-
ufactured by D-Wave Systems Incorporation for solv-
form computation, instead of electronic devices fol-
ingcombinatorialoptimizationproblems(Boixoetal.,
lowing the laws of classical physics and used by clas-
2014a, 2014b, Boixo et al., 2015a, 2015b, Johnson
sicalcomputers(NielsenandChuang,2000andWang,
et al., 2011, Jones, 2013, McGeoch, 2014, Rønnow
2012). Here, classical computers mean today’s elec-
et al., 2014 and Vinci et al., 2014, Hen et al., 2015,
tronic based computers. Two quantum computing ap-
Venturelli et al., 2014, and Martin-Mayor and Hen,
proachesarelogic-gatebasedquantumcomputingand
2015).
adiabatic quantum computing (Aharonov et al., 2007,
The D-Wave quantum computer is an analog com-
Deutsch, 1985, DiVincenzo, 1995, Browne, 2014,
puting hardware device that is designed and built to
Farhietal.,2000,2001,2002andJohnsonetal.,2011).
physically implement quantum annealing. It consists
Logic-gate based quantum computing has as its pur-
ofaquantumprocessorbasedonsuperconductingflux
pose the development of a quantum version of classic
quantum bits (where quantum bits, called qubits for
logic gate operations and the construction of general
short, are the quantum analog of classic bits ones and
purpose (or universal) quantum computers. Adiabatic
zeros, see more details about qubits in Section 2) and
quantum computing is based on quantum annealing to
surrounding system such as a cooling apparatus and
build special purpose quantum computers (i.e., quan-
a magnetic shielded box (Boixo et al., 2014a, 2014b,
tum annealers) for solving tough combinatorial opti-
ClarkeandWilhelm,2008,Devoret,WallraandMarti-
mization problems, where quantum annealing is the
nis, 2004, DiCarlo et al., 2009, Johnson et al., 2011,
quantum analog of classical annealing such as simu-
Lanting et al., 2014 and Pudenz, Albash and Lidar,
lated annealing with thermodynamics in classical an-
2014).Thehardwaredeviceisaquantumannealer,and
nealing replaced by quantum dynamics, and quantum
thequbitsareanarrayofchilledsuperconducting nio-
annealers refer to physical hardware implementations
biumloopsthatcanbeengineeredtoveryquicklyfind
of quantum annealing (see more details about classi-
the lowest point in an energy “landscape” of hills and
cal annealing and quantum annealing later in this sec-
valleysassociatedwithaquantumsystem.Theanneal-
tionandSection3).SinceitsintroductionbyFeynman ing idea is to define an optimization problem through
(1981/82), quantum computation has been proclaimed the energy landscape and represent the lowest energy
to be the future of computing, and intensive research point as the solution of the problem posed. Such op-
effortsareunderwayaroundtheglobetoinvestigatea timization problems turn up in everything from travel
numberoftechnologiesthatcouldleadtomorepower- salesman problem to integer factoring, from software
fulandmoreprevalentquantumdevicesforbettercom- verification and validation to object recognition and
putation,communicationandcryptography.Ontheone classification,andfromgenomesequenceanalysisand
hand, in spite of tremendous progresses made in the proteinfoldingtoportfoliooptimizationandriskanal-
past two decades, general purpose quantum comput- ysis(Brookeetal.,1999,Farhietal.,2000,2001,2002,
ers of practical scale still have a long way to go. On McGeoch, 2014, and Shor, 1994). For example, D-
theotherhand,thedevelopmentofquantum technolo- Wave devices have been tested on simple application
giesisatthecriticalpointwherequantumcommunica- problems in graphs and networks, machine learning,
tion devices and special purpose quantum computers, artificial intelligence and computational biology (Bian
such as quantum annealers, quantum simulators and et al., 2012, O’Gorman et al., 2014, Perdomo-Ortiz
quantum crypto devices, can be built with capabilities et al., 2012, Perdomo-Ortiz et al., 2014, andRieffel
exceeding classical computer based devices (Aspuru- et al., 2014). As a case in point, the lowest-energy ar-
Guzik et al., 2005, Britton et al., 2012, Browne, 2014, rangement of a protein is thought to be its preferred
Brumfiel, 2012, Nielsen and Chuang, 2000, Neumann state, and protein folding is to find the lowest-energy
et al., 2008, and Wang, 2012). For these quantum de- point in its energy landscape; a D-Wave device has
vices, it becomes increasingly important to test and/or been arranged to manipulate the qubits to reach their
lowest-energy state and solve the problem of folding
DepartmentofMathematicalSciences,Worcester a simple protein (McGeoch, 2014 and Perdomo-Ortiz
PolytechnicInstitute,Worcester,MA01609,USA(e-mail: et al., 2012). For a physical system, its lowest-energy
[email protected]). point is referred to as a ground state in physics, with
364 Y.WANG,S.WUANDJ.ZOU
itshigherenergypointsasexcitedstates.Analogousto further tested on a 512-qubit D-Wave device to dis-
classicalannealingsuchassimulatedannealing,quan- tinguish (or associate) its input–output behavior from
tum annealing finds a ground state by allowing the (or with) those obtained from the classical and quan-
qubitstoexploitaquantumeffectcalledquantumtun- tum annealing models in terms of the population ra-
nelinginthesensethattheygothroughtheenergyhills, tiooftheisolatedandclusteredgroundstates.Rønnow
rather than climbing over the energy hills in classi- et al. (2014) reported to find no evidence of quantum
calthermalannealing(seemoredetailsaboutquantum speedup for a 512-qubit D-Wave device on the prob-
tunnelinginSection3).Itmaybemoreefficienttofind lem of random spin glass instances, despite of some
a ground state using quantum annealing via quantum otherspeedupclaims(Browne,2014,Henetal.,2015,
tunneling than classical annealing via thermal jump Katzgraber,HamzeandAndrist,2014,McGeoch,2014
(Boixoetal.,2015a,2015b, Brookeetal.,1999,Farhi andVenturellietal.,2014).
et al., 2001, Jones, 2013, McGeoch, 2014, Perdomo- AttemptstoquantifythequantumnatureofD-Wave
devices have been not only met with excitement but
Ortizetal.,2014,Pudenz,AlbashandLidar,2014,and
alsoconfrontedwithsuspicion.Thesestudiesposefun-
Santoroetal.,2002).
damentalquestionsregardingthedistinguishabilitybe-
In contrast to the fact that research laboratories can
tweenquantumannealersandclassicalthermalanneal-
usually manage quantum computers with up to about
ers.Wemayboilthedistinguishabilityquestionsdown
a dozen of qubits, D-Wave devices utilize a solid state
to statistical analysis of annealing methods and the
architecture with over a thousand of interlaced super-
associated statistical inference problems. The consis-
conducting flux qubits. The manufacturing methods
tenceorinconsistenceclaimsintheliteratureregarding
andcomputingtechnologiesoftheD-Wavedevicesare
the studies of the D-Wave devices along with MCMC
well documented, yet it is challenging to understand
based annealing methods are largely based on causal
their computational power. Whether the D-Wave de-
informal inspection of correlations and histogram pat-
vicesdisplayalarge-scalequantumbehaviororcanbe
terns. The first part of the paper reviews quantum an-
described by a classical model has attracted tremen-
nealing and its implementations by D-Wave devices
dous attention, and it remains controversial to deter-
and MCMC based annealing approaches. We illus-
mine whether quantum or classical effects play a cru-
trate computing experiments for quantum annealing
cial role in exhibiting the computational input–output
andcarryouttheMCMCsimulationsbyclassicalcom-
behaviors of the D-Wave devices. Boixo et al. (2014a,
puters to generate data and compare quantum anneal-
2014b)employedrandomIsingmodelinstancestotest
ing with classical annealing. The second part of the
on a 128-qubit D-Wave device and compared the out-
paper proposes statistical methodologies to analyze
puts to the data generated from MCMC based anneal-
theinput–outputbehaviorsoftheD-Wavedevicesand
ingmodels.Theyfoundthatthecomputingexperimen-
MCMC based annealing models. We have established
tal data exhibit similar bimodal histogram patterns for
theasymptotictheorytojustifytheproposedstatistical
the D-Wave device and simulate quantum annealing
methodologies, and conducted numerical simulations
butunimodalhistogrampatternsforsimulatedanneal-
to check their finite sample performances. The statis-
ing and a classical spin model. The results are used to
tical shape tests support bimodal but not U-shape his-
argue that the D-Wave performance is consistent with
togram patterns exhibited in the data from a D-Wave
simulatedquantumannealingbutinconsistentwithbe-
device, and confirm U-shape histogram patterns in the
haviorsoftheclassicalsimulatedannealingandclassi- datafromthesimulatedquantumannealingandSSSV
calspinmodels.Shinetal.(2014)proposedaMCMC modelsandunimodal(includingmonotone)histogram
model,whichissimplyreferredtoastheSSSVmodel, patterns in the data from the classical simulated an-
to yield at least as strong correlation with the input– nealing model. We further explore the energy gap be-
output data of the D-Wave device as simulated quan- tweenthegroundstatesandtheexcitedstates,andun-
tum annealing and disputed the implication in Boixo cover some possible sources for the U-shape patterns.
etal.(2014a,2014b)ofquantumeffectsintheD-Wave Ouranalysisresultsprovidestatisticalevidencetosug-
device. Albash et al. (2014) adopted a measure based gestthattheinput–outputdatafromtheD-Wavedevice
on the total variation distance between the probabil- arenotconsistentwiththestatisticalbehaviorsofthese
ity distributions over excited states and demonstrated MCMCbasedannealingmodels.
thattheD-WavedataarenotwellcorrelatedwithSQA The rest of the paper proceeds as follows. Section 2
or SSSV data by the measure. Vinci et al. (2014) provides a brief introduction to quantum mechanics
QUANTUMANNEALINGANDMCMCSIMULATIONS 365
and quantum computation. The main review and re- described by its quantum state, often given in terms
searchpartsofthepaperaregivenbySections3and4, of a unit vector |ψ(cid:2) in Cd. To study the quantum sys-
respectively. Section 3 introduces quantum annealing tem,weperformmeasurementsonthesystemtoobtain
anditsimplementationsbyD-Wavedevicesandbythe data,whichmaybecarriedoutintermsoftheso-called
MCMCbasedannealing methodsinthecontextofthe observables.AnobservableMisdefinedtobeaHermi-
Ising model. We describe D-Wave computing experi- tianmatrixonCd.Assumethattheeigen-ecomposition
ments and carry out MCMC simulations by classical ofMisasfollows:
computers to generate data and compare quantum an- (cid:3)r
nealing with classical annealing. Section 4 proposes M= λ Q ,
a a
statistical methodologies to analyze input–output data a=1
fromaD-WavedeviceandfromvariousMCMCbased
whereλ ,...,λ aretherealeigenvaluesofM,andQ
1 r a
annealingmodels,andpresentsourfindingsonthesta-
istheprojectionontotheeigen-spacecorrespondingto
tisticalbehaviorsoftheD-WavedeviceandtheMCMC
theeigenvalueλ .Accordingtoquantumtheory,when
a
basedannealingmodels.Section5featuresconcluding
we measure the quantum system in terms of M under
remarkswithalistoffuturestatistical researchtopics.
the state |ψ(cid:2), the measurement outcome (cid:4) is a ran-
AllproofsarerelegatedintheAppendix.
domvariablethattakesvaluesin{λ ,λ ,...,λ },with
1 2 r
probability distribution P((cid:4)=λ )=tr(Q |ψ(cid:2)(cid:3)ψ|)=
a a
2. ABRIEFQUANTUMBACKGROUND
(cid:3)ψ|Q |ψ(cid:2), a =1,2,...,r. Statistically, we may per-
a
Quantum physics is counter intuitive and hard to formmeasurementsonMforthequantumsystemmul-
grasp.Itisfundamentallydifferentfromwhatwemight tiple times to obtain measurement data and infer the
expectonthebasisofeverydayexperiences.Asamat- quantum state from the data. For a quantum system at
ter of fact, Neils Bohr, a founder of quantum theory, a given time point, its state vector comprises all infor-
once said, “Anyone who thinks he understands quan- mation about the system in the sense that it can yield
tum theory doesn’t really understand it.” As this pa- theprobabilitydistributionsofmeasurementoutcomes
per involves only quantum computing based on quan- andmayencapsulateallthingsofimportanceaboutthe
tum annealing, we will bypass high-level models of quantumsystem.
quantum computation and keep our review on quan- Todescribethetimeevolutionofaquantumsystem,
tum physics and quantum computation at the minimal we denote by |ψ(t)(cid:2) the state of the quantum system
level. at time t, which is also referred to as a wave func-
tion. The states |ψ(t )(cid:2) and |ψ(t )(cid:2) at times t and
2.1 Notation 1 2 1
t are connected through |ψ(t )(cid:2) = U(t ,t )|ψ(t )(cid:2),
2 2 1 2 1
For the purpose of this paper, we consider only the where U(t ,t )=exp[−iH(t −t )] is a unitary ma-
1 2 2 1
finitedimensioncase.DenotebyCd thed-dimensional trix, and H is a Hermitian matrix on Cd. In fact, the
complexspace.Givenavectorψ inCd,wefollowthe continuoustimeevolutionof|ψ(t)(cid:2)isgovernedbythe
convention in quantum mechanics and quantum com-
Schrödingerequation
putation to use Dirac notations ket |·(cid:2) and bra (cid:3)·| to
indicate that |ψ(cid:2) and (cid:3)ψ| arecolumn and row vectors, √−1∂|ψ(t)(cid:2) =H(cid:4)(cid:4)ψ(t)(cid:5) orequivalently
respectively.Denotebysuperscripts∗,(cid:5)and†thecon- ∂t
(2.1) (cid:4) (cid:5) √ (cid:4) (cid:5)
jugate of a complex number, the transpose of a vec- (cid:4)ψ(t) =e− −1Ht(cid:4)ψ(0),
tor or matrix, and the conjugate transpose operation,
respectively. A natural inner product in Cd is given where H is a possibly time-dependent Hermitian ma-
(cid:2)
by (cid:3)u|v(cid:2) = d u∗v = (u∗,...,u∗)(v ,...,v )(cid:5), trix on Cd, which is known as the Hamiltonian of
j=1 j j 1 d 1 d
where (cid:3)u|=(u1,..√.,ud) and |v(cid:2)=(v1,...,vd)(cid:5), and the quantum system. The Schrödinger equation shows
the modulus |u|= (cid:3)u|u(cid:2). We say a matrix A is Her- thatforthequantumsystem,itsHamiltoniancancom-
mitian if A=A†, and a matrix U is said to be unitary pletely describe the dynamic evolution of its quantum
ifUU†=U†U=I,whereIisanidentitymatrix. states. It should be stressed that besides the approach
by the Schrödinger equation to describe the quantum
2.2 QuantumPhysics
evolution,thereareotherformulationsofquantumme-
A quantum system is completely characterized by chanics such as the so-called matrix mechanics cre-
its state and the time evolution of the state. A d- atedbyWernerHeisenberg,MaxBornandPascualJor-
dimensional quantum system at a given time can be dan and the path integral formulation due to Richard
366 Y.WANG,S.WUANDJ.ZOU
Feynman. See Holevo (1982), Sakurai and Napolitano d = 2b, with computational basis states of the form
(2010),Shankar(1994),Wang(2012,2013),Wangand |x x ···x (cid:2), x =0 or 1,j =1,...,b. For example,
1 2 b j
Xu(2015),Xu(2015)andCaietal.(2016). the states of two qubits are unit vectors in C4, with
the superposition states |ψ(cid:2) = α |00(cid:2) + α |01(cid:2) +
2.3 QubitsandTheirQuantumProperties 00 01
α |10(cid:2)+α |11(cid:2), where amplitudes α’s are complex
10 11
Quantum computation grapples with understanding numberssatisfying|α |2+|α |2+|α |2+|α |2=
00 01 10 11
how to take advantage of the enormous information 1. Two qubits can be made to interact and stay in an
hidden in the quantum systems and to harness the im- entangledstatesothattheymaintainpersistent,instant
mense potential of atoms and photons for the purpose influence on each other no matter how far apart they
of computing. It utilizes strange quantum phenomena
become and regardless of the nature of the interven-
suchasquantumsuperposition,quantumentanglement
ing medium. A consequence of entanglement is that
and quantum tunneling to do the trick of performing
measuring one of the two entangled qubits will make
computation and processing information. It intends to
the state of the other entangled qubit to be definite,
develop quantum computing devices for solving cer-
and thus completely predictable, even though the two
tain tough computational problems faster and/or more
qubitsmaybeveryfaraway.
efficientthanclassicalcomputers. The quantum system of b qubits is described by Cd
Any computers must utilize some states of physical
with each superposition state specified by 2b ampli-
systems to store digits. Classical computers use volt-
tudes. As 2b increases exponentially in b, it is very
agelevelstoencodebits 0 and 1.Analogtobits 0 and
easy for such a system to grow with an enormously
1 in classic computation, quantum computation uses
big space. As a result, we encounter a key signature
qubits |0(cid:2) and |1(cid:2). However, unlike classical compu-
inhandling aquantum systemthatitscomplexityusu-
tation where bits are either 0’s and 1’s and transistors
allygrowsexponentiallywithitssize.Ittakesanexpo-
can crunch the ones and zeroes individually, quantum
nential number of bits of memory on a classical com-
computation allows qubits to encode ones and zeroes
puter to store the state of a quantum system, and sim-
simultaneouslythroughwhatisknownasquantumsu-
ulations of quantum systems via classical computers
perposition.Aqubitcanberealizedbyvariousphysical
face great computational challenge. For example, be-
systems, such as the quantum spin of a particle where
|0(cid:2)and|1(cid:2)correspondtothespinupstate|↑(cid:2)andspin sides the√memory storage difficulty, we need to eval-
down state |↓(cid:2) of the particle, respectively. According uate e− −1Ht in (2.1) for the study of the quantum
to quantum physics, besides states |0(cid:2) and |1(cid:2) the par- evolution; as the matrix size of Hamiltonian H grows
ticlemayalsoexistinasuperpositionstate,whichisa exponentially with the size of the system, it is ex-
blendofspinupanddownstatessimultaneously.Thus, tremelydifficulttoexponentiatethequantumHamilto-
whileatanyinstantaclassicalbitcanbeeither 0or1, nian using mathematical analysis or classical comput-
aqubit canbeasuperposition of both |0(cid:2) and |1(cid:2), that ers.Ontheotherhand,sincequantumsystemsareable
is,asuperpositionqubitisinaquantumstatethatmay tostoreandkeeptrackanexponentialnumberofcom-
beviewedasbothoneandzeroatthesametime.Math- plexnumbersandperformdatamanipulationsandcal-
ematically, asuperposition qubit may take the form of culations asthesystemsevolve,quantum computation
|ψ(cid:2)=α |0(cid:2)+α |1(cid:2),whereα andα aretwocomplex andquantuminformationsearchforwaystoaccessthe
0 1 0 1
numbers satisfying |α0|2 +|α1|2 =1. In other words, exponential information reservoir hidden in the quan-
a qubit represents unit vectors in a two-dimensional tum systems and utilize the tremendous potential for
complexvectorspaceC2,and|0(cid:2)and|1(cid:2)consistofan computing.SeeNielsenandChuang(2000)andWang
orthonormal basis for the space and are often referred (2011,2012).
toasthecomputationalbasis.Aclassicalbitcanbeex-
amined to determine whether it is in the states 0 or 1,
3. QUANTUMANNEALINGANDMCMC
but for a qubit we can not determine its state and find
SIMULATIONS
the values of α and α by examining it. By quantum
0 1
physics, we can perform a measurement on qubit |ψ(cid:2) Annealing a material by heating and slowly cool-
andobtaineithertheresult0,withprobability|α |2,or ing it to enhance its quality is an ancient technique
0
theresult1,withprobability|α |2. that has been used for materials like glasses and met-
1
Like classic bits, we can define multiple qubits. alsoverseventhousandyears.Mimickingthisprocess
The states of b qubits are unit vectors in Cd, where by computer simulations creates simulated annealing
QUANTUMANNEALINGANDMCMCSIMULATIONS 367
as an optimization method. It is based on the anal- optimization problem with some probability. The key
ogy between the behavior of a complex physical sys- idea behind quantum annealing is to replace thermal
tem with many degrees of freedom and an optimiza- fluctuations in SA by quantum fluctuations via quan-
tionproblemoffindingtheglobalminimumofagiven tum tunneling so that the system is kept close to some
objective(orcost)functiondependingonmanyparam- instantaneousgroundstateofthequantumsystemdur-
eters. By viewing the objective function of the opti- ing the quantum annealing evolution, analog to some
mization problem as the energy of the physical sys- quasi-equilibriumstatetobekeptduringthetimeevo-
tem, we naturally formulate the optimization problem lutionofSA.
asthe(possibleNP-hard)problemoffindingminimum Bothclassicalandquantumannealingtechniquesare
energy configurations (or ground states) of the many- powerfultoolsforsolvinghardoptimizationproblems,
body physical system, and develop computer simu- whether they are utilized as physical devices or sim-
lation algorithms to mimic the system behavior and ulation methods. The physical scheme is to employ a
search for the minimum energy configurations of the natural system or build a device to engineer a physi-
correspondingphysicalmodel. cal system so that ground states of the system repre-
A classical approach is simulated annealing (SA), sent the sought-after solution of an optimization prob-
which takes into account the relative configuration lem (McGeoch, 2014). The simulation approach is to
energies and a fictitious time-dependent temperature apply “escape” rules in computer simulations to pre-
when probabilistically exploring the immense search vent the system from getting trapped in local minima
space. For a given optimization problem, its objective of an energy or cost function, and eventually reach
function to be minimized is identified with the energy theglobalminimumwithsomeprobability(Riegerand
ofaphysicalsystem,andwethengivethephysicalsys- Kawashima, 1999 and Martonˇák, Santoro and Tosatti,
tem a temperature as an artificially-introduced control 2002).Inbothsituations,thesystemisallowedtoprob-
parameter.Byreducingthetemperaturegraduallyfrom abilistically explore its immense configuration space
a high value to zero, we wish to drive the system to and ultimately “freeze” in the global minimum with
the state with the lowest value of the energy (objec- certain probability, and by enough repeated tries we
tive function), and thus reach the solution of the op- can find the global minimum and solve the optimiza-
timization problem. The initial temperature is usually tionproblem.
set high relative to the system energy scale in order
3.1 ClassicalIsingModelandSimulatedAnnealing
to induce thermal fluctuations and sample its config-
urationsusingMCMCsimulationsliketheMetropolis The Ising model is often used to describe natural
algorithm, and the annealing process escapes from lo- systems in physics, and many optimization problems
cal minima by thermal fluctuations to climb over en- can be mapped into physical systems described by the
ergy barriers and searches for lower energy configura- Isingmodelwhosegroundstatesprovidethesolutions
tions. As the system evolves with the change of tem- totheoptimizationproblems.Examplesincludetravel-
perature sufficiently slow, the system is expected to ing salesman problem, portfolio optimization, integer
stay close to thermal equilibrium during time evolu- factoring, social economics network, protein folding,
tion, and thus in the end we lead the SA system to the protein modeling and statistical genetics. See Irback,
zero-temperature equilibrium state, the lowest-energy Peterson and Potthast (1996), Majewski, Li and Ott
state.SeeBertsimasandTsitsiklis(1992),Kirkpatrick, (2001),McGeoch(2014)andStauffer(2008).
GelattandVecchi(1983)andWinker(2001). Consider the Ising model described by a graph G =
Quantumannealingisbasedonthephysicalprocess (V(G),E(G)),where V(G) and E(G) standforthever-
of a quantum system whose lowest energy, or equiv- texandedgesetsofG,respectively.Eachvertexisoc-
alently, a ground state of the system, represents the cupiedbyarandomvariabletakingvaluesin{+1,−1},
solution to an optimization problem posed. It starts and each edge represents the coupling (or interaction)
with building a simple quantum system initialized in between the two vertex variables connected by the
its ground state, and then moves the simple system edge. A configuration s={s ,j ∈V(G)} is defined to
j
graduallytowardthetargetcomplexsystem.According be an assignment of a set of values to all vertex vari-
to the quantum adiabatic theorem (Farhi et al., 2000, abless ,j ∈V(G).Verticesarealsoreferredtoassites,
j
2001, 2002), as the system slowly evolves, it tends to and vertex variables as spins in physics, with +1 for
remaininagroundstate,andhencemeasuringthestate spinupand−1forspindown.Forexample,considera
of the final system will yield an answer to the original graph corresponding to a two-dimensional lattice with
368 Y.WANG,S.WUANDJ.ZOU
a magnet placed at each lattice site pointing either up of Hc(s). SA finds the minimum with repeated tries
I
or down. Suppose that there are b lattice sites labelled by using MCMC methods such as the Metropolis–
by j = 1,...,b. Site variable s stands for a binary Hastings algorithm to generate configuration samples
j
random variable indicating the position of the magnet from the Boltzmann distribution with temperature de-
atsite j,where s =±1 maybeinterpretedasthe jth creasing slowly. It starts with a random initial spin
j
magnetpointingupordown,respectively. configuration and randomly flips spins at each time
TheHamiltonianoftheclassicalIsingmodelisgiven step; we always accept a new spin configuration if it
by lowers the energy and accept it probabilistically us-
(cid:3) (cid:3)
ing the Metropolis rule otherwise; and we gradually
(3.1) Hc(s)=− J s s − h s ,
I ij i j j j lowerthetemperaturetoreducetheescapeprobability
(i,j)∈E(G) j∈V(G)
oftrappinginlocalminima.Specifically,theSAalgo-
where (i,j) stands for the edge between sites i and rithm initializes spins with −1 and+1 randomly and
j, with the first sum over all pairs of vertices with independently from each other to obtain initial spins
edge (i,j) ∈ E(G), Jij stands for the interaction (or s(0) = {s(0)}. We update spins one by one, and each
j
coupling) between sites i and j associated with edge completeupdatingoverallspinsconstitutesonesweep.
(i,j) ∈ E(G), and hj describes an external magnetic Atthekthsweep,forspini,weattempttoflipitsstate
fieldonvertexj ∈V(G).Asetoffixedvalues{Jij,hj} s(k−1) to new state s(k) = −s(k−1) while keeping all
i i i
is referred to as one instance of the Ising model. For a
other spins unchanged, and calculate energy change
given configuration s, the energy of the Ising model is (k−1)
between its original state s and the newly flipped
equaltoHc(s). i
I (k)
states ,
According to Boltzmann’s law, the probability of a i
(cid:6) (cid:7)
given configuration s is described by the Boltzmann (cid:8)E(k)=−h s(k)−s(k−1)
(orGibbs)distribution i i i i
Pβ(s)= e−βHcI(s), Zβ =(cid:3)e−βHcI(s), −(cid:3)i−1Jijsj(k)(cid:6)si(k)−si(k−1)(cid:7)
Zβ s j=1
whereβ =(kBT)−1 isaninversetemperature,withkB − (cid:3)b J s(k−1)(cid:6)s(k)−s(k−1)(cid:7).
a generic physical constant called the Boltzmann con- ij j i i
stant and T the absolute temperature, and the normal- j=i+1
ization constant Zβ is called the partition function. If The new state s(k) is accepted with probability min{1,
we take kB = 1, T is interpreted as the fundamental exp(−(cid:8)E(k)/Ti)}, that is, we change spin i’s state
temperature of the system with units of energy, and β i k
from s(k−1) to new state s(k) if (cid:8)E(k) ≤0 and other-
istreatedasthereciprocaltothefundamentaltempera- i i i
ture.TheconfigurationprobabilityPβ(s)representsthe wiseupdateitsstatewithprobabilityexp(−(cid:8)Ei(k)/Tk),
probability that the physical system is in a state with whereTk istheannealingscheduletolowerthetemper-
configurationsinequilibrium. ature. Typical annealing schedules for lowering tem-
When a combinatorial optimization is represented perature include Tk proportional to k1 or lo1gk. See
by the Ising model, the objective is to find a ground Bertsimas and Tsitsiklis (1992), Geman and Geman
state of the Ising model, that is, we need to find a (1984)andHajek(1988).
state whose configuration minimizes the energy func-
3.2 QuantumAnnealing
tion Hc(s). If the Ising model has b sites, the config-
I
uration space is {−1,+1}b with the total number of WedescribeaquantumIsingsystemassociatedwith
theconfigurationsequalto2b.Forasystemwithmany quantumannealingthroughthesamegraphG usedfor
sites,becauseoftheexponentialcomplexity,itisanex- the classical Ising model, where similar to the clas-
tremelydifficulttasktonumericallyfindgroundstates sical case, the vertex set V(G) now represents quan-
and solve the minimization problem. In fact, it is pro- tum spins, with edge set E(G) for couplings (or inter-
hibitive fordeterministic exhaustive searchalgorithms actions) between two quantum spins. As illustrated in
to solve the minimization problem with such an ex- Section 2.3, each quantum spin may realize a qubit,
ponential growth search space. Annealing approaches and thus each vertex represents a qubit, with the to-
suchasSAareoftenemployedtoprobabilisticallyex- tal number of qubits in the system equal to the total
plore the search space and find the global minimum number, b, of the vertices in G. The quantum system
QUANTUMANNEALINGANDMCMCSIMULATIONS 369
is characterized by complex space Cd, where d =2b, identicalmatricesandtensorproductsigns.Becauseof
withitsquantumstatedescribedbyaunitvectorinCd thetensorproducts,alltermsinHq arediagonalmatri-
I
anditsdynamicevolutiongovernedbytheSchrödinger cesofsize2b,andsodoesHq.
I
equation defined in (2.1) via a quantum Hamiltonian, The goal of finding a ground state of quantum
q q
whichisaHermitianmatrixofsize d.Theenergiesof Hamiltonian H is to search for an eigenvector of H
I I
the quantum system correspond to the eigenvalues of corresponding to its smallest eigenvalue, or equiva-
thequantumHamiltonian,withagroundstatebeingan lently, a quantum spin configuration with the minimal
eigenvector corresponding to the smallest eigenvalue. energy. Because Hq involves only tensor products of
I
TospecifythequantumHamiltonianweneedtodefine commutingdiagonalmatricesI andσ ,Hq isadiag-
j j I
(cid:8) (cid:9) (cid:8) (cid:9)
onal matrix with eigenvalues equal to its diagonal en-
1 0 0 1
Ij = 0 1 , σxj = 1 0 , tries.Moreover,itiseasytoshowfromthesamestruc-
q
ture of (3.1) and (3.2) that the eigenvalues of H are
(cid:8) (cid:9) I
σz = 1 0 , j =1,...,b, actually the 2b values of classical Hamiltonian HcI(s),
j 0 −1 s∈{−1,+1}b. Thus, the system governed by Hq be-
I
where σx and σz are called Pauli matrices in x and z haves essentially like a classical system, and finding
j j
the minimal energy of the quantum Ising Hamiltonian
axes, respectively (here we do not need Pauli matrix
q
H is equivalent to finding the minimal energy of the
in y axis, see Nielsen and Chuang, 2000 and Wang, I
classicalIsingHamiltonianHc.
2012). To describe the quantum system, we replace I
each classical vertex variable s =±1 in (3.1) by σz Sincetheoriginaloptimizationproblemdescribedin
j j
forquantumspin j.Paulimatrix σz hastwoeigenval- Section 3.1 can be formulated in the quantum frame-
j work,thecomputationaltaskforsolvingtheoptimiza-
ues ±1 whose eigen-states |+1(cid:2) and |−1(cid:2) correspond
to spin up state |↑(cid:2) and spin down state |↓(cid:2), respec- tion problem remains the same as in the classical case
so far. The key in quantum annealing is to engineer
tively, so that a qubit at vertex j is realized by the
a magnetic field orthogonal to the Ising axis and ob-
quantum spin. Quantum configurations consist of the
2b possible combinations of 2b eigen-states |±1(cid:2) of tain the Hamiltonian of the Ising model in the trans-
q
the Pauli matrices {σz}b . Replacing s in classical verse field. With HI representing a potential energy,
j j=1 j thetransversefieldstandsforakineticenergythatdoes
Ising Hamiltonian HcI(s) by σzj, we obtain the follow- not commute with Hq, thus it induces transitions be-
ingquantumHamiltonianofthequantumIsingmodel: I
(cid:3) (cid:3) tween the up and down states of each single spin, and
(3.2) Hq =− J σzσz − h σz, turns the model behavior from classical to quantum.
I ij i j j j
(i,j)∈E(G) j∈V(G) Specifically, assume that the transverse magnetic field
where J is the Ising coupling along the edge (i,j)∈ isgovernedbyaquantumHamiltonian
ij (cid:3)
E(G), and hj is the local field on the vertex j ∈V(G). (3.3) HX =− σxj,
Here, we use the convention in the quantum literature j∈V(G)
that σz and σzσz in (3.2) denote their tensor products
j i j whereweagainfollowthequantumconventionthatσx
alongwithidenticalmatricesinsuchawaythat j
stands for tensor products of b matrices of size 2, that
σziσzj ≡I1⊗···⊗Ii−1 is,
⊗σ(cid:10)zi ⊗Ii+1⊗·(cid:11)·(cid:12)·⊗Ij−1⊗σzj(cid:13) σxj ≡I1⊗···⊗Ij−1⊗(cid:10) σ(cid:11)(cid:12)xj⊗(cid:13)Ij+1⊗···⊗Ib,
verticesi andj vertexj
⊗Ij+1⊗···⊗Ib, whichdoesnotcommutewithσz inHq.Sinceσx has
j I j
σzj ≡I1⊗···⊗Ij−1⊗(cid:10) σ(cid:11)(cid:12)zj⊗(cid:13)Ij+1⊗···⊗Ib. |twvjo,±e1i(cid:2)g=env(1a,lu±e1s)±†,1HwXithhascosirmrepspleoenidgienngveecigtoernv|ve+ct(cid:2)o=rs
vertexj |v1,+1(cid:2)⊗|v2,+1(cid:2)⊗···⊗|vb,+1(cid:2) corresponding to its
That is, we have one matrix, either a Pauli matrix σz smallest eigenvalue, namely, |v+(cid:2) is the ground state
j
or an identity matrix I , to act on the jth qubit, each ofH .
j X
term in (3.2) is a tensor product of b matrices of size Thequantumannealingprocedurestartswithaquan-
two,and thequantum convention simplyidentifies the tumsystemdrivenbythetransversemagneticfieldH
X
qubitswithPaulimatricesforrealactionsbutomitsthe and initialized in its ground state. During the process
370 Y.WANG,S.WUANDJ.ZOU
ofquantumannealing,thesystemisevolvedgradually degrees of freedom of quantum nature, and then inge-
from initial Hamiltonian H to the final target Hamil- niouslycontrolsthestrengthofthesequantumfluctua-
X
tonian Hq. For such Hamiltonian changes, the quan- tionsviaannealingschedulesA(t)andB(t)inorderto
I q
tum adiabatic theorem indicates that the system tends makethesystemfinallyreachagroundstateofH ,just
I
toremainingroundstatesoftheinstantaneousHamil- like the slow reduction of temperature in classical an-
tonian through quantum tunneling (Farhi et al., 2000, nealing. Initially, the strength of quantum fluctuations
2001,2002,andMcGeoch,2014).Attheendofquan- is set to be relatively high for the system to search
for the global structure of the configuration space,
tum annealing, if the system is in a ground state of
q corresponding to the initial high-temperature situation
the final Hamiltonian H , an optimal solution is ob-
I
in SA, and then the strength is slowly decreased to
tainedbymeasuringthesystem.Specifically,quantum
finally vanish to recover the original target system
annealing is realized by an instantaneous Hamiltonian
hopefullyinthelowest-energystate.Quantumanneal-
fortheIsingmodelinthetransversefieldasfollows:
ing differs from the classical annealing in terms of
(3.4) H (t)=A(t)H +B(t)Hq, t ∈[0,t ], quantum-mechanical fluctuations in quantum anneal-
D X I f
ing instead of thermal fluctuations in SA. While clas-
whereA(t)andB(t)aretime-dependentsmoothfunc-
sical annealing relies on thermal fluctuations to make
tions controlling the annealing schedules, and tf is thesystemhopfromstatetostateoverintermediateen-
the total annealing time. Typically, A(tf)=B(0)=0, ergy barriers and search for the desired lowest-energy
A(t) is decreasing and B(t) is increasing. The grad- state, quantum annealing replaces thermal hopping by
ualmoveoftheHamiltonianfromHD(0)toHD(tf)is quantum-mechanical fluctuations for state transitions.
achieved through controlling the annealing schedules Thequantumfluctuationsinquantumannealingarere-
A(t)andB(t).Atinitialtimet =0,HD(0)=A(0)HX, alized via quantum tunneling that allows the anneal-
at the final time t =t , H (t )=B(t )Hq, and A(0) ing process to explore the different states by traveling
f D f f I
andB(t )areknownscalars,thenH (t)hasthesame directly through energy barriers, rather than climbing
f D
q
eigenvectors as H attheinitial timeandas H atthe over them thermally. Quantum tunneling refers to the
X I
final time, where the corresponding eigenvalues differ quantum phenomenon where particles tunnel through
by factors of A(0) and B(t ), respectively. Therefore, abarrierinthecircumstancethatcouldnotbepossible
f
the quantum annealing evolution driven by H (t) es- by classical physics. The tunneling process cannot be
D
sentiallymovesthesystemH initializedatitsground directly observed nor adequately explained by classi-
X
statetothefinalsystemHq.Accordingtothequantum cal physics. Quantum tunneling is often explained us-
I
ingtheHeisenberguncertaintyprincipleandthewave-
adiabatic theorem (Aharonov et al., 2007, Born and
particle duality of matter in quantum physics (Sakurai
Fock, 1928, Farhi et al., 2000, 2001, 2002, McGeoch,
andNapolitano,2010andShankar,1994).Here,wetry
2014, and Morita and Nishimori, 2008), for appropri-
toprovidesomeheuristicexplanation.Tofacilitateour
ately chosen A(t) and B(t) we can measure the quan-
understanding of the phenomenon, we may compare
tum system at the final annealing time t to find a
f
q particles attempting to travel between barriers to balls
ground state of H with sufficiently high probability.
I trying to roll over a hill or balls trying to penetrate
Thus, with certain probability the quantum annealing
a wall. Classical physics indicates that balls without
procedure driven by (3.4) can obtain the global min-
sufficient energy to surmount the hill would roll back
imum of the objective function Hc(s) and solve the
I down, and balls which do not have enough energy to
minimizationproblem.
penetrate the wall would either bounce back or bury
Inquantumannealing,themovefromtheinitialsim-
themselvesinsidethewall.Inbothcases,ballswillnot
q
plesystemH tothefinalcomplexsystemH isoften
X I be able to reach the other sides of the hill and wall.
realized through engineering magnetic fields. For ex-
Whiletheseparticleswithoutenoughenergytoclassi-
ample, we may manufacture the move from the initial
cally surmount a barrier cannot get to the other side
system to the final system by turning magnetic fields of the barrier by classical physics, quantum physics
on and then off adiabatically, like the physical imple- predicts that the particles can, with some probability,
mentation by D-Wave devices that initializes its sys- cross the barrier and reach the other side, thus tunnel-
tem in the ground state of HX by chilling the sys- ing through the barrier. Figure 1 provides some car-
tem to a near absolute zero temperature and turning toon illustrations of thermal hopping versus quantum
on the transverse magnetic field. Quantum annealing tunneling.Theillustrationsmayindicatethatifenergy
induces quantum fluctuations by introducing artificial barriers are high and thin, quantum tunneling may be
QUANTUMANNEALINGANDMCMCSIMULATIONS 371
FIG.1. Acartoonillustrationofquantumtunnelingvs.thermalclimbingonthetoppanelwithannealingelucidationsofquantumtunneling
ontheleftbottompanelandthermalclimbingontherightbottompanel.
more efficient than thermal hopping, and quantum an- 3.3 PhysicalImplementationofQuantum
nealing may be faster than classical thermal annealing AnnealingbytheD-WaveQuantumComputer
to reach the equilibrium distribution. In fact, research
The D-Wave quantum computer is a commercially
has shown that it can be more efficient to explore the
available computing hardware device designed to im-
state space quantum mechanically, and quantum an-
plement quantum annealing. As illustrated in Fig-
nealing may have advantages over thermal annealing
ure 2, it includes a quantum processor, a magnetically
(Brookeetal.,1999,Denchevetal.,2016,Farhietal.,
shielded box and a cooling equipment housed inside
2000, 2001, 2002, Kechedzhi and Smelyanskiy, 2015
a ten square meter shielded room. The quantum pro-
andMcGeoch,2014).
cessor consists of superconducting flux qubits based
Quantum annealing may be implemented by phys- onniobiumsuperconductingquantuminterferencede-
ical devices or MCMC simulations. We will illus- vices,andtheshieldedboxandthecoolingequipment
trate in next two sections D-Wave physical devices protect the system from coupling to the environment
to implement the unitary dynamic evolution of quan- and maintain the system at near absolute zero temper-
tum annealing by natural Schrödingier dynamics and ature to keep superconducting flux qubits in quantum
MCMC methods to approximate the dynamic evolu- states.D-WaveOne“Rainier”deviceof128qubitsand
tionofquantumannealingbyartificialtimeevolutions D-WaveTwo“Vesuvius”deviceof512qubitswerere-
ofMonteCarlodynamics. leasedinMay2011andMay2013,respectively,while