Increasing the Hardness of Posiform Planting Using Random QUBOs for Programmable Quantum Annealer Benchmarking

Roughly, our algorithm generates a set of smaller non-overlapping random QUBOs and one larger posiform planted QUBO that covers the hardware graph. The smaller random QUBOs are then fused with the large posiform planted QUBO covering the hardware graph. The aim of combining the larger (posiform planted) QUBO with smaller QUBOs is to alter the coefficients in such a way that the QUBO becomes harder to solve. At the same time, it is possible to add on the smaller QUBOs to the posiform planted QUBO without changing its optimal solution.

The disjoint graphs corresponding to these QUBOs are generated using the Networkx [76] implementation of the Kernighan–Lin bisection algorithm [77]. Specifically, the hardware graph is recursively bisected into equally sized partitions until a threshold is reached. The total number of partitioned subgraphs is increased starting at a maximum number of variables (for all subgraphs) of $50$. This recursive bisection technique ensures that the final partitioned subgraphs will be identically sized, unless there are any odd divisions in which case the number of variables in the partitioned subgraphs will be different from each other by at most $1$ variable (this is only true because the number of partitions is always a power of 2).

Note that the construction of these disjoint hardware subgraphs is similar to the tiled embeddings used in parallel quantum annealing [29, 31, 32, 78, 79]. Figure 14 in the appendix shows an example of the partitioned Advantage_system4.1 hardware graph and Figure 15 shows the same for the Advantage2_prototype1.1 device.

The bisection partitioning algorithm is stochastic, and so different runs can produce different partitions. When generating the random QUBO instances, a new different disjoint bisection partition is constructed for every instance.

The steps used to construct the hardware-native QUBO problems are defined as follows:

1. 1.

   Generate the bisected subgraphs that cover the entire hardware graph applying the Kernighan–Lin bisection algorithm recursively.

2. 2.

   Choose linear and quadratic random coefficients for all of the edges and nodes within the disjoint induced subgraphs of the hardware graph computed in the previous step. Here, we generate two distinct types of random QUBOs using two different discrete coefficient sets drawn from either $\{-1,-0.9,\dots,0.9,1\}$ (with $20$ different values) in increments of 0.1, or $\{-1,1\}$ (with $2$ values). We will denote these two different types of QUBOs using the shorthand of lin₂₀ and lin₂, respectively.

3. 3.

   Use an exact classical solver to compute an optimal variable assignment for each of the random subproblems. In this case, we use CPLEX [80], which can deterministically find an optimal variable assignment (given sufficient compute time). To this end, the optimization problem is formulated as a binary Mixed Integer Quadratic Program (MIQP).

4. 4.

   Concatenate together each of the single optimal variable assignments for each random problem found in the previous step (do this based on some consistent bit-ordering, in this case we used consistent node indexing of the hardware graph). This concatenated bitstring thus gives a single variable state for every qubit (spin) in the hardware graph.

5. 5.

   Use the posiform planting algorithm introduced in [59] to generate a QUBO whose planted ground-state matches exactly the concatenated bitstring computed in the previous step. There are various valid choices for generating this posiform QUBO, although some may affect the coefficient range of the resulting QUBO problem, potentially leading to suboptimal performance when encoded onto current D-Wave quantum annealers. In our case, we use the 2-SAT solver Minisat [81] (as in [59]). All posiform coefficients corresponding to a 2-SAT clause are chosen as $1$. To conserve compute time, we only attempt to solve the 2-SAT instance being generated in batches. This may cause the generated 2-SAT instance to be (slightly) larger than is needed to ensure the uniqueness of the planted solution.

6. 6.

   ”Glue” together all of the random QUBO problems and the posiform planted QUBO, which is done by simply summing all of the terms together. The obtained final QUBO inherits from posiform planting that it has a unique and known optimal solution. The posiform planted QUBO spans the entirety of the chip and connects together the small random QUBOs, making this final QUBO connected (and in particular, it covers a majority of the edges in the hardware graph and all of the nodes in the hardware graph). A parameter that we add into this QUBO combination step is a constant scale factor which we will refer to as the posiform scaling coefficient. This coefficient multiplies all of the terms in the posiform QUBO prior to their addition to the random QUBO problems. We evaluated this coefficient experimentally at $0.1$ and $0.01$. The reason for this is empirical. We observed in experiments that scaling the posiform QUBO coefficients by a factor prior to combining it with the smaller QUBOs resulted in harder problems for simulated annealing.

The process outlined in Section 2.1, whereby smaller QUBOs (solved to optimality with CPLEX) are combined with the larger posiform planted QUBO (having a unique solution), preserves the uniqueness of the solution. This section gives a formal proof of this statement.

Although some of the random QUBOs being generated may have multiple solutions (i.e., multiple minima), the combined final QUBO has a unique one. To see that, consider the random QUBOs $R_{1},\dots,R_{k}$ and the posiform QUBO $P$. Denote the variables of $P$ by $X=\{x_{1},\dots,x_{n}\}$, and let $\mathrm{sub}_{i}(X)$ represent the subset of $X$ used as variables in each $R_{i}$, for $i=1,\dots,k$. Let $X^{*}=\{x_{1}^{*},\dots,x_{n}^{*}\}$ be the unique minimum of $P$. We claim that $X^{*}$ is also a unique minimum of the final QUBO $Q=\sum_{i}R_{i}+\alpha P$, where $\alpha$ is the posiform coefficients scale factor.

To prove that claim, let $\hat{X}=\{\hat{x}_{1},\dots,\hat{x}_{n}\}$ be a binary assignment of the variables such that $\hat{X}\neq X^{*}$. Since, by construction, $\mathrm{sub}_{i}(X^{*})$ is a minimum of $R_{i}$ for each $i=1,\dots,k$, it follows that

Furthermore, $X^{*}$ is a unique minimum of $P$ and $\hat{X}\neq X^{*}$. Then

Combining eq. (3) and eq. (4), we get

Hence, $Q(X^{*})<Q(\hat{X})$, and since this holds for any $\hat{X}\neq X^{*}$, $X^{*}$ is a unique minimum of $Q$.

The compute time required to sample an optimal solution with high confidence, in particular for heuristic probabilistic solvers, is given by the Time-to-solution (TTS) measure [82], which is defined as

where $T_{\text{anneal}}$ is the average QPU time used per anneal-readout cycle and $p$ is the proportion of samples that correspond to an optimal solution. $T_{\text{anneal}}$ is measured by dividing the total QPU time (total QPU time is defined as the QPU access time from the D-Wave system, which includes programming time, readout time, and anneal time) by the number of anneals measured in the given dataset. This TTS measure quantifies the expected compute time required to sample an optimal solution at least once, with 99% probability. Note that the only compute time used to measure the TTS is the QPU access time; this does not include local embedding and data processing CPU time.

The generated QUBO problems are also solved using simulated annealing [75] to compare the performance of quantum annealers against a classical general-purpose solver. Simulated annealing performance is evaluated as a function of the number of Metropolis-Hastings update sweeps, ranging from 1 to 10,000. The comparison uses the Python 3 package dwave-neal²²2https://github.com/dwavesystems/dwave-neal, a C++ implementation with a Python wrapper. The default geometric annealing schedule is used, generating $10,000$ samples for each parameter setting and problem instance. In this implementation, a sweep of variable updates is performed in a fixed order for each step of $\beta$, where $\beta$ corresponds to the number of Metropolis-Hastings updates in the simulated annealing schedule. The simulated annealing CPU time is measured as process time.

The Gurobi optimization software [83] is used to solve the same QUBO instances that are solved using simulated annealing and quantum annealing (up to problem sizes where it is still feasible). Gurobi provides an optimality guarantee of found solutions, and can deterministically find an optimal variable assignment of the optimization problem. We use Gurobi to solve the QUBO problems as binary Quadratic Programs. Gurobi version 11.0.3 is used for all simulations, and the compute platform is a Red Hat Linux node with Intel(R) Xeon(R) CPU E5-2695 v4 2.10GHz. Gurobi settings are all default except the following; a single thread is used, the compute time limit is $4,000,000$ seconds, and the MIP gap is $1\mathrm{e}{-8}$. The compute time reported as the classical CPU time time.

The annealing times used to sample the problem instances is varied over $\{1,2,3,4,\ldots 19,20,100,1000,2000\}$ microseconds. In the case of Advantage_system4.1 and Advantage2_prototype2.3, annealing time of $500$ nanoseconds is also used. The measured sample distributions are computed by executing $100$ anneal-readout cycles for each job, which is then repeated multiple times for each parameter configuration so as to build a large sample distribution for each QUBO problem instance. Some characteristics of the three D-Wave systems are given in Table 1. The D-Wave QPUs used in this study typically operate at or near $16$ milliKelvin.

The Advantage2_prototype1.1 with posiform scaling coefficient $0.01$ hardware compatible QUBO models (for both lin₂₀ and lin₂ cases) are sampled on the D-Wave hardware using $10000$ samples for each of the random $100$ QUBO models and annealing times. All other QUBO models are sampled with $1000$ samples per parameter combination. The Advantage2_prototype1.1 QUBOs with posiform scaling $0.01$ are solved with a larger number of samples because the sampling success rate was found to be especially low ($0$ for many problem instances), and thus we increased the sample count with the goal of being able to quantify the TTS.

The largest random QUBOs that are exactly solved using CPLEX is limited by the compute time required to solve the models. We limit the compute time to be a maximum of 2 days on a single compute node, using a single thread. For any larger Zephyr graphs, CPLEX reached this compute time limit, as did larger Pegasus graph instances. Therefore, larger random QUBOs could be solved, but would require either using a faster optimization software, or using much more compute time.

In this section we investigate the success rate with which the optimal (that is, the planted) solution is being sampled on different D-Wave devices. The devices being considered are Advantage_system4.1, Advantage2_prototype1.1, and Advantage_prototype2.3. To this end, we generate $100$ random posiform planted QUBOs which always cover the entire hardware graph of the respective architecture being used, and glue on a varying number of smaller random-coefficient QUBOs having appropriately chosen sizes.

We start with Advantage_system4.1. When gluing on smaller QUBOs to the posiform planted QUBO, we employ a posiform QUBO scaling coefficient of $0.1$. Figure 1 reports the success rate for Advantage_system4.1 in finding the unique planted optimal solution. All proportions being reported are with respect to the generated ensemble of $100$ posiform planted QUBOs. Figure 1 presents success rates as a function of annealing time. We observe that the optimal solution sampling success rate is generally low, with a maximum of approximately $0.25$. Although the shapes of the plots in the two columns for lin₂ and lin₂₀ look similar, the scales on the left are much lower, indicating smaller success probabilities. We observe that smaller annealing times lead to a higher success rate on the Advantage_system4.1, and that especially for coefficients chosen from lin₂₀ (right column), success rates flatten off closer to $0$ with longer annealing times.

We repeat the same experiment on the Advantage2_prototype1.1 device using a posiform QUBO scaling coefficient of $0.1$ when combining the posiform planted QUBO with smaller QUBOs to disguise the solution. Similarly to the previous experiment, Figure 2 again shows the the ground-state sampling success rates. We observe in Figure 2 that in contrast to the previous experiment on Advantage_system4.1, the success probabilities for Advantage2_prototype1.1 are much higher for low annealing times, with many generated QUBOs being solvable with probability $1.0$. Moreover, there is a clear decrease to zero in probability as the annealing times increase. This pattern in consistent across the two coefficient choices (left and right columns), as well as the number and sizes of the smaller QUBOs being glued onto the posiform planted QUBO.

This experiment is repeated for Advantage_prototype2.3 using a posiform QUBO scaling coefficient is $0.1$ in Figure 3. The results in Figure 3 for the Advantage_prototype2.3 device confirm the ones seen for Advantage2_prototype1.1 in that the achieved success probabilities are close to $1$ for low annealing times and decrease to closer to $0$ as the annealing times increase. As before, all subplots show a qualitative similar behavior for the different coefficient choices (columns) or the number and sizes of the smaller QUBOs being glued onto the posiform planted QUBO (rows).

Figure 4 shows the quantum annealing ground-state sampling success rate for QUBOs tailored to the Advantage2_prototype1.1 hardware graph, but now using a posiform scaling coefficient of $0.01$. This shows a clear decrease of ground-state sampling success as compared to Figure 2 with a larger posiform scaling coefficient. We observe that by making the posiform scaling coefficient small and making the random QUBO coefficient discrete in $\{+1,-1\}$, we can make these solution planted QUBOs much harder for the D-Wave hardware to solve optimally.

Figure 5 shows quantum annealing ground-state sampling success rates for the QUBOs defined on the Advantage2_prototype2.3 hardware graph, again with a posiform scaling coefficient of $0.01$ (as opposed to Figure 3). Here we see a significant decrease in optimal solution sampling rate, where for the lin₂ QUBO instances (left column) the unique optimal is rarely sampled. Thus we again observe that the posiform scaling coefficient attenuates the computational hardness of these planted QUBOs. The overall consistent finding is that smaller posiform QUBO coefficients with respect to the fused random QUBOs make the problems harder for the D-Wave quantum annealers to solve.

Figures 1, 2, 3 all show a consistent trend of higher ground-state sampling rates at smaller annealing times, and lower ground-state sampling rates at longer annealing times. This is a very interesting property of quantum annealing and these particular optimization problems – for many other optimization problems solved with quantum annealing in other studies, this trend is exactly reversed; usually better optimal solution sampling occurs at longer annealing times [30, 84, 15, 85, 17, 16, 56].

Figures 1, 2, 3, 4, 5 all show a very consistent trend that the sampling success rate, although it does depend on annealing time, does not change noticeably as a function of the varying random QUBO sizes that were fused to the posiform planted QUBO. This can be seen by looking at the differences between the rows of all of these figures – each row denotes a different random QUBO size threshold. This is a counter-intuitive finding since larger random QUBOs generally require more compute time to solve to optimality, and thereby we would expect to see lower sampling success rates.

For all problem instances generated on Advantage_system4.1 using posiform scale factor of $0.01$, the D-Wave hardware never sampled any of the optimal solutions out of the $1000$ samples obtained for each individual instance and QUBO setting configuration. For this reason, there are no ground-state sampling probability figures for this particular QUBO configuration, as there is for posiform coefficient scaling of $0.1$ on the Advantage_system4.1 instances (Figure 1). However, there is a natural question of how close the D-Wave the hardware was to sampling the planted solution. Figure 6 answers this question by plotting the energy difference between the planted solution energy and the lowest energy sample found on D-Wave hardware. The more negative this quantity is, the further away the quantum annealing hardware sample was from the planted solution, and the closer to $0$ this quantity is, the closer the hardware sampling was to the optimal solution.

In order to have a base for comparison we repeat the experiments of the previous section and solve the same QUBOs with the classical simulated annealing algorithm [75]. In particular, we solve the same set of QUBOs as in Section 3.2.

Figure 7 shows simulation results when solving the same set of QUBOs generated for Figure 5, that is QUBOs tailored to the Advantage2_prototype2.3 hardware graph. Two observations are noteworthy. First, the displayed plots are qualitatively similar, in that they all show that the success probability increases with the number of Metropolis-Hastings updates, which is expected. In particular, the plots are similar for the QUBOs generated for QUBOs with coefficients from lin₂ (left column) and lin₂₀ (right column). Second, we observe that for a low number of Metropolis-Hastings updates, the success probability is (close to) zero, showing that the generated QUBOs are not trivial, and that moreover, more Metropolis-Hastings updates are needed for a non-zero success probability when the coefficients come lin₂₀ (right column).

Figure 8 shows qualitatively similar results to Figure 7, where a small number of spin updates yields low ground-state success rate and a larger number of spin updates yeilds ground-state success rates closer to $1$. However, Figure 8 with posiform scaling of $0.1$ has a much higher ground-state success probability compared to Figure 7 with posiform scaling of $0.01$. This shows that the posiform scale factor can attenuate the hardness of these QUBO problems for simulated annealing. This same quality of the larger posiform scaling coefficients being easier was seen in the quantum annealing sampling results.

Appendix C shows three additional simulated annealing sampling figures (Figures 18, 19, 20). All of the simulated annealing sampling plots show that the optimal solution is able to be reliably sampled for all QUBO configurations if a sufficiently large number of variable spin updates is utilized.

An interesting property of the ground-state sampling curves from simulated annealing for the posiform scaling $0.01$ and lin₂ coefficient distributions (see the left columns of Figures 5 and 7) is a slight down-turn of success probability just before $10$ Metropolis-Hastings spin updates, before the success probability continues to increase. This non-monotonic behavior seems to be a characteristic of these specific problem instances as it is not seen in the other posiform solution planted QUBO configurations, and the cause is not known.

In this section we record the TTS metric for solving the previously considered ensembles of QUBOs. The TTS is computed as in Section 2.3. Here TTS is computed for each individual hardware tailored QUBO instance using the entirety of the parameter sweep and all samples obtained for each QUBO, the results for which are described in Section 3.2. In particular, this means that the sample success rate is measured from the entirety of the samples across all annealing times used, and therefore the QPU time per anneal is the average QPU time across different annealing times. This approach to measure TTS was taken because the sampling success rate can vary quite significantly as a function of the annealing time. In particular, for some annealing times the measured optimal solution sampling success rate is $0$.

Figure 9 shows TTS measurements for the full ensemble of posiform planted QUBOs with posiform scale of $0.1$, sampled on the three D-Wave QPUs. Figure 10 shows the same for the posiform scalng of $0.01$ (not including Advantage_system4.1 because the optimal solution sampling rate was always zero). Each subplot shows a scatterplot of the TTS measurements as a function of the maximum random QUBO sizes that were glued to the whole-chip posiform QUBO. As the TTS cannot be measured if the optimal solution is never found, there are missing data points in these TTS measures.

Figures 9 and 10 show that the measured TTS largely stays in the same range irrespective of the size (and number) of the QUBOs glued onto the whole-chip posiform QUBO.

Table 2 gives exact TTS quantities, and range of TTS quantities, for each QUBO parameter type.

Figure 11 shows the Gurobi classical compute time required to deterministically find the unique optimal solution of the posiform planted instances, defined on Advantage2_prototype1.1. Figures 12 and 13 shows the same Gurobi runtime scaling for the QUBOs defined on the other two D-Wave QPUs. However, in these cases, the cumulative compute time was too large to compute in a reasonable amount of time for the larger random QUBO subproblem sizes, and therefore these two figures are limited to at most the $200$ variable size cutoff.

Fourth observations are noteworthy. First, for Gurobi there is a clear dependence on the size of the random QUBO instances for the required runtime to find the optimal solution. This is consistent with previous studies that have used integer programming tools such as Gurobi to solve random spin glass instances defined on D-Wave hardware graphs [15]. In general, these random spin glasses get significantly harder to solve optimally as the problem size increases. This is significant because it shows that despite the gluing of the relatively easy-to-solve posiform planted QUBOs, the random QUBOs do make the overall problem very computationally challenging for state of the art optimization software such as Gurobi. Notably, this characteristic of random QUBO size dependence is not shared by either quantum annealing or simulated annealing ground-state sampling results.

Second, we see that the more discrete coefficient problem instances are harder to solve using Gurobi, and the finer coefficient precision problem instances are easier to solve using Gurobi. This property is also seen in the simulated annealing and quantum annealing results.

Third, the QUBOs with posiform scaling coefficient $0.01$ and $\pm 1$ coefficients are clearly the hardest instances to solve using Gurobi. This scaling property is consistent with both the quantum annealing results and the simulated annealing results.

Fourth, the exact runtimes used to solve each instance vary quite significantly for each random QUBO size. For some instances, the range of solution times spans several orders of magnitude.