Emulating Quantum Interference with Generalized Ising Machines
EEmulating Quantum Interference with Generalized Ising Machines
Shuvro Chowdhury, Kerem Y. Camsari,
1, 2 and Supriyo Datta School of Electrical and Computer Engineering, Purdue University, IN 47907, USA Currently at the Department of Electrical and Computer Engineering,University of California Santa Barbara, Santa Barbara, CA 93106, USA (Dated: July 16, 2020)The recent groundbreaking demonstration of quantum supremacy in noisy intermediate scalequantum (NISQ) computing era has led to an intense activity in establishing finer boundaries be-tween classical and quantum computing. In this paper, we use established techniques based onquantum Monte Carlo (QMC) to formulate a systematic procedure for translating any sequence of d quantum gates acting on n q-bits into a Boltzmann machine (BM) having n + g ( d ) classical spinsor p-bits with two values “0” and “1”, but with a complex energy function E . Using this procedurewe emulate Shor’s algorithm with up to 36 q-bits using 90 p-bits, on an ordinary laptop computer inless than a day, while a naive Schr¨odinger implementation would require multiplying matrices with ≈ elements. Even larger problems should be accessible on dedicated Ising Machines. However,we also identify clear limitations of the probabilistic approach by introducing a quantitative metric S Total for its inefficiency relative to a quantum computer. For example, a straightforward proba-bilistic implementation of Shor’s algorithm with n q-bits leads to an S Total ∼ exp ( − n/ n/ instead ofthe polynomial scaling expected for true quantum computers. This is because quantum algorithmslike Grover’s search and Shor’s factorization lead to BM’s where specific outputs are suppressednot by the usual statistical factor exp ( −(cid:60) ( E )), but by a cancellation of paths through the phaseexp ( − i (cid:61) ( E )), a manifestation of the well-known “sign problem” in QMC. But it may be possible to“tame” this problem with appropriate transformations that introduce real parts into the energy andincrease S Total . Finally, we present an example featuring a standard optimization algorithm basedon a purely real energy function to which we add an imaginary part (cid:61) ( E ), thereby augmentingthe statistical suppression of Feynman paths with quantum-like phase cancellation. This exampleillustrates how the sign problem encountered in classical annealers can be turned into a computa-tional resource for quantum annealers as the augmented algorithm can potentially be implementedin hardware quantum annealers. I. INTRODUCTION
Quantum computing is based on the use of quantumgates to perform d successive unitary transformations(gates) U (1) , U (2) , · · · , U ( d − , U ( d ) on a set of n q-bitsso that their wavefunctions evolve from an initial | ψ (0) (cid:105) to a final | ψ ( d ) (cid:105) (Fig. 1a) [1]. Each of these wavefunc-tions | ψ ( d ) (cid:105) has 2 n components, coming from a tensorproduct of n single q-bit wavefunctions with two com-plex components each. In classical computing, a directdeterministic calculation requires us to multiply 2 n × n transformation matrices with exponentially large mem-ory requirements as n increases. By contrast, a quantumcomputer requires only n q-bits which naturally live in2 n × n Fock space.Powerful algorithms such as Shor’s period-finding andGrover’s Search Algorithm [2, 3] that creatively makeuse of quantum interference in the 2 n × n dimensionalFock space, led to an intense activity in recent years todevelop quantum computers [4–6], along with parallel ef-forts that have pushed the boundaries of classical simula-tion of quantum systems [7–13]. However, all such effortsof classical emulation are limited to special cases and it iswidely believed that the in the most general cases, quan-tum circuits can only be emulated by exponential costsin either memory, time or both [4].In this paper we first establish a common platform link- ing quantum and probabilistic circuits to highlight theirdifferences and facilitate the cross-fertilization of ideas.In the spirit of earlier works [4, 14–20], we use a Feynmanpaths approach to establish a systematic methodology fortranslating any sequence of d operations acting on n q-bits into a Boltzmann machine (BM) having ( n + g ( d ))classical spins or p-bits having two values “0” and “1”.The first n p-bits represent the input q-bits, while theother p-bits represent the q-bits after the application ofsuccessive gating operations (Fig. 1b).There is, however, a key difference with the typicalreal valued energy functions in traditional BMs and theenergy functions we use in this work. Although the p-bits are classical with values that can only be “0” or “1”,the energy function E = E (1) + E (2) + · · · + E ( d ) canhave complex values. The time evolution of the BM isgoverned by the real part (cid:60) ( E ) and can be emulated withstandard sampling techniques like Gibbs sampling. Butunlike standard BM, individual samples make complexcontributions of unit magnitude given by exp ( − i (cid:61) ( E )),which on summing, approaches the exact wavefunction ofthe corresponding q-circuit, similar in spirit to existingQuantum Monte Carlo (QMC) based approaches [21].Similar concepts have been discussed in the context ofadiabatic quantum computing (AQC) which is based onan evolution operator of the form exp ( − βH ), H beingthe full Hamiltonian matrix. These ideas have also been a r X i v : . [ c ond - m a t . d i s - nn ] J u l FIG. 1 . Mapping between quantum circuits andBoltzmann machines: (a) In quantum computing, a se-quence of unitary gates U (1) , U (2) , · · · , U ( d ) is applied to aninitial wavefunction | ψ (0) (cid:105) and a final wavefunction | ψ ( d ) (cid:105) isobtained. (b) the quantum circuit in (a) is mapped intoa Boltzmann machine with p-bits. Note that although thequantum gates in Fig. 1a are time-ordered, the p-bit networkin Fig. 1b is a reciprocal BM described with an energy func-tion which is obtained by modeling each k -th quantum gatevia a corresponding complex energy function, E ( k ) . Never-theless the time sequence of gates is reflected in the BM’senergy function. While in quantum case the number of q-bitsremains the same, in the Boltzmann analog, the number ofp-bit required increases with the number of applied gates d as n + g ( d ) (see text for a discussion on the nature of g ( d )). applied to gated quantum computing (GQC) [4] basedon unitary transformations U ∝ exp ( − iH ) implementedwith man-made quantum gates.Recently, Boltzmann machines with complex weightshave also been employed with machine learning tech-niques to various classes of quantum Hamiltonians. It hasbeen shown that certain representations of Boltzmannmachines can be trained to obtain the ground state of aHamiltonian, and can even be used to represent quantumstates exactly [22–25]. A key difference between thiswork and the methods presented in this paper is thatthere is no training of weights in the complex valued BMsin our approach and the weights are obtained what mightbe called “one-shot” learning as we describe in Section II. (b) (c) H H H d = 1 d = 2 +1+1 +1-1 + + (constructive interference)(destructive interference) (a) Simple double slit experiment
000 001 010 011 100 101 110 111 +ve-ve C o n t r i b u t i o n p e r s a m p l e p-bit path configuration S Output q-bit configuration, at d = 2
FIG. 2 . Quantum interference with probabilistic sam-pling: (a) One q-bit with a string of Hadamard (H) gatesapplied in sequence. The input is clamped to | (cid:105) which canevolve into different output states | (cid:105) and | (cid:105) through differ-ent intermediate paths. (b) The contribution of each sampleto four possible paths at d = 2 (the other four paths are notvisited because those require the input to be set at | (cid:105) ). Allpaths contribute equally in terms of magnitude but the 011path has a phase (hatch filled) which is opposite to that of theother three paths. (c) Average sign (see Eq. (4) in the text)plotted against different outputs which reflects that half ofthe total samples are wasted to get the “null” at | (cid:105) . Double slit interference:
The basic ideas behindour approach can be appreciated with a simple example.Consider one q-bit driven by a sequence of Hadamardgates (Fig. 2a), each represented by a transformation con-necting q-bits in planes d − d : U Hadamard = (cid:18) √ (cid:19) (cid:20) | ( d − (cid:105) | ( d − (cid:105)(cid:104) ( d ) | +1 +1 (cid:104) ( d ) | +1 − (cid:21) (1)Two applications of the gate result in an identity trans-formation:[ U Hadamard ] = (cid:20) | ( d − (cid:105) | ( d − (cid:105)(cid:104) ( d ) | (cid:104) ( d ) | (cid:21) (2)If the q-bit is initialized to | (cid:105) , then after one gate it hasequal probability of being either in | (cid:105) or | (cid:105) , but aftertwo gates the q-bit is back to the | (cid:105) state with 100%probability.This is a very simple illustration of the classic doubleslit interference that Feynman used in his lectures to il-lustrate the difference between quantum and classical, orbetween electrons and bullets as he put it [26]. As wewill see, our Boltzmann machine emulates this quantuminterference using a complex energy function given by E = iπ ( s s + s s ) + constant (3)where s , s , s represent the initial state, the state af-ter one gate and the state after two gates respectively.Starting from | (cid:105) , the q-bit will also end up in | (cid:105) afterthe two Hadamard gates. But the bullet starting froma “0” (corresponding to s = 0) can end up in each ofthe final states “0” (corresponding to s = 0) and “1”(corresponding to s = 1) via two paths (each path isdenoted by s → s → s ) as follows0 → → → → → → → → (cid:60) ( E ) and appear with equalprobability in the sampling process giving the same resultfor both outputs if the imaginary part is ignored as shownin Fig. 2b. But the complex Boltzmann machine weightseach path according to exp ( − i (cid:61) ( E )) so that the totalcontribution isFinal state in 0: → e − iπ (0+0) + e − iπ (0+0) = 2Final state in 1: → e − iπ (0+0) + e − iπ (0+1) = 0 . Our complex Boltzmann machine thus emulates quan-tum interference using classical bullets by associatingcomplex energies with the paths taken by the bullet andhence in principle provides results that are exact whenall possible paths taken by a bullet (between a given in-put state to a given output state) are considered intoaccount.
Sign problem:
But the p-bits with complex phasesare still fundamentally very different from q-bits whichwould all end up in the correct state after two idealHadamard gates. By contrast our p-bits can take fourpaths each with probability 1 /
4. Two of the paths addup to give 1 /
2, while the other two cancel to give zeroand are in a sense wasted which is shown in Fig. 2c. Thisis a problem that becomes more acute as the depth of thecircuit is increased. With two Hadamard gates the heightof the correct peak is 1 /
2; with d ( d being even) gatesit is 1 / d/ so that an acceptable peak to null ratio willrequire a larger number of samples. This is essentiallythe sign problem well-known in QMC [27].Quantitatively this problem can be characterized bythe average sign S α defined for a given output configura-tion α as
2 8 14 20
Number of gates, d S = = 0= /6= /8= /5 0 /4 /20.0 0.250.50 FIG. 3 . Improvement in average sign with rotatedHadamard gates:
Average sign for the correct peak ( α = 0)with input “0”) is plotted against the number of gates forvarious rotation angles, θ . Note that the number of gatesare even so that for a given number of gates the output ofthe ordinary and the rotated Hadamard gates are exactly thesame. Each of the dashed lines follows a 2 − ζd fit. Inset showsthe dependence of ζ (periodic with a period of π/
2) withrespect to the rotation angle θ . S α = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:80) paths to α e −(cid:60) ( E ) e − i (cid:61) ( E ) (cid:80) γ (cid:80) paths to γ e −(cid:60) ( E ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (4)One way to characterize the efficiency of the prob-abilistic emulation of quantum circuits is in terms ofthe total sign by summing over average signs S α forall α , i.e., S total = (cid:80) S α with the number of samples N T needed for a given accuracy scaling inversely as N T ∝ /S total . At one extreme we have stoquasticproblems where the BM energy is purely real giving S total ≈
1, and probabilistic emulation can be quiteefficient. At the other extreme we have classic quantumalgorithms like Grover’s search and Shor’s algorithmwhere the BM energy is purely imaginary, and 1 /S total increases exponentially with the size of the problem.Such problems can be basis transformed to reduce1 /S total , a process that is commonly referred to astaming the sign problem in the QMC literature. A toyexample of this is presented next. Taming the sign problem:
It is apparent fromEq. (4) that if (cid:61) ( E ) = 0 the average signs of all outputstates should add up to 1, indicating the absence of anysign problem, while the present example represents theother extreme where (cid:60) ( E ) = constant, so that all pathsare equally likely and any discrimination between differ-ent outputs occurs through their average signs. Muchwork has gone into taming [28–30] the sign problem soas to transfer the path discrimination at least partially to (cid:60) ( E ). In the present case for example, with d = 12, thesign for the peak equals 1 / = 0 . V U
Hadamard V † ) torotate the Hadamard gate by an angle of θ around the y -axis ( V = R y ( θ )) to transform it into U H,rotated = (cid:20) | ( d − (cid:105) | ( d − (cid:105)(cid:104) ( d ) | cos ( θ + π/
4) sin ( θ + π/ (cid:104) ( d ) | sin ( θ + π/ − cos ( θ + π/ (cid:21) (5)then it is possible to somewhat improve the averagesign by choosing an appropriate θ as it is shown in Fig. 3. Organization of the paper:
In Section II we de-scribe how we obtain our rules for translating one q-bitand two q-bit gates into a complex energy function E gov-erning the corresponding complex Boltzmann machine,using the Feynman path approach [4]. In Section IIIwe present results for a probabilistic implementation ofa two q-bit Grover search circuit illustrating how thenumber of paths reaching each output is essentially thesame, and nulls arise due to a cancellation rather than a suppression of separate paths. Like the example above,here too (cid:60) ( E ) is trivial so that probabilistic samplingdoes not provide any output selectivity, and the Groverpeak arises entirely from constructive addition which isreflected in its large average sign compared to the sup-pressed outputs. In Section IV we present results for aprobabilistic implementation of Shor’s period-finding cir-cuit, which too gives rise to peaks through cancellationrather than suppression of paths. We present scalingresults for circuits with n = 20 , , , ,
36 q-bits em-ulated with 50 , , , ,
90 p-bits respectively showingthat the number of samples needed to establish a fixedpeak to null ratio increases as 2 n/ unlike the polyno-mial scaling expected for a q-bit implementation. The36 q-bits system is emulated with 90 p-bits and has beenperformed on an ordinary laptop in less than a day, acalculation that would require us to multiply matriceswith ≈ elements in a deterministic approach. Fi-nally in Section IV we present an example featuring astandard optimization algorithm based on a purely realenergy function to which we add an imaginary part (cid:61) ( E ),thereby augmenting the statistical suppression of Feyn-man paths with quantum-like phase cancellation. II. COMPLEX ENERGY FUNCTIONS FORQUANTUM GATES
The basic rule for writing a complex energy functionfor a sequence of gate operations can be obtained as fol-lows. First we note that the elements of the overall trans-formation matrix are given by a matrix product of theindividual matrices U α,β = (cid:88) p , ··· ,p d − U ( d ) α,p d − U ( d − p d − ,p d − · · · U (2) p ,p U (1) p ,β (6) where | β (cid:105) and | α (cid:105) corresponds to the initial and the fi-nal state of the q-bits respectively. We express each k -th element of the transformation matrices in terms of acorresponding energy function, which in general can becomplex: U ( k ) p,q = e − E ( k ) p,q ⇒ E ( k ) p,q = − ln (cid:16) U ( k ) p,q (cid:17) (7)Using Eq. (7) in Eq. (6) , we can write U α,β = (cid:88) p , ··· ,p d − e − (cid:0) E ( d ) α,pd − + E ( d − pd − ,pd − + ··· + E (2) p ,p + E (1) p ,β (cid:1) = (cid:88) p , ·· ,p d − e −(cid:60) ( E ) e − i (cid:61) ( E ) (8)where (cid:60) ( E ) and (cid:61) ( E ) represent the real and imaginaryparts of the total energy function E obtained by summingthe individual ones E ( p , · · · , p d − ) = E ( d ) α,p d − + E ( d − p d − ,p d − + · · · + E (1) p ,β (9)Eq. (8) is an exact result, but it represents a sum over alarge number of paths β → p → p → · · · → p d − → α (from an initial state α to a final state β ) which growsexponentially in n .Usually the energies are real so that e − E can beinterpreted as a probability and probabilistic approachesallow us to sample the most important paths based onpowerful algorithms like Metropolis or Gibbs sampling.For complex energies we can still sample the pathsbased on the real part of E while the imaginary partcan be interpreted as the complex contribution of unitmagnitude contributed by a particular path. In the nextwe will show how to obtain energy functions for one andtwo q-bit gates and then we will present the algorithmwhich allows us to include complex contribution ofvarious paths with the usual sampling procedure. One q-bit gates:
Any one q-bit gate is in generaldescribed by a transformation matrix of the form U ( d ) = (cid:20) | ( d − (cid:105) | ( d − (cid:105)(cid:104) ( d ) | a b (cid:104) ( d ) | c A (cid:21) (10)Using Eq. (7) we can immediately write the energyfunction in tabular form E ( d ) = (cid:20) s ( d − =0 s ( d − =1 s ( d ) =0 − ln ( a ) − ln ( b ) s ( d ) =1 − ln ( c ) − ln ( A ) (cid:21) It is straightforward to convert this tabular result intoa Boolean sum of products expression [31] E ( d ) = − (1 − s ( d ) )(1 − s ( d − ) ln ( a ) − (1 − s ( d ) ) s ( d − ln ( b ) − s ( d ) (1 − s ( d − ) ln ( c ) − s ( d − s ( d ) ln ( A ) (11)which simplifies to E ( d ) = − ln ( a ) + s ( d ) ln ( a /c ) + s ( d − ln ( a /b )+ s ( d − s ( d ) ln ( b c /a A ) (12)Note that this energy function has linear and quadraticterms corresponding one-body and two-body interactionsin an Ising model, which require a Boltzmann machinewith a linear synaptic function derived from the gradientof the energy function. It is straightforward to check thatwith a = b = c = − A = 1 in Eq. (11), we obtain theresult stated earlier in Eq. (3). Two q-bit gates:
Any two q-bit gate is in generaldescribed by a transformation matrix of the form U ( d ) = | ( d − (cid:105) | ( d − (cid:105) | ( d − (cid:105) | ( d − (cid:105)(cid:104) ( d ) | a b a b (cid:104) ( d ) | c A c A (cid:104) ( d ) | a b a b (cid:104) ( d ) | c A c A (13) Using Eq. (7) we can immediately write the energy func-tion in tabular form E ( d ) = s ( d − s ( d − s ( d ) s ( d )
00 10 01 1100 − ln ( a ) − ln ( b ) − ln ( a ) − ln ( b ) − ln ( c ) − ln ( A ) − ln ( c ) − ln ( A ) − ln ( a ) − ln ( b ) − ln ( a ) − ln ( b ) − ln ( c ) − ln ( A ) − ln ( c ) − ln ( A ) (14) Once again this tabular result can be translated into aBoolean function like Eq. (12), but the energy functionwill have three-body and four-body interactions whosegradient leads to non-linear synaptic terms. Eq. 14 withsuch three-body and four-body terms represents a higher-order Ising model that has been discussed in the learningcontext by Ref. [32]. Such “generalized” Ising models canbe solved much like ordinary Ising models with two-bodyinteractions, provided that the synaptic feedback can becomputed, for example by an FPGA [33]. Naturally, theresistive crossbar arrays that accelerate linear synapticoperations [34] would not be suitable for this purpose.It is possible to eliminate the three-body and four-bodyterms at the expense of additional p-bits by decomposingthe gates in terms of a sequence of a rotation operations[14, 35], making the synaptic function linear. Alterna-tively, the three and four-body interactions can be re-duced to standard 2-body interactions through the useof auxiliary variables [36, 37].
Ising Machines:
In recent years, various incarna-tions of special-purpose hardware accelerators known as“Ising Machines” have emerged [33, 38–43] to simulatethe statistical mechanics of Ising models, onto whichmany known combinatorial optimization problems havebeen mapped [44]. These special purpose machines canobtain a very large number of samples per second [45]that is far more than what can be obtained using or-dinary computers. In principle, our methodology pre-sented in this paper should be directly implementable insuch dedicated Ising Machines to simulate larger quan-tum systems by computing the phase exp ( − i (cid:61) ( E )) for each path as a post-processing step to compute outputstates of interest. Sampling with complex contributions:
Algo-rithm 1 assumes that the mapping from quantum circuitto p-bit has already been carried out and available for theprocedure to use. In the current setting, the algorithmalso assumes that the entries of the input wavefunction, ψ in can only be “0” and “1”, though any other non-trivial input wavefunction can be encoded by insertingadditional gates in the quantum circuit.Certain gates like controlled-NOT (CNOT),controlled-controlled-NOT (CCNOT or Toffoli) orthe modular exponentiation part of Shor’s order findingcircuit are deterministic. In a true quantum computerthese would need to be implemented with quantumgates. But in our probabilistic emulation we canimplement them simply as ordinary logic gates. Sincelogic gates are directed we would not have a single BMany more; instead we would have separate BMs withdirected connections. Consequently, the input p-bitsfor a deterministic gate should be updated before theoutput p-bits although the p-bits within each BM canbe updated in random order. Algorithm 1
Estimating the output wavefunction ψ out from U ψ in using probabilistic sampling N T : number of time samples N m : number of p-bits n in : number of input p-bits n out : number of output p-bits procedure EstimatePsi ( ψ in ) (cid:46) Output: ψ out for j ← n out do ψ out,j ← (cid:46) Initialize count for ψ out for j ← n in do m j ← ψ in,j (cid:46) Initialize input p-bits for k ← N T do for j ← ( n in + 1) to N m do m j ← compute E from the mapping with { m } m j ← compute E from the mapping with { m } I j ← (cid:60) ( E − E ) draw r from uniform distribution [0 , if σ ( I j ) > r then (cid:46) σ ( x ) = (1 + exp ( − x )) − m j ← s ← binary2decimal( { m } ) ψ out,s ← ψ out,s + exp ( − i (cid:61) ( E )) else m j ← s ← binary2decimal( { m } ) ψ out,s ← ψ out,s + exp ( − i (cid:61) ( E )) Z ← (cid:80) nout j =1 | ψ out,j | for j ← → n out do ψ out,j ← ψ out,j /Z (cid:46) Normalize ψ out Stoquastic Hamiltonians:
In this paper we fo-cus on GQC which is based on unitary transformations U ∝ exp ( − iH ). In this case the energy functions arerarely real since that requires the elements of the U ma-trix to be all real and positive, which is seldom the case.However, our approach is also applicable to adiabaticquantum computing (AQC) based on exp ( − βH ) whichcan often have purely real and positive elements leadingto purely real energy functions. Such Hamiltonians areclassified as stoquastic (see for example, [46]) and can beemulated with standard BM’s. Non-stoquastic Hamilto-nians requiring complex BM’s form a special subset of allproblems of interest in AQC, while in GQC many prob-lems belong to this category.Interestingly, in one respect GQC is simpler than AQCbecause the GQC evolution operators are naturally builtout of the product of few (typically 1 and 2) q-bit oper-ations. U = exp ( − iH ) exp ( − iH ) · · · By contrast, it is not straightforward to do the re-verse, namely to break up the evolution operatorexp ( − βH ) for AQC into a product of separate termscorresponding to the components H = H + H + · · · , exp ( − βH ) (cid:54) = exp ( − βH ) exp ( − βH ) · · · unless thecomponents H , H , · · · commute. The standard ap-proach is the Suzuki-Trotter transformation [47] whichbreaks up H into r replicas by writingexp ( − βH ) = [exp ( − βH/r )] r ≈ (cid:0) exp ( − βH /r ) exp ( − βH /r ) · · · (cid:1) r assuming that r is large enough to make the commutatorsof H /r , H /r, · · · negligible. This paper focuses on GQCwhere replicas are not needed as explained above. III. TWO Q-BIT GROVER SEARCH
The three q-bit circuit in Fig. 4 represents a standardGrover search circuit [1] using ten Hadamard (H) gates,four Pauli-X gates, one CNOT gate and one CCNOTgate. Each of these gates adds one p-bit, so that theoverall p-bit network requires 19 p-bits, three of whichrepresent the input which is clamped to | (cid:105) . At plane1, all eight possible output states are equally probablebut half the wavefunctions are positive while the otherhalf are negative. All results are described well by theprobabilistic method with the correct signs. At plane4, a single peak rises, again correctly emulated by thep-network.Once again it is instructive to look at the total sign, S total (= (cid:80) S α ) at various planes: d = 1 → S total ≈ d = 2 → S total ≈ . d = 3 → S total ≈ . d = 4 → S total ≈ . d = 1(a) HHH HH Hd = 2 H Hd = 3 HHd = 4d = 4d = 3 d = 2d = 1(b)
Output configura�on S Output configura�on S Output configura�on S Output configura�on S FIG. 4 . Grover search circuit using three q-bits: (a)Circuit schematic (from [1]). Dashed box around the CC-NOT gate is the oracle that can identify | (cid:105) ), (b) Outputsat planes 1, 2, 3 and 4 showing the results obtained from theprobabilistic emulation with 19 p-bits. These outputs are gen-erated from 10 samples. See text for a discussion regardinghow the total sign ( (cid:80) S α ) progressively decreases as measuredin increasing depth, d. IV. SHOR’S ALGORITHM
The n = t + m q-bit circuit in Fig. 5 represents thegeneral outline of the Shor’s algorithm designed to findthe order, r of the function f ( x ) = a x mod N . If wechoose N = 143 and a = 43, then the algorithm requires t = 16 (satisfying 143 < t < × ) in the first reg-ister and m = 8 (satisfying m = (cid:100) log N (cid:101) ) in the second HHH HH H H H fi r s t r e g i s t e r t q - b i t ss e c o n d r e g i s t e r m q - b i t s R − R − R − R − H inverse quantum Fourier transform, IQFT t 1 U − t 2 U − U U modular exponen�a�on, t 3 U − x a mod N t 11 R −− t 11 R −− R − FIG. 5 . Shor’s order-finding circuit:
The circuit consists of two quantum registers: the first register is of size t q-bits andinitially set to | (cid:105) state, while the second register is of size m q-bits and initially set to | (cid:105) state. First, a bank of Hadamardgates are applied to the top register to create a superposition of all classically possible states, (cid:80) t − x =0 | x (cid:105) . Then a modularexponentiation operation, f ( x ) = a x mod N is performed on the bottom register for each | x (cid:105) in the superposed state. The U transformation in the modular exponentiation block performs U | u (cid:105) = | au mod N (cid:105) . After that a t q-bit inverse quantumFourier transform is performed which generates r distinct peaks (which is the order of f ( x )) in the probability distribution ofthe q-bits in the top register. register, totaling n = 24 q-bits. The modular exponenti-ation seems to be the detailed component to implementwith quantum computers and different approaches havebeen suggested to address its concrete implementation[48–52].The probabilistically mapped circuit corresponding tothe quantum circuit in Fig. 5 is shown in Fig. 6. Themodular exponentiation operation is computed determin-istically such that the m q-bit second register requiresonly m output p-bits that are loaded with the classicallycomputed value of a x mod N at each time step duringsampling. The mapping presented in this work leadsto all imaginary coupling and biases and therefore themapped circuit basically acts as a random number gen-erator with phase calculations performed within the post-processing of the samples. But as we have noted earlier,the advantages of probabilistic sampling can be exploitedby appropriately taming the circuit that introduces somereal parts (that vary with the state of the circuit) in thecoupling and biases.For the quantum circuit in Fig. 5 with n = 16 + 8 q-bits, Fig. 7 shows the probability distribution for each ofthe 2 = 65 ,
536 outputs obtained using 10 samples.The 6 peaks are clearly evident, giving r = 6. The factors are now obtained fromgcd( a r/ ± , N ) ∈ { , } In our probabilistic implementation of small Shor cir-cuits, we construct a full joint probability distribution ofp-bits that correspond to the q-bits of the first registerby taking enough samples and obtain the period “r” bycounting the total number of peaks, for simplicity. Wenote however that this approach becomes impractical forlarge periods and the standard method of deducing theperiod “r” from a single peak measurement by a contin-ued fraction algorithm needs to be adopted in such cases.Once again, it is interesting to stress the key differencewith standard Boltzmann machines where we engineerthe real energy landscape so that most paths endpreferentially on the peaks. By contrast, the pathsin this complex Boltzmann machine reach all possibleoutputs almost equally. Nevertheless only six peaksstand out because the paths reaching them all have thesame phase and add coherently, while for other outputsthe paths cancel out. The Boltzmann machine withcomplex energies allows us to emulate the other q-bitnetwork accurately, but many more samples are neededto get the correct distribution, since a large number iswasted canceling each other out at the nulls. x a mod N x t i n p u t p - b i t s m i n p u t p - b i t s t o u t p u t p - b i t s m o u t p u t p - b i t s FIG. 6 . Shor’s order-finding circuit mapped with p-bits:
The circuit consists of 3 t + 2 m p-bits; corresponding tothe t q-bits in the first register in the quantum circuit, thereare t input p-bits and all these p-bits are initially set to 0.Each Hadamard gate in the quantum circuit then adds a p-bit each but the rotation operations do not add any. m q-bitsin the second register correspond to additional m input p-bitsand these are initially set to 0 except the least significant onewhich is set to 1. x is derived from the states of the p-bitsinside the dashed box as shown and a x mod N is computedclassically and the result is set to another m p-bits. The biasinputs to each individual p-bits are not shown explicitly. Scaling properties:
To investigate the scaling prop-erties we use a simpler version of the circuit by first ap-plying an additional transformation, T such that T : a x mod N (cid:55)→ x mod r and if r (the order of a x mod N ) is allowed only to be apower of 2, then this allows us to perform the modularexponentiation with simple CNOT circuits. Althoughthis uses the explicit knowledge of r but is not uncommon[53]. The simplified Shor’s circuit is shown in Fig. 8. Notethat the U transformation is now defined differently thanthe original one, i.e., we use U | u (cid:105) = | u + 1 mod r (cid:105)
0 21845 43691 65535
Output configuration P r o b a b i l i t y FIG. 7 . Shor’s general order-finding circuit appliedto a = 43 , N = 143 : The output obtained from 10 samplesshows 6 peaks (corresponding to r = 6).) HHH n / q - b i t s n / q - b i t s H t 1 U − t 2 U − U Umodular exponentiationx mod r IQFT
FIG. 8 . Simplified Shor’s circuit: both registers are ofthe size n/ | (cid:105) state. The U transformation is now defined as U | u (cid:105) = | u + 1 mod r (cid:105) . Thelimitation on r is that it can only be a power of 2. This redef-inition helps to simplify the implementations of U j matriceswith the use of simple CNOT gates. The arrows on someCNOT gates indicate that depending on the actual order r ifless than maximum order ( r max ), would need to be replacedwith identity (no gate). so that U j | u (cid:105) = U | u + j mod r (cid:105) . and note that j can be written as j n − n − + j n − n − + · · · + j where { j n − j n − · · · j } denotes the binary rep-resentation of j . The use of this new definition of U re-quires the second register now to be set with | (cid:105) stateinitially instead of | (cid:105) state used in the general Shor’scircuit. Also note that now half the q-bits ( n/ l max )in the first register represent the maximum order, r max =2 l max of f ( x ) for all possible a values. The other half areauxiliary q-bits required for the modified modular expo-nentiation. In the examples presented we use just twoCNOTs in the modular exponentiation block leading tofour peaks in the output indicating that the order of f ( x )is 4.Fig. 9 shows the peaks on a logarithmic scale for asimulation with 36 q-bits emulated with 90 p-bits using N T = 10 samples. Note how the four peaks stand outfrom the rest making them clearly identifiable. l o g ( P r o b a b i l i t y ) Output q-bit configuration 1050 0.5 1.0 1.5 2.0 2.5-10-8-6-4-2 n = 36
FIG. 9 . Shor’s period-finding circuit with n = 36 : Alloutput peaks on a logarithmic scale obtained from a proba-bilistic simulation with N T = 10 samples. Note that out of2 = 262 ,
144 possible outputs, four peaks are clearly sepa-rated from the rest.
Fig. 10 shows the ratio (
P eak /P eak
4) plottedagainst the number of samples N T . At the end of the sim-ulation, the probabilities are sorted in a descending orderof their magnitude and the 4-th (we expect to see fourpeaks) and the 5-th highest probabilities (this peak isnon-desired, ideally should be zero) are labeled as P eak
P eak n/ suggesting that N T ∝ n/ .This is of course inferior to the polynomial scaling in n expected from a q-bit implementation that has madeShor’s algorithm legendary. But we observe that the scal-ing depends on the number of q-bits and not the numberof p-bits which in this example equals 5 n/ × non-sparsetransformation matrices with over 10 elements. Bycontrast, the scaling curves shown in Fig. 9 required nosignificant memory and were carried out on an ordinary (a)(b) -5 0 5 10 log (N T /2 (n/2) ) -8-6-4-2 0 l o g ( P e a k / P e a k ) n = 20n = 24n = 28n = 32n = 36
12 14 16 18 20 22 log (N T ) l o g ( P e a k / P e a k ) n = 20n = 24n = 28n = 32n = 36 FIG. 10 . Computational scaling for a probabilisticsimulation of Shor’s period-finding circuit: (a) Resultsshown for n = 20, 24, 28, 32 and 36 suggest that (b) thenumber of samples needed for a specified peak suppressionratio scales as N T ∝ n/ . See text for discussion. laptop.Finally we note that this is a straightforward imple-mentation leading to a Boltzmann machine with a trivialreal part that distributes paths equally among all possi-ble outputs, relying on sign cancellation to make peaksemerge. In this mode, the probabilistic sampling methodis no different from random guessing, representing anexample of an extremely severe sign problem. It maybe possible to tame this sign problem with proper basistransformations so that the real part plays a more sig-nificant role, thereby making use of the powerful proba-bilistic sampling techniques. V. AUGMENTED PROBABILISTICALGORITHMS
Finally we present an example of a standard classicaloptimization algorithm implemented using a real energyfunction within the same framework, making it straight-forward to augment it with an imaginary component in-troducing quantum interference. Consider for example0 (a)
Exact Boltzmann probabilitiesK = 0/N K = 10/N K = 50/N K = 100/N(b)
Sampled probabilities (without phase):(c)
Sampled probabilities (with phase):
0 13 23 300 132330 Y
0 13 23 300 132330 -3 Y
0 13 23 30 X -3 Y
0 13 23 300 132330 Y
0 13 23 300 132330 -2 Y
0 13 23 30 X -2 Y
0 13 23 300 132330 Y
0 13 23 300 132330 Y
0 13 23 30 X -1 Y
0 13 23 300 132330 Y
0 13 23 300 132330 Y
0 13 23 30 X -3 Y FIG. 11 . Non-stoquasticity to augment classical optimization algorithms: (a) Heatmap of exp ( − ε ) (see Eq. 15)calculated exactly as a function of ( X, Y ). (b) Heatmap of probabilities calculated from 10 samples without including thephase. (c) Same as (b), but including the phase exp ( − i (cid:61) ( E )) (see Eq. 16). For K = 0 /N (random guessing) and K = 10 /N ,both minima (13 ,
23) and (23 ,
13) are located correctly. For K = 50 /N only one minimum (13 ,
23) is located. For K = 100 /N the solution is stuck in the wrong minima (fourth columns of second and third row). the problem of factorizing a number N into X and Y bydefining a cost function C = (cid:12)(cid:12) XY − N (cid:12)(cid:12) and using a BMto minimize it by defining an energy ε proportional to C by a constant K : ε = K × C = K (cid:12)(cid:12) XY − N (cid:12)(cid:12) (15)This is not a particularly attractive optimization algo-rithm [54], we choose it only for illustrative purposes.The point is that we can implement it using the sameapproach that we have described in this paper, the onlydifference being that the energy is not obtained from U-matrix elements defined by quantum gates using Eq. (7).Instead it is given simply by Eq. (15).The difficulty with this algorithm can be appreciatedby looking at the top two rows of Fig. 11. The top rowshows the statistical factor exp ( − ε ) for increasing valuesof K . For the smallest K = 1 /N the minima are shal-low, progressively breaking up into deep valleys as K isincreased. Clearly we would like a large value of K so thecorrect answers are located in deep valleys well separatedfrom the rest.But the top row represents the exact result which iseasy to calculate for this toy problem with N = 13 × K , sampling easily approaches the exact distribution, but for large K it gets stuck inlocal minima far from the two global minima at (13 , , (cid:61) ( E ) into the energythat selects the exact solutions through interference. Thelast row in Fig. 11 shows the results obtained by samplingusing the complex energy function E = ε + i πrC (16)Here r is a random number between 0 and 1 and can beimplemented in hardware by including auxiliary p-bitswhose collective state determines r so that it fluctuatesrandomly.The imaginary part of the energy 2 πrC makes the con-tributions of different paths fluctuate wildly for all out-puts except for the correct solution with C = 0. As a re-sult when we sample with the phases taken into account,the correct solutions (13 ,
23) and (23 ,
13) are singled outas shown in the last row of Fig. 11. With K = 50 /N the samples get stuck on a fraction of the valleys, andthe phase helps select out only one minimum (13,23).With an even larger K = 100 /N the samples are stuckin wrong local minima, and the phase has no chance offinding the right answers.1This approach does not really enhance the classicalprobabilistic algorithm, since one could identify the cor-rect answer simply by asking if C = 0 rather than invokeinterference through a complex energy function. Even ifthe imaginary phase is added to the paths while samplingis performed, these phases do not influence the samplingprobabilities, as such, there would be no computationalbenefit for classical samplers, as they would suffer fromthe sign problem. Quantum annealers however could en-gineer the implementation of such random phases thatwould influence sampling dynamics and the algorithmcould lead to a speed-up in identifying the correct an-swer like the Grover’s algorithm [3], but that remains tobe seen. Our purpose here is to show that the frameworkpresented in this paper is general enough to encompassquantum gates as well as classical optimization, thus fa-cilitating a cross-fertilization of ideas. V. CONCLUSIONS
In summary, we have presented what could be viewedas an extension of sign-corrected QMC to GQC algo-rithms, leading to a systematic procedure for translat-ing any given quantum circuit into a network of classi-cal probabilistic p-bits described by a complex energyfunction E . The real part of E governs the time evolu-tion just like standard Boltzmann machines. But unlikeBoltzmann machines, each sample makes a complex con-tribution of unit magnitude given by exp ( − i (cid:61) ( E )). Theprobabilistic method provides significant memory advan-tages allowing us to perform calculations on an ordinarylaptop, that would have required the multiplication ofmatrices with ∼ elements in a deterministic ap-proach. Even larger problems should be accessible ondedicated Ising computers.However, we note that the probabilistic emulation suf-fers from the sign problem. A straightforward transla-tion of quantum algorithms like Grover search and Shor’speriod-finding leads to energy functions whose real partsare essentially constant and provide no discriminationamongst different outputs. Nulls in the output do notarise from a suppression of paths, but from a cancellationof paths through the phase exp ( − i (cid:61) ( E )). Consequentlymany samples are wasted on low probability outputs, un-like a true quantum computer. We introduce a quanti-tative metric for this inefficiency of a probabilistic em-ulation relative to a quantum computer in terms of thetotal sign which measures the severity of the well-known“sign problem” in QMC algorithms. The inverse of S Total indicates the factor by which the computation time isincreased. For example, a straightforward probabilisticimplementation of Shors algorithm for n qbits leads toan S Total ∼ exp ( − n/ n/ −(cid:60) ( E )) with path cancellationthrough the phase factor exp ( − i (cid:61) ( E )). This augmentedoptimization algorithm could potentially be implementedin hardware quantum annealers that can attach the imag-inary phase factor to different samples through a collec-tion of hidden variables [55]. VI. APPENDIX:ZERO ELEMENTS IN U-MATRIX
The procedure laid out in Sec. II for obtaining energyfunctions for given transformation matrices is straight-forward. One point to note, however, is that often thereare zero elements in the matrix which will lead to singularvalues for the logarithmic functions in Eqs. (12) or (14).A general but approximate way to deal with zero ele-ments is to replace them with e − J , J being a suitablylarge positive number. However, in special cases an ex-act approach is possible, which is best illustrated withexamples. One q-bit examples:
Consider the transformation R z ( γ ) representing a one q-bit rotation about the z -axis: R z ( γ ) = e − iγσ z / = (cid:20) | ( d − (cid:105) | ( d − (cid:105)(cid:104) ( d ) | e − iγ/ (cid:104) ( d ) | e + iγ/ (cid:21) (17)Replacing the zero off-diagonal elements with e − J , J be-ing a large positive number, we have E ( d ) = − i γ (cid:0) − s ( d − − s ( d ) (cid:1) + lim J →∞ J (cid:0) s ( d − + s ( d ) − s ( d − s ( d ) (cid:1) (18)This is the straightforward approximate approach de-scribed above. Alternatively, we could eliminate s ( d ) asan independent variable by setting it equal to s ( d − , sothat s ( d ) = s ( d − E ( d ) = − i γ (cid:0) − s ( d − (cid:1) (19)This approach can also be used for zeros on the diag-onal. Consider the Pauli-X gate: X = (cid:20) (cid:21) (20)We can eliminate s ( d ) as an independent variable by set-ting it equal to 1 − s ( d − , so that s ( d ) = 1 − s ( d − E ( d ) = 0 (21)Note that this is essentially a deterministic NOT gate.2 Two q-bit exchange operator:
Consider the twoq-bit exchange operator acting between q-bit 1 and 2defined as follows: J = | ( d − (cid:105) | ( d − (cid:105) | ( d − (cid:105) | ( d − (cid:105)(cid:104) ( d ) | e iγ/ (cid:104) ( d ) | e − iγ/ (cid:104) ( d ) | e − iγ/ (cid:104) ( d ) | e iγ/ (22)Since the matrix is diagonal, we can eliminate s ( d )1 , s ( d )2 by writing s ( d )1 = s ( d − , s ( d )2 = s ( d − E ( d ) = iγ (cid:0) − s ( d − (cid:1)(cid:0) − s ( d − (cid:1) (23)For two q-bit operations with off-diagonal elements, thegeneral approach would be to use exp ( − J ) to replacezero elements if present. Instead we may be able toreduce the number of independent variables for specificcases as illustrated below. CNOT Gate:
Consider the CNOT gate described bythe transformation matrix: U CNOT = | ( d − (cid:105) | ( d − (cid:105) | ( d − (cid:105) | ( d − (cid:105)(cid:104) ( d ) | (cid:104) ( d ) | (cid:104) ( d ) | (cid:104) ( d ) | In this case, we can write s ( d )1 = s ( d − + s ( d − − s ( d − s ( d − s ( d )2 = s ( d − E ( d ) = 0 (24)It can be seen that this is essentially a deterministic XORgate: s ( d )1 = XOR (cid:0) s ( d − , s ( d − (cid:1) (25) CCNOT Gate:
Similarly for the CCNOT gate, wecan write s ( d )1 = s ( d − + s ( d − s ( d − − s ( d − s ( d − s ( d − s ( d )2 = s ( d − , s ( d )3 = s ( d − E ( d ) = 0 (26)This too is essentially a deterministic gate: s ( d )1 = XOR (cid:18) s ( d − , (cid:0) s ( d − AND s ( d − (cid:1)(cid:19) (27) ACKNOWLEDGMENTS
This work was supported in part by ASCENT, oneof six centers in JUMP, a Semiconductor Research Cor-poration (SRC) program sponsored by DARPA. KYCgratefully acknowledges support from Center for Scienceof Information (CSoI), an NSF Science and TechnologyCenter, under grant CCF-0939370. The authors thankBrian Sutton and Marc Cahay for extensive commentsand suggestions on an earlier version of this manuscript.The authors are also grateful to Diptiman Sen for usefuldiscussions related to non-stoquastic Hamiltonians. [1] Michael A. Nielsen and Isaac L. Chuang,
Quantum Com-putation and Quantum Information: 10th AnniversaryEdition , 10th ed. (Cambridge University Press, USA,2011).[2] Peter W. Shor, “Polynomial-time algorithms for primefactorization and discrete logarithms on a quantum com-puter,” SIAM Review , 303–332 (1999).[3] Lov K. Grover, “A fast quantum mechanical algorithmfor database search,” in Annual ACM Symposium onTheory of Computing (ACM, 1996) pp. 212–219.[4] Sergio Boixo, Sergei V. Isakov, Vadim N. Smelyanskiy,Ryan Babbush, Nan Ding, Zhang Jiang, Michael J.Bremner, John M. Martinis, and Hartmut Neven, “Char-acterizing quantum supremacy in near-term devices,”Nature Physics , 595–600 (2018).[5] Sergey Bravyi, David Gosset, and Robert K¨onig, “Quan-tum advantage with shallow circuits,” Science , 308–311 (2018).[6] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, Rupak Biswas, SergioBoixo, Fernando G. S. L. Brandao, David A. Buell, BrianBurkett, Yu Chen, Zijun Chen, Ben Chiaro, RobertoCollins, William Courtney, Andrew Dunsworth, Ed-ward Farhi, Brooks Foxen, Austin Fowler, Craig Gidney,Marissa Giustina, Rob Graff, Keith Guerin, Steve Habeg-ger, Matthew P. Harrigan, Michael J. Hartmann, AlanHo, Markus Hoffmann, Trent Huang, Travis S. Humble,Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Dvir Kafri,Kostyantyn Kechedzhi, Julian Kelly, Paul V. Klimov,Sergey Knysh, Alexander Korotkov, Fedor Kostritsa,David Landhuis, Mike Lindmark, Erik Lucero, DmitryLyakh, Salvatore Mandr`a, Jarrod R. McClean, MatthewMcEwen, Anthony Megrant, Xiao Mi, Kristel Michielsen,Masoud Mohseni, Josh Mutus, Ofer Naaman, MatthewNeeley, Charles Neill, Murphy Yuezhen Niu, Eric Os-tby, Andre Petukhov, John C. Platt, Chris Quintana,Eleanor G. Rieffel, Pedram Roushan, Nicholas C. Ru-bin, Daniel Sank, Kevin J. Satzinger, Vadim Smelyan- skiy, Kevin J. Sung, Matthew D. Trevithick, AmitVainsencher, Benjamin Villalonga, Theodore White,Z. Jamie Yao, Ping Yeh, Adam Zalcman, Hartmut Neven,and John M. Martinis, “Quantum supremacy using aprogrammable superconducting processor,” Nature ,505–510 (2019).[7] Daniel Gottesman, “The heisenberg representation ofquantum computers,” (1998).[8] Scott Aaronson and Daniel Gottesman, “Improved sim-ulation of stabilizer circuits,” Phys. Rev. A , 052328(2004).[9] Adam Smith, M. S. Kim, Frank Pollmann, and JohannesKnolle, “Simulating quantum many-body dynamics on acurrent digital quantum computer,” npj Quantum Infor-mation , 106 (2019).[10] Yiqing Zhou, E. Miles Stoudenmire, and Xavier Wain-tal, “What limits the simulation of quantum computers?”(2020).[11] Juan Carrasquilla, Di Luo, Felipe P´erez, Ashley Milsted,Bryan K. Clark, Maksims Volkovs, and Leandro Aolita,“Probabilistic simulation of quantum circuits with thetransformer,” (2019).[12] John Napp, Rolando L. La Placa, Alexander M. Dalzell,Fernando G. S. L. Brandao, and Aram W. Harrow, “Effi-cient classical simulation of random shallow 2d quantumcircuits,” (2019).[13] Kyungjoo Noh, Liang Jiang, and Bill Fefferman, “Effi-cient classical simulation of noisy random quantum cir-cuits in one dimension,” arXiv preprint arXiv:2003.13163(2020).[14] M. Van den Nest, W. D¨ur, R. Raussendorf, and H. J.Briegel, “Quantum algorithms for spin models and simu-lable gate sets for quantum computation,” Phys. Rev. A , 052334 (2009).[15] Joseph Geraci and Daniel A Lidar, “Classical ising modeltest for quantum circuits,” New Journal of Physics ,075026 (2010).[16] G De las Cuevas, W D¨ur, M Van den Nest, and M AMartin-Delgado, “Quantum algorithms for classical lat-tice models,” New Journal of Physics , 093021 (2011).[17] Michael J. Bremner, Ashley Montanaro, and Dan J.Shepherd, “Average-case complexity versus approximatesimulation of commuting quantum computations,” Phys.Rev. Lett. , 080501 (2016).[18] Keisuke Fujii and Tomoyuki Morimae, “Commutingquantum circuits and complexity of ising partition func-tions,” New Journal of Physics , 033003 (2017).[19] Bjarni Jnsson, Bela Bauer, and Giuseppe Carleo,“Neural-network states for the classical simulation ofquantum computing,” (2018), arXiv:1808.05232 [quant-ph].[20] Christian Pehle and Christof Wetterich, “Neuromorphicquantum computing,” (2020), arXiv:2005.01533 [cond-mat.dis-nn].[21] M. Suzuki, Quantum Monte Carlo methods in equilibriumand nonequilibrium systems: proceedings of the NinthTaniguchi International Symposium, Susono, Japan,November 14-18, 1986 , Springer series in solid-state sci-ences (Springer-Verlag, 1987).[22] Giuseppe Carleo and Matthias Troyer, “Solving thequantum many-body problem with artificial neural net-works,” Science , 602–606 (2017).[23] Giuseppe Carleo, Yusuke Nomura, and MasatoshiImada, “Constructing exact representations of quantum many-body systems with deep neural networks,” NatureCommunications , 5322 (2018).[24] Xun Gao and Lu-Ming Duan, “Efficient representation ofquantum many-body states with deep neural networks,”Nature Communications , 662 (2017).[25] Juan Carrasquilla, Giacomo Torlai, Roger G. Melko, andLeandro Aolita, “Reconstructing quantum states withgenerative models,” Nature Machine Intelligence , 155–161 (2019).[26] The Feynmann Lectures on Physics, Volume III .[27] Matthias Troyer and Uwe-Jens Wiese, “Computationalcomplexity and fundamental limitations to fermionicquantum monte carlo simulations,” Phys. Rev. Lett. ,170201 (2005).[28] Milad Marvian, Daniel A. Lidar, and Itay Hen, “Onthe computational complexity of curing non-stoquastichamiltonians,” Nature Communications , 1571 (2019).[29] Dominik Hangleiter, Ingo Roth, Daniel Nagaj, and JensEisert, “Easing the monte carlo sign problem,” (2019).[30] Erez Berg, Max A. Metlitski, and Subir Sachdev, “Sign-problem–free quantum monte carlo of the onset of antifer-romagnetism in metals,” Science , 1606–1609 (2012).[31] John F Wakerly, Digital Design Principles and Practices (Prentice Hall, Englewood Cliffs, NJ, 2016).[32] T J Sejnowski, “Higher-order boltzmann machines,” in
AIP Conference Proceedings 151 on Neural Networks forComputing (American Institute of Physics Inc., USA,1987) p. 398403.[33] Peter L McMahon, Alireza Marandi, Yoshitaka Harib-ara, Ryan Hamerly, Carsten Langrock, Shuhei Tamate,Takahiro Inagaki, Hiroki Takesue, Shoko Utsunomiya,Kazuyuki Aihara, et al. , “A fully programmable 100-spincoherent ising machine with all-to-all connections,” Sci-ence , 614–617 (2016).[34] M. Hu, J. P. Strachan, Z. Li, E. M. Grafals, N. Davila,C. Graves, S. Lam, N. Ge, J. J. Yang, and R. S. Williams,“Dot-product engine for neuromorphic computing: Pro-gramming 1t1m crossbar to accelerate matrix-vector mul-tiplication,” in (2016) pp. 1–6.[35] Jaehyun Kim, Jae-Seung Lee, and Soonchil Lee, “Im-plementing unitary operators in quantum computation,”Phys. Rev. A , 032312 (2000).[36] Shuxian Jiang, Keith A Britt, Alexander J McCaskey,Travis S Humble, and Sabre Kais, “Quantum annealingfor prime factorization,” Scientific reports , 1–9 (2018).[37] JD Biamonte, “Nonperturbative k-body to two-bodycommuting conversion hamiltonians and embeddingproblem instances into ising spins,” Physical Review A , 052331 (2008).[38] Masanao Yamaoka, Chihiro Yoshimura, Masato Hayashi,Takuya Okuyama, Hidetaka Aoki, and Hiroyuki Mizuno,“A 20k-spin ising chip to solve combinatorial optimiza-tion problems with cmos annealing,” IEEE Journal ofSolid-State Circuits , 303–309 (2015).[39] Takahiro Inagaki, Yoshitaka Haribara, Koji Igarashi,Tomohiro Sonobe, Shuhei Tamate, Toshimori Honjo,Alireza Marandi, Peter L McMahon, Takeshi Umeki,Koji Enbutsu, et al. , “A coherent ising machine for2000-node optimization problems,” Science , 603–606(2016).[40] Tianshi Wang and Jaijeet Roychowdhury, “Oscillator-based ising machine,” arXiv preprint arXiv:1709.08102(2017). [41] Jeffrey Chou, Suraj Bramhavar, Siddhartha Ghosh,and William Herzog, “Analog coupled oscillator basedweighted ising machine,” Scientific reports , 1–10(2019).[42] William A Borders, Ahmed Z Pervaiz, Shunsuke Fukami,Kerem Y Camsari, Hideo Ohno, and Supriyo Datta, “In-teger factorization using stochastic magnetic tunnel junc-tions,” Nature , 390–393 (2019).[43] S Dutta, A Khanna, J Gomez, K Ni, Z Toroczkai, andS Datta, “Experimental demonstration of phase transi-tion nano-oscillator based ising machine,” in (IEEE,2019) pp. 37–8.[44] Andrew Lucas, “Ising formulations of many np prob-lems,” Frontiers in Physics , 5 (2014).[45] Brian Sutton, Rafatul Faria, Lakshmi A Ghantasala,Kerem Y Camsari, and Supriyo Datta, “Autonomousprobabilistic coprocessing with petaflips per second,”arXiv preprint arXiv:1907.09664 (2019).[46] Sergey Bravyi, David P Divincenzo, Roberto I Oliveira,and Barbara M Terhal, “The complexity of stoquas-tic local hamiltonian problems,” arXiv preprint quant-ph/0606140 (2006).[47] Masuo Suzuki, “Relationship between d-dimensionalquantal spin systems and (d+1)-dimensional ising sys-tems equivalence, critical exponents and systematic ap-proximants of the partition function and spin correla-tions,” Progress of Theoretical Physics , 1454–1469(1976).[48] Aidan Dang, Charles D. Hill, and Lloyd C. L. Hollen-berg, “Optimising Matrix Product State Simulations ofShor’s Algorithm,” Quantum , 116 (2019).[49] Archana Tankasala and Hesameddin Ilatikhameneh,“Quantum-kit: Simulating shor’s factorization of 24-bit number on desktop,” (2019), arXiv:1908.07187 [quant-ph].[50] Ben Lanyon, Till Weinhold, Nathan Langford, M Barbi-eri, Daniel James, Alexei Gilchrist, and Andrew White,“Experimental demonstration of a compiled version ofshor’s algorithm with quantum entanglement,” Physicalreview letters , 250505 (2008).[51] Thomas Monz, Daniel Nigg, Esteban A. Martinez,Matthias F. Brandl, Philipp Schindler, Richard Rines,Shannon X. Wang, Isaac L. Chuang, and Rainer Blatt,“Realization of a scalable shor algorithm,” Science ,1068–1070 (2016).[52] Lieven M. K. Vandersypen, Matthias Steffen, GregoryBreyta, Costantino S. Yannoni, Mark H. Sherwood, andIsaac L. Chuang, “Experimental realization of shor’squantum factoring algorithm using nuclear magnetic res-onance,” Nature , 883–887 (2001).[53] Michael R. Geller and Zhongyuan Zhou, “Factoring 51and 85 with 8 qubits,” Scientific Reports , 3023 (2013).[54] P. Henelius and S. Girvin, “A statistical mechanics ap-proach to the factorization problem,” (2011).[55] I. Ozfidan, C. Deng, A.Y. Smirnov, T. Lanting, R. Har-ris, L. Swenson, J. Whittaker, F. Altomare, M. Bab-cock, C. Baron, A.J. Berkley, K. Boothby, H. Chris-tiani, P. Bunyk, C. Enderud, B. Evert, M. Hager, A. Ha-jda, J. Hilton, S. Huang, E. Hoskinson, M.W. Johnson,K. Jooya, E. Ladizinsky, N. Ladizinsky, R. Li, A. Mac-Donald, D. Marsden, G. Marsden, T. Medina, R. Molavi,R. Neufeld, M. Nissen, M. Norouzpour, T. Oh, I. Pavlov,I. Perminov, G. Poulin-Lamarre, M. Reis, T. Prescott,C. Rich, Y. Sato, G. Sterling, N. Tsai, M. Volkmann,W. Wilkinson, J. Yao, and M.H. Amin, “Demonstrationof a nonstoquastic hamiltonian in coupled superconduct-ing flux qubits,” Phys. Rev. Applied13