QUANTUM CONSTRAINED HAMILTONIAN OPTIMIZATION

Abstract

A device includes applying coupling and transformation operations to quantum states according to a Hamiltonian specification. Information is received at a digital computer based in part on measurements of the quantum states. The digital computer provides information for preparing quantum states associated with quantum processing elements based in part on the information. A control module applies coupling and transformation operations based on interaction with the digital computer for processing the constrained optimization problem. The processing includes preparing quantum states associated with quantum processing elements characterized by a summation of a constraint Hamiltonian and an objective Hamiltonian. The processing further includes operating the control module to evolve a time-dependent Hamiltonian by forming a sum of a first term having the constraint Hamiltonian and a second term. The second term includes a product that is initially equal to the objective Hamiltonian and is evolved into a negative of the objective Hamiltonian.

Description

BACKGROUND OF THE INVENTION

Constrained combinatorial optimization problems represent a class of computational challenges where the goal is to find the optimal solution from a finite set of feasible solutions, subject to a set of constraints. These problems are characterized by discrete decision variables and a multitude of constraints that limit the feasible solution space. The objective is to maximize or minimize an objective function while adhering to the specified constraints. Common examples include the traveling salesman problem, job scheduling, and network design. The complexity of these problems arises from the discrete nature of the decision variables, leading to a combinatorial explosion of possible solutions and requiring sophisticated algorithms for efficient exploration of the solution space.

Addressing constrained combinatorial optimization problems involves the development and application of various optimization techniques. Metaheuristic algorithms, such as genetic algorithms, simulated annealing, and ant colony optimization, are commonly employed due to their ability to navigate large solution spaces and find near-optimal solutions. These algorithms iteratively explore and exploit the solution space, balancing the trade-off between exploration for diverse solutions and exploitation to refine the current best-known solution. Additionally, mathematical programming techniques, such as integer programming and constraint programming, offer formal frameworks for modeling and solving constrained combinatorial optimization problems. These methods provide a systematic approach to formulate the problem, express constraints, and leverage optimization solvers to find optimal or near-optimal solutions.

Despite the advancements in optimization algorithms, constrained combinatorial optimization problems remain challenging due to their NP-hard nature. Efficiently solving large instances of these problems requires a combination of algorithmic innovations, problem-specific heuristics, and computational resources. The complexity and ubiquity of constrained optimization problems motivates the development of quantum algorithms as a new approach to tackle them.

BRIEF DESCRIPTION OF THE DRAWINGS

For a more complete understanding of the invention, reference is made to the following description and accompanying drawings, in which:

FIG. 1 is a diagram of an example hybrid quantum-classical computing system in accordance with respective examples.

FIG. 2 illustrates summary data of MIS simulation results in accordance with respective examples.

FIG. 3 illustrates summary data of DMDS simulation results with random directed graphs in accordance with respective examples.

FIG. 4 illustrates summary data of knapsack simulation results in accordance with respective examples.

FIG. 5 illustrates histograms results for a specific example of the knapsack problem in accordance with respective examples.

FIG. 6 illustrates summary data of combinatorial auction simulation results in accordance with respective examples.

FIG. 7 illustrates histograms of combinatorial simulation results for a specific example of a combinatorial auction in accordance with respective examples.

FIG. 8 is a flowchart for processing a constrained optimization problem in accordance with respective examples.

DETAILED DESCRIPTION

The invention can be implemented in numerous ways, including as a process; an apparatus; a system; a composition of matter; a computer program product embodied on a computer readable storage medium; and/or a processor, such as a processor configured to execute instructions stored on and/or provided by a memory coupled to the processor. In this specification, these implementations, or any other form that the invention may take, may be referred to as techniques. In general, the order of the steps of disclosed processes may be altered within the scope of the invention. Unless stated otherwise, a component such as a processor or a memory described as being configured to perform a task may be implemented as a general component that is temporarily configured to perform the task at a given time or a specific component that is manufactured to perform the task. As used herein, the term ‘processor’ refers to one or more devices, circuits, and/or processing cores configured to process data, such as computer program instructions.

A detailed description of one or more examples of the invention is provided below along with accompanying figures that illustrate various features of the invention. The invention is described in connection with such examples, but the invention is not limited to any example. The scope of the invention is limited only by the claims and the invention encompasses numerous alternatives, modifications and equivalents. Numerous specific details are set forth in the following description in order to provide a thorough understanding of the invention. These details are provided for the purpose of example and the invention may be practiced according to the claims without some or all of these specific details. For the purpose of clarity, technical material that is known in the technical fields related to the invention has not been described in detail so that the invention is not unnecessarily obscured.

Described herein is an adiabatic quantum algorithm, quantum constrained Hamiltonian optimization (Q-CHOP), that may be used for solving constrained optimization problems. In various implementations, Q-CHOP enforces a Hamiltonian constraint at all times, and constructs an adiabatic path directly out of the problem's objective function to convert the worst feasible state into the best feasible state. In contrast to the standard approaches to adiabatic quantum optimization, Q-CHOP constructs the entire adiabatic schedule out of optimization problem data, without invoking an ad-hoc driver Hamiltonian or initial state that have no relation to the problem at hand. Q-CHOP further comes equipped with a natural choice for the constraint-enforcing “penalty factor”, leaving only the quantum runtime as a hyperparameter that may need careful tuning for performance. Q-CHOP consistently outperforms the standard, penalty-based adiabatic algorithm in numerical experiments across a diverse set of problems, suggesting that Q-CHOP is a promising general-purpose method for constrained combinatorial optimization.

Referring to FIG. 1, an example hybrid quantum-classical computing system 100 includes a quantum processor 102 comprising a plurality of quantum processing elements associated with respective quantum states. In some examples, the quantum processor 102 may be a physical quantum processor. In other examples, the quantum processor 102 may be a simulation of a physical quantum processor. The quantum processor 102 is configured to apply coupling operations (e.g., multi-qubit operations) and transformation operations (e.g., single-qubit operations) to a plurality of the quantum states according to a quantum circuit specification (e.g., QC specification 108 stored on a storage medium 110) defining a plurality of quantum gate operations. A digital computer 104 comprises at least one central processing unit, and is configured to: receive information based at least in part on measurements of one or more quantum states associated with respective quantum processing elements of the quantum processor 102, and provide information for preparing one or more quantum states associated with respective quantum processing elements of the quantum processor 102 based at least in part on the received information. A control module 106 is configured to control the applied coupling and transformation operations (e.g., using radio frequency electromagnetic fields coupled to a coupled array of superconducting circuits, or using optical electromagnetic fields coupled to atoms) based on interaction with the digital computer 104 for processing a constrained optimization problem. The processing comprises: preparing quantum states associated with a plurality of the quantum processing elements characterized by a first Hamiltonian representing a constraint of the constrained optimization problem, preparing quantum states associated with a plurality of the quantum processing elements characterized by a second Hamiltonian representing an objective function of the constrained optimization problem, and operating the control module to evolve a time-dependent Hamiltonian according to an evolution that includes: (1) forming an expression that is initially equal to a product of the first Hamiltonian and at least one time-dependent function that includes the second Hamiltonian as an argument of an exponential function, or (2) forming a sum of a first term comprising the first Hamiltonian and a second term, where the second term comprises a product of (A) a time-dependent scalar function and (B) a time-dependent operator that is initially equal to the second Hamiltonian and is evolved into the negative of the second Hamiltonian.

In some implementations, the quantum computing algorithm represented by the quantum circuit specification is used for evolving the time-dependent Hamiltonian over a first iteration comprising a sequence of discrete time intervals. These discrete time intervals can be variationally optimized using a processor (e.g., a classical processor such as the processor core(s) of the digital computer 104, or a quantum processor, or a hybrid quantum-classical computing system that includes both a classical processor and a quantum processor). For example, after evolving the time-dependent Hamiltonian over the first iteration, a variational optimization of the second Hamiltonian may be performed to determine a second iteration comprising a second sequence of discrete time intervals. The variational optimization may include performing gradient descent, Adam, etc. The same sequence of discrete time intervals may be evolved repeatedly in order to generate an expectation value of the second Hamiltonian for a given sequence of discrete time intervals. In such a case, the variational optimization would optimize the expectation value of the second Hamiltonian.

As way of background, consider a constrained combinatorial optimization problem of the form:

$\begin{array}{cc}\mathrm{minimize}:f\left(x\right)& \left(1\right)\end{array}$ $\mathrm{subject}\mathrm{to}:C_i\left(x\right)=0\mathrm{for}\mathrm{all}i.$ $x\in {\left\{0,1\right\}}^N\mathrm{minimize}:f\left(x\right)$ $\mathrm{subject}\mathrm{to}:C_i\left(x\right)=0\mathrm{for}\mathrm{all}i.$ $x\in {\left\{0,1\right\}}^N$

Here f is the objective function and C_i are constraints. It is assume without loss of generality that the all constraint functions are nonnegative, which implies that they are minimized by constraint-satisfying (feasible) bitstrings; otherwise, C_i(x) can be replaced in Eq. (1) by {tilde over (C)}_i(x)=C_i(x)². The treatment of inequality constraints is discussed further below.

Under the assumption of nonnegative constraints, the optimization problem in Eq. (1) can be encoded into a system of qubits by the objective and constraint Hamiltonians:

$\begin{array}{cc}{H}_{\mathrm{obj}}=\sum _{x}f\left(x\right)❘x\rangle \langle x❘,H_{\mathrm{con}}=\sum _{i,x}{C}_i\left(x\right)❘x\rangle \langle x❘,& \left(2\right)\end{array}$

thereby associating “good” bitstrings of low objective-function value with low-energy states of H_obj, and identifying the set of feasible bitstrings with the ground-state subspace of H_con. Thusly encoded, the optimization problem in Eq. (1) can be solved via the adiabatic algorithm.

To this end, the standard procedure is to first construct a problem Hamiltonian H_prob=H_obj+λH_con with a penalty factor for constraint violation, A, that is sufficiently large to ensure that the ground state of H_prob is feasible. One then chooses a driver Hamiltonian H_driver that satisfies two properties: (i) it has an easy-to-prepare ground state, and (ii) it does not commute with H_prob. Finally, the quantum algorithm prepares the ground state of H_driver, and slowly changes the Hamiltonian of the system from H_driver to H_prob over the course of a total time T, for example via the interpolating Hamiltonian:

$\begin{array}{cc}{H}_{\mathrm{interp}}\left(s\right)=\left(1-s\right)H_{\mathrm{driver}}+{\mathrm{sH}}_{\mathrm{prob}},& \left(3\right)\end{array}$

sweeping s: 0→1 as the time t: 0→T.

T is referred to as the “quantum runtime” of the algorithm. Property (i) is necessary to ensure that the algorithm solves the optimization problem by adiabatically tracking the ground state of H_interp, while property (ii) ensures that H_interp has off-diagonal terms that induce nontrivial dynamics.

Choosing an interpolating Hamiltonian of the form in Eq. (3) is not required, and indeed sometimes different paths from H_driver to H_prob can improve performance, in the extreme case turning algorithmic failure into success. For simplicity, however, the interpolating Hamiltonian is considered. This approach is referred to as the standard adiabatic algorithm (SAA).

The adiabatic theorem ensures that the SAA prepares an optimal solution in time T that scales with an inverse power of the smallest spectral gap Δ_interp^min of H_interp(minimized over s∈[0, 1]), provided that there are no symmetries that fragment Hilbert space into uncoupled sectors. However, the gap Δ_interp^min is usually not known in advance, so in practice the SAA requires making an ad hoc choice of some “long” time T, and hoping that T is sufficiently large. If the selected time is too small, the SAA will not only fail to provide an optimal solution, it may even fail to provide a feasible solution.

There are a few additional shortcomings of the SAA that are worth noting. First, the SAA requires making a careful choice of the penalty factor λ. Increasing λ makes it more likely that the SAA finds feasible solutions to the optimization problem at hand. However, doing so also increases the runtime of the algorithm, or equivalently the depth of a quantum circuit that implements this algorithm with fixed error. Making an optimal choice of λ has been shown to be a difficult task in practice and NP-hard in general, though known heuristic strategies exist that work well in certain cases.

Second, the SAA relies on an ad hoc initial state and driver Hamiltonian H_driver that has no connection to the problem of interest. A common choice for the initial state is the uniform superposition over all bitstrings (representing the uniform prior), and a transverse magnetic field for which the uniform superposition is a unique ground state. However, even here one must choose an overall normalization (prefactor) for H_driver, analogous to the choice of penalty factor λ. Finally, there is an intuitive sense in which the SAA does “too much work”: it searches the configuration space of all bitstrings, feasible or otherwise.

Described below are various examples of quantum constrained Hamiltonian optimization (Q-CHOP) and its application to problems with various objectives and constraints. The Pauli-Z operator for qubit j∈{1, 2, . . . , N} is denoted by Z_j=|00|_j−|11|_j, the Pauli-X operator by X_j=|01|_j+|10|_j, and the Pauli-Y operator by Y_j=−iZ_jX_j. The generator of a global rotation about axis α∈{x, y, z} on the Bloch sphere is denoted by S_α, e.g.

$S_z=\frac{1}{2}{\sum }_j{Z}_j,$

such that R_α(θ)=e^−iθSα rotates all qubits about axis α by the angle θ.

The general strategy for Q-CHOP is as follows:

• • 1. Identify an initial Hamiltonian H_init whose energy within the space _feas of feasible states is minimized by an easy-to-prepare state |x_init∈_feas. Prepare |x_init.
  • 2. Construct a ramp H_ramp(s) that traces a continuous path from H_init to the objective Hamiltonian H_obj, with H_ramp(0)=H_init and H_ramp(1)=H_obj.
  • 3. Pick a time T and a positive real number A, and evolve from time t=0 to t=T under the total Hamiltonian

$\begin{array}{cc}{H}_{\mathrm{tot}}\left(s\right)=H_{\mathrm{con}}+{\lambda }^{-1}{H}_{\mathrm{ramp}}\left(s\right),& \left(4\right)\end{array}$

• •  with s: 0→1 as t: 0→T.In various examples below, explicit choices of H_init and H_ramp that are constructed out of the objective Hamiltonian H_obj are provided. However, general features of the strategy in steps 1-3 are discussed first.

Intuitively the constraint Hamiltonian H_con will restrict dynamics to the feasible subspace _feas, while the ramp H_ramp will adiabatically transfer the initial state |x_init to an optimal (or approximate) solution to the optimization problem. In this general form, Q-CHOP shares many features with the SAA. The scalar A is analogous to the penalty factor in the SAA, and indeed the Hamiltonians of Q-CHOP and the SAA at s=1 are equal (up to an overall normalization factor). Similarly, the choices of an initial Hamiltonian H_init and a ramp H_ramp for Q-CHOP are analogous to choices of a driver Hamiltonian H_driver and a path from H_driver to H_prob in the SAA.

Though this general form of Q-CHOP does not practically reduce the number of choices necessary to solve the constrained optimization problem, it effectively reduces the size of the classical search space of the quantum algorithm, which can reduce requisite runtimes. Moreover, Q-CHOP provides a theoretical framework for making a suitable choice of λ. Denote the operator norm of M by ∥M∥, let ∥H_ramp∥=max_s∈[0,1] ∥H_ramp(s)∥, and denote the spectral gap of H_con by Δ_con. If

$\begin{array}{cc}H_{\mathrm{ramp}}/\lambda <{\Delta }_{\mathrm{con}}/2,\mathrm{and}& \left(a\right)\end{array}$ $\begin{array}{cc}{H}_{\mathrm{ramp}}\mathrm{changes}\mathrm{slowly}\mathrm{with}\mathrm{respect}\mathrm{to}{\Delta }_{\mathrm{con}}^{-1},& \left(b\right)\end{array}$

then second-order perturbation theory guarantees that Q-CHOP will yield feasible states with probability

$\begin{array}{cc}{p}_{\mathrm{feas}}=1-O\left({\delta }^2\right),\mathrm{where}\delta =\frac{H_{\mathrm{ramp}}}{{\mathrm{\lambda \Delta }}_{\mathrm{con}}}.& \left(5\right)\end{array}$

Intuitively, condition (a) ensures that the H_ramp is not “strong enough” to violate constraints, while condition (b) precludes rapid (e.g., Floquet) driving schemes that resonantly pump energy into the system. The spectral gap Δ_con can typically be determined (or at least lower-bounded) by inspection of the constraint Hamiltonian H_con.

Larger values of penalty factor λ provide better guarantees of feasible solutions, while smaller values of λ are desirable to increase the weight of the ramp H_ramp that drives evolution within the feasible subspace. Condition (a) and the guarantee in Eq. (5) provide an upper bound on “reasonable” choices of λ, namely λ_max≈2∥H_ramp∥/Δ_con. In practice, even the upper bound λ≤λ_max is too stringent, as the above discussion still holds if the norm ∥H_ramp(s)∥ is replaced by ∥(1−P_feas)H_ramp(S)P_feas∥, where P_feas is a projector onto the feasible subspace. One should therefore expect smaller values of λ to still yield predominantly feasible states, although the limits on λ may need to be found empirically, since the norm ∥(1−P_feas)H_rampP_feas∥ is generally not known a priori. Moreover, the ramp that constructed and described below has the feature that ∥H_ramp∥≤∥H_obj∥. The fact that a suitable choice of λ can be determined from ∥H_obj∥ and Δ_con directly, rather than from a significantly more complex interpolating Hamiltonian in the case of the SAA, is a notable advantage of Q-CHOP, since these quantities are typically easy to estimate or bound appropriately. A uniform linear objective

$H_{\mathrm{obj}}=S_z=\frac{1}{2}{\sum }_j{Z}_j\simeq -{\sum }_j{❘1\rangle \langle 1❘}_j$

(where ≃ denotes equality up to a constant), for example, has norm ∥H_obj∥=N/2, so if Δ_con=1 then λ=N is a reasonable choice. However, in the worst case, determining the optimal choice of λ may be NP-hard as is the case for the SAA.

It is often the case that the objective function of an optimization problem is linear in the decision variables of the problem, in which case the corresponding objective Hamiltonian can be written (up to an irrelevant constant shift) as a weighted sum of single-qubit Pauli-Z operators. This Hamiltonian is then odd with respect to qubit inversion, which is to say that x|H_obj|x=−x|H_obj|x, where x=(x₁, x₂, . . . )∈{0, 1}^N is a bitstring and x=(1−x₁, 1−x₂, . . . ) is its bitwise complement. It is also sometimes the case that finding good solutions to a constrained optimization problem is difficult, but finding the worst feasible solution is easy, either by inspection or with an efficient classical algorithm. This is the case for the problem of finding a maximum independent set (MIS), whose worst feasible state is the all-0 bitstring. When the objective is odd, the asymmetry between best and worst feasible states motivates the choice of initial Hamiltonian H_init=−H_obj, and a ramp that continuously flips all qubits: H_ramp(θ)=−H_obj(θ), with

$\begin{array}{cc}{H}_{\mathrm{obj}}\left(\theta \right)=R_y\left(\theta \right)H_{\mathrm{obj}}{R}_y{\left(\theta \right)}^{†}& \left(6\right)\end{array}$

where parameterising the ramp is by the angle θ=sπ for brevity, and R_y(θ)=e^−iθSy rotates all qubits about the y axis by θ.

When these two conditions hold, namely objective Hamiltonian is odd worst feasible state |x_worst is easy-to-prepare, the Q-CHOP can be done as follows. Prepare |x_worst, pick a quantum runtime T and positive number λ≈2∥H_obj∥/Δ_con, and evolve from time t=0 to t=T under the total Hamiltonian

$\begin{array}{cc}{H}_{\mathrm{tot}}\left(\theta \right)=H_{\mathrm{con}}-{\lambda }^{-1}{H}_{\mathrm{obj}}\left(\theta \right),& \left(7\right)\end{array}$

sweeping θ: 0→π as t: 0→T, for example with the linear ramp θ=πt/T. Here ∥H_obj∥ is the operator norm of H_obj, and Δ_con can be replaced by an estimate or lower bound on the spectral gap of H_con. The prescription for choosing λ is only a guideline, and smaller values of λ may work well in practice. λ was set equal to the qubit number N in all examples described herein.

In the case of MIS, Q-CHOP may reduce to the specialized algorithm based on non-Abelian adiabatic mixing, up to an adiabatically vanishing addition of ∂_tθS_y to H_tot(θ) in Eq. (7). However, the general formalism of Q-CHOP extends to more general constrained optimization problems as well. Intriguingly, the addition of ∂_tθS_y can be identified with approximate counterdiabatic driving, and its form depends only on the linearity of the MIS objective function. An added term of ∂_tθS_y is therefore expected to marginally improve the performance of Q-CHOP for any optimization problem with linear objectives, and similar approximate counterdiabatic driving terms can likely be derived for non-linear objectives as well.

In various examples, Q-CHOP may be sped up for problems with single-body objective Hamiltonians of the form:

$\begin{array}{cc}{H}_{\mathrm{obj}}=\frac{1}{2}{\sum }_j{c}_j{Z}_j,& \left(8\right)\end{array}$

where {c_j} are arbitrary real coefficients. For convenience and clarity, H_obj,z=H_obj, as well as rotated versions H_obj,x and H_obj,y in which the Pauli-Z operators in H_obj,z are respectively replaced by Pauli-X and Pauli-Y operators. The restriction of Hamiltonians to the feasible subspace is denoted as _feas by a superscript, such that H_obj,α^feas=P_feasH_obj,αP_feas, where P_feas is a projector onto _feas.

Q-CHOP may rely on the adiabatic theorem to track the ground state of a Hamiltonian H(θ) by slowly varying a parameter θ. Evolving too quickly introduces errors due non-adiabatic effects. These effects can be mitigated by the introduction of a counterdiabatic drive, in which the Hamiltonian is replaced according to

$\begin{array}{cc}H\to H_{\mathrm{CD}}=H+{\partial }_t\theta G,& \left(9\right)\end{array}$

where the gauge field G can be found by minimizing the Hibert-Schmidt norm of the operator

$\begin{array}{cc}𝒪\left(G,H\right)={\partial }_{\theta }H+i\left[G,H\right].& \left(10\right)\end{array}$

Here [G, H]=GH−HG, and the (squared) Hilbert-Schmidt norm of is ∥∥_HS²=tr().

When time evolution is sufficiently slow to avoid exciting states outside the feasible subspace, Q-CHOP effectively evolves from time t=0 to t=T under the Hamiltonian

$\begin{array}{cc}{H}_{\mathrm{tot}}^{\mathrm{feas}}\left(\theta \right)=-{\lambda }^{-1}{H}_{\mathrm{obj}}^{\mathrm{feas}}\left(\theta \right),& \left(11\right)\end{array}$ $\mathrm{with}$ $\begin{array}{cc}{H}_{\mathrm{obj}}^{\mathrm{feas}}\left(\theta \right)=P_{\mathrm{feas}}{R}_y\left(\theta \right)H_{\mathrm{obj}}{R_y\left(\theta \right)}^{†}{P}_{\mathrm{feas}},& \left(12\right)\end{array}$

where the angle θ: 0→π as t: 0→T. As H_tot^feas(θ) is the Hamiltonian to augment with a counterdiabatic drive, next is to find a gauge field G that minimizes the norm of (G, H_tot^feas) which is equivalent to minimizing the norm of (G, H_obj^feas). To this end, the following expansion may be used:

$\begin{array}{cc}{H}_{\mathrm{obj}}^{\mathrm{feas}}\left(\theta \right)=\cos \left(\theta \right)H_{\mathrm{obj},z}^{\mathrm{feas}}+\sin \left(\theta \right)H_{\mathrm{obj},x}^{\mathrm{feas}}& \left(13\right)\end{array}$

and in turn

$\begin{array}{cc}{\partial }_{\theta }{H}_{\mathrm{obj}}^{\mathrm{feas}}\left(\theta \right)=-\sin \left(\theta \right)H_{\mathrm{obj},z}^{\mathrm{feas}}+\cos \left(\theta \right)H_{\mathrm{obj},x}^{\mathrm{feas}}.& \left(14\right)\end{array}$

It is unclear how to minimize the norm of (G, H_obj^feas) exactly, so minimizing it variationally with respect to the ansatz may be used:

$\begin{array}{cc}{G}_g=g_x{S}_x^{\mathrm{feas}}+g_y{S}_y^{\mathrm{feas}}+g_z{S}_z^{\mathrm{feas}},& \left(15\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{S}_{\alpha }^{\mathrm{feas}}=P_{\mathrm{feas}}{S}_{\alpha }{P}_{\mathrm{feas}},& \left(16\right)\end{array}$

for which

$\begin{array}{cc}𝒪\left(G_g,H_{\mathrm{obj}}^{\mathrm{feas}}\right)={\partial }_{\theta }{H}_{\mathrm{obj}}^{\mathrm{feas}}\left(\theta \right)+i\cos \left(\theta \right)\left[G_g,H_{\mathrm{obj},z}^{\mathrm{feas}}\right]+i\sin \left(\theta \right)\left[G_g,H_{\mathrm{obj},x}^{\mathrm{feas}}\right].& \left(17\right)\end{array}$

Nonzero values of g_x and g_z only contribute to imaginary components to (G_g, H_obj^feas), and conversely all imaginary components of (G_g, H_obj^feas) involve factors of g_x or g_z. To minimize the Hilbert-Schmidt norm of (G_g, H_obj^feas), set g_x=g_z=0 and define G_y=g_yS_y^feas, leaving:

$\begin{array}{cc}𝒪\left(G_y,H_{\mathrm{obj}}^{\mathrm{feas}}\right)={\partial }_{\theta }{H}_{\mathrm{obj}}^{\mathrm{feas}}\left(\theta \right)+{\mathrm{ig}}_y\cos \left(\theta \right)\left[S_y^{\mathrm{feas}},H_{\mathrm{obj},z}^{\mathrm{feas}}\right]+{\mathrm{ig}}_y\sin \left(\theta \right)\left[S_y^{\mathrm{feas}},H_{\mathrm{obj},x}^{\mathrm{feas}}\right].& \left(18\right)\end{array}$

The projector P_feas commutes with H_obj,z (both are diagonal in the computational basis), which allows a simplification:

$\begin{array}{cc}\left[S_y^{\mathrm{feas}},H_{\mathrm{obj},z}^{\mathrm{feas}}\right]=P_{\mathrm{feas}}\left[S_y,H_{\mathrm{obj},z}\right]P_{\mathrm{feas}}={\mathrm{iH}}_{\mathrm{obj},x}^{\mathrm{feas}}.& \left(19\right)\end{array}$

Further substituting ∂_θH_obj^feas(θ) from Eq. (14) into Eq. (18), provides:

$\begin{array}{cc}𝒪\left(G_y,H_{\mathrm{obj}}^{\mathrm{feas}}\right)=\cos \left(\theta \right)\left(1-g_y\right)H_{\mathrm{obj},x}^{\mathrm{feas}}-\sin \left(\theta \right)\left(H_{\mathrm{obj},z}^{\mathrm{feas}}-{\mathrm{ig}}_y\left[S_y^{\mathrm{feas}},H_{\mathrm{obj},x}^{\mathrm{feas}}\right]\right).& \left(20\right)\end{array}$

The sin(θ) term vanishes as θ→0 or π, at which point (g_yS_y^feas, H_obj^feas) is minimized by g_y=1. Near θ=π/2, the sin(θ) term dominates, but here no definite statements about the norm of (G_y, H_obj^feas) can be made without knowledge of the feasible subspace. Nonetheless, given that i [S_y, H_obj,x]=H_obj,z, i [S_y^feas, H_obj,x^feas] is generally expected to be “H_obj,z^feas-like” so a minimum near g_y˜1 is expected. In some examples, the gauge field G=S_y is a reasonable (albeit sub-optimal) choice to mitigate non-adiabatic effects. In the case of linear objectives, the total Hamiltonian for Q-CHOP with approximate counterdiabatic driving thus becomes

$\begin{array}{cc}{H}_{\mathrm{tot}}\to H_{\mathrm{tot},\mathrm{CD}}=H_{\mathrm{con}}-{\lambda }^{-1}{H}_{\mathrm{obj}}\left(\theta \right)+{\partial }_t\theta S_y.& \left(21\right)\end{array}$

The small factor of ∂_tθ˜1/T prevents ∂_tθ S_y from generating a population transfer to non-feasible states, so ∂_tθ S_y≈∂_tθ S_f^feas.

In some examples, the requirement of an odd objective used in Sec. ?? is relaxed and a modification of the algorithm that no longer requires the worst feasible state to be easy-to-prepare. In the case of a non-odd objective, the prescription in Eq. (6) needs modification to ensure that the rotated objective Hamiltonian changes sign when θ: 0→π. To this end, the objective can be decomposed into Pauli strings as

$\begin{array}{cc}{H}_{\mathrm{obj}}=\frac{1}{2}\sum _{S\subset {\mathbb{Z}}_N}{c}_S\prod _{j\in S}{Z}_j,& \left(22\right)\end{array}$

where _N={1, 2, . . . , N} is the set of qubit indices. Next one qubit in each term is rotated, and average over each choice of qubits, altogether defining

$\begin{array}{cc}{H}_{\mathrm{obj}}\left(\theta \right)=\frac{1}{2}\sum _{S\subset {\mathbb{Z}}_N}\frac{c_S}{❘S❘}\sum _{j\in S}{R}_y\left(\theta \right)Z_j{R_y\left(\theta \right)}^{†}\prod _{k\in S\backslash \left\{j\right\}}{Z}_k,& \left(23\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{R}_y\left(\theta \right)Z_j{R_y\left(\theta \right)}^{†}=\cos \left(\theta \right)Z_j+\sin \left(\theta \right)X_j.& \left(24\right)\end{array}$

In some examples, the objective is decomposed into even and odd terms, rotate the odd terms globally as in Eq. (6), and rotate the even terms “locally” as in Eq. (23). These examples may have the benefit of reduced computational overheads when discretizing Q-CHOP for gate-based quantum computation. In some examples, the terms of a Hamiltonian may be made up of a sum of components.

In some examples of Q-CHOP, the worst feasible state of a constrained optimization problem is easy to prepare is required. In other examples, however, this requirement is relaxed to that of having at least one feasible state that is easy to prepare. For example, let |x_★ be an easy-to-prepare eigenstate of H_obj, and let E_★=x_★|H_obj|x_★ be the corresponding objective energy. Define a subproblem with the same constraint Hamiltonian as the original optimization problem, H_con, but a new objective Hamiltonian

$\begin{array}{cc}{H}_{\mathrm{obj}}^{★}=-{\left(H_{\mathrm{obj}}-E_{★}\right)}^2,& \left(25\right)\end{array}$

for which |x_★ is the worst feasible state. By construction, H_obj^★ is minimized within the feasible subspace by either the best feasible state |x_best or worst feasible state |x_worst of the original optimization problem, depending on which corresponding objective energy (with respect to H_obj) is farther from E_★. If the optimization subproblem defined by (H_obj^★, H_con) is solved by |x_best, then after solving the subproblem the process is complete. Otherwise, the worst feasible state |x_worst identified by the subproblem to solve the original optimization problem (H_obj, H_con) can be used. The relaxation of Q-CHOP from needing the worst feasible state to needing any feasible state comes at the cost of having to solve a subproblem with an objective H_obj^★ that has quadratically larger norm, and can therefore be expected to require a larger penalty factor λ^★ and a larger quantum runtime T*. There may additionally be concerns if the energy E_★ has degeneracy d_★>>1, in which case a population transfer into the ground state of the final Q-CHOP Hamiltonian can be expected to be suppressed by a factor of 1/d_★.

Next, the treatment of inequality constraints in Q-CHOP and the SAA is discussed. An inequality constraint of the form D(x)≥0 can be converted into an equality constraint by the introduction of a slack variable: D(x)=n, where n≥0 is a new decision variable for the optimization problem. This equality constraint can then be rearranged and squared to arrive at the constraint

$\begin{array}{cc}{C}_D\left(x,n\right)={\left(D\left(x\right)-n\right)}^2=0,& \left(26\right)\end{array}$

at the cost of expanding the optimization domain as {0,1}^N→{0, 1}^N×{0, 1, 2, . . . }. The Hamiltonian enforcing this constraint takes the form

$\begin{array}{cc}{H_{\mathrm{con}}^D=(\sum _{x}D\left(x\right)❘x\rangle \langle x❘-{\hat{n}}_D)}^2,& \left(27\right)\end{array}$

where {circumflex over (n)}_D=Σ_n≥0 n|nn|_D is a number operator that can be identified with an auxiliary bosonic mode indexed by the constraint function D. In this sense, bosonic modes may be resources for enforcing inequality constraints.

In practice, the dimension of a slack variable can be truncated at D_max=max_x D(x) to encode a problem into a finite-dimensional Hilbert space. A discrete-variable encoding of the slack variable then requires └log₂(D_max+1)┐ qubits, though it may be desirable to encode the slack variable directly into a qudit or oscillator. If the maximum of D is not efficiently computable, D_max can be replaced by an upper bound on D(x). Expanding

$\begin{array}{cc}D\left(x\right)=\sum _{S\subset {\mathbb{Z}}_N}{c}_S^D\prod _{j\in S}{x}_j,& \left(28\right)\end{array}$

where each x_j∈{0, 1}, denotes the constant contribution to D(x) by c₀^D=c_{{ }}^D, and the set of remaining coefficients by C_D={c_S^D: S⊂_N, |S|>0}. A suitable upper bound for D(x) is then given by the sum of positive coefficients in C_D,

$\begin{array}{cc}{\tilde{D}}_{\max }=\mathrm{sum}\left({𝒞}_D\bigcap {\mathbb{Z}}_{+}\right)+c_0^D,& \left(29\right)\end{array}$

where _± denotes the set of positive (+) or negative (−) integers. If the coefficients in C_D have a nontrivial greatest common divisor gcd(C_D)>1, the dimensionality of the slack variable can be further reduced by pruning its allowed values to the set

$\begin{array}{cc}{𝕊}_D=\left(\gcd \left({𝒞}_D\right)\mathbb{Z}+c_0^D\right)\bigcap \left[0,{\tilde{D}}_{\max }\right].& \left(30\right)\end{array}$

For the SAA, the introduction of a slack variable may force a choice of initial state and driver Hamiltonian for the corresponding slack qudit. A natural generalization of the qubit-only case is to interpret the initial state |+^⊗N and transverse-field driver Hamiltonian

${H_{\mathrm{driver}}=-\frac{1}{2}{\sum }_j{X}_j\simeq -{\sum }_j❘+\rangle \langle +❘}_j$

respectively as a uniform prior and a sum of projectors onto the uniform prior for each decision variable. Here ≃ denotes equality up to an overall constant. This interpretation motivates initializing the slack qudit into a uniform superposition of its domain, and adding the (negative) projector onto this state to the driver Hamiltonian. In some examples, this generalization of initial state and driver Hamiltonian is used.

In some examples of Q-CHOP, there are two subtleties to address for the treatment of inequality constraints. First, once an initial worst feasible state |x_worst is identified, the slack variable for inequality constraint D may be initialized to |n_worst=|D(x_worst). Second, merely adding a slack qudit with the constraint in Eq. (27) to the Q-CHOP Hamiltonian fragments Hilbert space into uncoupled sectors indexed by the value of n=D(x). Initializing the slack variable in |n_worst thereby enforces, in effect, the much stronger constraint that D(x)=n_worst. To remedy this situation, a strategy of engineering an alternative adiabatic path that adds thoughtfully chosen mixing terms to the Q-CHOP Hamiltonian may be used.

In order to satisfy the constraint enforced by Eq. (27), added mixing terms may be used to couple the decision variables and the slack variables, changing them in a correlated manner. A Hamiltonian term that flips qubit j, for example, should simultaneously change the value of all slack variables associated with inequality constraints that involve variable x_j, so as to ensure that inequality constraints remain satisfied. However, engineering a strictly constraint-satisfying mixing Hamiltonian would involve complex multi-control logic that may generally be inefficient to implement. Accordingly, in some examples, slack variables are allowed to change to any value, and rely on the always-on constraint Hamiltonian of Q-CHOP to make constraint-violating terms off-resonant. For problems with inequality constraints, the objective Hamiltonian of Q-CHOP may be modified using a slack variable mixing operator S(θ). As a specific example of S(θ), the objective Hamiltonian of Q-CHOP may be modified as H_obj(θ)→H_obj(θ)×S(θ), with

$\begin{array}{cc}𝒮\left(\theta \right)=1+\sin \left(\theta \right)\prod _D❘{𝕊}_D❘P_{+}^D.& \left(31\right)\end{array}$

Here |_D| is the dimension of the slack qudit for constraint D and P₊^D is a projector onto is uniform superposition |+_D=Σ_n∈D|n/√{square root over (|_D|)}, such that |_D|P₊^D=|mn|_D. The factor of sin θ ensures that the added mixing terms vanish at θ=0 and π, which is necessary to preserve the initial and final ground states of the Q-CHOP Hamiltonian.

If the polynomial D(x) has degree d and coefficients with magnitude O(1), then |_D|=O(N^d). It is therefore worth noting that the slack operator in Eq. (31) grows with the product of the slack variable dimensions, Π_D|_D|. Though it is not surprising for the cost of solving an optimization problem to grow with the degree of its inequality constraints, it is highly desirable to have a cost that does not grow multiplicatively with the number of inequality constraints.

Numerical Experiments and Analysis

Described below are specific examples that compare the performance of Q-CHOP numerically using exact simulation and to SAA. A comprehensive performance comparison between Q-CHOP and the SAA requires a detailed, problem-aware assessment of the numerous choices that are made to instantiate simulation Hamiltonians, including choices of normalization factors, driver Hamiltonians (for the SAA), penalty factor, adiabatic schedules, and total evolution times.

To simulate Q-CHOP, an initial state is evolved with the Hamiltonian in Eq. (7) using a linear ramp, θ=πt/T. In the case of the SAA, the initial state |+^⊗N∝(|0+|1)^⊗N is evolved with the Hamiltonian

$\begin{array}{cc}{H}_{\mathrm{SAA}}\left(t\right)=-\left(1-\frac{t}{T}\right)S_x+\frac{t}{T}\left(H_{\mathrm{con}}+{\lambda }^{-1}{H}_{\mathrm{obj}}\right),& \left(32\right)\end{array}$

where

${S_x=\frac{1}{2}{\sum }_j{X}_j\simeq {\sum }_j❘+\rangle \langle +❘}_j$

is a transverse field normalized to have a spectral range equal to the qubit number N and the penalty factor is λ=N. For optimization problems with slack qudits, the qudits are initialized to a uniform superposition and projectors are added onto these initial states to S_x. In all the described experiments, the objective Hamiltonian is normalized to make the optimization problems scale-invariant. Specifically, H_obj is set to H_obj/(H_obj), where the normalization factor (H_obj) defined in terms of the nonzero coefficients of H_obj in Eq. (22) by

$\begin{array}{cc}{𝒩\left(H_{\mathrm{obj}}\right)}^2=\underset{\underset{S\ne \left\{\right\}}{c_S\ne 0}}{𝔼}{c}_S^2.& \left(33\right)\end{array}$

Here denotes an average, and the empty set is excluded to make the norm (H_obj) independent of constant shifts to H_obj. This normalization is equivalent to that of up to a constant factor of four. The constant factor is chosen so that (S_z)=1 for the uniform linear objective S_z≃−Σ_j|11|_j with spectral range N.

Unlike the objective Hamiltonians, constraint Hamiltonians typically have a natural expression as a sum of projectors, each of which has a spectral range of 1. This ensures that the constraint Hamiltonians have a spectral gap of at least 1, and typically equal (or otherwise close) to 1. Therefore, no additional preconditioning of the constraint Hamiltonians are performed unless specified otherwise. For problems with inequality constraints, all constraint coefficients are divided by their greatest common divisor by taking D(x)→D(x)/gcd(C_D∪c₀^D), prune the set of allowed slack variable values, and modify the SAA and Q-CHOP Hamiltonians to address slack variables as discussed above.

To highlight the features and results and various implementations, the following metrics are used to measure the performance of an optimization algorithm. The first is the in-constraint approximation ratio r, which is defined so that the best and worst feasible states have values of 1 and 0, respectively. The operator whose expectation is equal to the in-constraint approximation ration is given by

$\begin{array}{cc}r=P_{\mathrm{feas}}\left(\frac{H_{\mathrm{obj}}-E_{\mathrm{worst}}}{E_{\mathrm{best}}-E_{\mathrm{worst}}}\right)P_{\mathrm{feas}},& \left(34\right)\end{array}$

where P_feas is a projector onto the feasible-state manifold, and =|H_obj| for ∈{best, worst} are the objective energies of the best and worst feasible states.

TABLE 0.1
Problems simulated below and features
that are relevant for Q-CHOP.
		odd	no inequality
	problem	objective	constraints
	MIS	yes	yes
	DMDS	yes	yes
	knapsack	yes	no
	comb. auction	yes	no

The second quantity is the overlap between the state produced by the quantum algorithm and the state encoding the optimal solution to the optimization problem. The corresponding operator is the projector onto the optimal solutions of an optimization problem:

$\begin{array}{cc}{P}_{\mathrm{opt}}=\sum _{x:\langle x❘r❘x\rangle =1}❘x\rangle \langle x❘.& \left(35\right)\end{array}$

For optimization problems whose objective function values may lie anywhere on the real line, finding the exact solution is too strict of a criterion for success. In this case, the overlap with E-optimal states may be used, with the operator given by

$\begin{array}{cc}{P}_{ϵ}=\sum _{x:\langle x❘r❘x\rangle \ge 1-ϵ}❘x\rangle \langle x❘.& \left(36\right)\end{array}$

Table 0.1 summarizes the problems that are simulate in the numerical experiments described below, and some features that are relevant for Q-CHOP. All simulations were performed by numerically integrating equations of motion using the DOP853 method of scipy.integrate.solve_ivp, with both absolute and relative error tolerances set to 10⁻⁸.

Maximum Independent Set

The problem of finding a maximum independent set (MIS) is considered first. Given a graph G=(V, E), a set S⊂V is an independent set if none of its vertices are neighbors; that is, {v, w}∉E for all v, w∈S. A maximum independent set is an independent set containing the largest possible number of vertices. Candidate solutions S⊂V to MIS can be identified by bitstrings x={x_v: v∈V}∈{0, 1}^|V| for which x_v=1 if and only if v∈S. The identification of subsets is referred to as the membership encoding into |V|-bit strings. MIS can thus be defined by the following integer program:

$\begin{array}{cc}\mathrm{maximize}:\sum _{v\in V}{x}_v& \left(37\right)\end{array}$ $\mathrm{subject}\mathrm{to}:x_v{x}_w=0\mathrm{for}\mathrm{all}\left\{v,w\right\}\in E,$ $x\in {\left\{0,1\right\}}^{❘V❘}.$

The constraints of MIS are sometimes equivalently written as x_v+x_w≤1, but the formulation with equality constraints in Eq. (37) avoids the need for slack variables. The worst feasible solution to MIS is the empty set, identified by the all-0 bitstring x=(0, 0, . . . ). Many constrained optimization problems, including as maximum clique, minimum vertex cover, and maximum set packing, are straightforwardly reducible to MIS with no significant overhead.

Solving MIS with an adiabatic quantum algorithm was previously considered by Yu, Wilczek, and Wu. Specifically, solving MIS with Q-CHOP is equivalent to running the Yu, Wilzcek, and Wu algorithm with λ=−4 and T=4N². However, the previously known algorithm provides little insight into the choice λ. As a specific example, previous work did not discuss the fact that making |∂_tφ| too large (|λ| too small) makes the algorithm for MIS fail due to population leakage outside of the feasible subspace. In contrast, the formalism of Q-CHOP makes clear from Eq. (5) that one should choose |λ|˜N. The Q-CHOP formalism additionally clarifies that choosing ∂_tφ>0 (λ<0) makes these algorithms track the state of maximal energy within the feasible subspace, making them more vulnerable to mixing with non-feasible states.

FIG. 2 illustrates the summary of MIS simulation results 200. Expectation values of (top row) the approximation ratio r and (bottom row) the optimal-state projector P_opt in simulations of MIS with Q-CHOP 202 and the SAA 204. (Left-most panels) time-series data in simulations of 10 random graphs with N=10 vertices and quantum runtime T=2πN². The inset shows the in-constraint probability P_feas over time for the same simulations. (Middle panels) values at the end of simulations with the same random graphs, but different quantum runtimes T. (Right-most panels) values at the end of simulations with different qubit numbers (vertices). Lines in the right-most panels track average values for each qubit number N.

To evaluate Q-CHOP on MIS and to compare its performance to that of SAA, random (Erdös-Rényi) graphs G(N, p) with N vertices and edge-creation probability p=0.3 were considered. FIG. 2 summarizes simulation results from solving MIS with Q-CHOP and the SAA on random graphs for different qubit numbers N and quantum runtimes T. The primary takeaway from these results is that Q-CHOP not only outperforms the SAA in finding high-quality solutions on average (measured by the approximation ratio r), but also finds the optimal solution with significantly higher probability (measured by the projector P_opt). The probability of finding an optimal solution also grows much faster with respect to the quantum runtime T when solving MIS with Q-CHOP as opposed to the SAA. Moreover, for Q-CHOP the probability of finding an optimal solution remains roughly constant with respect to N when T˜N², whereas for the SAA this probability decays with N. For the SAA, a linear regression test rejects the null hypothesis of constant P_opt with respect to N (at T=2πN²) with a p-value of 6×10⁻⁴ (in favor of a negative slope), whereas for Q-CHOP, the null hypothesis is rejected with p-value of 0.01 in favor of a positive slope. While the utility of a linear regression test is limited for a quantity that is bounded on [0, 1], it illustrates the statistical significance of the trend observed in FIG. 2.

Directed Minimum Dominating Set

Another example is the problem of finding a directed minimum dominating set (DMDS). Given a directed graph G=(V, E), a vertex v∈V is dominated by a vertex u∈V if the directed edge (u, v)∈E. A dominating set of G is then a subset S⊂V for which every vertex of the graph is either contained in S, or dominated by a member of S. A minimum dominating set of a graph is a dominating set containing the fewest possible vertices. DMDS can be defined by the following integer program:

$\begin{array}{cc}\mathrm{minimize}:\sum _{v\in V}{x}_v& \left(38\right)\end{array}$ $\mathrm{subject}\mathrm{to}:{\stackrel{\_}{x}}_v\prod _{\left(u,v\right)\in E}{\stackrel{\_}{x}}_u=0\mathrm{for}\mathrm{all}v\in V,$ $x\in {\left\{0,1\right\}}^{❘V❘},$

where the binary variable x denotes membership in the dominating set and x_v=1−x_v for shorthand. The constraint for vertex v ensures that either v∈S (in which case x_v=0), or one of its neighbors u: (u, v)∈E is in S (in which case x_u=0). The worst feasible solution to DMDS is the all-inclusive set S=V, identified by the all-1 bitstring x=(1, 1, . . . ). The integer program formulation of DMDS in Eq. (38) superficially resembles that of MIS in Eq. (37), and indeed the problem of finding a minimum dominating set on an undirected graph can be reduced to solving MIS on a hypergraph. The definition of DMDS on a directed graph, however, prevents its straightforward reduction to a variant of MIS.

FIG. 3 illustrates a summary of DMDS simulation results with random directed graphs, presented in a format identical to FIG. 2. FIG. 3 provides a summary of using Q-CHOP 302 and the SAA 304 to solve DMDS on directed Erdös-Rényi graphs with edge-creation probability p=0.3, and existing edges are assigned a random orientation 300. Altogether, the results and conclusions from FIG. 3 are strikingly similar to those for MIS, likely due to the similarity of the underlying ensemble of graphs that define the MIS and DMDS problem instances. As with MIS, Q-CHOP yields high-quality solutions to DMDS on average, and has a higher probability of finding an optimal solution. The optimal state probability for Q-CHOP again grows faster with quantum runtime T, and does not decay with problem size N (Q-CHOP data in the is consistent with the null hypothesis of no decay, while for the SAA the null hypothesis is rejected with p-value 5×10⁻⁴), suggesting a potential advantage in scaling.

Knapsack

Q-CHOP is useful for problems beyond graph problems. For example, Q-CHOP may be used on the knapsack problem, a widely studied combinatorial optimization problem with applications in logistics and finance. The difficulty of solving the knapsack problem is the basis of knapsack cryptosystems, which are candidates for post-quantum cryptography, and closely-related lattice-based cryptosystems are under active consideration for the development of post-quantum cryptographic standards. A knapsack problem consists of a collection of items I to be placed in a “knapsack”. Each item i∈I has a weight w_i and value v_i, and the knapsack has weight capacity W. The goal is to choose a subset of items that maximizes value without exceeding the capacity of the knapsack. Again identifying subsets of I with the membership encoding, the knapsack problem can be defined by the following integer program:

$\begin{array}{cc}\mathrm{minimize}:\sum _{i\in I}{v}_i{x}_i& \left(39\right)\end{array}$ $\mathrm{subject}\mathrm{to}:\sum _{i\in I}{w}_i{x}_i\le W,$ $x\in {\left\{0,1\right\}}^{❘I❘}.$

The knapsack problem thereby introduces inequality constraints for Q-CHOP and the SAA, so simulating a knapsack problem instance requires a total Hilbert space dimension of 2^|I|(W+1). The worst feasible state of knapsack problem is the all-0 state.

The task of generating classically hard knapsack problem instances was previously studied by Jooken, Leyman, and De Causmaecker, with an implementation of their instance generator provided at: https://github.com/JorikJooken/knapsackProblemlnstances.

In addition to a number of items N and knapsack capacity W, an ensemble of random instances generated by the above method is defined by a number of groups g that determine item weight and value statistics, a fraction of items f that are assigned to the last group, an upper bound s on the value and weight of items in the last group, and a noise parameter ∈. As an experiment, to produce random knapsack problem instances with N items, the above generator implementation with (W, g, f, s, E)=(2N, ┌N/2┐, 1/N, N, 0.1) can be used. Any knapsack problem instances in which any item has weight greater than W are rejected.

FIG. 4 illustrates a summary of knapsack simulation results 400, presented in a format identical to FIG. 1. N as shown in FIG. 4 is the number of items in a knapsack instance. FIG. 5 illustrates histograms results for a specific example of the knapsack problem 500. In FIG. 5, the (Top panels) histograms are of final approximation ratios r and optimal state probabilities P_opt from knapsack simulations with N=8 items and quantum runtime T=2πN². The data in the top panels here is a coarse-grained slice of the data at time t=T in the left-most panels of FIG. 4. Darker bars 508 indicate that Q-CHOP 502 and SAA 504 values lie within the same histogram bin. In the (Bottom panels), the histograms of the performance differences (as measured by r and P_opt) between Q-CHOP and the SAA 506 in the same simulations as the top panels are shown. Q-CHOP outperforms the SAA in every simulated knapsack problem instance with N=8 items.

FIG. 4 shows the results of knapsack problem simulations with Q-CHOP 402 and the SAA 404, and additional data about results at N=8 and T=2πN² is provided in FIG. 5. The performance of Q-CHOP and the SAA for the knapsack problem is more variable than that for MIS and DMDS, which can be attributed to a richer ensemble of problem instances. The SAA generally fails to produce final states that lie within the feasible-state manifold, and in some cases essentially fails to find an optimal state altogether. Interestingly both Q-CHOP and the SAA exhibit plateau behavior in the growth of r and P_opt with progress t/T throughout the solution of individual problem instances. Even so, Q-CHOP achieves reliably high approximation ratios, and finds optimal states with higher probability than the SAA in every simulated instance (see FIG. 5). Q-CHOP also achieves significantly higher approximation ratios and optimal-state probabilities with shorter quantum runtimes T than the SAA, though there is high variability in the growth of P_opt with T. The performance gap between Q-CHOP and the SAA persists for different problem sizes N.

Combinatorial Auction

As another example, the problem of optimizing a combinatorial auction (CA), which is commonly used mechanism for market-based resource allocation and procurement, is considered. A combinatorial auction consists of a collection I of items to be sold, and a set B of bids that offer to purchase baskets of items. Each item i∈I has multiplicity m_i. Each bid b=(P_b, q_b)∈B offers to pay p_b for the basket of items specified by the vector q_b=(q_bi: i∈I), where q_bi is the quantity of item i in the basket requested by bid b. The goal of a combinatorial auction is for the seller to accept a subset of the bids that maximizes payments without exceeding the available inventory of items. Identifying the subset of accepted bids by the membership encoding, a combinatorial can be defined by the following integer program:

$\begin{array}{cc}\mathrm{maximize}:\sum _{b\in B}{p}_b{x}_b& \left(40\right)\end{array}$ $\mathrm{subject}\mathrm{to}:\sum _{b\in B}{q}_{\mathrm{bi}}{x}_b\le m_i\mathrm{for}\mathrm{all}i\in I,$ $x\in {\left\{0,1\right\}}^{❘B❘}.$

Due to the inequality constraints and associated slack variables, a simulating a combinatorial auction problem instance requires a total Hilbert space dimension of 2^|I|Π_i∈I(m_i+1). The worst feasible solution of the combinatorial auction is the all-0 state.

As can be seen by comparing Eq. (39) to Eq. (40), the knapsack problem can be interpreted as a combinatorial auction consisting of one item with multiplicity W. A combinatorial auction of unique items (in which all multiplicities m_i=1) can be also be reduced to an instance of weighted MIS on a graph (V, E) in which nodes are identified with bids, V=B, and weighted by the payments p_b in the objective function, as in Eq. (40). The edges E are obtained by inserting, for each item i∈I, a complete graph on the bids that request item i. In other words, E=∪_i∈I E_i with E_i={{a, b}: q_ai=q_bi=1}.

FIG. 6 illustrations a summary of combinatorial auction simulation results 600, presented in a format identical to FIG. 1. In FIG. 6, N is the number of bids in the combinatorial auction. FIG. 7 illustrates histograms of combinatorial simulation results for N=8 bids and quantum runtime T=2πN² 700, presented in a format identical to Figure FIG. 5. Darker bars 708 indicate that Q-CHOP 702 and SAA 704 values lie within the same histogram bin. In the (Bottom panels), the histograms of the performance differences (as measured by r and P_opt) between Q-CHOP and the SAA 706 in the same simulations as the top panels are shown. Q-CHOP outperforms the SAA in every simulated combinatorial auction problem instance with N=8 bids.

To generate combinatorial auction problem instances, the implementation of the Combinatorial Auction Test Suite provided by Gasse, Chelelat, Ferroni, Charlin, and Lodi was used. Instances of combinatorial auctions with |I|=3 items and N bids for varying N were generated. All other parameters defining an ensemble of combinatorial auctions are set to the defaults as in the provided implementation. These defaults set all multiplicities m_i=1, so the combinatorial auction instances are equivalent to instances of weighted MIS drawn from an appropriately defined ensemble of weighted graphs. For the sake of testing Q-CHOP and the SAA with inequality constraints, however, the formulation of the combinatorial auction in Eq. (40) were kept, which is simulated with auxiliary slack variables.

FIG. 6 and FIG. 7 summarize the combinatorial auction simulation results. Similarly to the results for the knapsack problem, both Q-CHOP 602, 702 and SAA 604, 704 performance is variable across different problem instances, though for the SAA the growth of r and P_opt tapers off with progress t/T rather than hitting a plateau. Despite exhibiting plateau behavior, Q-CHOP still achieves more reliably high approximation ratios, higher optimal-state probabilities, and a faster growth of these performance metrics with quantum runtime T. The performance gap between Q-CHOP and the SAA persists for different problem sizes N.

FIG. 8 is a flowchart of an example process for processing a constrained optimization problem 800. In some implementations, one or more process blocks of FIG. 8 may be performed by a quantum processor, a classical processor, or a control system.

At 810 a quantum processor applies coupling and transformation operations to a plurality of quantum states according to a Hamiltonian specification. The quantum processor may include a plurality of quantum processing elements associated with respective quantum states. At 820 a digital computer receives information based at least in part on measurements of one or more quantum states associated with respective quantum processing elements of the quantum processor. At 830, the digital computer provides information for preparing one or more quantum states associated with respective quantum processing elements of the quantum processor based at least in part on the information. At 830, a control module applies the applied coupling and transformation operations based on interaction with the digital computer for processing the constrained optimization problem. At 840, the processing includes preparing quantum states associated with a plurality of the quantum processing elements characterized by a summation of a constraint Hamiltonian representing a constraint of the constrained optimization problem and an objective Hamiltonian representing an objective function of the constrained optimization problem. At 850, the control module evolves a time-dependent Hamiltonian according to an evolution that includes forming a sum of a first term having the constraint Hamiltonian and a second term. The second term includes a product of (A) a time-dependent scalar function and (B) a time-dependent operator that is initially equal to the objective Hamiltonian and is evolved into a negative of the objective Hamiltonian.

Although FIG. 8 shows example blocks of process 800, in some implementations, process 800 may include additional blocks, fewer blocks, different blocks, or differently arranged blocks than those depicted in FIG. 8. Additionally, or alternatively, two or more of the blocks of process 800 may be performed in parallel.

It will thus be seen that the objects set forth above, among those made apparent from the preceding description, are efficiently attained and, because certain changes may be made in carrying out the above method and in the construction(s) set forth without departing from the spirit and scope of the invention, it is intended that all matter contained in the above description and shown in the accompanying drawings shall be interpreted as illustrative and not in a limiting sense.

Although the foregoing examples have been described in some detail for purposes of clarity of understanding, the invention is not limited to the details provided. There are many alternative ways of implementing the invention. The disclosed examples are illustrative and not restrictive.

Claims

1. An apparatus for processing a constrained optimization problem, the apparatus comprising: a quantum processor comprising a plurality of quantum processing elements associated with respective quantum states, and configured to apply coupling and transformation operations to a plurality of the quantum states according to a Hamiltonian specification;a digital computer comprising at least one central processing unit, the digital computer configured to: receive information based at least in part on measurements of one or more quantum states associated with respective quantum processing elements of the quantum processor; andprovide information for preparing one or more quantum states associated with respective quantum processing elements of the quantum processor based at least in part on the received information; anda control module configured to control the applied coupling and transformation operations based on interaction with the digital computer for processing the constrained optimization problem, the processing comprising: preparing quantum states associated with a plurality of the quantum processing elements characterized by a summation of a constraint Hamiltonian representing a constraint of the constrained optimization problem and an objective Hamiltonian representing an objective function of the constrained optimization problem; andoperating the control module to evolve a time-dependent Hamiltonian according to an evolution that includes forming a sum of a first term comprising the constraint Hamiltonian and a second term, where the second term comprises a product of (A) a time-dependent scalar function and (B) a time-dependent operator that is initially equal to the objective Hamiltonian and is evolved into a negative of the objective Hamiltonian.

2. The apparatus of claim 1, wherein the objective Hamiltonian comprises an even term, wherein the time-dependent operator locally rotates the even term.

3. The apparatus of claim 1, wherein the objective Hamiltonian comprises a plurality of even terms, wherein for each of the plurality of even terms the time-dependent operator comprises a second sum of local rotations of components of the even term divided by a number of components.

4. The apparatus of claim 3, wherein the objective Hamiltonian comprises a plurality of odd terms, wherein for each of the plurality of odd terms the time-dependent operator comprises a sum of local rotations of components of the odd term divided by a number of components.

5. The apparatus of claim 3, wherein the objective Hamiltonian comprises a plurality of odd terms, further comprising: partitioning the objective Hamiltonian into an odd objective Hamiltonian and an even objective Hamiltonian; andglobally rotating the odd objective Hamiltonian.

6. The apparatus of claim 1, wherein the objective Hamiltonian comprises a sum of weighted terms that represent the constrained optimization problem, and at least two of the weighted terms have different weights from each other.

7. The apparatus of claim 6, wherein the weighted terms correspond to respective vertices, edges, or hyperedges of a graph or hypergraph.

8. The apparatus of claim 1, wherein the constrained optimization problem comprises a weighted constrained optimization problem.

9. The apparatus of claim 8, wherein the weighted constrained optimization problem comprises a problem selected from the group consisting of: weighted maximum independent set, weighted maximal clique, weighted minimum vertex cover, weighted maximum set packing, weighted minimum dominating set, weighted minimum set cover, and weighted minimum dominating set on a directed graph.

10. The apparatus of claim 1, wherein the constrained optimization problem comprises an inequality constraint and wherein the objective Hamiltonian is modified by a slack variable mixing operator.

11. The apparatus of claim 10, the slack variable mixing operator is an identity plus a term that mixes slack variable amongst themselves.

12. The apparatus of claim 1, wherein the constrained optimization problem comprises a knapsack problem.

13. The apparatus of claim 1, wherein the constrained optimization problem comprises a combinatorial auction problem.

14. A method for processing a constrained optimization problem, the method comprising: applying, using a quantum processor, coupling and transformation operations to a plurality of quantum states according to a Hamiltonian specification, wherein the quantum processor comprises a plurality of quantum processing elements associated with respective quantum states;receiving, at a digital computer, information based at least in part on measurements of one or more quantum states associated with respective quantum processing elements of the quantum processor;providing, from the digital computer, information for preparing one or more quantum states associated with respective quantum processing elements of the quantum processor based at least in part on the information; andapplying, from a control module, the applied coupling and transformation operations based on interaction with the digital computer for processing the constrained optimization problem, the processing comprising: preparing quantum states associated with a plurality of the quantum processing elements characterized by a summation of a constraint Hamiltonian representing a constraint of the constrained optimization problem and an objective Hamiltonian representing an objective function of the constrained optimization problem; andoperating the control module to evolve a time-dependent Hamiltonian according to an evolution that includes forming a sum of a first term comprising the constraint Hamiltonian and a second term, where the second term comprises a product of (A) a time-dependent scalar function and (B) a time-dependent operator that is initially equal to the objective Hamiltonian and is evolved into a negative of the objective Hamiltonian.

15. The method of claim 14, wherein the objective Hamiltonian comprises an even term, wherein the time-dependent operator locally rotates the even term.

16. The method of claim 14, wherein the objective Hamiltonian comprises a plurality of even terms, and wherein for each of the plurality of even terms the time-dependent operator comprises a sum of local rotations of components of the even term divided by a number of components.

17. The method of claim 16, wherein the objective Hamiltonian comprises a plurality of odd terms, wherein for each of the plurality of odd terms the time-dependent operator comprises a second sum of local rotations of components of the odd term divided by a number of components.

18. The method of claim 16, wherein the objective Hamiltonian comprises a plurality of odd terms, further comprising: partitioning the objective Hamiltonian into an odd objective Hamiltonian and an even objective Hamiltonian; andglobally rotating the odd objective Hamiltonian.

19. The method of claim 14, wherein the constrained optimization problem comprises an inequality constraint and wherein the objective Hamiltonian is modified by a slack variable mixing operator.

20. The method of claim 19, the slack variable mixing operator is an identity plus a term that mixes slack variable amongst themselves.