Spin-Variable Reduction Method for Handling Linear Equality Constraints in Ising Machines

We propose a spin-variable reduction method for Ising machines to handle linear equality constraints in a combinatorial optimization problem. Ising machines including quantum-annealing machines can effectively solve combinatorial optimization problems. They are designed to find the lowest-energy solution of a quadratic unconstrained binary optimization (QUBO), which is mapped from the combinatorial optimization problem. The proposed method reduces the number of binary variables to formulate the QUBO compared to the conventional penalty method. We demonstrate a sufficient condition to obtain the optimum of the combinatorial optimization problem in the spin-variable reduction method and its general applicability. We apply it to typical combinatorial optimization problems, such as the graph <inline-formula><tex-math notation="LaTeX">$k$</tex-math><alternatives><mml:math><mml:mi>k</mml:mi></mml:math><inline-graphic xlink:href="shirai-ieq1-3239539.gif"/></alternatives></inline-formula>-partitioning problem and the quadratic assignment problem. Experiments using simulated-annealing and quantum-annealing based Ising machines demonstrate that the spin-variable reduction method outperforms the penalty method. The proposed method extends the application of Ising machines to larger-size combinatorial optimization problems with linear equality constraints.


I. INTRODUCTION
C OMBINATORIAL optimization problems find the optimal combination of decision variables to minimize or maximize an objective function under a set of given constraints [1]. Solving a combinatorial optimization problem with many decision variables is challenging because the number of solution candidates increases exponentially as the number of decision variables increases. Although typical examples found in textbooks are the satisfiability problem, quadratic assignment problem (QAP), and graph k-partitioning problem (GkPP), combinatorial optimization problems are ubiquitous in daily life. Examples include the logistics optimization, traffic route optimization, and drug discovery.
Ising machines are specialized computers for combinatorial optimization problems. Most Ising machine hardware operates on the basis of simulated annealing (SA) [2], [3], [4], [5], [6], Manuscript [7] and quantum annealing (QA) [8], [9] (see the review of Ising machines in [10]). The quadratic unconstrained binary optimization (QUBO) is used to solve a combinatorial optimization problem in Ising machines [11], [12], [13], [14], [15], [16]. The QUBO is described by binary variables called spins, which take values of 0 or 1. 1 Then, the solution space of a combinatorial optimization problem is transformed into the spin configuration space {0, 1} ⊗N , where N is the number of spins. The QUBO problem finds the spin configuration to minimize the energy function defined by the QUBO. Ising machines search for the lowest-energy state of a QUBO, which is called the ground state. Efficient methods have been proposed for improving the performance of Ising machines [17], [18], [19], [20]. The optimal solution of a combinatorial optimization problem is obtained from the ground state. This paper develops a method for Ising machines.
The energy function of a QUBO is generally given as the sum of the terms for an objective function and the constraints. Feasible solutions (FSs) satisfy the constraints. The conventional method, which is called the penalty method, gives the constraint term by increasing the energy of infeasible solutions [12]. The spin configuration space is separated into many FS subspaces with lower energies by the infeasible solution space with higher energies. Fig. 1(a) applies the penalty method to a combinatorial optimization problem with two binary variables and a single linear equality constraint. The QUBO is expressed by H penalty = x 1 − x 2 + 5(x 1 + x 2 − 1) 2 , where the first term x 1 − x 2 is the objective function and the second term (x 1 + x 2 − 1) 2 is the constraint term. The penalty 5 is given to the constraint term to increase the energy of infeasible solutions. Then, the two FSs (x 1 , x 2 ) = (0, 1) and (1, 0) are separated by infeasible solutions with large energies. The transition probability between the FSs is low when the thermal or quantum fluctuation is small. In this situation, searching for the optimum using an Ising machine is difficult. SA-based or QA-based Ising machine studies have reported that the penalty method falls into a local minimum solution and fails to reach the optimum of the QUBO (e.g., Refs. [21], [22]). Therefore, developing an efficient constrainthandling method is important.
In this study, we propose a new constraint-handling method called the spin-variable reduction (SVR) method. The proposed method, which is based on the idea of variable reduction, expresses a spin variable using other spin variables in terms Fig. 1. Examples of the penalty method and the spin-variable reduction method in a simple combinatorial optimization problem with two binary variables and a single linear equality constraint. Both methods formulate the combinatorial optimization problem into QUBOs. (a) The penalty method gives a QUBO H penalty as a function of (x 1 , x 2 ). Two infeasible solutions have a large energy. (b) The spin-variable reduction method gives a QUBO H svr as a function of x 1 . H svr is independent of x 2 and x 2 is obtained by x 1 using the variable relationship in the linear equality constraint (i.e., x 2 = 1 − x 1 ). of the variable relationships in linear equality constraints. The SVR method reduces the number of spin variables compared to the penalty method and decreases the number of the infeasible solutions. Fig. 1(b) applies the SVR method to the combinatorial optimization problem given in (a). The QUBO is expressed by The SVR method has one less spin than the penalty method, and thus H svr depends on a single spin x 1 . The reduced spin value x 2 is obtained by x 1 using the variable relationship in the linear equality constraint (i.e., x 2 = 1 − x 1 ). The SVR method induces a transition between the FSs (x 1 ) = (0) and (1) without passing large-energy infeasible solutions. In this way, reducing the number of the infeasible solutions leads to an efficient search for the optimum.
The main contributions of this paper are as follows: r We propose a new constraint-handling method called the SVR method for QUBO problems. The SVR method reduces the number of spins compared to the penalty method.
r We give a sufficient condition to obtain the optimum of the combinatorial optimization problem in the SVR method. This condition covers a wide range of constraints, which typically appear in QUBO problems.
r We apply the SVR method to two NP-hard combinatorial optimization problems: GkPP and QAP. Then we experimentally solve QUBO problems using an SA-based Ising machine and a QA-based Ising machine. The SVR method outperforms the penalty method. The rest of this paper is organized as follows. Section II defines the constrained combinatorial optimization problem and reviews the penalty method. Section III proposes the SVR method. Section IV formulates the SVR method in GkPP and QAP. Section V details the experimental results using SA-based and QA-based Ising machines. Section VI discusses the experimental results. Section VII summarizes this paper. Appendices, which can be found on the Computer Society Digital Library at http://doi. ieeecomputersociety.org/10.1109/TC.2023.3239539, give supplements on the application of the SVR method to d-dimensional systems, the domain-wall method for 2-dimensional systems, the results for extended version of GkPP, and the ideal QA simulation.

A. Definitions
This paper considers the following combinatorial optimization problem: Here, Q({x i } i∈V ) and C k ({x i } i∈V ) = 0 denote the objective function and the k-th constraint, respectively. M is the set of the linear equality constraints.
Here, the equivalence of combinatorial optimization problems is defined as: Definition 1: Given the set of the optimal solutions of problem a and problem b as

B. QUBO and Penalty Method
QUBO [23] is a common input format for Ising machines and is defined on an undirected graph G qubo = (V qubo , E qubo ), where V qubo and E qubo are the vertex set and edge set, respectively. They are given as 1} is a binary variable called the spin and Q i,j , Q i,i , Q 0 ∈ R. The QUBO problem finds the ground state of H({x i } i∈V qubo ).
(3) For a sufficiently large λ, C k ({x i } i∈V ) is zero for k ∈ M . Thus, the following theorem is obtained.
is expressed in the QUBO and C k ({x i } i∈V ) is a linear equality constraint in (1). It should be noted that even if Q({x i } i∈V ) contains higher order terms (e.g., x 1 x 2 x 3 ), introducing auxiliary variables always reduces these to second order terms [24]. To the best of our knowledge, only the penalty method is known to

III. SPIN-VARIABLE REDUCTION METHOD
Here, a new constraint-handling method called the SVR method is formulated (Section III-A). The SVR method formulates Problem 1 as a QUBO problem. The QUBOs are explicitly given for combinatorial optimization problems with typical constraints (Section III-B). Table I lists the notations used in this section.

A. Spin-Variable Reduction Method
Our proposed method is based on the idea of variable reduction. First, we introduce subsets V 1 , V 2 ⊆ V and M 1 , M 2 ⊆ M , where V and M were defined in Problem 1 (i.e., V is a set of the spin variables and M is a set of the linear equality constraints).
Independent Spins: i ∈ V 1 when a spin variable is used to represent other spin variables in terms of the variable relationships in linear equality constraints. The spin variables in V 1 are called independent spins.
Dependent spins i ∈ V 2 when a spin variable is expressed and calculated by independent spins. V 2 is given as V 2 = V \ V 1 . The spin variables in V 2 are called dependent spins.
Eliminated equality constraints k ∈ M 2 when the linear equality constraints C k ({x i } i∈V ) are eliminated along with the variable reduction. This is called the eliminated equality constraint.
The complement set of M 2 is defined as M 1 := M \ M 2 . A one-to-one correspondence exists between V 2 and M 2 (i.e., V 2 M 2 ) since the one equality constraint is used to eliminate one spin.
Dependent spins can be expressed in terms of independent spins. For k ∈ M 2 , the equality constraint gives where {ã i,k } for i ∈ V 2 and k ∈ M 2 can be regarded as a matrix with dimension |V 2 |. When the constraints in k ∈ M 2 are linearly independent, the matrix is given by a regular square matrix.
Here, the inverse matrix is introduced with a matrix elementã −1 where δ i,j is the Kronecker's delta. Then, (4) gives for i ∈ V 2 where a j,i = − k∈M 2ã j,kã Next, a new objective function and constraints are redefined as functions of independent spins. They are obtained by replacing the dependent spins by (6) and for k ∈ M 1 It should be noted that C k ({x i } i∈V 1 ) = 0 maintains a linear equality constraint. Then, the combinatorial optimization problem can be reformulated as and the constraints of dependent spins that for i ∈ V 2 Equation (9) gives the constraint for independent spins since R i ({x j } j∈V 1 ) can take values other than 0 or 1.
Proof of Lemma 2 The equivalence of Problems 1 and 3 are proven in terms of Definition 1. We consider sets of FSs for Problems 1 and 3, which are denoted by V (1) and For Let the sets of optimal solutions of Problems 1 and 3 be V (1) * and V (3) * , respectively. Then, Thus, Problems 1 and 3 are equivalent. Problem 3 is reformulated. We use the fact that It gives the following Problem 4 and Lemma 3.
Lemma 3: Problem 3 is equivalent to Problem 4 under the conditions that Proof of Lemma 3 For a sufficiently large λ, 1: under the conditions in (16). This is because the conditions restrict the range of , which is nothing but Problem 3. Lemma 2 shows an equivalence of Problems 1 and 3 and Lemma 3 shows an equivalence of Problems 3 and 4 under the conditions in (16). Thus, we obtain the following theorem.
Theorem 4: Problems 1 and 4 are equivalent under the conditions in (16).
The linearity of R i ({x j } j∈V 1 ) (see (5)) leads to the following proposition.
Proposition 5: Theorem 4 and Proposition 5 give the corollary: Corollary 6: (Optimality) The optimal solutions of Problem 1 are obtained by solving the QUBO problem of H svr ({x i } i∈V 1 ) in Problem 4 when the conditions in (16) are satisfied and Since V 1 ⊆ V , the SVR method requires fewer spins to formulate the QUBO than the penalty method.
Algorithm 1 practically formulates Problem 4 from Then, iteratively update the subsets and the constraints. Each update finds k ∈ M 1 and i ∈ V 1 that satisfy the condition where a j ∈ Z for j ∈ V 1 ∪ {0} (Line 3). Then, move k from M 1 to M 2 and i from V 1 to V 2 (Line 4). Accordingly, . The algorithm ends when k ∈ M 1 and i ∈ V 1 satisfying the condition in (17) does not exist.  (15)). Note that R i ({x m } m∈V 1 ) satisfies the condition in (16). In Algorithm 1, the condition in (17) determines which spin variable can take as a dependent spin. The condition is easily checked by looking at the coefficients of a linear equality constraint. The choice of dependent spins is arbitrary as long as the conditions are satisfied and may affect the performances of Ising machines. However, as will be shown in experiments (see Table IX in Section V), the performances are almost independent of the choice of dependent spins.
Below, the SVR method is demonstrated in simple examples of Problem 1.
The optimal solution is given by ( To solve the constrained combinatorial optimization problem in the SVR method, let x 1 be a dependent spin. The dependent spin is expressed by independent spins x 2 and x 3 as This equation meets the condition in (16) since a 0,1 = 2, a 2,1 = −1, and a 3,1 = −2 are integers. Then, the SVR method gives the QUBO as H svr is independent of dependent spin x 1 . Problem 4 finds the ground state of H svr .
Next, using this example, we show that the SVR method is inapplicable when (16) is not satisfied. Let x 3 be a dependent spin instead of x 1 . Then, does not meet the condition in (16) since a 1,3 = −1/2 and a 2,3 = −1/2 are not integers. The QUBO in the SVR method reads  1, 0). However, the constraint x 3 = R 3 (x 1 , x 2 ) indicates that x 3 takes an inappropriate value of 1/2. Thus, (16) is generally necessary to obtain the optimal solution in the SVR method. The condition in (16) is hereafter referred to as the SVR condition. Example 2: For x 1 , x 2 ∈ {0, 1}, minimize x 1 subject to a linear equality constraint x 1 + x 2 = 1.
The optimal solution is given by (x 1 , x 2 ) = (0, 1). Let x 1 be a dependent spin. Then, R 1 (x 2 ) = 1 − x 2 meets the SVR condition. The SVR method gives the QUBO without the constraint term as H svr = R 1 (x 2 ) since Here, we use x 2 2 = x 2 . The constraint term disappears because the domain of R 1 (x 2 ) matches the domain of a spin variable. Problem 4 minimizes the value of H svr . The ground state of H svr is x 2 = 1, and then x 1 = R 1 (x 2 ) gives the optimum (x 1 , x 2 ) = (0, 1).

B. Applications
Here, the SVR method is applied to combinatorial optimization problems with typical constraints in one-dimensional and two-dimensional systems. Application to the higherdimensional systems is straightforward (see Appendix A, available in the online supplemental material). Note that Algorithm 1 explicitly gives a QUBO for the SVR method in other types of constrained combinatorial optimization problems of Problem 1.
1) One-Dimensional System: First, consider Problem 1 with Q({x i } i∈V ) in the QUBO form under a linear equality constraint given by where k ∈ N. The constraint is called the k-hot constraint. It appears in the G2PP, the job sequencing problem, and graphcoloring problem. If the r-th spin is set as a dependent spin (i.e., r ∈ V 2 ), then Here, R r ({x i } i∈V 1 ) satisfies the SVR condition. The SVR method and the penalty method respectively give the QUBO H (1) svr and H (1) penalty (see Problem 4 and Problem 2) as where τ (1) and λ (1) are the constraint coefficients. The SVR method has one less spin than the penalty method. It should be noted that a different QUBO formulation called the domainwall method gives a QUBO with |V | − 1 spins for k = 1 [21]. Regardless, our method is applicable to arbitrary k.
2) Two-Dimensional System: Next, consider Problem 1 with Q({x i,s } (i,s)∈L (r) ⊗L (c) =V ) in the QUBO form under constraints given by where k (r) ∈ N and k (c) ∈ N. Here, the number of 1-valued spins in the FSs are set to N 1 = k (r) |L (r) | = k (c) |L (c) |. The constraints appear in the QAP, GkPP (k ≥ 3), and the traveling salesman problem. The QUBO obtained by the SVR method depends on the selection of dependent spins. Using a simple example where the spins in the p-th row or in the q-th column are set as dependent spins (i.e., V 2 = {(p, s)} s∈L (c) ∪ {(i, q)} i∈L (r) ), the dependent spins are represented by the independent spins as and Equations (27) and (28) satisfy the SVR condition.
The SVR method and the penalty method give the QUBO H (2) svr and H (2) penalty (see Problems 4 and 2), respectively, as where τ (2) and λ (2) are the constraint coefficients. The SVR method gives the QUBO with (|L (r) | − 1)(|L (c) | − 1) spins. By contrast, the penalty method gives the QUBO with |L (r) ||L (c) | spins. The domain-wall method formulates the QUBO with (|L (r) | − 1)|L (c) | spins when k (r) = 1 [25] (see also Appendix B, available in the online supplemental material). The SVR method at |L (r) | = |L (c) | and k (r) = k (c) = 1 is equivalent to the inserted method in [26]. Hence, our proposed method is superior to other known methods in terms of the number of spins.

IV. QUBO FORMULATION BY SPIN-VARIABLE REDUCTION METHOD
This section formulates the QUBO of combinatorial optimization problems by the SVR method. We consider the G2PP as examples of the one-dimensional problems and the GkPP (k ≥ 3) and the QAP as examples of the two-dimensional problems.

A. G2PP
The G2PP is specified by sets of vertices and edges, V g2p and E g2p with even |V g2p |. It partitions the vertex set into two subsets of equal size such that the number of edges connecting the two subsets is minimized. It is an NP-hard problem [27].
Here, the G2PP is transformed into Problem 1. First, we place a spin x i ∈ {0, 1} for each vertex i ∈ V g2p and set x i = 0 when vertex i is in the first subset and x i = 1 when vertex i is in the second subset. The objective function (i.e., the number of edges connecting the two subsets) is given as [12], [28] where k i denotes the degree of vertex i. The constraint is expressed as The constraint is satisfied when half of the spins are 0 and the other half are 1. Equations (31) and (23) with k = |V g2p |/2 have the same form. Hence, (25) gives the QUBOs for the SVR method and the penalty method. Next, we demonstrate the SVR method for an undirected graph with four spins (see Fig. 2). The optimal solutions are given by (x 1 , x 2 , x 3 , x 4 ) = (0, 0, 1, 1) or (1,1,0,0). We set x 1 as a dependent spin. Namely, Then the objective function for the SVR method is obtained as The objective function is independent of dependent spin x 1 . The constraint is described by Overall, the SVR method gives a QUBO as

B. GkPP
The GkPP (k ≥ 3) is an extension of the G2PP. It partitions the vertex set of graph G gkp = (V gkp , E gkp ) into k subsets of equal size such that the number of edges connecting the different subsets is minimized. Here, we assume |V gkp |/k to an integer.
The GkPP is transformed into Problem 1. First, we place spin x i,s ∈ {0, 1} for i ∈ V gkp and s ∈ {1, . . . , k} =: K. The spins denote the subset of each vertex. When x i,s = 1, vertex i belongs to the s-th subset. Then the objective function is given as The constraints are expressed as The constraints are satisfied if each vertex belongs to a subset and each subset size is equal to |V gkp |/k. Equation (37) satisfies the form in (26) with k (r) = 1 and k (c) = |V gkp |/k. Thus, (29) gives the QUBOs for the SVR method and the penalty method. See Appendix C, available in the online supplemental material, for an extended version of GkPP with inequality constraints.

C. QAP
The QAP assigns facilities to different locations to minimize the total transport cost [29]. The total transport cost is defined as (38) where |L qap | =: n fac , π denotes a permutation, and L qap ⊗ L qap =: V qap represents the set of facilities and locations. Here, f i,j indicates the flow amount between facilities i and j, while d s,t denotes the distance between locations s and t. The QAP is an NP-hard problem [30].
Here, the QAP is transformed into Problem 1. First, we place spin x i,s ∈ {0, 1} for each (i, s) ∈ V qap . The spins denote the location of each facility. When x i,s = 1, facility i is placed  i,j for i, j ∈ {1, 2, 3} and (b) d s,t for s, t ∈ {a, b, c}. (c) Optimal assignment for the location of the facilities. at location s. Then the objective function (i.e., S(π)) is given as [12], [14] The constraints are expressed as If all the constraints are satisfied, the permutation is described as x i,s = 1 when s = π(i). Otherwise x i,s = 0. Equation (40) adopts the form in (26) with k (r) = k (c) = 1. Consequently, (29) gives the QUBOs for the SVR method and the penalty method.
As an example, we demonstrate the SVR method in a QAP with three facilities and locations. We set (f 1,2 , f 2,3 , f 1,3 ) = (1, 0, 0) and (d a,b , d b,c , d a,c ) = (0, 1, 1) (see Fig. 3). We assumed that f i,j and d s,t are symmetric (i.e., f i,j = f j,i and d s,t = d t,s ) and the diagonals are zero (i.e., f i,i = 0 and d s,s = 0). There are two optimal solutions. One allocates facility 1 to location a, facility 2 to location b, and facility 3 to location c and the other allocates facility 1 to location b, facility 2 to location a, and facility 3 to location c. The total transport cost in the optimal solutions is 0. Here, we set the dependent spins as V 2 = {(1, a), (1, b), (1, c), (2, a), (3, a)}. Namely, The objective function is obtained as The constraint term is described by gives the lowest value of H qap = 0. The lowest-value solutions give the optimal solutions by calculating the dependent-spin values using

V. EXPERIMENTAL EVALUATION USING ISING MACHINES
This section compares the performances of the SVR method and the penalty method for GkPPs and QAPs when an SA-based Ising machine or a QA-based Ising machine is used. Additionally, we show the performance of the domain-wall method for GkPPs (k ≥ 3) and QAPs (see the QUBO formulation in Appendix B, available in the online supplemental material).

A. Set Up of Ising Machines
SA [31], [32], [33] and QA [34], [35] are meta-heuristic algorithms that address combinatorial optimization problems. Ising machines based on SA or QA are designed to solve QUBO problems. Hereafter, the SA-based Ising machine [7] and the QA-based Ising machine [36] are referred to as classical IM and quantum IM, respectively. The classical and quantum IM hardware embed a maximum of 130,000 spins 2 and 180 spins on a complete graph, respectively [7], [37]. The IM experiments run on a MacBook Pro with a 2.8GHz quad core Intel Core i7 (16GB RAM) using Python 3.7.6 as the implementation language. Table V lists the parameter sets for classical IM and quantum IM. The annealing time was set to 1 second for the classical IM and 20μ second for the quantum IM. The number of runs was set to 10 for the classical IM and 1,000 for the quantum IM. The remaining values were set to the default for each IM.

B. Methods
We adopted GkPPs and QAPs to compare the performances of the SVR method, the penalty method, and the domain-wall method. Note that the domain-wall method is inapplicable to the G2PP. The G2PP instances are specified by the undirected graph G g2p = (V g2p , E g2p ). For each problem size |V g2p |, 20 undirected graphs with edge density of 0.5 were generated. Connecting the vertices i and j in V g2p with half probability created the undirected graphs. The QUBO for a G2PP has the constraint coefficient τ (1) in the SVR method and λ (1) in the penalty method ( (25)). The argument in [12] gives a sufficient condition to satisfy the constraint in the ground state of the QUBO as The same argument holds for the SVR method. Thus, τ (1) and λ (1) are set as the value of the right-hand side of (45) in both the classical and quantum IM experiments. GkPP (k ≥ 3) generated 5 undirected graphs with edge density of 0.5 for each graph size and QAP used instances in Refs. [38], [39]. The QAP instances are respectively called nug-n fac and lipa-n fac -a. The constraint coefficient in the QUBO for a GkPP and a QAP is τ (2) in the SVR method, λ (2) in the penalty method (29), and κ in the domain-wall method (Appendix B, available in the online supplemental material). The constraint coefficients were determined using the procedure in [40] in classical IM experiments. We calculated the FS rate and the cost function. The FS rate represents the ratio of the number of obtained FSs to the number of runs. The cost function was calculated for the FSs. A smaller value of the cost function indicates a better performance. The FS rate and cost function typically increase with the constraint coefficients [15]. The dependences indicate that the constraint coefficient has an optimal value. To systematically determine the optimal value of the constraint coefficient, we set the threshold value of the FS rate as r th = 0.8 and repeatedly solved each problem instance while varying the constraint coefficient. The precision threshold was set to 10. The quantum IM experiments used the optimal value of the constraint coefficient obtained in classical IM experiments.
The classical IM hardware has physical spins on a complete graph architecture [7]. Hence, the number of physical spins on the architecture is identical to the number of logical spins in a QUBO (i.e., |V qubo |). On the other hand, the quantum IM hardware has physical spins on a sparse graph architecture [37]. Therefore, minor embedding is generally required to map G qubo onto the architecture [41], [42]. The minor-embedding process constructs chains of physical spins to represent logical spin states. The physical spins in a chain interact via ferromagnetic coupling. The quantum IM experiments calculated the chain break fraction, which is the ratio of chains whose spins have different values. When the chain-break fraction averaged over the 1,000 runs is greater than 0.1, the ferromagnetic coupling strength in a chain was increased. Here, the number of logical spins and physical spins are denoted by N L and N P , respectively.
Below, the average cost function, the minimum cost function, the FS rate, and N P are used. The average cost function is the value averaged over the FSs, while the minimum cost function is the lowest value among the FSs. In GkPPs, each value was averaged over problem instances. The average cost function and the minimal cost function are denoted by C ave and C min , respectively. Table VI gives the results of G2PP in a classical IM for different |V g2p | using the SVR method and the penalty method. The SVR method selects x 1 as a dependent spin. Both methods find C ave = C min for small sized problems, indicating that the classical IM obtains the optimum for all 20 problem instances. By contrast, the SVR method outperforms the penalty method for large-sized problems. The parenthesis denotes the relative value of C ave in the SVR method when C ave in the penalty method is set to 1. The decrease in the relative value with |V g2p | indicates the superior performance of the SVR method for larger-sized problems. The FS rate is always 1 in both methods. Table VII shows the results of GkPP (k ≥ 3) for various k. The SVR method selects the first-raw spins {x 1,s } s∈K and the first-column spins {x i,1 } i∈V gkp as dependent spins. The smallsized problem (i.e., |V gkp | = 12) finds C ave = C min for all the methods. On the other hand, for a large-sized problem (i.e., |V gkp | = 360), the penalty method has worst performance among them. The SVR method is best for large k, and the domain-wall method is best for small k. The constraint coefficients in the penalty method are larger than the others. Table VIII shows the results of QAP in a classical IM with different n fac . The SVR method selects the first-raw spins {x 1,s } s∈L qap and the first-column spins {x i,1 } i∈L qap as dependent spins. The SVR method finds the optimum when n fac ≤ 25, whereas the penalty method finds the optimum when n fac ≤ 20 and the domain-wall method finds the optimum when n fac ≤ 15. For larger n fac , the SVR method outperforms the other methods. The domain-wall method fails to yield a FS rate above the threshold value in lipa80a even at κ = 10 4 . The constraint coefficients in the SVR method τ (2) are smaller than those in the penalty method λ (2) , especially for lipa-n fac -a problem instances. Table IX shows the G2PP results of |V g2p | = 1024, the G8PP results of |V g8p | = 360, and the QAP results of lipa50a for 10 different choices of dependent spins. The G2PP randomly selects dependent spin r from the vertex set V g2p and the G8PP and the QAP randomly sets p-th row spins and q-th column spins to dependent spins. In all cases, the SVR method produces smaller average cost functions and minimum cost functions than those of the penalty method and the domain-wall method, indicating that the results shown in Tables VI, VII, and VIII are independent of the choice of dependent spins. In G8PP and QAP, the optimal values of τ (2) are almost independent of the choice of dependent spins. Table X gives the results of G2PP in a quantum IM for different |V g2p | using the SVR method and the penalty method. Both methods provide the same value of C min as that in a classical IM for |V g2p | = 8. For larger |V g2p |, the quantum IM gives larger values of C ave and C min than the classical IM. C ave and C min are almost the same in both methods. By contrast, the FS rate in the SVR method is higher than that in the penalty method. The SVR method has smaller N P than the penalty method.   Table XI shows the results of GkPP (k ≥ 3) for k ∈ {3, 4, 6} and |V gkp | = 12. When k = 6, any method does not find a FS of all problem instances. When k is 3 or 4, all methods produce approximately the same values of C ave and C min , but the SVR method gives the highest FS rate. The SVR method has the lowest N P of the three methods for all k. Table XII shows the results of QAP in a quantum IM for different n fac . When n fac is 10 or 12, any method fails to find an FS with the quantum IM. For a smaller number of facilities, all methods generate almost the same values of C ave and C min , but the SVR method gives the highest FS rate. The SVR method has the smallest N P among the three methods.

E. Preprocessing-Time Comparison
Table XIII shows the preprocessing times for generating QUBOs of GkPP and QAP using the SVR method and the penalty method. We used a software package from [7]. For all problems, the SVR method takes longer to preprocess than the penalty method. For G2PPs, the increase rate is approximately 5% regardless of problem size |V g2p |, and preprocessing time scales as |V g2p | 2 . The increase rate of GkPPs (k ≥ 4) increases with k. The processing time scales as k 2 for the SVR method and k for the penalty method. These are explained by QUBO edge number scalings (i.e., |E qubo | ∼ k 2 for H svr and |E qubo | ∼ k for H penalty ). For QAPs, the increase rate is approximately 30% regardless of problem size n fac , and preprocessing time scales as n 4 fac .

VI. DISCUSSION
In a classical IM, the SVR method generates smaller cost functions in both GkPP and QAP than the penalty method. We explain the results in terms of the energy landscape structure. We define the set of FS subspaces {S α } such that S α is a set of FSs and two FSs are in the same S α if the single-spin flips can make a transition between the two FSs without passing an infeasible solution. A single-spin flip denotes the change of a spin value from 0 to 1 or vice versa. For a G2PP, the penalty method needs at least two-spin flips to induce a transition between the FSs. Each FS belongs to the different S α , and thus |S α | exponentially increases with the number of spins |V g2p |  TABLE XII  RESULTS OF THE SVR METHOD, THE PENALTY METHOD, AND THE DOMAIN-WALL METHOD IN QAPS USING A QUANTUM IM. OPT. DENOTES THE OPTIMUM   TABLE XIII  PREPROCESSING TIME FOR QUBO GENERATION IN THE SVR METHOD AND  THE PENALTY METHOD as |S α | = |V g2p | |V g2p |/2 ∼ |V g2p | −1/2 2 |V g2p | . By contrast, the SVR method gives |S α | = 1. For example, in the four-spin G2PP (Table III), the transition from FS (x 2 , x 3 , x 4 ) = (0, 0, 1) to FS (0, 1, 1) is induced by a single-spin flip of x 3 . For a GkPP (k ≥ 3) and QAP, the penalty method needs at least four-spin flips to induce a transition between the FSs. The SVR method reduces |S α | compared to the penalty method. For example, in a QAP with n fac = 3 (Table IV), (x 2,b , x 3,b , x 2,c , x 3,c ) = (0, 0, 0, 1), (1, 0, 0, 0), and (1, 0, 0, 1) are within the same FS subspace. The reduction of |S α | efficiently searches lower-energy solutions. In a quantum IM, the SVR method gives a higher FS rate in both GkPP and QAP than the other methods. This is attributed to the smaller number of infeasible solutions.
Due to the minor embedding, the SVR method and the penalty method in the quantum IM have the same number of FS subspaces because a multi-spin flip of physical spins is necessary to flip a single logical spin. Therefore, the results for an ideal QA simulation without minor embedding (see the simulation details in Appendix D, available in the online supplemental material) should be compared with the quantum IM results. Fig. 4 shows the annealing-time dependences of the average cost function and the FS rate in the G2PPs with |V g2p | = 16. For the ideal QA, the average cost function approaches the optimum and the FS rate approaches 1 with the annealing time. For a long annealing time, the average cost function is smaller in the SVR method than the penalty method, implying that reducing the number of FS subspaces improves the performance of the ideal QA. Moreover, the SVR method has a higher FS rate than the penalty method at a short annealing time, but at a long annealing time they have the same FS rate. By contrast, the quantum IM results weakly depend on the annealing time. The average cost function is smaller and the FS rate is higher as the annealing time increases. Note that ideal QA finds better solutions in less time than the quantum IM. We expect this to be due to the minor embedding because the minor embedding reduces quantum fluctuation effects [17]. An interesting future work is to theoretically and experimentally investigate the effect of minor embedding on the performance.

VII. CONCLUSION
We propose the SVR method to handle the linear equality constraints for combinatorial optimization problems in the QUBO form. The SVR method reduces the number of spin variables in the QUBO problem compared to the penalty method and the other known methods. Here, we determine the sufficient condition to obtain the optimum of the combinatorial optimization problem in the SVR method and provide QUBO formulations for combinatorial optimization problems with typical constraints in one-dimensional and two-dimensional systems. Therefore, the SVR method extends the application of Ising machines to larger-size combinatorial optimization problems with linear equality constraints. Additionally, we compare the performances of the SVR method, the penalty method, and the domain-wall method in GkPP and QAP using an SA-based Ising machine and a QA-based Ising machine. The SVR method outperforms the other methods in terms of the cost function in an SA-based Ising machine and the FS rate in a QA-based Ising machine.

ACKNOWLEDGMENTS
This work was supported in part by JST CREST under Grant JPMJCR19K4, Japan and in part by KAKENHI under Grant 21K03391. The authors thank Shu Tanaka for his fruitful comments. T. S. thanks the Supercomputer Center, the Institute for Solid State Physics, The University of Tokyo for use of the facilities. Tatsuhiko Shirai (Member, IEEE) received the BSci, MSci, and DrSci degrees from the University of Tokyo, in 2011, 2013, and 2016, respectively. He is presently an assistant professor with the Department of Computer Science and Communications Engineering, Waseda University. His research interests include quantum dynamics, statistical mechanics, and computational science. He is a member of the JPS.
Nozomu Togawa (Member, IEEE) received the BEng, MEng, and DrEng degrees from Waseda University, in 1992, 1994, and 1997, respectively, all in electrical engineering. He is presently a professor with the Department of Computer Science and Communications Engineering, Waseda University. His research interests include VLSI design, graph theory, and computational geometry. He is a member of the IEICE and IPSJ.