Analysis and Experimental Demonstration of Amplitude Amplification for Combinatorial Optimization

1  Introduction

Quantum Amplitude Amplification (QAA) is the generalization of Grover’s algorithm, first proposed nearly three decades ago [1] as a means of searching for a marked subset of an unstructured database. This problem is framed as having access to a black-box boolean oracle function f:{1,2,3,...,N}→{0,1} mapping elements of the database to either 0 or 1, for which the solver must find the instance(s) f(x)=1 in the fewest calls to f. An optimal classical algorithm requires O(N) calls to the oracle function. Encoding f as a unitary operation, Grover’s algorithm solves this problem in only O(N) calls to f using a quantum computer. It was shown shortly thereafter that Grover’s algorithm is optimal [2,3], and later that a deterministic (fixed-point) quantum search for marked states is possible with the same speedup [4–7].

Fast-forward nearly 30 years, and the mathematical framework of Amplitude Amplification has been analyzed and generalized in numerous ways [8–15], used to solve various problems beyond database searching [16–23], and incorporated into subroutines for other quantum algorithms [24–27]. Many of these theoretical advancements were developed in the absence of physical hardware, but simultaneously there have been several successful implementations of QAA across multiple quantum technology platforms [28–35].

In this study we focus on a particular variant of QAA which solves combinatorial optimization problems [36–42]. This formulation of QAA uses an oracle operator similar to the phase-separator operator in QAOA [43–49], applying phases proportional to every solution of a discrete cost function (Hamiltonian), which we call cost oracles. Although we shall use the name oracle in this study to keep with normal convention, we stress that these operators contain no black-box element to them. The motivation for these oracle operations is their quantum circuit efficiency, addressing a common point of criticism of implementing Grover’s [50–52], namely the requirement of full N-Toffoli gate operators or quantum dictionaries [21,53] for the oracle.

QAA as defined in this study consists of applying two alternating operations: oracle and diffusion, each containing a free parameter which can in principle vary with each application. One way to determine these values is through measurement feedback and a classical optimizer, effectively QAOA using diffusion as the mixer [47–49]. Alternatively, it has been shown that cost oracle QAA can succeed in the same manner as Grover’s, using a single parameter value for each operator [37–42]. As in this study, using π as the diffusion free parameter setting has been shown to be capable of achieving high probabilities of measurement for optimal solutions [38,39,41,42] (90%+).

The challenge of QAA as described above is that every combinatorial optimization problem results in a unique oracle operation, which in turn requires determining the optimal free parameter setting(s) for the oracle, diffusion, and total iterations k. For both regular [1] and deterministic (fixed-point) Grover’s [5–7] optimal values for each are exactly computable, while here we show that unconstrained linear cost functions (no quadratic or higher terms) produce a symmetry of solutions which can be used to derive an exact formula for predicting oracle parameter settings. Using simulations of problem sizes up to 40 qubits we show peak probabilities and iterations k for finding all possible solutions, highlighting that algorithmic performance for globally optimal solution identification becomes increasingly Grover-like with problem size.

We conclude with experimental results showing a single iteration of generalized QAA up to 5 qubits, run on both IBMQ and IonQ, with and without noise mitigation techniques implemented by the respective vendors [54,55]. Each experiment is designed around varying the free parameter in either the oracle or diffusion operator, comparing the measured probabilities of each basis state with their predicted theoretical values. Our results show the first ever experimental demonstration of generalized QAA using cost oracles.

2  Amplitude Amplification

We begin by defining the QAA algorithm, given in Algorithm 1.

Algorithm 1 applies to both parametrized variations of QAA using Grover’s oracle UG(ϕ) and cost oracles UC(ps). Starting from the N-qubit equal superposition state |s⟩ defined as

|s⟩=12N∑i2N|Zi⟩,(1)

where Zi is a length-N binary string, we perform k iterations of alternating oracle (UG or UC) and diffusion (Us) operations until concluding with a measurement on |Ψ⟩ in the computational basis (the Bloch sphere z-axis of each qubit). Ideally, the result of this measurement yields the target state(s) |Zi⟩ solving the problem encoded by the oracle. The optimal runtime of Algorithm 1 using UG(π) together with Us(π) is k≈π42N/Nm for Nm solutions [1–3].

2.1 Oracle Operations

The Grover oracle UG(ϕ) is a quantum operation defined by the unitary operator

UG(ϕ)|Zi⟩={eiϕ|Zi⟩if |Zi⟩∈marked,|Zi⟩otherwise,(2)

which applies a phase of eiϕ to marked target states. The number of marked states, which we denote Nm, can be as few as one and typically no more than 2N/4 [13]. Realizing this fully entangling operation on quantum hardware is a known challenge [50–52].

We shall define M and N as sets containing all marked and non-marked Zi, respectively. The Grover oracle subdivides the state |Ψ⟩ into a two-dimensional subspace spanned by the orthogonal basis states in M and N as defined by

|m⟩≡1Nm∑Zi∈M|Zi⟩,|n⟩≡1Nn∑Zi∈N|Zi⟩.(3)

For an N-qubit system one typically has Nn+Nm=2N, but more generally this summation can be any integer when using qudits [56] or qubits [57]. Writing |Ψ⟩ in terms of the collective states |n⟩ and |m⟩ allows the oracle operation to be rewritten as

UG(ϕ)=|n⟩⟨n|+eiϕ|m⟩⟨m|.(4)

Now we shall extend the formalism of Eqs. (3) and (4) to oracles which apply more than two unique phases, which is how we define cost oracles in this study. Given a cost function C(Z), we define Ci as the value obtained from evaluating the bit string Zi (C(Zi)=Ci). The result of applying a cost oracle to any basis state is given by

UC(ps)|Zi⟩=eiC(Zi)⋅ps|Zi⟩.(5)

Next we can generalize the notion of the collective states |n⟩ and |m⟩ to |Ci⟩, defined as

|Ci⟩≡1Ni∑Zi∈Ci|Zi⟩.(6)

Analogous to |n⟩ and |m⟩, which are defined by the sets of binary strings N and M in Eq. (3), here each collective state |Ci⟩ is defined by the corresponding set Ci. Specifically, each Ci is the set of all Zj which evaluate to the same cost function value C(Zj)=Ci. We can then rewrite the cost oracle as

UC(ps)=∑iD|Ci⟩⟨Ci|eiCi⋅ps.(7)

The unitary operator UC has one free parameter ps, which we use to stand for the phase scale [41,42] (equivalent to γ in the QAOA literature [47–49]). Analogous to UG, applying UC subdivides |Ψ⟩ into a D-dimensional subspace spanned by the orthogonal collective states |Ci⟩, where D is the number of unique evaluations obtainable from C(Z). Importantly, the quantum circuit construction of UC does not require any knowledge of these collective states or the |Zi⟩ contained within them (see [41,42] and Appendix E), which are assumed to be unknown to the solver.

To conclude, D=2 is by definition a Grover oracle, while D∈[3,2N] defines cost oracles, where D is problem dependent on C(Z). Note that Eqs. (6) and (7) reduce to the form of (3) and (4) for the case D=2, ps=ϕ, {|C1⟩,|C2⟩}={|n⟩,|m⟩}, {C1,C2}={0,1}, {N1,N2}={Nn,Nm}. Thus, our equations for cost oracles are a generalization of Grover’s to higher dimensions D>2.

2.2 Diffusion

We now turn to the diffusion operator Us(θ), defined by

Us(θ)=I−(1−eiθ)|s⟩⟨s|,(8)

corresponding to the quantum circuit construction shown in Fig. 1. This is the standard construction of Us(θ) as originally proposed by Grover [1], also referred to as the Grover mixer in QAOA [46,47]. This multiqubit operation is the N-qubit equivalent to P(θ) shown in Eq. (9), applying a phase of eiθ to the basis state |1⟩⊗N, and no phase to all other |Zi⟩.

P(θ)=[100eiθ](9)

Figure 1: Quantum circuit for the diffusion operator Us(θ) acting on N qubits.

Because oracle operations only apply phases, in QAA the diffusion operator is solely responsible for increasing and decreasing amplitude magnitudes (probabilities). The manner in which it does so is proportional to α¯, the arithmetic mean of all amplitudes αi (amplitudes of each basis state: αi=⟨Zi|Ψ⟩), given by

α¯=12N⟨s|Ψ⟩=12N∑i=12N⟨Zi|Ψ⟩(10)

=12N∑i=1DNi⋅αi,(11)

while

⟨Zj|Us(θ)|Ψ⟩=αj−(1−eiθ)α¯.(12)

Eq. (12) is the change in amplitude experienced by each basis state |Zi⟩ resulting from diffusion. The typical choice for diffusion is θ=π, which from Eq. (12) results in Δαi=αi−2α¯. This change in amplitude is then maximal when αi and α¯ are π phase different, which is the situation created by UG(π) for the marked state(s) [1]. Next in Section 3, we show that this same π phase matching condition can be achieved for UC encoding unconstrained linear cost functions.

3  Solving Linear Cost Functions

For UG, the optimal ϕ, θ, and k values for Algorithm 1 are determined by the size of the unmarked and marked subspaces Nn and Nm [1,6,7]. For cost oracles UC, these parameters are determined by Ci and Ni. In practice however, these values are unknown except through evaluations C(Z), so a realistic implementation of QAA for combinatorial optimization needs to be able to determine ps, θ, and k values without excessive classical precalculation. Here we demonstrate one such example for unconstrained linear cost functions and discuss the mathematical properties which make it possible.

3.1 Problem Symmetry and ᾱ

Consider the case of a set of N randomly selected, real-valued weights Wi with W={W1,W2,...,WN}. We define linear cost functions in this study as

C(Zi)=∑i=1NWi⋅zi,(13)

where each binary variable zi is assigned to one weight Wi. Implementing C(Z) according to Eq. (13) as a cost oracle UC(ps) is achieved through single qubit phase gates P(Wi⋅ps) on each qubit [41,42] (Appendix E). The key feature of this C(Z) that allows one to predict optimal ps values is the symmetric distribution of the solution space about the arithmetic mean C¯, given by

C¯=12N∑i2NCi(14)

=C(Zj)+C(¬Zj)2.(15)

Eq. (14) is generic to all C(Z), while (15) is a property of linear C(Z) as defined by Eq. (13) (proof in Appendix C). For linear C(Z) the value of C¯ can be obtained by evaluating any Zj and its bitwise inverse ¬Zj. The significance of this symmetry of Ci solutions about C¯ is that the resulting mean amplitude after the first oracle operation UC(ps)|s⟩ is

α¯=|α¯|eiC¯⋅ps,(16)

whereby the phase of α¯ is proportional to C¯ (Appendix C). Even though the value of α¯ itself requires complete information of all |Ci⟩, for linear C(Z) knowing the phase of α¯ does not. Thus, using Eq. (16), we can derive a formula for ps such that the phase of α¯ is made to be exactly π different from the phase applied to any target collective state |Ci⟩ of one’s choosing.

C¯⋅ps−Ci⋅ps=±π,(17)

ps=±πC¯−CiforCi≠C¯(18)

Eq. (18) ensures a π phase difference between α¯ and the target state |Ci⟩ after the first oracle operation, while in Appendix C we show that this phase difference holds for subsequent iterations when using Us(π). A visualization of this property is shown in Fig. 2, displaying a series of complex plane plots for |Ψ⟩ after applying UC(ps) according to Eq. (18) for the first five iterations, corresponding to the 10-qubit C(Z) given in Appendix A, Fig. A1 (W1) as well as Fig. A2. As indicated by the axes scale below each plot, the probability of measuring the target state is increasing with each subsequent iteration in a manner exactly analogous to Grover iterations Us(π)UG(π) (π phase between |m⟩ and α¯), but for an oracle encoding a cost function Us(π)UC(ps).

Figure 2: Five complex plane plots of |Ψ⟩ after the application of the operators shown above each panel, corresponding to the 10-qubit C(Z) given in Appendix A. The value of ps used in UC is from Eq. (18), for Ci=2. Each colored point in the plots represents the amplitudes contained within one collective state |Ci⟩, with Ni indicated by the accompanying color scale. The extrema states |Cmin⟩ and |Cmax⟩ are plotted as square and triangle markers, respectively. In each panel the value of α¯ is depicted by the red ×, and the origin of the complex plane with a black +. A red-dotted line is drawn in each panel from α¯ to the collective state |Ci=2⟩, illustrating their π phase difference via the line’s intersection with the origin.

3.2 Algorithmic Performance

The success of QAA for solving a combinatorial optimization problem encoded as UC is dependent on correctly choosing the three parameters ps, θ, and k. A reliable means of determining these parameters for problems such as QUBO or harder [37–42] remains an open research question, while here shall show that Eq. (18) for ps together with Us(π) is sufficient for reliable algorithmic performance.

For the 20-qubit linear C(Z) composed of the integer weights W2 given in Appendix A, Fig. 3 shows the joint peak probabilities (solid-colored lines) obtainable for the ten most optimal |Ci⟩ states as a function of ps. Specifically, each colored curve represents the combined probability of measuring |Ci=T⟩ or its inverse |Ci=2C¯−T⟩ (for example, |Ci=−223⟩ and |Ci=194⟩ shown in blue). Over the range of ps values shown, each peak probability was obtained by simulating Algorithm 1 up to the iteration k where it peaks for the first time (see Fig. 4). The resonance-like shape of these curves is a well-understood feature of Grover’s algorithm [8,58–60], also observed for harder combinatorial problems such as QUBO [41,42], and emerges from the condition for constructive interference for the target state (see Appendix D, Fig. A3 for more details).

Figure 3: (solid-color lines) Combined peak probability as a function of ps for the ten highest (lowest) |Ci⟩ states, for the N=20 linear C(Z) composed of the weights W2 from Appendix A. (gray-dashed lines) The ps values produced from Eq. (18) for the integers T∈[−223,−209], with Ni values reported atop each line.

Figure 4: Simulations of QAA showing probability vs. iterations k for qubit sizes N=10,20,30,40 (top right corner), comparing the joint probability of |Ci=2⟩ and its inverse (solid-black lines) to |m⟩ for standard Grover’s (dashed-red lines) for Nm=2. The ps value used for each UC is from Eq. (18), while the weights of each linear C(Z) are W1 from Appendix A.

Accompanying the curves in the figure are vertical gray lines corresponding to ps values obtained from Eq. (18) for the integers T∈[−223,−209]. Above each line is the corresponding Ni value for each |Ci=T⟩ state (or zero if no such solution exists). As evidenced by the closeness of the lines to each of the ten probability curve peaks, linear C(Z) with integer weights is the ideal scenario for QAA using Eq. (18). Specifically, a solver knows to only check integer values for Eq. (18). Although it is a simple problem instance, it demonstrates that at 20 qubits an exact equation for ps is capable of producing measurement probabilities of 60%–75%+ for |Ci⟩ near the extrema. However, these probabilities are low compared to larger qubit problem sizes, such as those shown in Fig. 4.

Fig. 4, together with Fig. 5 illustrate two key features of QAA algorithmic performance. Firstly, as qubit size increases the achievable probabilities of solutions near the extrema, and the iterations k to reach these probabilities both become increasingly Grover-like. Secondly, within a single problem instance the peak achievable probabilities of each |Ci⟩ decrease for solutions further away from the extrema. This first feature is an important strength for cost oracle QAA, illustrating UC(ps) is capable of producing comparable probabilities to UG(ϕ) without the need for expensive oracle circuits [21,50–53], especially for extrema Ci which are typically the solutions of interest for combinatorial optimization problems.

Figure 5: (top) Plot of peak achievable probabilities for each |Ci⟩ state (colored circles) as a function of σ(ps) using ps from Eq. (18), for the N=40 cases W1 and W3 given in Appendix A. (bottom) Plot of iterations k corresponding to the probabilities shown in the top plot. For comparison, lines indicating the number of iterations for Grover’s (UG(π)) using various Nm are plotted. Beside each line is Δk, showing a 5% increase in k from Grover’s.

Focusing now on Fig. 5, the top plot shows the complete spectrum of achievable probabilities for all |Ci⟩ for two 40-qubit linear C(Z) using ps from Eq. (18), while the bottom plot shows the corresponding number of iterations k needed to achieve these probabilities. For completeness, given below in Eq. (19) is σ(ps), the standard deviation of all Ci solutions after scaling by ps, which is a convenient way to plot all Ci solutions for both C(Z) on the same axis.

σ(ps)=∑i2N(Ci⋅ps−C¯⋅ps)22N19)

The significance of Fig. 5 is the bottom plot, particularly the black lines which show the value of k for standard Grover’s for various values of Nm. Beside each black line is a 5% increase in k, labeled Δk, illustrating that for these 40-qubit problem instances the optimal number of iterations for cost oracle QAA is consistently around 5% more than standard Grover’s. Thus, the results shown in Figs. 3–5 illustrate a special case of QAA, oracles encoding linear C(Z) composed of integer weights, whereby a solver can reliably determine near optimal parameter settings for ps, θ, and k. These results are particularly encouraging because the same underlying features have been observed for QAA using oracles encoding QUBO [37–42], suggesting that these harder problem instances may also have reliable strategies for approximating optimal parameter values as well. In particular, the important condition which led to Eq. (18) is solution space symmetry, which may be applicable to some NP-hard problem instances.

4  Experimental Demonstration

To conclude this study here we present experimental results demonstrating current hardware progress towards the realization of generalized QAA [12] on both IBMQ’s superconducting qubits as well as IonQ’s trapped ion qubits. Each experiment is a single iteration of Us(θ)UC(ps) or Us(θ)UG(ϕ), varying either θ or ps across 100 different values and tracking the measured probabilities of each basis state |Zi⟩. At each parameter value we ran the quantum circuit 10,000 times, totaling one million circuit executions per experiment per qubit size. Our results confirm that these probabilities match the theoretically predicted values from our equations, to our knowledge the first ever such experimental demonstration of generalized QAA using cost oracles.

In total we conducted three unique experiments, which we refer to by number, with quantum circuits for each experiment given in Appendix E, Figs. A4 and A5 (for circuit depth and further details on circuit decomposition for UG and Us see Figs. A6–A8). Each experiment was implemented for qubit sizes N=2,3,4,5, run on five unique devices: IonQ’s aria1 ion trap, IBMQ’s fez and torino (heron), and IBMQ’s kyiv, and brisbane (eagle) superconducting processors, with and without error mitigation techniques provided by the respective hardware vendors [54,55]. The implementation of the multicontrolled phase gate CN-P(θ) used for both Us(θ) and UG(π) is adapted from [61,62], discussed in detail in Appendix E

4.1 The First Iteration

Here we present equations for the theoretical quantum states of experiments 1–3, corresponding to the first iteration of QAA using Grover and cost oracles, analogous to [7,12]. We begin with Eqs. (20) and (21) below.

α¯(ϕ)=(12N)32⋅(Nn+eiϕNm)(20)

α¯(ps)=(12N)32⋅∑jDeiCj⋅psNj,(21)

which are the average amplitudes α¯ (Eq. (11)) as a function of ϕ and ps for the states UG(ϕ)|s⟩ and UC(ps)|s⟩, respectively. Using these equations in combination with (12), we obtain the quantum states after the first iteration for both oracle types given by

Us(θ)UG(ϕ)|s⟩=Nn(12N−α¯(ϕ)(1−eiθ))|n⟩+Nm(eiϕ2N−α¯(ϕ)(1−eiθ))|m⟩(22)

Us(θ)UC(ps)|s⟩=∑jDNj(eiCj⋅ps2N−α¯(ps)(1−eiθ))|Cj⟩.(23)

Experiments 1 and 2 correspond to Eq. (23), whereby in exp. 1 we prepared and measured 100 quantum circuits for the state Us(π)UC(ps) (varying the parameter ps) and likewise Us(θ)UC(1) for exp. 2 (varying θ). Experiment 3 corresponds to states of Eq. (22) Us(θ)UG(π), Nm=1, varying θ over the full range of 2π. Taking the squared magnitude of the amplitudes given in Eqs. (22) and (23) yields the expected theoretical measurement probabilities for each collective state, which we experimentally verify in the coming subsections.

4.2 Performance Metric

Because each experiment is a composition of 100 different quantum circuits, here we shall briefly describe the quantity f, which is a metric that pools the results from all 100 circuits into a single value comparing experimentally measured probabilities vs. theory, full mathematical description given in Appendix B. The quantity f is a ratio comparing experimentally observed measurement counts of each |Zi⟩ to probabilities predicted by Eqs. (22) and (23): f=1 means that the basis state was observed with the same probability as theoretically expected, while f=0 means an observed probability of 1/2N (the equal superposition state |s⟩), i.e., a completely decohered quantum state. Values of f<0 occur when a basis state’s measured probability was found to be on the opposite side of 1/2N as predicted from theory, signaling that the qubits had reached an unintended final |Ψ⟩ not attributable to decoherence. Shown in Fig. 6 is an example plot illustrating three theoretical values of f (no experimental data) for experiment 2, the quantum state |Ψ⟩=Us(θ)UC(1) over the full 2π range of θ (see Fig. A1 for the |Ci⟩ corresponding to each color).

Figure 6: (solid-colored) Probabilities representing f=1 for each basis state as a function of diffusion angle θ for experiment 2, N=2, as predicted from Eq. (23) over the range θ∈[0,2π]. For comparison, probabilities corresponding to f=0.7 (dashed-colored) and 0 (black-dash-dotted) are also plotted.

4.3 Improvement from Error Mitigation

A recent advancement in both IBMQ and IonQ’s commercially available quantum computing platforms is the addition of error mitigation as an option for users [54,55]. More specifically each service offers both Pauli gate twirling and dynamical decoupling. In this study we report on qubit performance for error mitigation ‘on’ and ‘off’, which refers to the use of both techniques together by the respective hardware vendors, or neither. An example of difference in performance for on vs. off is illustrated in Fig. 7. In each plot the solid-colored lines represent f=1 (Eq. (23)), while the colored data points show the measured probability counts for each |Zi⟩ state (colored according to the |Ci⟩ each basis state belongs to), with fexp for the entire experiment reported in the top right (fexp given in Appendix B). Examples of f<0 for individual |Zi⟩ can be seen in the top plot (error mitigation off), illustrating noisy |Ψ⟩ states with clear underlying structure (ripe for error mitigation/correction techniques) as opposed to complete decoherence.

Figure 7: A comparison of measurement results for exp. 1, N=3, on IBMQ’s brisbane with (bottom) and without (top) error mitigation. The solid-colored lines represent probabilities for f=1 (Eq. (23)), while the colored data points are the experimental measurement probabilities obtained for each basis state. Reported in the top right corner is the cumulative f for the experiment, fexp from Eq. (A9) in Appendix B.

4.4 Experimental Results

Given in Fig. 8 below are the full results for exps. 1–3, with error mitigation techniques turned on (top) and off (bottom). The numerical values shown in each table are fexp from Eq. (A9) (Appendix B), the average f value from all 2N basis states. The color of each cell indicates the ranked performance of each device within a single experiment, from best to worst: [blue, green, yellow, orange, red]. Cells that are instead colored gray are for values fexp<0.15.

Figure 8: Two tables summarizing the experimental results of exps. 1–3 for error mitigation on (top) and off (bottom). The numerical value of each cell is fexp from Eq. (A9), with colors indicating ranked performance. Gray cells indicate values less than 0.15.

A summary of the results shown in Fig. 8 is as follows. The effect of error mitigation provided by both vendors improved fexp in 55 out of the 60 experiment runs, producing an overall average increase in fexp by 0.39. Regarding the performance of the individual devices, IBMQ’s heron processors fez and torino ranked first and second in 13 out of 18 experiments for qubit sizes N≥3. However, with error mitigation IonQ’s aria1 trap was first in all three 2-qubit experiments, including the highest overall fexp value of 0.959 as well as the highest 3-qubit value of 0.922. Shown in Fig. 9 are the highest fexp achieved for all experiments and qubit sizes.

Figure 9: Experiment results showing the highest fexp achieved for each combination of experiment and qubit size. The data points in each plot are Meas.(|Zi⟩) from Eq. (A4), colored according to |Ci⟩ states (see Appendix A). The solid-colored lines in each plot represent f=1 for each collective state, obtained from Eqs. (22) and (23). Reported in the top right corner is the cumulative f for the experiment, fexp from Eq. (A9) in Appendix B, as well as fm for exp. 3 (f for only the blue data points corresponding to |m⟩, the all |1⟩ state).

5  Conclusion

In this study, we have generalized the common |n⟩ and |m⟩ collective state formalism of Grover’s algorithm to cost oracles encoding combinatorial optimization problems [37–42], and shown that for the special case of unconstrained linear C(Z) an exact equation exists for determining the free parameter value of the oracle operator. In Section 3, we used simulations of QAA up to 40 qubits to demonstrate the closeness in algorithmic performance of cost oracle QAA using our equation for ps with standard Grover’s. And finally, in Section 4, we verify our derived equations of generalized QAA [7,12] on two different qubit technologies, showcasing progress in state-of-the-art commercial quantum devices and their respective error mitigation capabilities.

Outlook and Future Research

The motivation for studying cost oracle QAA [37–42] stems from the quantum circuit efficiency of UC [63,64] as compared to the full N-Toffoli gates needed for UG [21,23,31,32,50–53], especially when considering the additional overhead needed for error correction [65,66]. For both oracle cases however, implementing the diffusion operator still remains as a technological bottleneck. Progress in circuit depth efficiency [67–70] is one solution towards achieving large N-Toffoli gates, as well as going beyond 2-qubit gates [71–74]. Additionally there is also the potential for realizing QAA outside of the gate-based model, such as quantum programmable processors [75,76] for implementing fixed unitary operators such as diffusion.

In this work we focus on linear cost functions composed of integer weights to demonstrate ideal conditions for cost oracle QAA. Specifically, integer weights limit the possible values obtainable from the cost function, which in turn are used for selecting the free parameter setting of the oracle via our derived equation. The same algorithmic performance as demonstrated in this study has been shown for QUBO and more complex problems [38–42], but determining the free parameter setting of the oracle for these harder problem instances remains an open research question. Alternatively, one way to sidestep this issue at the cost of increased runtime is the Grover-Mixer QAOA [47–49] style of determining optimal parameter settings for each iteration via measurement results and a classical optimizer. The advantage of QAA as demonstrated in Section 3 is efficiency, requiring ideally only a few measurements to find the desired solution of a cost function.

Lastly, the experimental results of Section 4 show current commercial hardware’s ability to achieve QAA up to N=5 qubits, in essence demonstrating their ability to implement the N-Toffoli gate operation. Based on experimental trends of the last 5+ years [29–35] and progress in circuit efficiency [67–70], it does not appear likely that this number N will dramatically increase in the next several years without a significant technological breakthrough (extending N into the 10s or 100s of qubits). By contrast, current hardware (especially superconducting qubits which support parallel gates) is already powerful enough to implement cost oracle operations up to hundreds of qubits, making diffusion the technological bottleneck limiting QAA. Therefore, another avenue for future research is exploring QAA performance using less than full N-qubit diffusion, instead replaced by multiple low qubit-sized diffusions in parallel. This would significantly reduce circuit depth, but also introduces more θ degrees of freedom which complicate the algorithm. If small diffusions in parallel can produce probabilities comparable to full N-qubit diffusion, then QAA could potentially be a viable near term quantum algorithm alongside QAOA for solving combinatorial optimization problems.

Acknowledgement: Any opinions, findings, conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of AFRL. This project was supported in part by an appointment to the NRC Research Associateship Program at AFRL, administered by the Fellowships Office of the National Academies of Sciences, Engineering, and Medicine.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Daniel Koch. Methodology, Daniel Koch, Brian Pardo. Formal Analysis, Daniel Koch, Kip Nieman. Experimental Data Generation, Daniel Koch, Brian Pardo. Writing, Revewing and Editting, Daniel Koch, Brian Pardo, Kip Nieman. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data and code files that support the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

Appendix A Cost Functions & Distributions

Found in Sections 3 and 4 are references to a particular problem cases of linear cost functions which can be found here. Given in Eq. (A1) is W1, the set of positive integers from 1 to N, which are the sets of weights for Figs. 2, 4, and 5, as well as experiments 1 and 2 in Section 4. Also shown in Fig. 5 is a second 40-qubit cost function corresponding to the weights W3 given in Eq. (A3).

W1={1,2,3,...,N}(A1)

W2={−44,−35,−33,−32,−23,−20,−11,−11,−10,−4,2,6,9,11,11,17,21,34,40,43}(A2)

W3={−731,−722,−676,−668,−663,−662,−564,−563,−555,−409,−209,−189,−135,−43,1,3,28,48,73,127,139,156,160,286,307,308,427,461,490,512,548,551,568,583,589,642,776,917,929,948,949}(A3)

For the results of experiments 1 and 2 found in Section 4, which use UC(ps) encoding linear C(Z) according to Eq. (A1), Fig. A1 below provides a table all of the collective states |Ci⟩ along with the basis states |Zi⟩ contained within them.

Figure A1: A table displaying all of the collective states |Ci⟩ for experiments 1 and 2, corresponding to a linear C(Z) composed of the weights W1 from Eq. (A1) up to N=5, as well as the basis states |Zi⟩ contained within them.

Appendix B Fidelity Metric

Here we discuss in full mathematical detail the calculation of f, the metric of performance used in Section 4. We begin with Eq. (A4) for Meas.(|Zi⟩), which is simply the experimentally observed probability of a basis state. ‘Shots’ refers to the number of times a quantum circuit was run, producing that number of measurement results, which in this study is 10,000. Given in Eq. (A5) is the the quantity ΔPi, which is the difference in probability between measured counts of a |Zi⟩ basis state (Eq. (A4)) and its theoretical expected probability given by Eq. (23) or (22).

Meas.(|Zi⟩)=Counts(|Zi⟩)shots(A4)

ΔPi=|Meas.(|Zi⟩)−|⟨Zi|Ψ⟩|2|(A5)

The problem with using raw probability difference to evaluate performance is that good or bad is relative for each basis state. This was the motivation for introducing f, which compares each ΔPi against probabilities corresponding to a completely decohered equal superposition state, given below in Eq. (A6) as ΔP~i.

ΔP~i=|12N−|⟨Zi|Ψ⟩|2|(A6)

The computation of f is based on the ratio between root mean square calculations of ΔPi and ΔP~i values from all 100 circuits, as shown in Eq. (A7) (note that i refers to the basis state while j refers to one of the 100 quantum circuits). For each basis state |Zi⟩ we compute fi according to Eq. (A8) (where RMS~ is obtained from using ΔP~i values in Eq. (A7)). Averaging all 2Nfi together yields fexp, given in Eq. (A9).

RMSi=1100∑j100ΔPij2(A7)

fi=1−RMSiRMS~i(A8)

fexp=12N∑i2Nfi(A9)

Figure A2: Histogram of all Ci solutions to a 10-qubit C(Z) using the weights W given in Eq. (A1). (red line) The value C¯ from Eq. (14).

Appendix C Symmetry Proofs for Linear C(Z)

To assist with the forthcoming derivations, given below in Fig. A2 is a histogram of all Ci solutions to the 10-qubit linear C(Z) composed of the set of weights W1 from Eq. (A1). Each black circle shown in the figure corresponds to one collective state |Ci⟩ containing Ni basis states.

Eq. (16) from Section 3 states that for any UC(ps) encoding a linear C(Z) the resulting phase of α¯ is equal to the C¯⋅ps. Here we shall derive the equations leading to this relation, showing explicitly the computation of α¯ after the first two applications of UC(ps). The results which follow are specifically for C(Z) according to Eq. (13) and do not generalize to more complex problems such as QUBO. We begin with Eq. (A10) below which defines a pair of inverse solutions Zi and ¬Zi, two binary strings that have equal and opposite 0 and 1 values on every bit such that their binary sum equals 2N−1 (base-10). For linear C(Z) the evaluation of any pair of inverse Zi is equal to the full sum of W.

Zi+¬Zi=111...1=2N−1(A10)

C(Zi)+C(¬Zi)=∑iNWi=Wsum(A11)

Using Eq. (A11) we can now prove Eq. (11), the mean cost function value C¯ according to Eq. (14) is always equal to half of the sum of any pair of inverse solutions, given by Eqs. (A12)–(A15) below.

C¯=12N∑i2NC(Zi)(A12)

=12N∑i2N−1C(Zi)+C(¬Zi)(A13)

=2N−1⋅Wsum2(A14)

=C(Zj)+C(¬Zj)2(A15)

For the next derivation let us define ΔCi, given in Eq. (A16), which is simply the difference between the mean cost function value and any C(Zi).

ΔCi=C¯−C(Zi)(A16)

Illustrated in Fig. A2 is the symmetry property of linear C(Z) that all solutions come in equal and opposite pairs that are equidistant from C¯. Given in Eqs. (A17)–(A19) below is the proof that these symmetric solutions are in fact pairs of inverse Zi.

C(Zi)+C(¬Zi)=2C¯(A17)

C¯−ΔCi+C(¬Zi)=2C¯(A18)

C(¬Zi)=C¯+ΔCi(A19)

Using Eqs. (A16) and (A19), we can now prove the relation between C¯ and α¯ given in Eq. (16) of Section 3. The following equations use the amplitudes αi corresponding to the state |Ψ⟩=UC(ps)|s⟩, the first oracle application acting on the equal superposition state. By rewriting the summation for α¯ first in terms of inverse Zi pairs and then substituting in ΔCi and C¯ equivalents, we arrive at Eq. (A24).

α¯1=12N∑j2Nαj(A20)

=12N2N∑j2NeiC(Zj)⋅ps(A21)

=1(2N)3/2∑j2N−1(eiC(Zj)⋅ps+eiC(¬Zj)⋅ps)(A22)

=1(2N)3/2∑j2N−1eiC¯⋅ps(e−iΔCj⋅ps+eiΔCj⋅ps)(A23)

=eiC¯⋅ps(2N)3/2∑j2N−12cos⁡(ΔCj⋅ps)(A24)

To summarize, the equations above show that the value of α¯1 after the first application of UC is equal to eiC¯⋅ps times a constant, confirming the relation given by Eq. (16), which in turn is used to derive Eq. (17) for finding ps values which create a π phase difference between α¯ and any |Ci⟩. Here we have shown that the phase of α¯ is proportional to C¯, but the full computation of α¯ still requires complete knowledge of all Ci.

So far we have shown that Eq. (16) holds for the first oracle application, so next we shall prove that it holds for the second as well. Given in Eq. (A25) below are the amplitudes αj′ corresponding to the state |Ψ⟩=Us(π)UC(ps)|s⟩.

αj′=eiC(Zj)⋅ps2N−2|α¯1|eiC¯⋅ps(A25)

Eqs. (A26)–(A29) below show the calculation of α¯2, the mean amplitude following the second oracle application.

α¯2=12N∑j2Nei2C(Zj)⋅ps2N−2|α¯1|ei(C¯+C(Zj))⋅ps(A26)

=12N∑j2N−1(ei2C(Zj)⋅ps+ei2C(¬Zj)⋅ps)2N−2|α¯1|eiC¯⋅ps(eiC(Zj)⋅ps+eiC(¬Zj)⋅ps)(A27)

=ei2C¯⋅ps2N∑j2N−1(e−i2ΔCj⋅ps+ei2ΔCj⋅ps)2N−2|α¯1|(e−iΔCj⋅ps+eiΔCj⋅ps)(A28)

=ei2C¯⋅ps2N−1∑j2N−1cos(2ΔCj⋅ps)2N−2|α¯1|cos(ΔCj⋅ps)(A29)

In Eq. (A29) above we see that the phase of α¯2 is equal to 2C¯⋅ps, but once again the full value of α¯2 requires complete knowledge of every amplitude value. To conclude, we note that the derivation of two iterations given above is general, showing that phase relation given in Eq. (16) holds, also supported by the simulation results of Section 3, but does not constitute a complete proof by induction.

Appendix D Resonance Width in Grover’s

The QAA algorithm using iterations of Us(θ)UG(ϕ) in Algorithm 1 is optimal when both free parameters are [θ,ϕ]=[π,π] [1–3]. However, it is important to understand the behavior of the algorithm when these phases differ, particularly when considering noisy gates which lead to imperfect implementations of the oracle [8,58]. In this appendix we show the relation between peak achievable probabilities as a function of ϕ around π following closely the approach of [59]. The resonance phenomena demonstrated here for Grover’s at ϕ=π is the same as those shown in Fig. 3 for each of the individual |Ci⟩ states [38,41,42]. Eq. (18) from Section 3 produces near optimal ps values for linear C(Z), but in general determining ps values which align with |Ci⟩ resonance peaks for QUBO and harder problem instances is an open research question [40]. We begin with Eq. (A30) below, the equal superposition state |s⟩ expressed as a column vector. For simplicity the initial amplitudes of |m⟩ and |n⟩ are sin⁡β=Nm/2N and cos⁡β=Nn/2N, respectively.

|s⟩=sin⁡β|m⟩+cos⁡β|n⟩=[sin⁡βcos⁡β].(A30)

Next we need the 2 × 2 matrix form of Us(θ)UG(ϕ), which we call G in Eq. (A31) below.

G≡Us(θ)UG(ϕ)(A31)

=[eiϕ(1+(eiθ−1)sin2⁡(β))(eiθ−1)sin⁡(β)cos⁡(β)eiϕ(eiθ−1)sin⁡(β)cos⁡(β)1+(eiθ−1)cos2⁡(β)](A32)

Since G given in Eq. (A32) is a unitary matrix, it can be represented generically by the matrix decomposition given in Eq. (A33), where U is some unitary 2 × 2 matrix and D is a diagonal matrix with values λ+ and λ−. The matrix U has columns given by the normalized eigenvectors of Eq. (A34), with eigenvalues given in Eqs. (A35) and (A36).

G=U∗DU,(A33)

u+=[e−iϕ2cos⁡xsin⁡x]u−=[−sin⁡xeiϕ2cos⁡x](A34)

λ±=−ei(θ+ϕ)2±w(A35)

cos⁡w=cos⁡(ϕ−θ2)−2sin⁡ϕ2sin⁡θ2sin2⁡β.(A36)

Using Eqs. (A33)–(A36), we can write the general form of the Grover operator G after t iterations, given in Eq. (A37).

Gt=[eiwtcos2⁡x+e−iwtsin2⁡xie−iϕ2sin⁡(wt)sin⁡(2x)ieiϕ2sin⁡(wt)sin⁡(2x)eiwtsin2⁡x+e−iwtcos2⁡x](A37)

Next we need expressions for the angle x, which can be found by using the equation Gu+=λ+u+, or equivalently by setting the off-diagonal elements of U∗GU to 0. Solving this first expression yields the equations for sin⁡x and cos⁡x given in Eqs. (A38)–(A40) below.

sin⁡x=lm−12sin⁡θ2sin⁡(2β)(A38)

cos⁡x=lm−12[sin⁡w+sin⁡(ϕ−θ2)+2cos⁡ϕ2sin⁡θ2sin2⁡β](A39)

lm=(sin⁡θ2sin⁡(2β))2+[sin⁡w+sin⁡(ϕ−θ2)+2cos⁡ϕ2sin⁡θ2sin2⁡β]2(A40)

Using the equations above for Gt, we can compute the amplitude of |m⟩ at any iteration t using Eq. (A41) below.

⟨m|Gt|s⟩=sin⁡βcos⁡(wt)+i(eiϕ/2cos⁡βsin⁡(2x)+sin⁡βcos⁡(2x))sin⁡(wt)(A41)

Taking the absolute value squared of Eq. (A41) yields Pm(t)=|⟨m|Gt|s⟩|2, the probability of measuring |m⟩ after t iterations. In the large 2N limit, this probability simplifies significantly to Eq. (A42) below.

Pm(t)≈sin2⁡(2x)sin2⁡(wt)(A42)

The probability Pm(t) is maximal at t=π2w, which after substituting into the equation above yields Eq. (A43) for Pmax, the maximum achievable probability for |m⟩ as a function of ϕ and θ.

Pmax≈4sin2⁡θ22Nsin2(θ−ϕ2)+4sin⁡θ2sin⁡ϕ2cos(θ−ϕ2)(A43)

And finally, we are interested in the case of θ=π, yielding Eq. (A44) for Pmax as a function of UG oracle angle ϕ. Shown in Fig. A3 are four plots for various problem sizes N, illustrating the resonance behavior sharply peaked around ϕ=π, becoming narrower with increasing problem size N.

Pmax=1sin2⁡ϕ2+2N4cos2⁡ϕ2.(A44)

Figure A3: Peak achievable probability Pmax from Eq. (A44) as a function of oracle phase ϕ for various N-qubit problem sizes. Each plot shows the resonance behavior around ϕ=π for Grover’s algorithm using Us(π) for diffusion searching for a single marked state Nm=1.

Fig. A3 illustrates the degree of precision in ϕ necessary for Grover’s algorithm as problem size N increases, also representing the required accuracy in ps for cost oracles [38,42]. We can quantify this required precision with Eq. (A45) below, the full width half maximum (FWHM) for problem sizes N>2.

δϕ=2cos−1⁡(22N−4)(A45)

.

Appendix E Quantum Circuits

Here we present the quantum circuits for experiments 1–3 from Section 4, with further details on the decomposition of the CN-P(θ) gate given in the next subsection. Shown in Fig. A4 below is the quantum circuit for experiments 1 and 2, corresponding to the state Us(θ)UC(ps)|s⟩ given in Eq. (23). The oracle UC here encodes W1 from Eq. (A1). Experiment 1 uses θ=π while varying ps, while conversely experiment 2 uses ps=1 while varying θ. The value of N′ is given by Eq. (A46) below.

N′=∑i=1Ni=12(N2+N)(A46)

Figure A4: Quantum circuit for experiments 1 and 2.

The quantum circuit for experiment 3 is given in Fig. A5, using UG(π) as the oracle while varying θ in the diffusion operator. The CN-Z operator shown in the figure was implemented using CN-P(π) as discussed in the next subsection. The quantum state produced from this circuit is Eq. (22). Given in Fig. A6 is a table detailing the circuit depth and 2-qubit gate count for all three experiments on IBMQ’s four processors, where circuit depth is defined as parallel gate layers [63].

Figure A5: Quantum circuit for experiment 3.

Figure A6: Circuit depth and number of 2 qubit gates across all experiments and devices.

The most resource-intensive component of both the Grover oracle and the diffusion operator is the fully-entangling CN-P(θ) gate. This multiqubit gate must be transpiled into 1 & 2-qubit instructions, for which there are several approaches. In this work, we opt for a simple decomposition that is not optimized for any particular architecture. We send the same quantum circuit to both IBMQ and IonQ, decomposing CNP(θ) into 2-qubit CX and single qubit P(θ) operations. For IBMQ we transpiled directly to the backend’s qubits, while for IonQ qubit selection and gate transpilation were left to the hardware provider.

The decomposition procedure used in this study is adapted from [61,62]. Consider a sequence of binary variables zi for i∈[1,N]. We want to apply a 2×2 unitary gate U, in our case the phase gate P(θ) from Eq. (9), to the Nth qubit if and only if zi=1 for all i∈[1,N−1]. This operation can be expressed by the identity given in Eq. (A47).

2m−1(z1∧z2∧⋯∧zm)=∑k1zk1−∑k1<k2(zk1⊕zk2)+⋯+(−1)m−1(z1⊕⋯⊕zm)(A47)

The operation above can be implemented using a Gray sequence to determine the control and target qubits. An N-bit Gray code is the sequence of 2N binary numbers with adjacent elements in the sequence differing by a single bit flip. The particular Gray sequence that we use can be determined recursively as follows. We denote a Gray sequence by GN=(g1,g2,…,g2N), with the base case Gray sequence given by G1=(0,1). We can recursively define a Gray sequence given by GN+1=(0GN,1G~N), where G~ is the reverse ordering of G. For example, this produces G2=(00,01,11,10) followed by G3=(000,001,011,010,110,111,101,100). To implement this code as gates, we note that the leading 1 corresponds to the target qubit and the remaining 1’s correspond to the XOR combination for which we apply the gate V or V∗. The parity of the bits determine which is applied: V if even or V∗ if odd. Given in Fig. A7 below is the quantum circuit used in this study for implementing a CN-P(θ) gate, with further details given in Fig. A8.

Figure A7: Quantum circuit for implementing the gray code.

Figure A8: Quantum circuits for each of the unitaries Un as shown in Fig. A7. For achieving an N-qubit CN-P(θ) phase gate, the operators V and V∗ are the single qubit phases gates P(θ′) and P(−θ′), respectively, where θ′=θ/2N−1.

We note that the approach used in this study requires O(2N) CNOT gates, and thus is only practical for small N instances. However, our approach did give shallower circuit depth and fewer gate counts than the native CN-P(θ) implementation by Qiskit at the time of running (pre-Qiskit 2.0). For example, at N=5 our transpilation for diffusion was 48 CNOTs as compared to Qiskit’s MCPhase gate which transpiled to 70.

References

1. Grover LK. A fast quantum mechanical algorithm for database search. In: Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing; 1996 May 22–24; Philadelphia, PA, USA. p. 212–29. [Google Scholar]

2. Boyer M, Brassard G, Høyer P, Tapp A. Tight bounds on quantum searching. Fortschritte Phys. 1998;46:493–506. doi:10.1002/(sici)1521-3978(199806)46:4/5<493::aid-prop493>3.0.co;2-p. [Google Scholar] [CrossRef]

3. Zalka C. Grover’s quantum searching algorithm is optimal. Phys Rev A. 1999;60:2746–51. [Google Scholar]

4. Long GL. Grover algorithm with zero theoretical failure rate. Phys Rev A. 2001;64:022307. doi:10.1103/physreva.64.022307. [Google Scholar] [CrossRef]

5. Grover LK. Fixed-point quantum search. Phys Rev Lett. 2005;95:150501. doi:10.1103/physrevlett.95.150501. [Google Scholar] [PubMed] [CrossRef]

6. Yoder TJ, Low GH, Chuang IL. Fixed-point quantum search with an optimal number of queries. Phys Rev Lett. 2014;113(21):210501. doi:10.1103/physrevlett.113.210501. [Google Scholar] [PubMed] [CrossRef]

7. Roy T, Jiang L, Schuster DI. Deterministic Grover with a restricted oracle. Phys Rev Res. 2022;4:L022013. [Google Scholar]

8. Long GL, Li YS, Zhang WL, Niu L. Phase matching in quantum searching. Phys Lett A. 1999;262(1):27–34. doi:10.1016/s0375-9601(99)00631-3. [Google Scholar] [CrossRef]

9. Høyer P. Arbitrary phases in quantum amplitude amplification. Phys Rev A. 2000;62(5):052304. doi:10.1103/physreva.62.052304. [Google Scholar] [CrossRef]

10. Biham E, Biham O, Biron D, Grassl M, Lidar DA, Shapira D. Generalized Grover’s quantum search algorithms using recursion equations. Phys Rev A. 2001;63(1):012310. doi:10.1103/physreva.63.012310. [Google Scholar] [CrossRef]

11. Byrnes T, Forster G, Tessler L. Generalized Grover’s algorithm for multiple phase inversion states. Phys Rev Lett. 2018;120(6):060501. doi:10.1103/physrevlett.120.060501. [Google Scholar] [CrossRef]

12. Kwon H, Bae J. Quantum amplitude amplification operators. Phys Rev A. 2021;104(6):062438. doi:10.1103/physreva.104.062438. [Google Scholar] [CrossRef]

13. Szablowski PJ. Understanding the mathematics of Grover’s algorithm. Quantum Inf Process. 2021;20(5):191. doi:10.1007/s11128-021-03125-w. [Google Scholar] [CrossRef]

14. Sun Y, Wu LA. Quantum search algorithm on weighted databases. Sci Rep. 2024;14(1):30169. doi:10.1038/s41598-024-81701-7. [Google Scholar] [PubMed] [CrossRef]

15. Suzuki Y, Gluza M, Son J, Tiang BH, Ng NHY, Holmes Z. Grover’s algorithm is an approximation of imaginary-time evolution. arXiv:2507.15065. 2025. [Google Scholar]

16. Brassard G, Høyer P, Mosca M, Tapp A. Quantum amplitude amplification and estimation. In: Contemporary mathematics. Vol. 305. London, UK: AMS; 2002. p. 53–74. [Google Scholar]

17. Childs AM, Goldstone J. Spatial search by quantum walk. Phys Rev A. 2004;70:022314. [Google Scholar]

18. Aaronson S, Ambainis A. Quantum search of spatial regions. Theory Comput. 2005;1:47–79. doi:10.1109/sfcs.2003.1238194. [Google Scholar] [CrossRef]

19. Janmark J, Meyer DA, Wong TG. Global symmetry is unnecessary for fast quantum search. Phys Rev Lett. 2014;112(21):210502. doi:10.1103/physrevlett.112.210502. [Google Scholar] [CrossRef]

20. Anikeeva G, Markovic O, Borish V, Hines JA, Rajagopal SV, Cooper ES, et al. Number partitioning with Grover’s algorithm in central spin systems. PRX Quantum. 2021;2(2):020319. doi:10.1103/prxquantum.2.020319. [Google Scholar] [CrossRef]

21. Gilliam A, Woerner S, Gonciulea C. Grover adaptive search for constrained polynomial binary optimization. Quantum. 2021;5:428. [Google Scholar]

22. Nieman K, Durand H, Patel S, Koch D, Alsing PM. Investigating an amplitude amplification-based optimization algorithm for model predictive control. Digit Chem Eng. 2024;10:100134. [Google Scholar]

23. Zhang Z, Paredes R, Sundar B, Quiroga D, Kyrillidis A, Duenas-Osorio L, et al. Grover-QAOA for 3-SAT: quadratic speedup, fair-sampling, and parameter clustering. Quantum Sci Technol. 2025;10:015022. [Google Scholar]

24. Harrow AW, Hassidim A, Lloyd S. Quantum algorithm for linear systems of equations. Phys Rev Lett. 2009;103(15):150502. doi:10.1103/physrevlett.103.150502. [Google Scholar] [PubMed] [CrossRef]

25. Brassard G, Høyer P, Tapp A. Quantum counting. Vol. 1443. Berlin/Heidelberg, Germany: Springer; 1998. p. 820–31. [Google Scholar]

26. Grinko D, Gacon J, Zoufal C, Woerner S. Iterative quantum amplitude estimation. npj Quant Inf. 2021;7:52. [Google Scholar]

27. Muser T, Zapusek E, Bellis V, Reiter F. Provable advantages of kernel-based quantum learners and quantum preprocessing based on Grover’s algorithm. Phys Rev A. 2024;110(3):032434. doi:10.1103/physreva.110.032434. [Google Scholar] [CrossRef]

28. Chen K, Li CM, Zhang Q, Chen YA, Goebel A, Chen S, et al. Experimental realization of one-way quantum computing with two-photon four-qubit cluster states. Phys Rev Lett. 2007;99(12):120503. doi:10.1103/physrevlett.99.120503. [Google Scholar] [PubMed] [CrossRef]

29. Figgatt C, Maslov D, Landsman KA, Linke NM, Debnath S, Monroe C. Complete 3-qubit Grover search on a programmable quantum computer. Nat Commun. 2017;8:1918. doi:10.1038/s41467-017-01904-7. [Google Scholar] [PubMed] [CrossRef]

30. Mandviwalla A, Ohshiro K, Ji B. Implementing Grover’s algorithm on the IBM quantum computers. In: Proceedings of the 2018 IEEE International Conference on Big Data; 2018 Dec 10–13; Seattle, WA, USA. p. 2531–7. [Google Scholar]

31. Zhang K, Rao P, Yu K, Lim H, Korepin V. Implementation of efficient quantum search algorithms on NISQ computers. Quantum Inf Process. 2021;20(7):233. doi:10.1007/s11128-021-03165-2. [Google Scholar] [CrossRef]

32. Zhang K, Yu K, Korepin V. Quantum search on noisy intermediate-scale quantum devices. Europhys Lett. 2022;140(1):18002. doi:10.1209/0295-5075/ac90e6. [Google Scholar] [CrossRef]

33. Pokharel B, Lidar DA. Better-than-classical Grover search via quantum error detection and suppression. npj Quantum Inf. 2024;10(1):23. doi:10.1038/s41534-023-00794-6. [Google Scholar] [CrossRef]

34. Thorvaldson I, Poulos D, Moehle CM, Misha SH, Edlbauer H, Reiner J, et al. Grover’s algorithm in a four-qubit silicon processor above the fault-tolerant threshold. Nat Nanotechnol. 2025;20(4):472–7. doi:10.1038/s41565-024-01853-5. [Google Scholar] [CrossRef]

35. AbuGhanem M. Characterizing Grover search algorithm on large-scale superconducting quantum computers. Sci Rep. 2025;15(1):1281. doi:10.1038/s41598-024-80188-6. [Google Scholar] [PubMed] [CrossRef]

36. Shyamsundar P. Non-Boolean quantum amplitude amplification and quantum mean estimation. Quantum Inf Process. 2023;22(12):423. doi:10.1007/s11128-023-04146-3. [Google Scholar] [CrossRef]

37. Satoh T, Ohkura Y, Meter RV. Subdivided phase oracle for NISQ search algorithms. IEEE Trans Quantum Eng. 2020;1:3100815. doi:10.1109/tqe.2020.3012068. [Google Scholar] [CrossRef]

38. Benchasattabuse N, Satoh T, Hajdušek M, Van Meter R. Amplitude amplification for optimization via subdivided phase oracle. arXiv:2205.00602. 2022. [Google Scholar]

39. Tani K, Tsuchiya S, Tani S, Takeuchi Y. Quantum algorithm for unstructured search of ranked targets. Phys Scr. 2025;100(7):075114. doi:10.1088/1402-4896/ade377. [Google Scholar] [CrossRef]

40. Zhukov A, Lebedev A, Pogosov W. Grover’s search meets Ising models: a quantum algorithm for finding low-energy states. Comput Phys Commun. 2025;313:109627. [Google Scholar]

41. Koch D, Cutugno M, Karlson S, Patel S, Wessing L, Alsing PM. Gaussian amplitude amplification for quantum pathfinding. Entropy. 2022;24(7):963. doi:10.3390/e24070963. [Google Scholar] [PubMed] [CrossRef]

42. Koch D, Cutugno M, Patel S, Wessing L, Alsing PM. Variational amplification for solving QUBO problems. Quantum Rep. 2023;5(4):625–58. doi:10.3390/quantum5040041. [Google Scholar] [CrossRef]

43. Farhi E, Goldstone J, Gutmann S. A quantum approximate optimization algorithm. arXiv:1411.4028. 2014. [Google Scholar]

44. Hadfield S, Wang Z, O’Gorman B, Rieffel EG, Venturelli D, Biswas R. From the quantum approximate optimization algorithm to a quantum alternating operator ansatz. Algorithms. 2019;12(2):34. doi:10.3390/a12020034. [Google Scholar] [CrossRef]

45. Maciejewski FB, Bach BG, Dupont M, Lott PA, Sundar B, Bernal Neira DE, et al. A multilevel approach for solving large-scale QUBO problems with noisy hybrid quantum approximate optimization. In: Proceedings of the 2024 IEEE High Performance Extreme Computing Conference (HPEC); 2024 Sep 23–27; Wakefield, MA, USA. p. 1–10. [Google Scholar]

46. Blekos K, Brand D, Ceschini A, Chou C-H, Li R-H, Pandya K, et al. A review on quantum approximate optimization algorithm and its variants. Phys Rep. 2024;1068(2):1–66. doi:10.1016/j.physrep.2024.03.002. [Google Scholar] [CrossRef]

47. Bärtschi A, Eidenbenz S. Grover mixers for QAOA: shifting complexity from mixer design to state preparation. arXiv:2006.00354. 2020. [Google Scholar]

48. Bridi GA, Marquezino FdL. Analytical results for the quantum alternating operator ansatz with Grover mixer. Phys Rev A. 2024;110(5):052409. doi:10.1103/physreva.110.052409. [Google Scholar] [CrossRef]

49. Xie N, Xu J, Chen T, Lee X, Saito Y, Asai N, et al. Performance upper bound of a Grover-mixer quantum alternating operator ansatz. Phys Rev A. 2025;111(1):012401. doi:10.1103/physreva.111.012401. [Google Scholar] [CrossRef]

50. Stoudenmire EM, Waintal X. Opening the black box inside Grover’s algorithm. Phys Rev X. 2024;14(4):041029. doi:10.1103/physrevx.14.041029. [Google Scholar] [CrossRef]

51. Babbush R, McClean JR, Newman M, Gidney C, Boixo S, Neven H. Focus beyond quadratic speedups for error-corrected quantum advantage. PRX Quantum. 2021;2(1):010103. doi:10.1103/prxquantum.2.010103. [Google Scholar] [CrossRef]

52. Hoefler T, Haener T, Troyer M. Disentangling hype from practicality: on realistically achieving quantum advantage. arXiv:2307.00523. 2023. [Google Scholar]

53. Gilliam A, Venci C, Muralidharan S, Dorum V, May E, Narasimhan R, et al. Foundational patterns for efficient quantum computing. arXiv:1907.11513v3. 2021. [Google Scholar]

54. IBMQ. Error suppression and mitigation techniques. [cited 2024 Nov 5]. Available from: https://docs.quantum.ibm.com/guides/error-mitigation-and-suppression-techniques. [Google Scholar]

55. IonQ. Debiasing and sharpening. [cited 2024 Nov 11]. Available from: https://ionq.com/resources/debiasing-and-sharpening. [Google Scholar]

56. Wang Y, Hu Z, Sanders BC, Kais S. Qudits and high-dimensional quantum computing. Front Phys. 2020;8:589504. doi:10.3389/fphy.2020.589504. [Google Scholar] [CrossRef]

57. Shukla A, Vedula P. An efficient implementation of a quantum search algorithm for arbitrary N. Eur Phys J Plus. 2025;140(6):575. doi:10.1140/epjp/s13360-025-06518-3. [Google Scholar] [CrossRef]

58. Long GL, Li YS, Zhang WL, Tu CC. Dominant gate imperfection in Grover’s quantum search algorithm. Phys Rev A. 2000;61(4):042305. doi:10.1103/physreva.61.042305. [Google Scholar] [CrossRef]

59. Hsieh JY, Li CM, Chuu DS. An improved phase error tolerance in a quantum search algorithm. Chin J Phys. 2004;42(5):585. doi:10.5772/14391. [Google Scholar] [CrossRef]

60. Romanelli A, Auyuanet A, Donangelo R. Quantum search with resonances. Phys A Stat Mech Its Appl. 2006;360(2):274–84. doi:10.1016/j.physa.2005.05.101. [Google Scholar] [CrossRef]

61. Barenco A, Bennett CH, Cleve R, DiVincenzo DP, Margolus N, Shor P, et al. Elementary gates for quantum computation. Phys Rev A. 1995;52(5):3457–67. doi:10.1103/physreva.52.3457. [Google Scholar] [PubMed] [CrossRef]

62. Schuch N. Implementation of quantum algorithms with Josephson charge qubits. Regensburg, Germany: Universität Regensburg; 2002. [Google Scholar]

63. McKay DC, Hincks I, Pritchett EJ, Carroll M, Govia LCG, Merkel ST. Benchmarking quantum processor performance at scale. arXiv:2311.05933. 2023. [Google Scholar]

64. Pelofske E, Bärtschi A, Eidenbenz S. Short-depth QAOA circuits and quantum annealing on higher-order ising models. npj Quantum Inf. 2024;10:30. doi:10.21203/rs.3.rs-3238348/v2. [Google Scholar] [CrossRef]

65. Google Quantum AI. Suppressing quantum errors by scaling a surface code logical qubit. Nature. 2023;614(7949):676–81. [Google Scholar]

66. Google Quantum AI and Collaborators. Quantum error correction below the surface code threshold. Nature. 2025;638(8052):920–6. [Google Scholar]

67. He Y, Luo M, Zhang E, Wang HK, Wang XF. Decompositions of n-qubit Toffoli gates with linear circuit complexity. Int J Theor Phys. 2017;56(7):2350–61. doi:10.1007/s10773-017-3389-4. [Google Scholar] [CrossRef]

68. Dreier F, Fleckenstein C, Aigner G, Fellner M, Stahn R, Lanthaler M, et al. Connectivity-aware synthesis of quantum algorithms. In: Proceedings of the 2025 IEEE International Conference on Quantum Computing and Engineering (QCE); 2025 Aug 30–Sep 5; Albuquerque, NM, USA. p. 394–95. [Google Scholar]

69. Nie J, Zi W, Sun X. Quantum circuit for multi-qubit Toffoli gate with optimal resource. arXiv:2402.05053. 2024. [Google Scholar]

70. Silva JDS, da Silva AJ. Logarithmic depth decomposition of approximate multi-controlled single-qubit gates without ancilla qubits. arXiv:2507.00400. 2025. [Google Scholar]

71. Maslov D. Advantages of using relative-phase Toffoli gates with an application to multiple control Toffoli optimization. Phys Rev A. 2016;93(2):022311. doi:10.1103/physreva.93.022311. [Google Scholar] [CrossRef]

72. Kim Y, Morvan A, Nguyen LB, Naik RK, Jünger C, Chen L, et al. High-fidelity three-qubit iToffoli gate for fixed-frequency superconducting qubits. Nat Phys. 2022;18(7):783–8. doi:10.1038/s41567-022-01590-3. [Google Scholar] [CrossRef]

73. Warren CW, Fernández-Pendás J, Ahmed S, Abad T, Bengtsson A, Biznárová J, et al. Extensive characterization and implementation of a family of three-qubit gates at the coherence limit. npj Quantum Inf. 2023;9(1):44. doi:10.1038/s41534-023-00711-x. [Google Scholar] [CrossRef]

74. Katz O, Cetina M, Monroe C. Programmable N-body interactions with trapped ions. PRX Quantum. 2023;4:030311. [Google Scholar]

75. Harris NC, Carolan J, Bunandar D, Prabhu M, Hochberg M, Baehr-Jones T, et al. Linear programmable nanophotonic processors. Optica. 2018;5(12):1623–31. doi:10.1364/optica.5.001623. [Google Scholar] [CrossRef]

76. Mohit F, Guanzon J, McKinlay J, Weinhold TJ, Myers CR, Almeida MP, et al. Quantum mechanics can find a needle in a haystack every time. arXiv:2506.06435. 2025. [Google Scholar]