Hybrid Optimization Method Using Simulated-Annealing-Based Ising Machine and Quantum Annealer

To demonstrate the effectiveness of our hybrid optimization method, we considered a fully connected random Ising model, which was represented as a complete graph. The choice of an Ising model was motivated by the fact that once the spin is fixed on a complete graph, all spins show a uniform edge density in all sub-Ising models. Using other Ising models may lead to various graph structures or the generation of isolated spins in each sub-Ising model, making it difficult to consistently evaluate the hybrid method due to variations in the hardware's influence on D-Wave Advantage. To mitigate this, we adopted a complete graph.

In addition, when the magnetic field coefficients, \(h_{i}=0\), Ising models display a twofold degeneracy for all states, rendering the sub-Ising model generation of the hybrid method ineffective. As a countermeasure, it is proposed to set one randomly chosen spin to either −1 or +1, as described in the previous study.^43,44) However, in this study, we set the magnetic field coefficient \(h_{i}\neq 0\) and assumed no trivial degeneracy. The coefficients of the magnetic field and interaction were randomly selected according to a Gaussian distribution with a mean of zero and a standard deviation of unity. Notably, zero was excluded for both the interactions and magnetic fields.

In this study, we performed simulations using SA as a non-quantum-type Ising machine and D-Wave Advantage as a quantum annealer. Most non-quantum-type Ising machines employ SA as the basis of their internal algorithms.^{2,3,7,49–51)} The SA parameters included a randomly set initial state and an initial temperature \(T_{\text{initial}}\) of \(\lceil 2v_{\text{max}}\rceil\). Here, \(v_{\text{max}}\) was the maximum value among the \(v_{i}\) defined by \(v_{i}=|h_{i}+\sum_{j}J_{ij}|\), which is the absolute value of the sum of the magnetic fields and interactions for i-th spin in the original Ising model. The initial temperature was set sufficiently high to accommodate for the transition between arbitrary states in the beginning of SA. The temperature schedule was set to the power-law decay for every outer loop, which is given by \(T(t)=T_{\text{initial}}\times r^{t}\), where r is the cooling rate, and t is the t-th outer loop. The outer loop was set for each simulation. The cooling rate was set such that the final temperature was equal to 0.1, which is sufficiently low on the energy scale of the original Ising model. The inner loop was set at 160, which is the size of the original Ising model.

When using D-Wave Advantage, an original Ising model was annealed 100 times, and the best solution was selected. The other parameters of D-Wave Advantage were set to their default values,⁵²⁾ i.e., the annealing time was set to 20 µs, and the problem was rescaled such that the magnetic fields were in the range \([-2, 2]\), and the interactions were in the range \([-1, 1]\). To compare the performance of D-Wave Advantage, we set the number of spins in the original Ising model was set to the maximum number that can be embedded in D-Wave Advantage. D-Wave Advantage used in this study had 5,627 qubits and a Pegasus graph topology,⁵³⁾ allowing a complete graph of 177 spins to be embedded using miner embedding.⁵⁴⁾ However, several qubits could not be used due to defects. After conducting multiple tests, a stable embeddable size of 160 spins was set as the original Ising model size.

First, we performed the hybrid optimization method using several original Ising models to evaluate its performance. Table I summarizes the parameters for the outer loop, sub-Ising model size (m), \(N_{\text{I}}\), \(N_{\text{S}}\), \(N_{\text{E}}\), and \(N_{\text{L}}\). D-Wave Advantage and preprocessing SA were considered for the reference solution accuracy. The parameters for the D-Wave Advantage were utilized consistently in both as a part of the hybrid method and when used independently, as described in Sect. 3.1. Figure 2 shows the energy densities of the solutions obtained by D-Wave Advantage (DW), preprocessing SA (SA), and hybrid optimization method (HM). The SA data use the \(X_{\text{best}}\) obtained from each solution pool as preprocessing, and the energy density is the internal energy per spin (i.e., \(H/n\)). Figure 2 also shows the number of loops spent until the completion of the hybrid method. The number of loops represents the count of returning to the solution pool creation from the evaluation step. As an example, when \(N_{\text{L}} = 3\) and \(X_{\text{best}}\) was not updated at all, the number of loops was two. These data were obtained from the average and standard deviation for ten simulations. The hybrid method achieved a higher solution accuracy than either DW or SA.

»View larger version

Figure 2. (Color online) (Top) Energy densities by different solvers in multiple original Ising models. White circles are the average of ten runs. Top and bottom bars denote the highest and lowest energy densities, respectively. Black diamonds denote the outliers. (Bottom) Number of loops by the hybrid optimization method. Error bars are the standard deviations in ten runs.

»View table Table I. Parameters of SA and hybrid optimization methods. \(N_{\text{I}}\) is the number of solution samples to be obtained at the sampling step, \(N_{\text{S}}\) is the number of the original Ising models (reference) for generating sub-Ising models at the conversion step, \(N_{\text{E}}\) is the number of newly obtained solutions through the repetition of the conversion and execution, and \(N_{\text{L}}\) is the number of repeats of the same \(X_{\text{best}}\) to end the hybrid method. These parameters are described in Sect. 2 and Fig. 1.

Second, we evaluated the solution accuracy of the preprocessing dependency. This is because in practical scenarios, it can be challenging to set suitable SA parameters for a specific problem, and thus, the accuracy of SA is not guaranteed.^37,38,55) The hybrid method was performed using four different preprocessings: random spin states (Random), performing SA with the outer loop set to 10, 50, or 100. The parameters of the hybrid method are the same as those listed in Table I. Figure 3 shows the energy densities and number of loops with different preprocessings. Regardless of the solution accuracy of the preprocessing, the hybrid method improved the solutions from that of preprocessing. Even for the preprocessing with an extremely low solution accuracy (e.g., Random), the hybrid method demonstrated a certain level of solution accuracy. However, a lower preprocessing solution accuracy required more loops to improve the solution. Conversely, when the outer loop = 100, SA and HM indicated almost the same energy density. It is possibly because of the high accuracy of SA, where several of the ten simulation trials resulted in values close to the optimal solution for the problem.

»View larger version

Figure 3. (Color online) Performance of hybrid optimization methods in preprocessing with different solution accuracies. (Top) Energy densities of preprocessing and the hybrid optimization method. White circles are the average of ten runs. Top and bottom bars denote the highest and lowest energy densities, respectively. Black diamonds denote outliers. Scales of energy densities differ between the upper and lower across the omitted portions. (Bottom) Number of loops by the hybrid optimization method. Error bars are standard deviations in ten runs.

Next, we performed the hybrid method with different parameter patterns to evaluate the effect of the parameters on solution accuracy. Figure 4 shows the energy density and number of loops for \(N_{\text{L}}\) set from 1 to 5, while the outer loop was set to 50. All other parameters are the same as those in Table I. Note that the solution pool sets obtained by the preprocessing SA were the same for all conditions of \(N_{\text{L}}\). In addition, Fig. 4 shows the improvement percentage in the \(X_{\text{best}}\) of the preprocessing SA among the ten simulations of the hybrid method. The solutions were improved for \(N_{\text{L}}\geq 2\), although the solution with the highest energy corresponding to the top bar or an outlier in the box plot was not improved. As \(N_{\text{L}}\) increased, both the percentage of improved solutions and number of loops increased. There are two possible reasons for the increased number of loops. First, the hybrid method requires additional loops prior to termination. Second, the solutions improve with increasing \(N_{\text{L}}\).

»View larger version

Figure 4. (Color online) Performance of hybrid optimization methods with different \(N_{\text{L}}\). (Top) Energy densities of preprocessing SA and the hybrid method. White circles are the average of ten runs. Top and bottom bars denote the highest and lowest energy densities, respectively. Black diamonds denote outliers. Gray squares denote the percentage of improved solutions from the preprocessing SA in the ten simulations of the hybrid optimization method. (Bottom) Number of loops by the hybrid optimization method. Error bars are standard deviations in ten runs.

Figure 5 shows the energy density and number of loops with some sets of \(N_{\text{E}}\) and \(N_{\text{S}}\). Note that the sets of solution pools obtained by the preprocessing SA were the same for all conditions of \(N_{\text{S}}\) and \(N_{\text{E}}\), and the outer loop was set to 50. All other parameters are the same as those listed in Table I. The hybrid method better improved the solution more when \(N_{\text{E}}\) was 10 or 20 compared to when it was 5. This is because a larger number of \(N_{\text{E}}\) resulted in more solutions generated per loop, and there was a higher possibility of solution improvement prior to the loop termination. Additionally, the effect of \(N_{\text{S}}\) was not significant when \(N_{\text{E}}\) was the same. The number of loops was higher for \(N_{\text{E}}\) of 10 and 20, which is where an improvement is observed.

»View larger version

Figure 5. (Color online) Performance of hybrid optimization methods with different sets of \(N_{\text{E}}\) and \(N_{\text{S}}\). (Top) Energy densities of preprocessing SA and the hybrid optimization method. White circles are the average of ten runs. Top and bottom bars denote the highest and lowest energy densities, respectively. Black diamonds denote outliers. (Bottom) Number of loops by the hybrid optimization method. Error bars are standard deviations.

Finally, we performed the hybrid method with different sub-Ising model sizes to evaluate the sub-Ising model size dependency. Figure 6 shows the energy density and number of loops for the sub-Ising model sizes (m) of 16, 40, 80, 120, and 144. Notably, the sets of solution pools obtained by the preprocessing SA were the same for all sub-Ising model sizes. The outer loop was set to 50. All other parameters are the same as those mentioned in Table I. Constructing sub-Ising models with extremely small or large numbers of spins (e.g., a sub-Ising model size = 16 or 144) relative to the size of the original Ising model (\(n=160\)) do not improve the solution using the hybrid method (Fig. 6). When the sub-Ising model size was small, the fixed spins significantly influenced the sub-Ising model, and there was little freedom in the obtained solutions. By contrast, when the sub-Ising model size was large, the influence of the solution accuracy of D-Wave Advantage became more significant. If the solution accuracy of D-Wave Advantage for the original Ising model of size 160 is low (shown in Fig. 2), the sub-Ising model size increases, the solution accuracy decreases, and the solution does not improve. These results suggest that the solution accuracy of the hybrid method depends on the sub-Ising model size.

»View larger version

Figure 6. (Color online) Performance of hybrid optimization methods with different sub-Ising model size. (Top) Energy densities of preprocessing SA and the hybrid method. White circles are the average of ten runs. Top and bottom bars denote the highest and lowest energy densities, respectively. (Bottom) Number of loops by the hybrid optimization method. Error bars are standard deviations.