Probing for quantum speedup in spin-glass problems with planted solutions

The availability of quantum annealing devices with hundreds of qubits has made the experimental demonstration of a quantum speedup for optimization problems a coveted, albeit elusive goal. Going beyond earlier studies of random Ising problems, here we introduce a method to construct a set of frustrated Ising-model optimization problems with tunable hardness. We study the performance of a D-Wave Two device (DW2) with up to 503 qubits on these problems and compare it to a suite of classical algorithms, including a highly optimized algorithm designed to compete directly with the DW2. The problems are generated around predetermined ground-state configurations, called planted solutions, which makes them particularly suitable for benchmarking purposes. The problem set exhibits properties familiar from constraint satisfaction (SAT) problems, such as a peak in the typical hardness of the problems, determined by a tunable clause density parameter. We bound the hardness regime where the DW2 device either does not or might exhibit a quantum speedup for our problem set. While we do not find evidence for a speedup for the hardest and most frustrated problems in our problem set, we cannot rule out that a speedup might exist for some of the easier, less frustrated problems. Our empirical findings pertain to the specific D-Wave processor and problem set we studied and leave open the possibility that future processors might exhibit a quantum speedup on the same problem set.

Figure 1

Examples of randomly generated loops and couplings on the DW2 Chimera graph. (Top) Qubits and couplings participating in the loops are highlighted in green and purple, respectively. Only even-length loops are embeddable on the Chimera graph. (Bottom) Distribution of $J$ values for a sample problem instance with $N=126$ spins and edges and 101 loops. It is virtually impossible to recover the loop Hamiltonians $H_j$ from a given $H_{\mathrm{Ising}}$. The couplings are all eventually rescaled to lie in $[-1,1]$. We always set the local fields $h_i$ to zero as nonzero fields tended to make the problems easier.

Figure 2

Time to solution as a function of clause density. Shown is $\mathrm{TTS}(L,\alpha ,0.5)$ (log scale) for (a) DW2 and (b) HFS as a function of the clause density. The different colors represent the different Chimera subgraph sizes, which continue to $L=12$ in the HFS case. In both cases there is a clear peak. From the HFS results we can identify the peak position as being at $\alpha =0.17\pm 0.01$, which is consistent with the peak position in the DW2 results. Error bars represent $2\sigma$ confidence intervals.

Figure 3

Frustration fraction. Shown is the fraction of frustrated couplings (the number of frustrated couplings divided by the total number of couplings, where a frustrated coupling is defined with respect to the planted solution) as a function of clause density for different Chimera subgraphs $C_L$, in the case of loops of length $\ge 8$, averaged over the 100 instances for each given $\alpha$ and $N$. There is a broad peak at $\alpha \approx 0.25$. This is the clause density at which there is the largest fraction of frustrated couplings; it is near where we expect the hardest instances to occur, in good agreement with Fig. 2.

Figure 4

Sketch of the relation between the $ln\left(\mathrm{TTS}\right)$ curves for optimal and suboptimal annealing times. The annealing time needs to be optimized for each problem size. Blue represents the TTS with a size-independent annealing time $t_a$. Red represents the optimal TTS corresponding to having an optimal size-dependent annealing time $t_a^{\mathrm{opt}}\left(L\right)$, i.e., the lower envelope of the full series of fixed annealing time TTS curves. This curve need not be linear as depicted, though we expect it to be linear for NP-hard problems. The blue line upper bounds the red line since by definition ${\mathrm{TTS}}_{\mathrm{DW}2}(\lambda ,t_a^{\mathrm{opt}}\left(L\right))\le {\mathrm{TTS}}_{\mathrm{DW}2}(\lambda ,t_a)$. The vertical dotted line represents the problem size $L^{*}$ at which $t_a=t_a^{\mathrm{opt}}\left(L^{*}\right)$. To the left of this line $t_a>t_a^{\mathrm{opt}}$ and the slope of the $\text{fixed-}{t}_a$ TTS curve lower bounds the slope of the optimal TTS curve, since for very small problem sizes a large $t_a$ results in insensitivity to problem size, and the success probability is essentially constant. The opposite happens to the right of this line, where $t_a<t_a^{\mathrm{opt}}$, and where the success probability rapidly drops with $L$ at fixed $t_a$.

Figure 5

Median time to solution over instances. Plotted is $\mathrm{TTS}(L,\alpha ,0.5)$ (log scale) as a function of the Chimera subgraph size $C_L$, for a range of clause densities and for all solvers we tested. Note that only the scaling matters and not the actual TTS, since it is determined by constant factors that vary from processor to processor, compiler options, etc. All algorithms' timing reflects the result after accounting for parallelism, as described in Appendix pp1. Error bars represent $2\sigma$ confidence intervals. The DW2 annealing time is $t_a=20\mu \text{s}$.

Figure 6

Speedup ratio. Plotted is the median speedup ratio $S_X(L,\alpha ,0.5)$ (log scale) as defined in Eq. (3) for all algorithms tested. A negative slope indicates a definite slowdown for the DW2. A positive slope indicates the possibility of an advantage for the DW2 over the corresponding classical algorithm. This is observed for $\alpha >0.4$ in the comparison to SAS, SQA, SSSV, and HFS (see Fig. 10 for a more detailed analysis). Error bars represent $2\sigma$ confidence intervals.

Figure 7

Success probability correlations. The results for all instances at $\alpha =0.35$ are shown. Each data point is the success probability for the same instance, (a) for DW2 at $t_a=20\mu \text{s}$ and $t_a=40\mu \text{s}$, (b) for DW2 at $t_a=20\mu \text{s}$ and SAA at $50000$ sweeps and ${\beta }_f=5$ [in dimensionless units, such that $max\left(J_{ij}\right)=1$]. Perfect correlation means that all data points would fall on the diagonal, and a strong correlation is observed in both cases. The data is colored by the problem size $L$ and shows a clear progression from high success probabilities at small $L$ to large success probabilities at small $L$. Qualitatively similar results are seen for $t_a=20\mu \text{s}$ vs $t_a=40\mu \text{s}$ at all $\alpha$ values and for DW2 vs SAA at intermediate $\alpha$ values (see Figs. 22 and 23 in Appendix pp2).

Figure 8

Correlation between the DW2 data for two different annealing times and SAA. Plotted is the normalized Euclidean distance $D({\vec{p}}_1,{\vec{p}}_2)$ for ${\vec{p}}_1$ and ${\vec{p}}_2$ being, respectively, the ordered success probability for DW2 at $t_a=20\mu \text{s}$ and $t_a=40\mu \text{s}$ (blue circles), DW2 at $t_a=20\mu \text{s}$ and SAA (yellow squares), DW2 at $t_a=40\mu \text{s}$ and SAA (green diamonds). For comparison, the Euclidean distance between two random vectors with elements $\in [0,1]$ is $\sim 0.4$. SAA data is for $50000$ sweeps and ${\beta }_f=5$. The correlation with SAA degrades slightly for $t_a=40\mu \text{s}$. Error bars represent $2\sigma$ confidence intervals and were computed using bootstrapping (see Appendix pp2 for details). In each comparison, to construct ${\vec{p}}_1$ and ${\vec{p}}_2$ we fixed $\alpha$ and used half the instances (for bootstrapping purposes) for $L\in [2,8]$.

Figure 9

Scaling coefficients of the number of runs. Plotted here is $b\left(\alpha \right)$ [Eq. (8)] for the DW2 at all annealing times (overlapping solid lines and large symbols), for the HFS algorithm, for SAS with an optimal number of sweeps, for SAS with noise and an optimal number of sweeps, and for SAA with a large-enough number of sweeps that the asymptotic distribution has been reached at ${\beta }_f=5$. The scaling coefficients of HFS and of optimized SAS each set an upper bound for a DW2 speedup against that particular algorithm. In terms of the scaling coefficient the DW2 result is statistically indistinguishable (except at $\alpha =0.1$) from the result of running SAA at $S=50000$ and ${\beta }_f=5$. The coefficients shown here are extracted from fits with $L\ge 4$ (see Fig. 26 in Appendix pp2). Error bars represent $2\sigma$ confidence intervals.

Figure 10

Difference between the scaling coefficients. Plotted here is the difference between the scaling coefficients in Fig. 9, $b_{\mathrm{X}}\left(\alpha \right)-b_{\mathrm{DW}2}\left(\alpha \right)$, where $X$ denotes the HFS algorithm, SAS with an optimal number of sweeps, SAS with noise and an optimal number of sweeps, or SAA with $S=50000$ and ${\beta }_f=5$. When the difference is nonpositive there can be no speedup since optimizing $t_a$ can only increase $b_{\mathrm{DW}2}\left(\alpha \right)$; conversely, when the difference is positive a speedup is still possible, i.e., not accounting for the error bars, for $\alpha \le 0.05$ and $\alpha \ge 0.4$ for HFS and for $\alpha \le 0.15$ and $\alpha \ge 0.35$ for SAS without noise. These ranges shrink if the error bars are accounted for, but notably, for most $\alpha$ values SAS with noise does not disallow a limited speedup, suggesting that control noise may play an important factor in masking a DW2 speedup. Error bars represent $2\sigma$ confidence intervals.

Figure 11

Scaling of SAA at different final temperatures. SAA is run at $S=50000$ and various final inverse temperatures. The peak at $\alpha \sim 0.2$ remains a robust feature. The data marked “Noisy” is with $5\%$ random noise added to the couplings $J_{ij}$.

Figure 12

Annealing schedule of the DW2. The annealing curves $A\left(t\right)$ and $B\left(t\right)$ are calculated using rf–superconducting-quantum-interference-device models with independently calibrated qubit parameters. Units of $\hslash =1$. The operating temperature of 17 mK is also shown.

Figure 13

The DW2 Chimera graph. The qubits or spin variables occupy the vertices (circles) and the couplings $J_{ij}$ are along the edges. Of the 512 qubits, 503 were operative in our experiments (green circles) and 9 were not (red circles). We utilized subgraphs comprising $L\times L$ unit cells, denoted $C_L$, indicated by the solid black lines. There were $31,70,126,198,284,385,503$ qubits in our $C_2,\cdots ,C_8$ graphs, respectively.

Figure 14

An example of the view of the Chimera graph in the HFS algorithm. Each vertex is one half of a unit cell and is $2^4\text{-dimensional}$. The algorithm repeatedly finds the minimum energy configuration of the graph over trees (like the one highlighted) given the state of the remaining vertices (in gray).

Figure 15

Ground-state degeneracy. Number of unique solutions as a function of clause density for different Chimera subgraph sizes $C_L$. As the clause density is increased, the number of unique solutions found decreases to one (up to the global bit flip symmetry). Shown is the median degeneracy; i.e., we sort the degeneracies of the 100 instances for each value of $L$ and $\alpha$, and find the median. Our procedure counts the degenerate solutions and stops when it reaches ${10}^5$ solutions. If the median has ${10}^5$ solutions, then we assume that not all solutions were found and hence the degeneracy for that value of $L$ and $\alpha$ is not plotted. These are solutions on the used qubits (e.g., there are many instances for each $\alpha$ at $L=8$ that use $<503$ qubits); to account for the $n_{uq}$ unused qubits, we multiply the degeneracy by $2^{n_{uq}}$.

Figure 16

Scatter plot of the TTS for the HFS algorithm and the degeneracy at $L=8$ and $\alpha =0.4$ (100 instances total). Even though there is a wide range of degeneracy over several orders of magnitude, we do not observe any trend in the TTS. The degeneracy accounts for the fact that some qubits are not coupled into the problem (e.g., if $n$ qubits are not specified for that particular problem, then the degeneracy is $2^n$ times the directly counted degeneracy). The Pearson correlation coefficient is $-0.046$.

Figure 17

Comparison of the 25th (left) and 75th (right) percentiles of the TTS (log scale) for all algorithms as a function of clause density $\alpha$. The different colors represent the different Chimera sizes tested. All solvers show a peak at the same density value of $\alpha \approx 0.2$.

Figure 18

Suboptimal annealing time and optimal sweeps for $\alpha =0.2$. Plotted is the TTS (log scale) as a function of size $L$ for (a) the DW2, with all available annealing times, (b) SQA, (c) SA, and (d) SSSV, for many different sweep numbers. The lower envelope gives the scaling curves shown in Fig. 5 for $\alpha =0.2$. The TTS curves flatten at high $L$ for the following reason: Each classical annealer was run $N_X$ times ($N_{\mathrm{SA}}=N_{\mathrm{SSSV}}={10}^4, N_{\mathrm{SQA}}={10}^3$), and our $\beta$ distribution is $\beta (0.5,N_X+0.5)$ for 0 successes, which has an average value of $\sim 1/\left(2N_X\right)$. This reflects the (Bayesian) information acquired after $N_X$ runs with 0 successes (one would not expect the probability to be 0). The flattening has no impact on the scale of the optimal number of sweeps.

Figure 19

TTS (log scale) for $L=8$ as a function of number of sweeps for (a) DW2, (b) SQA, (c) SA, and (d) SSSV used to identify the optimal number of sweeps.

Figure 20

Scaling of the SAA optimal number of sweeps. The optimal number of sweeps is extracted for each $L$ from Fig. 19. The scaling is roughly exponential for the smaller $\alpha$ values and appears to be close to exponential for the larger $\alpha$ values. Lines are guides for the eye.

Figure 21

The speedup ratios (log scale) for (a) the 25th and (b) the 75th percentile of time to solution as a function of system size for various $\alpha$. The different colors denote a representative sample of clause densities.

Figure 22

Success probability correlations. The results for all instances at all $\alpha$ values are shown for the DW2 data at $t_a=20\mu \text{s}$ and $t_a=40\mu \text{s}$. This complements Fig. 7; the color scheme corresponds to different sizes $L$, as in that figure.

Figure 23

Success probability correlations. The results for all instances at all $\alpha$ values are shown for the DW2 data at $t_a=20\mu \text{s}$ and SAA with $S=50000$ and ${\beta }_f=5$. This complements Fig. 7; the color scheme corresponds to different sizes $L$, as in that figure. The correlation gradually improves from poor for the lowest clause densities to strong at $\alpha =0.35$, then deteriorates again.

Figure 24

Success probability correlations. The results for all instances at all $\alpha$ values are shown for the DW2 data at $t_a=20\mu \text{s}$ and SQAA with $S=10000$ and $\beta =5$. This complements Fig. 7; the color scheme corresponds to different sizes $L$, as in that figure. The correlation gradually improves from poor for the lowest clause densities to strong at $\alpha =0.35$, then deteriorates again.

Figure 25

Success probability correlations. The results for all instances at all $\alpha$ values are shown for the SQAA data at $S=10000$ and $\beta =5$ and SAA with $S=50000$ and ${\beta }_f=5$. This complements Fig. 7; the color scheme corresponds to different sizes $L$, as in that figure. The correlation gradually improves from poor for the lowest clause densities to strong at $\alpha =0.35$, then deteriorates again. Note that we did not attempt to optimize the correlations between the two methods.

Figure 26

Exponential fits to the DW2 number of runs. In accordance with Eq. (8), the least-squares linear fits to $ln\left[r\right(L,\alpha ,0.5\left)\right]$ are plotted in (a) for $t_a=20\mu \text{s}$. To reduce finite-size scaling effects we exclude $L=2,3$ and perform the fit for $L\ge 4$. The intercept of the linear fits is the the DW2 scaling constant $a\left(\alpha \right)$ from Eq. (8), and is shown in (b). Error bars represent $2\sigma$ confidence intervals.

Figure 27

(a) The DW2 scaling coefficients $b\left(\alpha \right)$ from Eq. (8) for SAA at various sweep numbers and ${\beta }_f=5$, along with the DW2 scaling coefficients for all $t_a$'s we tried. (b) The SQAS and SAS scaling coefficient $b\left(\alpha \right)$ at the optimal number of sweeps for each. SAS had a final temperature of ${\beta }_f=5$ and SQAS was operated at $\beta =5$.

Figure 28

Histograms of the scale factor of the instances as a function of system size for various values of $\alpha$. All problems passed to DW2 must have all couplers in the range $[-1,1]$, so all couplers in a problem are rescaled down by a factor equal to the maximum absolute value of the couplers in the problem (and hence this quantity is called the scale factor). Since internal control error (ICE) is largely instance-independent, the larger the scaling factor of an instance, the worse the relative impact of ICE will be. We see a drift to larger scaling factors for increasing size $L$ and increasing clause density $\alpha$. Larger values of $\alpha$ obviously will have larger scale factors as there are, on average, more loops per qubit (and thus, per edge) and thus a larger maximum potential coupler strength. Since the edges included in a loop are generated randomly, the more edges available at a fixed clause density, the more opportunities for a single edge to be included in many loops by chance, causing the average scale factor to drift upward as a function of problem size.