Post-processing variationally scheduled quantum algorithm for constrained combinatorial optimization problems

We propose a post-processing variationally scheduled quantum algorithm (pVSQA) for solving constrained combinatorial optimization problems (COPs). COPs are typically transformed into ground-state search problems of the Ising model on a quantum annealer or gate-type quantum device. Variational methods are used to find an optimal schedule function that leads to high-quality solutions in a short amount of time. Post-processing techniques convert the output solutions of the quantum devices to satisfy the constraints of the COPs. pVSQA combines the variational methods and the post-processing technique. We obtain a sufficient condition for constrained COPs to apply pVSQA based on a greedy post-processing algorithm. We apply the proposed method to two constrained NP-hard COPs: the graph partitioning problem and the quadratic knapsack problem. pVSQA on a simulator shows that a small number of variational parameters is sufficient to achieve a (near-)optimal performance within a predetermined operation time. Then building upon the simulator results, we implement pVSQA on a quantum annealer and a gate-type quantum device. The experimental results demonstrate the effectiveness of our proposed method.


I. INTRODUCTION
Q UANTUM algorithms for solving combinatorial opti- mization problems (COPs) have attracted attention in both academia and industry.COPs are typically transformed into the ground-state search problems of the Ising models on quantum devices.Adiabatic quantum annealing is a representative quantum algorithm to search for the ground state by following an adiabatic path [1], [2].The computational time is approximately proportional to the square of the inverse of the spectral gap [3].Because the spectral gap generally decreases with problem size, adiabatic quantum annealing is difficult to execute for large-scale problems.
Recent advances in quantum technology have realized quantum devices such as quantum annealers [4], [5] and gatebased quantum devices [6].These devices provide diverse platforms for executing quantum computation, but they are susceptible to various forms of noise, including coherent and incoherent noise.Noise limits their applicability to quantum algorithms with low computational costs.For example, an experiment on a quantum annealer with 2, 000 spins re-vealed noise effects within a short timescale of approximately 100ns [7].Additionally, the computational accuracy of gatebased quantum devices decreases as the number of spins increases [6].
Given the severe limitations of the adiabatic quantum annealing and the near-term quantum devices, the development of quantum algorithms that obtain the low-energy states (i.e., high-quality solutions of COPs) within short operation time is crucial.Several ideas have been proposed toward this direction.One is to variationally optimize a schedule function to drive the system into the lower-energy state within a short amount of time (see Fig. 1 for an example of the schedule function s (m) (t) in this work).Examples include variationally scheduled quantum simulation for quantum annealers [8] and quantum approximate optimization algorithm (QAOA) for gate-based quantum devices [9].These methods can find a schedule function that leads to high-quality solutions within short operation time [10]- [13].
Another idea is to employ post-processing techniques.When constrained COPs are considered, post-processing en-Quantum computation on quantum annealer or gate-type quantum device s (m) (t) : schedule function FIGURE 1: Illustration of pVSQA for a two-spin COP with a constraint σ 1 + σ 2 = 0. Iterating the quantum computation and classical optimization m times gives a schedule function s (m) (t).p ′ (σ) is the probability distribution of the spin configurations obtained by quantum measurements.When a post-processing is performed, p(σ) becomes populated only with feasible solutions.The classical optimizer calculates the energy expectation value of the cost function C(σ) using p(σ) and updates the schedule function s (m) (t) (dashed line) to s (m+1) (t) (solid line).
hances the solution quality.The spin configurations for a constrained COP are divided into two solutions: feasible and infeasible solutions [14], [15].Feasible solutions are those that satisfy the constraints.Post-processing refers to a transformation of infeasible solutions into feasible ones.Postprocessing techniques have been applied to slot-placement problems and quadratic assignment problems [16], [17].This study combines a variational method and a postprocessing technique to efficiently obtain high-quality solutions of constrained COPs.Fig. 1 provides a flowchart of our proposed algorithm, which we call the post-processing variationally scheduled quantum algorithm (pVSQA).As an illustration, consider a COP with two spins σ = (σ 1 , σ 2 ) ∈ {−1, 1} ⊗2 subject to a constraint σ 1 + σ 2 = 0.The spin configurations (σ 1 , σ 2 ) = (1, −1) and (−1, 1) are feasible solutions, and all others are infeasible solutions.First, the quantum device generates a variational quantum state via quantum computation.Then quantum measurements of the quantum state provide the probability distribution of the spin configuration denoted by p ′ (σ).The probability distribution generally has a population over the infeasible solutions (i.e., (σ 1 , σ 2 ) = (−1, −1) and (1, 1)).Post-processing transforms the infeasible solutions into feasible ones.Hence, the probability distribution p(σ) is populated only with feasible solutions.The classical optimizer calculates the energy expectation value of the cost function using p(σ) and updates the schedule function.Repeating this process finds an optimized schedule function that gives low-energy states.
The three main contributions of this paper are given as follows: • We propose pVSQA to efficiently provide high-quality solutions of constrained COPs.We obtain a sufficient condition for constrained COPs to apply pVSQA based on a greedy post-processing algorithm.We apply it to the graph partitioning problem (GPP) and the quadratic knapsack problem (QKP).• We execute pVSQA on a simulator.pVSQA requires a small number of variational parameters to achieve a (near-)optimal performance within a predetermined operation time.
• The experiments on a quantum annealer and a gatebased quantum device demonstrates that pVSQA outperforms conventional quantum annealing with and without a post-processing and QAOA without a postprocessing.The rest of this paper is organized as follows.Section II provides the settings and preliminaries.Section III proposes pVSQA.Section IV exhibits the simulator results.Section V shows the experimental results on a quantum annealer and a gate-based quantum device.Section VI discusses the results.Section VII summarizes this paper.

II. PRELIMINARIES
Ising models generically map various computationally hard COPs [14], [18], [19].Ising models are defined on undirected graph G = (V, E) as where σx,y,z i are the Pauli spins acting on site i, J ij ∈ R denotes an interaction between spin i and spin j, h i ∈ R denotes a magnetic field on spin i, and H 0 ∈ R. The Ising model for a constrained COP is given as where Ĥobj and Ĥcst are Ising models for the objective function and the constraints of the COP, respectively.The cost function C(σ) is given as Here Quadratic unconstrained binary optimization (QUBO) is equivalent to the Ising model [20] and is defined on an undirected graph G q = (V q , E q ) as where x i ∈ {0, 1}, q ij ∈ R, and Q 0 ∈ R. Replacing x i in the QUBO by (σ z i + 1)/2 gives the Ising model.This study considers the quantum computation of the following Hamiltonian dynamics on a simulator or quantum devices, based on the quantum annealing algorithm and QAOA.The quantum state |ψ(t)⟩ obeys the Schrödinger equation: where the Planck constant ℏ = 1.Ĥq denotes the source of quantum fluctuations, and herein the transverse-field Hamiltonian Ĥq = − i∈V σx i is adopted.The initial state |ψ(0)⟩ is set to the ground state of Ĥq .The schedule function is characterized by s(t) ∈ [0, 1] for t ∈ [0, T ], where the operation time is denoted by T .s(t) is variationally determined in pVSQA, whereas it is fixed in conventional quantum annealing.

III. PVSQA
This section details pVSQA.First, pVSQA is formulated.Then we give three types of schedule functions since the function form depends on the type of quantum devices.Finally, we give a condition for the applicability of pVSQA to constrained COPs.

A. FORMULATION OF PVSQA
Fig. 1 illustrates the pVSQA.The algorithm first gives the Ising model corresponding to a COP and the initial schedule function s (0) (t) for t ∈ [0, T ] as the input.Then iterating the following process m times updates the schedule function to s (m) (t).
The quantum device performs the quantum computation along the schedule function s (m) (t) to generate |ψ(T )⟩.The probability distribution of the spin configuration is given as When using a real quantum device, the probability distribution is approximately obtained by repeating the quantum projection measurements of |ψ(T )⟩ on the computational basis.
Then a post-processing is performed.The post-processing is defined by mapping P : P σ → σ, which transforms all infeasible solutions into feasible ones (see Sec. III-C for the explicit mapping of P ).Furthermore, the probability distribution can be renormalized by using the top r% of solutions in terms of the cost function [21].The probability distribution is recalculated as where θ(•) is the Heaviside step function and C r is chosen to satisfy the normalization of p r (σ) (i.e., σ p r (σ) = 1).
Here, p r (σ) is populated only with feasible solutions.The energy expectation value of the cost function is given as A classical optimizer modifies the schedule function s (m) (t) to lower the energy expectation value and updates s (m) (t) to s (m+1) (t).It relies on an optimization algorithm such as the gradient descent method or the Powell method to construct a progressively optimized sequence of variational parameters [22].Finally, when the schedule function is optimized, the performance of pVSQA is evaluated using probability distribution p r (σ).
As a comparison to pVSQA, a variational quantum algorithm without a post-processing called VSQA is introduced.The quantum device performs the quantum computation in (5).The probability distribution can be renormalized as where 0 < r ≤ 100, C ′ r is a threshold to satisfy σ p ′ r (σ) = 1, and the cost function C ′ (σ) is given as A classical computer calculates the energy E ′ as Since the constraint Hamiltonian H cst of the cost function increases the energy of infeasible solutions, feasible solutions are easily obtained.Then a classical optimizer lowers E ′ to update the schedule function s (m) (t).
Here, it is noted that pVSQA converges to the ground state of ĤIsing as T tends to infinity.The reasoning is as follows.Let us recall that the adiabatic theorem ensures the convergence of quantum annealing [3].Since the schedule function of quantum annealing is a specific case of pVSQA, the optimal value of the cost function for pVSQA is no larger than that of quantum annealing.Given that the cost function takes its minimum value at the ground state, it follows that pVSQA attains this state as well.

B. TYPE OF SCHEDULE FUNCTIONS
We consider three types of schedule functions: continuous, linear, and QAOA schedules.The continuous schedule represents s(t) as a continuous function, the linear schedule represents s(t) as a linear function, and the QAOA schedule has s(t) that takes values of either zero or one.Although the continuous schedule shows the best performance among the three, its implementation on quantum devices remains challenging.On the other hand, the linear schedule can be performed on quantum annealers, whereas the QAOA schedule is easily implemented on gate-based quantum devices.First, we discuss the optimal schedule function s * (t) in the continuous schedule.According to optimal control theory, for a short operation time T , s * (0) = 1 and s * (T ) = 0, and the schedule function between the two can be a smooth path [23].The results can be extended straightforwardly to pVSQA with a post-processing.The linear schedule has two variational parameters s 1 and s 2 , which represent the initial and final values of s(t), respectively (see Fig. 3 (a)).The schedule function is given as It should be noted that previous studies fixed s(0) = 0 and s(T ) = 1, and then used the values of s(t) during quantum simulation as variational parameters [8], [12], [13].Unlike previous studies, we adopt s 1 and s 2 as variational parameters since the optimal schedule function in the continuous schedule does not always start with s(0) = 0 or end with s(T ) = 1 (see Fig. 2).The linear schedule is approximately implemented on quantum annealers (see Sec. V-B1 for details).
The QAOA schedule has p layers.Each layer consists of two variational parameters (see Fig. 3 (b)).Thus, there are 2p variational parameters denoted by {s i } 2p i=1 .The time duration of s(t) = 1 and s(t) = 0 for layer ℓ ∈ {1, . . ., p} is respectively determined by where s 0 = 0. We impose the following conditions on the variational parameters: s i−1 ≤ s i for i ∈ {1, . . ., 2p} and s 2p ≤ T .The quantum state is not time-evolved between t ∈ [s 2p , T ].These conditions allow the performance of the QAOA schedule to be compared with the other types of schedules.The QAOA schedule is easily implemented on real gate-based quantum devices (see Sec. V-B2 for details).

C. APPLICABILITY OF PVSQA
We consider following COPs.Problem 1.For binary variables x i ∈ {0, 1} for i ∈ V , find a configuration that minimizes the objective function Q obj under the following constraints where Here, b m min , b m max , and a m i for i ∈ V m are integers and M is the number of constraints.
Table 1 classifies the constraints into two types: independent constraint and dependent constraint.Constraints are called independent when V m ∩ V n = ∅ for all pairs of m and n (m ̸ = n).Otherwise, they are called dependent.
pVSQA is applicable when post-processing transforms all infeasible solutions into feasible ones (i.e., mapping P (•) in (7) exists).Here, we provide a sufficient condition to have a greedy post-processing algorithm.The conditions hold for COPs with independent constraints and cover a wide range of COPs including GPPs, knapsack problems, and graphcoloring problems [14].On the other hand, a post-processing algorithm should be considered for each COP with dependent constraints.For example, an adhoc post-processing was found for quadratic assignment problems [16].It is herein noted that a post-processing algorithm with a low computational cost is not found for all types of the constraints.This situation occurs when finding a feasible solution is NPcomplete.Extending the applicability of pVSQA remains a challenge.
Below, we first prove a theorem for COPs with independent constraints.Then we provide a greedy post-processing algorithm and demonstrate its applicability in two constrained COPs: GPP and QKP.

1) Theorem
First, in Problem 1, we consider an objective function Q obj with a single constraint (i.e., M = 1): where b 1 min , b 1 max , and a 1 i for i ∈ V 1 ⊆ V are integers.Then we prove a following lemma.Lemma 1.For Problem 1 with a constraint in (15), consider a Hamiltonian Q ′ as where A ′ denotes a constraint coefficient and Then all local minimum solutions of Q ′ are feasible solutions under the following conditions: 1) A ′ > ∆Q 1 obj where ∆Q 1 obj is the maximum energy change of Q obj when flipping a binary variable in Here, a flip indicates a change in the value of a binary variable from 0 to 1 or vice versa.Then the local minimum solution is defined as a solution that does not lower the energy of Q ′ by flipping any binary variable.
Proof of Lemma 1.To prove the lemma, we show that all infeasible solutions are not local minimum solutions under the three conditions.Suppose that {x i } i∈V is an infeasible solution.Then is satisfied.When ( 18) holds, we find from the second condition that for any This indicates that the configurations obtained by flipping a single binary variable do not satisfy (19).
Next, we show that the infeasible solution is not a local minimum solution using a proof by contradiction.Assume that the infeasible solution is a local minimum solution.Equation (20) implies that if a flip increases the value of i∈V 1 a 1 i x i , it lowers the value of Q ′ when A ′ > ∆Q 1 obj .Therefore, a flip should not increase the value of i∈V 1 a 1 i x i , implying that a 1 i ≥ 0 when x i = 1 and a 1 i ≤ 0 when x i = 0. Thus, for any configuration {x ′ i }.This contradicts the third condition that feasible solutions exist and indicates that {x i } is not a local minimum solution.For infeasible solutions that satisfy (19), similar arguments show that they are not local minimum solutions.Consequently, when the three conditions are satisfied, all local minimum solutions are feasible solutions.
Then Lemma 1 gives a following theorem.
Theorem 2. For Problem 1 with constraints in (14), consider a Hamiltonian Q ′ as where A ′ denotes a constraint coefficient and Then all local minimum solutions of Q ′ are feasible solutions under the following conditions: 1) A ′ > ∆Q obj where ∆Q obj is the maximum energy change of Q obj when flipping a binary variable in 3) Constraints are independent (i.e., V m ∩ V n = ∅ for all m and n (m ̸ = n)).4) Feasible solutions exist.
Proof of Theorem 2. We first show that all local minimum solutions of Q ′ satisfy the mth constraint in (14) under the four conditions of Theorem 2. Let us rewrite (22) in the form of (16) as where Thus, the first condition of Lemma 1 is satisfied.The remaining conditions of Lemma 1 trivially hold from the conditions of Theorem 2. Lemma 1 implies that all local minimum solutions of Q ′ satisfy the mth constraint in (14).
The above argument holds for any m ∈ {1, . . ., M }.Thus, all local minimum solutions of Q ′ are feasible solutions under the four conditions of Theorem 2.
We give some remarks on Theorem 2. First, it includes the equality constraint called k-hot constraint.An inequality constraint in (14) with b m min = b m max = k, and a m i = 1 for i ∈ V m is reduced to the k-hot constraint: The Second, we comment on the conditions of Theorem 2. As for the first condition, when Q obj is formulated as a QUBO (i.e., Q obj = (i,j)∈Eq q ij x i x j + i∈Vq q i x i ), an upper bound of ∆Q obj is given as A ′ can be set to a value larger than the right-hand side of the equation.The second and third conditions are easily checked by looking into {a m i } i∈V m , b m min , and b m max in (14).The fourth condition is checked by performing a local optimization of Q ′ (e.g.Algorithm 1) for a spin configuration.Under Algorithm 1 Greedy post-processing algorithm: P (σ ′ ) = σ in ( 7) for i ∈ V do if ∆Q j < 0 then 9: x j ← 1 − x j 10: end if 13: end while 14: σ i ← 2x i − 1 for i ∈ V the first, second, and third conditions, the spin configuration always reaches a feasible solution by repeating a local optimization when the fourth condition is satisfied.Otherwise, the feasible solutions do not exist.Post-processing is given by a local optimization method using the property that all local minimum solutions are feasible solutions.The local optimization method flips a binary variable to lower energy Q ′ .Since infeasible solutions are not local minimum solutions, repeating the local optimization method can give a feasible solution.
In this study, we adopt a greedy method as a local optimization method.Algorithm 1 gives a post-processing algorithm based on a greedy method.The algorithm provides a mapping P in pVSQA (see (7)).The greedy method flips a binary variable resulting in the lowest energy of Q ′ when it is flipped (line 4-12).We perform the greedy post-processing algorithm for all spin configurations obtained via the quantum computation.Namely, if the obtained feasible solution is not a local minimum solution, post-processing transforms it into a feasible solution with a lower energy.
Fig. 4 shows the problem-size dependence of the number of flips to convert a spin configuration into a local minimum solution in the greedy post-processing algorithm for GPPs with a |V |/2-hot constraint.The data represent the mean and the standard deviation over 1000 random spin configurations for each of 10 different GPPs.The number of flips increases with the problem size, but the rate of increase slows down with the problem size.

3) Examples of constrained COPs: GPP and QKP
We consider two types of constrained COPs: GPP and QKP.Both problems are categorized as NP-hard [24], [25].
GPP divides the node set of a given graph G g = (V g , E g ) into two subsets with equal size, minimizing the number of cuts.Here, it is assumed that the graph has an even number of nodes.The cut denotes the edge connecting the two subsets.The binary variable x i ∈ {0, 1} is assigned for each node such that x i = 0 (x i = 1) indicates that node i belongs to the one (the other) subset.Then this problem has a |V g |/2-hot constraint, which is given as pVSQA is applicable to GPP since the constraint satisfies the conditions of Theorem 2. QUBO for the GPP is given as [26] Q where k i denotes the degree of node i.Here, Q obj and Q cst are QUBOs for the objective function and the constraints, respectively.A ≥ 0 is the constraint coefficient.QKP finds a set of items to yield the highest profit within the capacity of the knapsack.The set of n items with weights {w i } n i=1 and profits {p ij } n i≤j , and the capacity of the knapsack C are given.Here, p ii denotes the profit of item i and p ij denotes the additional profit when item i and item j are chosen.The binary variable x i ∈ {0, 1} is assigned for each item such that x i = 1 when item i is chosen and x i = 0 when item i is not chosen.This problem has an inequality constraint, which is given as pVSQA is applicable to QKP since the constraint satisfies the conditions of Theorem 2. QUBO for the QKP is given as where Q obj and Q cst are QUBOs for the objective function and the constraints, respectively.A ≥ 0 is the constraint coefficient.Here, we adopt a mapping in [27] for an inequality constraint.Unlike other mappings [28], it does not require additional auxiliary binary variables and is suitable for nearterm quantum devices with severe limitations on the input number of spins.

IV. SIMULATOR RESULTS
The performance of pVSQA is initially evaluated using a simulator.The simulator employs the fourth-order RungeKutta method to solve the Schrödinger equation in (5).We use Python 3.7.6 integrated with intel Fortran as the implementation language.

A. INSTANCES OF COPS
We generate 10 instances for GPPs and 10 instances for QKPs for each problem size.The GPP instances are specified by the undirected graph G g = (V g , E g ).Random graphs are created by connecting each pair of vertices with a probability of 0.5 to form edges.The QKP instances with n items are generated by benchmarking instances with the 100 items listed in [29].They are labeled 100_100_i where i ∈ {1, . . ., 10} denotes the instance.Each instance is specified by a profit matrix {p  [30].
For each problem instance, we calculate the Ising model using the QUBOs in (30) and (32).The coefficients of the Ising model are rescaled so that the maximal absolute value of {J ij } (i,j)∈E and {h i } i∈V is one.This rescaling determines the timescale for quantum computation.

B. METHOD
We compare the performances of pVSQA, VSQA, pQA, and QA.Here, pQA represents quantum annealing with a postprocessing, and QA represents quantum annealing without a post-processing.In pQA and QA, we adopt a linearly interpolated schedule function from s(0) = 0 to s(T ) = 1 (i.e., s(t) = t/T ).pVSQA and VSQA optimize the schedule function using a classical optimizer.For continuous schedules, we follow the method in [23] based on optimal control theory and use a gradient descent method as a classical optimizer with an initial condition of s(t) = 0.5 for t ∈ [0, T ].We set r = 100 in ( 7) and ( 9).For linear and QAOA schedules, we adopt the Powell optimizer implemented in SciPy as a classical optimizer [31].For the linear schedule, the maximum number of iterations is set to 10 with initial values of s 1 = s 2 = 0.5.For the QAOA schedule, the maximum number of iterations is set to 10p and initially s i = 0 for i ∈ {1, . . ., 2p} [32].The default values are used for the remaining input parameters.
The constraint coefficient A is determined using the procedure in [28].Hereafter, p 100 (σ) and p ′ 100 (σ) respectively denote the probability distributions in ( 7) and ( 9) obtained for either the optimized schedule function in pVSQA and VSQA or the given schedule function in pQA and QA.We calculate the feasible solution rate and the average cost function.The feasible solution rate is 1 for pVSQA and pQA due to the post-processing.The feasible solution rate for VQA and QA is calculated as where F is the feasible solution subspace.The average cost function is given as for VSQA and QA. ( C ave is calculated for the feasible solutions for VSQA and QA.A smaller value of the average cost function indicates a higher performance.To systematically determine the optimum value of the constraint coefficient, we set the threshold value of the feasible solution rate to 0.1.We repeatedly solve each problem instance while varying the constraint coefficient.Then we determine the parameter regime that yields a feasible solution rate above the threshold value.Finally, we search for the optimal value of A to minimize the average cost function.The precision threshold of the constraint coefficient is set to 1 for GPPs and 200 for QKPs.The performances of pVSQA, VSQA, pQA, and QA are compared using the success probability averaged over the 10 problem instances for each problem size.The success probability p suc is given as where S is the set of optimum solutions.the continuous schedule for all operation times (Fig. 5 (a)).

Fig
For the QAOA schedule, p = 1 is sufficient for the bangbang regime (T ∼ 10 −1 ), whereas p = 3 is appropriate for the bang-anneal-bang regime (T ∼ 1) to achieve the success probability of the continuous schedule (Fig. 5 (b)).However, the performance of the QAOA schedule degrades in the anneal regime (T ∼ 10) because the time steps in the QAOA schedule are too large to accurately approximate the continuous schedule.Fig. 5 (c) compares the success probability of pVSQA with VSQA, pQA, and QA when the operation time is T = 1 and the number of layers for the QAOA schedule is p = 3. Filled symbols connected by solid lines denote the results For all the problem sizes, the linear and QAOA schedules exhibit almost the same performances.pVSQA shows the best performance among them and significantly outperforms the conventional methods of QA and QAOA (i.e., VSQA for the QAOA schedule).
The simulator results for the QKPs are qualitatively similar to those for the GPPs (Fig. 6).The linear schedule approximates the continuous schedule well for all operation times (Fig. 6 (a)).The QAOA schedule with a small value of p provides a good approximation of the continuous schedule in the bang-bang regime and the bang-anneal-bang regime, but its performance degrades in the anneal regime (Fig. 6 (b)).Fig. 6 (c) compares the success probability of pVSQA with the other methods.Note that we set p = 3 for the QAOA schedule and T = 1.pVSQA outperforms the other methods.
Table 2 shows the optimized values of variational parameters s 1 and s 2 for the linear schedule.GPPs and QKPs display the similar dependences of s 1 and s 2 on the operation time.For a short operation time, s 1 ≃ 1 and s 2 ≃ 0. In contrast, s 1 ≃ 0 and s 2 ≃ 1 for a long operation time.As the operation time increases, s 1 gradually decreases, whereas s 2 gradually increases.These operation time dependences of s 1 and s 2 are consistent with the optimal schedule functions in the continuous schedule (see Fig. 2).

V. EXPERIMENTS WITH QUANTUM DEVICES
We conducted experiments of pVSQA on a quantum annealer and a gate-based quantum device.We use a D-Wave annealer as a quantum annealer.The quantum annealing hardware embeds a maximum of 180 spins on a complete graph [33].As a gate-based quantum device, we use ibmq_brisbane and ibmq_osaka, which are IBM Quantum Canary Processors [34].The number of spins available is 127.Both experiments used Python 3.7.6 as the implementation language.with 50 items and 2 problems with 100 items are generated.We use the benchmarking instances 100_100_i to generate the problem instances [29].Problem instances with 50 items are generated in the same process as Sec.IV.The problem instances are named as QKP_n_i, where n is the number of items and i is the label of the benchmarking instance.The gate-based quantum device uses a 4-node GPP (Fig. 7).The optimal cost function for each problem is obtained by a heuristic Ising-model classical solver [30] or reference [29].

A. INSTANCES OF COPS
For each problem instance, we generate the Ising model using the QUBOs in (30) and (32).For the experiment on the quantum annealer, auto-scaling of the solver is used to adjust the coefficients of the Ising model.For the experiments on the gate-based quantum device, the coefficients of the Ising model are rescaled as shown in Sec.IV.

B. METHOD 1) Implementation on the quantum annealer
We compare the performances of pVSQA, pQA, and QA on the quantum annealer.Here, VSQA is excluded due to its poor performance compared to pVSQA and pQA (Sec.IV).We set the operation time to 5µs.pVSQA adopts the schedule function close to that of the linear schedule.The schedule function is controlled by a parameter called anneal fraction s, which takes a value between 0 and 1.The strength of the transverse fields decreases with s.The values of s at initial time s 0 and at end time s 3 are fixed to 0 and 1, respectively.We use the values of s at 0.5µs and 4.5µs as variational parameters.They are denoted by s 1 and s 2 , respectively.The set of parameters {s i } 3 i=0 produces a piecewise linear schedule function.Two types of classical optimizers are adopted to determine s 1 and s 2 .One uses the Powell optimizer implemented in SciPy, where the initial values of s 1 and s 2 are 0.5 and the maximum number of iterations is 10.The other uses the grid-search optimizer to find the optimal values of s 1 and s 2 .The precision thresholds of s 1 and s 2 are set to 0.1.The probability distribution p 50 (σ) in ( 7) is generated by the top 50 high-quality solutions in terms of the cost function (i.e., C(σ)) out of the 100 solutions obtained by the quantum measurements.pQA and QA use the default schedule function.Namely, the schedule function begins with s = 0 and ends with s = 1 with a linear interpolation between them.
Then we optimize the constraint coefficient A and the chain strength for each problem.The chain strength is a parameter for minor embedding [35], [36].Here, the parameters are determined as follows.We obtain 100 solutions by performing quantum computation along the given schedule function for pQA and QA.We use the top 50 high-quality solutions in terms of the cost function (i.e., C(σ) for pQA and C ′ (σ) for QA) to compute both the feasible solution rate and the average cost function.The feasible solution rate is given by the ratio to obtain the feasible solutions in the 50 solutions.Then a parameter regime is selected to yield a feasible solution rate above the threshold value of 0.1.Afterwards, the optimal values of the constraint coefficient and the chain strength are searched to minimize the average cost function.The average cost function is given by the average of C(σ) and C ′ (σ) over the feasible solutions in the 50 solutions.The precision threshold of the constraint coefficient A is 1 for GPPs and 1000 for QKPs, whereas that of the chain strength is 10.pVSQA uses the same values of A and chain strength as pQA.
Finally, the performance is measured by the success probability and the residual energy.Success probability p suc is the ratio to obtain the optimal solutions in the 50 solutions.Residual energy R ave 50 (R min ) is the difference between the average (minimum) cost function and the optimum one.Smaller values of the residual energies indicate better performances.
In experiments, each algorithm is executed 5 times for each problem instance.In Section V-C, we show the average and the standard deviation of the success probability and the residual energy based on the 5 samples.

2) Implementation on the gate-based quantum device
We compare the performances between pVSQA and VSQA because pQA and QA are difficult to execute on the gatebased quantum device.The quantum circuit for the QAOA schedule with p layers is given as where Ûq (s) := e −i Ĥqs , ÛIsing (s) := e −i ĤIsings , and |ψ(0)⟩ is the ground state of Ĥq .The operators are time-ordered (i.e., the operators with i = 1 are on the right-most, whereas the operators with i = p are on the left-most).It should be noted that various mixer terms, tailored to explore within the feasible solution subspace [37], can be selected as Ûq (s).However, for the purpose of this study, the transverse field (i.e., Ûq (s) = exp (i i∈V σx i s)) is specifically chosen to illustrate the effectiveness of a post-processing approach.We use Qiskit library to build the quantum circuit [38].Here, we adopt the QAOA schedule with T = 1 and p = 1.We set the constraint coefficient A optimized by maximizing the success probability in simulator experiments with the precision threshold of 1.We use the Powell optimizer implemented in SciPy as a classical optimizer with a maximum of 10 iterations.The initial values of the variational parameters (i.e., s 1 and s 2 ) are set to zero.The probability distribution p 100 (σ) for pVSQA in (7) and p ′ 100 (σ) for VQA in (9) are generated by the 50 solutions obtained by the quantum measurements.Finally, the performance on the gate-based quantum device is measured by the success probability from 50 solutions.
In experiments, each algorithm is executed 3 times.Below, we show the average and the standard deviation of the success probability based on the 3 samples.

C. RESULTS
Table 3 displays the performances of pVSQA, pQA, and QA on the quantum annealer.Post-processing significantly improves the performance.The choice of the classical optimizer influences the pVSQA performance.The grid-search optimizer outperforms the Powell optimizer.pVSQA with the grid-search optimizer yields better results than pQA for all problem instances, but pVSQA with the Powell optimizer exhibits performance similar to pQA, particularly for QKPs.pVSQA with the grid-search optimizer finds optimal solutions, except for GPPs with 64 nodes.Even for problems with low (or no) success probability, pVSQA with the gridsearch optimizer produces smaller residual energies (R ave 50 and R min ) than pQA and QA.
Table 4 shows the performances of pVSQA and VSQA on the gate-based quantum device.pVSQA outperforms VSQA.The success probability of pVSQA on the quantum device is slightly worse than that on the simulator.This difference is attributed to the inherent noise in the quantum device.

VI. DISCUSSION
The simulator result shows that the linear schedule almost reproduces the success probability of the continuous schedule in both GPPs and QKPs.This is expected with short or long operation times.Because the quantum state does not change significantly from the initial one when the operation time is short, the success probability is close to the random sampling with a post-processing.This situation is almost independent of the schedule function.When the operation time is long, the conventional quantum annealing schedule (i.e., linear schedule with s 0 = 0 and s 1 = 1) achieves the ground state.Consequently, a small modification of the schedule function does not affect the performance.However, this result is nontrivial for an intermediate operation time.Although the optimum schedule function for the continuous schedule is a complicated bang-anneal-bang type (see Fig. 2 (b)), a simple linear schedule gives a good approximation.The QAOA schedule also provides a suitable approximation of the continuous schedule in the bang-bang regime and the bang-anneal-bang regime.These results indicate that a small number of variational parameters is sufficient to implement pVSQA.
The experimental results on a quantum annealer demonstrate the effectiveness of our proposed algorithm in both GPP and QKP.To investigate these results in more detail, we illustrate the grid-search results of R ave 50 for both GPP and QKP in Fig. 8.They are rescaled in units of their respective σ values, the standard deviations of R ave 50 averaged over s 1 and s 2 .Circles and squares represent the optimal values of s 1 and s 2 obtained using the grid-search optimizer and the Powell optimizer, respectively.Triangles correspond to s 1 and s 2 values in pQA (Note that the default schedule function is used in pQA experiments and s 1 and s 2 are not explicitly specified).For GPPs, the upper-left regime (i.e., small s 1 and large s 2 ) yields lower energies for 32 nodes but higher energies for 64 nodes, which is consistent with the simulator result.As the problem size increases, the crossover between the bang-bang regime and the anneal regime generally occurs at longer operation times.In other words, for a given operation time, the crossover from the anneal regime to the bang-bang regime appears as the system size increases.This observation indicates that the pVSQA is more effective than pQA when the problem size is large.For QKPs, the Powell optimizer fails to produce higherquality solutions than pQA.Fig. 8 reveals that the ranges of R ave 50 , representing the difference between the maximum and the minimum values, in the rescaled unit for QKPs are smaller than those for GPPs.We hypothesize that the smaller ranges contribute to the relatively poor performances of the Powell optimizer for QKPs.Even in this situation, the grid-search optimizer consistently obtains lower-energy solutions.It is imperative to develop a classical optimizer that efficiently finds the optimal variational parameters on quantum annealers.
The gate-based quantum device demonstrates that pVSQA outperforms the conventional QAOA.However, the experiments in this study are limited to a small-sized problem due to our restricted execution time on quantum devices.In the future, the pVSQA performance should be evaluated on large-sized problems.

VII. CONCLUSION
We proposed pVSQA to efficiently solve constrained COPs.We provided a sufficient condition for constrained COPs to apply pVSQA based on a greedy post-processing algorithm.To demonstrate its feasibility for the state-of-the-art quantum annealer and gate-based quantum device, we applied it to two NP-hard COPs: GPP and QKP.pVSQA on a simulator achieves a (near-)optimal performance within a predetermined operation time using a small number of variational parameters.Then we experimentally verified the superior performances of pVSQA on real quantum devices compared to conventional quantum annealing with and without a postprocessing and QAOA without a post-processing.

FIGURE 2 :FIGURE 3 :
FIGURE 2: (a)-(c): Optimal schedule functions in the continuous schedule for different operation times of T = 0.1, 1, and 10, respectively.Each line corresponds to the optimal schedule function for a different GPP.(d): Success probability as a function of operation time.Points in Fig.(a)-(c) are marked.

Fig. 2 (
a)-(c) shows the three types of optimal schedule functions for 10 GPPs with 8 nodes for different values of the operation time.Fig. 2 (d) plots the success probability as a function of the operation time.Here, the success probability denotes the probability of obtaining the optimum solution of a GPP when quantum computation is executed along the optimal schedule function and a subsequent post-processing without a renormalization (i.e., r = 100 in (7)) is performed (see Sec. III-C for details).The success probability increases with operation time.Short operation times produce a bangbang-type schedule function.Quantum computation starts with the Ising Hamiltonian and quickly transitions to the transverse-field Hamiltonian.The transition occurs at the midpoint of the operation time.Intermediate operation times result in a bang-anneal-bang-type schedule function.Quantum computation begins with the Ising Hamiltonian, which follows a continuous function that ends with the transversefield Hamiltonian.Long operation times generate schedule functions that closely resemble the conventional quantum annealing, which start with the transverse-field Hamiltonian and end with the Ising Hamiltonian.Below, we call the short, intermediate, and long operation time regimes the bangbang, bang-anneal-bang, and anneal regimes, respectively.Reference [23] reported similar observation.
obj and ∆Q m obj are respectively introduced as the maximum energy change of Qobj and Q obj when flipping a binary variable in V m .The independent constraints indicate ∆ Qm obj = ∆Q m obj , and then the first condition of Theorem 2 gives second condition (i.e., b m max −b m min ≥ max i∈V m |a m i |−1) always holds and the fourth condition implies k ≤ |V m |.The constraint Hamiltonian for the k-hot constraint is given as

5 :j
∆Q i ← Energy change of Q ′ when flipping x i ← min(arg min i∈V ∆Q i ) 8:

FIGURE 4 :
FIGURE 4: Number of flips to obtain a local minimum solution in the greedy post-processing algorithm as a function of |V | for GPPs.The solid line and shaded region represent the mean and the standard deviation over 1000 random spin configurations for each of 10 different GPPs, respectively.
generate QKPs with n items, we use a profit matrix {p .Here, ⌊x⌋ := max{m ∈ Z|m ≤ x} is the usual floor function.The optimum solution for each problem is obtained by a state-ofthe-art heuristic Ising-model classical solver

FIGURE 5 :
Fig.5shows the results for the GPPs.The success probability for the linear schedule increases with the operation time.The linear schedule almost reproduces the success probability of

FIGURE 6 :
FIGURE 6: Results for QKPs.Success probability of 8 item QKPs as a function of operation time for (a) the linear schedule and (b) QAOA schedule.Continuous schedule is depicted by pentagons (black) for comparison.(c): Dependences of the success probability on the number of items for pVSQA, VSQA, pQA, and QA.In pVSQA and VSQA, data for the linear schedule and QAOA schedule are represented by filled and open symbols, respectively.

TABLE 2 :
Optimized values of variational parameters s 1 and s 2 for the linear schedule.Values in parentheses denote the standard deviation of the average.

4 FIGURE 7 :
FIGURE 7: Four-node GPP used in experiments on the gatebased quantum device.Dashed line denotes the optimum solution.

TABLE 3 :
Experimental results of pVSQA, pQA, and QA for GPPs and QKPs on the quantum annealer.Average residual energy, minimum residual energy, success probability, and feasible solution rate are denoted by R ave 50 , R min , p suc , and p FS , respectively.Values in parentheses denote the standard deviation of the average.Average indicates the average of R ave 50 , R min , and p suc over the GPPs and QKPs with the same problem size.Opt denotes the cost function of the optimal solution.

TABLE 1 :
Type of constraints.

TABLE 4 :
Success probability of pVSQA and VSQA for the GPP on the gate-based quantum device (qpu) and the simulator (cpu).Values in parentheses denote the standard deviation of the average.