Quantum Adiabatic Algorithms, Small Gaps, and Different Paths

MIT-CTP 4076 CERN-PH-TH-2009/175 Quan tum A diabatic Algorithms, Small Gaps, and Differen t P aths Edw ard F arhi, 1 Jeffrey Goldstone, 1 Da vid Gosset, 1 Sam Gutmann, 2 Harv ey B. Meyer, 1, 3 and P eter Shor 1, 4 1 Center for The or etic al Physics, Massachusetts Institute of T e chnolo gy, Cambridge, MA 02139 2 Dep artment of Mathematics, Northe astern University, Boston, MA 02115 3 Physics Dep artment, CERN, 1211 Geneva 23, Switzerland 4 Dep artment of Mathematics, Massachusetts Institute of T e chnolo gy, Cambridge, MA 02139 Abstract W e construct a set of instances of 3SA T which are not solv ed efficiently using the simplest quan tum adi- abatic algorithm. These instances are obtained by picking random clauses all consisten t with tw o disparate plan ted solutions and then penalizing one of them with a single additional clause. W e argue that by ran- domly mo difying the b eginning Hamiltonian, one obtains (with substantial probability) an adiabatic path that remov es this difficulty . This suggests that the quantum adiabatic algorithm should in general b e run on each instance with many different random paths leading to the problem Hamiltonian. W e do not know whether this tric k will help for a random instance of 3SA T (as opp osed to an instance from the particular set w e consider), esp ecially if the instance has an exp onen tial num b er of disparate assignments that violate few clauses. W e use a con tinuous imaginary time Quan tum Monte Carlo algorithm in a no vel wa y to numerically in vestigate the ground state as well as the first excited state of our system. Our arguments are supplemen ted b y Quantum Monte Carlo data from simulations with up to 150 spins. 1 I. INTR ODUCTION Quan tum adiabatic algorithms are designed for classical combinatorial optimization problems [8]. In the simplest case, suc h algorithms work by adiabatically evolving in the ground state of a system with Hamiltonian H ( s ) = (1 − s ) H B + sH P that is a function of a parameter s which is increased from 0 to 1 as a function of time. H B is called the b eginning Hamiltonian and H P , whic h is instance dep endent, is called the problem Hamiltonian. The minim um (for s ∈ [0 , 1] ) eigen v alue gap b etw een the ground state and first excited state of H ( s ) is related to the runtime of the adiabatic algorithm. If H ( s ) has an exponentially small minim um gap then the corresp onding algorithm is inefficient, whereas a minimum gap which scales in verse p olynomially corresp onds to an efficien t quantum adiabatic algorithm. Whether or not quan tum adiabatic algorithms can b e used to solv e classically difficult optimiza- tion problems efficiently remains to b e seen. Some numerical studies hav e examined the quantum adiabatic algorithm on random sets of instances of optimization problems where these sets are though t to b e difficult for classical algorithms. These studies ha v e rep orted p olynomial scaling of the minimum gap out to ab out 100 bits [6, 12, 20]. Whether this scaling p ersists at high bit num b er has recen tly b een called in to question [21]. Mean while there has b een no rigorous analytical result that characterizes the p erformance of the quan tum adiabatic algorithms on random instances of NP-complete problems. Ov er the years there ha ve b een a num b er of prop osed examples whic h w ere mean t to demonstrate failures of the adiabatic algorithm on sp ecific problems. In reference [18], v an Dam et al constructed examples intended to foil the adiabatic algorithm, but these examples used a nonlo cal cost function. Related 3SA T examples of v an Dam and V azirani [19] indeed cannot b e solv ed efficiently using the quan tum adiabatic algorithm. How ev er it w as sho wn in [5] that suc h 3SA T instances do not p ose a problem for the quantum adiabatic algorithm if, ha ving fixed a sp ecific problem Hamiltonian, one randomly c ho oses multiple interpolating paths b et w een the initial and final Hamiltonians and runs the adiabatic algorithm once for eac h random path. Fisher [9] has constructed an interesting but sp ecialized example on which the quantum adiabatic algorithm is inefficient. (A later example of Reic hardt [17] is based on this.) W e do not kno w if random path change helps here. (It is in teresting to note that in this case the run time scales like c √ n where n is the n umber of bits.) The authors of reference [22] p oin ted out that a certain adiabatic algorithm for 3SA T is not efficient. How ever, this w as due to a p ervers e and av oidable nonlo cal choice of b eginning Hamiltonian H B [7]. One purp ose of this pap er is to discuss a differen t t yp e of challenge to adiabatic optimization 2 Figure 1: Energy levels of the Hamiltonian H 0 ( s ) b efore adding the last term h to the problem Hamiltonian. The lo wer curve coincides with the upp er curve only at s = 1 . whic h has recently come to light. (See references [2 – 4] and in a different context [15].) It seems to us that in the history of challenges to adiabatic optimization, this may b e the most serious. The tak eaw ay message of this work is that a small minimum gap for H ( s ) can arise when the Hamiltonian H P has features which are seen in the following construction. First supp ose that w e start with a problem Hamiltonian H 0 P whic h has tw o degenerate ground states | z 1 i and | z 2 i corresp onding to bit strings of length n that differ in order n bits. W e then exp ect that the t wo lo west eigenv alues of H 0 ( s ) = (1 − s ) H B + sH 0 P will lo ok something lik e figure 1. Note that one curve is alw ays b elow the other except at s = 1 even though in the figure they app ear to meld b ecause they ha ve the same slop e at s = 1 . Supp ose that the upp er curv e (the first excited state for s near 1 ) approaches the state | z 2 i and the low er curve approaches the state | z 1 i as s → 1 . W e now form the problem Hamiltonian H P = H 0 P + h , where h is a term that penalizes the state | z 1 i but not the state | z 2 i . Then, for the tw o lo west eigenv alues of H ( s ) = (1 − s ) H B + sH P , we exp ect to hav e the situation pictured in figure 2, where there is a small gap near s = 1 . In this pap er we discuss a simple wa y of generating random instances of 3SA T where the cor- resp onding adiabatic Hamiltonian H ( s ) has the difficulty discussed ab ov e. W e argue in section I I I that this problem can b e o v ercome (for the set of instances we consider) by randomizing the choice of b eginning Hamiltonian. W e present a contin uous imaginary time Quantum Monte Carlo algorithm, which is a mo difica- 3 Figure 2: Energy levels of the Hamiltonian H ( s ) . There is a tiny gap at s ? . tion of the Heat Bath algorithm of Krzak ala et al [14]. (“Quantum Mon te Carlo” is a completely classical n umerical technique for finding prop erties of quantum systems and is not a quantum al- gorithm.) W e use this numerical technique in a no v el manner which allows us to in v estigate both the ground state and the first excited state of our Hamiltonian and thereby detect the presence or absence of a small gap b et ween them. This differs from the standard application of the Quan tum Mon te Carlo metho d in that w e are able to obtain information ab out the first excited state using a simple procedure which we hav e v alidated at low bit num b er by comparing our results to exact n umerical diagonalization. These Quantum Monte Carlo simulations supp ort our claim that the problem w e describ e can b e o vercome using random path change. I I. PR OBLEMA TIC INST ANCES W e no w describ e the metho d we use to generate n bit random instances of 3SA T that lead to quan tum adiabatic Hamiltonians with small minimum gaps. T o do this w e first generate an instance with exactly tw o satisfying assignments giv en by the bit strings 111 ... 1 and 000 ... 0 . Each clause c of the 3SA T instance specifies a subset of 3 bits i 1 ( c ) , i 2 ( c ) , i 3 ( c ) ∈ { 1 , .., n } and a particular assignmen t w 1 ( c ) w 2 ( c ) w 3 ( c ) to those three bits whic h is disallo wed. In order to only generate instances which are consistent with the bit strings 111 ... 1 and 000 ... 0 , we only use clauses for which w 1 ( c ) w 2 ( c ) w 3 ( c ) ∈ { 100 , 010 , 001 , 110 , 101 , 011 } . 4 W e add suc h clauses one at a time uniformly at random and stop as so on as these t w o bit strings are the only bit strings which satisfy all of the clauses that hav e b een added. (In practice to c heck whether or not this is the case we use a classical 3SA T solver.) W e write m for the total n umber of clauses in the instance. W e note that the num b er of clauses m obtained using this pro cedure scales lik e n log n . W e need this man y clauses in order to ensure that eac h bit is in volv ed in some clause. On the other hand, when the n umber of clauses is 5 n log n the probability of additional satisfying assignmen ts go es to zero as n → ∞ . W e now consider the problem Hamiltonian H 0 P corresp onding to this instance, which w e define to b e H 0 P = m X c =1 1 + ( − 1) w 1 ( c ) σ i 1 ( c ) z 2 ! 1 + ( − 1) w 2 ( c ) σ i 2 ( c ) z 2 ! 1 + ( − 1) w 3 ( c ) σ i 3 ( c ) z 2 ! . Eac h term in this sum is 1 if i 1 ( c ) i 2 ( c ) i 3 ( c ) violates the clause and 0 otherwise. F or the b eginning Hamiltonian w e choose H B = n X i =1  1 − σ i x 2  . (1) The t w o low est eigenv alues of the Hamiltonian H 0 ( s ) = (1 − s ) H B + sH 0 P will then b oth approach 0 as s → 1 , as in figure 1. The ground state for v alues of s whic h are sufficiently close to 1 will approac h either | 000 ... 0 i or | 111 ... 1 i as s → 1 . Supp ose for v alues of s close to 1 the state that approac hes | 000 ... 0 i has low est energy . Then w e add an extra term h 0 whic h acts on bits 1 , 2 and 3 and whic h p enalizes this state but not | 111 ... 1 i h 0 = 1 2  1 + σ 1 z 2   1 + σ 2 z 2   1 + σ 3 z 2  and pic k H P = H 0 P + h 0 . (Note that this extra term has a m ultiplicativ e factor of 1 2 , to a void any degeneracy of the first excited state at s = 1 .) In the case where the low est energy state near s = 1 approaches | 111 ... 1 i w e instead add h 1 = 1 2  1 − σ 1 z 2   1 − σ 2 z 2   1 − σ 3 z 2  . The Hamiltonian H ( s ) = (1 − s ) H B + sH P is then exp ected to hav e a small gap near s = 1 as depicted in figure 2. W e exp ect the lo cation s ? of the minim um gap to approach s = 1 as n → ∞ . 5 Lo cation of the Minimum Gap In order to determine the dep endence of the lo cation s ? of the av oided crossing on the num b er of spins n and the num b er of clauses m , we consider the p erturbative corrections to the energies of the states | 000 ... 0 i and | 111 ... 1 i around s = 1 . W e will show that the low order terms in the p erturbation series reliably predict a crossing at 1 s ? = 1 − Θ( 1 n 1 / 4  m n  3 4 ) . With the Hamiltonian constructed in the previous section, one of these states has energy 0 at s = 1 (call this the low er state | z L i ) and the other state has energy 1 2 at s = 1 (w e call this the upp er state | z U i ). W e write our Hamiltonian as H ( s ) = (1 − s ) n 2 + s " −  1 − s s  n X i =1 σ i x 2 + H P # (2) and we then consider the term −  1 − s s  P n i =1 σ i x 2 as a p erturbation to H P expanding around s = 1 . T o understand why we trust p erturbation theory to predict the lo cation of the near crossing, consider a system whic h is comp osed of t w o disconnected sectors, A and B , so the corresp onding Hamiltonian is of the form   H A ( s ) 0 0 H B ( s )   . In this situation the generic rule that lev els do not cross do es not apply and w e can easily imagine that the t wo low est levels lo ok like what we show in figure 3, where the levels actually cross at s ? . Imagine that low order p erturbation theory around s = 1 can b e used to get go o d appro ximations to the ground state energies of H A ( s ) and H B ( s ) for s near s ? , even for s somewhat to the left of s ? . Then it is p ossible to accurately predict s ? . Our situation is very close to this. W e can think of A as consisting of the states close to z L in Hamming weigh t, and B as states close to z U in Hamming w eight. Similarly , w e view H A and H B as the restrictions of H to these sectors. Note that it tak es n p o wers of the p erturbation to connect | 000 ... 0 i and | 111 ... 1 i and this is why we view A and B as essen tially disconnected. Although figure 3 lo oks like figure 2, it is figure 2 that depicts the actual situation, where the t wo levels av oid crossing. This true near cross means that the p erturbation series in the actual theory will diverge very close to s ? . How ev er this div ergence will only b e seen at high order, in fact at an order which is prop ortional to n . The low order terms of the p erturbation series in figure 3 1 The notation f ( n ) = Θ( g ( n )) means that, for n sufficiently large, bg ( n ) ≤ f ( n ) ≤ cg ( n ) for some constants b and c . 6 Figure 3: A true energy level crossing can arise from t wo disconnected sectors. are the same as the low order terms of the p erturbation series in figure 2, so w e can trust lo w order p erturbation theory to lo cate s ? . (W e argue b elo w that as a function of the num b er of bits, n , s ? go es to 1 as n go es to infinit y . This implies that the radius of conv ergence of the p erturbation theory for the full H ( s ) , expanded ab out s = 1 , go es to 0 as n go es to infinity . This fact has no b earing on our argument that low order p erturbation theory can b e used to accurately predict s ? . ) F or small v alues of the parameter 1 − s s , the energies of the tw o states under consideration can b e expanded as E L ( s ) = (1 − s ) n 2 + s " 0 +  1 − s s  2 e (2) L +  1 − s s  4 e (4) L + ... # (3) E U ( s ) = (1 − s ) n 2 + s " 1 2 +  1 − s s  2 e (2) U +  1 − s s  4 e (4) U + ... # . (4) It is easy to see that each expansion (inside the square brac k ets) only con tains even p o wers. Note that e (2) L is guaranteed to b e negative as is alwa ys the case for 2nd order p erturbation theory of the ground state. In addition, since we added a term to p enalize the state whic h had smaller energy near s = 1 b efore adding the clause, we exp ect that e (2) U < e (2) L . W e are interested in the b eha viour of the difference E L ( s ) − E U ( s ) for randomly generated instances (as generated using the prescription of the previous section) as a function of the num b er of spins n in the limit n → ∞ . This requires us to further in vestigate the behaviour of eac h coefficient e ( k ) L and e ( k ) U as a function of n and m . 7 In fact to lo cate s ? w e need only to go to second order in p erturbation theory . F rom equation 2 we view H P as the unperturb ed Hamiltonian, and − 1 2 P n i =1 σ i x as the perturbation, with 1 − s s as the expansion parameter. No w σ i x | z i = | z ⊕ ˆ e i i so at second order w e get e (2) L = 1 4 n X i =1 |h z L ⊕ ˆ e i | σ i x | z L i| 2 h z L | H P | z L i − h z L ⊕ ˆ e i | H P | z L ⊕ ˆ e i i = − 1 4 n X i =1 1 h z L ⊕ ˆ e i | H P | z L ⊕ ˆ e i i . Similarly , e (2) U = 1 4 n X i =1 1 1 2 − h z U ⊕ ˆ e i | H P | z U ⊕ ˆ e i i where the 1 2 in the denominator is h z U | H P | z U i . The exp ected energy p enalty incurred when flipping a bit of either z U or z L is of order m n since eac h bit is typically inv olved in Θ( m n ) clauses. So the co efficien ts e (2) L and e (2) U are of order n  n m  since the energy denominators in volv ed are Θ( m n ) . W e now show that their difference is of order √ n  n m  3 2 . W rite e (2) U − e (2) L = n X i =1 d i where for eac h i , d i = 1 4  1 1 2 − h z U ⊕ ˆ e i | H P | z U ⊕ ˆ e i i + 1 h z L ⊕ ˆ e i | H P | z L ⊕ ˆ e i i  . (5) Recall that H P = H 0 P + h where h is the p enalty term from the final clause which acts only on the first 3 bits. Therefore, for i = 4 , 5 , ..., n d i = 1 4  − 1 h z U ⊕ ˆ e i | H 0 P | z U ⊕ ˆ e i i + 1 h z L ⊕ ˆ e i | H 0 P | z L ⊕ ˆ e i i  . Our pro cedure for generating instances is symmetric b etw een the strings 000 ... 0 and 111 ... 1 so a veraging ov er instances it is clear that the mean of d i for i = 4 , 5 ..., n is 0 . Th us we exp ect P n i =4 d i to b e (appro ximately) Gaussian with mean 0 and standard deviation prop ortional to √ nσ ( d ) , where σ ( d ) is the standard deviation of each d i for i ∈ { 4 , 5 , ..., n } . T o compute σ ( d ) we note that d i = 1 4  h z U ⊕ ˆ e i | H 0 P | z U ⊕ ˆ e i i − h z L ⊕ ˆ e i | H 0 P | z L ⊕ ˆ e i i h z L ⊕ ˆ e i | H 0 P | z L ⊕ ˆ e i ih z U ⊕ ˆ e i | H 0 P | z U ⊕ ˆ e i i  . Again using the symmetry b etw een all zeros and all ones, we conclude that the n umerator is of order p m n and the denominator is of order  m n  2 . Hence we exp ect σ ( d ) to b e Θ   n m  3 2  . So e (2) U − e (2) L 8 is of order √ n  n m  3 2 . W e will no w lo cate s ? using second order p erturbation theory and afterw ards argue that higher orders do not change the result. Returning to equations 3 and 4, equating the t wo energies at second order we hav e 0 = s ? " 1 2 +  1 − s ? s ?  2  e (2) U − e (2) L  # so s ? = 1 − Θ( 1 n 1 / 4  m n  3 4 ) . (6) The 4th order correction to the energy of the lo wer state is given by e (4) L = h z L | V  φ L H P  2 V | z L ih z L | V φ L H P V | z L i − h z L | V φ L H P V φ L H P V φ L H P V | z L i where V = − 1 2 n X i =1 σ i x and φ L = 1 − | z L ih z L | . W riting H P ( z ) = h z | H P | z i , this can b e expressed as e (4) L = 1 16 n X i =1 n X j =1 1 ( H P ( z L ⊕ ˆ e i )) 2 H P ( z L ⊕ ˆ e j ) − 1 16 X i 6 = j 1 H P ( z L ⊕ ˆ e i ) H P ( z L ⊕ ˆ e j ) H P ( z L ⊕ ˆ e i ⊕ ˆ e j ) − 1 16 X i 6 = j 1 ( H P ( z L ⊕ ˆ e i )) 2 H P ( z L ⊕ ˆ e i ⊕ ˆ e j ) = n X i =1 1 16 1 ( H P ( z L ⊕ ˆ e i )) 3 + X i 6 = j 1 16 H P ( z L ⊕ ˆ e i ⊕ ˆ e j ) − H P ( z L ⊕ ˆ e i ) − H P ( z L ⊕ ˆ e j ) ( H P ( z L ⊕ ˆ e i )) 2 H P ( z L ⊕ ˆ e j ) ( H P ( z L ⊕ ˆ e i ⊕ ˆ e j )) . No w consider the terms in this expression corresp onding to indices i, j for whic h i 6 = j and the bits i and j do not app ear in a clause together. Under these conditions w e hav e H P ( z L ⊕ ˆ e i ⊕ ˆ e j ) = H P ( z L ⊕ ˆ e i ) + H P ( z L ⊕ ˆ e j ) . So w e can write e (4) L = n X i =1 1 16 1 ( H P ( z L ⊕ ˆ e i )) 3 (7) + X i 6 = j clausemates 1 16 H P ( z L ⊕ ˆ e i ⊕ ˆ e j ) − H P ( z L ⊕ ˆ e i ) − H P ( z L ⊕ ˆ e j ) ( H P ( z L ⊕ ˆ e i )) 2 H P ( z L ⊕ ˆ e j ) ( H P ( z L ⊕ ˆ e i ⊕ ˆ e j )) . 9 Here the subscript “clausemates” indicates that we sum only ov er pairs of indices whic h app ear together in at least one clause of the 3SA T instance corresp onding to H P . F or e (4) U , since the unp erturb ed energy is 1 2 w e obtain e (4) U = n X i =1 1 16 1  H P ( z U ⊕ ˆ e i ) − 1 2  3 (8) + X i 6 = j clausemates 1 16  H P ( z U ⊕ ˆ e i ⊕ ˆ e j ) − 1 2  −  H P ( z U ⊕ ˆ e i ) − 1 2  −  H P ( z U ⊕ ˆ e j ) − 1 2   H P ( z U ⊕ ˆ e i ) − 1 2  2  H P ( z U ⊕ ˆ e j ) − 1 2   H P ( z U ⊕ ˆ e i ⊕ ˆ e j ) − 1 2  . Let’s lo ok at the first sum in each of equations 7 and 8 where i go es from 1 to n . Eac h term as i go es from 4 to n is of order  n m  3 and so the difference of the tw o sums is Θ( √ n  n m  3 ) . The second sums (those which are restricted to clausemates) con tain of order m terms. Each denominator is of order  m n  4 and the n umerators are Θ(1) since the only contribution to the numerator is from clauses in H P whic h in v olve bits i and j together. Separating out the terms where i and/or j are 1 , 2 or 3 , we conclude that the con tribution to e (4) U − e (4) L from the clausemate terms is Θ( √ n  n m  7 2 ) . F or our instance generation m grows lik e n log n so this clausemate con tribution is asymptotically dominated by the first term whic h scales lik e Θ( √ n  n m  3 ) . So the fourth order con tribution to the difference of energies E U ( s ) − E L ( s ) is Θ( s h √ n  n m  3  1 − s s  4 i ) . At s ? whic h we determined at second order to b e 1 − Θ( 1 n 1 / 4  m n  3 4 ) , the fourth order con tribution to the energy difference is Θ  1 √ n  . The fourth order corrections can therefore b e neglected in determining the lo cation of s ? . Sixth order and higher con tributions to the difference are even smaller. I I I. FIXING THE PR OBLEM BY P A TH CHANGE In our instance generation we manufactured a small gap by p enalizing the planted assignment corresp onding to the energy eigenstate with the smallest energy near s = 1 . Since the slop es of the t wo curves in figure 1 are the same at s = 1 , the second deriv atives determine which eigen v alue is smaller near s = 1 . After p enalization we hav e e (2) U < e (2) L whic h is consistent with the near crossing in figure 2. Supp ose instead that we p enalize the assignment corresp onding to larger energy in figure 1. Then w e exp ect the situation depicted in figure 4 where no level crossing is induced. W e imagine that the instances that we manufacture with a small gap as in figure 2 are a mo del for what migh t b e encountered in running the quantum adiabatic algorithm on some instance a quan tum computer is actually trying to solve. There is a strategy for ov ercoming this problem. The idea is to pro duce figure 4, with reasonable probabilit y , b y randomly mo difying the adiabatic path whic h ends at H P , of course making no use of the prop erties of the particular instance. F or this 10 Figure 4: Energy levels with no a voided crossing near s = 1 . Here the second deriv ative of the upp er curve B is greater than the second deriv ative of the low er curv e A. purp ose one could use an y random ensemble of paths H ( s ) such that H (1) = H P and the ground state of H (0) is simple to prepare. How ever in this pap er w e only consider randomly changing the b eginning Hamiltonian. W e ha ve made this choice so that we are able to use our Quantum Monte Carlo metho d to numerically verify our arguments. Lik e other Quantum Mon te Carlo metho ds, the metho d we use do es not work when the Hamiltonian has nonzero off-diagonal matrix elements with p ositiv e sign. Starting with an instance where H ( s ) = (1 − s ) H B + sH P has a tiny gap due to the problem discussed ab o v e, w e now consider a differen t adiabatic path ˜ H ( s ) = (1 − s ) ˜ H B + sH P obtained b y k eeping the same problem Hamiltonian H P but c ho osing a different b eginning Hamiltonian ˜ H B in a random fashion whic h w e will prescrib e b elo w. W e argue that the small gap near s = 1 is then remo ved with substan tial probability , so that by rep eating this procedure a constant n umber of times it is p ossible to find an adiabatic path without a small gap near s = 1 . The wa y that we c ho ose a random Hamiltonian ˜ H B is to first dra w n random v ariables c i for i = 1 , 2 , ..., n , where each c i is c hosen to b e 1 2 or 3 2 with equal probabilit y . W e then take ˜ H B = n X i =1 c i (1 − σ i x ) 2 . (9) W e write ˜ e (2) U and ˜ e (2) L for the analogous quantities to e (2) L and e (2) U for the new Hamiltonian ˜ H ( s ) = (1 − s ) ˜ H B + sH P . The p oin t is that b y randomizing ˜ H B in the w ay we prescrib e, there is a 11 substan tial probability that one will obtain ˜ e (2) U − ˜ e (2) L > 0 , and in that case one exp ects no av oided crossing near s = 1 . W rite ˜ e (2) U − ˜ e (2) L = n X i =1 c 2 i d i where the { d i } are fixed b y the instance (and are defined in equation 5). Since w e hav e fixed the problem Hamiltonian H P , the only random v ariables app earing in the ab o ve equation are the c i . W e ha ve c 2 i = 5 4 so the mean v alue of ˜ e (2) U − ˜ e (2) L is then ˜ e (2) U − ˜ e (2) L = 5 4 ( e (2) U − e (2) L ) < 0 . But more imp ortan tly ˜ e (2) U − ˜ e (2) L = Θ( √ n  n m  3 2 ) . The v ariance of this difference is V ar  ˜ e (2) U − ˜ e (2) L  = n X i =1 d 2 i V ar ( c 2 i ) = n X i =1 d 2 i · 1 , whic h is Θ( n  n m  3 ) . F or a fixed instance with a corresp onding fixed set { d i } the random v ariable P i c 2 i d i is appro ximately Gaussian and from its mean and v ariance w e see that the probabilit y that ˜ e (2) U − ˜ e (2) L is positive and in fact greater than a √ n  n m  3 2 , for a > 0 , is b ounded aw a y from 0 indep enden t of n . This means that there is a go o d chance that randomizing H B turns the situation depicted in figure 2 in to the situation depicted in figure 4. In the case of tw o planted satisfying assignments with one p enalized to pro duce a small gap when the b eginning Hamiltonian is H B of equation 1, we ha ve shown that a random c hoice for the b eginning Hamiltonian ˜ H B of equation 9 can with substantial probabilit y remov e the small gap. This gives further weigh t to the idea that when running the quan tum adiabatic algorithm on a single instance of some optimization problem, the programmer should run the quan tum adiabatic algorithm rep eatedly with differen t paths ending at H P [5]. What if there are k > 2 satisfying assignmen ts and all but one are p enalized? W e ha ve considered instances of 3SA T for which the corresp onding problem Hamiltonian has a unique ground state and a nondegenerate first excited state that is far from the ground state in 12 Hamming weigh t. In this case w e hav e shown that the path change strategy succeeds in removing a tiny gap. What happ ens when w e ha ve k m utually disparate solutions and we p enalize all but one of them? W e sho w (under some assumptions) that for large n , path change will succeed after a num b er of tries polynomial in k . W e can therefore hop e for success if there are k disparate assignmen ts which violate few clauses and k scales p olynomially with n . On the other hand when k is sup erp olynomial in n we ha v e no reason to b e optimistic ab out the p erformance of the quantum adiabatic algorithm. The energies E r ( s ) for r = 0 , 1 , .., k − 1 of the k states under consideration can b e expanded as in equations 3 and 4. W rite e (2) r for the second order corrections in these expansions. If we add a small num b er of clauses to this instance which penalize all the solutions except for a randomly c hosen one which we call z 0 (forming a new problem Hamiltonian ˆ H p ) then the differences e (2) r − e (2) 0 will not c hange substantially . So we are interested in the differences e (2) r − e (2) 0 = n X i =1 d ri , where d ri = 1 4  − 1 h z r ⊕ ˆ e i | H p | z r ⊕ ˆ e i i + 1 h z 0 ⊕ ˆ e i | H P | z 0 ⊕ ˆ e i i  , (10) and where z 0 , z 1 , z 2 , ..., z k − 1 are the k ground states of H p . The adiabatic algorithm applied to this instance with problem Hamiltonian ˆ H p will succeed if all of the ab o v e differences are p ositive, for r = 1 , ..., k − 1 . There is a 1 k c hance of this o ccurring. If they are not all p ositiv e, then we ha ve encountered an instance whic h we exp ect to ha ve a small gap. In this case, w e randomize the b eginning Hamiltonian as in the previous subsection. This pro duces new v alues ˜ e (2) r for the second order energy corrections and the relev ant differences are given by ˜ e (2) r − ˜ e (2) 0 = n X i =1 c 2 i d ri , (11) where the random v ariables c i are 1 2 or 3 2 with equal probabilit y (and the d ri are fixed). With this fixed problem Hamiltonian, how many times do w e ha v e to randomize the b eginning Hamiltonian (b y dra wing a random set { c i } ) to pro duce algorithmic success? F or large n , the random v ariables ˜ e (2) r − ˜ e (2) 0 are appro ximately jointly Gaussian with exp ectation and v ariance given by ˜ e (2) r − ˜ e (2) 0 = 5 4 n X i =1 d ri V ar  ˜ e (2) r − ˜ e (2) 0  = n X i =1 ( d ri ) 2 . 13 The correlation (the cov ariance divided by the pro duct of the standard deviations) b et w een any t wo of these random v ariables is given by Corr  ˜ e (2) q − ˜ e (2) 0 , ˜ e (2) r − ˜ e (2) 0  = P i d q i d ri h P j ( d q j ) 2 P j ( d rj ) 2 i 1 2 . (12) W e now argue that the ab o ve correlation is typically equal to 1 2 when n is large. T o see this, we can expand the ab ov e expression using equation 10. W e imagine that the instances hav e b een generated so that the la w of large num b ers applies and that in the limit of large n the quantit y 1 n n X i =1  1 h z r ⊕ ˆ e i | H p | z r ⊕ ˆ e i i · 1 h z r ⊕ ˆ e i | H p | z r ⊕ ˆ e i i  approac hes an r indep enden t v alue and the quantit y 1 n n X i =1  1 h z q ⊕ ˆ e i | H p | z q ⊕ ˆ e i i · 1 h z r ⊕ ˆ e i | H p | z r ⊕ ˆ e i i  approac hes a v alue indep endent of r and q for r 6 = q . Using this in equation 12, we obtain that the correlation is equal to 1 2 (for r 6 = q and in the limit of large n ). Recall that the problem Hamiltonian has b een fixed and therefore the { d ri } are set. W e ha ve argued that for large n the differences n ˜ e (2) r − ˜ e (2) 0 o ha ve appro ximately a join t normal distribution with pairwise correlations equal to 1 2 . The probability that a random c hoice of the { c i } mak e the algorithm succeed is Pr h ˜ e (2) r − ˜ e (2) 0 > 0 , r = 1 , ..., k − 1 i = Pr " n X i =1 c 2 i d ri > 0 , r = 1 , ..., k − 1 # = Pr   P n i =1 c 2 i d ri − 5 4 P n i =1 d ri q P n i =1 ( d ri ) 2 > − 5 4 P n i =1 d ri q P n i =1 ( d ri ) 2 , r = 1 , ..., k − 1   . Note that in the last line we hav e subtracted the mean and divided b y the standard deviation. Using a standard tric k we can rewrite the ab ov e probability as (in the limit of large n), Pr   w r − w 0 > − 5 4 P n i =1 d ri q P n i =1 ( d ri ) 2 , r = 1 , ..., k − 1   , (13) where w 0 and { w r } are indep enden t normal random v ariables with mean zero and v ariance 1 2 . This follo ws from the fact that the joint normal distribution is sp ecified by its marginal distributions and pairwise correlations. Essentially w r and w 0 are the (normalized) changes in the quantities e (2) r and e (2) 0 pro vided by the randomization. A t this p oint, in order to ev aluate equation 13 we need to understand the magnitude of the RHS of the inequality . Consider typical v alues for the sum P n i =1 d ri . As in section I I we assume that 14 this sum is Gaussian in the limit of large n with mean 0 and standard deviation appro ximated by q P n i =1 d 2 ri . (In the argument b elo w, we use some prop erties of the Gaussian distribution whic h follow from the inequalities Z ∞ a e − x 2 dx ≤ Z ∞ a x a e − x 2 dx = 1 2 a e − a 2 Z ∞ a e − x 2 dx ≥ Z a +1 a e − x 2 dx ≥ e − ( a +1) 2 . ) F or the k − 1 approximately Gaussian random v ariables − P n i =1 d ri √ P n i =1 ( d ri ) 2 with mean 0 and v ariance 1 , for some constan t B it will hold with high probabilit y (in the limit of large n ) that − 5 4 P n i =1 d ri q P n i =1 ( d ri ) 2 < B p log k for r = 1 , ..., k − 1 . No w the random v ariable w 0 satisfies w 0 < − 2 B √ log k with probability whic h is only p olynomially small in k . F or eac h w r the probability that w r < − B √ log k is p olynomially small in k and we use the b ound Pr h w r − w 0 ≥ B p log k , r = 1 , ..., k − 1 i ≥ Pr h w 0 < − 2 B p log k i · k − 1 Y r =1  1 − Pr h w r < − B p log k i . The first term on the RHS is only p olynomially small, and eac h of the k − 1 terms in the pro duct are at least 1 − 1 k 2 if B is large enough. So the probability that the adiabatic algorithm succeeds is at w orst p olynomially small in k . By rep eating the randomization a p olynomial num b er of times w e exp ect to succeed with high probabilit y . IV. QUANTUM MONTE CARLO AND NUMERICAL RESUL TS Con tin uous Imaginary Time Quan tum Mon te Carlo This section is a review of con tinous imaginary time Quantum Monte Carlo (whic h is a classical path in tegral simulation technique for extracting prop erties of quantum systems [16]). In particular w e will show ho w this metho d can b e used to compute thermal exp ectation v alues of Hermitian op erators at in v erse temp erature β . W e start with a Hamiltonian H which we write as H = H 0 + V where H 0 is diagonal in some known basis {| z i} , and V is purely off diagonal in this basis. W e require that all the nonzero matrix elemen ts of V are negativ e. F or the Hamiltonian H ( s ) = (1 − s ) H B + sH P 15 with H B as in equation 1, w e hav e H 0 = sH P + (1 − s ) n 2 V = − (1 − s ) n X i =1 σ i x 2 . Here we include the factor of 1 − s in the definition of V since w e are not doing p erturbation theory in this quan tity . The partition function can b e expanded using the Dyson series as follo ws T r h e − β H i = T r  e − β H 0 ∞ X m =0 ( − 1) m Z β t m =0 dt m Z t m t m − 1 =0 dt m − 1 ... Z t 2 t 1 1 =0 dt 1 V I ( t m ) ...V I ( t 1 )  . Here we use the notation V I ( t ) = e tH 0 V e − tH 0 . W e can then insert complete sets of states and take the trace to obtain the path in tegral T r h e − β H i = ∞ X m =0 X { z 1 ,...,z m }  ( − 1) m h z 1 | V | z m ih z m | V | z m − 1 i ... h z 2 | V | z 1 i Z β t m =0 dt m Z t m t m − 1 =0 dt m − 1 ... Z t 2 t 1 =0 dt 1 e − R β t =0 H 0 ( z ( t )) dt  . (14) In this form ula we hav e used the notation H 0 ( z ( t )) = h z ( t ) | H 0 | z ( t ) i , where the function z ( t ) is defined b y z ( t ) =                            z 1 , 0 ≤ t < t 1 z 2 , t 1 ≤ t < t 2 . . . z m , t m − 1 ≤ t < t m z 1 , t m ≤ t ≤ β . So in particular Z β t =0 H 0 ( z ( t )) dt = h z 1 | H 0 | z 1 i ( t 1 + β − t m ) + h z 2 | H 0 | z 2 i ( t 2 − t 1 ) + ... h z m | H 0 | z m i ( t m − t m − 1 ) . W e view the function z ( t ) as a path in imaginary time which b egins at t = 0 and ends at t = β . Then equation 14 is a sum ov er paths, where ev ery path is assigned a p ositive weigh t according to a measure ˜ ρ ˜ ρ ( P ) = ( − 1) m h z 1 | V | z m ih z m | V | z m − 1 i ... h z 2 | V | z 1 i dt 1 ...dt m e − R β t =0 H 0 ( z ( t )) dt . 16 Note that the fact that ˜ ρ is p ositiv e semidefinite follo ws from our assumption that all matrix elemen ts of V are negative. W e write ρ = ˜ ρ Z ( β ) (15) for the normalized distribution o ver paths. W e now discuss how one can obtain prop erties of the quantum system as exp ectation v alues with resp ect to the classical measure ρ. W e are interested in ground state prop erties, but in practice w e select a large v alue of β and compute thermal exp ectation v alues at this inv erse temp erature. If β is sufficien tly large then these expectation v alues will agree with the corresp onding exp ectation v alues in the ground state. It is straightforw ard to show that for an y op erator A whic h is diagonal in the | z i basis, the thermal exp ectation v alue can b e written as T r [ Ae − β H ] T r [ e − β H ] = h 1 β Z β 0 A ( z ( t )) dt i ρ (16) where A ( z ( t )) = h z ( t ) | A | z ( t ) i . The notation hi ρ means av erage with resp ect to the classical proba- bilit y distribution ρ ov er paths. This immediately allows us to estimate quan tities suc h as the diagonal part of the Hamiltonian T r [ H 0 e − β H ] T r [ e − β H ] = h 1 β Z β 0 H 0 ( z ( t )) dt i ρ . (17) Another quan tity that we find useful in our study is the Hamming weigh t op erator defined by W = n X i =1  1 − σ i z 2  . (18) whic h can also b e estimated using equation 16. In order to estimate the thermal exp ectation v alue of the full Hamiltonian H 0 + V we use equation 17 as well as the expression T r [ V e − β H ] T r [ e − β H ] = −h m β i ρ (19) where m is the num b er of transitions in the path. (Equation 19 is not as simple to derive as equation 16.) So b y generating paths from the distribution ρ and then computing av erages with respect to ρ , we can ev aluate the thermal exp ectation v alue of the energy at a given inv erse temp erature β and b y taking β sufficien tly large we can approximate the ground state energy . Generating paths from the distribution ρ is itself a challenging task. W e use a mo dified version of the heat bath algorithm of Krzak ala et al [14] whic h we describ e in the app endix. As with other Quan tum Monte Carlo metho ds, this algorithm is a Mark o v Chain Monte Carlo metho d. In order to sample from the distribution ρ , one defines a Mark ov Chain ov er the state space consisting of 17 all paths, where the limiting distribution of the c hain is ρ . T o obtain samples from ρ one starts in an initial path z init ( t ) and then applies some num b er N eq uil of iterations of the Mark o v Chain. If N eq uil is sufficiently large then the distribution of the paths found by further iteration will b e close to ρ. W e use these subsequen t paths to compute av erages with resp ect to ρ . Equilibration of the Quan tum Monte Carlo and Iden tification of Level Crossings The discussion in the previous section demonstrates how one can estimate exp ectation v alues of v arious quantities in the ground state of a quantum system. The metho d is “contin uous time” so there is no discretization error and for fixed β quan tities can b e arbitrarily well appro ximated with enough statistics. W e use the Quantum Monte Carlo metho d in a nonstandard wa y in order to b e able to study the low est tw o eigenstates of our Hamiltonian (as opp osed to just the ground state). As describ ed in the previous section, the standard pro cedure for generating configurations is to equilibrate to the distribution ρ from some initial path z init ( t ) (we call this the seeded path) b y applying the Mon te Carlo update N eq uil times. In order to ensure that one has equilibrated after N eq uil Mon te Carlo up dates, one could for example do sim ulations with tw o or more differen t initial paths (seeded paths) and c hec k that the v alues of Mon te Carlo observ ables app ear to conv erge to the same seed indep enden t v alues. F or eac h instance we consider, we run tw o differen t Monte Carlo simulations that are seeded with t wo different paths z 0 init ( t ) and z 1 init ( t ) . These seeds are paths with no flips in them, corresp onding to the t wo states 000 ... 0 and 111 ... 1 , z 0 init ( t ) = 000 ... 0 for all t ∈ [0 , β ] z 1 init ( t ) = 111 ... 1 for all t ∈ [0 , β ] . With eac h seed, we run the Mon te Carlo sim ulation for some total n um b er N total of Mon te Carlo sw eeps, taking data ev ery k sw eeps. Here a single sweep is defined to b e n iterations of our Monte Carlo up date rule as describ ed in the app endix (where n is the n umber of spins). W e then remov e the first N eq uil Mon te Carlo samples, and use the remaining Monte Carlo samples to estimate the thermal a verages h H i = T r [ H e − β H ] T r [ e − β H ] h W i = T r [ W e − β H ] T r [ e − β H ] using equations 17,19 and 18. 18 In order to con vince the reader that the Mon te Carlo algorithm works correctly , and that we can obtain information ab out the first excited state, we show data at 16 bits where exact numerical diagonalization is p ossible. In figures 5 and 6 w e show the result of this pro cedure for a 3SA T instance with 16 bits, with problem Hamiltonian H 0 P (b efore adding the p enalt y term h , so that the lev els are degenerate at s = 1 ). This instance has 122 clauses. The inv erse temp erature is β = 150 . The total num b er of Monte Carlo sweeps at each v alue of s is N total = 200000 , with data tak en after ev ery fifth sweep (this gav e us 40000 data samples). W e use N eq uil = 2500 samples solely for equilibration at each v alue of s . Note that the Mon te Carlo sim ulations with the t w o different seeds (corresp onding to the circles and crosses in the figures) only agree for v alues of s less than roughly 0 . 4 . W e interpret this to mean that at these v alues of s the Quantum Monte Carlo has equilibrated to the prop er limiting distribution ρ regardless of the seed. As s increases past 0 . 4 , the tw o sim ulations abruptly b egin to giv e differen t results. In this case the sim ulation seeded with z 1 init ( t ) finds a metastable equilibrium of the Mark ov Chain (i.e not the limiting distribution ρ ), whic h corresp onds in this case to the first excited state of the Hamiltonian. F or s ab ov e 0 . 4 we can see from the comparison with exact diagonalization that the tw o differently seeded v alues allow us to compute the energies (figure 5) and Hamming weigh ts (figure 6) of the low est tw o energy levels of our Hamiltonian. What is happ ening here is that for large s , the quan tum system can b e thought of as consisting of tw o disconnected sectors. One sector consists of states in the z basis with low Hamming w eight and the other with Hamming w eigh t near n . The Quan tum Monte Carlo “equilibrates” in eac h sector dep ending on the initial seed as can b e seen b y the smo oth data in eac h sector. (Note that data is taken indep endently at eac h v alue of s.) The low er of the cross and circle at eac h v alue of s in figure 5 is the ground state energy . Of course if the classical algorithm ran for long enough then the circles in figure 5 w ould lie on the crosses for every v alue of s . In figures 7 and 8 we show data taken for the same instance after the p enalty clause is added that remo ves the degeneracy at s = 1 . In this case the tw o lev els hav e an av oided crossing as exp ected and whic h can b e seen by exact n umerical diagonalization in figure 7. How ever in the Monte Carlo simulation with tw o seeds there are tw o essentially disconnected sectors and the tw o levels b eing track ed do cross. When we see this b eha viour in the Monte Carlo simulation, we in terpret it as comp elling evidence that there is a tiny gap in the actual system, whic h o ccurs where the curv es cross. In figure 8 first lo ok at the curve s coming from the exact numerical diagonalization. The Hamming weigh t of the ground state decreases and then abruptly rises at the lo cation of the minim um gap. The Hamming weigh t of the first excited state also undergo es an abrupt change 19 at the location of the minimum gap. Not surprisingly the exact diagonalization clearly sho ws the b eha viour w e exp ect with the manufactured tiny gap. Now lo ok at the Monte Carlo data which is in terpreted by first lo oking at figure 7. In figure 7 the true ground state is alwa ys the low er of the circles and the crosses. F or s b elow s ? ≈ 0 . 423 it is the crosses and for s larger than s ? it is the circles. A ccordingly in figure 8 the Hamming weigh t of the ground state is trac ked by the crosses for s b elow s ? and by the circles for s ab o v e s ? . W e can therefore conclude, based only on the Mon te Carlo data, that the Hamming w eight of the ground state changes abruptly . A t higher bit n umber w e do not ha ve exact numerical diagonalization a v ailable but we can still use the Quan tum Mon te Carlo to extract key information. Randomizing the Beginning Hamiltonian F or the 16 spin Hamiltonian depicted ab ov e, we now consider randomizing the b eginning Hamil- tonian H B as describ ed in section I I I. In order to minimize the num b er of times we run the Mon te Carlo algorithm (this will b e more of an issue at high bit n um b er where simulations are v ery time consuming), w e generated many differen t sets of co efficients { c i } and calculated the differences ˜ e (2) U − ˜ e (2) L with fixed problem Hamiltonian H 0 P and b eginning Hamiltonians ˜ H B = n X i =1 c i  1 − σ i x 2  . A ccording to our discussion in section I I I, w e exp ect the a v oided crossing near s = 1 to b e remo ved for choices of co efficients c i suc h that the difference ˜ e (2) U − ˜ e (2) L > 0 . W e made a histogram of the ˜ e (2) U − ˜ e (2) L (sho wn in figure 9) and randomly chose three sets of co efficients { c i } suc h that ˜ e (2) U − ˜ e (2) L > 1 2 for each. Our analysis predicts that in these cases the crossing will b e remov ed. As exp ected, each of these three sets of co efficien ts resulted in an adiabatic Hamiltonian with the small gap remov ed, although in one there was another small gap at a smaller v alue of s . W e plot the Mon te Carlo data for one of these sets in figures 10 and 11. In section I I w e argued that w e could manufacture a small gap in a 3SA T instance as sketc hed in figure 2. Our 16 bit instance is a concrete example of this and the tiny gap can b e seen in figure 7. In section I I I w e argued that by randomizing the b eginning Hamiltonian we should be able to pro duce figure 4. This is what w e see concretely at 16 bits in figure 10. W e now tell the same story at 150 bits using only Mon te Carlo data since exact numerical diagonalization is not p ossible. 20 Data for an Instance of 3SA T with 150 bits In addition to v alidating our metho d at 16 bits, we studied 3SA T instances with 25 , 75 and 150 bits using our Quan tum Mon te Carlo simulator. The data from these simulations for the most part supp orted our arguments. In this section we presen t Monte Carlo data taken for a double plant instance of 3SA T with 150 spins and 1783 clauses. In this case the inv erse temp erature β = 300 , and the total n umber of Monte Carlo sweeps at eac h v alue of s is N total = 100000 , with data tak en ev ery fifth sw eep (giving 20000 data samples). The first 2500 data samples at each v alue of s are remo ved for equilibration. W e ran our simulations on a 648 pro cessor SiCortex computer cluster in em barrassingly parallel fashion. W e used a differen t pro cessor for eac h v alue of s and each v alue of the seed. Data taken at the low er v alues of s to ok the longest to accum ulate, in some cases more than 10 da ys on a single pro cessor for a single data p oin t. Figures 12 and 13 show the energy and the Hamming weigh t b efore the p enalt y clause is added, running with t wo different seeds. In figure 12 for large v alues of s , the crosses are alw ays b elow the circles and track the ground state which ends at 000 ... 0 as can b e seen in figure 13. W e ha ve also plotted the second order p erturbation theory energies for these tw o levels, expanding around s = 1 . Note the go o d agreement b et w een the second order p erturbation theory and the Monte Carlo data all the w ay down to s = 0 . 35 . Since the lo wer curv e corresp onds to 000 ... 0 , we p enalize this assignment by the addition of a single clause to form H P attempting to man ufacture a near crossing. In figure 14, the circle data is b elo w the cross data for s near 1 but the t w o curv es cross near s = . 49 . This is seen clearly in figure 15 where the energy difference is plotted. F rom the Monte Carlo data we conclude that H ( s ) has a tiny gap. The lo cation of the av oided crossing is w ell predicted by second order p erturbation theory as can b e seen in figure 15. The Mon te Carlo data at 150 bits shows that w e can make an instance of 3SA T with a tiny gap following the pro cedure outlined in section I I. W e now use the Monte Carlo sim ulator to show that a randomly chosen ˜ H B can alter the Hamiltonian H ( s ) so that this small gap b ecomes large. F or the instance at hand we first compute the difference ˜ e (2) U − ˜ e (2) L for 100000 randomly chosen sets of coefficients. The histogram of these differences is plotted in figure 17. After doing this we randomly selected (from this set) tw o sets of co efficien ts such that ˜ e (2) U − ˜ e (2) L > 1 2 . F or b oth of these sets of co efficien ts we saw that the crossing at s ≈ 0 . 49 was no longer present, although in one case there app eared to b e a new crossing at a m uc h low er v alue of s . Figures 18 and 19 show the Monte Carlo data for the choice of ˜ H B whic h do es not app ear to hav e an y crossing. At 150 bits w e see 21 comp elling evidence that the story outlined in sections I I and II I is true. Is s ? near 1 ? W e argued that s ? should go to 1 for large enough n . Ho w ev er our 16 bit example has s ? ≈ 0 . 42 and the 150 bit example has s ? ≈ 0 . 49 . Recall from equation 6 that s ? = 1 − Θ( 1 n 1 / 4  m n  3 4 ) . A t 16 bits with m = 122 we hav e 1 n 1 / 4  m n  3 4 = 2 . 29 and at 150 bits with m = 1783 we ha ve 1 n 1 / 4  m n  3 4 = 1 . 83 . Although asymptotically m is of order n log n , for these v alues of n w e are not y et in the regime where 1 n 1 / 4  m n  3 4  1 . Ev en though s ? is not near 1 for our instance at 150 bits, w e see from figure 15 that second order p erturbation theory can b e used to predict the lo cation of s ? . This is b ecause the fourth order con tribution to the energy difference is quite small. So already at 150 bits w e can predict the presence or absence of an a voided crossing using second order p erturbation theory . V. CONCLUSIONS W e ha ve introduced a new Quan tum Mon te Carlo tec hnique to analyze the p erformance of quan tum adiabatic algorithms for instances of satisfiability . Using seeded configurations (in the Mon te Carlo simulation) corresp onding to disparate low lying states, our technique exp osed the presence or absence of an exp onen tially small gap without ever actually computing the gap. W e used this metho d to n umerically inv estigate a set of random instances of 3SA T which w ere designed to exp ose a w eakness of the adiabatic algorithm. W e confirmed that this w eakness can b e o vercome for our set of instances b y using path change. Our n umerical work and the main part of our analysis p ertains to instances of 3SA T with 2 plan ted satisfying assignments. W e ha v e also considered the scenario where an instance has k satisfying assignments and then all but one are penalized with a small num b er of clauses. Ho w ev er, in the k > 2 case w e hav e made certain assumptions which mak e our analysis p ossible and we do not know if they apply to randomly generated instances. In our scenario the adiabatic algorithm with path change will succeed in a num b er of tries whic h is p olynomially large in k . Here k is fixed and n go es to infinity but we take the calculation in section II I as an indication that the algorithm will succeed when k gro ws p olynomially with n . Throughout this pap er we hav e assumed that our instances hav e a unique satisfying assignmen t. In this case w e b eliev e that the crucial distinction is whether there are p olynomially man y or 22 exp onen tially many low lying m utually disparate states abov e the unique satisfying assignmen t. If there are p olynomially many , we hav e argued that the adiabatic algorithm with path c hange will succeed, but if there are exp onentially many we ha ve no reason to b e optimistic ab out the p erformance of the algorithm. When there are exp onentially many satisfying assignments our analysis do es not apply . This situation w as considered in recent work by Knysh and Smelyanskiy [13]. Our results giv e further evidence that path change m ust b e considered an in tegral part of the quan tum adiabatic algorithm. F or an y giv en instance, the algorithm should b e run with many differen t randomly selected paths which end at the problem Hamiltonian. As long as the algorithm succeeds on at least a p olynomially small fraction of the trials, it can b e used to solv e decision problems. VI. A CKNO WLEDGEMENTS W e would like to thank Mohammad Amin, Y ale F an, Floren t Krzak ala, Jeremie Roland, and P eter Y oung for interesting discussions. W e also thank Alan Edelman for generously offering us access to his SiCortex computer cluster and Andy Lutomirski for fixing it. This work was sup- p orted in part by funds provided by the U.S. Departmen t of Energy under co op erativ e research agreemen t DE-FG02-94ER40818, the W. M. Keck F oundation Center for Extreme Quan tum Infor- mation Theory , the U.S Army Researc h Lab oratory’s Army Researc h Office through grant num b er W911NF-09-1-0438, the National Science F oundation through grant n umber CCF-0829421, and the Natural Sciences and Engineering Researc h Council of Canada. [1] Dorit Aharonov, Wim v an Dam, Julia Kemp e, Zeph Landau, Seth Lloyd, and Oded Regev. Adia- batic quantum computation is equiv alent to standard quan tum computation. SIAM JOURNAL OF COMPUTING , 37:166, 2007. arXiv:quant-ph/0405098. [2] Boris Altshuler, Hari Kro vi, and Jeremie Roland. Adiabatic quan tum optimization fails for random instances of NP-complete problems, 2009. [3] M. H. S. Amin and V. Choi. First order quantum phase transition in adiabatic quantum computation, 2009. [4] Y ale F an. Adiabatic quantum algorithms for b o olean satisfiability , 2009. (RSI summer pro ject at MIT). [5] Edward F arhi, Jeffrey Goldstone, and Sam Gutmann. Quantum adiabatic evolution algorithms with differen t paths, 2002. arXiv:quant-ph/0208135. 23 [6] Edward F arhi, Jeffrey Goldstone, Sam Gutmann, Joshua Lapan, Andrew Lundgren, and Daniel Preda. A quan tum adiabatic evolution algorithm applied to random instances of an NP-complete problem. Scienc e , 292:472–475, 2001. arXiv:quant-ph/0104129. [7] Edward F arhi, Jeffrey Goldstone, Sam Gutmann, and Daniel Naga j. Ho w to make the quan tum adiabatic algorithm fail. International Journal of Quantum Information , 6:503, 200 8. arXiv:quant- ph/0512159. [8] Edward F arhi, Jeffrey Goldstone, Sam Gutmann, and Michael Sipser. Quantum computation b y adia- batic ev olution, 2000. arXiv:quant-ph/0001106. [9] Daniel Fisher. Priv ate communication based on references [10] and [11]. [10] Daniel S. Fisher. Critical behavior of random transv erse-field ising spin c hains. Phys. R ev. B , 51(10):6411–6461, Mar 1995. [11] Daniel S. Fisher and A. P . Y oung. Distributions of gaps and end-to-end correlations in random transv erse-field ising spin chains. Phys. R ev. B , 58(14):9131–9141, Oct 1998. [12] T ad Hogg. A diabatic quan tum computing for random satisfiability problems. Physic al R eview A , 67:022314, 2003. quant-ph/0206059. [13] S. Knysh and V. Smely anskiy. On the relev ance of av oided crossings a wa y from quan tum critical p oint to the complexit y of quantum adiabatic algorithm. May 2010. [14] Florent Krzak ala, Alb erto Rosso, Guilhem Semerjian, and F rancesco Zamp oni. On the path integral represen tation for quantum spin mo dels and its application to the quan tum cavit y metho d and to monte carlo sim ulations. Physic al R eview B , 78:134428, 2008. [15] Y oshiki Matsuda, Hidetoshi Nishimori, and Helm ut G. Katzgrab er. Ground-state statistics from annealing algorithms: Quan tum vs classical approaches. New Journal of Physics , 11:073021, 2009. [16] N. V. Prokof ’ev, B. V. Svistunov, and I. S. T upitsyn. Exact quantum monte carlo pro cess for the statistics of discrete systems. ZH.EKS.TEOR.FIZ. , 64:853, 1996. arXiv:cond-mat/9612091. [17] Ben W. Reic hardt. The quantum adiabatic optimization algorithm and lo cal minima. In STOC ’04: Pr o c e e dings of the thirty-sixth annual A CM symp osium on The ory of c omputing , pages 502–510, New Y ork, NY, USA, 2004. ACM. [18] Wim v an Dam, Michele Mosca, and Umesh V azirani. How p o werful is adiabatic quantum computation?, 2002. arXiv:quant-ph/0206003. [19] Wim v an Dam and Umesh V azirani. Limits on quantum adiabatic optimization (unpublished man uscript). [20] A. P . Y oung, S. Kn ysh, and V. N. Smelyanskiy . Size dep endence of the minim um excitation gap in the quan tum adiabatic algorithm. Physic al R eview L etters , 101:170503, 2008. [21] A.P Y oung. Priv ate communication, 2009. [22] Marko Znidaric and Martin Horv at. Exp onen tial complexity of an adiabatic algorithm for an NP- complete problem. Physic al R eview A , 73:022329, 2006. arXiv:quant-ph/0509162. 24 App endix: A mo dified V ersion of the Heat Bath algorithm of Krzak ala et al The authors of [14] give a Quantum Monte Carlo algorithm for spin systems in a transv erse field. The algorithm w e use, which is describ ed in this section, is a mo dified version of that algorithm. The mo dification whic h we ha ve made is describ ed in subsection 2; everything else in this section constitutes a review of reference [14]. Lik e other w orldline Quan tum Mon te Carlo tec hniques, the algorithm samples the appropriate probability distribution ρ (see section IV) ov er paths in imaginary time via Mark o v Chain Mon te Carlo. Ho wev er this algorithm is only applicable to the case where the Hamiltonian is of the form H = H 0 + V , where H 0 is diagonal in the computational basis | z i and V = − n X i =1 c i σ i x for some set of co efficien ts { c i } whic h are all p ositiv e. F or a Hamiltonian of this form, the distribution ρ o ver paths (from equation 15) is given by ρ ( P ) = 1 Z ( β ) (Π m r =1 c i r ) dt 1 ...dt m e − R β t =0 H 0 ( z ( t )) dt where in this case a path is sp ecified b y an n bit string z 1 (call this the starting state) at time t = 0 and a sequence of flips which o ccur in bits lab eled i 1 , ..., i m at times t 1 , ..., t m (whic h are ordered), where each i r ∈ { 1 , ..., n } and t r ∈ [0 , β ] for r ∈ { 1 , ..., m } . Another w ay of sp ecifying a path is to sp ecify the path P j of each spin j ∈ { 1 , ..., n } . So a path P of the n spin system can b e written P = ( P 1 , P 2 , ..., P n ) . F or each j ∈ { 1 , .., n } , P j sp ecifies the j th bit of the starting state z 1 as w ell as the times at whic h bit flips o ccur in bit j. Note that we only need to consider paths which flip each bit an even num b er of times, since only these paths o ccur with nonzero probability . An example of a path for a system with 2 spins is giv en in figure 20. 25 Figure 20: An example of a path P for 2 spins. The red dots indicate bit flips. In this example the starting state z 1 = 10 and there are four flips in the first spin and 2 in the second (i.e i 1 , i 2 , i 4 , i 6 = 1 and i 3 , i 5 = 2 ). The times of these flips are lab eled t 1 through t 6 . W e now define in detail the Marko v Chain which has limiting distribution ρ. It will b e useful to define p j ( P j | P 1 , ..., P j − 1 , P j +1 , ..., P n ) (for each j ∈ { 1 , ..., n } ) to b e the conditional probability distribution of the path of the j th bit, conditioned on the remainder of the path b eing fixed. As in [14], the up date rule for the Mon te Carlo algorithm consists of the following 3 steps: 1. Randomly and uniformly choose a spin j (where j ∈ { 1 , ..., n } ). 2. Remov e all flips in the path which o ccur in bit j . Also remov e the j th bit from the starting state z 1 . This step corresp onds to wholly remo ving the path P j of spin j . 3. Draw a new path P j (starting v alue and flip times) for spin j from the conditional distribution p j ( P j | P 1 , ..., P j − 1 , P j +1 , ..., P n ) . Of course, the nontrivial part of this algorithm is in sp ecifying a pro cedure which executes step (3) in the ab ov e. T o do this, the authors of reference [14] note that for any index j ∈ { 1 , ..., n } , the diagonal part of the Hamiltonian can b e written as H 0 = g j + f j σ j z (20) where g j and f j are op erator v alued functions of all { σ k z } except for σ j z . F or a given path P , we define F j ( z ( t )) = h z ( t ) | f j | z ( t ) i for t ∈ [0 , β ] . The function F j ( z ( t ) ) is then piecewise constant, and 26 for a giv en path can b e written as F j ( z ( t )) =                            h z 1 | f j | z 1 i , 0 ≤ t < t 1 h z 2 | f j | z 2 i , t 1 ≤ t < t 2 . . . h z m | f j | z m i , t m − 1 ≤ t < t m h z 1 | f j | z 1 i , t m ≤ t ≤ β . Although the ab o ve expression in volv es all the bit flips in the path, the function F j ( z ( t )) can actually only change v alue at times where bit flips o ccur in bits other than the j th bit, since the op erator f j do es not in v olve σ j z . Let us then write F j ( z ( t )) =                      h 0 , 0 = ˜ t 0 ≤ t ≤ ˜ t 1 h 1 , ˜ t 1 ≤ t ≤ ˜ t 2 . . . h q , ˜ t q ≤ t ≤ β = ˜ t q +1 . In this expression the times ˜ t s corresp ond to times at which bit flips o ccur in bits other than bit j (also, h q = h 0 ). With this notation, the pro cedure of reference [14] that generates a new path for bit j consists of the follo wing: 1. Compute the v alue of F j ( z ( t )) as a function of imaginary time along the path. 2. In this step we generate b oundary conditions for the path of bit j at times ˜ t 0 , ˜ t 1 , ... ˜ t q . In other words w e c ho ose q + 1 v alues s 0 , ..., s q ∈ { 0 , 1 } suc h that at time ˜ t r the bit j will b e set to the v alue s r in the new path that w e are generating. T o do this, w e sample the v alues s 0 , ..., s q ∈ { 0 , 1 } for spin j at the times ˜ t 0 , ..., ˜ t q from their joint distribution, which is giv en b y Z ( s 0 , s 1 , ..., s q |{ h 0 , ..., h q } , { ˜ t 1 , ..., ˜ t q } ) = h s 0 | A q | s q ih s q | A q − 1 | s q − 1 i ... h s 1 | A 0 | s 0 i T r [ A q A q − 1 ...A 0 ] where A i = e − λ i [ h i σ j z − c j σ j x ] and λ i = ˜ t i +1 − ˜ t i . W e also define s q +1 = s 0 . 27 3. Having chosen b oundary conditions, w e now generate subpaths for bit j on each interv al [ ˜ t i , ˜ t i +1 ] of length λ i . Suc h a subpath is sp ecified by a num b er of flips w and the time offsets τ 1 , τ 2 , ..., τ w ∈ [0 , λ i ] at whic h flips o ccur (note that the starting v alue of the bit is determined b y the b oundary conditions). The n umber of flips is restricte d to be either ev en or odd dep ending on the b oundary conditions that were chosen for this interv al in the previous step. In this step, the subpath for eac h interv al [ ˜ t i , ˜ t i +1 ] is dra wn from the distribution g i ( τ , ..., τ w ) = 1 h s i +1 | A i | s i i c j w e − s i h i [( τ 1 − 0) − ( τ 2 − τ 1 )+( τ 3 − τ 2 ) − ... +( λ i − τ w )] dτ 1 ...dτ w (21) 4. Put all the subpaths together to form a new path P j for bit j on [0 , β ] . In section 1 we review the metho d outlined in [14] for sampling from the distribution Z ( s 0 , s 1 , ..., s q | ( h 0 , ..., h q ) , ( ˜ t 1 , ..., ˜ t q )) in step (2) of the ab ov e. In section 2 we outline our metho d for sampling from the distribution g i ( τ 1 , ..., τ w ) , whic h differs from the metho d suggested in reference [14]. 1. Generating Boundary Conditions for A Single Spin P ath The prescription outlined in [14] for generating a set s 0 , s 1 , ..., s q from the distribution Z ( s 0 , s 1 , ..., s q | ( h 0 , ..., h q ) , ( ˜ t 1 , ..., ˜ t q )) = h s 0 | A q | s q ih s q | A q − 1 | s q − 1 i ... h s 1 | A 0 | s 0 i T r [ A q A q − 1 ...A 0 ] is as follo ws. First generate s 0 ∈ { 0 , 1 } according to the distribution p ( s 0 ) = 1 T r [ A q A q − 1 ...A 0 ] h s 0 | A q A q − 1 ...A 0 | s 0 i . (computing these probabilities in volv es m ultiplying q tw o b y tw o matrices). Then generate s 1 ∈ { 0 , 1 } according to p ( s 1 | s 0 ) = 1 h s 0 | A q A q − 1 ...A 1 A 0 | s 0 i h s 0 | A q A q − 1 ...A 1 | s 1 ih s 1 | A 0 | s 0 i . Then generate s 2 ∈ { 0 , 1 } from the distribution p ( s 2 | s 1 , s 0 ) = 1 h s 0 | A q A q − 1 ...A 1 | s 1 i h s 0 | A q A q − 1 ...A 2 | s 2 ih s 2 | A 1 | s 1 i and so on. Note that this generates the correct distribution since p ( s 0 ) p ( s 1 | s 0 ) p ( s 2 | s 1 , s 0 ) ...p ( s q | s q − 1 , ..., s 0 ) = Z ( s 0 , s 1 , ..., s q | ( h 0 , ..., h q ) , ( ˜ t 1 , ..., ˜ t q )) . 28 2. Algorithm for Sampling from the Single Spin P ath In tegral with Fixed Boundary Con- ditions W e now present an algorithm whic h samples from the normalized probabilit y distribution ov er single spin paths S ( t ) for t ∈ [0 , λ ] ( S ( t ) takes v alues in { 0 , 1 } ). Here we only present the case of paths with b oundary conditions S (0) = 0 and S ( λ ) = 1 . The other three cases are completely analogous. W e can parameterize suc h a path with fixed b oundary conditions by the times { τ 1 , ..., τ w } at which S ( t ) changes v alue. Note that the n um b er w of suc h flips is o dd due to our choice of b oundary conditions. The distribution o v er paths that we aim to sample from is given by c w e − h [( τ 1 − 0) − ( τ 2 − τ 1 )+( τ 3 − τ 2 ) − ... +( λ − τ w )] dτ 1 ...dτ w h 1 | e − λ [ hσ z − cσ x ] | 0 i . (22) It will also b e useful for us to mak e the change of v ariables from the times ( τ 1 , ..., τ w ) to the w aiting times ( u 1 , ..., u w ) defined b y u 1 = τ 1 u j = τ j − τ j − 1 j ≥ 2 Then the w eight assigned to each path is c w e − h [ u 1 − u 2 + u 3 − ... + λ − P w k =1 u k ] du 1 ...du w h 1 | e − λ [ hσ z − cσ x ] | 0 i whic h we can also write as c w e − R λ t =0 h (1 − 2 S ( t )) dt du 1 ...du w h 1 | e − λ [ hσ z − cσ x ] | 0 i . The algorithm is as follo ws: 1. Start at t=0 in state S (0) defined b y the b oundary conditions. Define B 1 = 1 − 2 S (0) . Set i=1. 2. Draw the waiting time u i un til the next flip from the distribution f ( u i ) = [ p h 2 + c 2 + B i h ] e − u i [ √ h 2 + c 2 + B i h ] If P i j =1 u j > λ then go to step 3. Otherwise define B i +1 = − B i and set i → i + 1 and rep eat step 2. 29 3. T ake the path you hav e generated (which will in general b e longer than λ ), and lo ok at the segmen t [0 , λ ] . If this path satisfies the b oundary condition at t = λ then take this to b e the generated path. Otherwise, thro w aw ay the path and rep eat from step (1). W e no w show that this algorithm generates paths from the distribution (22). Before conditioning on the b oundary conditions b eing satisfied, the probability of generating a sequence of waiting times in ( u 1 , u 1 + du 1 ) , ( u 2 , u 2 + du 2 ) , ... ( u w , u w + du w ) follo wed by any waiting time u w +1 suc h that u w +1 > λ − P w j =1 u i is giv en by f ( u 1 ) f ( u 2 ) ...f ( u w ) du 1 du 2 ...du w Prob  u w +1 > λ − w X j =1 u i  = f ( u 1 ) f ( u 2 ) ...f ( u w ) du 1 du 2 ...du w e − ( λ − P w j =1 u i )[ √ h 2 + c 2 + B w +1 h ] =  w Y i =1 [ p h 2 + c 2 + B i h ]  e − P w i =1 u i √ h 2 + c 2 e − P w i =1 u i B i h du 1 ...du w e − ( λ − P w j =1 u j )[ √ h 2 + c 2 + S B w +1 h ]  =        c w e − λ √ c 2 + h 2 e − R λ t =0 (1 − 2 S ( t )) h du 1 ...du w , if w is even c w  q 1 + ( h c ) 2 + B 1  h c   e − λ √ c 2 + h 2 e − R λ t =0 (1 − 2 S ( t )) h du 1 ...du w , if w is o dd . In the last line w e ha ve used the difference of squares form ula to simplify consecutiv e terms: [ √ c 2 + h 2 + h ][ √ c 2 + h 2 − h ] = c 2 . When we condition on fixed b oundary conditions (whatev er they ma y b e), this generates the correct distribution ov er paths. 30 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.5 2 2.5 3 3.5 4 Energy of a 16 Bit Instance with 2 Planted Solutions s Energy Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Exact Ground State Energy Exact First Excited State Energy Figure 5: The discontin uit y in the circle data that o ccurs near s = 0 . 4 is a Monte Carlo effect that we understand. As can b e seen from the exact numerical diagonalization there is no true discontin uit y in either the ground state energy or the first excited state energy . F or s greater than 0 . 4 the Monte Carlo sim ulation is in a metastable equilibrium that corresp onds to the first excited state. The true ground state energy at eac h s is alwa ys the low er of the circle and the cross at that v alue of s. 31 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 2 4 6 8 10 12 14 16 18 20 Hamming Weight of a 16 Bit Instance with 2 Planted Solutions s Hamming Weight Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Exact Ground State Hamming Weight Exact First Excited State Hamming Weight Figure 6: T ogether with figure 5, we see that the discontin uity in the circles appears in data for both the energy and the Hamming weigh t. This is not indicative of a phase transition in the physical system (as evidenced by the smo oth curves computed by exact numerical diagonalization), but is purely a result of the w ay in which we use the Monte Carlo metho d. 32 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1.5 2 2.5 3 3.5 4 Energy of a 16 Bit Instance After Adding the Penalty Clause s Energy Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Exact Ground State Energy Exact First Excited State Energy 0.422 0.4225 0.423 0.4235 3.59 3.591 3.592 3.593 3.594 3.595 Figure 7: After adding the p enalty clause we see that the energy lev els ha ve an av oided crossing at s ≈ 0 . 423 . The inset sho ws exact numerical diagonalization near the av oided crossing where w e can resolve a tiny gap. The ground state energy is w ell approximated by the low er of the circle and cross at eac h v alue of s . 33 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 2 4 6 8 10 12 14 16 18 20 Hamming Weight of a 16 bit Instance After Adding the Penalty Clause s Hamming Weight Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Exact Ground State Hamming Weight Exact First Excited State Hamming Weight Figure 8: In this figure there is a phase transition which occurs near s ≈ 0 . 423 . W e see from the exact n umerical diagonalization that the Hamming weigh ts of the first excited state and ground state undergo abrupt transitions at the p oint where there is a tin y av oided crossing in figure 7. There is also a jump in the Mon te Carlo data plotted with circles that o ccurs before the av oided crossing: this is a Monte Carlo effect as discussed earlier and has no physical significance. If we lo ok at the Monte Carlo data in figure 7 we conclude that b elow s ≈ 0 . 423 the crosses represent the ground state and after this p oint the circles represen t the ground state. F rom the curren t figure, along with figure 7, we conclude that the Hamming w eight of the ground state jumps abruptly at s ? ≈ 0 . 423 . 34 −2 −1.5 −1 −0.5 0 0.5 1 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Histogram of Curvature Differences (16 spins) Number of Samples ˜ e ( 2 ) U − ˜ e ( 2 ) L Figure 9: The histogram of ˜ e (2) U − ˜ e (2) L for 1 million differen t choices of co efficients c i sho ws a substantial tail for which ˜ e (2) U − ˜ e (2) L > 0 . These sets of co efficien ts corresp ond to b eginning Hamiltonians ˜ H B for which we exp ect the small gap in figure 7 to be remov ed. 35 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1 1.5 2 2.5 3 3.5 Energy with Choice of HB that Removes a Cross (16 spins) s Energy Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Exact Ground State Energy Exact First Excited State Energy Figure 10: A random b eginning Hamiltonian remov es the crossing seen in figure 7. The problem Hamiltonian is the same as that in figure 7. The circles are alw ays b elo w (or equal to) the crosses for all v alu es of s . This means that the circles track the ground state for all s and we see no sign of a small gap in the Monte Carlo data. This is consistent with the display ed exact numerical diagonalization. 36 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 2 4 6 8 10 12 14 16 18 20 Hamming Weight with Choice of HB that Removes a Cross (16 spins) s Hamming Weight Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Exact Ground State Hamming Weight Exact First Excited State Hamming Weight Figure 11: F rom figure 10 w e see that the ground state of the Hamiltonian corresp onds to the circles for all v alues of s and the Hamming weigh t of the circles here go es smo othly to the Hamming w eight of the unique satisfying assignmen t. (The jump in the Monte Carlo data corresp onding to the crosses is due to the Mon te Carlo effect discussed earlier.) 37 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 10 15 20 25 30 35 40 45 Energy of a 150 Bit Instance with 2 Planted Solutions s Energy Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Second Order Perturbation Theory Second Order Perturbation Theory 0.4 0.5 0.6 0.7 0.8 0 0.1 0.2 0.3 0.4 0.5 ∆ E s Monte Carlo Data Second Order Perturbation Theory Figure 12: The crosses, which represent the Monte Carlo data seeded with 111 ... 1 , are b elow (or equal to) the circle data for all v alues of s . (This is seen more clearly in the inset which shows the p ositive difference b et ween the circle and cross v alues). W e conclude that the crosses track the ground state energy which is smo othly v arying. The jump in the circle data is a Monte Carlo effect and for s ab ov e 0 . 2 the circles track the first excited state. 38 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 20 40 60 80 100 120 140 160 180 Hamming Weight of a 150 Bit Instance with 2 Planted Solutions s Hamming Weight Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Figure 13: Comparing with figure 12 w e see that the Hamming weigh t for the first excited state and the ground state are continous functions of s . W e only obtain data for the first excited state for v alues of s larger than s ≈ 0 . 2 . 39 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 10 15 20 25 30 35 40 45 Energy of a 150 Bit Instance After Adding the Penalty Clause s Energy Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Second Order Perturbation Theory Second Order Perturbation Theory Figure 14: Adding the p enalty clause makes the cross data go ab o ve the circle data at s ≈ 0 . 49 . This is sho wn in more detail in figure 15 where we plot the energy difference b etw een the first tw o levels as a function of s . W e interpret this to mean that the Hamiltonian H ( s ) has a tiny gap at s ? ≈ 0 . 49 . 40 0.4 0.45 0.5 0.55 0.6 0.65 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 Energy Difference between the Two Lowest Levels (150 spins) s Energy Difference Monte Carlo Data Second Order Perturbation Theory Figure 15: The energy difference, circles min us crosses, from figure 14 near the v alue of s where the difference is 0 . Note that second order p erturbation theory do es quite well in predicting where the difference go es through zero. 41 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 20 40 60 80 100 120 140 160 180 Hamming Weight of a 150 Bit Instance After Adding the Penalty Clause s Hamming Weight Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Figure 16: Looking at figure 14 we see that the ground state is represented by the crosses to the left of s ≈ 0 . 49 and is represen ted by the circles after this v alue of s . T racking the Hamming weigh t of the ground state, w e conclude that it changes abruptly at s ≈ 0 . 49 . 42 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 0 50 100 150 200 250 300 350 Histogram of Curvature Differences (150 spins) Number of Samples ˜ e ( 2 ) U − ˜ e ( 2 ) L Figure 17: Histogram of ˜ e (2) U − ˜ e (2) L for 100000 c hoices of co efficien ts c i for our 150 spin instance. Note that a go od fraction hav e ˜ e (2) U − ˜ e (2) L > 0 . 43 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 10 15 20 25 30 35 40 Energy with Choice of HB that Removes a Cross (150 spins) s Energy Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Second Order Perturbation Theory Second Order Perturbation Theory 0.4 0.5 0.6 0.7 0.8 0.4 0.5 0.6 0.7 0.8 ∆ E s Monte Carlo Data Second Order Perturbation Theory Figure 18: A random choice of co efficients such that ˜ e (2) U − ˜ e (2) L > 1 2 giv es rise to an H ( s ) where there is no longer an av oided crossing. The circles here corresp ond to the ground state for all s since the cross data is alw ays ab o ve (or equal to) the circle data for all s. This can b e seen in the inset where we hav e plotted the energy difference, crosses minus circles. The crosses hav e a Monte Carlo discontin uity near s ≈ 0 . 2 , after whic h they corresp ond to the first excited state. 44 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 20 40 60 80 100 120 140 160 180 Hamming Weight with Choice of HB that Removes a Cross (150 spins) s Hamming Weight Monte Carlo Seeded with 000...0 Monte Carlo Seeded with 111...1 Figure 19: Lo oking at figure 18 we see that the ground state corresp onds to the circles for all v alues of s so we see here that the Hamming weigh t of the ground state go es smo othly to its final v alue as s is increased. W e tak e this as further evidence that this choice of ˜ H B w ould corresp ond to success for the quantum adiabatic algorithm for this instance. 45