Multi-Spin-Flip Engineering in an Ising Machine

A merge process is proposed to engineer a multi-spin flip in an Ising machine. The merge process deforms the Hamiltonian (energy function) of the Ising model. We prove a theorem for the merge process and show that a single-spin flip in the deformed Hamiltonian is equivalent to a multi-spin flip in the original Hamiltonian. A merge process induces a transition within the subspace of feasible solutions. We propose a hybrid simulated annealing (SA) algorithm with the merge process. The hybrid algorithm outperforms the conventional SA algorithm, genetic algorithm, and tabu search in the binary quadratic knapsack problems (QKP) and the quadratic assignment problems (QAP). Finally, the hybrid merge process is used in a real Ising machine. The performance is improved in QKP and QAP instances. The merge process is generally applicable to existing Ising machine hardware because the deformed Hamiltonian keeps the format of the Ising model.


INTRODUCTION
I SING machines have attracted attention both in academia and industry as efficient solvers for combinatorial optimization problems. A combinatorial optimization problem finds a set of decision variables to minimize or maximize an objective function subject to some constraints [1]. Typical examples found in textbooks are the travel salesman problem, knapsack problem, and graph partitioning. Additionally, combinatorial optimization problems are ubiquitous in daily life. Examples include the job scheduling problem, shift-planning problem, and traffic route optimization. When an Ising machine is used, a combinatorial optimization problem is mapped to an Ising problem. Ising problems adopt the format of the Ising model [2]. The Ising model has a variable called a spin, which takes a value of 1 or À1. The solution space in a combinatorial optimization problem is mapped to a spin configuration space fÀ1; 1g N where N is the number of spins. An Ising problem finds the spin configuration to minimize an energy function called a Hamiltonian, which is defined in the Ising model. Efficient methods have been proposed to map the Ising problem for diverse situations such as typical NP-complete problems [3], [4], machine learning [5], portfolio optimization [6], route optimization [7], [8], and material design [9]. Most Isingmachine hardware has been developed based on the operational principle of simulated annealing (SA) [10], [11], [12], [13], [14], [15].
It is still challenging to enhance the performance of an Ising machine to obtain the optimal solution of an Ising problem in a feasible time. The performance deterioration is attributed to an energy landscape of the Hamiltonian for an Ising problem (see Fig. 1a). When a combinatorial optimization problem with constraints is mapped to an Ising problem, the Hamiltonian H is given as H ¼ H obj þ H c , where H obj and H c encode the objective function and the constraints, respectively. H c ¼ 0 when the solution satisfies the constraints. Otherwise, H c > 0. The solutions satisfying the constraints are called feasible solutions (FSs). The FSs have the lower values of H (energies) than the infeasible solutions. Most of Ising-machine hardware adopts a single-spin flip as an operational principle to search for the optimum solution [10], [11], [12], [15]. Here, a flip is denoted by changing a spin value from À1 to 1 or from 1 to À1. A multi-spin flip is usually necessary to induce a transition among the FSs (see Fig. 1a). However, existing Isingmachine hardware has not implemented a multi-spin flip as an operational principle owing to the huge hardware cost. Thus, Ising machines should pass high-energy infeasible solutions to obtain the optimal solution. The transition to a high-energy solution is difficult to occur at low temperature, which deteriorates the performance of Ising machines. The problem is more serious for the daily-life combinatorial optimization problems because they have various kinds of constraints.
There is a previous study on engineering a multi-spin flip [16]. The SA performance was enhanced because simultaneous flips of multi spins can induce transitions within the subspace of FSs [16]. However, implementing the multispin flips into existing real Ising-machine hardware is quite difficult because it changes the operational principle of the Ising machine. This study overcomes this difficulty by providing a merge process that deforms the Hamiltonian of the Ising model. As far as we know, no work has done on deforming the energy landscape of the Ising model.
In this study, we propose a merge process to solve the problem of Ising machines. A merge process engineers a multi-spin flip by deforming the Hamiltonian of the Ising model. An appropriate choice of a merge process can induce a transition within the subspace of FSs by a single-spin flip. Fig. 1b provides a schematic picture of a merge process for a two-spin flip. The merge process encodes the two-spin values ðs 1 ; s 2 Þ ¼ ð1; À1Þ and ðÀ1; 1Þ into a spin value t ¼ 1 and À1, respectively. A single-spin flip from t ¼ 1 to À1 in the deformed Hamiltonian corresponds to a two-spin flip from ðs 1 ; s 2 Þ ¼ ð1; À1Þ to ðÀ1; 1Þ in the original Hamiltonian. The deformation keeps the energies the same in a spin-configuration subspace (see theorem 1 and corollary 2 in Sec. 3), and thus the optimization for a deformed Hamiltonian efficiently searches the optimal solution for the original Hamiltonian. Here, we provide an SA algorithm hybridized with a merge process. The algorithm is applied to two combinatorial optimization problems with different types of constraints: the binary quadratic knapsack problem (QKP) and the quadratic assignment problem (QAP), both of which are categorized into NP-hard problems [17], [18]. The hybrid algorithm outperforms the conventional SA algorithm, genetic algorithm (GA), and tabu search (TS). Additionally, we hybridize a merge process into an algorithm of an existing Ising machine to solve QKP and QAP. The merge process is applicable to existing Ising-machine hardware since the deformed Hamiltonian keeps the format of the Ising model. The hybrid algorithm shows superior performances.
The main contributions of this paper are summarized as follows.
We propose a merge process that engineers a multispin flip by deforming the Hamiltonian of an Ising model. A single-spin flip in the deformed Hamiltonian is equivalent to a multi-spin flip in the original Hamiltonian. We prove a strong theorem for the merge process. We propose hybrid merge processes with SA or an algorithm of Ising-machine hardware. Engineering a multi-spin flip efficiently searches an optimum in Ising-machine hardware, which is based on singlespin-flip operational principle. We apply the hybrid merge processes to QKP and QAP as examples of difficult combinatorial optimization problems. The hybrid merge process shows superior performances in all problem instances and reduces a residual energy by an average of 86% and a maximum of 93% compared to the conventional SA. When applying to Ising-machine hardware, it reduces a residual energy by an average of 50% and a maximum of 71% (see Tables 3, 4, 5, and 7 in detail). The rest of this paper is organized as follows. Section 2 introduces SA, GA, and TS as metaheuristics on Ising models. Herein, we briefly explain SA as an operational principle of Ising machines. Section 3 introduces the merge process. Section 4 hybridizes the merge process into SA or an algorithm of an Ising machine. Section 5 shows the experimental results. Section 6 concludes this paper and provides future outlooks.  [19], [20], [21], GA [22], ant colony optimization [23], TS [24], and iterated local search [25]. Among them, SA, GA, and TS have been used to efficiently search the lowest-energy solution of the Ising models [26], [27], [28], [29]. Metaheuristics search the optimal solution of an objective function in combinatorial optimization problems. Here, an Ising model is adopted as the objective function since it is the common format for Ising machines. The Hamiltonian of an Ising model is defined as where s i 2 fÀ1; 1g is the spin, N is the number of spins, and H 0 is a constant. Here, J i;j 2 R and h i 2 R denote the interaction between spins i and j and the bias on spin i, respectively. The interaction coefficients are symmetric (i.e., J i;j ¼ J j;i ) and the diagonal elements are set to 0 (i.e., J i;i ¼ 0). The Ising problem is formulated to find the spin configuration fs i g NÀ1 i¼0 that minimizes the value of H.

Simulated Annealing
SA is a single-solution-based metaheuristic algorithm which improves a single candidate solution by using a local search [26]. Algorithm 1 shows the SA algorithm implemented by the single-spin-flip Markov Chain Monte Carlo (MCMC) method. Most Ising-machine hardware is based on this algorithm. The initial solution has a random spin configuration. Then the spin configuration is repeatedly updated. Each iteration randomly selects a spin and generates a candidate solution where the chosen spin is flipped from the current solution. Next, consider the transition from the current solution to the candidate solution. The probability of making the transition is specified by a transition  probability W ðDE; T Þ, where T is the temperature and DE is the difference between the current-solution energy E and the candidate-solution energy E 0 . That is Two well-known choices of the transition probability are the heat-bath method and the Metropolis method [30]. SA gradually decreases the temperature. The convergence theorem ensures that SA obtains the lowest-energy solution of H when the temperature decreases at a sufficiently slow rate [31].

Genetic Algorithm
GA is a population-based metaheuristic algorithm which utilizes multiple candidate solutions during the search process [27], [29]. The population in each generation is fixed to n p . The first generation consists of random spin configurations. The population is repeatedly replaced by using genetic operators on the basis of their energies in Eq. (1). The genetic operators are crossover and mutation. The crossover operator splits two solutions (parents) in the same manner and combines different regions to generate a new solution (child). The mutation operator randomly selects spins in a solution (parent) and generates a new solution (child) by randomizing the selected spin values. The number of new solutions which are generated by the crossover and mutation operators are denoted by n crossover and n mutation , respectively.

Tabu Search
TS is a single-solution-based mataheuristic algorithm which relies on the use of memory structures called tabu list [26], [28]. As in SA, TS iteratively updates a solution using a local search. When considering a transition to a candidate solution, tabu list forbids or penalizes certain transitions that would return to a recently encountered solution. Tabu list is specified by a set of attributes that change in going from the current solution to the next (e.g., the spin configuration and the value of a flipped spin). Tabu tenure n t denotes the size of tabu list (i.e., the number of recently encountered solutions to be referred).

MERGE
To enhance the performance of an Ising machine, a mechanism of inducing a transition among the FSs without passing an infeasible solution is necessary. Herein, we introduce a merge process to provide the mechanism. The merge process deforms the Hamiltonian of an Ising model.

Theorem
We newly prove a strong theorem for the merge process.
The following theorem gives a deformed Hamiltonian H m ðfm i g; fs i g; HÞ.
where d i;j denotes the Kronecker delta. For a given spin configuration fs i g NÀ1 i¼0 , consider the spin configuration fs i g NÀ1 i¼0 wherẽ s i ¼ s i s m i . Then, H m ðfm i g; fs i g; HÞ in the spin configuration fs i g NÀ1 i¼0 has the same energy as H in the spin configuration fs i g NÀ1 i¼0 . H m ðfm i g; fs i g; HÞ is independent of s i when m i 6 ¼ i. H m ðfm i g; fs i g; HÞ The interaction term is rewritten as s k s ' J k;' s m k s m ' ; s k s ' J k;' d m k ;i d m ' ;j s i s j ; In the derivation of Eq. (6), from the first line to the second line and from the fourth line to the fifth line, we have used the symmetry of J i;j (i.e., J i;j ¼ J j;i ). The bias term is rewritten as Substituting Eqs. (6) and (7) into Eq. (5) and comparing it with Eq. (3) gives fJ m i;j g, fh m i g, and H m 0 in Eq (4).
Theorem 1 indicates that the energies for a deformed Hamiltonian are lower bounded by the lowest energy of H. The relationss i ¼ s i s m i give a correspondence between a solution in H and a solution in H m ðfm i g; fs i g; HÞ. Thus, searching an optimum for a deformed Hamiltonian provides a high-quality solution for H. Existing Ising machines can be used to search the optimum, since the deformed Hamiltonian has the format of an Ising model. When the energy in a deformed Hamiltonian is identical to the lowest energy of H in a spin configuration of fs i g NÀ1 i¼0 ,s i ¼ s i s m i gives an optimum for H.
In the following, spin i is called a merged spin when m i 6 ¼ i, while an unmerged spin when m i ¼ i. When m i 6 ¼ i, spin m i is called the merging spin of spin i. The merging spin is always an unmerged spin because m m i ¼ m i holds. Theorem 1 indicates that a deformed Hamiltonian is independent of the merged spins. If Then, we obtain the following corollary to theorem 1.
Corollary 2. H m ðfm j g; fs j g; HÞ has the same energy spectrum as H in the spin configuration subspace satisfying s i ¼ s i s m i for all i 2 f0; 1; . . . ; N À 1g. Namely

Merge
The merge process engineers a multi-spin flip based on theorem 1 and corollary 2. Suppose that a current solution satisfies s i ¼ s i s m i for all spins. Note that s i ¼ 1 for unmerged spins. Then theorem 1 gives a deformed Hamiltonian H m ðfm i g; fs i g; HÞ and corollary 2 implies H m ðfm i g; Below, simple Ising models demonstrate that merge processes for a two-spin flip and a four-spin flip induce a transition within the subspace of FSs.

Merge Process for a Two-Spin Flip
We consider a two-spin Ising model, which is expressed as This model has a one-hot constraint for fs 0 ; s 1 g (i.e., one spin value is 1 and the other is À1 in the set fs 0 ; s 1 g). Onehot constraint appears in problems like graph coloring and clique cover [3], [4]. Fig. 2a shows the energy landscape of H 1 . The FSs are the optimum at ðs 0 ; s 1 Þ ¼ ðÀ1; 1Þ with an energy of À2 and a local minimum at ðs 0 ; s 1 Þ ¼ ð1; À1Þ with an energy of 2. Here, the local minimum is defined as the solution where a single-spin flip does not lower the energy. The energy required to move from a local minimum to another one by single-spin flips is called energy barrier. The energy barrier between a local minimum ðs 0 ; s 1 Þ ¼ ð1; À1Þ and the optimum ðs 0 ; s 1 Þ ¼ ðÀ1; 1Þ is 38. A two-spin flip is necessary to induce a transition between the FSs. Suppose that the solution is in a local minimum of H 1 . To engineer a two-spin flip, we set spin 1 as a merged spin and spin 0 as the merging spin. Then, Here, the deformed Hamiltonian H 0 1 is independent of s 1 since spin 1 is a merged spin. Fig. 2b shows the energy landscape of H 0 1 . This model has optimums at ðs 0 ; s 1 Þ ¼ ðÀ1; 1Þ and ðÀ1; À1Þ with an energy of À2 and degenerate solutions with an energy of 2 at ðs 0 ; s 1 Þ ¼ ð1; 1Þ and ð1; À1Þ. The energy barrier in H 1 disappears in H 0 1 . Consequently, the solution ðs 0 ; s 1 Þ ¼ ð1; À1Þ is no longer a local minimum in H 0 1 and a transition to an optimum ðs 0 ; s 1 Þ ¼ ðÀ1; À1Þ can occur by a single-spin flip of spin 0. The energy of ðs 0 ; s 1 Þ ¼ ðÀ1; À1Þ in H 0 1 is À2, which is identical to the lowest energy in H 1 . Then, theorem 1 gives the optimum in H 1 by updating the merged-spin value s 1 by s 1 s m 1 ¼ Às 0 . Here, spin 0 in an optimum of H 0 1 takes a value of À1. Thus, spin 1 is updated as Overall, the merge process induces the transition between the FSs of H 1 with a single-spin flip of s 0 . Namely, a single-spin flip in H 0 1 is equivalent to a two-spin flip in H 1 . The merge process for a two-spin flip is effective when k-hot constraint (k 2 N) exists. One-hot constraint is a special case of k-hot constraint with k ¼ 1. k-hot constraint appears in problems like graph partitioning [3]. The constraint Hamiltonian is given by H 1 contains H 2 with n ¼ 2 and k ¼ 1. The solutions where k spins take a value of 1 are the FSs. At least a two-spin flip is necessary to induce a transition among the FSs. Suppose solution A two with s i ¼ 1 and s j ¼ À1 and consider solution B two where spin i is flipped from solution A two and solution C two where both spins i and j are flipped from solution A two . When solution A two is a FS, solution B two is an infeasible solution but solution C two is a FS. To engineer a two-spin flip from solution A two to solution C two , we set spin j as a merged spin and spin i as the merging spin of spin j. Namely, m j ¼ i and s j ¼ s j s m j ¼ À1 according to corollary 2. The other spins are set as unmerged spins (i.e., m k ¼ k and s k ¼ 1). Then, theorem 1 gives the deformed Hamiltonian H 0 Hence, a transition from solution A two to solution B two is possible without increasing the energy in H 0 2 . Then, updating the merged spin j as s j ¼ s j s m j ¼ 1 transforms solution B two into solution C two . In this way, the merge process induces a transition within the subspace of FSs.

Merge Process for a Four-Spin Flip
We consider a four-spin Ising model, which is given as This model has one-hot constraints for fs 0 ; s 1 g, fs 2 ; s 3 g, fs 0 ; s 2 g, and fs 1 ; s 3 g. The FSs are the optimal at ðs 0 ; s 1 ; s 2 ; s 3 Þ ¼ ð1; À1; À1; 1Þ with an energy of À1 and a local minimum at ðs 0 ; s 1 ; s 2 ; s 3 Þ ¼ ðÀ1; 1; 1; À1Þ with an energy of 1. The remaining infeasible solutions have larger energies. A four-spin flip is necessary to induce a transition between the FSs. Suppose that the solution is in a local minimum of H 3 and consider a merge process to engineer a four-spin flip. Let us set spins 1, 2, and 3 as merged spins and spin 0 as their merging spin. Then, The solution ðs 0 ; s 1 ; s 2 ; s 3 Þ ¼ ðÀ1; 1; 1; À1Þ is no longer a local minimum in H 0 3 . A single-spin flip of spin 0 can induce a transition to an optimum in H 0 3 (i.e., ðs 0 ; s 1 ; s 2 ; s 3 Þ ¼ ð1; 1; 1; À1Þ). The energy of ðs 0 ; s 1 ; s 2 ; s 3 Þ ¼ ð1; 1; 1; À1Þ in H 0 3 is À1, which is identical to the lowest energy in H 3 . Then, theorem 1 gives the optimum in H 3 by updating the merged-spin values as s i s m i (i.e., s 1 ¼ s 1 s m 1 ¼ À1, s 2 ¼ s 2 s m 2 ¼ À1, and s 3 ¼ s 3 s m 3 ¼ 1). The solution ðs 0 ; s 1 ; s 2 ; s 3 Þ ¼ ð1; À1; À1; 1Þ is the optimum in H 3 . Thus, the merge process induces a transition between the FSs of H 3 with a single-spin flip of s 0 . Namely, a single-spin flip in H 0 3 is equivalent to a four-spin flip in H 3 . The merge process for a four-spin flip is effective when k-hot constraints (k 2 N) interact with each other. The constraint Hamiltonian is given by H 4 appears in problems like travel salesman problem and QAP [3], [4].

HYBRID ALGORITHMS
We hybridize a merge process into SA or the algorithm for an Ising machine, which is abbreviated as IM. The hybrid algorithms are called Merge SA and Merge IM, respectively. It is noted that the merge process cannot be easily hybridized into GA and TS. The algorithm that incorporates a merge process should meet two properties: The first is a single-solution-based algorithm. This is necessary to give a single deformed Hamiltonian. The other is the Markov property. This is because the information in the past search history is meaningless to the newly generated deformed Hamiltonian. SA and IM have both properties, but GA and TS do not (that is, GA is population based and TS includes a non-Markovian memory).
Algorithm 2 shows Merge SA and Merge IM. The spin values are repeatedly updated according to single-spin-flip SA or the guiding principle of Ising-machine hardware for the deformed Hamiltonian H m ðfm j g; fs j g; HÞ. Merge interval gives the frequency of recreating the deformed Hamiltonian. Before the recreation, the merged-spin values are updated according to theorem 1 (i.e., s j ¼ s j s m j ). The deformed Hamiltonian is specified by fm j g NÀ1 j¼0 and fs j g NÀ1 j¼0 . The process to assign them should depend on combinatorial optimization problems. The assignment processes for QKP and QAP are discussed in the following section (see Sec. 5.2). When implementing Merge IM, repeated communications between a central processing unit (CPU) and an Ising machine are necessary. In Algorithm 2, lines 3-5 and lines 12-16 are performed on the CPU, while lines 7-10 are performed on an Ising machine. In each communication, the CPU sends the QUBO matrix for the deformed Hamiltonian, the spin configuration in the current solution, and the current temperature to the Ising machine. Then, the CPU receives an updated spin configuration from the Ising machine. Table 1  The overall time complexity for a merge process is given by OðN 2 Þ. The time complexity for creating the deformed Hamiltonian is OðN 2 Þ (lines 4-5 or lines 15-16 in Algorithm 2) since the deformed Hamiltonian has NðN À 1Þ=2 interactions and N biases. The time complexity for updating the merged-spin values (lines 12-14 in Algorithm 2) is proportional to the number of merged spins, and thus it is OðNÞ. Developing a massive parallel computing design using the graphics processing unit (GPU) or field programmable gate array (FPGA) is of great interest, but it is out of the scope of this study.

EXPERIMENTAL RESULTS
This section compares the algorithm performance with and without a merge process. Section 5.1 gives Ising problem formulations of QKP and QAP. They are NP-hard problems [17], [18]. It has been reported that existing Isingmachine hardware is unable to provide even a near-optimal solution of them [32], [33]. Section 5.2 shows the assignment process of fm i g and fs i g for each combinatorial optimization problem. Section 5.3 describes the simulation. Section 5.4 provides the comparison methods. Section 5.5 compares SA, GA, TS, and Merge SA, while Section 5.6 compares IM and Merge IM.

Ising Problem Formulations of QKP and QAP
Here, QKP and QAP are formulated as Ising problems. To obtain the Ising model for each combinatorial optimization problem, we introduce its equivalent, which is referred to as the quadratic unconstrained binary optimization (QUBO). QUBO is defined as where x i 2 f0; 1g is a binary variable and Q 0 is a constant. Here, Q i;j 2 R is regarded as the ði; jÞ-element of a matrix, which is called the QUBO matrix.

Binary Quadratic Knapsack Problem
QKP finds a set of items to yield the highest profit within the capacity of the knapsack [34]. QKP has applications, for example, in very-large scale integrated circuit design [35] and compiler design [36]. A binary variable x i is introduced to show whether item i is chosen. Item i is chosen if x i ¼ 1.
If not, then x i ¼ 0. The overall profit within the knapsack is given by where n item is the number of items. The matrix fp i;j g is referred to as the profit matrix [32]. The overall weight of items within the knapsack is given by where w i is the weight of item i. QKP is formulated to maximize P ðfx i gÞ subject to an inequality constraint W ðfx i gÞ C, where C is the capacity of the knapsack. This study uses QKP instances in Refs. [37], [38]. The profit matrix density r is chosen from r 2 f0:25; 0:5; 0:75; 1g for n item ¼ 200 and from r 2 f0:25; 0:5g for n item ¼ 300. For each n item and r, ten instances are randomly generated. Problems are named as n d i where n ¼ n item , d ¼ 100r, and i specifies the seed of an instance [39]. QUBO for a QKP has the following form [40] A QKP > 0 is the constraint coefficient and fy i g mÀ2 i¼0 are the auxiliary binary variables. The constraint QUBO in Eq. (19) is 0 only when the overall weight is smaller than C but larger than C À m þ 1. We set m to the largest weight of items according to [32]. The Ising model for a QKP is obtained via the transformation of the binary variables fx i g n item À1 i¼0 and fy i g mÀ2 i¼0 into spins [6], which gives the Hamiltonian H QKP with N ¼ n item þ m À 1.

Quadratic Assignment Problem
QAP assigns facilities to different locations such that the total transport cost is minimized [41]. The total transport cost is defined as where p gives a permutation and n fac are the numbers of facilities and locations. Here, f i;j denotes the flow amount between facilities i and j, while d a;b denotes the distance between locations a and b. d a;b is symmetric (i.e., d a;b ¼ d b;a ). QAP is formulated to find the permutation that minimizes SðpÞ. This study uses QAP instances in Ref. [42] for n fac 30 and in Ref. [43] for n fac ! 40. f i;j is symmetric for n fac 30 and asymmetric for n fac ! 40.
QUBO for a QAP is given as where Q QAP obj ðfx i;a gÞ and Q QAP c ðfx i;a gÞ encode the objective function and the constraint for the permutation, respectively. A QAP > 0 is the constraint coefficient. A binary variable x i;a is introduced to show the locations of the facilities. When x i;a ¼ 1, facility i is placed at location a. Then, QUBO for the objective function and constraint is respectively given by [3], [4] Q QAP obj ðfx i;a gÞ ¼ Q QAP c ðfx i;a gÞ ¼ X The constraint QUBO is 0 if and only if each facility is placed at one location and each location has one facility. When the constraint is satisfied, the permutation is described as x i;a ¼ 1 when a ¼ pðiÞ. Otherwise, x i;a ¼ 0.
The Ising model for a QAP is obtained via the transformation of the binary variables fx i;a g into spins [6], which gives the Hamiltonian H QAP with N ¼ n fac 2 . Below, the spin corresponding to x i;a is represented by s in fac þa .

Assignment
Processes of fm i g and fs i g for QKP and QAP Here, processes to assign fm i g and fs i g for QKP and QAP are provided. The assignment gives a deformed Hamiltonian in Merge SA and Merge IM (see Algorithm 2). In QKP, a two-spin flip is effective to replace item i inside the knapsack by item j outside the knapsack when both items cannot be selected within the capacity of the knapsack. To implement the process by using single-spin flips, removing item i from the knapsack or adding item j into the knapsack is necessary. The solutions where both items are selected or neither item is selected have larger energies than the solutions where either item is selected. Thus, a transition to a high-energy solution is necessary. On the other hand, a two-spin flip allows to replace the items without passing a high-energy solution. Algorithm 3 provides the pseudo code to assign fm i g and fs i g in QKP. Spins are initially set as unmerged spins (i.e., m i ¼ i and s i ¼ 1). Then fm i g and fs i g are repeatedly updated. A probability p QKP determines whether spin i is chosen as a merged spin. Then, the merging spin j is randomly selected from unmerged spins (i.e., m i ¼ j). Corollary 2 gives s i ¼ s i s m i ¼ s i s j . This assignment engineers a two-spin flip of spins i and j by a single-spin flip of spin j in the deformed Hamiltonian. The merged spins and merging spins are selected from the spins for items.
In QAP, the Hamiltonian H QAP contains H 4 in Eq. (14) with n ¼ n fac and k ¼ 1. Thus, engineering a four-spin flip induces a transition among the FSs of H QAP without passing an infeasible solution. Algorithm 4 provides the pseudo code to assign fm i g and fs i g in QAP. Spins are initially set as unmerged spins, and then fm i g and fs i g are repeatedly updated. A probability p QAP determines whether a spin in the set of spins fs in fac þa g n fac À1 a¼0 is chosen as a merged spin. Then, randomly select four unmerged spins that satisfy ðs in fac þa ; s in fac þb ; s jn fac þa ; s jn fac þb Þ ¼ ð1; À1; À1; 1Þ for i 6 ¼ j and a 6 ¼ b. Here, spin jn fac þ b is chosen as the merging spin and the other three spins are chosen as the merged spins. Corollary 2 gives fs k g for the merged spins (e.g., s in fac þa ¼ s in fac þa s m in fac þa ¼ 1). This assignment engineers a four-spin flip of spins in fac þ a, in fac þ b, jn fac þ a and jn fac þ b by a single-spin flip of spin jn fac þ b in the deformed Hamiltonian. Both algorithms require a post process to satisfy the condition of m i that m m i ¼ m i . If spin m i is a merged spin (i.e., m m i 6 ¼ m i ), then m i and s i are updated by m m i and s m i s i , respectively. This update keeps the relation of s i ¼ s i s m i , which is a condition for corollary 2. m i and s i are repeatedly updated until m m i ¼ m i holds.

Simulation Details
Here, simulation details for SA, GA, TS, Merge SA, IM and Merge IM are provided. The SA, GA, and Merge SA experiments were conducted on the ISSP supercomputer at Kashiwa city in Japan using Fortran 90 as the implementation language. The TS experiments used a solver in a software library. The IM and Merge IM experiments were performed on real Ising-machine hardware using Python 3.7.6 as the implementation language. Table 2 gives the parameter set for each algorithm. We applied SA to the Ising model by adopting the single-spin-flip MCMC method (see Algorithm 1). The number of iterations is set as 10 5 N.
Here, the heat-bath method is adopted as W ðDE; T Þ. The temperature interval is given by N. The temperature at the kth update is given as T k ¼ T i e Àk , where and T i are the decay rate and the initial temperature, respectively. The decay rate is determined so that the temperature at the end of the iterations (i.e., T 10 5 ) is the final temperature. The initial temperature is larger than the maximum of DE. Here, E h :¼ max i ½ð P NÀ1 j¼0 jJ ij jÞ þ jh i j is used since it gives an upper bound of DE. The initial temperature is given by  a h E h , where a h takes a value of 1, 10, or 100. On the other hand, the final temperature is smaller than the energy scale of each combinatorial optimization problem. Here, E ' ¼ 1 is used for QKP and E ' ¼ min i;jði6 ¼jÞ f i;j min a;bða6 ¼bÞ d a;b is used for QAP. The final temperature is given by a ' E ' , where a ' takes a value of 1, 0.1, or 0.01. The number of runs is set to 20. We applied GA to the Ising model following Ref. [29]. The number of iterations for replacing the population is fixed to 10 3 1 . n p is set to 200 and the new population is generated through genetic operators, ðn crossover ; n mutation Þ ¼ ð50; 50Þ. In addition, we introduce 50 random solutions and retain the 50 lowestenergy solutions from the previous population. The uniform crossover operator [22] is used to generate a new solution (child) from two solutions (parents). Mutation randomly selects 0.1 N spins in the solution (parent) and randomizes their values according to probability p GA . p GA is the probability that a selected spin value of the solution (child) is 1. We set p GA ¼ 0:5 in QKP and p GA ¼ ðn fac Þ À1 in QAP 2 . p GA is also used to generate the 50 random solutions (i.e., the number of ðþ1Þ-spins in a random solution is p GA N on average). The crossover and mutation select parents according to the roulett-wheel selection given in Ref. [29]. The number of runs is set to 20.
TS experiments used a solver in [44], which implements a TS algorithm in Refs. [26], [28]. The solver runs up to approximately 1000 spins. The input parameters such as tabu tenure are set to the default values (e.g., n t ¼ 20). The number of runs is set to 20.
We used Digital Annealer as an Ising machine [10]. The solver is based on the single-spin-flip SA. The Ising-machine hardware has a maximum of 8192 spins on a complete graph. The hardware constraint gives the number of iterations. The parameters of the Ising machine are set similar to those for the SA (see Table 2). In Merge SA and Merge IM, merge interval is given as MN (M 2 N). The probabilities p QKP and p QAP assume values of 0.1, 0.2, or 0.3. Merge SA and Merge IM have 20 and 5 runs, respectively. The other parameters are set similar to those for the SA.

Comparison Methods
We used the procedure detailed in Ref. [40] to compare the performance among SA, GA, TS, and Merge SA in QKP and QAP. We calculated the FS rate and the residual energy for  Number of runs  20 20 20 20  20  5 1. GA calculates the energy to generate a population, while SA calculates the energy difference to update a solution. The calculation of energy is computationally much more expensive than that of energy difference [26]. Here, the number of iterations in GA is smaller than that in SA so that the total complexities of both algorithms are the same.
2. In QAP, the number of ðþ1Þ-spins in FSs is n fac , and thus p GA should be set as the density of ðþ1Þ-spins (i.e., p GA ¼ n fac =N). In QKP, the number of ðþ1Þ-spins in FSs is not known in advance, and thus p GA ¼ 0:5 is adopted. each algorithm. The FS rate denotes the ratio of the number of obtained FSs to the number of runs. FSs are given as solutions that satisfy Q QKP c ðfx i g; fy i gÞ ¼ 0 in QKP and Q QAP c ðfx i;a gÞ ¼ 0 in QAP. The residual energy was calculated for the FSs. The residual energies in QKP and QAP are respectively defined as E QKP res ¼ P Ã À P ðfx i gÞ; where P Ã denotes the maximum of P ðfx i gÞ obtained in Ref. [39] and S Ã denotes the minimum of SðpÞ found in Ref. [45]. The smaller the value of the residual energy, the better the performance. The FS rates and residual energies typically increase with the constraint coefficients (i.e., A QKP and A QAP ) [6]. The dependences indicate that the constraint coefficient has an optimal value, where the FS rate is high and the residual energy is small. To systematically determine the optimum value of the constraint coefficient, we set the threshold value of the FS rate as r th ¼ 0:9. We repeatedly solved each problem instance while varying the hyperparameters: A QKP or A QAP , a h , and a ' in SA, A QKP or A QAP in GA and TS, and A QKP or A QAP , p QKP or p QAP , a h , and a ' in Merge SA. We study the parameter regimes of 1 A QKP 500 and 10 A QAP 5000, and set the precision threshold of A QKP and A QAP to 1 and 10, respectively. Then we determined a parameter set that yields a FS rate above the threshold value. Finally, we searched for the optimal set to minimize the residual energy averaged over the FSs. Below, we denote the residual energies averaged over FSs and the minimal residual energies at the optimal values of hyperparameters in QKP and QAP by E QKP res;ave , E QAP res;ave and E QKP res;min , E QAP res;min , respectively. The performances of SA, GA, TS, and Merge SA were compared using the average residual energies and the minimal residual energies.
We used the following method to compare the performances of IM and Merge IM in QKP and QAP. The hyperparameters for IM (i.e., A QKP or A QAP , a h , and a ' ) were determined using the same method as those for SA, GA, TS, and Merge SA. Namely, we found a parameter set that yields a FS rate above r th and searched for the optimal set to minimize E QKP res;ave or E QAP res;ave . Here, r th is set to 0.4 in QKP and 0.9 in QAP. Merge IM used the optimal set of hyperparameters in IM and searched for the optimal value of p QKP or p QAP that minimizes E QKP res;ave or E QAP res;ave . When the FS rate is less than r th for any value of p QKP or p QAP , the constraint coefficient is increased until a parameter set that yields a FS rate above r th is found.
Merge SA uses M ¼ 1 to compare the performance with SA, GA, and TS since a small M shows a better performance (see Fig. 3). On the other hand, M is set to 10 3 in Merge IM. It is difficult to set M to a smaller value in Merge IM because each communication takes the order of one minute (see Table 1).
Tables 3, 4, and 6 give the optimal set of hyperparameters for each QKP instance. A QKP increases with n item and r. For the given values of n item and r, A QKP in TS is the largest, and A QKP in Merge SA is typically larger than those in SA and GA. a h and a ' are problem dependent. By contrast, p QKP in Merge SA and Merge IM is nearly independent of n item and r on average. It can be approximated as 0.2.
Tables 5 and 7 provide the optimal set of hyperparameters for each QAP instance. A QAP increases with n fac . For each n fac , A QAP are almost the same for all the algorithms except for GA and TS. A QAP in GA and TS are larger than the others. a h and a ' are problem dependent in SA, while a h =a ' takes a small value (i.e., 1 or 10) in Merge SA and IM, indicating a better performance as the temperature slowly decreases. p QAP takes a value of 0.2 or 0.3 in both Merge SA and Merge IM.
GA and TS do not yield a FS rate above the threshold value within the parameter regime of A QKP in some QKP instances (e.g., instance 300 25 8 in GA). GA does not yield a high FS rate in QAP instances with large n fac . TS experiments of QAP were performed up to n fac ¼ 30 due to the size limitation of the solver [44].  residual energy in Merge SA is approximately 12% of that in SA. SA, GA, and TS give larger E QKP res;min than Merge SA for all 60 instances, except for instance 200 25 3 of SA. Table 5   Merge IM for all instances, except for instance 300 25 5. In QAP, IM provides larger values of E QAP res;ave and E QAP res;min than Merge IM in all instances. The enhanced performances indicate that merge processes engineer two-spin flips for QKP and four-spin flips for QAP in real Ising-machine hardware, which is based on a single-spin-flip operational principle. The M-dependences in Merge SA (see Fig. 3) imply that Merge IM should show a better performance when M takes a smaller value.

CONCLUSION
This study proposed a merge process that engineers a multi-spin flip in Ising-machine hardware. We proved a theorem for the merge process and provided two types of merge processes for a two-spin flip and a four-spin flip. Simple Ising models demonstrated that a merge process induces a transition within the subspace of FSs. We also provided an SA algorithm hybridized with a merge process. The hybrid algorithm outperforms the conventional SA algorithm, GA, and TS in QKP and QAP in terms of a residual energy. Finally, an enhanced performance of the algorithm of an Ising machine hybridized with a merge process is demonstrated in QKP and QAP.
A merge process may enhance the performance of various types of existing Ising machines. In the future, we plan to incorporate a merge process into them and confirm its effectiveness.