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

Ising machines have the potential to realize fast and highly accurate solvers for combinatorial optimization problems. They are classified based on their internal algorithms. Examples include simulated-annealing-based Ising machines (non-quantum-type Ising machines) and quantum-annealing-based Ising machines (quantum annealers). Herein we propose a hybrid optimization method, which utilizes the advantages of both types. In this hybrid optimization method, the preprocessing step is performed by solving the non-quantum-annealing Ising machine multiple times. Then sub-Ising models with a reduced size by spin fixing are solved using a quantum annealer. The performance of the hybrid optimization method is evaluated via simulations using Simulated Annealing (SA) as a non-quantum-type Ising machine and D-Wave Advantage as a quantum annealer. Additionally, we investigate the parameter dependence of the proposed hybrid optimization method. The hybrid optimization method outperforms the preprocessing SA and the quantum annealing machine alone in fully connected random Ising models.

An Ising model is defined on an undirected graph G = (V, E), where V and E are sets of vertices and edges, respectively.It consists of spins, interactions, and magnetic fields.The Hamiltonian (or energy function) H of the Ising model is defined as where σ i is the spin on the vertex i ∈ V and takes a value of either −1 or +1.The coefficients of the RHS in Eq. ( 1), h i and J i j , are the magnetic field on the vertex i ∈ V and the interaction on the edge (i, j) ∈ E, respectively.Various internal algorithms operate Ising machines. 34)This * skikuchi@keio.jpstudy focuses on quantum-annealing-based Ising machines (quantum annealers) and simulated-annealing-based Ising machines (non-quantum-type Ising machines).Quantum annealing (QA) is a heuristic algorithm that introduces quantum fluctuations. 35)Quantum annealers show quantum tunneling due to their quantum fluctuations.Consequently, quantum annealers should be advantageous compared to classical methods.][38] The current advantage of non-quantum-type Ising machines over quantum annealers is that they can find lower-energy states of larger-scale Ising models.Previous studies have compiled the specifications of various Ising machines. 39,40) ][43][44][45][46][47][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 with an Ising machine.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) owever, no studies have investigated hybrid methods that utilize different types of Ising machines.
In this study, we propose a hybrid optimization method that combines the advantages of a non-quantum-type Ising machine and a quantum annealer, examine the parameter dependence, and evaluate its performance via simulations.The rest of this paper is organized as follows.Section 2 introduces the hybrid method.Section 3 performs the simulations using SA as a non-quantum-type Ising machine and D-Wave Advantage 4.1 (D-Wave Advantage) as a quantum annealing machine, with a fully connected Ising model that can be embedded in D-Wave Advantage.The simulations investigate the influence of the parameters on the performance and evaluate the per-formance of the hybrid method compared to existing ones.Finally, Sect. 4 concludes our study and provides future research.

Hybrid optimization method
Here we detail our hybrid optimization method inspired by the previous study. 46)Figure 1 depicts the proposed method.The proposed method is composed of four steps.
In the sampling step, N I solutions are sampled using a nonquantum-type Ising machine and registered in the solution pool.Because a non-quantum-type Ising machine shows a stochastic behavior, solutions with different spin states can be sampled.The solution pool may contain duplicate solutions.This step ends by temporarily selecting the solution with the lowest energy in the solution pool as the best solution X best .
In the conversion step, N S solutions (X 1 , X 2 , •••, X N S ) are randomly selected from the solution pool, and using an approach called "sample persistence". 43,44) he spins are fixed to obtain a sub-Ising model.Duplicate solutions may be included when performing N S random selections.Sample persistence assumes that spins with the same value across independently obtained samples are candidates to be fixed at a persistent value.The remaining spins are a group of difficult spins to solve.
Then sub-Ising models are obtained as follows.Let n and m be the number of spins in the original Ising model and that in the sub-Ising model, respectively.In the original Ising model, which is represented by Eq. ( 1), let S = { σ 1 , σ 2 , ..., σ n } be a set of spins.First, d i is calculated for each spin of the N S solutions to determine the spins not fixed at −1 or +1 into the sub-Ising model.The value of d i is given as where σ i,k is the i-th spin in the k-th solution.If N S is even (odd), the minimum value of d i is 0 (1).The minimum value of d i is the most unstable spin, whereas the maximum value of d i is N S , which is a persistent spin.To obtain the more unstable spins, collect m spins of the sub-Ising model in order of the smallest d i .Let S ′ be a set of spins in the sub-Ising model.
Next, randomly select a tentative solution from the N S solutions.Let σi = ( σ1 , σ2 , ••• , σn ) be a spin state of the tentative solution, where σi is −1 or 1.Then, the Hamiltonian of the sub-Ising model is generated from S ′ and σi as where In the execution step, a quantum annealer searches for lower-energy states of the sub-Ising model.The spins contained in S ′ of the tentative solution are updated with the spin state of the sub-Ising model.This becomes the new solution.Then, the conversion step and the execution step are repeated N E times to obtain N E new solutions are obtained.These new solutions are added to the solution pool, constructing an expanded solution pool.
In the evaluation step, N I low energy solutions are collected from the expanded solution pool to create a new solution pool.The lowest energy solution in the new solution pool is updated as the new X best .Then the flow from the solution pool creation to the evaluation step is iterated using the new solution pool.This is repeated until X best remains the same for N L consecutive times.Finally, the hybrid optimization method outputs the lowest energy in the solution pool X best at the end of the iteration as the solution.

Setup of experiment
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 all spins show a uniform edge density in all sub-Ising models once the spin is fixed on a complete graph.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 a sub-Ising model generation of the proposed hybrid method ineffective.As a countermeasure, we set one randomly chosen spin to either −1 or +1 as described in the previous study. 43,44) owever, in this study, we set the magnetic field coefficient h i 0 and assumed no trivial degeneracy.The coefficients of the magnetic field and the interactions were randomly selected according to a Gaussian distribution with a mean of zero and a standard deviation of unity.It should be noted that zero was excluded for both the interactions and the 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.51][52] The SA preprocessing parameters were a randomly set initial state and an initial temperature T initial of ⌈2v max ⌉.Here, v max waas the maximum value among the v i defined by v i = h i + j J i j , which is the absolute value of the sum of the magnetic fields and the 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 powerlaw decay for every outer loop, which is given by T (t) = T initial × r t , where r is the cooling rate and t is the t-th outer loop.The outer loop was set for each original Ising model.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, one original Ising model was annealed 100 times, and the best solution was selected.

Spins
Repeat the loop until the solution with the lowest energy in the solution pool ( $%&' ) remains the same for   consecutive times.The other parameters of D-Wave Advantage were set to their default values. 53)To compare the performance of D-Wave Advantage, 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, 54) allowing a complete graph of 177 spins to be embedded using miner embedding. 55)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.

Numerical experiment
First, we performed the proposed 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 I , N S , N E , and N L .D-Wave Advantage and preprocessing SA were considered for the reference solution accuracy.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 X best obtained from each solution pool as the 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 is complete.The number of loops represents the count of returning to the solution pool creation from the execution step.As an example, when N L = 3 and X best was not updated at all, the number of loops is 2.These data were obtained from the average and standard deviation for ten simulations.The hybrid method achieved a higher solution accuracy than either of DW or SA.
Second, we evaluated the solution accuracy of the preprocessing dependency.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 in Table I. Figure 3 shows the energy densities and the number of loops with the different preprocessings.Regardless of the solution accuracy of the preprocessing, the hybrid method improved the solutions from that of preprocessing.Moreover, even for the preprocessing with an extremely low solution accuracy (e.g., Random), the hybrid method gave a certain level of solution accuracy.However, a lower preprocessing solution accuracy required more loops to improve the solution.
Next, we performed the hybrid method with different parameter patterns to evaluate the effect of the parameters on Table I.Parameters of SA and hybrid optimization methods.N I is the number of solution samples to be obtained at the sampling step, N S is the number of reference original Ising models for generating sub-Ising models at the conversion step, N E is the number of newly obtained solutions through repetition of the conversion and execution steps, and N L is the number of repeats of the same X best to end the hybrid method.These parameters are described in Sect.  the solution accuracy.Figure 4 shows the energy density and the number of loops for N L set from 1 to 5, while the outer loop was set to 50.All other parameters were the same as those in Table I.Note that the solution pool sets obtained by the preprocessing SA were for all conditions of N L .In addition, Fig. 4 shows the improvement percentage in X best of the preprocessing SA among the ten simulations of the hybrid method.The solutions were improved for N L ≥ 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 L increased, both the percentage of improved solutions and the number of loops increased.There are two possible reasons for the increased number of loops.First, the hybrid method requires more loops prior to termination.Second, the solutions improve with an increase in N L .
Figure 5 shows the energy density and the number of loops with some sets of N E and N S .Note that the sets of solution pools obtained by the preprocessing SA were the same for all conditions of N S and N E , and the outer loop was set to 50.All other parameters were the same as those in Table I.The hybrid method improved the solution more when N E was 10 or 20 compared to when N E was 5.This is because a larger number of N 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 S was not significant when N E was the same.The number of loops was higher for N E of 10 and 20, which is where an improvement is observed.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 the number of loops for the sub-Ising model sizes (m) of 16, 40, 80, 120, and 144.Note that the sets of solution pools obtained by the preprocessing SA were for all sub-Ising model sizes.The outer loop was set to 50.All other parameters were the same as those in Table I.Constructing sub-Ising models with extremely small or large numbers of spins (e.g., a sub-Ising model size = 40 or 144) relative to the size of the original  Ising model (n = 160) did 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 obtained.By contrast, when the sub-Ising model size was large, the influence of the solution accuracy of D-Wave Advantage became more significant.Assuming that the solution accuracy of D-Wave Advantage for the original Ising model of size 160 was low (shown in Fig. 2), the sub-Ising model size became larger, the solution accuracy decreased, and the solution did not improve.These results suggest that the solution accuracy of the hybrid method depends on the sub-Ising model size.

Conclusion and future work
We propose a hybrid optimization method, which utilizes the advantages of two types of Ising machines (non-quantumtype Ising machines and quantum annealers) Then we evaluated the performance of the hybrid method via simulations using SA as a non-quantum-type Ising machine and D-Wave Advantage as a quantum annealer.The hybrid method achieves a higher solution accuracy than the D-Wave Advantage and preprocessing SA alone.Regardless of the preprocessing accuracy, the hybrid method improves the solution.However, there is a trade-off between the preprocessing accuracy and the number of loops for solution improvement.
We also evaluated the effect of the hybrid method parameters on the solution accuracy.Increasing the N L value, which determines the termination condition of the hybrid method, improves the solution accuracy while simultaneously increasing the number of loops.Next, we assessed the solution accuracy of the hybrid method with multiple combinations of N E and N S .The solution is enhanced when N E is 10 or 20 compared to when N E is 5.However, the appropriate sub-Ising model spins should be selected because the hybrid method has a sub-Ising model dependency.
This method first constructs a solution pool with a nonquantum-type Ising machine.Then it solves N E sub-Ising models using a quantum annealer.Using the hybrid method, it is more beneficial to solve multiple of the reduced-spin-size sub-Ising models embedding the original Ising model size into D-Wave Advantage.If the sub-Ising models can obtaine solutions in parallel, the hybrid method should be accelerated using the concept of parallel quantum annealing. 56)This approach should be considered in the future.
We also evaluated the performance of the hybrid method on the original Ising model with a problem size that can be embedded in D-Wave Advantage.However, several methods 42,[45][46][47] have been proposed to solve larger sizes of problems because the problem size that can be embedded in the quantum annealer is limited. Terefore, the performance of the hybrid method should be evaluated in the future using an original Ising model that cannot be embedded in a quantum annealer.Although the performance evaluation of the hybrid method was conducted using simulations with SA and D-Wave Advantage, the hybrid method should be evaluated using an actual Ising machine in the future.Furthermore, a theoretical investigation on the origin of the improved solution accuracy by the proposed method is insufficient in this paper.Since a theoretical investigation should contribute to parameter setting and performance improvement, it remains as a future work.

Fig. 1 .
Fig. 1. (Color online) Scheme of hybrid optimization method.The procedure consists of four steps: the sampling step, conversion step, execution step, and evaluation step.Gray, blue and red solid circles indicate that the spin has undetermined value, −1 and +1.Black arrows represent the flow.The sub-Ising model size (m), N I , N S , N E , and N L are the parameters of the method.

Fig. 2 .
Fig. 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 method.Error bars are the standard deviations obtained in ten runs.

Fig. 3 .
Fig. 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 obtained in ten runs.

Fig. 4 .
Fig. 4. (Color online) Performance of hybrid optimization methods with different N 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 obtained in ten runs.

Fig. 5 .
Fig. 5. (Color online) Performance of hybrid optimization methods with different sets of N E and N S .(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.(Bottom) Number of loops by the hybrid optimization method.Error bars are standard deviations.

Fig. 6 .
Fig. 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 method.Error bars are standard deviations.