Grover Adaptive Search for Constrained Polynomial Binary Optimization

In this paper we discuss Grover Adaptive Search (GAS) for Constrained Polynomial Binary Optimization (CPBO) problems, and in particular, Quadratic Unconstrained Binary Optimization (QUBO) problems, as a special case. GAS can provide a quadratic speed-up for combinatorial optimization problems compared to brute force search. However, this requires the development of efficient oracles to represent problems and flag states that satisfy certain search criteria. In general, this can be achieved using quantum arithmetic, however, this is expensive in terms of Toffoli gates as well as required ancilla qubits, which can be prohibitive in the near-term. Within this work, we develop a way to construct efficient oracles to solve CPBO problems using GAS algorithms. We demonstrate this approach and the potential speed-up for the portfolio optimization problem, i.e. a QUBO, using simulation and experimental results obtained on real quantum hardware. However, our approach applies to higher-degree polynomial objective functions as well as constrained optimization problems.


Introduction
Using the laws of quantum mechanics, quantum computers offer novel solutions for resource-intensive problems.
A commonly-studied class of combinatorial optimization problems are Quantum Unconstrained Binary Optimization (QUBO) problems, with applications in resource allocation, finance, machine learning, and partitioning. There are multiple approaches to solve QUBO problems on a quantum computer, discussed in the following paragraphs.
First, quantum annealing [15,16] is a metaheuristic for adiabatic quantum computers. Similar to simulate annealing, it can be used to approximate optimal solutions of QUBO problems.
Further, there exist variational quantum optimization heuristics, such as Variational Quantum Eigensolver (VQE) and Quantum Approximate Optimization Algorithm (QAOA) [10][11][12][13]. VQE and QAOA are heuristics designed for near-term, noisy quantum computers without performance guarantees. However, for QAOA, it is known that in the infinite depth limit, the algorithm recovers adiabatic evolution and would converge to the optimal solution.
Last, there are Grover-based [1] optimization algorithms, such as Grover Adaptive Search (GAS) [17][18][19]. GAS iteratively applies Grover Search to find the optimum value of an objective function, by using the best-known value as a threshold to flag all values smaller than the threshold in order to find a better solution. The algorithmic framework comes with a quadratic speed-up, however it likely requires an error-corrected fault-tolerant quantum computer due to the depth of the resulting circuits. One of the challenges inherent in GAS is the creation of efficient oracles.
In this paper, we provide a framework for automatically generating efficient oracles for solving Constrained Polynomial Binary Optimization (CPBO)a generalization of QUBO-with GAS. The objective function and constraints need to be efficiently encoded, for which we use a Quantum Dictionary [20], a pattern for representing key-value pairs as entangled quantum registers, that turns out to be efficient for polynomial functions -in particular for quadratics representing QUBO problems. The approach relies on the addition of classical numbers to a quantum register in superposition, conditioned on the state of another quantum register. It is similar to the method used in Quantum Fourier Transform (QFT) adders [21]. Given a boolean polynomial, the coefficient of each monomial is added to the value register conditioned on the qubits in the key register corresponding to the variables present in the monomial.
We test our algorithm on the portfolio optimization problem [22,23]. Multiple variants of this problem have been studied in the quantum optimization literature, ranging from convex continuous formulations [24] to QUBOs [13,25,26]. Here we investigate a QUBO formulation as well as a formulation with an inequality budget constraint, the latter not being compatible with other approaches like quantum annealing, VQE, or QAOA-as those approaches can only handle linear equality constraints through quadratic penalty terms.
The remainder of this paper is organized as follows. Sec. 2 introduces GAS in general. Sec. 3 introduces QUBO problems, and shows how we can efficiently generate oracles to solve them using GAS algorithms, as well as how this approach extends to more general CPBO problems. In Sec. 4, we apply the developed technique to a concrete test case -portfolio optimization -and demonstrate it via simulation using Qiskit [27]. Sec. 5 concludes this paper and discusses possible directions of future research.

Grover Adaptive Search
Optimization problems are often solved by sequential approximation methods. In many cases, such methods are the only choice, but they may be computationally more efficient even when a solution to a problem can be expressed in a closed form. GAS works in a similar way, as it repeatedly uses Grover Search to randomly sample from all solutions better than the current one.
Grover Search is often described as a search algorithm, because it was initially formulated in the context of finding a single state of interest in a superposition of n-qubit quantum states. The algorithm has been generalized to the case of multiple states of interest, in which case it is better interpreted as a sampling algorithm. It amplifies the amplitudes of the states of interest within a larger search space, thus, increasing the probability of measuring one of the target states.
Grover Search -the core element of GAS -needs three ingredients: 2. An oracle operator O, that recognizes the states of interest and multiplies their amplitudes by -1. For instance, suppose I ⊂ {0, . . . , 2 n −1} denotes the set of target states and A = H ⊗n , then 3. The Grover diffusion operator D, that multiplies the amplitude of the |0 n state (or, equivalently, all states except |0 n ) by -1. if y < y i then 7 x i+1 = x, y i+1 = y, and k = 1 8 else 9 x i+1 = x i , y i+1 = y i , and k = λk 10 i = i + 1; 11 until a termination condition is met; The diffusion operator has the net effect of inverting all amplitudes in the quantum state about their mean. This causes all the amplitudes of the states of interest to be magnified, while the amplitudes of all other states are decreased. More precisely, applying the Grover operator G = ADA † O the right number of times to state A |0 n -i.e. evaluating G r A |0 n for an integer r ≥ 0 -will maximally amplify the amplitudes of the states of interest. The optimal number of applications r depends on the number N = 2 n of all states and the number s of states of interest, and is equal to π 4 N s . This implies a probability of sampling a target state of at least 1/2, which corresponds to a quadratic speed-up compared to classical search. Since s is in general unknown, we can either use Quantum Counting algorithms [28][29][30][31] to find s, or apply a randomized strategy.
The latter is the essence of [32], where an algorithm for applying Grover Search for unknown s is presented. This was then used to create a minimumfinding algorithm [17], which we refer to as GAS. In the following we outline GAS, which is formally given in Alg. 1.
Consider a function f : X → R for n binary variables, where for ease of presentation assume X = {0, 1} n , for which we are interested in finding min x∈X f (x). The main idea of GAS is to construct A y and O y for a given threshold y such that they flag all states x ∈ X satisfying f (x) < y, such that we can use Grover Search to find a solutionx with a function value better than y. Then we set y = f (x) and repeat until some formal termination criteria is met, e.g. based on the number of iterations, time, or progress in y.
While implementations of GAS vary around the specific use case [19,33], the general framework still loosely follows the steps described in [17]. In the following, we will show how operator A and oracle O can be efficiently constructed for QUBO as well as CPBO problems.

QUBO and CPBO Oracles
A QUBO problem with n binary variables is specified by a quadratic operator represented by a matrix Q ∈ R n×n , vector b ∈ R n , and constant c ∈ R, defined as or more compactly as min x∈{0,1} n ( In the following, we show how to efficiently construct GAS oracles for QUBO problems. We will construct A y such that it prepares a n-qubit input register to represent the equal superposition of all |x n and a m-qubit output register to (approximately) represent the corresponding |f (x) − y m . Then, the oracle O y should flag the states with a negative value in the output register. Note that in the implementation discussed, the oracle operator is actually independent of y, but this is not a requirement. For clarity, we will refer to the oracle as O when the oracle is independent of y. More formally, we show how to construct the oracles such that where |x is the binary encoding of the integer x. Furthermore, we will show how the developed technique can be used to extend GAS to higher-degree polynomials of binary variables, as well as to constrained optimization.

Construction of operator A
To construct A, we will use a Quantum Dictionary, as introduced in [20], and summarize the construction in the following subsections.

Encoding a Single Integer Value
Given an m-qubit register and an angle θ ∈ [−π, π), we wish to prepare a quantum state whose state vector represents a "periodic signal" equivalent to a geometric sequence of length 2 m . This can be implemented using a unitary operator defined by Fig. 1.
The simplest implementation of U G (θ) uses the phase gate R(θ) that, when applied to a qubit, rotates the phase of the amplitudes of the states having 1 in the position corresponding to that qubit. In Qiskit, this gate is the U 1 (θ) operator [27]. The circuit for U G (θ) is shown in Fig. 2, and consists of applying the gate R(2 i θ) to the qubit m − 1 − i in the m-qubit register prepared in the state of equal superposition in equation (1).
, applied to an m-qubit register. R denotes the phase gate, which rotates the amplitudes of states having 1 in the position corresponding to the qubit it was applied to.
In [20] an alternative way to implement U G was introduced by applying controlled R y rotations to an ancillary register containing the encoding of an eigenstate of R y , as explained in Appendix A.
Note that QFT applied to a register containing the binary encoding of a non-negative integer also creates a geometric sequence of amplitudes in that register. U G can be seen as a shortcut for the QFT when the encoded numbers are known classically, as we avoid multi-qubit interactions.
Given an integer −2 m−1 ≤ k < 2 m−1 , if we apply U G ( 2π 2 m k), followed by the inverse QFT to an m-qubit register prepared in the state of equal superposition in equation (1), we end up with k (mod 2 m ) being encoded in the register, as shown in Fig. 3. This representation is called the binary Two's Complement of k, which just adds 2 m to negative values k, similar to the way we can represent negative angles with their complement, e.g. equating −π/4 with 7π/4. The reason this representation occurs naturally in this context is due to the fact that rotation composition behaves like modular arithmetic. The same method can be used to encode non-integers, as discussed in Appendix B.

Encoding a Superposition of Polynomial Values
Next we will discuss the application of U G (θ) to a register |z m representing the values of a function, controlled on a register |x n representing the inputs of the same function. In general, the application of a unitary operator U to |z m controlled by a set J ⊆ {0, . . . , n − 1} of qubits of |x n can be expressed as and an example is shown in Fig. 4. The general form of a polynomial of n variables is: Each subset J ⊆ {0, . . . , n − 1} has a corresponding monomial j∈J x j . The free term is represented by a ∅ .

Construction of Operator A
Now we are ready to define the operator A, as shown in Fig. 5. It consists of applying a controlled geometric sequence transformation C J (U G ( 2π 2 m a J )) for each subset J ⊆ {0, . . . , n − 1} with a non-zero coefficient a J , followed by a single application of the inverse QFT. Note that QUBO polynomials have only monomials of degree less than or equal to 2, i.e. |J| ≤ 2, so we only need to control on single and pairs of qubits. An example of encoding a single monomial is shown in Fig. 6. Figure 6: An example of encoding the monomial 2x1x3, with input register |x 4 and output register |z 3 . R denotes the phase gate, which rotates the amplitudes of states having 1 in the position corresponding to the qubit it was applied to. The operator A prepares a state where the |x n register holds all 2 n inputs in equal superposition, entangled with the corresponding values P (x) encoded in the |z m register:

Gate Count
where we are assuming that an appropriate size m for the value register is known. The desired A y operator is obtained by adding the −y constant to the free term of the original polynomial. This construction works for polynomials of arbitrary degree.
A detailed summary of the number of gates required to construct A for a QUBO is given in Table 1. In general, the number of gates scales as the number of monomials in the polynomial to be loaded. This matches the scaling of the description of the problem under the assumption that the input weights and output are represented using roughly the same number of bits, thus, the asymptotic scaling of our approach is optimal.
An alternative approach to construct A would be to use quantum arithmetic. However, even uncontrolled in-place addition of a classically given m-bit number to a m-bit quantum register requires 2m qubits and 2m − 1 Toffoli gates [35][36][37]. This then would have to be done O(n 2 ) times in a 2-controlled way, which leads to significantly larger circuits and m more qubits than Figure 7: A constrained optimization circuit, which encodes the functions corresponding to a given set of constraints into an m-qubit register, and flips a related indicator qubit for each condition that is satisfied. In this example, given a function f : X → R of n binary variables, the global indicator qubit (shown at the bottom of the circuit) is set to |1 if and only if f (x) + y-where y is the threshold parameter from GAS-and the cost of the input (denoted C(x)) are both less than 0.
Note that the encoding method is the same as described in Sec. 3.1, meaning that the binary representation of the values is in Two's Complement, and thus we only need to control on a single qubit (the most significant bit) to determine if a value is negative. There is a trade-off between reusing the m qubits and reversing the computation (shown here), or adding additional value qubits for each constraint, allowing one to skip the uncompute. In either case, we can replace the oracle that flags "good states" by a multi-controlled Z gate, with a number of controls equal to one more than the number of constraints. Here, we add one qubit and use a Toffoli gate to get back to the setting introduced before.
required by our approach, and also requires explicit encoding of negative integers.

Construction of oracle O
At each step of the algorithm, we are adding a constant to the polynomial, and searching for remaining negative values. This means the oracle just needs to recognize negative integers. Since values are represented in Two's Complement, where the most significant (left-most) bit designates the sign of the number, a single qubit in the value register can be used to recognize negative integers. The typical oracle that multiplies target amplitudes by −1 can be applied. Note that the oracle stays unchanged between iterations, because we add a constant to the polynomial, which may lead to overflow in the value register. In order to avoid that, we may need to increase the number of qubits in the value register by 1. Alternatively, threshold-based oracles (which are potentially more expensive) can be used, that will reduce the search space by filtering the numbers above a given threshold.

Constrained Optimization
It is common for optimization problems to impose additional constraints-e.g. the total number or cost of assets in a portfolio may be subjected to an upper bound. Such constraints translate into further reductions of the search space based on the key register. For example, the number of assets in a portfolio corresponds to the number of 1s in the binary representation of the input of the objective function, called the Hamming weight.
It is straightforward to use the GAS oracles introduced in Sec. 3.1 and 3.2 to take into account polynomial equality and inequality constraints on the key register. Therefore, we can add additional registers to evaluate other polynomials. Whether an inequality constraint is satisfied or not can again be mapped to the sign qubit by applying an appropriate shift to the polynomial. Equality constraints are a bit more expensive, as they require the detection of a particular state, which essentially has the same complexity as the Grover diffusion operator D.
This leads to a set of qubits flagging target states: one qubit identifying the states that correspond to objective values below the current threshold, and one qubit for each constraint. Applying a logical ANDoperation to all of them essentially acts as an intersection of the individually flagged states and allows to construct oracles for CPBO. An example is shown in Fig. 7, which we demonstrate in Sec. 4.1.

Test Cases
In the remainder of this paper we demonstrate the proposed technique on the portfolio optimization problem, and then show a simple example on quantum hardware.

Portfolio Optimization
Suppose an investment universe consisting of n assets, denoted by i = 1, . . . , n, their corresponding expected returns µ ∈ R n and covariance matrix Σ ∈ R n×n . Furthermore, we consider a given risk factor q ≥ 0, which determines the considered risk appetite. The resulting objective function is In other words, we want to minimize the weighted variance minus the expected portfolio return. Setting q = 0 implies a risk neutral investor, while increasing q increases its risk averseness.
In the presented form, portfolio optimization is a QUBO problem. We can extend it by imposing a budget constraint of the form where B ∈ {0, . . . , n} denotes the number of assets to be chosen. In general, equality constraints can be recast as penalty terms and added to the objective (since we minimize), where λ > 0 is a large number to enforce the constraint to be satisfied. This results in a quadratic term that will again lead to a QUBO problem. However, with the methodology introduced in Sec. 3.3, we can also model more complex constraints, e.g. a budget inequality constraint of the form where c i ∈ R denote the asset prices, which does usually make more sense in practice than (10). In the following we consider two examples, one with an objection function that we want to minimize, and then the same problem with added constraints.

Finding a Minimum Value
Consider a portfolio of n = 3 assets, risk factor q = 0.5, and returns described by: which leads to the formulation The objective function with added constant −y (where y is the current threshold) has an associated A y operator and the oracle O recognizes negative values, as introduced in Sec. 3. To perform the experiment we need 7 qubits split into two registers, n = 3 input qubits and 4 output qubits. While we only need 3 qubits in the output register, we add 1 extra to accommodate for the threshold shift. We set the initial threshold y 1 = 0, and stop searching if an improvement has not been found in three consecutive iterations of the algorithm.
For each iteration of GAS, we determine the number of applications of the Grover iterate as defined in Alg. 1. We apply A y for the current threshold, and Iteration 2 (y = 5, r = 1) Figure 8: The output probabilities of GAS for three iterations, searching for a minimum value of a function (Eq. 14), each with thresholds y and r applications of the Grover operator.
then apply the Grover iterate G y = A y DA † y O for the predetermined number of applications. If the measured value is less than y, we update the threshold. This process repeats until we have seen no improvement for three consecutive iterations.
Classically, we can determine the original minimum value by keeping track of the total shift, or by calculating the value of the objective function for the minimum key. The results of this simulated experiment are shown in Fig. 8.

Additional Constraints
We can impose a budget inequality constraint of the form discussed in Eq. 12 to the previous problem, where B < 2 and the price of each asset is 1. As shown in Sec. 3.3, we can implement this constraint by adding an additional register to the existing quantum circuit, and then encode the Hamming weight of the binary representation of the keys in that register. Note that we only need to control on the most significant qubit of the constraint register to determine if B < 2. Otherwise, the procedure for applying GAS is the same as in the original problem. The results of this simulated experiment are shown in Fig. 9.

Trials on Real Hardware
The experiments discussed in this section were run on the IBM ibmq_toronto device, with Quantum Volume 32 [38].
The configuration details for ibmq_toronto at the time of our experiments is given in Appendix C. Readout error-mitigation techniques were applied to the results of each circuit [27,39]. Iteration 2 (y = 1, r = 1) Figure 9: The output probabilities of three iterations of GAS, with respective thresholds y and r applications of the Grover operator. We want to find a minimum value of Eq. 14, with the additional constraint that the binary representation of the corresponding key has a Hamming weight less than 2.
Let us consider a simple example of a polynomial minimization which can be run on current quantum hardware: min x∈{0,1} 2 (−2 + x 1 + x 2 ). (15) As in the previous subsection, we set the initial threshold y 1 = 0, and stop searching after three iterations with no improvement. Note that due to the probabilistic nature of the algorithm there is a nonzero probability that an invalid key-value pair will be measured. In addition, the noise inherent in the present era of quantum hardware further impacts the results and increases the probability of wrong results. We repeat the computation several times, and take the outcome with the maximum probability as the measured result. As long as the noise is not too strong, this still achieves a good key-value pair.
The results of the experiment are shown in Fig. 10. Note that in this example, the valid measurement outcomes are 0 → −2, 1 → −1, 2 → −1, and 3 → 0. In the first iteration the four outcomes are shown in an approximately-equal superposition, and sampled randomly (we do not apply the Grover iterate). We measure 1 → −1, and thus we update the threshold to y 2 = −1 and the number of iterations to r = 1.
In the second iteration of Grover Adaptive Search the results of the amplification are shown, where 0 → −2 (the minimum) has the highest probability.  Figure 10: The output probabilities of GAS for two iterations run on real quantum hardware, with respective thresholds y and r applications of the Grover operator. The valid measurement outcomes are shown in blue, while the invalid measurement outcomes are shown in grey.

Conclusion
In this paper we introduced an efficient way to implement the oracles required for solving Constrained Polynomial Binary Optimization problems using Grover Adaptive Search. This problem class is very general and contains for instance QUBO problems. Our approach significantly reduces the number of gates required compared to standard quantum arithmetic approaches, i.e. it lowers the requirements to apply GAS on real quantum hardware for practically relevant problems. We demonstrated our algorithm on the portfolio optimization problem, i.e. a QUBO, where we could reliably find the optimal solution, and on real quantum hardware. Within this manuscript we focused mainly on problems with integer coefficients. The handling of non-integers is discussed in Appendix B.

Code Availability
All code associated with this publication can be found in the IBM Qiskit Optimization module, currently hosted at https://github.com/Qiskit/qiskitoptimization.

Acknowledgments
This material is for informational purposes only and is not the product of JPMorgan Chase & Co.'s Research Department. This material is not intended as research, a recommendation, advice, offer or solicitation for the purchase or sale of any financial product or service, and is not a research report and is not intended as such. This material is not intended to represent any position or opinion of JPMorgan Chase & Figure 11: An example of an integer encoding (left), accompanied by two examples of a non-integer encoding (middle and right). Note that in the integer case only one outcome is possible, whereas the encoding of a non-integer leads to a superposition of approximations.
Co. JPMorgan Chase & Co. disclaims any responsibility or liability whatsoever for the quality, accuracy or completeness of the information herein, and for any reliance on, or use of this material in any way. ©2021 JPMorgan Chase & Co.
IBM, the IBM logo, and ibm.com are trademarks of International Business Machines Corp., registered in many jurisdictions worldwide. Other product and service names might be trademarks of IBM or other companies. The current list of IBM trademarks is available at https://www.ibm.com/legal/copytrade.

A Quantum Dictionary
The Quantum Dictionary was introduced in [20] as a quantum computing pattern for encoding functions, in particular polynomials, into a quantum state using geometric sequences. The paper shows how quantum algorithms like search and counting applied to a quantum dictionary allow to solve combinatorial optimization and QUBO problems more efficiently than using classical methods. Encoding a geometric sequence can be done using the phase gate as we explained earlier, but it is worth mentioning an alternative method described in [20], which uses the R y family of gates. The R y gate is applied to an ancillary register containing one of its eigenstates prepared by an operator E(R y ), conditioned on both the key and value registers. The rotation angle θ is different for each application, representing a number that contributes to the values corresponding to keys that have the conditioned key as a subset. In particular, when encoding a polynomial, a rotation will be applied for each of its coefficients.
In [20] we used the R y operator, with the eigenstate 1/ √ 2(i |0 + |1 ), independent of the rotation angle. This state can be prepared by the circuit in Fig. 13. The corresponding eigenvalue of R y (2θ) is e iθ .

B Handling Non-Integers
When handling non-integers, we have two choices. The first is to approximate them with fractions with a common denominator and encode the numerator into quantum registers before performing computations, cf. Appendix B.1. The second, cf. Appendix B.2, is to phase-encode real numbers and let the inverse QFT convert the result of the computation into a superposition of approximations.

B.1 Approximating Real Coefficients by Fractions
If we relax the assumption that all coefficients are integers, we can approximate non-integers by dividing all values in µ and Σ by the largest (scaling the range of the coefficients to [−1, 1)), and approximating each value k by a fraction k 2 m with −2 m−1 ≤ k < 2 m−1 , where m is the number of value qubits. As an example, suppose n = 3, m = 5, q = 0.5, and which results in the optimization problem min x∈{0,1} 3 −(−16x 1 + 5x 2 + 10x 3 ). (20)

B.2 Encoding Real Coefficients as Fejér Distributions
Recall the unitary operator U G (θ) from Fig. 1, where θ ∈ [−π, π). As discussed in Sec. 3.1, the process for encoding an integer −2 m−1 ≤ k < 2 m−1 is to apply U G ( 2π 2 m k) to an m-qubit register in equal superposition, followed by a single application of the inverse Quantum Fourier Transform. For the general case of a real number, shown in Fig. 14, where angle θ ∈ [−π, π). The application of the same sequence of gates from the integer case results in a quantum state whose state vector consists of the inner product between G(θ) and the Fourier bases G( 2π 2 m j), representing a similarity measure between θ and 2π 2 m j, for 0 ≤ j ≤ 2 m − 1, where G(θ) denotes the geometric sequence vector (1, e iθ , . . . , e i(2 m −1)θ ). We will call the operator preparing this state U Fejér because the outcome probability distribution is the Fejér distribution [40]: (21) If θ = 2π 2 m a, for a real number −2 m−1 ≤ a < 2 m−1 , then U Fejér (θ) |0 m prepares a state whose two most likely measurement outcomes are the closest two integers to a. The probability of measuring one of them is at least 81% [41]. Fig. 11 shows the probability distribution of the outcomes for multiple values of a where m = 4.

C Hardware Specifications
The error map for ibmq_toronto in Fig. 15 was generated on the day of experimentation using Qiskit's visualization library [27].