Parallel Quantum Annealing

Ising spin glass models have been extensively studied in the literature, for instance with respect to the stability of the system, its phase transitions, and its magnetization distribution [10]. To find the ground state of such spin glass models, as well as for other problem Hamiltonians, quantum annealing has been shown to find ground states with higher probability than classical (thermal) simulated annealing [11].

The fact that a quantum system will stay near its ground state if its Hamiltonian varies slowly enough led to adiabatic quantum algorithms [12], for which it was hypothesized early that they might be able to outperform classical algorithms on instances of NP-complete problems. In particular, quantum annealing can be shown to converge to the ground state energy at a faster rate than its classical counterpart [13]. A comprehensive review on combinatorial optimization problems (such as the ones studied in the present article), spin glasses, and quantum annealing via adiabatic evolution is available in the literature [14].

There are very few published results yet concerning parallelism and quantum annealing, and none of them discusses the type of parallelism we propose. Some papers consider pairing quantum computing systems with modern HPC infrastructure [15]. The focus there lays on the design of macroarchitecture, microarchitecture, and programming models needed to integrate near-term quantum computers with HPC systems, and the near-term feasibility of such systems. They do not consider, however, parallelism at the level of a quantum device.

In a different context [16], the authors study the problem of simulating dynamical systems on a quantum annealer. Such systems are intrinsically sequential in the time component, which makes them difficult to simulate on a quantum computer. They propose to use the parallel in time idea from classical computing to reformulate a problem so that it can be solved as the task of finding a ground state of an Ising model. Such an Ising model is then solved on a quantum annealer. However, unlike our work, there is only one Ising problem solved at a time.

Another approach considers a framework built upon Field Programmable Gate Array chips (FPGA) in connection with probabilistic bits to simulate Gibbs samplers[17]. The authors use their algorithm on blocks of conditionally independent nodes of the input graph, which are updated in parallel using a quantum annealer. The authors remark that technically, all blocks have to be completed before some synchronization step is carried out; however, if not all blocks are completed, the network is still able to find exact ground states. However, the parallel architecture proposed by the authors [17] is not quantum.

To the best of our knowledge, we are the first to propose solving multiple independent problems in parallel on a single quantum device. In a previous contribution [18], the parallel quantum annealing methodology has been used in a specific manner as a way to solve $255$ small QUBOs (with a maximum of $4$ variables each) simultaneously on the D-Wave 2000Q at Los Alamos National Laboratory with a relatively small number of anneals. In contrast to the previous work [18], the present article expands on the parallel quantum annealing idea by investigating larger problem sizes, considering the NP-hard Maximum Clique problem, and applying it to the new DW Advantage generation that uses a Pegasus qubit connection topology [19].

The D-Wave annealers offer a variety of tuning parameters that can be chosen by the user before each anneal. Those encompass the annealing time, the chain strength, anneal offsets, etc. We perform a grid search on the DW 2000Q in order to find reasonable values for both the chain strength and the annealing time, with the goal to minimize the QPU time while maximizing the probability of finding the ground state of each of the embedded problems.

In all experiments, we set the annealing time to $50$ microseconds, and compute the chain strength using the uniform torque compensation method included in the D-Wave Ocean SDK [21] with a UTC prefactor of $0.2$. Lastly, We set the programming thermalization time to $0$ microseconds, and the readout thermalization time to $0$ microseconds. We use the same parameters for the experiments with DW Advantage as we do for the ones with DW 2000Q. However, the performance of DW Advantage might be potentially improved by additional parameter optimization.

Since quantum annealers are stochastic (and heuristic) solvers, they are not guaranteed to find the ground state (the global minimum) of a problem of the type of eq. (1) in each anneal. Therefore, to have a meaningful metric of accuracy, Time-To-Solution is used to quantify the expected time it takes to reach an optimum solution with a $99$ percent confidence [22, 23, 24, 25]. It is defined, in the sequential case, as

$\displaystyle\text{TTS}_{seq}=\frac{1}{A}(T_{\text{QPU}}+U)\frac{\log(0.01)}{\log(1-p)},$

(3)

where $p$ is the probability of finding the ground state within a single anneal, $T_{\text{QPU}}$ is the total D-Wave QPU time (specifically the qpu-access-time), and $U$ is the total CPU process time used to compute the unembedded solutions in seconds. Note that in eq. (3), we do not include the embedding time because we compute full clique embeddings, and then save and re-use those embeddings.

In the context of the present work, a generalization of the standard TTS measure is required since we solve $K$ problems of size $N$ simultaneously on the same D-Wave chip, using a total QPU processing time of $T_{\text{QPU}}$ and $A$ anneals.

For each problem $i\in\{1,\ldots,K\}$ we solve simultaneously, we record the proportion of correct solutions $p_{i}$ among the $A$ anneals. We weigh those using the formula

$$p_{K}=\frac{1}{K}\sum_{i=1}^{K}p_{i},$$

where every $p_{i}$ must be non-zero. We generalize eq. (3) as

$\displaystyle\text{TTS}_{ens}(K)=\frac{1}{A}\left(\frac{T_{\text{QPU}}}{K}+U\right)\frac{\log(0.01)}{\log(1-p_{K})},$

(4)

where $U$ is the total CPU unembedding time used to unembed the solutions of each of the $K$ problems. We refer to $\text{TTS}_{ens}$ as the ensemble Time-To-Solution, since it captures the time to reach an optimal solution for each problem in a group of problems that are solved simultaneously.

This section studies the accuracy of parallel quantum annealing per se. Using the parameters of Section 2.1, we investigate the probability with which a ground state is found as a function of the problem size. Figure 2 shows results for both DW 2000Q (a) and for DW Advantage (b). We observe that, as expected, DW 2000Q finds the ground state probability (GSP) more reliably for smaller problems, though the probability varies quite considerably with the graph density. Denser problem graphs are typically easier than sparser ones. For cliques of roughly size $40$ onward, DW 2000Q is unable to find the ground state. For DW Advantage we observe a similar picture. The GSP is higher for smaller problems than for larger ones, with an even bigger spread by graph density. As before, denser problem graphs are easier to solve reliably than sparser ones. As expected, DW Advantage is able to solve larger inputs, in particular it finds the ground state with appreciable probability of up to $40\%$ even for inputs of size $40$.

Something we noted in these results is that the parallel quantum annealing method never found all ground state solutions in a single anneal. Instead, different sets of ground state solutions were found with each anneal - and over a sequence of anneals we eventually find the solutions to all of the problems. Interestingly, the rate at which the ground state solutions were found was not the same across all problems - indeed it was heavily biased. This bias is likely due to differences in the embeddings.

Similarly to Figure 2, Figure 3 shows the ensemble TTS measure for parallel quantum annealing as a function of the problem size, and for various graph densities. For DW 2000Q (a), we observe that this TTS measure increases as the problem size increases, which is as expected. The reason for the ”jump” at around problem size $40$ is unknown, but the behavior of D-Wave could be related to an internal hardware issue (during the time frame in which some of this data was taken on DW 2000Q, around problem sizes 38-46, there was an initially undetected temperature increase from the typical $0.015$ Kelvin of the hardware, likely leading to a poorer solution quality). In accordance with Figure 2, dense problem inputs seem to be easier for D-Wave, and thus incur a lower TTS value. For DW Advantage (b), we observe a similar picture. Interestingly, for almost all problem sizes that can be solved on DW 2000Q (roughly up to size $40$), the DW 2000Q is actually the better solver with respect to the TTS metric. For sizes above 40, DW Advantage is the better solver with respect to TTS.

To compare the benefits and trade-offs of solving several problems in parallel on the D-Wave quantum annealers, we first investigate the difference between the GSP found by both sequential and parallel quantum annealing. Here, sequential quantum annealing refers to solving the problems separately one after the other on the D-Wave architecture. The GSP in the parallel case is simply defined per problem as the probability of finding its individual ground state. Figure 4 shows results for DW 2000Q (a) and for DW Advantage (b). We observe that, on average across the graph densities considered, solving problems simultaneously causes D-Wave to find a ground state slightly less frequently than when solving them sequentially. For DW Advantage, this phenomenon is even more pronounced.

The reduced GSP we observe for parallel quantum annealing can be explained as follows. First, the problem complexity (i.e., the number of variables and quadratic terms) increases when solving multiple problems as opposed to one, potentially causing difficulties for the D-Wave annealer. Second, multiple embeddings on the chip mean that there is a closer physical proximity between the qubits used per embedding, potentially causing increased leakage and interaction between the involved qubits. However, unlike TTS, the GSP metric does not account for the fact that we are finding solutions to multiple problems.

Next, Figure 5 compares the sequential and parallel quantum annealing efficiencies with respect to the TTS metric. The TTS speedup is computed as the quotient of the sequential and parallel TTS values. For DW 2000Q (a), we observe that the proposed parallel quantum annealing approach yields substantial speedups (an average of $152$-fold over all problem sizes considered) compared to solving the same problems sequentially on the D-Wave chip. A stratification by graph density is not observable in this case. As expected, the speedup is more pronounced for smaller problem sizes, as in this case more problems can be solved in parallel. For DW Advantage, a considerable speedup (on average around $20$-fold) is observed over all problem sizes considered.

Note that the CPU unembedding time increases when unembedding larger datasets; meaning that for larger backends (and more variables used per anneal) the total unembedding time will increase. Therefore, the differences in TTS between the two backends seen in Figure 5 are not only caused by the differences in GSP seen in Figure 4.

We compare the performance of the proposed parallel quantum annealing approach with a classical solver for the Maximum Clique problem. We employ the fast maximum clique (FMC) solver [20]. FMC is run in exact mode throughout this experiment.

Figure 6 shows the speedup of parallel quantum annealing over FMC. For the probabilistic parallel QA, we measure ensemble TTS (see Section 2.2), whereas for the classical FMC solver, we measure CPU times (this is justified since TTS reduces to CPU time in the classical case). The speedup is computed as the quotient of the parallel TTS metric and the CPU process time used by FMC. For DW 2000Q (a), we observe that for both small problems (of any density) and for high densities up to problem size $40$, using parallel quantum annealing yields better TTS values than using FMC. In turn, FMC is superior to DW 2000Q for low density graphs. For problems exceeding size $40$, DW 2000Q is not the preferred choice (as seen in previous experiments), making FMC the faster choice for all densities. Interestingly, for DW Advantage (b), we observe that D-Wave is not superior to FMC apart from very high densities. In turn, the DW Advantage architecture does not suffer from a reduced accuracy for large problem sizes (compared to the DW 2000Q).

Finally, we consider how parallel quantum annealing can be used to increase the GSP, and therefore TTS, of a single hard QUBO problem. Such problems need to be run sequentially many times on the quantum annealer to get even a single optimal solution. Our idea is to run multiple replicas (identical copies of the QUBO) simultaneously using parallel quantum annealing.

We compare sequential quantum annealing and parallel quantum annealing on DW Advantage for a fixed set of $16$ random graphs with $N=40$ nodes and $p=0.5$. Each problem is first solved sequentially, and afterwards solved in parallel by embedding it $16$ times onto the hardware. Results are shown in Figure 7, where we see that parallel quantum annealing can reduce TTS due to improved GSP. Over all problems considered, the average TTS speedup (the quotient of sequential QA TTS and parallel QA TTS) is $7$, and the average GSP increase (the quotient of parallel QA GSP and sequential QA GSP) is $10$. Note that problem $0$ had worse TTS in the parallel quantum annealing case. This is due to an increased cost in the CPU unembedding time in the parallel case.

To compute the TTS measure of the parallel implementation, we use eq. (3), and not eq. (4). This is due to the fact that there are not $K$ distinct problem being solved in parallel; instead, we are solving the same problem $K$ times during the same annealing cycle. Therefore, $p$ is simply the proportion of anneals that found the optimal solution. Note, however, that we can have cases where the ground state solution is found multiple times in the same anneal (we do not distinguish between finding the ground state once or more than once in the embedded problems).

The Maximum Clique QUBO we use is given by [5, 26]

$\displaystyle H_{MC}=-A\sum_{v\in V}x_{v}+B\sum_{(u,v)\in\overline{E}}x_{u}x_{v},$

(5)

where the constants can be chosen as $A=1$, $B=2$ (see [26]).

Consider two QUBOs $H_{1}(x_{1},\ldots,x_{n})$ and $H_{2}(x_{n+1},\ldots,x_{m})$ with respective minimum solutions $x_{1}^{\ast},\ldots,x_{n}^{\ast}$ and $x_{n+1}^{\ast},\ldots,x_{m}^{\ast}$. Since both QUBOs do not share any variable, the ground state of $H_{1}+H_{2}$, which is now a function in $x_{1},\ldots,x_{n},x_{n+1},\ldots,x_{m}$, is precisely the sum of the ground states of $H_{1}$ and $H_{2}$, and the minimum bitstring yielding the ground state will be $x_{1}^{\ast},\ldots,x_{n}^{\ast},x_{n+1}^{\ast},\ldots,x_{m}^{\ast}$. This idea lays at the heart of parallel quantum annealing. Given we have an embedding that allows us to place $H_{1}$ and $H_{2}$ simultaneously on the quantum chip, we can therefore embed and solve $H_{1}$ and $H_{2}$ in one D-Wave call. Naturally, $H_{1}$ and $H_{2}$ do not need to be QUBOs of the same type of optimization problem, and neither do they need to be of equal size. Generalization to a larger number of QUBOs is straightforward.

The D-Wave hardware has a hardware precision limit [27, 28, 29] for the provided problem coefficients. Therefore, when constructing these independent QUBOs, another point to consider is how the minimum and maximum range of these independent QUBOs compare to each other. If one of these QUBOs’ coefficients dominate all of the other QUBOs, then most of the QUBO coefficients will effectively be lost in noise. Therefore, it might be necessary to normalize the minimum and maximum QUBO coefficients across all of the QUBOs that are being solved at the same time. In the experiments shown in Section 2, all of the problems are Maximum Clique QUBOs, which always have the same minimum and maximum coefficient range - therefore no coefficient normalization is used.

We aim to create multiple disjoint minor-embeddings of a complete graph of size $N$ onto the DW 2000Q and the DW Advantage hardware. This allows us to embed arbitrary QUBOs of size $N$ onto the hardware connectivity graph, meaning that we can solve multiple problems of arbitrary structure (up to size $N$) in parallel in a single backend call.

We compute separate embeddings for clique sizes $N$ in the range of $N\in\{8,\ldots,46\}$ for DW 2000Q, and $N\in\{10,\ldots,60\}$ for DW Advantage. As shown in Figure 1, ideally the embeddings should be chosen such that both the number of cliques (complete graphs) embedded on the quantum hardware, as well as the physical distance between the embeddings is maximized. The latter is conjectured to decrease spurious interactions between qubits of different QUBOs on the chip, thus helping to increase the solution quality.

The disjoint embeddings for both DW 2000Q and DW Advantage are computed using the method minorminer [30, 31]. In particular, we supply the graph to be embedded to minorminer as simply the union of the disjoint cliques we wish to embed. This is done for all clique sizes $N$ for which an embedding can be found by minorminer. Figure 8 shows both the number of embeddings we are able to place on the D-Wave quantum chip, as well as the chain lengths occurring in those embeddings for both DW 2000Q and DW Advantage.

Supplying minorminer with a graph consisting of the disjoint cliques we wish to embed is sub-optimal. In particular, when computing the embeddings, we would ideally be optimizing for the separation between the embeddings on the hardware. This remains for future work, thus close physical proximity of our embeddings may cause leakage between qubits.

As described in Section 1, after annealing, we need to unembed all chains. For this, we tried both the majority vote method and the random weighted method. All results shown in this paper use the random weighted unembedding method.

For each of the number of disjoint embeddings we fixed in Section 4.3, we generate an Erdős–Rényi random graph $G(n,p)$ of size $n=N$, varying edge probability $p$ from $0.10$ to $0.95$ in increments of $0.05$ such that each random graph is connected, does not have any degree $0$ nodes, and is not a full clique itself. Then, we compute the exact Maximum Clique solution using FMC [20] and Networkx [32, 33, 34, 35] (the Networkx solver is used to find any degenerate Maximum Clique solutions). This step is necessary in order to find the CPU time used by FMC, as well as determine the exact solution(s) for the purpose of computing ground state probability (GSP), which is used to compute the quantum annealer Time-to-Solution (TTS) metric.

Next, we arbitrarily assign each of the random graphs to one of the disjoint embeddings. Then, the Maximum Clique QUBO is computed for each of these random graphs, and each QUBO is then embedded onto its assigned embedding. For example, Figure 1 (a) shows that on DW 2000Q for $N=20$, we generate $12$ random graphs of size $20$ and then embed those 12 Maximum Clique QUBOs onto each of their respective (independent) embeddings. Then, this is repeated for each density.

Once the problems are embedded, we call the D-Wave backend using the parameter setting of Section 2.1. In order to get reliable results for the TTS computation, we call the backend $100$ times for each problem, and each backend call requests $1,000$ anneals. This results in $100,000$ samples in total per random graph. From these results, we first unembed the samples using the random weighted method, and then we compute how many of the $100,000$ samples correctly found the Maximum Clique solution. Using this information, we then compute the ensemble TTS metric using eq. (4).

Next, we apply the sequential quantum annealing method; meaning that we solve each of those QUBOs using separate backend calls. Importantly, each QUBO still uses the same physical hardware embedding as it did in the parallel case, so that differences in the embeddings will not change solution quality. In the example of Figure 1 (a), we would perform the $100$ backend calls for each one of the $12$ distinct problems, without the other $11$ embedded into the chip. This procedure is then repeated for all graph densities as well. Once again, we unembed using the random weighted method, and we also compute GSP for each of the problems. This procedure allows us to directly compare the solution quality for each of these problems when solving them separately and in parallel. Using this data, we can compute the TTS using eq. 3 for each of the random graphs. Thus, the total TTS for a group of graphs that were solved sequentially is the sum of each of their individual TTS values (in contrast to using parallel quantum annealing where we compute the ensemble TTS using eq. 4).

It is important to note that, especially for larger clique sizes (e.g., cliques of size $40$ and larger for DW 2000Q, and clique sizes $50$ onward for DW Advantage), not all problems can be solved to optimality by the quantum annealer. In this case, we encounter missing information in the computation of the TTS metric. We mitigate this problem by attempting to solve another randomly generated problem instead until the TTS metric can be computed. However, in settings where particular graphs need to be solved (as opposed to the setting with random graphs we consider), it is likely that D-Wave fails for large problems solved using a fixed embedding [1]. These failures can be attributed to longer chain lengths leading to broken chains, and therefore worse solution quality. Our results in section 2.6 show a possible solution to this problem; that is using many different embeddings to solve the same problem in parallel might allow D-Wave to more consistently find optimal solutions for large problem sizes.