University of Southern California Dissertations and Theses

Understanding physical quantum annealers: an investigation of novel schedules

(USC Thesis Other)

Understanding physical quantum annealers: an investigation of novel schedules

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content UNDERSTANDING PHYSICAL QUANTUM ANNEALERS: AN INVESTIGATION OF NOVEL SCHEDULES by Zoe Gonzalez Izquierdo A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) May 2020 Copyright 2020 Zoe Gonzalez Izquierdo Contents List of Tables iv List of Figures v 1 Introduction 1 1.1 Quantum annealing . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 The D-Wave 2000Q quantum annealing processor . . . . . . . . . . 3 1.3 Annealing schedules . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3.1 Forward annealing . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.2 Reverse annealing . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 Testing a quantum annealer as a quantum thermal sampler 11 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Experimental methods . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 The case of pausing before the minimum gap . . . . . . . . . 16 2.3.2 The case of pausing after the minimum gap . . . . . . . . . 18 2.3.3 Quantifying the impact of the quench . . . . . . . . . . . . . 20 2.3.4 Other observables: magnetization . . . . . . . . . . . . . . . 22 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3 Discriminatingnon-isomorphicgraphswithanexperimentalquan- tum annealer 28 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 The graph isomorphism problem . . . . . . . . . . . . . . . . . . . . 29 3.2.1 Graph invariants . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3.1 Choice of graphs . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4.1 G13 and G13p . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4.2 G27 and G27p . . . . . . . . . . . . . . . . . . . . . . . . . . 40 ii 3.4.3 G17 and G17p . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.4.4 Set of four G25 graphs . . . . . . . . . . . . . . . . . . . . . 42 3.4.5 Set of four G33 graphs . . . . . . . . . . . . . . . . . . . . . 43 3.5 Conclusions and outlook . . . . . . . . . . . . . . . . . . . . . . . . 49 4 Ferromagnetically shifting the power of pausing 52 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.2 Problems, quantum annealing and setup . . . . . . . . . . . . . . . 55 4.2.1 Problem definition . . . . . . . . . . . . . . . . . . . . . . . 55 4.2.2 Solving on a quantum annealer . . . . . . . . . . . . . . . . 56 4.2.3 Hardware,annealingparameters,andprobleminstancespecifics 58 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.3.1 Annealing without pause, effect of|J F | . . . . . . . . . . . . 61 4.3.2 Improvement with pause . . . . . . . . . . . . . . . . . . . . 61 4.3.3 Optimal pause location shift with|J F | . . . . . . . . . . . . 64 4.3.4 Help of partial gauges . . . . . . . . . . . . . . . . . . . . . 65 4.3.5 Picture for how pausing and|J F | affect p success . . . . . . . . 67 4.4 Conclusions and outlook . . . . . . . . . . . . . . . . . . . . . . . . 78 Reference List 81 Appendices 93 A Testing a quantum annealer as a quantum thermal sampler 94 A.1 Reverse annealing protocol . . . . . . . . . . . . . . . . . . . . . . . 94 A.2 Frustrated chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 A.3 Effect of temperature . . . . . . . . . . . . . . . . . . . . . . . . . . 97 B Discriminatingnon-isomorphicgraphswithanexperimentalquan- tum annealer 99 B.1 Graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 C Ferromagnetically shifting the power of pausing 106 C.1 Level-based mapping for BD-MST problems . . . . . . . . . . . . . 106 C.2 BD-MST problem instances . . . . . . . . . . . . . . . . . . . . . . 111 C.3 Embedding statistics . . . . . . . . . . . . . . . . . . . . . . . . . . 113 iii List of Tables C.1 n=5 graphs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 C.2 Graph weights used to generate different instances. . . . . . . . . . 112 C.3 Mapped problem size for n = 4− 6 using the current Chimera architecture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 C.4 Mapped problem size forn = 4− 10 using the future Pegasus archi- tecture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 iv List of Figures 1.1 Chimera graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Energy scales of driver and problem Hamiltonian . . . . . . . . . . 5 1.3 Diagram of forward annealing schedule . . . . . . . . . . . . . . . . 8 1.4 Diagram of reverse annealing schedule . . . . . . . . . . . . . . . . 9 2.1 Comparison ofhH IM i between DW and QMC for n = 138 ands p = 0.2 and relative difference as a function of size . . . . . . . . . . . . 17 2.2 Comparison ofhH IM i between DW and QMC for n = 138 ands p = 0.5− 0.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3 Comparison ofhH IM i between DW and ME . . . . . . . . . . . . . 21 2.4 Comparison ofhM 2 z i andhH IM i between DW and QMC through full anneal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5 Relative difference inhM 2 z i between DW and QMC as a function of size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 Exact difference inhH p i throughout the anneal between G13 and G13p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2 Experimental results for G13 and G13p . . . . . . . . . . . . . . . . 44 3.3 Experimental results for G27 and G27p . . . . . . . . . . . . . . . . 45 3.4 Experimental results for G17 and G17p . . . . . . . . . . . . . . . . 46 v 3.5 Experimental results for set of G25 graphs . . . . . . . . . . . . . . 47 3.6 Experimental results for set of G33 graphs . . . . . . . . . . . . . . 48 4.1 p success and TTS as functions of|J F | . . . . . . . . . . . . . . . . . 62 4.2 Shift of optimal pause location with|J F | for a demo instance . . . . 63 4.3 Effect of pause duration on p success and TTS . . . . . . . . . . . . . 73 4.4 Improvement of TTS with a short pause . . . . . . . . . . . . . . . 74 4.5 Shift of optimal pause location with|J F | for ensemble . . . . . . . . 75 4.6 Improvement in p success with partial gauges . . . . . . . . . . . . . . 75 4.7 Diagram of different regimes in the anneal . . . . . . . . . . . . . . 76 4.8 Shift of minimal gap with|J F | for toy problem . . . . . . . . . . . . 77 4.9 Shift of energy levels with|J F | for toy problem . . . . . . . . . . . . 77 A.1 Comparison ofhH IM i between DW and QMC for n = 138 ands p = 0.2 for reverse annealing . . . . . . . . . . . . . . . . . . . . . . . . 94 A.2 Comparison ofhH IM i between DW and QMC for n = 138 ands p = 0.5− 0.8 for reverse annealing . . . . . . . . . . . . . . . . . . . . . 95 A.3 Comparison ofhH IM i between DW and QMC for n = 138 ands p = 0.2 for a frustrated chain . . . . . . . . . . . . . . . . . . . . . . . . 97 A.4 Effect of temperature onhH IM i . . . . . . . . . . . . . . . . . . . . 98 B.1 G13 and G13p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 B.2 G17 and G17p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 B.3 G25p1, G25p2, G25p3 and G25p4 . . . . . . . . . . . . . . . . . . . 101 B.4 G27 and G27p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 B.5 G33, G33p, G33p1, G33p2 . . . . . . . . . . . . . . . . . . . . . . . 104 C.1 Comparison of embedded problem size between Chimera and Pegasus115 vi Chapter 1 Introduction 1.1 Quantum annealing Quantum annealing [1–8] is a metaheuristic designed for solving discrete com- binatorial optimization problems. When the solution to a problem of interest can be encoded as the ground state (GS) of a problem Hamiltonian H p , and initial- izing the system in the known GS of a different Hamiltonian (the driver, H d ), we can make use of the adiabatic theorem of quantum mechanics [9] to provide a guarantee that the GS can be reached with high probability for a sufficiently slow interpolation, in which the instantaneous time-dependent HamiltonianH(t) starts off at t = 0 as H d , and by the end of the evolution at t =t a (the annealing time), it has fully morphed into H p . The adiabatic theorem posits that a system will remain in an instantaneous eigenstate if change in the system happens slowly enough. This condition on the rate of change is dependent on the spectral gap Δ, i.e. the energy differ- ence between the GS and the first excited state (FES) of the instantaneous time- dependent Hamiltonian H(t): t a max 0≤t≤ta |hFES(t)| dH(t) dt |GS(t)i| Δ 2 (t) . (1.1) Adiabaticity is typically defined for zero temperature (T = 0) and a pure state, but can be extended for a closed system at finite T [10], as well as for an open 1 system [11], in which case the mixed state remains close to a Gibbs thermal state described by a density matrix ρ G = e −βH Z , (1.2) where β is the inverse temperature β = 1/k B T with k B the Boltzmann constant, and Z the partition function. An adiabatic condition when the gap Δ is 0 also exists [12]. Quantum annealing was first conceived purely as an algorithm [1] to be run on classical computers, a quantum counterpart to simulated annealing [13]. In simu- lated annealing, by starting the evolution of the system at a high temperature and slowly decreasing it to 0, one aims to find the global minimum of a complicated energy landscape. Analogously, in quantum annealing, quantum fluctuations are high at the beginning of the evolution, favoring exploration of the full energy land- scape, and decrease as time goes on, by which point the system is expected to have found the state with minimum energy. Quantum tunneling can help getting through large energy barriers that might otherwise prevent the system from leav- ing a local minimum, which is especially useful for energy landscapes with many valleys. But a decade ago, experimental quantum annealers [14–16] irrupted into the quantum computing market, bringing along great promise—and just as much skep- ticism [17, 18]. Whether these machines provide any advantages over their classical analogues — or whether they could even be considered genuine quantum proces- sors [19–23] — is still a topic of much debate [24–27]. While they have so far failed to outperform state-of-the-art classical methods at optimization tasks [28, 29] potential uses as Boltzmann samplers [30, 31] and as effective quantum simulators [32–35] have been gaining traction. 2 In recent years, some of the technological improvements have led to the pos- sibility of modifying the annealing schedules, with new features such as pausing and reverse annealing [36–40]. We will explore potential uses for these novel capa- bilities, and the physics that underlies the behavior we observe when they are deployed. 1.2 The D-Wave 2000Q quantum annealing pro- cessor Our experimental results are generated by a D-Wave 2000Q quantum anneal- ing processor (DW), located at NASA Ames Research Center. It consists of a lattice of niobium superconducting quantum interference device (SQUID) operat- ing at a temperature of approximately 12mK. The two lower energy levels of these superconducting loops act as the qubits (quantum bits) of the processor. Just like classical bits, they can be in either of two states (0 and 1, or equivalently -1 and +1, corresponding to the two energy levels) but they can also be in a superposition of the two states. These qubits are connected according to a Chimera architecture [15, 41]. Chimera graphs are made up of smaller unit cells, which can be repeated in two dimensions to achieve graphs of different sizes. Each cell is a complete bipartite graphK 4,4 , where a qubit is connected to four others within its cell and two more in adjacent cells. The Chimera graph C L is obtained by arranging the unit cells in a square pattern with L cells per side. The processor we are using features a C 16 , with a total of 256 cells and 2048 qubits (although a few are inoperative due to fabrication issues). Its complete connectivity graph is shown in Fig. 1.1. 3 Figure 1.1: The quantum annealer hardware adjacency graph, a 16× 16 Chimera graph. The nonworking qubits are not shown here. DW implements a time-dependent, transverse field Ising Hamiltonian: H(s) =A(s)H d +B(s)H p , (1.3) wheres∈ [0, 1] is a rescaled, dimensionless time parameter andA(s) andB(s) are functions of s such that A(0)B(0)' 0 and B(1)A(1)' 0. These functions are fixed by the hardware, as shown in Fig. 1.2. H d is a transverse-field driver Hamiltonian: H d = n X i=1 σ x i , (1.4) 4 0.0 0.2 0.4 0.6 0.8 1.0 s 0 1 2 3 4 5 6 Energy (GHz) A(s) B(s) Figure 1.2: Strengths of the driver and problem Hamiltonians as functions of the dimensionless time parameter s, shown in units of h = 1. whose GS is an equal superposition of all possible states, |GS(H d )i = 2 n X i=1 |ii, (1.5) and H p an Ising type problem Hamiltonian: H p = X hi,ji J ij σ z i σ z j + n X i=1 h i σ z i . (1.6) hi,ji denotes the qubitsi andj being linked together using a coupler, implemented by a tunable inductance. The coupler can be assigned a strength J ij ∈ [−2, 1]. A negative J ij encourages the qubits to be aligned, while with a positive J ij they would tend to point in opposite directions. There is also an external longitudinal magnetic field acting on each qubit, h i ∈ [−1, 1]. The corresponding Ising energy can be expressed as: E(~ s) = X hi,ji J ij s i s j + n X i=1 h i s i , (1.7) 5 where~ s =|s 1 ,...s n irepresentsthestateofthequbits, withs i ∈{−1, +1}.Equiva- lently, this can also be expressed as a quadratic unconstrained binary optimization (QUBO) problem, where a cost function is written in terms of binary variables x i ∈{0, 1}: C(~ x) = X i,j Q ij x i x j + n X i=1 Q ii x i . (1.8) We can easily switch between the two formulations by setting s i = 1− 2x i . While both Ising and QUBO formulations only include terms of up to quadratic order, higher order expressions can be incorporated by using ancilla variables to express them as sums of quadratic terms, which results in needing to use additional qubits. Because the Chimera graph is not complete, not all problems that can be encoded by a general Ising Hamiltonian can be directly translated into the annealer’s problem Hamiltonian. We call them native problems when the direct encoding is possible. For non-native problems, minor-embedding [42, 43] can be used if we wish to run them on DW. To embed a problem, several physical qubits are coupled ferromagnetically (J ij < 0), and act as a single logical variable from the original problem. 1.3 Annealing schedules Even though A(s) and B(s) are fixed, recent technological advancements have made it possible to exercise some control over the annealing schedule. This means that instead of having the dimensionless time parameter s be simply related to the time by s(t) = t/t a , we are able to redefine s(t) as a piece-wise function of t. This allows us, for instance, to pause at specific points in the anneal, where s (and thus the Hamiltonian) remains fixed for some period of time. It has been 6 shown that introducing such a pause at an appropriate time (shortly after the minimum gap) can increase the success rate of the quantum annealer when used as an optimizer [38, 40]. This is due to the thermal state of the system at this point starting to have a larger contribution from the GS, and the system being able to thermalize during the pause. We can also use different annealing rates ds/dt for different portions of the anneal. One use of this feature is performing approximate quenches by abruptly turningdownthestrengthofthedriverHamiltonianfromagivenpointmid-anneal. It has also been proposed that, combining a pause with a sufficiently fast quench to the end of the anneal, such that the evolution becomes non-adiabatic, it could be possible to effectively measure the system at the time of the pause, allowing the annealer to be used as a quantum simulator, in the same role of Quantum Monte Carlo methods [32–34]. We will use two qualitatively different types of schedules, distinguished by the way in which the first part of the anneal is performed, which we detail below. 1.3.1 Forward annealing Here, the system is assumed to be initialized in the thermal state of the Hamil- tonian at s(0) = 0, where H(s = 0)≈ H d . At the low temperature that DW operates, this thermal state has very high weight on the GS of H(0). The system is annealed froms = 0 to some intermediates =s p at a rate of 0< ds dt i ≤ 1μs −1 . We then pause at s = s p for some time t p and finally anneal to s = 1 at a rate 0< ds dt f ≤ 1μs −1 . A diagram of this schedule is shown in Fig. 1.3. 7 t 0 s p 1 s t p ds dt i ds dt f Figure 1.3: Diagram of the forward annealing schedule s(t), showing the initial anneal to s p , performed at a rate ds dt i , a pause of length t p and the final quench, performed at a rate ds dt f . The durations of the three parts are not to scale. 1.3.2 Reverse annealing For our ‘reverse’ anneal protocol, the system is initialized in a classical state (each qubit is either in a computational one or a computational zero state) chosen by the user at s(0) = 1 instead of s(0) = 0. We then anneal backwards from s = 1 to some intermediate point s = s p at a rate of 0 < ds dt i ≤ 1μs −1 , pause for some time t p and anneal forward to s = 1 at a rate of 0 < ds dt f ≤ 1μs −1 . A diagram for this protocol is shown in Fig. 1.4. For each subsequent anneal, we can either reinitialize at the chosen classical state, or use the one returned from the previous anneal as our initial state. Presumably doing so can decrease the thermalization time [37], and help with exploitation of regions of interest in the solution space [44, 45]. A key difference between the forward and reverse annealing protocols is that for a sufficiently large s p the reverse annealing protocol avoids crossing the minimum 8 t 0 s p 1 s t p ds dt i ds dt f Figure 1.4: Diagram of the reverse annealing schedule s(t), showing the initial anneal to s p , performed at a rate ds dt i , a pause of length t p and the final quench, performed at a rate ds dt f . The durations of the three parts are not to scale. gap encountered during the forward annealing protocol. This can help reduce transitions to excited states, which are more likely to happen when the spectral gap is smaller. 1.4 Outline Therestofthisthesisisstructuredintothreeparts. InChapter2,weinvestigate the potential of the current hardware to reliably sample from the thermal state associated with a target quantum Hamiltonian, fulfilling the same role as Quantum Monte Carlo methods. This could in principle be accomplished with a pause- and-quench type schedule, which under the appropriate conditions could mimic a measurement in the middle [36, 37]. We find, however, that the quantum annealer fails to produce the correct expectation values predicted by Quantum Monte Carlo. 9 Comparing to master equation simulations, we are able to explain this discrepancy by the limitations in the annealing speed imposed by current hardware. In Chapter 3, and in light of the previous results, we consider an application for which the pause and quench protocol as currently implemented can be useful. Quantum annealing had previously been proposed as a discriminator for graphs whose classical Ising spectra are the same, but could be differentiated by their quantum spectra [46, 47], hence solving certain cases of the Graph Isomorphism problem. We find that using a reverse annealing schedule with a pause and quench we are able to tap into the quantum regime of the anneal and distinguish such sets of graphs. Finally, in Chapter 4, we test pausing on a class of optimization problems that are not native to the annealer’s architecture. An appropriately timed pause had been previously shown to drastically increase the probability of success of native problems [38], but embedding brings along differences that might affect those results. We study the effect of the ferromagnetic coupling between physical qubitsrepresentingthesamephysicalone, anddevelopanapproachtogaugetrans- formations that is compatible with an asymmetric range in the coupling strength. We find that while pausing still helps, its effect is less drastic than for native problems, which our physical understanding supports. 10 Chapter 2 Testing a quantum annealer as a quantum thermal sampler The contents of this chapter have been produced in collaboration with professors Tameem Albash and Itay Hen. 2.1 Introduction Describing the static and dynamical properties of a many-body quantum sys- tem remains a considerable challenge in physics. Strong correlations among the particles of a many-body quantum system can lead to highly complex entangled states, which exist in a vector space whose dimension grows exponentially with the size of the system. Due to this exponential scaling, affecting both the time and memory that a classical computer needs in order to perform the relevant compu- tations, even simulations of systems with only a few correlated particles quickly become intractable. Tackling this problem with a quantum computer or simulator continues to inspire [48] and drive the field of quantum simulation [49–52]. The adiabatic paradigm of quantum computing [53–55] naturally lends itself to tackling the problem of studying ground state (GS) properties of quantum sys- tems. By preparing the system in the GS of a ‘trivial’ Hamiltonian and performing 11 a sufficiently slow interpolation towards the target many-body quantum Hamilto- nian, the adiabatic theorem of quantum mechanics [9, 56, 57] provides a guarantee that the state at the end of the evolution will be close to the target GS. But what if we are interested in finite temperature properties of the same sys- tem? For decades, the work-horse for addressing this question has been Quantum Monte Carlo (QMC) methods [58–60]. Despite extensive work on developing more sophisticated methods [61–64], no universal method whose space and time com- plexity scale polynomially with problem size is yet known. It remains an open research topic how to efficiently encode a generic thermal state into the GS of a quantum Hamiltonian [65, 66]. Alternatively, one can ask whether the open-system dynamics [67] of a system naturally has the thermal state as its steady state, and if so whether the dissipative dynamics can be used in lieu of a true quantum algorithm to prepare such a state. Having the steady state of the dissipative dynamics be the standard Gibbs state associated with the target Hamiltonian is a non-trivial assumption; it is known to be the case of Markovian weak-coupling limit master equations [68, 69] satisfying the Kubo-Martin-Schwinger (KMS) condition [70]. An open-system adiabatic theorem [11, 71–74] provides a guarantee that in the long-time limit the state of the system will be close to the desired thermal state. Recently, two experiments [36, 37] carried out on D-Wave quantum annealing processors were reported to have successfully reproduced certain thermal proper- ties of complex quantum systems. In these studies, the processor is used as an analogue quantum simulator: the superconducting circuit hardware [75–78] imple- ments a transverse-field Ising model (TFIM) Hamiltonian in its low-lying spec- trum. The system is allowed to thermally relax at a particular realization of the 12 TFIM Hamiltonian, and diagonal thermodynamic observables are calculated using measurements performed after quenching the system to a purely Ising Hamilto- nian. In Ref. [37], the emergence of the finite order parameter associated with the Kosterlitz-Thouless phase transition of the TFIM on the triangular lattice [79] was observed, while Ref. [36] observed the finite-size precursors of the phase transi- tions associated with the TFIM on three-dimensional cubic lattices. In both cases, a minor-embedding procedure [42, 43] was required to map the target Hamiltonian connectivity graph onto the hardware physical connectivity. Both studies found good agreement between theory and experimental observations when measured quantities were restricted to the order parameter of the studied system [37] or to critical parameters indicating the location of phase transitions [36]. Motivatedbytheseresults,weexamineinthisstudytheextenttowhichonecan use such hardware to provide thermal samples for quantum Hamiltonians imple- mented by the hardware, especially at ever increasing system sizes. As the simu- lation of quantum systems on quantum information processors gains traction and moves on to problems which are classically difficult (or virtually impossible) to solve, we face the question of how to verify or trust the results produced by these devices. The one-dimensional (1D) TFIM, a canonical system in quantum ther- modynamics, provides a simple case study, which hopefully helps build confidence in the device as more difficult problems are tackled. To that aim, we study the performance of DW on the 1D nonfrustrated TFIM. We verify the accuracy of the expectation values calculated from the experimental samples by comparing them to data obtained from QMC simulations. If such devices are to be used more widely as quantum simulators or as thermal state 13 samplers, we believe this is a crucial first step in the process of validating the reliability of the device for this task. Wefindingeneralthatthethermalexpectationvaluescalculatedusingthesam- plesfromtheannealerarenotingoodagreementwiththeexpectedvaluesobtained by QMC, especially in the quantum paramagnetic region where the transverse field dominates [80]. Master equation simulations [81] on small size problems confirm that a reason for this discrepancy is the inability of current devices to perform measurements at the target Hamiltonian; instead measurements are performed on an Ising Hamiltonian with negligible transverse field, which requires a quench from the target Hamiltonian to the Ising Hamiltonian. Our simulations indicate that increasing the maximum annealing rate by a factor of 10 5 would reproduce the expected thermal expectation values, but this would in turn require annealing times on the order of a few picoseconds. Extensivity of the Hamiltonian suggests that this annealing rate would need to decrease with the inverse of the system size, further suggesting that such a measurement approach is not a scalable approach to address questions on the static properties of thermal quantum states. 2.2 Experimental methods We consider the one-dimensional transverse-field Ising model (1D TFIM) [80] with the n-qubit Hamiltonian H TFIM =−AH TF +BH IM =−A n X i=1 σ x i −B n X i=1 σ z i σ z i+1 , (2.1) where we choose periodic boundary conditions (σ z n+1 = σ z 1 ). The two-fold degen- erate GS of the ferromagnetic Ising chain H IM , corresponding to all-spins up and 14 all-spins down, satisfies all the Ising couplings (it is nonfrustrated) and has GS energy E 0 =−n. Our objective will be to measure the expectation value of H IM for different ratios of A to B, which is ideally given by: hH IM i = Tr h H IM e −βH TFIM i Tr [e −βH TFIM ] . (2.2) Measurements at different rations of A to B cannot directly be performed on the current hardware. We can, however, try to mimic them with an appropriately chosen annealing schedule. To do this, we progress towards the thermal state associated with the Hamiltonian at s = s p by performing anneals up to s = s p and pausing at s p for a period of time thereby allowing the system to thermally relax to its steady state. We then quench as rapidly as allowed by the hardware (ds/dt = 1μs −1 ) towards s = 1, where ideally a computational basis measurement is performed. The success of this method depends on two factors: the ability of the hardware to thermalize to the desired Gibbs state, and the annealing rate during the quench being fast enough to prevent any changes to the state. Werestrictourattentionheretothis‘forward’annealingprotocol. InAppendix A.1, we discuss an alternative protocol using ‘reverse’ annealing; while the details of the two protocols are different 1.3, we find no qualitative differences between the results of the two protocols. Weincorporategaugeaveraging[19]toaverageoutsystematicbiasesthatmight exist on the qubits and/or couplers, such as certain qubits more readily aligning in one direction. Gauge averaging is carried out by repeating each run of n a anneals 100 times, where for each run we apply a transformation of the form J ij →a i a j J ij where a i ∈{−1, +1} is chosen at random. This transformation corresponds to 15 applying a unitary transformation to the Hamiltonian, so the energy spectrum is unchanged but the states are relabeled accordingly. For example, the ungauged classical state (s 1 ,...,s n ) is mapped to the state (a 1 s 1 ,...,a n s n ). The results in different gauges can be readily mapped back to the states of the original problem. 2.3 Results In what follows, we separate the discussion to two cases which we argue exhibit qualitatively different behaviors, according to whether the location of the pause is before or after the minimum gap. For the processor we use, the minimum gap of the TFIM occurs at s ∗ ≈ 0.346, with A(s ∗ )/h≈ 1.05GHz. When the pause takes place before the minimum gap, that is, s p <s ∗ , the driver Hamiltonian still dominates, and if the system is in the instantaneous thermal state of H(s p ), the expectation value ofH IM —the diagonal energy—in this region will be closer to the GS energy of H TF than to that of H IM . On the other hand, with a pause after the minimum gap, s p > s ∗ , we find the system in the region dominated by the problem Hamiltonian, and hence the expectation value of H IM will be closer to its GS energy. 2.3.1 The case of pausing before the minimum gap We first consider the case of s p <s ∗ , i.e., where the pause takes place prior to the system reaching its minimum gap. In this region, the strength of A(s p )H TF is considerably greater than that of B(s p )H IM , so there is little overlap between the GS of H IM and the instantaneous GS of H(s p ). The QMC results reflect this. 16 − 20 − 15 − 10 0 400 800 1200 1600 2000 t p (μ s) − 140 − 135 − 130 hH IM i QMC GS DW 0 100 200 300 400 500 n 6.8 7.0 7.2 7.4 7.6 7.8 |hHIMiQMC−h HIMiDW| |hHIMiQMC| Figure 2.1: Top: Expectation value for H IM for n = 138 and s p = 0.2. Solid blue curve corresponds to QMC prediction, dashed orange curve corresponds to the Ising GS energy, and data points correspond to the experimental results using the forward annealing protocols. Bottom: The relative difference between the experimental and QMC results as a function of problem size for s p = 0.2, t p = 1.9ms and 1500 anneals per gauge. Error bars correspond to the 95% confidence interval calculated using a bootstrap over the 100 gauges. For instance, with pause location s p = 0.2 and problem sizen = 138, the expecta- tion value of H IM ishH IM i exact sp=0.2 ≈−16.5, while the Ising GS energy is−138. When we calculate this same expectation value using the output from the annealer, however, we findhH IM i DW sp=0.2 ≈−132, suggesting that the experimental mea- surement is failing to capture the state of the system in the region before the 17 minimum gap, with the relative difference between experimental and QMC results being 700%. This is depicted in Fig. 2.1 (top). This experimental behavior of being very far from the expected value and very close to the Ising GS value holds for all the problem sizes we studied, ranging from 4 to over 500 qubits. Figure 2.1 (bottom) shows the relative difference between experiment and theory as a function of size. While very small problems seem to do slightly worse, this difference remains fairly constant for a large range of sizes. We note that the experimental results do not change in any significant way as we increase the pause timet p , so we deduce that the system is very close to its steady state at the pause point. A possible explanation for this dramatic difference is the upper limit on the annealing rate, which constrains our quench to a maximum ds dt f = 1μs −1 . This rate is likely not sufficiently fast to prevent the state from evolving during the quench and allowing the dynamics to populate the Ising GS. We validate this conjecture using master equation simulations in Sec. 2.3.3. 2.3.2 The case of pausing after the minimum gap When we choose the pause to occur after the minimum gap, i.e., in the case where s p >s ∗ , the Hamiltonian B(s p )H Ising dominates over A(s p )H TF . We there- fore expecthH IM i to be significantly closer to E 0 =−n. The experimental results in this region turn out to depend sensitively on the value of s p , as we show in Fig. 2.2. We observe the best agreement with QMC around s p = 0.6. When we pause at later points, however,hH IM i obtained from the annealer moves away from the correct value. 18 0 400 800 1200 1600 2000 t p (μ s) − 138.0 − 137.5 − 137.0 − 136.5 − 136.0 − 135.5 hH IM i QMC GS DW (a) 0 400 800 1200 1600 2000 t p (μ s) − 138.0 − 137.5 − 137.0 − 136.5 − 136.0 − 135.5 hH IM i QMC GS DW (b) 0 400 800 1200 1600 2000 t p (μ s) − 138.0 − 137.5 − 137.0 − 136.5 − 136.0 − 135.5 hH IM i QMC GS DW (c) 0 400 800 1200 1600 2000 t p (μ s) − 138.0 − 137.5 − 137.0 − 136.5 − 136.0 − 135.5 hH IM i QMC GS DW (d) Figure 2.2: hH IM i comparison between experiment and simulations for n = 138 and s p = 0.5 (a), 0.6 (b), 0.7 (c) and 0.8 (d). Solid blue curve corresponds to QMC result, dashed orange curve corresponds to the Ising GS energy, and data points correspond to the DW results using the forward annealing protocols. Error bars correspond to the 95% confidence interval calculated using a bootstrap over the 100 gauges. The dependence on s p is likely again due to the quench rate. For s p < 0.6, the thermal state ats p still does not have a complete overlap with the Ising GS, so the quench from s = s p to s = 1 is sufficiently slow for the system to repopulate the Ising GS. Hence we find that the experimentalhH IM i is lower than predicted by QMC. Fors p > 0.6, the HamiltonianA(s p )H TF is very weak relative toB(s p )H IM , so the dynamics are expected to be extremely slow and effectively frozen [38, 82]. 19 In this case, the system likely does not have enough time to thermalize. The points arounds p = 0.6 represent a ‘sweet spot’ where the system is still able to thermalize and is only minimally affected by the slowness of the quench. 2.3.3 Quantifying the impact of the quench In the previous section we have seen how the experimental samples fail to reproduce the correct thermal expectation values associated with the point where we pause s p . We argued that this is likely due to the quench rate from s =s p to s = 1 not being sufficiently fast, with the system continuing to evolve during the quench. In order to confirm this hypothesis, we use the adiabatic master equation (ME) [81] to simulate the quantum annealing process using different quench rates. The key feature of this numerical method is that for any fixed s Hamiltonian, the fixed point of the dissipative dynamics is the Gibbs state of H(s). Therefore, for any sufficiently long pause and sufficiently fast quench, we expect the simulation results to agree with the theoritical prediction. Because of the computational cost of the simulations, we are only able to obtain exact results for the smallest systems with n = 4, but already at these sizes we are able to see the adverse effects of the quench rate for reproducing the correct thermal expectation values. We first run these simulations using the same schedule as the annealer, with a quench annealing rate of ds dt = 1μs −1 . These are shown as the ‘slow quench’ results in Fig. 2.3, where we find the simulatedhH IM i results are in agreement with experiment at all pause times s p during the anneal. Specifically, we findhH IM i to be close to the Ising GS energy regardless of the value of s p . In fact, for this small system size, even if we considered closed system dynamics, i.e. we decoupled the system from the thermal environment, we would get the same result, pointing to 20 the fact that the quench is so slow that the evolution across the minimum gap is effectively adiabatic. Next we survey a wide range of quench rates to find how fast the quench must be to reproduce the correct thermal expectation values. For the n = 4 system, we find that we need to go up to ds dt = 10 3 μs −1 before seeing any change in the behavior ofhH IM i. As the annealing rate ds dt becomes larger,hH IM i gets closer to the exact results, and they finally become in good agreement when the annealing rate is ds dt = 10 5 μs −1 (the ‘fast quench’ in Fig. 2.3). This implies that we need to increase the current fastest rate allowed by the hardware by a factor of 10 5 , with the minimum time for a full anneal decreasing from 1μs to around 10 ps. 0.0 0.2 0.4 0.6 0.8 1.0 s p − 4.0 − 3.5 − 3.0 − 2.5 − 2.0 − 1.5 − 1.0 − 0.5 0.0 hH IM i ME rate=10 5 μ s − 1 exact ME rate=10 4 μ s − 1 ME rate=10 3 μ s − 1 ME rate=1μ s − 1 DW Figure 2.3: Comparison of numerical (ME) simulations using several different quench rates. The slowest one, 1μs −1 , matches that of the experiments, while the fastest one, 10 5 μs −1 , agrees with the analytical results. Simulations are for a system size of n = 4. Also shown are the exact results and the experimental ones. While performing simulations at larger system sizes becomes computationally prohibitive, we can provide a simple argument that suggests that the annealing rate must be even faster at larger system sizes. For simplicity, let us assume a 21 constant quench rate and an evolution froms =s p tos = 1 that is purely unitary. In order for the unitary dynamics to not change the state significantly we must require T exp   −i ds dt ! −1Z 1 sp H(s 0 )ds 0   −1 ≤ (2.3) for some suitably small , which may also need to decrease with system size in order to reproduce the thermal expectation values to the desired accuracy. Here |·| denotes the operator norm, but any appropriate distance measure can be used. If we expand the time-ordered exponential for small (ds/dt) −1 , it follows from the extensivity of the Hamiltonian that to ensure each term inside the norm remains small, the inverse of the annealing rate (ds/dt) −1 must scale at least as 1/n. Therefore, we already find that simply to ensure that the unitary dynamics does not change the state significantly, the quench rate must become faster as the system size is increased. 2.3.4 Other observables: magnetization We have so far considered the thermal expectation value of H IM as our bench- mark for DW’s behavior. The question arises whether the agreement—or lack thereof—between experimental results and simulations is observable dependent. In Refs. [36, 37], certain observables showed good agreement, while others did not. This is an important point to consider when evaluating DW’s potential as a quan- tum thermal sampler, as consistent behavior across different observables would be required for a fully functional sampler. We choose the squared longitudinal magnetization M 2 z —due to the symmetry of the system,hM z i is always 0. Looking at the overall picture (Fig. 2.4), it is 22 apparent that the two observables follow different patterns of behavior when DW’s results are compared to QMC; whilehH IM i always stays close to the value it should attain at the end of the anneal,hM 2 z i follows a trend that is qualitatively more similar to the QMC data. However, thehM 2 z i results produced by DW do not match QMC anywhere in the anneal (except at one point where they cross), and in fact their relative difference is much larger than forhH IM i for most system sizes (Fig. 2.5). These observations confirm our suspicion that the degree of agreement between DW and QMC can vary across different observables. This adds another layer of complication if we wish to use the annealer to predict thermal expectation values, as its behavior regarding a particular observable is not indicative of what would happen for a different one, and each case would need to be considered individually. 23 0.0 0.2 0.4 0.6 0.8 1.0 s p 0.0 0.2 0.4 0.6 0.8 1.0 hM 2 z i DW QMC 0.0 0.2 0.4 0.6 0.8 1.0 s p − 140 − 120 − 100 − 80 − 60 − 40 − 20 0 hH IM i DW QMC Figure 2.4: Comparison of DW and QMC results for an n = 138 chain at a range of locations through the anneal, forhM 2 z i (top) andhH IM i (bottom). Blue data points correspond to QMC results and green data points to DW with t p = 1.9ms and1500annealspergauge. Errorbarscorrespondtothe95%bootstrapconfidence interval. 24 0 100 200 300 400 500 n 5 10 15 20 |hM 2 z iQMC−h M 2 z iDW| |hM 2 z iQMC| Figure 2.5: Relative difference between the experimental and QMC results for hM 2 z i as a function of problem size for s p = 0.2, t p = 1.9ms and 1500 anneals per gauge. Error bars correspond to the 95% bootstrap confidence interval. 25 2.4 Conclusion In this chapter, we asked whether the current ability to control the annealing scheduleoftheD-Wave2000Qquantumannealingprocessorallowsustoaccurately probethestateofthesysteminthemiddleoftheanneal. Theargumentputforthis that if the system is quenched sufficiently rapidly froms =s p towardss = 1 where the measurement is performed, we would be effectively performing a measurement on the state at s p . We used the nonfrustrated one-dimensional transverse-field Ising model as a case study and found that the D-Wave 2000Q quantum annealing processor cannot be used to reliably reproduce the correct thermal expectation values for the observableH IM . We show in Appendix A.2 that our findings remain unchanged when we consider frustrated chains. We identified the most likely culprit to be the quench rate: with a maximum annealing rate of ds dt = 1μs −1 offered by the device, the system continues to evolve before the measurement occurs. We also verify in Appendix A.3 that fluctuations in the temperature cannot explain the discrepancy between the QMC predictions and our experimental observations. Beyond the experimental difficulty of achieving the necessary fast quench rates, another fundamental problem arises. In our work, we have assumed an ideal qubit (2-level) Hamiltonian, but in the hardware the effective qubit Hamiltonian is real- ized by projecting onto the lowest two energy levels of a superconducting flux qubit [83]. In this effective description, the computational basis is defined in terms of symmetric and anti-symmetric combinations of the ground and first-excited states of the flux qubit Hamiltonian at zero flux, and this basis changes along the annealing schedule [83, 84]. In the unitary dynamics, the above approximation 26 manifests itself as geometric terms in the effective qubit Hamiltonian that con- tributes to the evolution even in the limit of a sudden quench [84]. These effects are not captured by our ideal qubit assumption in Eq. (2.3.3), and while the effect on a single qubit may be small, the error associated with it accumulates with system size. An additional difficulty will be added as we seek to study problems of increas- ing complexity. When embedding is required, its negative effect on the likelihood of obtaining correct solutions is particularly harmful for sampling problems [85], whereweareinterestedinstatesatallenergiesasopposed toonlyGSslikeforopti- mization problems. This effect is worsened with problem size and the complexity of the embedding. Our results, considered along those of Refs. [36, 37]—where experiments on the same platform were able to reproduce the correct thermal expectation val- ues of certain observables (but not others) in far more complicated systems that required embedding—suggest that the ability of the annealer to produce correct mid-anneal predictions is highly dependent on the observable and the problem considered, and their particular susceptibility to the quench. It is likely that for the systems in Refs. [36, 37], the quench did not change the samples in a signifi- cant way for the expectation values of observables that were calculated. However, as the simple example of the 1D TFIM shows, it is impossible to know a-priori whether a system will be susceptible to the slowness of the quench, and hence extreme care must be taken when interpreting results from the variable annealing schedule as ‘measurements-in-the-middle.’ We therefore conclude that currently available quantum annealing devices are not well posed to function as substitutes for quantum Monte Carlo simulators. 27 Chapter 3 Discriminating non-isomorphic graphs with an experimental quantum annealer The contents of this chapter have been produced in collaboration with professors Itay Hen and Klas Markström, and Ruilin Zhou. 3.1 Introduction Inthepreviouschapter, wesawthatusingapause-and-quenchannealingsched- ule is not equivalent to performing an actual measurement in the middle. However, results obtained using said protocol differ from the ones we would obtain through a standard anneal without a pause. This might make the method a viable can- didate to solve certain instances of the Graph Isomorphism problem—the task of deciding whether two graphs are isomorphic, or the same graph up to vertex relabeling. The possibility of using quantum annealers in this manner was first brought up almost a decade ago [46] and later attempted on the experimental annealer available at the time [47]. Graphs with the same classical Ising spectra but different quantum ones were tested, and some of them could be distinguished by looking at refinements of their classical spectra, while for others there was no appreciabledifference. The reason forthelackof success insomecasesisthatwhen 28 measurements were performed, at the end of the anneal, only the classical Ising component of the Hamiltonian remains, and with graphs having the same classical spectra, this measurement at the end is not a sufficient discriminator. But a mech- anism that allows tapping into the quantum spectra—such as a pause-and-qunech schedule—could in principle distinguish between such graphs. We put this hypothesis to the test by running sets of non-isomorphic, co-Ising graphs on a quantum annealer programmed with a reverse schedule with a pause and quench, and find that, even with measurements taken at the end of the anneal, this method succeeds in distinguishing the graphs we test. 3.2 The graph isomorphism problem A graphG = (V,E) is a pair of sets such thatE⊆ [V ] 2 . V is the set of vertices, pairs of which can form the edges belonging to E [86]. We denote the number of vertices|V| by n. Two graphs G and G 0 are said to be isomorphic, G ∼ =G 0 , when there exists a bijection between them that preserves vertex adjacency, that is, one can be converted into the other by a simple relabeling of its vertices. Determining whether two graphs are isomorphic is what we call the Graph Isomorphism (GI) problem. There are a plethora of methods for solving it, tai- lored to specific cases, and in the last few years a new quasi-polynomial general case algorithm [87] deposed the previously exponential one, which had remained unchallenged for 30 years [88]. This problem has also been explored at large from the quantum computing angle [46, 47, 89–93]. Note that the aim of this work is not to beat the performance of classical algorithms for the GI problem, but to prove that this method for solving it, first proposed in [46], is now technically possible to implement. Future work will be 29 necessary to optimize its performance and determine whether it can surpass that of other methods. 3.2.1 Graph invariants A graph invariant is a map that takes a graph as its argument and assigns equal values to isomorphic graphs. In other words, a parameter of the graph that remains unchanged under vertex relabelling. One possible approach to the GI problem is comparing certain invariants. However, most invariants that can be calculated efficiently can also assign the same value to non-isomorphic graphs, and this degeneracy renders them useless for solving the GI problem. An invariant that unequivocally determines the graph—one whose value is necessarily different for non-isomorphic graphs—is called complete. A complete invariant that would be useful for solving the GI problem is not known. Some exist that are incredibly complex [94], and others, while trivial in their definition, would incur too high a cost to calculate. For instance, we could enumerate every finite graph, assigning them a unique index. That index would obviously fully determine the graph, but the enumerating would be too costly. There are many matrix invariants associated to a graph, one of the most ele- mentary ones being the adjacency matrix A(G). A(G) ij = 1 if{ij} is an edge in G and A(G) ij = 0 otherwise. The spectrum of the adjacency matrix is also a common invariant, and, until a counter-example was presented in 1957 [95], it was believed that all graphs were determined by their spectrum—that it was a complete invariant. For the vast majority of small graphs (n≤ 11) this is in fact the case [96], which explains how long this misconception was held. It was later 30 found that cospectral non-isomorphic graphs are in fact fairly common, and many cospectral constructions are currently known [96–98]. Wecanusetheadjacencymatrixasabuildingblocktoconstructinvariants. By assigning to each edge{ij} in the graph the product of individual Pauli matrices σ z i σ z j , we obtain the Ising Hamiltonian, establishing a strong connection between graphs and physics: H I (G) =J X {i,j}∈E(G) σ z i σ z j =J X i,j A(G) i,j σ z i σ z j , (3.1) where J is a constant representing the strength of the interaction. The spectrum of H I (G) is a stronger graph invariant than that of A(G), in the sense that there are fewer non-isomorphic graphs which share this spectrum, as compared to simply co-spectral non-isomorphic graphs. However, the former still exist, and thus the spectrum ofH I (G) is not complete. We call two graphs co-Ising if their respective Ising Hamiltonians have the same spectra [47]. A yet stronger spectral invariant can be obtained by adding a longitudinal field to the previous Hamiltonian: H L (G) =J X i,j A(G) i,j σ z i σ z j +h n X i=1 σ z i , (3.2) with h being the strength of the longitudinal field. Two graphs are said to be longitudinal co-Ising if their respectiveH L have the same spectra. Non-isomorphic, longitudinal co-Ising graphs are also known. Adding a transverse field toH L results in the quantum Ising Hamiltonian [99]: H(G) =J X i,j A(G) i,j σ z i σ z j +h n X i=1 σ z i + Γ n X i=1 σ x i , (3.3) 31 where Γ represents the strength of the transverse field. Whether the spectrum of this Hamiltonian is a complete graph invariant remains an open question—non- isomorphic, quantum co-Ising graphs have not yet been found. This spectrum might not seem particularly helpful for distinguishing graphs, since calculating it entails the diagonalization of a 2 n × 2 n matrix, which quickly becomes intractable asn increases. However, the quantum Ising Hamiltonian is physically implemented by quantum annealers, which could in principle be used to distinguish graphs by their quantum spectra. 3.3 Methods Recall that DW implements a time-dependent, quantum Ising Hamiltonian H(s) =A(s)H d +B(s)H p , where H d = n P i=1 σ x i is the transverse-field driver Hamil- tonian and H p is the Ising type problem Hamiltonian: H p = X hi,ji J ij σ z i σ z j + n X i=1 h i σ z i . (3.4) When all the couplings have the same strength—J ij = J∀i,j—and the same happens with the longitudinal fields—h i = h∀i—, H p is the same as H L from Sec. 3.2, making H(s) a time dependent version of H(G): H(G,s) =A(s)   ˜ J X i,j [A(G)] i,j σ z i σ z j + ˜ h n X i=1 σ z i   + +B(s) n X i=1 σ x i , (3.5) where we have defined ˜ J≡ J/A(s) and ˜ h≡ h/A(s), and set Γ = B(s). In our runs, we always have ˜ J = ˜ h = 1. 32 Because the spectrum of the DW Hamiltonian is a potential complete invariant, it should in principle be possible to use DW for solving the GI problem by taking appropriate measurements. This approach was first proposed in [46], and later attempted in [47] using an older generation annealer. It was expanded upon in [90]. There are two challenges when implementing this method. The first one is that only diagonal measurements can be performed, precluding the direct measurement of the quantum spectrum. This is not too problematic, since expectation values of diagonal observables will still depend on the full quantum Hamiltonian: hOi = Tr[Oe −βH ] Tr[e −βH ] , (3.6) so differences in its spectrum will affect our results. We will measure the classical energy E≡hH p i, the z-magnetization M z , the spin-glass order parameter [46]: Q 2 = v u u t 1 N(N− 1) X i6=j hσ z i σ z j i 2 , (3.7) and the next-nearest neighbor interaction energy: Ω 2 = X i,j [A 2 (G)] i,j σ z i σ z j . (3.8) The second challenge, just like in the previous chapter, is only being able to take measurements at the end of the anneal, when H(s) is classical and only the H p component has non-zero strength. To illustrate why we would like to take measurements far from s = 1, consider Fig. 3.1. It shows the exactly calculated ΔE =E G13 −E G13p for one of the pairs 33 0.0 0.2 0.4 0.6 0.8 1.0 s − 0.03 − 0.02 − 0.01 0.00 E G13 − E G13p Figure 3.1: Exact difference in classical energy E =hH p i throughout the anneal at T= 12mK (D-Wave’s operating temperature) for the smallest pair of graphs we test, G13 and G13p. Their small size (n = 13) allows us to exactly calculate their spectrum by direct diagonalization. of graphs we will test, G13 and G13p, with n = 13 (small enough to directly diagonalize). These graphs look classically identical, as evidenced by the ΔE = 0 at both ends, and their differences only emerge during the quantum regime. A measurement near s = 0.3, where ΔE peaks, would be ideal for distinguishing them, but only measurements at s = 1 are currently possible. This does not necessarily mean that we will be unable to capture the differences in quantum spectra. In [47] the authors had hoped that the noisy, quantum evolution of the non- isomorphic graphs would lead to distinguishable statistics even when measuring only at the end. Some graphs—but not all—could be distinguished in this way, by looking at sorted individual values of a triplet of observables. We are now able to exercise a greater amount of control over how that evolution happens. Results from the previous chapter suggest that allowing the system to reach its thermal 34 state during a pause and then quenching to s = 1 is not equivalent to measuring during the pause, due to the quench not being, in general, sufficiently fast to fully prevent evolution. Measurements might not faithfully represent the thermal state at H(s p ). But for our current purpose, we do not require that the state remain unchanged during the quench, only that enough information is preserved so that the differences between quantum spectra are not completely washed out by the end of the anneal. We will use a reverse annealing schedule with a pause, like the one described in Sec. 1.3. A qualitatively depiction of our schedule is shown in Fig. 1.4. We prepare the system in the initial classical state with all spins up, and anneal back from s = 1 to some intermediate pause location s = s p in 1μs. We will test different pause locations s p ∈ [0.1, 0.2,..., 0.9], and pause there for a time t p = 100μs. We want the system to thermalize during the pause, and the information about that thermal state to be preserved as much as possible before we can measure ats = 1. To this end, we perform the final part of the anneal (from s = s p to s = 1) at the fastest rate allowed by the hardware, ds/dt = 1μs −1 , trying to prevent furhter evolution from the thermal state at s = s p . After the first anneal is completed, we take the final state for this anneal as the initial state for the next, and do the same for each subsequent anneal. To improve our statistics we perform gauge averaging, where a different unitary transformation is applied to the problem instance each time, in order to mitigate intrinsic biases of the physical components. For each run we perform 200 different random gauges, with 1000 anneals each. Bootstrapping is then performed over the gauges to obtain mean values. 35 Because our goal here is to prove that currently available quantum annealers can measure quantum differences in graphs that are classically equivalent, we have used a set of parameters that allowed us to do so. A focused study on the choice of the parameters that impact the time to solution—such as number of anneals, number of gauges, annealing time, pause time—will be needed to quantify and optimize the performance of this method. 3.3.1 Choice of graphs We want to test the ability of the annealer to distinguish graphs based on their quantumspectra. Tothisend, wechoosepairs(orlargertuples)ofnon-isomorphic, co-Ising graphs. Let us define the Ising polynomial, or partition function for the Ising model [100]: Z(G,x,y) = X ~ s∈Ω x E(~ s) y Mz (~ s) , (3.9) where ~ s is a spin configuration, Ω is the space of all possible configurations, and E(~ s) and M z (~ s) are the energy and magnetization of the configuration ~ s, that is, E(~ s) =− P <i,j> s i s j and M z (~ s) = P i s i . Co-Ising graphs have the same Ising polynomial. For smalln (≤ 10), a full classification of equivalence classes of co-Ising graphs is known. The first such examples were given in [100], where they were found via a complete search of all small graphs, and the search was extended in [101], where non-isomorphic graphs which have the same partition function for all two-state spin models were found. But due to the rapidly growing number of non-isomorphic graphs, even when restricted to e.g. regular graphs, and the rarity of co-Ising pairs any complete search will be limited to quite small sizes. 36 However, [100] also introduced a technique for constructing large families of non-ismorphic graphs with the same Ising-polynomial, and this method was extended in [102] to the random-cluster model. Using data from [100, 101] we have here used a version of this method to construct larger graphs to use as bench- marks. This method makes use of the rooted Ising-polynomial Z(G,v,x,y) of a graph G. Given a vertex v in the graph G this polynomial is defined as the usual Ising-polynomial except that we only sum over those states on G in which the spin of v is +1. Note that the full Ising-polynomial is given by Z(G,x,y) = Z(G,v,x,y)+Z(G,v,x, 1/y). So, if two graphsG 1 ,G 2 have the same rooted poly- nomials Z(G 1 ,v 1 ,x,y) and Z(G 2 ,v 2 ,x,y) for some vertex v 1 ∈ G 1 and v 2 ∈ G 2 then they will have the same Ising-polynomial, but two graphs may have the same Ising-polynomial even if no such pair of vertices exist. The rooted polynomial has the useful property that if we start with any two graphs, and root vertices v 1 and v 2 , and build a new graph G 3 by identifying the vertices v 1 and v 2 into a new vertex v 3 then Z(G 3 ,v 3 ,x,y) =Z(G 1 ,v 1 ,x,y)Z(G 2 ,v 2 ,x,y)/y. So, if we know the rooted polynomial for the smaller graphs G 1 and G 2 we can easily find the full Ising-polynomial ofG 3 , and if we start with two different pairs of non-isomorphic graphs with the same rooted polynomials we can build two larger graphs which are Co-Ising, but which may not be isomorphic. Similarly we can use more than two graphs in this product, as long as we divide by the right power of y. 37 Using the examples from [101] we search for pairs of graphs with the same rooted polynomials. Several such pairs exist, and we opted for a few pairs of tree- like starting graphs to keep them native to the Chimera architecture. While minor embedding [42] is a commonly used technique for dealing with non-native graphs, where several physical qubits linked by strong ferromagnetic couplings are used to represent a single node from the original graph, we avoid it for the present study. When we talk about embedding here, we refer only to the choice of physical qubits that will represent graph nodes with a one to one correspondence. Some of the pairs of graphs we study can be embedded using the same set of qubits. When this is the case, we also choose a third one, isomorphic to one of the other two, to ensure that any differences measured between the pair of non-isomorphic graphs is in fact due to their different structure, and not other factors such as noise. Other pairs or tuples have a structure such that they cannot be embedded using the same physical qubits. In this case, several different embeddings are used for each graph, and results averaged over them. For the differences between non- isomorphic graphs to be statistically significant, they must be larger than those found among the group of isomorphic graphs (i.e. the different embeddings). 3.4 Results 3.4.1 G13 and G13p The first pair of graphs{G13,G13p}, withn = 13, is small enough that we can perform exact diagonalization to calculate thermal averages. Fig. 3.1 shows the exact difference in energy E =hH p i between G13 and G13p, and similar results are obtained for other observables—there is a peak at some intermediate time in 38 the anneal but the difference is zero at the ends, confirming that the pair of graphs is classically (but not quantum) co-Ising. We also confirm that these thermal averages do not significantly change for a range of temperatures around 12 mK, D- Wave’s operating temperature, to ensure that any small temperature fluctuations that might occur during the annealing process will not have a large impact on our results. Whileexactresultstellusthat ΔE hasthelargestpeakvalueoftheobservables we survey, experimentally we find that the energy does not distinguish between G13 and G13p as clearly as the magnetization or the spin-glass order parameter (shown in Fig. 3.2) do. This is likely due to the fact that our pause-and-quench schedule is not able to mimic a true measurement in the middle, and that different observables are not equally susceptible to the quench, like we saw in the previous chapter. This pair of graphs was previously tested—and distinguished—in [47]. How- ever, to discriminate between the two graphs they needed to rely on refinements of the classical spectra—namely counting the individual occurrences of a triplet of observables—while, thanks to the introduction of a pause, we can do so by just considering the average value of a single observable. We embed G13 and G13p on DW using the same set of qubits, since their structure allows it. Alongside them, we run a third graph, G13i, isomorphic to G13, as a baseline for our results. Fig. 3.2 shows that DW can pick up the differ- ences in structure between G13 and G13p, while correctly finding G13 and G13i indistinguishable. 39 3.4.2 G27 and G27p The pair{G27,G27p} was also tested in [47], but found to be indistinguish- able even when looking at the individual occurrences for the triplet of observables {e,m, Ω 2 }—individual energy, magnetization and next-nearest neighbor interac- tion energy, respectively—which could tellG13 andG13p apart. The next-nearest neighbor interaction energy is defined as Ω 2 = P i,j [A 2 (G)] i,j σ z i σ z j . With the introduction of a pause, this triplet can now clearly tell the two graphs apart, as shown in Fig. 3.3 (top). In fact, we no longer need to look at individual occurrences, ashΩ 2 i suffices to distinguish between G27 and G27p (Fig. 3.3 bottom). Just like with the previous pair of graphs, these two can be embedded using the same qubits on DW, and we also run a third graph G27i isomorphic to G27. The average energy, magnetization and spin-glass order parameter cannot distinguish them, but Ω 2 can. This attests that the introduction of a pause can allow us to obtain information on the quantum spectrum that was previously inaccessible. 3.4.3 G17 and G17p After succeeding to distinguish the graphs used in [47], we tackle the new pairs and tuples that we have constructed, starting with the pair{G17,G17p}, with n = 17. Unlike the previously tested graphs, these two cannot be embedded using the same set of qubits. This provides an opportunity to start investigating how different embeddings affect our outcomes. While in the previous cases we added a third graph that was isomorphic to the first one in the pair, here we need to proceed with more caution to ensure that measured differences are not due to using different physical qubits and couplers. 40 We will run five different embeddings for each graph, and average results over them. This amounts to testing five isomorphic graphs, rather than just two. We also do it for both graphs in the pair, instead of just one. If the uncertainty for the group of isomorphic graphs is smaller than the difference between the non-isomorphic ones, then their differences, as measured by DW, are statistically significant. We start off with embeddings that use sets of qubits with the largest possible overlap. G17 and G17p cannot be embedded using the same set of qubits, but they can—in several different ways—using two sets that differ in only one qubit. We choose four different embeddings for each graph satisfying that condition: for embeddings number 1 through 4, G17 and G17p differ in only one qubit. And for each graph individually, embeddings with consecutive numbers differ in one qubit as well (as do embeddings number 4 and 1). Once we confirm that averaging over these four embeddings produces statis- tically different outcomes for G17 and G17p, we run a final embedding for each graph (number 5), using sets of qubits that are disjoint from the previous four, while still only differing in one qubit for G17 and G17p. The top panel in fig. 3.4 shows that G17 and G17p yield very different results for Q 2 , even when averaged over the five embeddings. The error bars are small compared to the difference between the two graphs. Q 2 also remains robust across the five embeddings (bottom two panels of fig. 3.4). There is no qualitative differ- ence between the more similar embeddings (1− 4) and the one with the disjoint qubit set (5). We also find that other observables—E,hM z i, Q 2 ,hΩ 2 i—are able to tell the two graphs apart in the same manner. 41 3.4.4 Set of four G25 graphs The next test case is a 4-tuple of graphs{G25p1,G25p2,G25p3,G25p4} with n = 25. Given our encouraging results on the previous pair when using several embeddings, we take it one step further here and find, for each of the four graphs, five fully different embeddings: all the qubit sets are disjoint and they also have different structures. This is a necessary step towards scaling up the use of quantum annealers to solve the GI problem. Even though we are manually picking our embeddings for this study—since we want to have full control over their characteristics and our goal is not yet to optimize performance—that will not remain viable at larger scales, where manually finding embeddings will be at best highly inefficient, and we must require that automatically generated embeddings yield correct results. For instance, DW’s embedding heuristic will usually return different embed- dings when run multiple times for the same graph. Although it is possible to choose a target graph—a subgraph of the full Chimera where we wish to place the logical graph—we do not expect in general to have enough information to be able to determine it. And if we want to run several isomorphic graphs to obtain trustworthy statistics, using several embeddings will likely still be necessary. We find that DW steps up to the task, and despite the added difficulty of distinguishing four similar graphs rather than only two, both the energy and the spin-glass order parameter can clearly tell them apart when pausing at an inter- mediate point in the anneal, as evidenced by Fig. 3.5. 42 3.4.5 Set of four G33 graphs The final set we test also contains four graphs{G33,G33p,G33p1,G33p2}. With n = 33, it is the largest that has been attempted so far. The first two graphs—G33 and G33p—can be embedded using the same qubits while the last two—G33p1 and G33p2—require different ones. Similarly to what we did for our first test cases, we embed G33 and G33p on the same qubits along with a third graph isomorphic to G33 (G33i). For G33p1 and G33p2, on the other hand, we run four completely different embeddings for each and average over them. All four graphs are distinguishable by the spin-glass order parameter with a pause located around s p = 0.7 (Fig. 3.6). 43 0.2 0.4 0.6 0.8 s p − 0.3 − 0.2 − 0.1 0.0 0.1 0.2 Δ hM z i hM z i G13 −h M z i G13i hM z i G13 −h M z i G13p 0.2 0.4 0.6 s p − 0.0050 − 0.0025 0.0000 0.0025 0.0050 0.0075 Δ Q 2 Q 2G13 − Q 2G13i Q 2G13 − Q 2G13p Figure 3.2: Experimental results for G13 and G13p. Top: Results for ΔhM z i for different pause locationss p . The red line shows the difference inhM z i between thetwoisomorphicgraphsG13 andG13i, whichremainsclosetozeroasweexpect. The green line shows the difference between the two non-isomorphic graphs G13 and G13p, which is non-zero, allowing us to tell them apart. Data points show the mean after performing a bootstrap over the 200 gauges, with the error bars at the 95% confidence interval. Bottom: Same for ΔQ 2 . These results confirm the annealer’s ability to distinguish this pair of graphs by their quantum spectra. The large error bars fors p close to 1 are likely due to our use of reverse annealing. Dynamics are slower in this region and we never reach s < s p . Because we reuse the final state at each anneal as the initial for the next, small random changes in the initial configuration are more likely to remain, and be different across gauges, leading to more disparate results when compared to annealing back to a smallers p . 44 m=-1 e=-34 Ω 2 =64 m=-1 e=-34 Ω 2 =68 m=-1 e=-34 Ω 2 =72 m=-1 e=-36 Ω 2 =68 m=-1 e=-36 Ω 2 =80 m=-3 e=-34 Ω 2 =60 m=-3 e=-34 Ω 2 =64 m=1 e=-34 Ω 2 =68 m=1 e=-34 Ω 2 =80 0 100 200 300 400 500 600 s p = 0.1 G27 G27i G27p 0.2 0.4 0.6 0.8 s p − 10.0 − 7.5 − 5.0 − 2.5 0.0 2.5 Ω 2 Ω 2 G27 − Ω 2 G27i Ω 2 G27 − Ω 2 G27p Figure 3.3: Top: Occurrences for the triplet of observables{e,m, Ω 2 } when a pause is introduced at s p = 0.1. G27 and G27i show the same results while G27p is clearly different. Bottom:hΩ 2 i suffices to discriminate betweenG27 andG27p. The red line shows the difference between the two isomorphic graphs G27 and G27i, staying relatively constant around 0, while the green line represents the difference between the non-isomorphic graphs G27 andG27p, which remains non- zero throughout. In both cases the mean is calculated from bootstrapping over the 200 gauges and error bars show the 95% confidence interval. 45 0.2 0.4 0.6 0.8 s p 0.35 0.40 0.45 0.50 0.55 Q 2 G17 G17p 0.2 0.4 0.6 0.8 s p 0.35 0.40 0.45 0.50 0.55 Q 2 G17 emb 1 G17 emb 2 G17 emb 3 G17 emb 4 G17 emb 5 0.2 0.4 0.6 0.8 s p 0.35 0.40 0.45 0.50 0.55 Q 2 G17p emb 1 G17p emb 2 G17p emb 3 G17p emb 4 G17p emb 5 Figure 3.4: Top: The spin-glass order parameterQ 2 clearly distinguishes between G17 and G17p when introducing a pause throughout most of the anneal. Results are averaged over 5 different embeddings for each graph. Middle and bottom: The spin-glass order parameter Q 2 remains robust across 5 different embeddings for G17 (middle) and G17p (bottom). For each embedding, the mean and 95% confidence interval are shown (obtained from bootstrap over gauges). Then the mean and standard error are calculated across the different embeddings. 46 0.2 0.4 0.6 0.8 s p − 28.0 − 27.5 − 27.0 − 26.5 − 26.0 − 25.5 hH p i G25p1 G25p2 G25p3 G25p4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 s p 0.3 0.4 0.5 0.6 0.7 Q 2 G25p1 G25p2 G25p3 G25p4 Figure3.5: TopTheclassicalenergyhH p iandBottomthespin-glassorderparam- eter Q 2 can distinguish between all four G25 graphs, with results averaged over five different embeddings for each graph, when pausing in the highlighted region. Error bars show the standard error of the mean, and individual embedding results were first calculated through a bootstrap over the 200 gauges. 47 0.2 0.4 0.6 0.8 s p 0.3 0.4 0.5 0.6 0.7 Q 2 G33 G33i G33p2 0.2 0.4 0.6 0.8 s p 0.3 0.4 0.5 0.6 0.7 Q 2 G33 G33p G33p1 G33p2 Figure3.6: TopThespin-glassorderparameterQ 2 issimilarforthetwoisomorphic graphs (red and green lines) throughout the anneal, while the non-isomorphic one has a different result for a pause at s p = 0.7. Bottom Q 2 can distinguish all four G33 graphs with a pause at s p = 0.7. Error bars show the standard error of the mean, and individual embedding results were first calculated through a bootstrap over the 200 gauges. 48 3.5 Conclusions and outlook In this work we have investigated how the use of a reverse annealing schedule with a pause boosts the ability of a quantum annealer to act as a discriminator for graphs with the same classical spectra but different quantum ones. We used the quantum annealer as a sampler, rather than its more standard role as an optimizer, since our goal is not to get a specific result, but instead to implement a test that can bring out the differences in quantum spectra between classically co-spectral but non-isomorphic graphs. The observables we measure would have led to equal outcomes in a purely classical regime. By introducing a mid-anneal pause followed by a quench we were able to access information about the full quantum Hamiltonian that typically remains out of reach for standard annealing schedules. This pause is usually necessary to access the region where the differences take place. Even though pausing and quenching as currently implemented is not a good substituteformeasuringinthemiddle—andinfactwhenexactresultsareavailable mid-anneal they do not match those calculated from the annealer’s output—the method works as intended by preserving enough information about the differences in quantum spectra that were present at the time of the pause. In fact, we were able to distinguish all the pairs and 4-tuples of graphs that we tested by looking at averages for one or several observables, where at least one of them—and usually more—showed statistically significant differences between the non-isomorphic graphs. Wehavethusextendedtheworkofpreviousstudiesonthepotentialofquantum annealers for solving certain hard instances of the graph isomorphism problem, 49 bringing to fruition the idea that classically co-spectral graphs could be distin- guished by their quantum spectra. We have also investigated how the use of different native embeddings affects these results. To make a fair comparison between two graphs, we also test at least another isomorphic one (i.e. the same graph with a different embedding) to establish a baseline. We find that the differences between non-isomoprhic graphs remain statistically significant even when averaged over several embeddings, and whether those embeddings are similar or not among themselves does not appear to have a discernible impact on the results. The structure of the graphs prevails over potential biases from the components and noise. This is a crucial step in the quest to expand the scope of this method. With the advent of new architectures with increased connectivity that will soon succeed Chimera, the pool of native graphs will expand significantly. To solve the graph isomorphism problem on a quantum annealer for larger, more connected graphs, and to do it more efficiently, being able to rely on heuristic methods for embedding- finding will likely be key. Robustness across embeddings will be even more essential for extending the method to non-native graphs requiring minor-embedding. When the logical graph does not have a one-to-one correspondence with the physical one, and instead a group of physical qubits must be coerced into acting together as a single logical variable, the potential for different outcomes from logically isomorphic graphs will greatly increase. Establishing a strong foundation for native graphs will bring us closer to successfully tackling the non-native case. Future work should also address performance. While the present study focuses on feasibility, a rigorous investigation of performance and its scaling with the size 50 of the graphs will be necessary to test this method against classical algorithms. Differences between reverse and forward annealing with a pause should be consid- ered, and the length of the pause, as well as other annealing parameters, adjusted to optimize the time to solution. 51 Chapter 4 Ferromagnetically shifting the power of pausing The contents of this chapter have been produced in collaboration with Drs. Shon Grabbe, Stuart Hadfield, Jeffrey Marshall, Zhihui Wang and Eleanor Rieffel. 4.1 Introduction It was previously observed that for a family of Ising problems of connectivity conforming to that of the hardware qubits, (native problems), a pause during the annealing can significantly boost the success probability in obtaining a ground state [38]. Strikingly, the location of the most effective pause is insensitive to the problem specifics and is shown to be universal for an ensemble of problem instances. In this chapter, we provide evidence that for a family of non-native problems—degree-bounded spanning tree problems—a universal pausing location still exists that can effectively boost both the success probability in finding the optimal solution, and the expected time-to-solution (TTS). For most application problems, on a hardware with restricted qubit connectiv- ity, minor embedding is required to enable coupling between two qubits that are not physically connected, where one logical qubit is represented by a set of physi- cal qubits with strong ferromagnetic coupling of magnitude|J F | 1 among them to 1 J F is always negative, so we always refer to its magnitude|J F |. 52 promote collective behavior. We call this set of physical qubits a vertex model. While|J F | can be set to a large value such that the embedded problem preserves the ground state of the logical problem, and analytical bounds on this value can be obtained [42], the energy spectrum in general varies with the value of|J F |, and its effect on the annealing often requires careful case-by-case consideration [103–105]. We also demonstrate that the annealing dynamics can be further influenced by varying the parameter|J F |. In particular, the optimal region for the pause location and the corresponding gain in success probability are both affected by the value of|J F |. A pause placed earlier in the anneal is more beneficial for a larger|J F | value. We extended the physical picture of the annealing developed in [38], where the interplay between adiabaticity and thermalization accommodates the boost to success probability by a pause, and we further describe how the shift in optimal pause location with|J F | is consistent in this picture. Furthermore, we performed perturbative analysis on the effect of|J F | on the minimalenergygapbetweenthegroundandthefirstexcitedstates,andournumer- ical simulations on small problems demonstrate the shift of minimal energy, both consistent with the shift observed in the annealing results of the degree-bounded minimum spanning tree problem set. More specifically, in this chapter, we show that: • Forafamilyofcombinatorialoptimizationproblemsthatrequiresembedding, specifically the bounded-degree minimal spanning tree (BD-MST) problems, an appropriate pause during the annealing can boost the success probability. The gain is dependent on the pause duration. • Choosing|J F | in the “extended" range, which is a newly added feature in DW, helps in the case of annealing without pausing. 53 • As we increase the value of|J F |, a systematic shift in the optimal pause location can be observed. • Along with the extension of implementable|J F | to an asymmetric range, a general strategy, gauge transformation, or, spin inverse transformation, that has been shown effective in suppressing noise, is yield invalid. We develop a partial gauge transformation techniquethatsignificantlyimprovesthesuccess probability, and enables the confirmation of the major claims in this paper. • A physical picture of the interplay between the minimal gap location, non- adiabaticity, and thermalization can be established, and it is consistent with the peak shift with|J F |. Our gap analysis and numerical simulations on small systems show that as|J F | increases, the minimal gap shifts to earlier and the gap size decreases, which backs up our understanding. • Withappropriatelengthandduration, apausecanleadtostatisticallysignif- icant improvement in the time to solution (TTS) for the set of spanning-tree problems we studied, compared with the best TTS achieved without pausing. BD-MST problems play a significant role within the family of network robust- ness problems. Our study—in particular the improvement demonstrated in the TTS—opens a pathway to leverage quantum technology for a broad spectrum of network related problems. Our physical picture and analysis of the shift in optimal pausing location shed light on further understanding and maneuvering quantum annealing for practical applications. The rest of the chapter is organized as follows: In Sec. 4.2, we define the problem and detail the parameters used in drawing the problem instances and in the annealing. Sec. 4.3 is devoted to results obtained from the annealer. Results 54 for annealing without pause are shown in Sec. 4.3.1, how pausing can be helpful is demonstrated in Sec. 4.3.2 and how pausing shifts with|J F | in Sec. 4.3.3. The technical treatment that enables the conclusive results, partial gauges, is discussed in Sec. 4.3.4, and we provide theoretical analysis and a physical picture for the shifting of optimal pause location with|J F | in Sec. 4.3.5. In Sec. 2.4 we summarize the results and discuss future work. 4.2 Problems, quantum annealing and setup 4.2.1 Problem definition Given an undirected graph G = (V,E) and set of weights w uv for each edge (uv)∈E, we seek a spanning tree T⊂E such that the tree weight P huvi∈T w uv is minimized. Definition 4.2.1.1. A spanning tree for a graph G is a subgraph of G that is a tree and contains all vertices of G. Spanning trees are important for several reasons. They play an important role in designing efficient routing algorithms. Some computationally hard problems, such as the Steiner tree problem and the traveling salesperson problem, can be solved approximately using spanning trees [106]. Spanning tree problems also find broad applications such as network design, bioinformatics, etc. For general graphs, determining if there exists a spanning tree of weight W can be decided in polynomial time, and different efficient algorithms exist to find a minimum weight tree; for example, using Kruskal’s algorithm which requires time O(|E| log|V|) [107]. (Special classes of graphs can be solved even faster.) 55 On the other hand, with the additional constraint that the maximum vertex degree of the spanning tree is at most Δ, even deciding whether there exists such a spanning tree becomes NP-complete for fixed Δ≥ 2 [108]. Problem I: Given an integer Δ≥ 2 and graph G = (V,E) with edge weights w uv , (uv) ∈ E, find a minimum weight spanning tree of maximum degree at most Δ. Finding a degree bounded spanning tree of cost at mostr times that of the opti- mum cost remains NP-hard for any r≥ 1 [109]. Hence, approximation algorithms aretypicallydesignedtoreturnaminimumspanningtreewithvertexdegrees‘close to’ Δ. In Refs. [110–112] a polynomial time algorithm is given which returns a spanning tree of degree at most Δ+1 and weight at mostOPT, whereOPT is the weight of the (unknown) optimal degree Δ spanning tree. Alternatively, heuristics exist which return valid degree bounded spanning trees but with suboptimal cost that may be difficult to quantify generally. A wide variety of approaches have been developed for this problem [113–118], including specific approximations for various special cases (e.g., geometric weights); see [119] for an extensive overview. 4.2.2 Solving on a quantum annealer Tosolveonaquantumannealer, anoptimizationproblemneedstobeexpressed in a quadratic unconstrained binary optimization (QUBO) form—or equivalently Ising—, and the binary variables therein can be readily translated into qubits, which are called logical qubits. A number of mappings of the BD-MST problem to QUBO can be found in [120]. We apply a resource-efficient mapping, the tight level-based mapping (see Appendix C.1 for details). 56 The resulting QUBO must further be embedded to conform with the hardware connectivity. In the embedding, we use the same coupling strength|J F | for all the couplings that connect two physical qubits that represent the same logical qubit (what we call a vertex model). Our study will have a focus on how|J F | affects the annealing performance. The shape of the annealing schedule, i.e., how the Hamiltonian is varied over time from the driving Hamiltonian H d to the problem (or cost) Hamiltonian H p , is fixed on the current hardware. However, one can insert a pause (i.e., holding the Hamiltonian constant for a specified time). The duration of this pause, as well as its location in the anneal, can be set by the user. In a previous study [38], it was shown that a pause at a location insensitive to the problem specifics boosts the probability of finding the ground state—the success probability—by orders of magnitude. This observation was attributed to the fact that, if introduces shortly after the minimum gap, a pause can give the system sufficient time to thermalize, at a point where the contribution of the GS to the thermal state is becoming significant, thus increasing the probability of finding the system in the GS when measurements are taken at the end. Associated with embedding, optimally setting the ferromagnetic coupling|J F | of the vertex model is a challenging task [104]. DW allows asymmetric extension of the pairwise qubit coupling strengthsJ i,j ∈ [−2, 1] (in addition to the canonical symmetric option J i,j ∈ [−1, 1]). One usage of this extension is to set|J F | in the extended range. We show how the extended values improve the success probability of our problems. On the other hand, the asymmetry yields invalid the general strategy of gauge transformation, which we normally use for its effectiveness in reducing noise effects 57 and obtaining higher quality output data. We designed a novel strategy, partial gauge transformation, that selectively applies the transformation only to couplings in the symmetric range [−1, 1]. For the case that only the embedding couplings are intheextendedrange, thisisequivalenttoapplyingageneralgaugetransformation to the logical problem prior to the embedding, and is simple to implement. We found that the partial gauge transformation significantly helped in both boosting the success probability and reducing the output variance. This strategy enabled emerging of the key features we observed: the positive role of an extended|J F | in the case of no pausing, and the shift of the optimal pause location with|J F |. 4.2.3 Hardware, annealing parameters, and problem instance specifics Each BD-MST problem instance consists of a weighted graph G = (V,E) and a degree bound Δ. The underlying graphs are chosen by exhausting all connected graphs with n =|V| = 5, which are listed in Table C.1 and have m =|E| ranging from 4 to 10. Weight sets are taken from the list in Table C.2. Graphs and weight sets can be combined to yield a large number of unique instances. Results are averaged over an ensemble of instances, the number of which will be specified along with the results in Section 4.3. For each instance, the objective Hamiltonian is calculated according to the tight level-based mapping described in Appendix C.1. For the degree bound we generally selected Δ = 2; we also tested Δ = 3 and our claims were found to hold for this value as well. The objective Hamiltonian is normalized to fit in the range [−1, 1], so that the extended J range is saved for the embedded qubits only. We chose a|J F | in the 58 range [1, 2]—initially exploring all values in that range at 0.1 intervals—and then used the annealer’s built-in heuristic algorithm for embedding. For the graphs we consider, the probability of finding an embedding is∼ 100%. The algorithm is run 30 times and the embedding with the smallest size (total number of physical qubits) is chosen. Detailed information about the typical size of the embedded problems for different graphs can be found in Appendix C.3, including the number of physical qubits and the size of vertex models. Embedding statistics for a future D-Wave architecture (Pegasus) are also shown. After the embedded problem is run on DW, outputs with any inconsistent values on physical qubits that represent the same logical qubit—or with violated penalty terms such that the output doesn’t encode a degree bounded spanning tree—are considered to be invalid answers, and discarded as failed runs. The retained valid answers are then verified against the exact solution of the problem, which are obtained through direct enumeration for the small problem sizes we consider. Several annealing timest a were initially tested for no-pause runs, and the short- est time allowed by DW, t a = 1μs, was found to be optimal in terms of TTS across problem instances for the instance ensembles. This agrees with previous studies [20, 24, 105]. Unless otherwise specified, runs are performed witht a = 1μs, 50000 anneals (or reads), and 100 partial gauges. For the pause runs, the whole range of possible pause locations (s p ∈ [0, 1]) was initially surveyed, confirming the existence of a peak in p success . Although the location of the peak is affected by|J F |, it always stays within the range [0.2, 0.5], so further runs were limited to this area of interest, withs p varied between 0.2 and 0.5 at 0.02 intervals. 59 A range of pause times t p was also studied, and shorter pause times found to yield better TTS. Most of our runs are performed with pause duration t p = 1μs. After optimal values for other parameters were found, additional t p values in the range [0.25, 2]μs were explored. 4.3 Results We use the probability of success (p success ) and time to solution (TTS) as our figures of merit for determining how likely a problem is to be solved, defined as: p success = # anneals with correct solution total # anneals (4.1) TTS = log(1− 0.99) log(1.0−p success ) t tot , (4.2) where the total timet tot takes into account both annealing and pause times, as well as the number of anneals. The annealer returns the solution with the minimum cost it has found. To ensure the validity of this solution, we first confirm that the resulting graph is in fact a spanning tree that satisfies the degree constraint, and also a true optimal solution by comparing with the minimum cost obtained by an exact classical algorithm. Any other outcome is weighted zero toward p success . TTS is a figure of merit reporting the projected time required to solve the problem with 99% confidence. These two measures are complementary to each other. While p success is directly determined by and hence provides a portal to understand the underlying physical process, TTS gives a more practical measure that is universal across different parameters and different solvers. Note that a higher success probability does not necessarily mean a higher TTS. For instance, 60 we might get a higher p success by using a longer annealing time t a = 100μs rather than a shorter one t a = 1μs, however, it also takes 100 times longer to do each run, and we may find that the shorter t a outperforms the longer one in TTS. 4.3.1 Annealing without pause, effect of|J F | We first establish the baseline results with no pause andt a = 1μs, which is the shortest that DW allows, and was chosen for consistently yielding the best TTS. By exploring the available range of|J F | values between 1 and 2, we confirm the advantage of using the extended|J F | range and identify its optimal value for the base case at|J F | ∗ = 1.6, as shown in Fig. 4.1 where the success probability and TTS are shown for a range of|J F | for the ensemble of instances with n = 5. This makes sense physically as stronger couplings within vertex models make it less likely for individual qubits to flip, which helps to avoid breaking up the model. Butatoolarge|J F |wouldmakeitincreasinglycostlyfortheclustertofliptogether, potentially preventing the system from leaving a non-optimal configuration. To boost the probability of success,|J F | must strike the right balance, leading to better chances of arriving at—and staying in—the correct configuration. 4.3.2 Improvement with pause Fig.4.2showsforademoinstance, at|J F | ∗ , thesuccessprobabilityasafunction of the pause location s p ∈ [0, 1]. We observe that the optimal pause location when|J F | = 1.6 is in a range consistent with that reported for the case of native-to-hardware Ising problems, (typicallyarounds p = 0.4in [38]). Weobserveanincreaseinsuccessprobabilityof about an order of magnitude for 1μs annealing with a pause of 100μs (see fig. 4.6, 61 − 2.0 − 1.9 − 1.8 − 1.7 − 1.6 − 1.5 − 1.4 J ferro 10 − 4 10 − 3 p success − 2.0 − 1.9 − 1.8 − 1.7 − 1.6 − 1.5 − 1.4 J ferro 10 3 10 4 10 5 TTS(μ s) Figure 4.1: Optimal|J F | for baseline. Top: Success probability for an ensemble of 45 instancesn = 5, 50K reads, 100 gauges, 1μs anneal and|J F | varied from 1.4 to 2.0 in increments of 0.1. The best performance is observed at|J F | ∗ = 1.6. Data points show the median obtained after bootstrapping over the set of instances. Error bars are at the 35 and 65 percentiles. Bottom: The corresponding TTS. while for shorter pausing like 1 to 2μs, a median gain of 3-4x is observed in 45 instances. For certain problems, the no-pause annealing failed to find a solution 62 0.20 0.25 0.30 0.35 0.40 0.45 0.50 s p 10 − 3 10 − 2 10 − 1 10 0 p success Jf =− 2 Jf =− 1.9 Jf =− 1.8 Jf =− 1.7 Jf =− 1.6 Jf =− 1.5 Jf =− 1.4 Jf =− 1.3 Jf =− 1.2 Jf =− 1.6 no pause Figure 4.2: Shift of optimal pause location. Left: Probability of success versus the annealing pause location for the demo (n = 5,m = 6ver1,w14,emb26) instance. The anneal was performed with 1 μs anneal time, 100 μs pause, 50K reads, and 100 partial gauges applied. A monotonic shift in the peak location with |J F | is observed. The horizontal curve corresponds to the no pause,|J F | = 1.6, 100 partial gauges, and anneal time of 1 μs. The reason for lower success probability for|J F | =1.2 is detailed in Sec. 4.3.3. even after 50K reads, while the annealing with a pause was able to find one. A physical picture supporting this phenomenon will be presented in Sec. 4.3.5. The probability of success grows monotonically with the pause duration in the range we studied (t a ∈ [0.25, 100]μs). However, the TTS not only depends on p success , but also scales with the total time spent on the annealing t a +t p , with t p being the pause duration. After sweeping through a range of pause duration values, we found that att p ∈ [0.25, 2]μs, we obtain an advantage in TTS compared to the best annealing results without pause. See Fig.4.3 where we show, for an ensemble of 45 instances with n = 5, that if the pause duration is set between 0.25μs and 2μs, the median TTS is lower (and thus better) than the case without pausing. 63 To better showcase how pausing helps, we define the TTS improvement ΔTTS = TTS(no pause)−TTS(pause), (4.3) with the two TTS values calculated for their respective optimal|J F | (|J F | ∗ = 1.6 for the no pause case and 1.8 for the pause case). A positive ΔTTS indicates that a pause improves upon the baseline results (i.e. reduces TTS). Fig 4.4 shows how ΔTTS remains consistently positive for t p ∈ [0.25, 2]μs and s p = 0.32. This improvement, while modest (about 30% for the optimal pause duration of 1μs), is statistically significant (except for t p = 0.5μs). 4.3.3 Optimal pause location shift with|J F | Oneinterestingnewavenuethatopensupwiththestudyofembeddedproblems is how the value of|J F | affects the benefits and effects of pausing. As previously discussed, thep success vss p curve typically shows a peak around an optimal pause location and is mostly flat far away from it. We have also seen that different|J F | lead to different p success . But when both parameters vary, rather than obtaining similarly positioned peaks for which the heights change with|J F |, we find that the peak position shifts toward earlier in the anneal, as can be seen in Fig. 4.2 for the demo instance. Such clear shifting is found in many instances, and results in a shift in the behavior of the whole instance ensemble, as shown in Fig. 4.5. Note that the success probability for smaller|J F | = 1.2, 1.3, even away from the peak, is clearly lower than for larger|J F | values. The reason is that when the ferromagnetic coupling is not very strong compared to the problem couplings, it is 64 more likely that the low-lying energy states are densely populated by states with an inconsistent vertex model (i.e. when not all the qubits are aligned and hence are no longer acting as a single variable). Accordingly, even when the annealer is doing well at finding the ground state or a low-lying state, such outcomes do not correspond to a valid solution of the original problem. Indeed, by applying simulated annealing to solve the embedded problem (which is too large to diagonalize exactly), we verified that in the range of|J F | we used, the ratio of inconsistent-vertex -model-state in the GS and low energy states is significantly higher for|J F | = 1.2, 1.3 than that for|J F | = 1.4 and above, indicating that between|J f | = 1.3 and 1.4, the ground state of the embedded problem is no longer a logical configuration. The mechanism for why the optimal pause location shifts toward earlier in the anneal with|J F | fits our theoretical understanding, which is laid out in sec- tion 4.3.5. 4.3.4 Help of partial gauges We developed a partial gauge transformation technique that significantly improved the success probability, and enabled the confirmation of the peak shift. WhenJ ij ∈ [−1, 1], it is straightforward to apply gauges. For our embedded prob- lems, however, we are making use of the extended J range, allowing J ij ∈ [−2, 1]. The extended range discourages the breaking of vertex models during annealing thanks to the stronger ferromagnetic couplings between physical qubits represent- ing the same logical variable, but it also impedes the use of standard gauges, since any couplings in the range [−2,−1) cannot change sign. 65 Our partial gauge method circumvents this issue by only applying the gauge transformation on the couplings within the interval [−1, 1]. Because the extended range is exclusively used on the vertex models in our problems, the partial gauge on the embedded problem is equivalent to applying a general gauge to the logical problem before embedding. Inapreviousstudy[121], theboostinp success bypausingisobservedforafamily of embedded problems, but no relation between the optimal pausing location and J F was observed. In our study, with the help of partial gauge transformation, the variance in the annealing output is significantly suppressed, resulted in the revelation of the shift of the peak in Sec. 4.3.3. Another benefit was a remarkable increase inp success , even with just 10 gauges. Ordinarily, we don’t expectp success to change much from gauge averaging. Solving the problem without gauge transformations is the same as picking one random gauge, which most times will be near average and not an outlier. Performing more gauges will lead to some of them being better and some worse, which in most cases would even out and not affect the initial result dramatically. But the problems we aresolvingherearedifficultforDW,andtheirtypicalp success isverycloseto0. The existence of such a lower bound explains the significant help in applying gauges: even if we get a bad gauge,p success cannot go below 0, so it is impossible to have an outlier in the negative direction. A good gauge, on the other hand, could yield a much higherp success . Doing as few as 10 gauges we are already likely to encounter at least one positive outlier, and that is the reason for the large improvement in p success . A higher number of gauges, while still able to benefit the results, will no longer have such a dramatic effect, as the spread in gauge quality will be similar to the 10 gauge case. 66 Fig 4.6 shows the increase in probability of success for an ensemble of 117 instances when 100 gauges are used. With the partial gauge transformation we are able to extend the benefits of general gauge averaging to embedded problems, while improving their results as difficulty increases. 4.3.5 Picture for how pausing and|J F | affect p success In this section we provide a physical picture explaining the shift of the optimal pausing location, numerical evidence on the change in minimal gap location and a perturbation analysis supporting these claims. How pausing helps We start with a review of the physical picture [38] that could account for the increase in success probability when introducing a pause in the middle of the annealing schedule, slightly after the minimal energy gap. We will use GS and FES to refer to the ground and the first excited states of the instantaneous quantum Hamiltonian. In the rest of the section we refer to the gap as the energy gap between the GS and the FES. At very early or late stages in the annealing, only one Hamiltonian—either the driver H d or the problem H p —dominates. Since both the problem Hamiltonian and the driving Hamiltonian are classical when acting alone, dynamics in these regions are almost classical. Because the temperature T is much lower than the energy scale, thermal relaxation rates remain slow. In the middle of the anneal, when the scales of H d and H p are comparable, the interplay between gap, adiabaticity (annealing speed relative to the gap) and 67 temperature will determine the system dynamics. This is where we expect to find significant population loss from the GS. In particular, when the gap is small enough, approximate instantaneous thermalization may occur, populating excited states. This region is also where non-adiabiatc transitions are expected to be most significant. We thus distinguish three different regimes in the anneal: Regime I: ||A(s)H d ||||B(s)H p ||. The instantaneous Hamiltonian is mainly H d , and its energy scale is much larger than the temperature, T. The system stays in the ground state of H d . Regime II:||A(s)H d ||∼||B(s)H p ||, and their energy scale is comparable to the temperature. Both thermal and quantum dynamics happen, and the minimal gap occurs in this region. Thermalization and non-adiabaticity are both contribut- ing to populating the FES (compared to the case of zero-temperature adiabatic evolution in which all population is in the GS). As the anneal goes on in this regime, it sequentially goes through the following stages: (a) Gap approaching temperature, system leaving adiabatic regime, but transi- tions (non-adiabatic and thermal) may still be relatively slow compared to the system evolution. (b) Gap is around its minimum and is much smaller than the temperature— thermalization happens almost instantaneously and the system is near the thermal equilibrium state. Quantum non-adiabatic effects could be strong enough to increase the population of the FES beyond its magnitude at ther- mal equilibrium. Whether this happens or not doesn’t affect the following picture of how pausing helps. 68 (c) Gap is larger than temperature, non-adiabaticity is weak. The system may still approximately equilibrate if given enough time (e.g. a pause), but will not fully thermalize during the standard evolution. At stage II (b), the instantaneous equilibration removes the state memory from the history. The system is simply in thermal equilibrium. Due to the closeness between the GS and FES, the FES is significantly populated. As the system enters stage II (c) from II (b), pausing would promote better thermalization, which could bring significant FES population in stage II (c) down to the GS, since relative to the gap the temperatureT is now lower, hence boosting the success probability. Regime III:||B(s)H p ||||A(s)H d ||,||B(s)H p || T, known as the frozen region in the literature [82]. Dynamics are slow and the state of the system will not change significantly. We provide a cartoon schematic of these different regions of interest during the anneal in Fig. 4.7. How|J F | shifts the optimal pausing location earlier An increase in|J F | would correspond to the shift of the minimal gap to earlier in the anneal, meaning that stage II. (b) occurs earlier in the anneal. The shift of the minimal gap can partly result from the increase in the relative norm of H p . In Figure. 4.8, for a small Ising problem embedded to 4 physical qubits with Chimera connectivity—which allows exact diagonalization of the instantaneous Hamiltonian—we show the change of minimal gap with|J F |. Proofsketchofgapscalingunder|J F |: Weapplyfirstordernon-degenerate perturbation theory. Let H(s) = H 0 (s) + B(s)λH F with H 0 (s) = A(s)H d + 69 B(s)H p +B(s)J F H F whereλ> 0,J F < 0andH F istheferromagneticHamiltonian for the vertex model. That is, we are considering the effect of weakening the vertex model infinitesi- mally by decreasing|J F |. Let us also take the large|J F |>|J ij | limit (and of course |J F |λ). We start with a simplified case, assuming the only vertex model is a chain of length 2 (the general case will be discussed below). Then H F = σ z k 1 σ z k 2 for two spins k 1 ,k 2 . Let|E i (s)i be the instantaneous ith eigenstate of H 0 (s). From now on, we consider s fixed at some value, and drop the explicit dependence. We can always decompose the instantaneous eigenstates in the computational basis, |E i i = X j a (i) j |z L j i + X k b (i) k |z B k i, (4.4) where|z L j i are logical states, and|z B k i has the chain broken. We compute the matrix elements hE i |H F |E i i = X j |a (i) j | 2 − X k |b (i) k | 2 . (4.5) Note, by normalization P k |b (i) k | 2 = 1− P j |a (i) j | 2 , and, denoting the logical proba- bility P (i) L := P j |a (i) j | 2 , hE i |H F |E i i = 2P (i) L − 1. (4.6) This tells us, to first order inλ> 0, that the low lying energy levels experience anincreaseinenergyupondecreasingtheferromagneticstrength(|J F |→|J F |−λ), i.e., E 0 i =E i +Bλ(2P (i) L − 1)>E i , assuming that P (i) L > 1/2 (which it should be for the GS and FES for s > 0). We see consistent behaviour with this picture in Fig. 4.9 (even though this figure is not in the perturbative limit). 70 Now, the gap Δ =E 1 −E 0 changes under λ, to first order, as Δ 0 = Δ +Bλ(hE 1 |H F |E 1 i−hE 0 |H F |E 0 i) (4.7) = Δ + 2Bλ(P (1) L −P (0) L ). (4.8) Note that at the start of the anneal,|E 0 (0)i =|+i ⊗N , and so P (0) L = 1/2 (in the specific case when the embedding contains just one additional qubit). At s = 0, the FESs are linear combinations containing one excitation in the x eigenbasis, i.e., a single|−i. Consider the symmetric FES, denoted|FE + i, where the state of the chain is 1 √ 2 (|− +i +| +−i) (and the other qubits are all|+i). This state is entirely in the logical subspace, due to the cancelling out of the|01i and|10i terms 2 . When the transverse field is ‘strong’, i.e., ‘near’ s = 0, by perturbation theory we indeed expect that Δ 0 > Δ. We see this in Fig. 4.8 (where the strongest chain,|J F | = 8, has the smallest gap, in the small s region). We also know that once the transverse field becomes weak relative to the prob- lem Hamiltonian (e.g. A(s)/B(s)< 1), then P (1) L −P (0) L → 0 as both the GS and FES become close to logical states. By interpolating between the two extremes, the above argument explains what we see in Fig. 4.8. Another argument that can come into play is that an increase in|J F | would yield “clusters” (physical qubits representing the same logical qubits) with stronger internal force. Changes of state in such clusters would take collective flipping of qubits to achieve (hence more demanding to quantum dynamics) and would incur a higher energy change (requiring a higher relative temperature). 2 A similar argument can be made in the case of longer chains: in the case of a length K chain, the symmetric state with singlex excitation hasP (1) L (s = 0) =K/2 K−1 >P (0) L (s = 0) = 1/2 K−1 . 71 Accordingly, the transition from stage II (b) to II (c) would happen earlier in the anneal. This effect might also be responsible for the less dramatic increase in success rate compared to the native Ising case: the associated energy barrier may require much higher relative temperature, while pausing earlier helps, the amount it can help is limited (because there is an interplay of the three influences which are correlated in a given annealing schedule and at a given temperature). 72 0.30 0.32 s p 10 − 3 6× 10 − 4 2× 10 − 3 3× 10 − 3 p success s p = 0.3 s p = 0.32 no pause t p = 0.25 μ s t p = 0.5 μ s t p = 0.75 μ s t p = 1.0 μ s t p = 2.0 μ s 0.30 0.32 s p 10 4 4× 10 3 6× 10 3 TTS (μ s) s p = 0.3 s p = 0.32 no pause t p = 0.25 μ s t p = 0.5 μ s t p = 0.75 μ s t p = 1.0 μ s t p = 2.0 μ s Figure 4.3: Effect of pause duration on p success and TTS. Left: With pause duration of{0.25, 0.5, 0.75, 1, 2} μs, and|J F |=1.8, the success probability for an ensemble of 45 instances for problem size n = 5 is shown for pause locations s p = 0.3 and s p = 0.32 (which we found to be optimal). The reference (horizontal line and band for median and 35 and 65 percentiles, respectively) is the no-pausing case with parameters optimal for TTS: t a = 1μs, and|J F |=1.6. Data points show the median, with error bars at the 35 and 65 percentiles, after performing 10 5 bootstraps over the set of instances. Right: the corresponding TTS. 73 0.28 0.30 0.32 0.34 0.36 s p − 3000 − 2000 − 1000 0 1000 2000 3000 ΔTTS ( μ s) s p = 0.3 s p = 0.32 t p = 0.25 μ s t p = 0.5 μ s t p = 0.75 μ s t p = 1.0 μ s t p = 2.0 μ s Figure 4.4: Improvement in TTS. ΔTTS represents the reduction in TTS accomplished by introducing a short pause at an optimal location s p . A posi- tive ΔTTS indicates that the TTS is reduced (improved) by the introduction of the pause. The 1μs pause is found to be optimal, improving TTS by 2· 10 3 μs (from 7· 10 3 μs for the no pause case). ΔTTS is calculated by subtracting the TTS for the pause case with with|J F | = 1.8 from that of the no pause case with |J F | = 1.6 (the optimal|J F | for each case). Data points are the median, error bars are 35 th and 65 th percentile obtained from 10 5 bootstraps over 45n = 5 instances. Data points are staggered for readability. 74 0.20 0.25 0.30 0.35 0.40 0.45 0.50 s p 10 − 3 10 − 2 p success |J f | = 1.2 |J f | = 1.4 |J f | = 1.6 |J f | = 1.8 |J f | = 2 Figure 4.5: Shift of optimal pause location (ensemble) with|J F |. Prob- ability of success for an ensemble of 9 instances with n = 5, t p = 100 μs, and t a = 1μs. As|J F | increases, the optimal pause location moves earlier in the anneal. Data points represent the median, obtained by bootstrapping over the different instances. Error bars are at the 35 and 65 percentiles. 0.20 0.25 0.30 0.35 0.40 0.45 0.50 s p 10 − 4 10 − 3 10 − 2 p success 100 gauges no pause no gauges no pause 100 gauges no gauges Figure 4.6: Effect of partial gauges. Partial gauges help boost median success probability for a set of 117 instances. |J F | = 1.6 is used for all data, and the baseline case without pause is also shown. For the pause results, t p = 100μs. The median is obtaied by bootstrapping over the different instances, and error bars (or bands) are at the 35 and 65 percentiles. 75 Figure 4.7: Cartoon diagram for the different regimes. Colored in purple the leftandrightmostregionsaretheadiabaticRegimesIandIII(whichfurtherextend to the left and right as indicated by the arrows). Regime 2 is subdivided into three regions, a,b and c, determined by the relative magnitudes of the instantaneous gap Δ and the temperature T. In IIa we expect the system to stop behaving strictly adiabatically, and in IIb instantaneous thermalization may occur if the relaxation time scale is small enough, as well as non-adiabatic transitions. In IIc, a pause may help repopulate the GS. This should be thought of as an approximate picture of what occurs, to aid the reader. In reality the transitions between these regions will not be sharp and defined by a single point during the evolution. 76 Figure 4.8: Shift of minimal gap with|J F |. Energy gap between the ground and the first excited states for the instantaneous quantum Hamiltonian during annealing for a toy problem of a complete graph of size 3, embedded to 4 physical qubits on Chimera connectivity. The gap is computed exactly by diagonalizing the instantaneous Hamiltonian. As|J F | increases, the instantaneous gap closes, and the minimum gap shifts to earlier in the anneal. Figure 4.9: Energy level shift with|J F |. The individual energy levels show an increase in energy upon decreasing the ferromagnetic strength (i.e. the case when λ> 0 in our perturbation theory). This is for a problem with a single chain of size 2 (complete graph of size 3, embedded size 4). 77 4.4 Conclusions and outlook In this work, we have studied the performance of a quantum annealer on the previouslyunconsideredclassofdegree-boundedminimumspanningtreeproblems, in one of the first applications of novel annealing schedules to embedded problems. We confirm that, like in the the previously reported case for native Ising prob- lems, anappropriatelytimedpausecanimprovetheprobabilityofsuccessandtime to solution for these problems. The region of the optimal pause location remains consistent across problem instances, as well as with the one obtained for the native case. We found that the boost provided by a pause is of lesser magnitude than for native problems. This is likely due to the more costly nature of changing configurations in an embedded problem; all qubits belonging to the same vertex model must flip together as a cluster, otherwise the model breaks and the solution is discarded. We also studied the effect of extending the range of coupling strength to allow for stronger ferromagnetic bonds, with J ∈ [−2, 1]. We showed that choosing |J F |—the coupling within a vertex model—in the extended range can further increase the probability of success, an observation which is supported by the phys- ical picture we present: stronger couplings within vertex models make it less likely for individual qubits to flip, which helps to avoid breaking up the model. But a too large|J F | would make it increasingly costly for the cluster to flip together, potentially preventing the system from leaving a non-optimal configuration. To boost the probability of success,|J F | must strike the right balance, leading to better chances of arriving at—and staying in—the correct configuration. 78 The choice of|J F | also impacts the optimal pause location, which shifts earlier in the anneal as|J F | increases. This is consistent with our theoretical picture: pausing increases the probability of success by promoting thermalization when done shortly after the minimum gap, making it easier for the population that migrated to an excited state to relax back to the ground state. Because a larger |J F | leads to an earlier minimum gap and a decrease in its size, the optimal pause location moves earlier as well. In order to use the extended J range without sacrificing the benefits provided by gauge transformations (which in their standard form are not compatible with an asymmetricalJ range) we developed a partial gauge transformation technique. This new approach can be applied to a general asymmetric J case. In partic- ular, when the extended range is only used for the embedded qubits, applying gauges to the rest amounts to gauge transforming the original unembedded prob- lem. We show that, as with the standard version, performing partial gauges helps curb uncertainty and, when applied to hard problems, the probability of success experiences a significant boost. We have studied the interplay between quantum annealing parameters in embedded problems, providing both deeper insights into the physics of these devices and pragmatic recommendations to improve performance on optimization and sampling problems. Even though the gain in TTS for this specific set of problem instances is not by orders of magnitude, the statistically significant gain is promising and the under- standingacquiredinthestudyopensavenueforfurthertailoringannealingparam- eters to improve the performance for practical optimization problems. 79 Our study offers a first step in understanding how to increase the solubility of problemsrequiringembeddingthroughtheoptimizationofannealingschedulesand embedding parameters. While we expect future hardware to vastly expand the set of native problems—such as the upcoming Pegasus architecture for D-Wave which will greatly improve upon Chimera’s connectivity—until an architecture featuring a complete graph or fully customizable couplers becomes available, embedding will remain a necessity as we strive to solve increasingly complex problems closer to real-world applications. Future work should tackle new classes of embedded problems, delve deeper into embedding mechanisms and features, and explore how these interact with the annealing parameters and the performance of quantum annealers on embedded problems in general. 80 Reference List [1] B.Apolloni, C.Carvalho, andD.deFalco. Quantumstochasticoptimization. Stochastic processes and their applications, 33, Dec 1989. [2] A. B. Finnila, M. A. Gomez, C. Sebenik, C. Stenson, and J. D. Doll. Quan- tum annealing: A new method for minimizing multidimensional functions. Chemical Physics Letters, 219(5–6):343–348, 3 1994. [3] J. Brooke, D. Bitko, T. F., Rosenbaum, and G. Aeppli. Quantum annealing of a disordered magnet. Science, 284(5415):779–781, 1999. [4] Tadashi Kadowaki and Hidetoshi Nishimori. Quantum annealing in the transverse Ising model. Phys. Rev. E, 58(5):5355, November 1998. [5] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, Joshua Lapan, Andrew Lundgren, and Daniel Preda. A quantum Adiabatic Evolution Algorithm Applied to Random Instances of an NP-Complete Problem. Science, 292(5516):472–475, April 2001. [6] Giuseppe E. Santoro, Roman Martoňák, Erio Tosatti, and Roberto Car. TheoryofquantumannealingofanIsingspinglass. Science, 295(5564):2427– 2430, 2002. [7] [8] Philipp Hauke, Helmut G. Katzgraber, Wolfgang Lechner, Hidetoshi Nishi- mori, and William David Oliver. Perspectives of quantum annealing: Meth- ods and implementations. 2019. [9] M. Born and V. Fock. Beweis des Adiabatensatzes. Z. Physik, 51(3-4):165– 180, March 1928. [10] Nikolai Il‘in, Anastasia Aristova, and Oleg Lychkovskiy. Adiabatic theorem for closed quantum systems initialized at finite temperature. arXiv e-prints, page arXiv:2002.02947, February 2020. 81 [11] J. E. Avron, M. Fraas, G. M. Graf, and P. Grech. Optimal time schedule for adiabatic evolution. Phys. Rev. A, 82:040304, Oct 2010. [12] A. Elgart and J.E. Avron. An adiabatic theorem without a gap condition: Two level system coupled to quantized radiation field. Physical Review A, 58, 03 1998. [13] Scott Kirkpatrick, C. Gelatt, and M. Vecchi. Optimization by simulated annealing. Science (New York, N.Y.), 220:671–80, 06 1983. [14] R. Harris, J. Johansson, A. J. Berkley, M. W. Johnson, T. Lanting, Siyuan Han, P. Bunyk, E. Ladizinsky, T. Oh, I. Perminov, E. Tolkacheva, S.Uchaikin,E.M.Chapple,C.Enderud,C.Rich,M.Thom,J.Wang,B.Wil- son, and G. Rose. Experimental demonstration of a robust and scalable flux qubit. Phys. Rev. B, 81(13):134510, April 2010. [15] R.Harris, M.W.Johnson, T.Lanting, A.J.Berkley, J.Johansson, P.Bunyk, E. Tolkacheva, E. Ladizinsky, N. Ladizinsky, T. Oh, F. Cioata, I. Perminov, P. Spear, C. Enderud, C. Rich, S. Uchaikin, M. C. Thom, E. M. Chapple, J. Wang, B. Wilson, M. H. S. Amin, N. Dickson, K. Karimi, B. Macready, C. J. S. Truncik, and G. Rose. Experimental investigation of an eight- qubit unit cell in a superconducting optimization processor. Phys. Rev. B, 82:024511, Jul 2010. [16] Mark Johnson, Mohammad Amin, S Gildert, Trevor Lanting, F Hamze, N Dickson, R Harris, Andrew Berkley, J Johansson, Paul Bunyk, E Chapple, C Enderud, Jeremy Hilton, Kamran Karimi, E Ladizinsky, Nicolas Ladizin- sky, T Oh, I Perminov, C Rich, and Geordie Rose. Quantum annealing with manufactured spins. Nature, 473:194–8, 05 2011. [17] Wim van Dam, Michele Mosca, and Umesh Vazirani. How Powerful is Adi- abatic Quantum Computation? arXiv e-prints, pages quant–ph/0206003, May 2002. [18] Marko Žnidarič and Martin Horvat. Exponential complexity of an adiabatic algorithm for an NP-complete problem. Phys. Rev. A, 73(2):022329, Febru- ary 2006. [19] Sergio Boixo, Tameem Albash, Federico M. Spedalieri, Nicholas Chancellor, and Daniel A. Lidar. Experimental signature of programmable quantum annealing. Nat. Commun., 4:2067, 06 2013. [20] Sergio Boixo, Troels Frimodt Rønnow, Sergei Isakov, Zhihui Wang, David Wecker, Daniel Lidar, John Martinis, and Matthias Troyer. Evidence for 82 quantum annealing with more than one hundred qubits. Nat Phys, 10, 04 2013. [21] John Smolin and Graeme Smith. Classical signature of quantum annealing. Frontiers in Physics, 2, 05 2013. [22] Seung Woo Shin, Graeme Smith, John A. Smolin, and Umesh Vazirani. How “Quantum” is the D-Wave Machine? arXiv e-prints, page arXiv:1401.7087, January 2014. [23] Tameem Albash, Walter Vinci, Anurag Mishra, Paul Warburton, and Daniel Lidar. Consistency tests of classical and quantum models for a quantum annealer. Physical Review A, 91:042314, 04 2015. [24] Troels F. Rønnow, Zhihui Wang, Joshua Job, Sergio Boixo, Sergei V. Isakov, David Wecker, John M. Martinis, Daniel A. Lidar, and Matthias Troyer. Defining and detecting quantum speedup. Science, 345(6195):420–424, 2014. [25] Victor Martin-Mayor and Itay Hen. Unraveling quantum annealers using classical hardness. Scientific Reports, 5, 02 2015. [26] Salvatore Mandrà, Helmut Katzgraber, and Creighton Thomas. The pitfalls of planar spin-glass benchmarks: Raising the bar for quantum annealers (again). Quantum Science and Technology, 2, 03 2017. [27] Salvatore Mandrà and Helmut Katzgraber. A deceptive step towards quan- tum speedup detection. Quantum Science and Technology, 3, 11 2017. [28] Itay Hen, Joshua Job, Tameem Albash, Troels Frimodt Rønnow, Matthias Troyer, and Daniel Lidar. Probing for quantum speedup in spin glass prob- lems with planted solutions. Physical Review A, 92, 02 2015. [29] HelmutKatzgraber,FirasHamze,ZhengZhu,AndrewOchoa,andHumberto Munoz-Bauza. Seeking quantum speedup through spin glasses: The good, the bad, and the ugly. Physical Review X, 5, 05 2015. [30] Mohammad H. Amin, Evgeny Andriyash, Jason Rolfe, Bohdan Kulchytskyy, and Roger Melko. Quantum boltzmann machine. Phys. Rev. X, 8:021050, May 2018. [31] Richard Y. Li, Tameem Albash, and Daniel A. Lidar. Improved Boltz- mann machines with error corrected quantum annealing. arXiv e-prints, page arXiv:1910.01283, October 2019. 83 [32] Andrew King, Juan Carrasquilla, Isil Ozfidan, Jack Raymond, Evgeny Andriyash, Andrew Berkley, Mauricio Reis, Trevor M. Lanting, Richard Har- ris, Gabriel Poulin-Lamarre, Anatoly Yu. Smirnov, Christopher Rich, Fabio Altomare, Paul Bunyk, Jed Whittaker, Loren Swenson, Emile Hoskinson, Yuki Sato, Mark H. Volkmann, and Mohammad Amin. Observation of topo- logical phenomena in a programmable lattice of 1,800 qubits. Nature, 560, 03 2018. [33] R. Harris, Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, T. Lanting, R. Li, T. Medina, R. Molavi, R. Neufeld, T. Oh, I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov, L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, and J. Yao. Phase transitions in a programmable quantum spin glass simulator. Science, 361(6398):162–165, 2018. [34] Andrew D. King, Jack Raymond, Trevor Lanting, Sergei V. Isakov, Masoud Mohseni, Gabriel Poulin-Lamarre, Sara Ejtemaee, William Bernoudy, Isil Ozfidan, Anatoly Yu. Smirnov, Mauricio Reis, Fabio Altomare, Michael Babcock, Catia Baron, Andrew J. Berkley, Kelly Boothby, Paul I. Bunyk, Holly Christiani, Colin Enderud, Bram Evert, Richard Harris, Emile Hoskin- son, Shuiyuan Huang, Kais Jooya, Ali Khodabandelou, Nicolas Ladizinsky, Ryan Li, P. Aaron Lott, Allison J. R. MacDonald, Danica Marsden, Gaelen Marsden, Teresa Medina, Reza Molavi, Richard Neufeld, Mana Norouzpour, Travis Oh, Igor Pavlov, Ilya Perminov, Thomas Prescott, Chris Rich, Yuki Sato, Benjamin Sheldan, George Sterling, Loren J. Swenson, Nicholas Tsai, Mark H. Volkmann, Jed D. Whittaker, Warren Wilkinson, Jason Yao, Hart- mut Neven, Jeremy P. Hilton, Eric Ladizinsky, Mark W. Johnson, and Mohammad H. Amin. Scaling advantage in quantum simulation of geometri- cally frustrated magnets. arXiv e-prints, page arXiv:1911.03446, November 2019. [35] Yuki Bando, Yuki Susa, Hiroki Oshiyama, Naokazu Shibata, Masayuki Ohzeki, Fernand o Javier Gómez-Ruiz, Daniel A. Lidar, Adolfo del Campo, Sei Suzuki, and Hidetoshi Nishimori. Probing the Universality of Topologi- cal Defect Formation in a Quantum Annealer: Kibble-Zurek Mechanism and Beyond. arXiv e-prints, page arXiv:2001.11637, January 2020. [36] R. Harris, Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, T. Lanting, R. Li, T. Medina, R. Molavi, R. Neufeld, T. Oh, I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov, L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, and 84 J. Yao. Phase transitions in a programmable quantum spin glass simulator. Science, 361(6398):162–165, 2018. [37] Andrew King, Juan Carrasquilla, Isil Ozfidan, Jack Raymond, Evgeny Andriyash, Andrew Berkley, Mauricio Reis, Trevor M. Lanting, Richard Har- ris, Gabriel Poulin-Lamarre, Anatoly Yu. Smirnov, Christopher Rich, Fabio Altomare, Paul Bunyk, Jed Whittaker, Loren Swenson, Emile Hoskinson, Yuki Sato, Mark H. Volkmann, and Mohammad Amin. Observation of topo- logical phenomena in a programmable lattice of 1,800 qubits. Nature, 560, 08 2018. [38] Jeffrey Marshall, Davide Venturelli, Itay Hen, and Eleanor G. Rieffel. Power of pausing: Advancing understanding of thermalization in experimental quantum annealers. Phys. Rev. Applied, 11:044083, Apr 2019. [39] Masaki Ohkuwa, Hidetoshi Nishimori, and Daniel A. Lidar. Reverse anneal- ing for the fully connected p -spin model. Physical Review A, 98:022314, August 2018. [40] G. Passarelli, V. Cataudella, and P. Lucignano. Improving quantum anneal- ing of the ferromagnetic p-spin model through pausing. Phys. Rev. B, 100:024302, Jul 2019. [41] P. I. Bunyk, E. M. Hoskinson, M. W. Johnson, E. Tolkacheva, F. Altomare, A. J. Berkley, R. Harris, J. P. Hilton, T. Lanting, A. J. Przybysz, and J.Whittaker. Architecturalconsiderationsinthedesignofasuperconducting quantum annealing processor. IEEE Transactions on Applied Superconduc- tivity, 24(4):1–10, Aug 2014. [42] Vicky Choi. Minor-embedding in adiabatic quantum computation: I. The parameter setting problem. Quant. Inf. Proc., 7(5):193–209, 2008. [43] VickyChoi. Minor-embeddinginadiabaticquantumcomputation: II.Minor- universal graph design. Quant. Inf. Proc., 10(3):343–353, 2011. [44] Davide Venturelli and Alexei Kondratyev. Reverse quantum annealing approachtoportfoliooptimizationproblems. Quantum Machine Intelligence, pages 17–30, 2019. [45] Masaki Ohkuwa, Hidetoshi Nishimori, and Daniel A. Lidar. Reverse anneal- ing for the fully connectedp-spin model. Phys. Rev. A, 98:022314, Aug 2018. [46] Itay Hen and A. P. Young. Solving the graph-isomorphism problem with a quantum annealer. Phys. Rev. A, 86:042310, Oct 2012. 85 [47] Walter Vinci, Klas Markström, Sergio Boixo, Aidan Roy, Federico M. Spedalieri, Paul A. Warburton, and Simone Severini. Hearing the Shape of the Ising Model with a Programmable Superconducting-Flux Annealer. Scientific Reports, 4:5703, July 2014. [48] Richard P. Feynman. Simulating physics with computers. Int. J. Theor. Phys., 21:467–488, 1982. [49] Seth Lloyd. Universal quantum simulators. Science, 273(5278):1073–1078, 08 1996. [50] J. Ignacio Cirac and Peter Zoller. Goals and opportunities in quantum sim- ulation. Nat. Phys., 8(4):264–266, 04 2012. [51] Andrew M. Childs, Dmitri Maslov, Yunseong Nam, Neil J. Ross, and Yuan Su. Towardthefirstquantumsimulationwithquantumspeedup. Proceedings of the National Academy of Sciences, 115(38):9456–9461, 2018. [52] Yudong Cao, Jonathan Romero, Jonathan P. Olson, Matthias Degroote, Peter D. Johnson, Mária Kieferová, Ian D. Kivlichan, Tim Menke, Borja Peropadre, Nicolas P. D. Sawaya, Sukin Sim, Libor Veis, and Alán Aspuru- Guzik. Quantum Chemistry in the Age of Quantum Computing. arXiv e-prints, page arXiv:1812.09976, Dec 2018. [53] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Michael Sipser. Quan- tum Computation by Adiabatic Evolution. arXiv:quant-ph/0001106, Jan- uary 2000. [54] Dorit Aharonov, Wim van Dam, Julia Kempe, Zeph Landau, Seth Lloyd, and Oded Regev. Adiabatic quantum computation is equivalent to standard quantum computation. SIAM J. Comput., 37(1):166–194, January 2007. [55] Dorit Aharonov and Amnon Ta-Shma. Adiabatic quantum state generation and statistical zero knowledge, 2003. [56] Tosio Kato. On the Adiabatic Theorem of Quantum Mechanics. J. Phys. Soc. Jpn., 5(6):435–439, November 1950. [57] SabineJansen,Mary-BethRuskai,andRuediSeiler. Boundsfortheadiabatic approximation with applications to quantum computation. J. Math. Phys., 48(10):102111, 2007. [58] David Ceperley and Berni Alder. Quantum monte carlo. Science, 231(4738):555–560, 1986. 86 [59] David Landau and Kurt Binder. A Guide to Monte Carlo Simulations in Statistical Physics. Cambridge University Press, New York, NY, USA, 2005. [60] M.E.J. Newman & G.T. Barkema. Monte Carlo Methods in Statistical Physics. Oxford Uinversity Press, 1999. [61] A. W. Sandvik. A generalization of Handscomb’s quantum Monte Carlo scheme — Application to the 1-D Hubbard model. J. Phys, A, 25:3667, 1992. [62] N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn. Exact, complete, and universal continuous-time worldline monte carlo approach to the statistics of discrete quantum systems. Journal of Experimental and Theoretical Physics, 87(2):310–321, 1998. [63] Anders W. Sandvik. Stochastic series expansion method with operator-loop update. Phys. Rev. B, 59:R14157–R14160, Jun 1999. [64] Tameem Albash, Gene Wagenbreth, and Itay Hen. Off-diagonal expansion quantum monte carlo. Phys. Rev. E, 96:063309, Dec 2017. [65] R. D. Somma, C. D. Batista, and G. Ortiz. Quantum approach to classical statistical mechanics. Phys. Rev. Lett., 99:030603, Jul 2007. [66] Hidetoshi Nishimori, Junichi Tsuda, and Sergey Knysh. Comparative study of the performance of quantum annealing and simulated annealing. Phys. Rev. E, 91(1):012104–, 01 2015. [67] Heinz-Peter Breuer and Francesco Petruccione. The Theory of Open Quan- tum Systems. Oxford University Press, 2002. [68] E.B.Davies. Markovianmasterequations. Communications in Mathematical Physics, 39(2):91–110, 1974. [69] Tameem Albash, Sergio Boixo, Daniel A Lidar, and Paolo Zanardi. Quantum adiabaticMarkovianmasterequations. New J. of Phys., 14(12):123016, 2012. [70] R. Haag, N. M. Hugenholtz, and M. Winnink. On the equilibrium states in quantum statistical mechanics. Comm. Math. Phys., 5:215–236, 1967. [71] Walid K. Abou Salem. On the Quasi-Static Evolution of Nonequilibrium Steady States. Ann. Henri Poincaré, 8(3):569–596, May 2007. [72] Alain Joye. General Adiabatic Evolution with a Gap Condition. Commun. Math. Phys., 275(1):139–162, October 2007. 87 [73] E. B. Davies and H. Spohn. Open quantum systems with time-dependent hamiltonians and their linear response. J. Stat. Phys., 19:511–523, 1978. [74] Lorenzo Campos Venuti, Tameem Albash, Daniel A. Lidar, and Paolo Zanardi. Adiabaticity in open quantum systems. Phys. Rev. A, 93:032118, Mar 2016. [75] R. Harris, J. Johansson, A. J. Berkley, M. W. Johnson, T. Lanting, Siyuan Han, P. Bunyk, E. Ladizinsky, T. Oh, I. Perminov, E. Tolkacheva, S.Uchaikin,E.M.Chapple,C.Enderud,C.Rich,M.Thom,J.Wang,B.Wil- son, and G. Rose. Experimental demonstration of a robust and scalable flux qubit. Phys. Rev. B, 81:134510, Apr 2010. [76] M W Johnson, P Bunyk, F Maibaum, E Tolkacheva, A J Berkley, E M Chapple, R Harris, J Johansson, T Lanting, I Perminov, E Ladizinsky, T Oh, and G Rose. A scalable control system for a superconducting adiabatic quantum optimization processor. Superconductor Science and Technology, 23(6):065004, 2010. [77] A J Berkley, M W Johnson, P Bunyk, R Harris, J Johansson, T Lanting, E Ladizinsky, E Tolkacheva, M H S Amin, and G Rose. A scalable read- out system for a superconducting adiabatic quantum optimization system. Superconductor Science and Technology, 23(10):105014, 2010. [78] P. I Bunyk, E. M. Hoskinson, M. W. Johnson, E. Tolkacheva, F. Altomare, AJ. Berkley, R. Harris, J. P. Hilton, T. Lanting, AJ. Przybysz, and J. Whit- taker. Architectural considerations in the design of a superconducting quan- tum annealing processor. IEEE Transactions on Applied Superconductivity, 24(4):1–10, Aug. 2014. [79] S.V.IsakovandR.Moessner. Interplayofquantumandthermalfluctuations in a frustrated magnet. Phys. Rev. B, 68:104409, Sep 2003. [80] Subir Sachdev. Quantum Phase Transitions. Cambridge University Press, 2 edition, 2011. [81] Tameem Albash, Sergio Boixo, Daniel A Lidar, and Paolo Zanardi. Quan- tum adiabatic markovian master equations. New Journal of Physics, 14(12):123016, dec 2012. [82] Mohammad H. Amin. Searching for quantum speedup in quasistatic quan- tum annealers. Phys. Rev. A, 92:052323, Nov 2015. [83] Sergio Boixo, Vadim N. Smelyanskiy, Alireza Shabani, Sergei V. Isakov, MarkDykman,VasilS.Denchev,MohammadH.Amin,AnatolyYuSmirnov, 88 Masoud Mohseni, and Hartmut Neven. Computational multiqubit tunnelling in programmable quantum annealers. Nat Commun, 7, 01 2016. [84] Walter Vinci and Daniel A. Lidar. Non-stoquastic hamiltonians in quantum annealing via geometric phases. npj Quantum Information, 3(1):38, 2017. [85] JeffreyMarshall, AndreaDiGioacchino, andEleanorG.Rieffel. ThePerilsof Embedding for Sampling Problems. arXiv e-prints, page arXiv:1909.12184, Sep 2019. [86] Reinhard Diestel. Graph Theory (Graduate Texts in Mathematics). Springer, 5 edition, 2016. [87] László Babai. Graph Isomorphism in Quasipolynomial Time. arXiv e-prints, page arXiv:1512.03547, December 2015. [88] L. Babai, W. M. Kantor, and E. M. Luks. Computational complexity and the classification of finite simple groups. In Proceedings of the 24th Annual Symposium on Foundations of Computer Science, SFCS ’83, pages 162–171, Washington, DC, USA, 1983. IEEE Computer Society. [89] Scott D. Berry and Jingbo B. Wang. Two-particle quantum walks: Entan- glement and graph isomorphism testing. Phys. Rev. A, 83:042317, Apr 2011. [90] Frank Gaitan and Lane Clark. Graph isomorphism and adiabatic quantum computing. Physical Review A, 89:022342, February 2014. [91] Kenneth M. Zick, Omar Shehab, and Matthew French. Experimental quan- tum annealing: case study involving the graph isomorphism problem. Sci- entific Reports, 5:11168 EP –, 06 2015. [92] P. W. Mills, R. P. Rundle, J. H. Samson, Simon J. Devitt, Todd Tilma, V. M. Dwyer, and Mark J. Everitt. On quantum invariants and the graph isomorphism problem. arXiv e-prints, page arXiv:1711.09842, November 2017. [93] Cristian S. Calude, Michael J. Dinneen, and Richard Hua. Qubo formula- tions for the graph isomorphism problem and related problems. Theoretical Computer Science, 701:54 – 69, 2017. At the intersection of computer science with biology, chemistry and physics - In Memory of Solomon Marcus. [94] Matthias Dehmer, Frank Emmert-Streib, and Martin Grabner. A computa- tional approach to construct a multivariate complete graph invariant. Infor- mation Sciences, 260:200–208, 03 2014. 89 [95] L. Von Collatz and U. Sinogowitz. Spektren endlicher grafen. Abhandlungen aus dem Mathematischen Seminar der Universitat Hamburg, 21(63-77), 1957. [96] Edwin R. van Dam and Willem H. Haemers. Which graphs are determined by their spectrum? Linear Algebra and its Applications, 373:241 – 272, 2003. Combinatorial Matrix Theory Conference (Pohang, 2002). [97] Chris Godsil and B D. McKay. Constructing cospectral graphs. Aequationes Mathematicae, 25:257–268, 12 1982. [98] Michael Langberg and Dan Vilenchik. Constructing cospectral graphs via a newformofgraphproduct. Linear and Multilinear Algebra, 66(9):1838–1852, 2018. [99] Elliott Lieb, Theodore Schultz, and Daniel Mattis. Two soluble models of an antiferromagnetic chain. Annals of Physics, 16(3):407 – 466, 1961. [100] Daniel Andrén and Klas Markström. The bivariate Ising polynomial of a graph. Discrete Appl. Math., 157(11):2515–2524, 2009. [101] Klas Markström. From the Ising and Potts models to the general graph homomorphism polynomial. In Graph polynomials, Discrete Math. Appl. (Boca Raton), pages 123–138. CRC Press, Boca Raton, FL, 2017. [102] Klas Markström and John C. Wierman. Aperiodic non-isomorphic lattices with equivalent percolation and random-cluster models. Electron. J. Com- bin., 17(1):Research Paper 48, 14, 2010. [103] Vicky Choi. The effects of the problem hamiltonian parameters on the mini- mum spectral gap in adiabatic quantum optimization. Quantum Information Processing, 19(3):90, 2020. [104] Yan-Long Fang and P. Warburton. Minimizing minor embedding energy: an application in quantum annealing. Arxiv preprint arXiv:1905.03291, 05 2019. [105] Davide Venturelli, Salvatore Mandrà, Sergey Knysh, Bryan O’Gorman, Rupak Biswas, and V. Smelyanskiy. Quantum optimization of fully- connected spin glasses. Physical Review X, 5, 06 2014. [106] Vijay V Vazirani. Approximation algorithms. Springer Science & Business Media, 2013. [107] ThomasHCormen, CharlesELeiserson, RonaldLRivest, andCliffordStein. Introduction to algorithms. MIT press, 2009. 90 [108] Michael R. Garey and David S. Johnson. Computers and Intractability: A Guide to the Theory of NP-Completeness. W. H. Freeman & Co., New York, NY, USA, 1979. [109] R Ravi, Madhav V Marathe, SS Ravi, Daniel J Rosenkrantz, and Harry B Hunt III. Many birds with one stone: Multi-objective approximation algo- rithms. In Proceedings of the twenty-fifth annual ACM symposium on Theory of computing, pages 438–447. ACM, 1993. [110] Michel X Goemans. Minimum bounded degree spanning trees. In Founda- tions of Computer Science, 2006. FOCS’06. 47th Annual IEEE Symposium on, pages 273–282. IEEE, 2006. [111] Martin Fürer and Balaji Raghavachari. Approximating the minimum degree spanning tree to within one from the optimal degree. In Proceedings of the third annual ACM-SIAM Symposium on Discrete Algorithms, pages 317–324. Society for Industrial and Applied Mathematics, 1992. [112] Mohit Singh and Lap Chi Lau. Approximating minimum bounded degree spanning trees to within one of optimal. In Proceedings of the thirty-ninth annual ACM symposium on Theory of computing, pages 661–670. ACM, 2007. [113] Jochen Könemann and R Ravi. A matter of degree: Improved approximation algorithms for degree-bounded minimum spanning trees. In Proceedings of the thirty-second annual ACM symposium on Theory of computing, pages 537–546. ACM, 2000. [114] MSZahrani,MartinJLoomes,JAMalcolm,andAndreasAAlbrecht. Alocal search heuristic for bounded-degree minimum spanning trees. Engineering Optimization, 40(12):1115–1135, 2008. [115] Samir Khuller, Balaji Raghavachari, and Neal Young. Low-degree spanning trees of small weight. SIAM Journal on Computing, 25(2):355–368, 1996. [116] Raja Jothi and Balaji Raghavachari. Degree-bounded minimum spanning trees. Discrete Applied Mathematics, 157(5):960–970, 2009. [117] Thang N Bui and Catherine M Zrncic. An ant-based algorithm for finding degree-constrained minimum spanning tree. In Proceedings of the 8th annual conference on Genetic and evolutionary computation, pages 11–18. ACM, 2006. [118] Thang N Bui, Xianghua Deng, and Catherine M Zrncic. An improved ant- based algorithm for the degree-constrained minimum spanning tree problem. IEEE Transactions on Evolutionary Computation, 16(2):266–278, 2011. 91 [119] Mohan Krishnamoorthy, Andreas T Ernst, and Yazid M Sharaiha. Compari- son of algorithms for the degree constrained minimum spanning tree. Journal of heuristics, 7(6):587–611, 2001. [120] AndrewLucas. Isingformulationsofmanynpproblems. Frontiers in Physics, 2:5, 2014. [121] Minsung Kim, Davide Venturelli, and Kyle Jamieson. Leveraging quantum annealing for large mimo processing in centralized radio access networks. In Proceedings of SIGCOMM ’19. ACM, 2019. [122] Jeffrey Marshall, Eleanor G. Rieffel, and Itay Hen. Thermalization, freeze- out, and noise: Deciphering experimental quantum annealers. Phys. Rev. Applied, 8:064025, Dec 2017. 92 Appendices 93 Appendix A Testing a quantum annealer as a quantum thermal sampler A.1 Reverse annealing protocol For our ‘reverse’ anneal protocol, the system is initialized in a classical state chosen at random at s = 1. We then anneal backwards from s = 1 to s = s p at a rate of 0 < ds dt i < 1μs −1 , pause for a time t p and quench to s = 1 at a − 20 − 15 − 10 0 500 1000 1500 2000 2500 3000 t tot (ms) − 140 − 135 − 130 hH IM i QMC GS DW rev Figure A.1: Expectation value for H IM for n = 138 and s p = 0.2. Solid blue curve corresponds to QMC prediction, dashed orange curve corresponds to the Ising GS energy, and data points correspond to the experimental results using the reverse annealing protocols. Error bars correspond to the 95% confidence interval calculated using a bootstrap over the 100 gauges. 94 0 500 1000 1500 2000 2500 3000 t tot (ms) − 138 − 137 − 136 − 135 − 134 − 133 − 132 hH IM i QMC GS DW rev (a) 0 500 1000 1500 2000 2500 3000 t tot (ms) − 138 − 137 − 136 − 135 − 134 − 133 − 132 hH IM i QMC GS DW rev (b) 0 500 1000 1500 2000 2500 3000 t tot (ms) − 138 − 137 − 136 − 135 − 134 − 133 − 132 hH IM i QMC GS DW rev (c) 0 500 1000 1500 2000 2500 3000 t tot (ms) − 138 − 137 − 136 − 135 − 134 − 133 − 132 hH IM i QMC GS DW rev (d) Figure A.2: hH IM i comparison between experiment and theory for n = 138 and s p = 0.5 (a), 0.6 (b), 0.7 (c) and 0.8 (d). Solid blue curve corresponds to QMC prediction, dashed orange curve corresponds to the Ising GS energy, and data points correspond to the DW results using the forward annealing protocols. Error bars correspond to the 95% confidence interval calculated using a bootstrap over the 100 gauges. rate of ds dt f = 1μs −1 . A diagram for this protocol is shown in Fig. 1.4. For each subsequent anneal, we use the configuration returned from the previous anneal as our initial state. Presumably doing so instead of initializing with random state decreases the thermalization time [37]. To allow the system to reach its steady state, we use increasing effective total time (t tot ), defined as the product of the number of anneals n a times the pause 95 time per anneal (t tot =n a t p ). There are upper limits for both the time per anneal (≤ 2 ms) as well as the overall total time (≤ 3000 ms). We choose a long pause time oft p = 1900μs, and only vary the number of anneals n a , in contrast with the forward annealing protocol where n a was fixed and a range of t p explored. Note that the maximum t tot is the same for both protocols. We show the results for s p < s ∗ in Fig. A.1 and the results for s p > s ∗ in Fig. A.2. We see little difference between the forward and reverse anneal protocols for s p <s ∗ , whichisconsistentwiththesystemthermalizingrapidlyinthistransverse- field-dominated regime [82], but we do find some quantitative differences for s p > s ∗ , which appear to become more pronounced away froms p = 0.6. The dependence on the annealing protocol is likely due to two different effects: i) the forward protocol must go through the minimum gap in order to reach s = s p while the reverse protocol never does, and ii) the varying initial condition of each anneal for the reverse anneal. We see that for s p ≤ 0.7, the reverse annealing protocol appears to reach lower average energies than the forward annealing protocol, but fors p = 0.8, the dynamics are likely too slow for it to have reached its steady state value. A.2 Frustrated chains We repeated the study reported on in the main text using frustrated Ising chains. The two models are identical except that while the nonfrustrated model can be mapped to couplers havingJ i,i+1 =−1∀i, in the frustrated case one of the couplings is antiferromagnetic and equals +1. In this case, the GS configuration cannot satisfy all edges; it is 2n-degenerate, with E 0 =−n + 2. We find that for 96 − 20 − 15 − 10 0 400 800 1200 1600 2000 t p (μ s) − 140 − 135 − 130 hH IM i QMC GS DW Figure A.3: At s p = 0.2, the n = 138 frustrated chain exhibits similar behavior to the non-frustrated case. While QMC calculations showhH IM i far from the GS of H IM , experimental results match what we would find at a much later point in the anneal. this model as well the results and conclusions are qualitatively very similar to the non-frustrated ones, withhH IM i very far from the value obtained through QMC before the minimum gap (see Fig. A.3). A.3 Effect of temperature The D-Wave 2000Q device operates at a temperature of approximately 12mK which is the value we used in our QMC simulations. However, temperature is not directly measured during the experiment, and fluctuations could have an effect on our measurements [122]. To eliminate fluctuations in temperature as a possible source of discrepancy, we repeat the simulations for a range of temperatures, and 97 0.2 0.4 0.6 0.8 s p − 140 − 120 − 100 − 80 − 60 − 40 − 20 0 hH IM i T = 25 mK T = 15 mK T = 5 mK T = 1 mK Figure A.4: QMC results for the n = 138 chain at different temperatures are shown, with very slight differences within the range of temperatures at which the annealer could potentially operate. confirm that the change inhH IM i is minimal and does not correspond to the behav- ior we observe experimentally, especially the lowhH IM i value early in the anneal (Fig. A.4). 98 Appendix B Discriminating non-isomorphic graphs with an experimental quantum annealer B.1 Graphs 1 2 3 4 5 6 7 8 9 10 11 12 13 1 2 3 4 5 6 7 8 9 10 11 12 13 Figure B.1: G13 and G13p G13 = [(1, 8), (1, 10), (1, 11), (1, 13), (2, 9), (2, 11), (2, 13), (3, 10), (3, 13), (4, 10), (5, 11), (6, 12), (7, 12), (9, 12), (12, 13)] G13p = [(1, 8), (1, 10), (1, 11), (1, 13), (2, 9), (2, 11), (2, 13), (3, 10), (3, 11), (4, 10), (5, 12), (6, 12), (7, 13), (8, 12), (12, 13)] 99 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Figure B.2: G17 and G17p G17 = [(1, 2), (1, 3), (1, 4), (1, 5), (4, 6), (4, 7), (5, 8), (5, 9), (5, 10), (6, 11), (10, 12), (10, 13), (10, 14), (11, 15), (11, 16), (12, 17)] G17p = [(1, 2), (1, 3), (1, 4), (1, 5), (4, 6), (5, 7), (5, 8), (5, 9), (6, 10), (6, 11), (9, 12), (9, 13), (9, 14), (10, 15), (12, 16), (12, 17)] 100 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 Figure B.3: G25p1, G25p2, G25p3 and G25p4 101 G25p1 = [(1, 6), (1, 7), (2, 6), (2, 14), (2, 22), (3, 7), (4, 8), (4, 9), (5, 8), (5, 9), (6, 9), (10, 14), (10, 15), (11, 15), (12, 16), (12, 17), (13, 16), (13, 17), (14, 17), (18, 22), (18, 23), (19, 23), (20, 24), (20, 25), (21, 24), (21, 25), (22, 25)] G25p2 = [(1, 6), (1, 7), (2, 6), (2, 14), (2, 18), (3, 7), (4, 8), (4, 9), (5, 8), (5, 9), (6, 9), (10, 14), (10, 15), (11, 15), (12, 16), (12, 17), (13, 16), (13, 17), (14, 17), (18, 23), (18, 25), (19, 22), (19, 24), (20, 23), (20, 25), (21, 24), (21, 25)] G25p3 = [(1, 6), (1, 7), (2, 6), (2, 10), (2, 18), (3, 7), (4, 8), (4, 9), (5, 8), (5, 9), (6, 9), (10, 15), (10, 17), (11, 14), (11, 16), (12, 15), (12, 17), (13, 16), (13, 17), (18, 23), (18, 25), (19, 22), (19, 24), (20, 23), (20, 25), (21, 24), (21, 25)] G25p4 = [(1, 7), (1, 9), (2, 6), (2, 8), (3, 7), (3, 9), (4, 8), (4, 9), (5, 1), (5, 10), (5, 18), (10, 15), (10, 17), (11, 14), (11, 16), (12, 15), (12, 17), (13, 16), (13, 17), (18, 23), (18, 25), (19, 22), (19, 24), (20, 23), (20, 25), (21, 24), (21, 25)] 102 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 Figure B.4: G27 and G27p G27 = [(1, 14), (1, 17), (2, 14), (2, 22), (3, 4), (3, 5), (4, 10), (4, 12), (5, 11), (5, 13), (6, 7), (6, 8), (6, 15), (7, 10), (7, 11), (8, 12), (8, 13), (9, 12), (9, 13), (9, 14), (10, 15), (11, 15), (14, 15), (16, 17), (16, 21), (17, 18), (18, 19), (19, 20), (20, 21), (22, 23), (22, 27), (23, 24), (24, 25), (25, 26), (26, 27)] G27p = [(1, 14), (1, 17), (2, 14), (2, 23), (3, 4), (3, 5), (4, 10), (4, 11), (5, 12), (5, 13), (6, 7), (6, 8), (6, 15), (7, 10), (7, 12), (8, 11), (8, 13), (9, 12), (9, 13), (9, 14), (10, 15), (11, 15), (14, 15), (16, 17), (16, 21), (17, 18), (18, 19), (19, 20), (20, 21), (22, 23), (22, 27), (23, 24), (24, 25), (25, 26), (26, 27)] 103 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 Figure B.5: G33, G33p, G33p1, G33p2 104 G33 = [(1, 6), (1, 7), (2, 6), (2, 14), (2, 22), (2, 30), (3, 7), (4, 8), (4, 9), (5, 8), (5, 9), (6, 9), (10, 14), (10, 15), (11, 15), (12, 16), (12, 17), (13, 16), (13, 17), (14, 17), (18, 22), (18, 23), (19, 23), (20, 24), (20, 25), (21, 24), (21, 25), (22, 25), (26, 30), (26, 31), (27, 31), (28, 32), (28, 33), (29, 32), (29, 33), (30, 33)] G33p = [(1, 6), (1, 7), (2, 6), (2, 14), (2, 22), (2, 26), (3, 7), (4, 8), (4, 9), (5, 8), (5, 9), (6, 9), (10, 14), (10, 15), (11, 15), (12, 16), (12, 17), (13, 16), (13, 17), (14, 17), (18, 22), (18, 23), (19, 23), (20, 24), (20, 25), (21, 24), (21, 25), (22, 25), (26, 31), (26, 33), (27, 30), (27, 32), (28, 31), (28, 33), (29, 32), (29, 33)] G33p1 = [(1, 6), (1, 7), (2, 6), (2, 10), (2, 18), (2, 26), (3, 7), (4, 8), (4, 9), (5, 8), (5, 9), (6, 9), (10, 15), (10, 17), (11, 14), (11, 16), (12, 15), (12, 17), (13, 16), (13, 17), (18, 23), (18, 25), (19, 22), (19, 24), (20, 23), (20, 25), (21, 24), (21, 25), (26, 31), (26, 33), (27, 30), (27, 32), (28, 31), (28, 33), (29, 32), (29, 33)] G33p2 = [(1, 6), (1, 7), (2, 6), (3, 7), (4, 8), (4, 9), (5, 8), (5, 9), (5, 10), (5, 18), (5, 26), (6, 9), (10, 15), (10, 17), (11, 14), (11, 16), (12, 15), (12, 17), (13, 16), (13, 17), (18, 23), (18, 25), (19, 22), (19, 24), (20, 23), (20, 25), (21, 24), (21, 25), (26, 31), (26, 33), (27, 30), (27, 32), (28, 31), (28, 33), (29, 32), (29, 33)] 105 Appendix C Ferromagnetically shifting the power of pausing C.1 Level-based mapping for BD-MST problems Consider a graph G = (V,E) with weights w(E) for each edge, from which we wish to obtain a minimal weighted spanning tree with maximum degree Δ. This will involve minimizing the sum of the weights of the tree edges, represented by the cost function C 0 = X p,v w pv x u,v (C.1) Several constraints will also be imposed to ensure that the graph is in fact a spanning tree and its degree is bounded by Δ. A root for the tree will be picked randomly or based on problem structure— generally, picking a high-degree vertex as the root will result in lower resource costs—and assigned to level 1. Its children will be at level 2, their children at level 3 and so on, leading to the ‘level-based’ designation. The variables x p,v appear in Eq. (C.1) represent the parent-child relationships in the tree; x p,v = 1 if p is the parent of v (and 0 if not). The indices p,v range over p = 1,...,n and v = 2,...,n, restricted to (intersected with) pairs (p,v) or (v,p) that occur in E. Thus there are 2 variables for every edge not containing 106 the root, and one for every root edge, giving 2m−d r total x p,v variables, with m being the number of edges in E and d r the degree of the root. Since our problem needs to be in QUBO form, the constraints will be expressed as penalty terms. The first penalty term enforces that every node (except the root) has exactly one parent, C pen1 = X v∈{2,...,n}   X p:(pv)∈E x p,v − 1   2 . (C.2) The number of terms in the sum is 2m−d r , i.e. equal to the number of variables x p,v . The second penalty term enforces that each vertex exists at exactly one level in the tree, C pen2 = X v∈{2,...,n} n X `=2 y v,` − 1 ! 2 . (C.3) It introduces the y v,` variables, with y v,` = 1 if v is at depth ` of the tree, v = 2,...,n, ` = 2,...,n. There are (n− 1) 2 such variables. However, since the number of variables will eventually determine how many logical qubits the problem requires, it is in our interest to reduce it as much as possible. By picking the root smartly the range of` can be reduced. We also carry out the following pre- processing: taking the original graph G = (V,E), the distance from each node to the one we have selected as the tree root is calculated. Given that it is impossible for a node to be at a level smaller than its distance to the root, we can avoid generating any y v,` for which that is the case, further bringing down the total number of y v,` variables. 107 The third penalty term enforces that the tree has degree at most Δ, C pen3 = v X p=2   X v:(pv)∈E x p,v − Δ−1 X j=1 z p,j   2 + ( X v:(1v)∈E x 1,v − Δ X j=1 z 1,j ) 2 . (C.4) It is separated into two terms to account for the fact that the root can have up to Δ children, while all other nodes cannot have more than (Δ− 1), since they have a parent. To enforce the inequality P v:(pv)∈E x p,v ≤ Δ− 1, the integer variable z p ∈ [0, Δ− 1] is introduced as a slack variable, and we can then enforce the inequality as P v:(pv)∈E x p,v = z p . The integer variable is further encoded into binary variablesz p,j . We use unary encoding, which can be applied straightforward to arbitrary values of Δ. The fourth and final penalty term enforces that the tree encoding is consistent, i.e., that if p is the parent of v then its level is one less than v’s, C pen4 = X p,v n X `=3 x p,v y v,` (1−y p,`−1 ) + dr X v=2 x 1,v (1−y v,2 ) + dr X v=2 y v,2 (1−x 1,v ). (C.5) The last two sums handle the edges connected to the root and their terms are quadratic, while the first sum deals with the remaining edges and produces cubic terms of the form x p,v y v,` (1−y p,`−1 ). While the original number of cubic terms would be (2m− 2d r )· (n− 2), thanks to the pre-processing of they v,` variables this number is reduced. Because cubic terms cannot be directly encoded in D-Wave, we must use ancilla variablesa p,v,` =x p,v y v,` . Substitutinga≡a p,v,` into the cubic terms we obtain quadratic ones, 4a−ay p,`−1 +x p,v y v,` − 2ax p,v − 2ay v,` . (C.6) 108 The total number of variables (and hence, logical qubits) without pre-processing is: 2m−d r + (n− 1) 2 +n(Δ− 1) + 1 + (2m− 2d r )· (n− 2)' 2mn +n 2 . (C.7) This would mean, for instance, that the complete graph K 5 with Δ = 3 would require between 86 and 100 logical qubits (depending ond r ). With pre-processing, we are able to bring this number down to 74. Finally, we can write the overall objective Hamiltonian as H p =C 0 +A(C pen1 +C pen2 +C pen3 +C pen4 ), (C.8) where we have defined the minimum penalty weight A to be larger than the maximum edge weight to guarantee that the optimal solution to the mapped prob- lem is in fact the desired optimal BD-MST: A =w max + = max (uv)∈E w uv +. (C.9) This is a sufficient condition∀ > 0. A solution with extra edges will not be returned, since we are minimizing the cost function and having an extra edge would always cost more compared to a tree. There is the possibility of a solution trying to attain a lower cost by removing an edge from a tree. But because A > w max , even if the BD-MST contained the maximum edge weight, and we removed that edge, the penalty would lead to a greater cost than leaving that edge. And more so if we removed any other edge instead. In our runs, we set = 0 (which means there could be a potential degeneracy for the tree containing the maximum weight edge and the non-tree 109 with this edge removed), and the solutions returned from the quantum annealer are checked to ensure they are in fact trees. 110 C.2 BD-MST problem instances All connected graphs of n=5, withm =|E| ranging from 4 to 10 are considered where an BD-MSTwith Δ = 2 exists. The edges in these graphs is provide in Table C.1. Additionally, the graph labeled ’m5ver5’ is included to demonstrate the BD- MSTwith Δ≥ 3. m label name edges 4 m4ver1 DhC (1,2), (2,3), (3,4), (4,5) 5 m5ver1 Dhc (1,2), (2,3), (3,4), (4,5), (1,5) 5 m5ver2 DiK (1,2), (2,3), (2,5), (3,4), (4,5) 5 m5ver3 DjC (1,2), (2,3), (2,4), (3,4), (4,5) 5 m5ver5 Dsc (1, 2),(1, 3), (1,4),(1,5),(4,5) 5 m5ver6 DKs (1, 2), (2, 3), (3, 4), (4, 5), (3, 5) 6 m6ver1 DyK (1,2), (1,5), (2,5), (2,3), (3,4), (4,5) 6 m6ver2 DjS (1,2), (2,3), (2,4), (2,5), (3,4), (4,5) 6 m6ver3 DjK (1,2), (2,3), (2,5), (3,5), (3,4), (4,5) 6 m6ver4 D{K (1,2), (1,5), (1,3), (2,3), (3,4), (4,5) 6 m6ver5 D{c (1, 2), (1, 3), (2, 3), (3, 5), (3, 4), (4, 5) 6 m6ver6 D]o (1, 2), (2, 3), (3, 4), (1, 4), (2, 5), (4, 5) 7 m7ver1 D|S (1,2), (1,5), (1,4), (2,5), (2,3), (3,4), (4,5) 7 m7ver2 DzW (1,2), (1,5), (2,5), (2,3), (2,4), (3,5), (4,5) 7 m7ver3 D|c (1,2), (1,3), (1,4), (1,5), (2,3), (3,4), (4,5) 7 m7ver4 D∼C (1,2), (1,3), (1,4), (2,4), (2,3), (3,4), (4,5) 7 m7ver5 D]w (1, 2), (2, 3), (3, 4), (1, 4), (4, 5), (2, 5), (3, 5) 7 m7ver6 Dh{ (1, 2), (1, 3), (1, 4), (1, 5), (2, 3), (3, 4), (4, 5) 8 m8ver1 D}k (1,2), (1,5), (1,3), (1,4), (2,5), (2,3), (3,4), (4,5) 8 m8ver2 Dz[ (1,2), (1,5), (2,5), (2,3), (2,4), (3,4), (3,5), (4,5) 9 m9ver1 D∼k (1, 2), (2, 3), (4, 5), (1, 5), (1, 4), (1, 3), (2, 5), (2, 4), (3, 5) 10 m10ver1 D∼{ (1, 2), (2, 3), (3, 4), (4, 5), (1, 5), (1, 4), (1, 3), (2, 5), (2, 4), (3, 5) Table C.1: n=5 graphs. For each graph, problem instances were generated by assigning a set of weights by sampling from one of the lists of weights appearing in Table C.2. The first m weights are used if the list is longer than the number of edges. 111 label weight list w2 [1, 2, 1, 2, 1, 2, 1, 2, 1, 2] w3 [1, 1, 2, 1, 1, 2, 1, 1, 2, 1] w4 [1, 1, 2, 2, 1, 1, 2, 2, 1, 1] w5 [1, 4, 1, 4, 1, 4, 1, 4, 1, 4] w6 [1, 3, 6, 1, 3, 6, 1, 3, 6, 1] w7 [1, 7, 1, 7, 1, 7, 1, 7, 1, 7] w8 [3, 2, 1, 3, 2, 1, 3, 2, 1, 3] w9 [4, 3, 2, 1, 4, 3, 2, 1, 4, 3] w10 [5, 4, 3, 2, 1, 5, 4, 3, 2, 1] w11 [6, 5, 4, 3, 2, 1, 6, 5, 4, 3] w12 [7, 6, 5, 4, 3, 2, 1, 7, 6, 5] w13 [1, 1, 3, 4, 2, 1, 2, 3, 4, 2] w14 [3, 2, 1, 1, 1, 1, 2, 4, 2, 2] w15 [2, 1, 2, 1, 4, 1, 1, 3, 3, 2] w16 [4, 3, 3, 4, 3, 3, 4, 3, 4 ] w17 [3, 4, 7, 5, 5, 5, 5] w18 [2, 1, 4, 1, 2, 1, 2] w19 [4, 6, 4, 7, 4, 7] w20 [1, 1, 2, 3, 2, 3] w21 [4, 5, 4, 5, 5] w22 [2, 2, 6, 2, 4] w23 [3, 3, 5, 2, 3, 2, 5, 2, 5] w24 [4, 3, 2, 2] w25 [2, 2, 6, 2, 4] w26 [4, 3, 3, 3] w27 [3, 4, 7, 5, 5, 5, 5] w28 [4, 6, 4, 7, 4, 7] w29 [6, 4, 2, 2] Table C.2: Graph weights used to generate different instances. 112 C.3 Embedding statistics Table C.3 contains mapped and embedded problem size for each graph on the current Chimera architecture. Analogous data for a future D-Wave architecture (Pegasus) is shown in table C.4. n m # logical variables # physical variables 4 6 35 108-150 4 5 29 65-121 4 4 25 60-116 4 4 23 47-82 4 3 20 40-76 5 4 32 83-140 5 5 42 121-215 5 5 43 138-169 5 5 44 149-205 5 5 47 151-220 5 5 39 112-179 5 6 50 169-205 5 6 53 194-255 5 6 49 181-246 5 6 46 148-272 5 6 50 170-249 5 6 50 164-217 5 7 54 193-247 5 7 58 229-300 5 7 50 171-260 5 7 55 226-281 5 7 56 227-273 5 7 50 162-224 5 8 58 219-284 5 8 64 287-362 5 9 66 299-413 5 10 74 380-485 6 15 137 1166-1293 Table C.3: Mapped problem size for n = 4− 6 using the current Chimera archi- tecture. 113 n m # logical variables # physical variables 4 6 35 54-71 5 10 74 164-207 6 15 137 1166-1293 7 21 230 1018-1291 8 28 359 2046-2712 9 36 530 3744-4454 10 45 749 6024-7889 Table C.4: Mapped problem size forn = 4− 10 using the future Pegasus architec- ture. The new Pegasus architecture will more than double the number of qubits with respect to the current D-Wave chip, and the connectivity will be greatly increased. Fig. C.1 shows the median embedded size for the fully connected graph of size n. The Pegasus architecture not only allows going up to n = 10 compared to Chimera’s n = 6, but also performs better than Chimera (smaller number of physical variables) for n = 4− 6. 114 0 200 400 600 # logical variables 10 2 10 3 10 4 # physical variables Chimera Pegasus Figure C.1: Comparison of embedded size between Chimera and Pegasus architec- tures for complete graphs withn between 4 and 10. Chimera embedding performed with D-Wave’s SAPI2 find_embedding routine with the D-Wave 2000Q hardware adjacency graph. Pegasus embedding performed with the new Ocean minorminer find_embedding routine. Median physical qubits as a function of the number of logical qubits. 115

Abstract (if available)

Abstract Experimental quantum annealers have made great progress in the last decade. In recent years, some of the technological improvements have led to the possibility of modifying the annealing schedules, with new features such as pausing and reverse annealing. We explore potential uses for these novel capabilities, and the physics that underlies the behavior we observe when they are deployed, through three different projects. In the first one, we demonstrate that using the currently implementable schedules one cannot in general reliably sample from the thermal state associated with a target quantum Hamiltonian. In the second one, we consider an application for which this limitation is not an obstacle, namely solving certain instances of the graph isomorphism problem. Finally, in the third one, we show how pausing can improve the performance of the quantum annealer on a class of optimization problems that are not native to the device’s architecture.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Error suppression in quantum annealing

PDF

Open-system modeling of quantum annealing: theory and applications

PDF

Explorations in the use of quantum annealers for machine learning

PDF

Error correction and quantumness testing of quantum annealing devices

PDF

The theory and practice of benchmarking quantum annealers

PDF

Dissipation as a resource for quantum information processing: harnessing the power of open quantum systems

PDF

Theory and simulation of Hamiltonian open quantum systems

PDF

Tunneling, cascades, and semiclassical methods in analog quantum optimization

PDF

Dynamic topology reconfiguration of Boltzmann machines on quantum annealers

PDF

Advancing the state of the art in quantum many-body physics simulations: Permutation Matrix Representation Quantum Monte Carlo and its Applications

PDF

Minor embedding for experimental investigation of a quantum annealer

PDF

Imposing classical symmetries on quantum operators with applications to optimization

PDF

Exploiting structure in the Boolean weighted constraint satisfaction problem: a constraint composite graph-based approach

PDF

Destructive decomposition of quantum measurements and continuous error detection and suppression using two-body local interactions

PDF

Optimization of the combinatoric closely spaced objects resolution algorithm with adiabatic quantum annealing

PDF

Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras

PDF

Demonstration of error suppression and algorithmic quantum speedup on noisy-intermediate scale quantum computers

PDF

Quantum information flow in steganography and the post Markovian master equation

PDF

Trainability, dynamics, and applications of quantum neural networks

PDF

DG structures on categorified quantum groups

Asset Metadata

Creator Gonzalez Izquierdo, Zoe (author)

Core Title Understanding physical quantum annealers: an investigation of novel schedules

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 05/07/2020

Defense Date 03/19/2020

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag OAI-PMH Harvest,quantum annealing,quantum computation

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Hen, Itay (committee chair)

Creator Email zgonzale@usc.edu,zoe.gonzalez.izquierdo@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-301138

Unique identifier UC11663616

Identifier etd-GonzalezIz-8457.pdf (filename),usctheses-c89-301138 (legacy record id)

Legacy Identifier etd-GonzalezIz-8457.pdf

Dmrecord 301138

Document Type Dissertation

Rights Gonzalez Izquierdo, Zoe

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA