Scaling Advantage in Approximate Optimization with Quantum Annealing

The pursuit of a quantum speedup using quantum processors is a major theme in modern physics and computer science. Two of the leading application areas are discrete optimization and quantum simulation. In the latter context, impressive recent quantum annealing (QA) advances have been made for fast, coherent anneals lasting on the order of the single superconducting flux qubit coherence time [20, 21]. While this diabatic approach is considered promIsing [22], it cannot be expected to scale up without the introduction of error correction or suppression, as decoherence and control errors pose insurmountable challenges to Hamiltonian quantum computation models [23, 24, 25, 26, 27, 28], just as they do for gate-model quantum computers. In the absence of a fault-tolerance threshold theorem [29] for QA, a variety of Hamiltonian error suppression techniques have been proposed and analyzed as ways to reduce the error rates of this computational model and the closely related model of adiabatic quantum computation [30, 31, 32, 33, 34, 35, 36], providing tools towards scalability.

However, despite these theoretical advances, there are currently practical limitations to the types and locality of programmable interactions in the Hamiltonians of quantum annealing hardware, even as new devices have started to emerge [37, 38, 39, 16]. To address these limitations, quantum annealing correction (QAC) [18] was developed as a realizable error suppression method for quantum annealing, targeting the available and restricted set of control operations in quantum annealers. QAC has been demonstrated to enhance the success probability of quantum annealing and mitigate the analog control problem  [18, 40, 41]. The QAC method is based on a repetition-code encoding of the Hamiltonian and does not fully realize a Hamiltonian stabilizer code. Despite this, it has been shown theoretically to increase the energy gap of the encoded Hamiltonian and reduce tunneling barriers, thus softening the onset of the associated critical dynamics as well as lowering the effective temperature [42].

Here, departing from the traditional paradigm of using QA for exact optimization, we demonstrate—by incorporating QAC—the first genuine scaling advantage in approximate optimization (low-energy sampling) using a quantum annealer. Even approximate optimization can be computationally hard unless $\mathrm{P=NP}$ [43, 44], so we do not expect the advantage we exhibit to amount to more than a polynomial speedup. However, whereas the scaling advantages reported in previous work were relative to simulated annealing [13, 16], the advantage we find here is over the best currently available general heuristic classical optimization method: parallel tempering with isoenergetic cluster moves (PT-ICM) [17]. This result is enabled by implementing QAC on the D-Wave Advantage quantum annealer for the Sidon-set spin glass instance class [9], embedded on the logical graph formed after the encoding step. The advantage of quantum annealing over PT-ICM is diminished without QAC, thus highlighting the crucial role of quantum error suppression.

Problems on the logical QAC graph can also be encoded directly by setting $J_{p}=0$, resulting in three uncoupled and unprotected, parallel classical copies of the problem instance. We then extract the energies of all three copies as independent annealing samples and denote this “unprotected” QA method by U3. We use the U3 method below to test whether QAC has a genuine advantage over simple classical repetition coding.

We define $[\operatorname{\text{TT}\varepsilon}]_{q}$ of an instance class as the $q$-th quantile of $\operatorname{\text{TT}\varepsilon}$ over the entire instance class. Here, we focus only on the median quantile, $q=0.5$, denoted $[\operatorname{\text{TT}\varepsilon}]_{\mathrm{Med}}$. For a given disorder, instance size, and $\varepsilon$-target, we find the annealing time $t_{f}$ (and penalty strength for QAC) that minimizes $[\operatorname{\text{TT}\varepsilon}]_{\mathrm{Med}}$. We restrict the penalty coupling strengths to the set $J_{p}\in\{0.1,0.2,0.3\}$ to reduce resource requirements for parameter optimization, as $J_{p}=0.2$ is the penalty strength that most frequently optimizes the success probabilities of individual instances, and the dependence on $J_{p}$ above $0.2$ becomes weak.

The shortest available annealing time on the D-Wave Advantage QPU accessed via Leap is $t_{f}^{\min}=0.5\,\mathrm{\mu s}$, and the bottom panels of Fig. 2 show that as the target residual energy density is increased, progressively larger problem sizes are needed to ensure that $t_{f}^{\text{opt}}(N)>t_{f}^{\min}$. We cannot rule out that with access to lower annealing times, one would find $t_{f}^{\text{opt}}(N)<t_{f}^{\min}$ for all $N$ values where we empirically find $t_{f}^{\text{opt}}(N)=t_{f}^{\min}$. We thus formulate a null hypothesis for each $N$ that $t_{f}^{\text{opt}}(N)\leq t_{f}^{\min}$ and compute a $P$ value as the empirical number of bootstrap samples whose $t_{f}^{\text{opt}}(N)=t_{f}^{\min}$, out of a total of $200$ samples (see Methods for details of our statistical analysis). To compute the $\operatorname{\text{TT}\varepsilon}$ scaling, i.e., the slope $\alpha$ in a fit to $\operatorname{\text{TT}\varepsilon}=cN^{\alpha}$, for each $\varepsilon$ we use only those $t_{f}^{\text{opt}}(N)$ values whose $P<0.05$ (filled circles in Fig. 2). We can thus be confident that the reported slopes reflect the true scaling of U3 and QAC.

Our first observation is that the QAC scaling is always better than the U3 scaling, which is consistent with previous studies concerning the effect of analog coupling errors (“$J$-chaos”) on the TTS for spin-glass instances [41]. Such errors are expected for the S28 instances due to the relatively high precision their specification requires.

Second, we observe that U3 and QAC reduce the absolute algorithmic runtime by four orders of magnitude compared to PT-ICM. However, this is not a scaling advantage, and since our PT-ICM calculations could be sped up by employing faster classical processors, we do not consider this a robust finding. Similarly, we exclude the programming and readout time used by the D-Wave QPU. The primary source of overhead is the readout time per sample, which scales with problem size and can reach $200\,\mathrm{\mu s}$ per sample for this QPU. This timing varies by hardware generation, and its inclusion obfuscates the dominant scaling source.

Third, and most significantly, we observe that as the target optimality gap increases QA’s scaling overtakes PT-ICM. Notably, at $\varepsilon=1\%$, QAC exhibits a scaling exponent of $1.69\pm 0.12$, compared to PT-ICM’s $1.93\pm 0.03$, and at $\varepsilon=1.25\%$, QAC and U3 exhibit scaling exponents of $1.15\pm 0.22$ and $1.76\pm 0.06$, respectively, compared to PT-ICM’s $1.87\pm 0.02$. The scaling of U3 at optimal annealing times continues to decrease as $\varepsilon$ increases; at $\varepsilon=1.5\%$ (not shown), the scaling exponents for U3 and PT-ICM become $1.60\pm 0.07$ and $1.86\pm 0.04$, respectively. This is robust evidence of a quantum annealing scaling advantage over the best available classical heuristic optimization algorithm.

We are unable to determine the scaling of QAC for $\varepsilon>1.25\%$, as we cannot confirm that $t_{f}^{\text{opt}}(N)>t_{f}^{\min}$ for any $N$; as can be seen in Fig. 2, already for $\varepsilon=1.25\%$ only the largest three $N$ values satisfy the $P<0.05$ criterion. However, given the consistently better scaling of QAC for lower values of $\varepsilon$, where $t_{f}^{\text{opt}}(N)>t_{f}^{\min}$ for QAC over a range of problem sizes, it is reasonable to conclude that the QAC scaling would be a further improvement over U3 if its true $t_{f}^{\text{opt}}(N)$ could be established for $\varepsilon>1.25\%$; this would require access to shorter annealing times or larger system sizes.

We note that it is unsurprIsing that, given a sufficiently large target optimality gap, the D-Wave QPU returns sample energies within that gap. Similarly, we can expect that for large enough $\varepsilon$, even simulated annealing or greedy descent will be nearly guaranteed success in polynomial time. The significance of our result is that QA reaches near-linear scaling at a smaller optimality gap target than PT-ICM. Thus, we refer to this result as an approximate optimization advantage for quantum annealing. We also note that Ref. [21] similarly reported a QA optimization advantage for the residual energy density (for 3D spin glasses), but this was done at a fixed problem size (of $N=5374$ physical qubits), and was instead concerned with the convergence of the residual energy density with the annealing time. We reemphasize that, in contrast, we are reporting a scaling advantage as a function of problem size, the proper context for quantum speedup claims. We explain in Methods that our results stand regardless of whether an optimality gap or residual energy density target is chosen for benchmarking.

We performed a finite-size scaling analysis and data collapse, and the collapsed Binder cumulant for the S28 instances is shown in Fig. 3. We find that $\mu_{\text{QAC}}=4.81\pm 0.22$ (at a penalty strength of $0.1$) compared to $\mu_{\text{U3}}=7.53\pm 0.47$. The reduction of $\mu$ is lost when $\lambda=0.2$ (not shown), suggesting that $\lambda=0.1$ is optimal in the sense of diabatic error suppression.

This effect indicates that QAC is much more effective at suppressing diabatic excitations. I.e., at equal annealing times, the dynamics are more adiabatic under QAC, in agreement with theory [42]. The significant reduction in $\mu$ suggests that in addition to $J$-chaos suppression, diabatic error suppression by QAC is responsible for the improved $\operatorname{\text{TT}\varepsilon}$ and shorter optimal annealing times. Additional context is provided in Methods.

There are a few limitations to the scope of our conclusion. First, our results do not imply that finding states within small, constant gaps (and indeed finding the ground state itself) is easy for quantum annealing, nor do they imply that all spin glass problems are amenable to an approximate optimization scaling advantage via quantum annealing. Second, being finite-range and two-dimensional limits the range of applications of the problem family we have studied here. To achieve an algorithmic quantum advantage in an application setting, the next challenge for quantum optimization is demonstrating a hardware-scalable advantage in densely connected problems at sufficiently small optimality gaps.

The scaling of PT-ICM is ideally evaluated using the parameters that best optimize the $\operatorname{\text{TT}\varepsilon}$ for each disorder realization, instance size, and target $\varepsilon$. However, a rigorous optimization of the number and choice of replica temperatures for all target $\varepsilon$’s and system sizes is computationally infeasible. To ensure our results hold for any choice of reasonably optimized temperatures, we repeated the $\operatorname{\text{TT}\varepsilon}$ evaluation with four temperature sets summarized in Table 1, which includes both logarithmically-spaced and feedback-optimized [1] temperatures. The $\operatorname{\text{TT}\varepsilon}$ for a given disorder, instance size, and energy target was chosen from the best $\operatorname{\text{TT}\varepsilon}$ out of the four temperature sets, illustrated in Fig. 4. At the optimality gap targets of interest, most temperature sets’ differences are not appreciable and are unlikely to affect our conclusions. A more comprehensive but computationally expensive optimization of PT-ICM would involve a grid search over a range of the parameters $\beta_{\mathrm{min}}$, $\beta_{\mathrm{max}}$, and $N_{\mathrm{icm}}$.

We also implemented and considered the performance of simulated annealing but do not present its results here as this algorithm was not competitive at large problem sizes. Our PT-ICM implementation is a general-purpose solver written in the Rust programming language, targeting CPU execution, and accepts instances with any connectivity. This is in contrast to previous studies that have used simulated annealing or the Hamze-de Freitas-Selby algorithm solvers that are specialized for problems defined on the Chimera graph [4, 5].

Next, we formulate a null hypothesis and compute $P$ values as follows. The null hypothesis is that the optimal annealing time $t_{f}$ is the minimum accessible $t_{f}^{\min}=0.5\,\mathrm{\mu s}$, i.e., that the true optimal annealing time is not above $t_{f}^{\min}$. The $P$ value is the empirical number of bootstrap samples whose optimal $t_{f}$ was $0.5\,\mathrm{\mu s}$, out of 200 samples. Filled circles in Fig. 2 mean $P<0.05$ for the probability that $t_{f}=0.5\mu s$ in the bootstrap sample, while open circles mean $P<0.20$. The filled circles show which points have the highest confidence that the optimal annealing time is not below $0.5\,\mathrm{\mu s}$.

While the collapse seen in Fig. 3 is visually convincing, the error estimates for $U$ are, unfortunately, too large to determine the KZ exponent or to extract $\nu$ and $z$. In addition, we do not observe that the extracted KZ exponent is universal: the estimate for $\mu$ that best collapses the Binder cumulant for binomial (as opposed to S28) spin glass instances is $\mu\approx 9\pm 1$ for U3, and is reduced to $\mu\approx 7.5\pm 0.6$ with QAC at $\lambda=0.1$. Thus, while the scaling ansatz yields a useful complementary way to quantify the advantage of QAC in terms of sampled quantities in addition to the $\operatorname{\text{TT}\varepsilon}$ metric, the lack of universality and the possibility that the spin-glass transition temperature is zero [6] do not clearly support a universal critical description of the annealing dynamics, quite unlike the conclusions of Refs. [7, 8].

This algorithm optimizes bulk energies throughout the volume of the spin-glass at the expense of large energy violations over regions scaling with the surface area of the spin-glass, i.e., at the patch boundaries. The boundary excitations become negligible for sufficiently large sizes compared to the bulk energies, with the latter eventually reaching the desired residual energy. More precisely, up to $4L_{0}$ boundary couplings may be violated per volume of $L_{0}^{2}$, so the residual energy density for this algorithm is upper-bounded by $4/L_{0}$. Thus, we estimate that the regime of system sizes where this algorithm applies for the $\rho$ targets we examine, e.g., 1%, would require a reliable, efficient solver for instances with at least a sub-problem side length of $L_{0}\approx 400$, resulting in a prefactor of $2^{400}\sim 10^{120}$, which is entirely impractical. Nevertheless, this algorithm could be a starting point for parallel and quantum-classical hybrid algorithms for massive, finite-dimensional problems. Furthermore, such an algorithm could reach a similar target $\varepsilon$ in linear time with increasing probability as the system size increases due to the decreasing variance of $E_{0}/(JN)$.

In an attempt to estimate the influence of precision and ground state degeneracy in the S28 disorder, we also studied range 6 (R6) instances for which the couplings were randomly drawn from $J\in\{\pm 1/6,\pm 2/6,\pm 3/6,\pm 4/6,\pm 5/6,\pm 1\}$. The main contrast between R6 and S28 is that the minimum non-zero local field a qubit may experience is much smaller in S28 disorder. Furthermore, in the absence of the Sidon property, the R6 disorder is more susceptible to ground-state degeneracies due to floppy qubits. However, it is less likely to occur than in a smaller range disorder case such as range 3.

The $\operatorname{\text{TT}\varepsilon}$ scaling for the R6 case is shown in Fig. 6. For a smaller $\varepsilon=0.75\%$, the scaling of QAC is better with R6 disorder than S28 disorder (Fig. 2), going from 2.26 to 2.02, then equalizes as $\varepsilon$ increases. Perhaps surprIsingly, the opposite is true of U3, which scales worse on R6 disorder at smaller epsilon (2.92 to 3.10) before scaling better than S28 disorder above $\varepsilon\approx 1.1\%$ (2.02 to 1.82). That is, while S28 is overall more challenging than R6 under QAC, low energy sampling of R6 is more challenging for unprotected QA than S28. This is likely caused by the worse relative precision of small coupling strengths in unprotected QA, which would not affect the S28 disorder.