 Full text:
 PDF (eReader) / PDF (Download) (864 kB)
Ising machines have been developed as fast and highly accurate solvers for combinatorial optimization problems. They are classified based on their internal algorithms, with examples including simulatedannealingbased Ising machines (nonquantumtype Ising machines) and quantumannealingbased Ising machines (quantum annealers). Herein, we have investigated the performance of a hybrid optimization method that capitalizes on the advantages of both types, utilizing a nonquantumtype Ising machine to enhance the performance of the quantum annealer. In this method, the nonquantumannealing Ising machine initially solves an original Ising model multiple times during preprocessing. Subsequently, reducedsize subIsing models, generated by spin fixing, are solved by a quantum annealer. Performance of the method is evaluated via simulations using Simulated Annealing (SA) as a nonquantumtype Ising machine and DWave Advantage as a quantum annealer. Additionally, we investigate the parameter dependence of the hybrid optimization method. The method outperforms the preprocessing SA and the quantum annealer alone in fully connected random Ising models.
Ising machines have gained immense attention as accurate and efficient solvers for combinatorial optimization problems.^{1}^{–}^{7}^{)} These problems find the optimal combination of decision variables to either minimize or maximize the objective function, while satisfying a set of constraints. Combinatorial optimizations are applied in diverse fields such as machine learning,^{8}^{–}^{11}^{)} materials science,^{12}^{–}^{18}^{)} portfolio optimization,^{19}^{,}^{20}^{)} protein folding,^{21}^{)} traffic optimization,^{22}^{–}^{26}^{)} quantum compiler,^{27}^{)} and blackbox optimization.^{12}^{,}^{28}^{,}^{29}^{)} To tackle combinatorial optimization problems using Ising machines, the problems are mapped as a mathematically constructed model in statistical mechanics called an Ising model or its equivalent model called a quadratic unconstrained binary optimization (QUBO) model.^{30}^{)} QUBO is also referred to as the unconstrained binary quadratic programming (UBQP) problem.^{31}^{)} Some combinatorial optimization problems were also formulated.^{32}^{,}^{33}^{)} Ising machines search for the ground state of an Ising model or a QUBO model, which represents the optimal solution to a problem.
An Ising model is defined on an undirected graph
Various internal algorithms operate Ising machines.^{34}^{)} This study focuses on quantumannealingbased Ising machines (quantum annealers) and simulatedannealingbased Ising machines (nonquantumtype Ising machines). Quantum annealing (QA) is a heuristic algorithm that introduces quantum fluctuations.^{35}^{)} QA can have an advantage over classical methods through quantum tunneling due to their quantum fluctuations. Nonquantumtype Ising machines implement simulated annealing (SA), a popular algorithm to solve combinatorial optimization problems.^{36}^{–}^{38}^{)} The current advantage of nonquantumtype Ising machines over quantum annealers is that they can find lowerenergy states of largerscale Ising models. Previous studies have compiled the specifications of various Ising machines.^{39}^{,}^{40}^{)}
Several methods have been proposed to enhance the performance of a quantum annealer.^{41}^{–}^{48}^{)} These methods involve two steps. First, smaller subproblems are generated using classical techniques or a quantum annealer through variable fixing or problem decomposition. Second, these subproblems are solved using a quantum annealer. Solving subproblems through this approach has produced improved metrics compared to directly solving the original problem using a quantum annealer in some studies.^{44}^{,}^{48}^{)} However, to the best of our knowledge, no studies have clearly investigated the hybrid methods utilizing a nonquantum type Ising machine to enhance the performance of a quantum annealer.
In this study, we investigated the performance of a hybrid optimization method that combines the advantages of a nonquantumtype Ising machine and a quantum annealer. Furthermore, we analyzed how various parameters influence the performance of the hybrid optimization method. The rest of this paper is organized as follows. Section 2 introduces the hybrid optimization method. Section 3 performs the simulations using SA as a nonquantumtype Ising machine and DWave Advantage 4.1 (DWave Advantage) as a quantum annealer, with a fully connected Ising model that can be embedded in DWave Advantage. The simulations evaluate the performance of the hybrid method compared to the existing ones and investigate the influence of the parameters on the performance. Finally, Sect. 4 concludes our study and provides future directions for further research.
Here, we describe in detail our hybrid optimization method inspired by the previous study.^{46}^{)} Figure 1 depicts the method, which comprises four steps.
Figure 1. (Color online) The scheme of the hybrid optimization method. The procedure consists of four steps: sampling, conversion, execution, and evaluation. Gray, blue, and red solid circles indicate that the spin has undetermined values, −1 and +1. Black arrows represent the flow. The subIsing model size (m),
In the sampling step,
In the conversion step,
Further, the subIsing models are obtained as follows. Let n and m be the number of spins in the original Ising model and that in the subIsing model, respectively. In the original Ising model, which is represented by Eq. (1), let
Next, randomly select a tentative solution from the
In the execution step, a quantum annealer searches for the lowerenergy states of the subIsing model. The spins contained in
In the evaluation step,
The differences between the algorithm described in our study and the one in the previous study^{46}^{)} are as follows:
•  Variable fixation method: The algorithm presented in the previous study^{46}^{)} targeted a QUBO where the values of the variables take either 0 or 1. Meanwhile, in this study, we modified the approach to develop a fixed spin determination algorithm suitable for the Ising model.  
•  Loop initialization point: The loop initialization point, following the evaluation step, was set before the sampling step by tabu search in the previous study.^{46}^{)} In other words, this method utilized the solutions included in the solution pool as initial inputs for tabu search. However, we anticipate that the preprocessing will be conducted with a nonquantumtype Ising machine. Given the conditions of SA used in our study, solutions can be obtained without being influenced by these initial solutions. Consequently, if we adopt SA in the same flow as the previous algorithm, the information from the solutions in the solution pool would be lost. Therefore, for our study, we adjusted the loop initialization point to be set after the sampling step.  
•  Termination condition: In the previous study,^{46}^{)} the termination condition was set when the average Hamming distance between pairs of solutions in the solution pool was less than or equal to the subproblem size. In our study, we investigate the influences of changes to this termination condition on the solution accuracy of the hybrid optimization method. Thus, to clarify the effect of parameters, we set the termination condition to when 
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 subIsing models. Using other Ising models may lead to various graph structures or the generation of isolated spins in each subIsing model, making it difficult to consistently evaluate the hybrid method due to variations in the hardware's influence on DWave Advantage. To mitigate this, we adopted a complete graph.
In addition, when the magnetic field coefficients,
In this study, we performed simulations using SA as a nonquantumtype Ising machine and DWave Advantage as a quantum annealer. Most nonquantumtype 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
When using DWave Advantage, an original Ising model was annealed 100 times, and the best solution was selected. The other parameters of DWave Advantage were set to their default values,^{52}^{)} i.e., the annealing time was set to 20 µs, and the problem was rescaled such that the magnetic fields were in the range
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, subIsing model size (m),
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.

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.
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
Figure 4. (Color online) Performance of hybrid optimization methods with different
Figure 5 shows the energy density and number of loops with some sets of
Figure 5. (Color online) Performance of hybrid optimization methods with different sets of
Finally, we performed the hybrid method with different subIsing model sizes to evaluate the subIsing model size dependency. Figure 6 shows the energy density and number of loops for the subIsing 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 subIsing model sizes. The outer loop was set to 50. All other parameters are the same as those mentioned in Table I. Constructing subIsing models with extremely small or large numbers of spins (e.g., a subIsing model size = 16 or 144) relative to the size of the original Ising model (
Figure 6. (Color online) Performance of hybrid optimization methods with different subIsing 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.
We evaluated the performance of the hybrid optimization method, which leverages the advantages of both nonquantumtype Ising machines and quantum annealers, through simulations using SA and DWave Advantage. Furthermore, we investigated the relationship between various parameters and the performance of the hybrid optimization method.
The hybrid method achieved a higher solution accuracy than the DWave Advantage and preprocessing SA alone. Regardless of the preprocessing accuracy, the hybrid method improved the solution. Even in practical scenarios where the solution accuracy of SA is not guaranteed, the possibility of improving the solution through the hybrid method has been suggested. However, there is a tradeoff between the preprocessing accuracy and number of loops for solution improvement. We also evaluated the effect of the hybrid method parameters on the solution accuracy. Increasing the
This method first constructs a solution pool with a nonquantumtype Ising machine. Then it solves
Acknowledgments
This article is based on the results obtained from a project, JPNP16007, commissioned by the New Energy and Industrial Technology Development Organization (NEDO). The computations in this work were partially performed using the facilities of the Supercomputer Center, the Institute for Solid State Physics, The University of Tokyo. S.T. was supported in part by JSPS KAKENHI (Grant Numbers JP21K03391, JP23H05447) and JST Grant Number JPMJPF2221. The Human BiologyMicrobiomeQuantum Research Center (Bio2Q) is supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan.
References
 1 M. W. Johnson, M. H. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, Nature 473, 194 (2011). 10.1038/nature10012 Crossref, Google Scholar
 2 M. Yamaoka, C. Yoshimura, M. Hayashi, T. Okuyama, H. Aoki, and H. Mizuno, IEEE J. SolidState Circuits 51, 303 (2016). 10.1109/JSSC.2015.2498601 Crossref, Google Scholar
 3 M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, and H. G. Katzgraber, Front. Phys. 7, 48 (2019). 10.3389/fphy.2019.00048 Crossref, Google Scholar
 4 T. Inagaki, Y. Haribara, K. Igarashi, T. Sonobe, S. Tamate, T. Honjo, A. Marandi, P. L. McMahon, T. Umeki, K. Enbutsu, O. Tadanaga, H. Takenouchi, K. Aihara, K.i. Kawarabayashi, K. Inoue, S. Utsunomiya, and H. Takesue, Science 354, 603 (2016). 10.1126/science.aah4243 Crossref, Google Scholar
 5 H. Goto, K. Tatsumura, and A. R. Dixon, Sci. Adv. 5, eaav2372 (2019). 10.1126/sciadv.aav2372 Crossref, Google Scholar
 6 M. Maezawa, G. Fujii, M. Hidaka, K. Imafuku, K. Kikuchi, H. Koike, K. Makise, S. Nagasawa, H. Nakagawa, M. Ukibe, and S. Kawabata, J. Phys. Soc. Jpn. 88, 061012 (2019). 10.7566/JPSJ.88.061012 Link, Google Scholar
 7 K. Yamamoto, K. Kawamura, K. Ando, N. Mertig, T. Takemoto, M. Yamaoka, H. Teramoto, A. Sakai, S. TakamaedaYamazaki, and M. Motomura, IEEE J. SolidState Circuits (2020). 10.1109/JSSC.2020.3027702 Crossref, Google Scholar
 8 H. Neven, V. S. Denchev, M. DrewBrook, J. Zhang, W. G. Macready, and G. Rose, Quantum 4 (2009). Google Scholar
 9 M. H. Amin, Phys. Rev. A 92, 052323 (2015). 10.1103/PhysRevA.92.052323 Crossref, Google Scholar
 10 M. H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and R. Melko, Phys. Rev. X 8, 021050 (2018). 10.1103/PhysRevX.8.021050 Crossref, Google Scholar
 11 D. O’Malley, V. V. Vesselinov, B. S. Alexandrov, and L. B. Alexandrov, PLOS ONE 13, e0206653 (2018). 10.1371/journal.pone.0206653 Crossref, Google Scholar
 12 K. Kitai, J. Guo, S. Ju, S. Tanaka, K. Tsuda, J. Shiomi, and R. Tamura, Phys. Rev. Res. 2, 013319 (2020). 10.1103/PhysRevResearch.2.013319 Crossref, Google Scholar
 13 K. Utimula, T. Ichibha, G. I. Prayogo, K. Hongo, K. Nakano, and R. Maezono, Sci. Rep. 11, 7261 (2021). 10.1038/s41598021862743 Crossref, Google Scholar
 14 T. Inoue, Y. Seki, S. Tanaka, N. Togawa, K. Ishizaki, and S. Noda, Opt. Express 30, 43503 (2022). 10.1364/OE.476839 Crossref, Google Scholar
 15 K. Endo, Y. Matsuda, S. Tanaka, and M. Muramatsu, Sci. Rep. 12, 10794 (2022). 10.1038/s41598022147354 Crossref, Google Scholar
 16 H. Sampei, K. Saegusa, K. Chishima, T. Higo, S. Tanaka, Y. Yayama, M. Nakamura, K. Kimura, and Y. Sekine, JACS 3, 991 (2023). 10.1021/jacsau.3c00018 Crossref, Google Scholar
 17 A. Tučs, F. Berenger, A. Yumoto, R. Tamura, T. Uzawa, and K. Tsuda, ACS Med. Chem. Lett. 14, 577 (2023). 10.1021/acsmedchemlett.2c00487 Crossref, Google Scholar
 18 K. HatakeyamaSato, Y. Uchima, T. Kashikawa, K. Kimura, and K. Oyaizu, RSC Adv. 13, 14651 (2023). 10.1039/D3RA01982A Crossref, Google Scholar
 19 G. Rosenberg, P. Haghnegahdar, P. Goddard, P. Carr, K. Wu, and M. L. De Prado, IEEE J. Sel. Top. Signal Process. 10, 1053 (2016). 10.1109/JSTSP.2016.2574703 Crossref, Google Scholar
 20 K. Tanahashi, S. Takayanagi, T. Motohashi, and S. Tanaka, J. Phys. Soc. Jpn. 88, 061010 (2019). 10.7566/JPSJ.88.061010 Link, Google Scholar
 21 A. PerdomoOrtiz, N. Dickson, M. DrewBrook, G. Rose, and A. AspuruGuzik, Sci. Rep. 2, 571 (2012). 10.1038/srep00571 Crossref, Google Scholar
 22 F. Neukart, G. Compostella, C. Seidel, D. von Dollen, S. Yarkoni, and B. Parney, Front. ICT 4, 29 (2017). 10.3389/fict.2017.00029 Crossref, Google Scholar
 23 H. Irie, G. Wongpaisarnsin, M. Terabe, A. Miki, and S. Taguchi, International Workshop on Quantum Technology and Optimization Problems, 2019, p. 145. 10.1007/9783030140823_13 Crossref, Google Scholar
 24 S. Bao, M. Tawada, S. Tanaka, and N. Togawa, International Symposium on VLSI Design, Automation and Test (VLSIDAT), 2021. 10.1109/VLSIDAT52063.2021.9427355 Crossref, Google Scholar
 25 S. Bao, M. Tawada, S. Tanaka, and N. Togawa, IEEE International Intelligent Transportation Systems Conference (ITSC), 2021, p. 3704. 10.1109/ITSC48978.2021.9564593 Crossref, Google Scholar
 26 Y. Mukasa, T. Wakaizumi, S. Tanaka, and N. Togawa, IEICE Trans. Inf. Syst. E104D, 1592 (2021). 10.1587/transinf.2020EDP7264 Crossref, Google Scholar
 27 S. Naito, Y. Hasegawa, Y. Matsuda, and S. Tanaka, arXiv:2303.02830. Google Scholar
 28 S. Izawa, K. Kitai, S. Tanaka, R. Tamura, and K. Tsuda, Phys. Rev. Res. 4, 023062 (2022). 10.1103/PhysRevResearch.4.023062 Crossref, Google Scholar
 29 Y. Seki, R. Tamura, and S. Tanaka, arXiv:2209.01016. Google Scholar
 30 S. Tanaka, Y. Matsuda, and N. Togawa, 25th Asia and South Pacific Design Automation Conference (ASPDAC), 2020, p. 659. 10.1109/ASPDAC47756.2020.9045126 Crossref, Google Scholar
 31 G. Kochenberger, J.K. Hao, F. Glover, M. Lewis, Z. Lü, H. Wang, and Y. Wang, J. Comb. Optim. 28, 58 (2014). 10.1007/s1087801497340 Crossref, Google Scholar
 32 A. Lucas, Front. Phys. 2, 5 (2014). 10.3389/fphy.2014.00005 Crossref, Google Scholar
 33 S. Tanaka, R. Tamura, and B. K. Chakrabarti, Quantum Spin Glasses, Annealing and Computation (Cambridge University Press, Cambridge, U.K., 2017). Google Scholar
 34 N. Mohseni, P. L. McMahon, and T. Byrnes, Nat. Rev. Phys. 4, 363 (2022). 10.1038/s42254022004408 Crossref, Google Scholar
 35 T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355 (1998). 10.1103/PhysRevE.58.5355 Crossref, Google Scholar
 36 S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Science 220, 671 (1983). 10.1126/science.220.4598.671 Crossref, Google Scholar
 37 D. S. Johnson, C. R. Aragon, L. A. McGeoch, and C. Schevon, Oper. Res. 37, 865 (1989). 10.1287/opre.37.6.865 Crossref, Google Scholar
 38 D. S. Johnson, C. R. Aragon, L. A. McGeoch, and C. Schevon, Oper. Res. 39, 378 (1991). 10.1287/opre.39.3.378 Crossref, Google Scholar
 39 D. Oku, M. Tawada, S. Tanaka, and N. Togawa, IEEE Trans. Comput. 71, 223 (2022). 10.1109/TC.2020.3045112 Crossref, Google Scholar
 40 M. Kowalsky, T. Albash, I. Hen, and D. A. Lidar, Quantum Sci. Technol. 7, 025008 (2022). 10.1088/20589565/ac4d1b Crossref, Google Scholar
 41 G. Rosenberg, M. Vazifeh, B. Woods, and E. Haber, Comput. Optim. Appl. 65, 845 (2016). 10.1007/s105890169844y Crossref, Google Scholar
 42 M. Booth, S. P. Reinhardt, and A. Roy, DWave Technical Report Series 141006AA, 2017. [https://docs.ocean.dwavesys.com/projects/qbsolv/en/latest/index.html]. Google Scholar
 43 H. Karimi, G. Rosenberg, and H. G. Katzgraber, Phys. Rev. E 96, 043312 (2017). 10.1103/PhysRevE.96.043312 Crossref, Google Scholar
 44 H. Karimi and G. Rosenberg, Quantum Inf. Process. 16, 166 (2017). 10.1007/s111280171615x Crossref, Google Scholar
 45 H. Irie, H. Liang, T. Doi, S. Gongyo, and T. Hatsuda, Sci. Rep. 11, 8426 (2021). 10.1038/s4159802187676z Crossref, Google Scholar
 46 Y. Atobe, M. Tawada, and N. Togawa, IEEE Trans. Comput. 71, 2606 (2022). 10.1109/TC.2021.3138629 Crossref, Google Scholar
 47 E. Pelofske, G. Hahn, and H. Djidjev, J. Signal Process. Syst. 93, 405 (2021). 10.1007/s11265020015501 Crossref, Google Scholar
 48 E. Pelofske, G. Hahn, and H. N. Djidjev, Quantum Inf. Process. 22, 219 (2023). 10.1007/s1112802303962x Crossref, Google Scholar
 49 C. Yoshimura, M. Hayashi, T. Okuyama, and M. Yamaoka, Int. J. Netw. Comput. 7, 154 (2017). 10.15803/ijnc.7.2_154 Crossref, Google Scholar
 50 T. Okuyama, M. Hayashi, and M. Yamaoka, IEEE International Conference on Rebooting Computing, ICRC 2017  Proceedings, 2017. 10.1109/ICRC.2017.8123652 Crossref, Google Scholar
 51 Fixstars Amplify Annealing Engine: Fixstars Amplify [https://amplify.fixstars.com/en/]. Google Scholar
 52 Solver Parameters  DWave System Documentation [https://docs.dwavesys.com/docs/latest/c_solver_parameters.html]. Google Scholar
 53 N. Dattani, S. Szalay, and N. Chancellor, arXiv:1901.07636. Google Scholar
 54 The Advantage System: Performance Update [https://www.dwavesys.com/media/kjtlcemb/141054aa_advantage_system_performance_update.pdf]. Google Scholar
 55 L. Ingber, Math. Comput. Modelling 18, 29 (1993). 10.1016/08957177(93)90204C Crossref, Google Scholar
 56 E. Pelofske, G. Hahn, and H. N. Djidjev, Sci. Rep. 12, 4499 (2022).10.1038/s41598022083948 Crossref, Google Scholar